RASSI-SO

4.5 Performing RASSI-Spin-Orbit calculations.

• In the previous MS-CASPT2 step, the matrix of an effective Hamiltonian in the CAS basis was diagonalized (per symmetry block). The eigenvectors form the CAS ′ basis; the eigenvalues are the MS-CASPT2 energies. Both are transferred to RASSI through theJobMixfile. RASSIwill add the spin-orbit couplings ⟨CAS′i|HˆSO|CAS′j⟩, HˆSO being the spin-orbit coupling contribution to the full Hamiltonian, to the input (diagonal) MS-CASPT2 effective Hamiltonian matrix in the CAS ′ basis and will diagonalize it to get the wavefunctions and energies which contain SOC effects in addition to electron correlation.

• MSIN option in RASSI-SO calculations with RASSI.

Due to residual symmetry breakings, some MS-CASPT2 eigenvalues (MS-CASPT2 en- ergies) that should be degenerate are not, or are only almost-degenerate. The MSINput option provides the energies to do a “manual shift” of the input eigenvalues so that the close-but-non degenerate eigenvalues are substituted by their mean value; this avoids a lot of problems of identification of SO results later on. This is why the tail.$RassiStates.$RassiType.rassi.input file (see below) to be used with the MSIN option specifies which blocks should be degenerate. Note that in case of non- degeneracy beyond a threshold the input preparer will issue a warning. The MSIN option is not available in standard releases of MOLCAS; slight modifications of the RASSI code are needed.

4.5.1 Prepare and submit all RASSI-SO calculations using the program RASSI.

• Spin-orbit coupling will not mix gerade and ungerade states and two independent SO calculations could be done. However, if one is interested in obtaining the electric dipole transition moments between gerade and ungerade states, both sets should be included in the RASSI-SO calculations at once. We will do this. This means that all the wave functions we have calculated so far should be now invoked. Note that spin-orbit coupling will mix different spins and different irrep subspecies. You will find this info in the files below. 

• Go to directory inputs/RASSI.

• To prepare a RASSI-SO calculation for gerade and ungerade states using the MSIN option ($RassiStates=g.and.u, $RassiType=msin), prepare the head and tail files

head.g.and.u.msin.rassi.input 

tail.g.and.u.msin.rassi.input.

• • The head file specifies the number of JobIph (and JobMix) files stored in the vectors/CASPT2 directory that will be used by RASSI and the number of roots to be used from each. Note that spin-orbit coupling will mix different spins and different irrep subspecies. You will find in the template file:

NROF JOBIPHS= 35 4 4 4 2 2 2 2 1 2 2 2 4 4 4 6 1 3 6 6 6 5 5 5 6 2 1 6 6 6 5 5 5 6 2 1

• • The tail file specifies the labels that identify the JobIph (and JobMix) files, as well as degeneracy information. All 35 $Block.$Irrep cases we have calculated should appear. 

  •  You will find in the template file:

      RASSI LIST

      DEGENERATE 3

      3Tg4.3T1g.ms-4

      3Tg6.3T1g.ms-4

      3Tg7.3T1g.ms-4

      DEGENERATE 3

      3Tg4.3T2g.ms-2

      3Tg6.3T2g.ms-2

      3Tg7.3T2g.ms-2

      3Ag.3Eg.ms-2

      3Ag.3A2g.ms-1

      DEGENERATE 3

      1Tg4.1T1g.ms-2

      1Tg6.1T1g.ms-2

      1Tg7.1T1g.ms-2

      DEGENERATE 3        

      1Tg4.1T2g.ms-4

      1Tg6.1T2g.ms-4

      1Tg7.1T2g.ms-4

      1Ag.1Eg.ms-6

      1Ag.1A2g.ms-1

      1Ag.1A1g.ms-3

      DEGENERATE 3

      1Tu2.1T1u.ms-6

      1Tu3.1T1u.ms-6

      1Tu5.1T1u.ms-6

      DEGENERATE 3

      1Tu2.1T2u.ms-5

      1Tu3.1T2u.ms-5

      1Tu5.1T2u.ms-5

      1Au.1Eu.ms-6

      1Au.1A2u.ms-2

      1Au.1A1u.ms-1

      DEGENERATE 3

      3Tu2.3T1u.ms-6

      3Tu3.3T1u.ms-6

      3Tu5.3T1u.ms-6

      DEGENERATE 3

      3Tu2.3T2u.ms-5

      3Tu3.3T2u.ms-5

      3Tu5.3T2u.ms-5

      3Au.3Eu.ms-6

      3Au.3A2u.ms-2

      3Au.3A1u.ms-1

      END OF RASSI LIST

Any order is possible in the RASSI LIST, but once a given input order is chosen, the info in the head and tail files must match!

• • Between the geometry independent head and tail info, the RASSI input preparer will insert the MSINput information to finish the preparation of an input file for each geometry.

• Edit and run do.prepare.rassi.

This will invoke the RASSI input preparer shells/pre-post/prepare.rassi.subblocks.ksh, which will use the files:

• • head.g.and.u.msin.rassi.input

  • printouts/CASPT2/*.msin.energies 

  • tail.g.and.u.msin.rassi.input

and prepare the files:

• • $GeomLab.$RassiStates.msin.rassi.input, which are geometry-dependend input files for RASSI.

The MSIN information is the set of MS-CASPT2 energies corresponding to the CAS ′ wave functions in the RASSI list. These energies are averaged as demanded by the DEGENARATE n keyword used in the tail file. Should the non-degeneracy be larger than a threshold, the shell will issue a warning. For this reason, the directory where the MS-CASPT2 energies are, should be given as a parameter to the input preparer shell. Check some input files. Note the average values obtained from actual quasi-degenerate energies by confronting the values in the *.msin.energies or *.caspt2.output files. 

• • Submit the RASSI-SO calculations for all points in the dPr−F grid.

  • For that, in the shells directory: Edit j.calculation and uncomment only the “ksh run.rassi.msin.sh g.and.u $RassiDir” line. Invoke prepare.jobs to prepare all the job files, update do.submit, and invoke it to submit all the jobs.

• The spin-orbit states of Pr4+ can be calculated adapting the procedure we have just described for Pr3+ to the simpler 4f1 manifold. $RassiStates=f1 is a reasonable choice now. 

4.5.2 Analyze and plot the results of RASSI.

• • Go to directory printouts/RASSI.

  • Grep the rc= value in all printouts. It should be 0. Check some printouts.

  • Use the RASSI analyzer to analyze the RASSI-SO calculations (assign Oh′ double group irreps) and to partially prepare input for EFIT program.

  • Update and invoke do.analyze.rassi.

This script shell links the $Block.$Irrep.ms-caspt2.assignments files in the CASPT2 printouts directory ../CASPT2 and calls the analyzer shells/pre-post/analyze.rassi.ksh.

• • The following files will be used by the analyzer: 

    • *.msin.rassi.output files

    • ../../inputs/RASSI/*.msin.rassi.input files

    • *.ms-caspt2.assignments

    • SO.REDUCTION.TABLE file, which is self explanatory 

    • *.msin.rassi.fixed files (see below)

  • The following files are produced by the analyzer:

    • *.msin.rassi.summary files with the results of the analyses of the SO wave- functions read from RASSI printouts

    • *.msin.efit.dat file for EFIT

    • A file rassi.table.tex.dat with SO wavef unction analyses made at all geometries, preformatted for a later use to produce a LATEX table.

  • The *.summary files must be checked very carefully. They contain the spin-orbit vectors (states) specifying their spin-orbit-free components. On the basis of this table and using SO.REDUCTION.TABLE, the analyzer tries to assign each SO state to a double-group Oh′ irrep in this case. Once it has done that it looks for broken degeneracies larger than a threshold (300 cm−1) and issues warnings. This must be checked! After that, it counts and prints out the total number of SO levels (called terms in the analyzer) obtained. Again, this must be checked because if it is different from what is expected this could reveal wrong assignments which would demand the user sets the assignment of particular states in the *.fixed file and rerun the analysis again. This happens for several geometries in this case: 

    • The total number of SO terms and states we should get for Pr3+ can easily be calculated with the excel file in printouts/RASSI/SO states.xlsx or similar. Here is a copy of the g.and.u case:

Spin-orbit states from the f(2) manifold

                      A1g A2g Eg  T1g T2g

SF Terms No. SF Terms

1A1g     3             1   0   0   0   0

1A2g     1             0   1   0   0   0

1Eg      3             0   0   1   0   0

1T1g     2             0   0   0   1   0

1T2g     4             0   0   0   0   1

3A1g     0             0   0   0   1   0

3A2g     1             0   0   0   0   1

3Eg      1             0   0   0   1   1

3T1g     4             1   0   1   1   1

3T2g     2             0   1   1   1   1

                       7   3   9   9  12 <-- to be checked in *.summary files

Total   49                            91

Spin-orbit states from the f(1)phi(1) manifolds:

                      A1u A2u Eu  T1u T2u

SF Terms No. Terms

1A1u     1             1   0   0   0   0

1A2u     2             0   1   0   0   0

1Eu      3             0   0   1   0   0

1T1u     6             0   0   0   1   0

1T2u     5             0   0   0   0   1

3A1u     1             0   0   0   1   0

3A2u     2             0   0   0   0   1

3Eu      3             0   0   0   1   1

3T1u     6             1   0   1   1   1

3T2u     5             0   1   1   1   1

                       7   7  14  21  21 <-- to be checked in *.summary files

Total   84                           168

g.and.u 133                          259 states <-- to be checked in *.summary files

                                     110 levels (or SO terms)

• • • After the first run of do.analyze.rassi, we see that the number of expected SO terms for each irrep is correct at the dPr−F distances 2.100, 2.150, and 2.200 A, but wrong at higher distances. Checking the *.summary file at 2.250 A, for instance, we see that SO states 73, 74, and 75 have almost pure 3T1g character (***** appears in the place of a 1.0000 because of a format issue, but this is not a problem). According to SO.REDUCTION.TABLE, four Oh′ double-group irreps could have contributions from 3T1g:

           3T1g (9 states)  -->  A1g, Eg, T1g, T2g  (9 states)

so that the analyzer cannot solve the assignment. Then, the analyzer writes down “unknown” as the first option and A1g as the second. It adopts the second option and the number of A1g states is higher than expected by 3, while the number of T1g states is short by 1. At shorter distances there is a tiny contribution from 1T1g which allows the analyzer to assign the three states as the components of the T1g term because of its reduction table: 

1T1g (3 states) --> T1g (3 states)

• • Having in mind that the three states appear to be degenerate (25127 cm−1), and that they seem to correspond clearly with the same states at shorter Pr-F distances, the correct assignment should be T1g. Then, we set the assignments accordingly in file *.2.250.*.fixed and similarly in the *.fixed files of other geometries, and run do.analyze.rassi again.

  • Something similar happens with states 39, 40, and 41 at distances 2.550 and 2.600 A. At 2.450 A, some broken symmetry at the crossing point between two T1u and T2u terms difficults the automatic identification of states 181, 182, 183 as T2u and 184, 185, 186 as T1u. At 2.600 A, some broken symmetry prevents the automatic identification of state 42 as A2g.

  • After all these “fixed” assignments, the numbers of SO terms provided by the analyzer are all right. Further checking of the *.summary files does not show signs of further troubles. 

  • The analysis of the RASSI-SO printouts for Pr4+ is much simpler and goes very smooth. The number of SO states and levels which should be obtained have been worked out for Ce3+ in Sec. 3.5.2. 

4.5.3 Plot and fit the RASSI-SO results.

• • Use EFIT to fit the RASSI-SO energy curves and to prepare input to plot them.

    • Go to the results directory.

    • Update msin.efit.inp.head and prep.msin.efit files and run the latter. Check file msin.efit.out.

    • To plot the RASSI-SO energy curves, run XMGRACE, open Pr.RAS.PT2.SO.agr, and import msin.curves.txt in the window of the RASSI-SO results using colors and order suggested in README.agr.

      • The plot with all CASSCF/MS-CASPT2/RASSI-SO results shows the effects of dynamic correlation and SO coupling progressively.

      • Another interesting plot is that of SO results with all the terms together in one graph and separated according to Oh′ double-group irreps in other graphs.

      • Use colors in the plots to show the parentage of states with either atomic SO terms or crystal field splitting of the 4f5d configuration using the analyses described in Sec. 4.5.6. 

    • Edit msin.efit.inp and run “efit.ksh msin” until the fitting of the energy curves has a good quality. The data in msin.efit.out immediately after the keyword SUMMARY will be later read to make a LATEX table in a report. 

It is important to check the shape of numerical energy curves: interaction between states of the same symmetry but different configurational character will produce avoided crossings which may result in very anharmonic (even doble-minimum) curves. This should guide the fitting and should serve to point out where the spectroscopic con- stants it produces are uncertain (as indicated usually by footnotes in the tables). 

4.5.4 Calculate absorption and emission electric dipole oscillator strengths and radiative lifetimes.

In the RASSI calculation for gerade and ungerade states together, the program calculates and prints out the components of the electric dipole transition moment between all states above a threshold, as well as the Einstein emission coefficients. The correct combination of this data gives the absorption and emission electric dipole oscillator strengths and the radiative lifetimes.

• • Go to the printouts/RASSI/OS directory.

  • Update un run do.oscillator.strength.rassi.

  • Given a geometry to calculate the oscillator strengths ($GeomLab), this script shell links the necessary

  • ../*.$GeomLab.*.msin.rassi.output and

  • ../*.$GeomLab.*.msin.rassi.summary files

  • and calls shells/pre-post/oscillator.strength.rassi.ksh.

  • The script shell oscillator.strength.rassi.ksh calculates the total absorption or emission oscillator strengths associated with the electronic transition from one state of the initial SO level to all states of the final SO level. The parameter $initL specifies the initial SO level; its sign sets a flag for:

         initL   = 0    total oscillator strengths are not calculated; useful to check the SO level list and input data

         initL   > 0    absorptions from $initL to all possible upper levels

         initL   < 0    emissions from $initL to all possible lower levels; 

                            emission lifetime of this state will also be calculated using the Einstein coefficients;

It produces the files

• • *.msin.rassi.OscStr.$initL, with absolute and relative oscillator strenghts (and radiative lifetime for emission), preceeded by the table of oscillator strength values between individual SO states read from the RASSI printout.

  • ABSORPTION.OS or EMISSION.OS files, which will be read later to make a LATEX table in a report.

4.5.5 Make a basic report with the RASSI-SO results.

• Make the body of a LATEX table that combines RASSI-SO spectroscopic constants calculated with EFIT with the analyses of the SO wave functions (made at a given geometry) and, optionally, with the absorption and/or emission oscillator strenghts (made at their own given geometries).

  • In the results directory, run do.SO.table. It will use the files:

    • msin.efit.out, with RASSI-SO spectroscopic constants

    • printouts/RASSI/rassi.table.tex.dat, with the SO wavefunction analyses 

    • printouts/RASSI/OS/ABSORPTION.OS and 

    • printouts/RASSI/OS/EMISSION.OS, with oscillator strengths

It will produce:

• SO.table.tex.dat, with the core of a LATEX table

• SO.banda.dat, with the basic data to plot absorption and emission spectra profiles with the program BANDA-LS.

• Make the report.

  • In the results/latex directory, issue make. It will make the report in report.pdf, which can be visualized.

4.5.6 Analysis of SO terms.

• File printouts/RASSI/rassi.table.tex.dat contains the analyses of the SO terms for each geometry analyzed and Table II of the report of the RASSI-SO results contains the analyses at one given geometry, usually chosen to be the closest to the ground state and emitting state minima. These analyses should help to find meaningul parentage of SO levels which eventually could serve to classify them beyond their Oh′ double-group irreps.

For the 4f2 manifold, consider the following:

• • The weights of MS-CASPT2 terms will take you to the table I of MS-CASPT2 results where the parentage with atomic SF terms should have been worked out.

  • The splitting of atomic J quantum number in Oh field (see table below).

An example: The analysis of the lowest four SO terms gives (dPr−F=2.600 ̊A):

         ----------- 3H --------------

         1.3T1g 2.3T1g 1.3T2g 1.3Eg

1.T2g    67                   26 (in%) 

1.A1g    27     70

1.Eg     27     27     43

1.T1g    16     37     32     13

which ultimately allows to associate them with SF atomic 3H. In addition, atomic 3H can lead to J= 4, 5, and 6: 3H4, 3H5, 3H6 spin-orbit terms which in an octahedral field transform as follows: (this auxiliary table is in SO states.xlsx file)

J  A1 A2 E  T1 T2

0  1

1           1

2        1     1

3     1     1  1

4  1  1     1  1   <<< which allows us to associate the four

5        1  2  1       states mentiond with the atomic J=4 of 3H: 3H4

6  1  1  1  1  2

Similarly, the next

         2Eg, 2T1g, 2T2g, 3T1g                           are associated with 3H, J=5: 3H5,

and    2A1g, 3T2g, 3Eg, 4T1g, 4T2g, 1A2g    are associated with 3H, J=6: 3H6, 

and so on. In some cases, association is not possible and some SO states correspond to mixings.

For the higher ungerade states, we can also play with

• the SF components whose 4f5deg , 4f5dt2g , or 4fφa1g character has been worked out above.

• the trends of bond lengths:

4f5deg    <    4f2                < 4f5dt2g

Pr4+         <   Pr3+ 4fφITE   < other states of Pr3+

• • the shape of the curves, which can be either regular (and common to all states of the same dominant configuration) or correspond to avoided crossings, which indicate that states of different dominant configurations and same symmetry avoid crossings and change configurational character from one side to the other of the crossing,

  • the match with the number of SO terms one can expect from each configuration:

Spin-orbit states for f(1)deg(1) states only

             A1u A2u Eu  T1u T2u

Terms No.

1A1u    0    1   0   0   0   0

1A2u    0    0   1   0   0   0

1Eu     1    0   0   1   0   0

1T1u    2    0   0   0   1   0

1T2u    2    0   0   0   0   1

3A1u    0    0   0   0   1   0

3A2u    0    0   0   0   0   1

3Eu     1    0   0   0   1   1

3T1u    2    1   0   1   1   1

3T2u    2    0   1   1   1   1

             2   2   5   7   7  <<<---

Spin-orbit states for f(1)dt2g(1) states only

             A1u A2u Eu  T1u T2u

Terms No.

1A1u    1    1   0   0   0   0

1A2u    1    0   1   0   0   0

1Eu     2    0   0   1   0   0

1T1u    3    0   0   0   1   0

1T2u    2    0   0   0   0   1

3A1u    1    0   0   0   1   0

3A2u    1    0   0   0   0   1

3Eu     2    0   0   0   1   1

3T1u    3    1   0   1   1   1

3T2u    2    0   1   1   1   1

             4   3   7  11  10  <<<---

Spin-orbit states for f(1)a1g(1) states only

             A1u A2u Eu  T1u T2u

Terms No.

1A1u    0    1   0   0   0   0

1A2u    1    0   1   0   0   0

1Eu     0    0   0   1   0   0

1T1u    1    0   0   0   1   0

1T2u    1    0   0   0   0   1

3A1u    0    0   0   0   1   0  

3A2u    1    0   0   0   0   1

3Eu     0    0   0   0   1   1

3T1u    1    1   0   1   1   1

3T2u    1    0   1   1   1   1

             1   2   2   3   4  <<<--- 

4.5.7 Plot absorption and emission spectra profiles.

• • Go to results/spectrum-profile.

  • Update abs*.inp and/or emi*.inp files by inserting results/SO.banda.dat and changing the parameters of the spectrum simulations.

  • Update and run do.spectrum-profile.

  • Run XMGRACE, read spectrum-profile.agr, and import the appropriate *.dat files.