7.39. The Multireference Correlation Module

7.39.1. General Description

A number of uncontracted multireference approaches are implemented in ORCA and reside in the orca_mrci module. All of these approaches start with a reference wavefunction that consists of multiple configurations (orbital occupation patterns). The reference wavefunction defined in the ref subblock can be a complete active space (CAS), restricted active space (RAS) or an arbitrary list of configurations. The total wavefunction is constructed by considering single and double excitations out of the reference configurations. These excited configurations are then used to generate configuration state functions (CSF) that have the proper spin and spatial symmetry. The number of wavefunction parameters rapidly grows with the number of reference functions. The orca_mrci module features a set of truncation criteria (TSel, TPre, TNat) that help to reduce the number of wavefunction parameters. Furthermore, by default, the program only considers reference configurations that already have the target spin and spatial symmetry. There are situations, where this is undesired and the restrictions can be lifted with the keyword rejectinvalidrefs false. For more information on the theory, the program module as well as its usage we recommend the review article by Neese et al.[628]. A tutorial type introduction to the subject is presented in chapter The Multireference Correlation Module of the manual and more examples in the CASSCF tutorial. The detailed documentation of all features of the MR-CI and MR-PT module is somewhat premature and at this point only a summary of keywords is given below. A thorough description of all technical and theoretical subtleties must wait for a later version of the manual.

The overall scaling of uncontracted approaches is steep. Hence, the methodology is restricted to small reference spaces and small molecules in general. Note that all integrals must be kept in memory! Internally contracted multireference approaches such as NEVPT2 do not share these bottlenecks. Aside from NEVPT2, ORCA features a fully internally contracted MRCI (FIC-MRCI) that resides in the orca_autoci module. For more details on the FIC-MRCI we refer to section CI methods using generated code.

%mrci
       # -----------------------------------------------------------
       # Orbital selection
       # NOTE: The orbitals are used as supplied. Thus, the ORDER of
       # orbitals is critical. Say you have
       #   nact    electrons in the active space
       #   nint    electrons in the internal space
       #   nfrozen electrons
       # * The first nfrozen/2 orbitals will not be included in the CI
       # * The next nint/2 orbitals will be doubly occupied in all
       #   references
       # * the nact electrons are distributed over the,say, mact
       #   orbitals according to the active space definitions.
       # The remaining orbitals are external.
       # IT IS YOUR RESPONSIBILITY THAT THE ORBITAL ORDERING MAKES
       # SENSE!
       # A sensible two-step procedure is:
       #  * generate some orbitals and LOOK AT THEM. Decide which ones
       #    to include in the CI.
       #  * re-read these orbitals with ! MORead NoIter. Perhaps use
       #    the "rotate" feature to reorder the MOs
       #    Then jump right into the CI which is defined in this se-
       #    cond job
       #
       # NOTE: the MRCI module respects the %method FrozenCore settings
       # -----------------------------------------------------------
       Loc  0,0,0
           # Localize orbitals in the internal (first flag), active
           # (second flag) and external space (third flag).
       UseIVOs  false
           # Use improved virtual orbitals in the CI

       # ---------------------------------
       # Method selection
       # ---------------------------------
       CIType   MRCI      # Multireference CI (default)
                MRDDCI1   # Difference dedicated CI 1-degree of freedom
                MRDDCI2   # Difference dedicated CI 2-degrees of freedom
                MRDDCI3   # Difference dedicated CI 3-degrees of freedom
                MRACPF    # Average coupled-pair functional
                MRACPF2   # Modified version of ACPF
                MRACPF2a  # A slightly modified version of ACPF-2a
                MRAQCC    # Average quadratic coupled-cluster
                MRCEPA_R  # Multireference CEPA due to Ruttink
                MRCEPA_0  # CEPA-0 approximation
                SORCI     # Spectroscopy oriented CI
                SORCP     # Spectroscopy oriented couplet pair approx.
                MRMP2     # Multireference Moeller-Plesset at second order
                MRMP3     # Multireference Moeller-Plesset at third order
                MRMP4     # Multireference Moeller-Plesset at fourth order
                          # but keeping only singles and doubles relative to
                          # the reference configurations.

       # ---------------------------------
       # Selection thresholds
       # ---------------------------------
       Tsel    1e-6   # Selection threshold for inclusion in the CI based
                      # 2nd order MP perturbation theory <0|H|I>/DE(MP)
       Tpre    1e-4   # Selection of configurations in the reference space
                      # after the initial diagonalization of the reference
                      # space only configurations with a weight large>Tpre
                      # to any root are included
       AllSingles   false
                      # include ALL SINGLES in the CI. Default is now TRUE!!!

       # perturbative estimate of the effect of the rejected configurations
       EunselOpt  0   # no correction
                  1   # based on the overlap with the 0th order wavefunction
                  2   # calculation with the relaxed reference space
                      # coefficients. This is the most accurate and only
                      # slightly more expensive

       # For CIType=MRCI,MRDDCI and SORCI the approximate correction for
       # higher excitations
       DavidsonOpt  Davidson1  # default
                    Davidson2  # modified version
                    Siegbahn   # Siegbahn's approximation
                    Pople      # Pople's approximation

       # For MRACPF,MRACPF2,MRAQCC and SORCP
       NelCorr      0
          # Number of electrons used for computing the average coupled-
          # pair correction.
          # =0 : set equal to ALL electrons in the CI
          # =-1: set equal to all ACTIVE SPACE electrons
          # =-2: set equal to ACTIVE SPACE electrons IF inactive doubles
          #      are excluded (as in MRDDCI)
          # >0 : set equal to user defined input value
       LinearResponse false
          # Use ground state correlation energy to compute the shift for
          # higher roots (not recommended)

       # ---------------------------------
       # Natural Orbital Iterations
       # ---------------------------------
       NatOrbIters    0 # default
          # number of average natural orbital iterations
       Tnat        1e-4
          # cutoff of natural orbitals. NOs with an occupation number less
          # then Tnat will not be included in the next iteration
          # Also, orbitals with occupation number closer than Tnat to 2.0
          # will be frozen in the next iteration
       Tnat2       -1
          # if chosen >0 then Tnat2 is the threshold for freezing the
          # almost doubly occupied orbitals. Otherwise it is set equal
          # to Tnat

       # ----------------------------------
       # Additional flags and algorithmic
       # details
       # ----------------------------------
       PrintLevel   2  # default. Values between 1 and 4 are possible

       DoDDCIMP2  false
          # for DDCI calculations: if set to true the program computes
          # a MP2 like correction for the effect of inactive double
          # excitations which are not explicitly included in the CI. This
          # is necessary if you compare molecules at different geometries
          # or compute potential energy surfaces.
       # ----------------------------------
       # The SORCP model
       # ----------------------------------
       CIType_in  # First step CIType
       CIType_fi  # Second step CIType
       Exc_in     # First step excitation scheme
       Exc_fi     # Second step excitation scheme
       Tsel_in    # First step Tsel
       Tsel_fi    # Second step Tsel
       Tpre_in    # First step Tsel
       Tpre_fi    # Second step Tpre

         # Thus, the SORCI model corresponds to CIType=SORCP with
         # CIType_in MRCI  CIType_fi MRCI
         # Exc_in    DDCI2 Cexc_fi   DDCI3
         # Tsel_in   1e-5  Tsel_fi   1e-5
         # Tpre_in   1e-2  Tpre_fi   1e-2

       # ----------------------------------
       # Multirerence perturbation theory
       # ----------------------------------
       MRPT_b  0.02   # Intruder state avoidance PT after Hirao (default 0.0)
                      # with this flag individual intruders are shifted away to
                      # to some extent from the reference space
       MRPT_shift 0.3 # Level shift introduced by Roos which shifts the entire
                      # excited manifold away in order to avoid intruder states.
                      # A correction is applied afterwards but results do depend
                      # on this (arbitrary) value to some extent.
       H0Opt  projected  # use an off-diagonal definition of H0
              Diagonal   # use a diagonal definition of H0 (much faster but maybe
                         # a little less reliable
       Partitioning MP   # Moeller plesset partitioning
                    EN   # Epstein-Nesbet partitioning (not recommended)
       Fopt  Standard    # Standard definition of MR Fock operators
             G3          # uses Anderson's g3 correction also used in CASPT2


       #---------------------------------------
       # restrict reference configurations
       #---------------------------------------
       RejectInvalidRefs true # by default reference CSFs are restricted
                              # to target spin and spatial symmetry

       # ======================================
       # Definitions of blocks of the CI Matrix
       # ======================================
       NewBlock 2 *    # generate a Block with doublet(=2) multiplicity
         Nroots  1     # number of roots to be generated
         Excitations cis   # CI with single excitations
                     cid   # CI with double excitations
                     cisd  # CI with single and double excitations
                     ddci1 # DDCI list with one degree of freedom
                     ddci2 # DDCI list with two degrees of freedom
                     ddci3 # DDCI list with three degrees of freedom
         Flags[_class_] 0 or 1
                  # Turn excitation classes on or off individually
                  # ``s'' stands for any SOMO, ``i'',``j'' for internal orbitals and
                  # ``a'',``b'' for external orbitals
                  #  Singles _class_ = ss, sa, is, ia
                  #  Doubles _class_ = ijss, ijsa, ijab,
                  #                    isss, issa, isab,
                  #                    ssss, sssa, ssab
                  # ``Flags'' takes priority over ``Excitations''. In fact ``Excitations''
                  # does nothing but to set ``Flags''. So, you can use ``Excitations''
                  # to provide initial values for ``Flags'' and then modify them
                  # with subsequent ``Flags'' assignments
         refs
               #
               # First choice - complete active space
               #
               CAS(nel,norb)  # CAS-CI reference with nel electrons in
                              # Norb orbitals
               #
               # Second choice - restricted active space
               #
               RAS(nel: m1 h/ m2 / m3 p)
                              # RAS-reference with nel electrons
                     # m1= number orbitals in RAS-1
                     # h = max. number of holes in RAS-1
                     # m2= number of orbitals in RAS-2 (any number of
                     #     electrons or holes)
                     # m3= number of orbitals in RAS-3
                     # p = max. number of particles in RAS-3
               #
               # Third choice - individually defined configurations
               #
                \{ 2 0 1 0\}
                \{ 1 1 1 0\}
                  etc.
                  # define as many configurations as you want. Doubly occupied MOs
                  # singly occupied MOs and empty MOs. Important notes:
                  #  a) the number of electrons must be the same in all references
                  #  b) the number of orbitals is determined from the number of
                  #     definitions. Thus, in the example above we have three active
                  #     electrons and four active orbitals despite the fact that the
                  #     highest orbital is not occupied in any reference.
                  #  The program determines the internal, active and external spaces
                  #  automatically from the number of active electrons and orbitals
               end
         end
         # there can be as many blocks as you want!!!

       # ----------------------------------
       # Density matrix generation flags
       #  First Key= State densities <I|D|I>
       #     =0: none
       #     =1: Ground state only (lowest root of all blocks; Electron only)
       #     =2: Ground state only (Electron and spin density)
       #     =3: Lowest root from each block (Electron density)
       #     =4: Lowest root from each block (Electron and spin density)
       #     =5: All states (Electron density)
       #     =6: All states (Electron and spin density)
       #  Second Key= Transition densities <I|D|J>
       #     needed for all transition intensities, g-tensor etc
       #     =0: none
       #     =1: from the ground state into all excited states (el)
       #     =2: from the ground state into all excited states (el+spin)
       #     =3: from all lowest states into all excited states (el)
       #     =4: from all lowest states into all excited states (el+spin)
       #     =5: all state pairs (el)
       #     =6: all state pairs (el+spin)
       # Note that for perturbation theory the density is computed as
       # an expectation value over the first (second) order wavefunction.
       # which is renormalized for this purpose
       # ----------------------------------
       Densities 1,1
       
       # ----------------------------------
       # Complete printing of the wavefunction 
       # ----------------------------------
       PrintWF    1      # CFG based printing (default)
                  det    # Determinant based wavefunction printing
       TPrintWF   3e-3   # Threshold for the printing of the CFGs/Dets

       # ----------------------------------
       # Algorithm for the solver
       # ----------------------------------
       Solver   Diag   # Davidson like solver
                DIIS   # DIIS like solver
         # both solvers have their pros and cons. The DIIS may converge
         # better or use less time since it only recomputes the vectors that
         # have not yet converged; The DIIS may be less sensitive to root flipping
         # effects but occasionally it converges poorly and states of the same
         # symmetry are occasionally a little problematic
         # For perturbation theory DIIS is always used.
       # For both solvers
       MaxIter   100  # the maximum number of iterations
       Etol      1e-6 # convergence tolerance for energies in Eh
       Rtol      1e-6 # convergence tolerance for residual

       # For Solver=Diag (Davidson solver)
       Otol      1e-16 # Orthogonality threshold for Schmidt process
       NGuessMat 512   # Dimension of the guess matrix 512x512
                       # be used to compute the initial guess of the actual MRCI calculation
       NGuessMatRefCI 512 # Dimension of the guess matrix 
                          # for the reference CI
        
       MaxDim       3  # Davidson expansion space = MaxDim * NRoots
       # For the Solver=DIIS. Particularly recommended for anything else but
       # straightforward CI and also for calculations in direct2 mode!
       MaxDIIS       5 # Maximum number of old vectors to be used in DIIS
       RelaxRefs  true # Relax reference space coefficients in the CI or
                       # freeze them to their zeroth order values
       LevelShift  0.4 # Level Shift for stabilizing the DIIS procedure
       
       # ----------------------------------
       # RI Approximation
       # ----------------------------------
       IntMode RITrafo #Use RI integrals
               FullTrafo #No RI (default)
       
       # ----------------------------------
       # Integral storage, memory and files
       # ----------------------------------
       IntStorage    FloatVals
                     DoubleVals  (default)
          # store integrals with float (4 byte) or double (8 byte)
          # accuracy in main memory
       FourIndexInts false (default)
                     True
          # Store ALL four index integrals over Mos in main memory
          # only possible for relatively small systems, perhaps up
          # to 150-200 MOs included in the CI
       MaxMemInt     256
          # Maximum amount of core memory devoted to the storage of
          # integrals. If NOT all three index integrals fit into main
          # memory the program fails
       MaxMemVec      16
          # Maximum amount of memory used for section of the trial and
          # sigma vectors. This is not a particularly critical variable
       KeepFiles    false
          # Keep integrals and CI program input file (.mrciinp). Then
          # you can manually edit the .mrciinp file which is a standard
          # ASCII file and run the MRCI program directly. The only thing
          # you cannot change is the orbital window.
       end

7.39.2. Properties Calculation Using the SOC Submodule

7.39.2.1. Zero-Field Splitting

The spin-orbit coupling (SOC) and spin-spin coupling (SSC) contributions to the zero-field splitting (ZFS) can be calculated very accurately using a wavefunction obtained from a multiconfigurational calculation of a multi-reference type such as CASSCF, MRCI, or MRPT as is described in QDPT Magnetic Properties Section Magnetic properties through Quasi Degenerate Perturbation Theory.

# In case that you want to run QDPT-SOC calculation with manually
#adjusted diagonal energies you can copy the following part into
#the %mrci soc block
#and modify it as needed(energies are given in
#wavenumbers relative to the lowest state)
# NOTE: It is YOUR responsibility to make sure that the CAS-CI state
#that you may want to dress with these energies correlate properly
#with the energies printed here. The order of states or even the
#identity of states may change with and without inclusion of
#dynamic correlation In the case that dynamic correlation strongly
#mixes different CAS-CI states there may not even be a proper
#correlation!
#
     EDiag[  0]        0.00  # root   0 of block 0
     EDiag[  1]    48328.40  # root   1 of block 0
     EDiag[  2]    48328.40  # root   2 of block 0
     EDiag[  3]    49334.96  # root   3 of block 0
     EDiag[  4]     7763.59  # root   0 of block 1
     EDiag[  5]     7763.59  # root   1 of block 1
     EDiag[  6]    11898.46  # root   2 of block 1
     EDiag[  7]    46754.23  # root   3 of block 1

Those transition energies can be substituted by a more accurate energies provided in the input file as follows:

%soc
dosoc true
dossc true
     EDiag[  0]        0.00  # root   0 of block 0
     EDiag[  1]    48328.40  # root   1 of block 0
     EDiag[  2]    48328.40  # root   2 of block 0
     EDiag[  3]    49334.96  # root   3 of block 0
     EDiag[  4]     7763.59  # root   0 of block 1
     EDiag[  5]     7763.59  # root   1 of block 1
     EDiag[  6]    11898.46  # root   2 of block 1
     EDiag[  7]    46754.23  # root   3 of block 1
end

Accurate diagonal energies generally improve the accuracy of the SOC and SSC splittings.

7.39.2.2. Local Zero-Field Splitting

The submodule can also be used to calculate the local ZFS splitting parameters of atomic centers. The method, referred to as local complete active space configuration interaction (L-CASCI), can be used to separate into atomic contributions the SOC part of the total ZFS tensor. The rational behind it and additional details are described in the original publication [716]; below are listed only the steps required to reproduce the calculation for the dimer complex presented there.

1. The first step consists in obtaining the molecular orbitals that are going to be used in the configuration interaction (CI) procedure. A good set of orbitals can be obtained from a restricted open-shell spin-averaged Hartree-Fock (SAHF) calculation. The relevant part of the input is listed below:

! def2-tzvp keepfock 

% scf
    hftyp rohf
    rohf_case sahf
    rohf_numop 2
    rohf_nel[1] 9
    rohf_norb[1] 10
    end

For the present Mn(II)Mn(III) dimer there are a total of 9 electrons distributed into 10 d-orbitals.

2. Next, the molecular orbitals are localized using one of the implemented localization schemes. Below is the orca_loc input used in this case:

sahf.gbw
sahf.loc
0
200 # first of the 10 d-orbitals
209 # last of the 10 d-orbitals
128
0.000001
0.75
0.65
2

3. Following this, the localized orbitals are made locally canonical by block diagonalizing the Fock matrix using the orca_blockf utility.

orca_blockf sahf.fsv sahf.loc 200 204 205 209

The first two numbers define the range of molecular orbitals localized on one center; the last two are for the second center.

4. The recanonicalized orbitals stored in the sahf.loc file can be then used to calculate the SOC contribution to the local ZFS of the Mn(III) center using the following MRCI input:

! zora-def2-tzvp def2-tzvp/c zora
! nomulliken noloewdin
! moread noiter allowrhf
! moread

% mrci
    citype mrci
    tsel 0
    tpre 0
    intmode ritrafo
    solver diis
    soc
        intmode ritrafo
        dosoc true
        end
    newblock 10 *
        nroots 5
        excitations none
        refs
            #  Mn(II)    Mn(III)
            {1 1 1 1 1  1 1 1 1 0}
            {1 1 1 1 1  1 1 1 0 1}
            {1 1 1 1 1  1 1 0 1 1}
            {1 1 1 1 1  1 0 1 1 1}
            {1 1 1 1 1  0 1 1 1 1}
            end
        end
    newblock 8 *
        nroots 45
        excitations none
        refs
            #  Mn(II)    Mn(III)
            {1 1 1 1 1  2 1 1 0 0}
            {1 1 1 1 1  2 1 0 1 0}
            {1 1 1 1 1  2 1 0 0 1}
            {1 1 1 1 1  2 0 1 1 0}
            {1 1 1 1 1  2 0 1 0 1}
            {1 1 1 1 1  2 0 0 1 1}
            {1 1 1 1 1  1 2 1 0 0}
            {1 1 1 1 1  1 2 0 1 0}
            {1 1 1 1 1  1 2 0 0 1}
            {1 1 1 1 1  1 1 2 0 0}
            {1 1 1 1 1  1 1 1 1 0}
            {1 1 1 1 1  1 1 1 0 1}
            {1 1 1 1 1  1 1 0 2 0}
            {1 1 1 1 1  1 1 0 1 1}
            {1 1 1 1 1  1 1 0 0 2}
            {1 1 1 1 1  1 0 2 1 0}
            {1 1 1 1 1  1 0 2 0 1}
            {1 1 1 1 1  1 0 1 2 0}
            {1 1 1 1 1  1 0 1 1 1}
            {1 1 1 1 1  1 0 1 0 2}
            {1 1 1 1 1  1 0 0 2 1}
            {1 1 1 1 1  1 0 0 1 2}
            {1 1 1 1 1  0 2 1 1 0}
            {1 1 1 1 1  0 2 1 0 1}
            {1 1 1 1 1  0 2 0 1 1}
            {1 1 1 1 1  0 1 2 1 0}
            {1 1 1 1 1  0 1 2 0 1}
            {1 1 1 1 1  0 1 1 2 0}
            {1 1 1 1 1  0 1 1 1 1}
            {1 1 1 1 1  0 1 1 0 2}
            {1 1 1 1 1  0 1 0 2 1}
            {1 1 1 1 1  0 1 0 1 2}
            {1 1 1 1 1  0 0 2 1 1}
            {1 1 1 1 1  0 0 1 2 1}
            {1 1 1 1 1  0 0 1 1 2}
            end
        end
    end

5. The three second order ZFS components printed at the end of the calculation (Second order D-tensor: component 0, etc.) are scaled using the S value for the complex, which in this case is 4.5 (9 electrons \(\times\) 0.5). In order to obtain the correct local value of the ZFS, the three matrices have to be rescaled using the S value for Mn(III), which is to 2. Note that the three matrices have different scaling prefactors, and the dependence on S is not the same:

\(\mathbf{D}^{SOC-(0) } \propto \frac{1}{S^2}\)

\(\mathbf{D}^{SOC-(-1) } \propto \frac{1}{S(2S-1) }\)

\(\mathbf{D}^{SOC-(+1) } \propto \frac{1}{(S+1)(2S+1) }\)

These equations can be used to calculate the required prefactors. For example in the case of the SOC-(0) the prefactor is equal to:

\(\mathbf{D}_{\mathrm{Mn(III) }}^{SOC-(0) } = \frac{4.5^2}{2^2}\cdot\mathbf{D}_{\mathrm{dimer} }^{SOC-(0) } = 5.0625 \cdot \mathbf{D}_{\mathrm{dimer} }^{SOC-(0) }\)

The final step is to scale the two remaining matrices using the appropriate prefactors, sum all three of them up, diagonalize the resulting the matrix, and use its eigenvalues to calculate the D and E parameters. These represent the local ZFS parameters of the Mn(III) center.

7.39.2.3. Zero-Field Splitting from an excited Multiplet

For an excited state Multiplet the Calculationof ZFS can be requested by

Lowest eigenvalue of the SOC matrix:    -149.86223277 Eh
Energy stabilization:    -2.54512 cm-1
Eigenvalues:     cm-1         eV      Boltzmann populations at T =  300.000 K
   0:            0.00        0.0000       3.36e-01
   1:            2.37        0.0003       3.32e-01
   2:            2.37        0.0003       3.32e-01
   3:         7757.65        0.9618       2.33e-17
   4:         7757.66        0.9618       2.33e-17
   5:        11913.81        1.4771       5.15e-26
soc
 DTensor true
 IStates 3,4,5
end
*****************************************
 EXCITED STATE ZERO-FIELD SPLITTING:
*****************************************

--------------------------------------------
Computing Excited State D Tensors of
Excited State Multiplet Consisting of States :  3  4  5
--------------------------------------------
          0         4     
          1         5     
          2         0     
          3         1     
          4         0     
          5         2     


--------------------------------------------
           ZERO-FIELD SPLITTING
(2ND ORDER SPIN-ORBIT COUPLING CONTRIBUTION)
--------------------------------------------

D   =   -2.668445  cm-1
E/D =    0.000103


...

--------------------------------------------------------
                  ZERO-FIELD SPLITTING
EFFECTIVE HAMILTONIAN SOC CONTRIBUTION
--------------------------------------------------------

D   =   -2.674495  cm-1
E/D =    0.009610

...

7.39.2.4. g-Tensor

The orca_mrci program contains an option to calculate g-tensors using MRCI wavefunctions. For a system with an odd number of electrons, the doubly degenerate eigenvalues obtained from the QDPT procedure represent Kramers pairs, which are used to build the matrix elements of the total spin operator and the total angular momentum operator from the Zeeman Hamiltonian. Denoting \(\Psi\) as a solution and \(\bar{\Psi}\) as its Kramers partner and using matrix element notations

(7.240)\[\Phi_{11}^{k} =\left\langle\Psi \right|\hat{{L} }_{k} +g_{e} \hat{{S} }_{k} \left| \Psi \right\rangle,\, \Phi_{12}^{k} =\left\langle\Psi \right|\hat{{L} }_{k} +g_{e} \hat{{S} }_{k} \left| \bar{\Psi} \right\rangle,\, k=x,y,z \]

The elements of g-matrix are obtained as:

(7.241)\[g_{kz} =2\Phi_{11}^{k} ,\, g_{ky} =-2\Im\left({ \Phi_{12}^{k} } \right) ,\, g_{kx} =2\Re\left({ \Phi_{12}^{k} } \right)\]

Then, the true tensor G is built from g-matrices:

(7.242)\[G=gg^{T} \]

G is subjected further to diagonalization yielding positive eigenvalues, the square roots of which give the principal values of g-matrix.

(7.243)\[g_{xx} =\sqrt{ G_{xx} } ,\, g_{yy} =\sqrt{ G_{yy} } ,\, g_{zz} =\sqrt{ G_{zz} } \]

A typical mrci block of the input file for a g-tensor calculation should (e.g. for a S=3/2 problem) look as the following:

%mrci   ewin -4,1000
        citype mrci
        cimode direct2
        intmode fulltrafo
        solver diis
        etol 1e-8
        rtol 1e-8
        tsel 1e-6
        tpre 1e-5
        soc
          PrintLevel 2
          GTensor true      # make g-tensor calculations
          NDoubGTensor 2    # number of Kramers doublets to account
                            # for every pair a separate
                            # calculation is performed
        end
        newblock 4 *
          excitations cisd
          nroots 10
          refs  cas(7,5) end
        end
end

The result for the first Kramers pair is printed as follows:

--------------
KRAMERS PAIR 1
--------------

Matrix elements Re<1|S|1>   -0.072128   0.024511  -2.998843 
Matrix elements Re<1|S|2>   -0.001088   0.000366  -0.002010 
Matrix elements Im<1|S|2>   -0.000354  -0.001037  -0.000173 
Matrix elements Re<1|L|1>   -0.027067   0.009209  -1.123531 
Matrix elements Re<1|L|2>   -0.000031   0.000010  -0.000763 
Matrix elements Im<1|L|2>   -0.000006  -0.000011  -0.000065 

-------------------
ELECTRONIC G-MATRIX
-------------------

g-matrix:
-0.002240     0.000754    -0.005551 
0.000720     0.002100     0.000477 
-0.198556     0.067498    -8.251703 

g-factors:
0.002220     0.002222     8.254370 iso =     2.752937

g-shifts:
-2.000100    -2.000098     6.252051 iso =     0.750618

Eigenvectors:
0.057426     0.998060     0.024055 
0.998327    -0.057244    -0.008177 
0.006784    -0.024484     0.999677

Here for the \(L\) and \(S\) matrix elements indices 1 and 2 are assumed to denote Kramers partners, and three numbers in the first row stand for \(x, y, z\) contributions.

In addition the g-tensor is calculated within the Effective Hamiltonian formalism.

----------------------------------------------
	ELECTRONIC G-MATRIX FROM EFFECTIVE HAMILTONIAN
	----------------------------------------------
	
	g-matrix:
	1.978874    -0.000345     0.018908 
	-0.000345     1.977899    -0.006433 
	0.018879    -0.006418     2.763402 
	
	g-factors:
	1.977789     1.978477     2.763909 iso =     2.240058
	
	g-shifts:
	-0.024530    -0.023843     0.761590 iso =     0.237739
	
	Eigenvectors:
	0.288884     0.957062     0.024060 
	0.957364    -0.288770    -0.008181 
	0.000882    -0.025397     0.999677 
	
	# The g-factors are square roots of the eigenvalues of gT*g
	# Orientations are the eigenvectors of gT*g

Finally and only within the MRCI module the g-tensor is evaluated by using the Sum Over States formalism[609]:

---------------------------------------------------------------------------
	SUM OVER STATES CALCULATION OF THE SPIN HAMILTONIAN (for g and HFC tensors)
	---------------------------------------------------------------------------
	
	
	Ground state index        =     0
	Ground state multiplicity =     4
	Ground state spin density = P[  1]
	State =   1  <0|P|I>=   2  <0|Q|I>=  19
	State =   2  <0|P|I>=   3  <0|Q|I>=  27
	State =   3  <0|P|I>=   4  <0|Q|I>=  34
	State =   4  <0|P|I>=   5  <0|Q|I>=  40
	State =   5  <0|P|I>=   6  <0|Q|I>=  45
	State =   6  <0|P|I>=   7  <0|Q|I>=  49
	State =   7  <0|P|I>=   8  <0|Q|I>=  52
	State =   8  <0|P|I>=   9  <0|Q|I>=  54
	State =   9  <0|P|I>=  10  <0|Q|I>=  55
	
	Origin for angular momentum        ... ( -0.0006, -0.0010,  0.0021)
	Kinetic Energy                     ... done
	Relativistic mass correction       ... done
	Gauge correction                   ... done
	Angular momentum integrals         ... done
	Reading Spin-Orbit Integrals       ... done
	-----------------------
	MATRIX ELEMENT PRINTING
	-----------------------
	
	Energy differences (DE=EI-E0) and spin-orbit matrix elements (SO=<I|HSO|0>) are 
	printed in cm**-1. Orbital Zeeman matrix elements (L=<I|L|0>) are printed in au.
	
	State     DE         LX        LY        LZ           SOX         SOY         SOZ
	
	1    1349.3     0.0464   -0.0158    1.9264     -23.432       7.965    -974.312
	2   13026.2    -0.6596    0.6888    0.0214     337.028    -351.116     -10.966
	3   13615.1    -0.6961   -0.6514    0.0113     354.225     332.219      -5.736
	4   56686.3    -0.0053    0.0077    0.0971       1.794      -1.696     -36.786
	5   56954.2    -0.0516   -0.0048   -0.0042      28.211       5.821       1.459
	6   56994.0    -0.0418    0.0233   -0.0025      15.185      -2.144       1.145
	7   63371.5    -0.0211    0.0226    0.0078       3.833      -2.948      -2.724
	8   64176.0    -0.0652    0.0032   -0.0002      32.779       6.146       0.063
	9   74309.9    -0.0007    0.0032    0.0380       0.183      -1.058     -13.517
	
	-------------------
	ELECTRONIC G-MATRIX
	-------------------
	
	raw-matrix : 
	2.025533    -0.000738     0.021755  
	-0.000738     2.024537    -0.007389  
	0.021755    -0.007389     2.928943  
	
	g-factors: 
	2.024122     2.025363     2.929527 iso =     2.326338
	
	g-shifts:
	0.021803     0.023044     0.927208 iso =     0.324018
	
	Eigenvectors: 
	0.533896    -0.845208     0.024064  
	0.845530     0.533866    -0.008182  
	-0.005932     0.024715     0.999677  
	
	
	Euler angles w.r.t. molecular frame (degrees):
	-76.5038       1.4564    -161.2223
	
	-----------------------------
	CONTRIBUTIONS TO THE G-MATRIX
	-----------------------------
	
	Term                                g1           g2           g3
	--------------------------------------------------------------------
	Relativistic mass correction:   -0.0008220   -0.0008220   -0.0008220
	Gauge correction            :    0.0000000    0.0000000    0.0000000
	g(OZ/SOC)                   :    0.0226250    0.0238662    0.9280297
	
	State   1                		:    0.0000000   -0.0000000    0.9279829
	State   2                		:    0.0013767    0.0223913    0.0000000
	State   3                		:    0.0212332    0.0014408    0.0000000
	State   4                		:    0.0000000    0.0000004    0.0000418
	State   5                		:    0.0000074    0.0000099    0.0000001
	State   6                		:    0.0000002    0.0000078    0.0000001
	State   7                		:    0.0000000    0.0000015    0.0000002
	State   8                		:    0.0000076    0.0000144    0.0000000
	State   9                		:    0.0000000    0.0000000    0.0000046
	-----------------------------------------
	Total g-shifts              :    0.0218030    0.0230442    0.9272077
	
	
	# The g-factors are square roots of the eigenvalues of gT*g
	# Orientations are the eigenvectors of gT*g

Note that within the SOS formalism in addition to the second order (SOC) contributions the bilinear to the field terms: Relativistic mass correction and diamagnetic spin-orbit term (Gauge) are evaluated. As can be seen these corrections are rather negligible in comparison to the second order SOC contributions and most of the time can be safely omitted. Moreover further insight is obtained by printing the individual contribution of each excited state to the g-tensor. In the example above the first excited state contributes to the \(g_z\) component while the next two to both the \(g_x\) and \(g_y\) components, respectively.

So to summarize the g-tensor calculations in the framework of wavefunction based methods like MRCI and/or CASSCF can be evaluated:

  • via the QDPT approach within an individual Kramers doublet. This is valid analysis only for non-integer spin cases. In particular for systems with well isolated Kramers doublets where the EPR spectrum originates only from one Kramers doublet defined within the pseudo spin 1/2 formalism. This analysis has been proven useful in determining the sign of the ZFS and the electronic structure of the system under investigation.[547]

  • within the effective Hamiltonian approach. This is a valid analysis for all spin cases as it provides the principal g-values of the system under investigation evaluated in the molecular axis frame. These g-values can be directly compared with the experimentally determined ones.[421]

  • within the sum over states formalism (SOS). As above this analysis is valid for all spin cases and is only available via the MRCI module.

7.39.2.5. Magnetization and Magnetic Susceptibility

The MRCI and CASSCF modules of ORCA allow for the calculation of magnetization and magnetic susceptibility curves at different fields and temperatures by differentiation of the QDPT Hamiltonian with respect to the magnetic field. For magnetic susceptibility, calculations are performed in two ways when a static field different from zero is defined: (i) as the second derivative of energy with respect to the magnetic field and (ii) as the magnetization divided by the magnetic field. Although the first method corresponds to the definition of magnetic susceptibility, the second approach is widely used in the experimental determination of \(\chi*T\) curves. If the static field is low, both formulas tend to provide similar values.

The full list of keywords is presented below.

%mrci
      citype mrci
      newblock 3 *
      excitations none
      refs cas(2,7) end
      end
      soc
        dosoc true
        domagnetization true       # Calculate magnetization (def: false)
        dosusceptibility true      # Calculate susceptiblity (def: false)
        LebedevPrec            5   # Precision of the grid for different field
                                   # directions (meaningful values range from 1
                                   # (smallest) to 10 (largest))
        nPointsFStep           5   # number of steps for numerical differentiation
                                   # (def: 5, meaningful values are 3, 5 7 and 9)
        MAGFieldStep       100.0   # Size of field step for numerical differentiation
                                   # (def: 100 Gauss)
        MAGTemperatureMIN    4.0   # minimum temperature (K) for magnetization
        MAGTemperatureMAX    4.0   # maximum temperature (K) for magnetization
        MAGTemperatureNPoints  1   # number of temperature points for magnetization
        MAGFieldMIN          0.0   # minimum field (Gauss) for magnetization
        MAGFieldMAX      70000.0   # maximum field (Gauss) for magnetization
        MAGNpoints            15   # number of field points for magnetization
        SUSTempMIN           1.0   # minimum temperature (K) for susceptibility
        SUSTempMAX         300.0   # maximum temperature (K) for susceptibility
        SUSNPoints           300   # number of temperature points for susceptibility
        SUSStatFieldMIN      0.0   # minimum static field (Gauss) for susceptibility
        SUSStatFieldMAX      0.0   # maximum static field (Gauss) for susceptibility
        SUSStatFieldNPoints    1   # number of static fields for susceptibility
      end
end

The same keywords apply for CASSCF calculations in rel block (instead of soc in MRCI). Although different aspects of integration and grid precision can be modified through keywords, default values should provide an accurate description of both properties. Calculated magnetization and susceptibility are printed in .sus and .mag files, respectively and also in the output file.

-------------------------------------------------------------------------------
FIELD DEPENDENT MAGNETIZATION AND MEAN SUSCEPTIBILITY (chi=M/B)
-------------------------------------------------------------------------------
TEMPERATURE (K)   M. FIELD (Gauss)   MAGNETIZATION (B.M.)   chi*T (cm3*K/mol)
-------------------------------------------------------------------------------

4.00               0.00               0.000000                  inf
4.00            5000.00               0.350759             1.567189
4.00           10000.00               0.688804             1.538788
4.00           15000.00               1.003466             1.494496
4.00           20000.00               1.287480             1.438115
4.00           25000.00               1.537346             1.373773
4.00           30000.00               1.752841             1.305282
4.00           35000.00               1.936067             1.235764
4.00           40000.00               2.090450             1.167516
4.00           45000.00               2.219920             1.102067
4.00           50000.00               2.328368             1.040315
4.00           55000.00               2.419335             0.982690
4.00           60000.00               2.495883             0.929301
4.00           65000.00               2.560582             0.880052
4.00           70000.00               2.615538             0.834730
-----------------------------------------------------------
-----------------------------------------------------------
	TEMPERATURE DEPENDENT MAGNETIC SUSCEPTIBILITY
	-----------------------------------------------------------
	STATIC  FIELD    TEMPERATURE       chi*T (cm3*K/mol)
	(Gauss)           (K)          M/B          d2E/dB2
	-----------------------------------------------------------
	
	0.00       1.00           ----       1.576836
	0.00       2.00           ----       1.576910
	0.00       3.00           ----       1.576951
	0.00       4.00           ----       1.576988
	0.00       5.00           ----       1.577023
	0.00       6.00           ----       1.577057
	0.00       7.00           ----       1.577091
	0.00       8.00           ----       1.577125
	0.00       9.00           ----       1.577159
	0.00      10.00           ----       1.577193
	0.00      11.00           ----       1.577227
	.....
	0.00     300.00           ----       1.586942
	1000.00       1.00       1.570517       1.558042
	1000.00       2.00       1.575324       1.572178
	1000.00       3.00       1.576246       1.574845
	1000.00       4.00       1.576590       1.575802
	1000.00       5.00       1.576768       1.576264
	1000.00       6.00       1.576880       1.576530
	1000.00       7.00       1.576961       1.576704
	1000.00       8.00       1.577026       1.576829
	.....

Note that the CASSCF module also supports the calculation of susceptibility tensors at non-zero user-defined magnetic fields. This is not yet possible with the MRCI module.

7.39.2.6. MCD and Absorption Spectra

The MRCI module of the ORCA program allows calculating MCD spectra and the SOC effects on absorption spectra. The formalism is described in detail by Ganyushin and Neese[283]. The approach is based on the direct calculation of the transition energies and transition probabilities between the magnetic levels. Namely, the differential absorption of LCP- and RCP photons for transitions from a manifold of initial states \(A\) to a manifold of final states \(J\). Using Fermi’s golden rule, the Franck-Condon approximation, assuming a pure electronic dipole mechanism and accounting for the Boltzmann populations of the energy levels, the basic equation of MCD spectroscopy may be written as (atomic units are used throughout):

(7.244)\[\frac{\Delta \varepsilon }{E}=\gamma \sum\limits_{a,j} { \left({N_{a} -N_{j} } \right)\left({ \left|{ \left\langle { \Psi_{a} \left|{m_{\text{LCP} } } \right|\Psi_{j} } \right\rangle} \right|^{2}-\left|{\left\langle { \Psi_{a} \left|{ m_{\text{RCP} } } \right|\Psi_{j} } \right\rangle} \right|^{2} } \right)f\left( E \right)} \]

Here \(a\) and \(j\) label members of the initial and state manifold probed in the experiments.

(7.245)\[N_{a} \left({ B,T} \right)=\frac{\exp \left({ -E_{a} /kT} \right)}{\sum\limits_i { \exp \left({ -E_{i} /kT} \right)} } \]

denotes the Boltzmann population and if the \(a\)-th ground state sublevel at energy \(E_{a}\), \(f\left( E \right)\) stands for a line shape function, and \(\gamma\) denotes a collection of constants. The electric dipole operators are given by:

(7.246)\[m_{\text{LCP} } \equiv m_{x} -im_{y} \]
(7.247)\[m_{\text{RCP} } \equiv m_{x} +im_{y} \]

They represent linear combinations of the dipole moment operator:

(7.248)\[\vec{{m} }=\sum\limits_N { Z_{N} \vec{{R} }_{N} } -\sum\limits_i { \vec{{r} }_{i} } \]

where \(N\) and \(i\) denotes summations of nuclei (at positions \(\vec{{R} }_{N}\) with charges \(Z_{N})\) and electrons (at positions \(\vec{{r} }_{i} )\) respectively. The calculated transition dipole moment are subjected to the space averaging over the Euler angles which is performed by a simple summation over three angular grids.

(7.249)\[\left({ \frac{\Delta \varepsilon }{E} } \right)_{ev} =\frac{1}{8\pi ^{2} }\int\limits_{\psi =0}^{2\pi } { \int\limits_{\phi =0}^{2\pi } {\int\limits_{\theta =0}^\pi{ \left({ \frac{\Delta \varepsilon }{E} } \right)\sin \theta d\theta d\phi d\psi } } } \approx \sum\limits_{\mu \eta \tau } { \left({ \frac{\Delta \varepsilon }{E} } \right)_{\mu \eta \tau } } \sin \theta_{\tau } \]

Finally, every transition is approximated by a Gaussian curve with a definite Gaussian shape width parameter. Hence, the final calculated MCD spectrum arises from the superposition of these curves.

As an illustration, consider calculation of a classical example of MCD spectrum of [Fe(CN) \(_{6}\)]\(^{3-}\). The mrci block of the input file is presented below.

%mrci   ewin -4,10000
        citype mrddci2
        intmode ritrafo
        Tsel  1e-6
        Tpre  1e-5
        etol 1e-8
        rtol 1e-8
        cimode direct2
        maxmemint 300
        solver diis
        davidsonopt 0
        nguessmat 150
        MaxIter 50
        LevelShift 0.5
        PrintLevel 3
        soc
        printlevel 3
        Domcd true        # perform the MCD calculation
        NInitStates 24    # number of SOC and SSC state to account
                          # Starts from the lowest state
        NPointsTheta 10   # number of integration point for
        NPointsPhi 10     # Euler angles
        NPointsPsi 10     #
        B 43500           # experimental magnetic field strength
                          # in Gauss
        Temperature 299.0 # experimental temperature (in K)
        end
        newblock 2 *
         nroots 12
         excitations cisd
         refs cas(23,12) end
         end
        end

The parameters B and Temperature can be assigned in pairs, i.e. B \(=\) 1000, 2000, 3000…, Temperature \(=\) 4, 10, 300…. The program calculates the MCD and absorption spectra for every pair. Now for every point of the integration grid the program prints out the Euler angles, the orientation of the magnetic field in the coordinate system of a molecule, and the energy levels.

Psi = 36.000 Phi = 72.000 Theta = 20.000

Bx = 8745.0 By = 12036.5 Bz = 40876.6

Energy levels (cm-1,eV):Boltzmann populations for T = 299.000   K
   0 :         0.000     0.0000         4.53e-01
   1 :         3.943     0.0005         4.45e-01
   2 :       454.228     0.0563         5.09e-02
   3 :       454.745     0.0564         5.08e-02
   4 :      1592.142     0.1974         2.13e-04
   5 :      1595.272     0.1978         2.10e-04
   6 :     25956.363     3.2182         2.59e-55
   7 :     25958.427     3.2184         2.56e-55
   8 :     25985.656     3.2218         2.25e-55
   9 :     25987.277     3.2220         2.23e-55
  10 :     26070.268     3.2323         1.49e-55
  11 :     26071.484     3.2325         1.49e-55
  12 :     31976.645     3.9646         6.78e-68
  13 :     31979.948     3.9650         6.67e-68
  14 :     32018.008     3.9697         5.56e-68
  15 :     32021.074     3.9701         5.48e-68
  16 :     32153.427     3.9865         2.90e-68
  17 :     32157.233     3.9870         2.84e-68
  18 :     42299.325     5.2444         1.81e-89
  19 :     42303.461     5.2450         1.78e-89
  20 :     42346.521     5.2503         1.45e-89
  21 :     42348.023     5.2505         1.44e-89
  22 :     42456.119     5.2639         8.53e-90
  23 :     42456.642     5.2640         8.51e-90

In the next lines, ORCA calculates the strength of LCP and RCP transitions and prints the transition energies, the difference between LCP and RCP transitions (denoted as C), and sum of LCP and RCP transitions (denoted as D), and C by D ratio.

dE     Na         C           D        C/D

  0 -> 1        3.943 4.53e-01  1.14e-13  8.13e-13   0.00e+00
  0 -> 2      454.228 4.53e-01  5.01e-09  9.90e-09   5.06e-01
  0 -> 3      454.745 4.53e-01 -4.65e-09  7.00e-09  -6.65e-01
  0 -> 4     1592.142 4.53e-01 -8.80e-08  1.02e-07  -8.67e-01
  0 -> 5     1595.272 4.53e-01 -2.29e-08  2.97e-08  -7.71e-01
  0 -> 6    25956.363 4.53e-01  1.22e+01  9.60e+01   1.27e-01
  0 -> 7    25958.427 4.53e-01  3.44e+01  3.52e+01   9.77e-01
  0 -> 8    25985.656 4.53e-01  3.83e+01  1.70e+02   2.25e-01
  0 -> 9    25987.277 4.53e-01 -7.73e+00  6.03e+01  -1.28e-01
  0 ->10    26070.268 4.53e-01 -6.11e+00  2.85e+01  -2.14e-01
  0 ->11    26071.484 4.53e-01  6.17e+00  9.21e+00   6.70e-01
  0 ->12    31976.645 4.53e-01  2.45e+01  6.21e+01   3.95e-01
  0 ->13    31979.948 4.53e-01 -6.58e+01  6.93e+01  -9.50e-01
  0 ->14    32018.008 4.53e-01  3.42e-01  1.07e+02   3.21e-03
  0 ->15    32021.074 4.53e-01 -6.16e+00  3.24e+01  -1.90e-01
  0 ->16    32153.427 4.53e-01 -4.73e+01  1.37e+02  -3.46e-01
  0 ->17    32157.233 4.53e-01 -1.02e+00  5.97e+01  -1.71e-02
  0 ->18    42299.325 4.53e-01  6.47e+00  2.11e+01   3.07e-01
  0 ->19    42303.461 4.53e-01 -2.59e+00  7.61e+00  -3.40e-01
  0 ->20    42346.521 4.53e-01  1.90e+01  8.99e+01   2.11e-01
  0 ->21    42348.023 4.53e-01  3.36e+00  3.55e+00   9.48e-01
  0 ->22    42456.119 4.53e-01  2.52e-01  4.86e-01   5.20e-01
  0 ->23    42456.642 4.53e-01 -2.01e+00  2.91e+00  -6.91e-01
  1 -> 2      450.285 4.45e-01  4.59e-09  6.87e-09   6.69e-01
  1 -> 3      450.802 4.45e-01 -4.96e-09  9.73e-09  -5.09e-01

All C and D values are copied additionally into the text files input.1.mcd, input.2.mcd…, for every pair of Temperature and B parameters. These files contain the energies and C and D values for every calculated transition. These files are used by the program orca_mapspc to calculate the spectra lines. The orca_mapspc program generates from the raw transitions data into spectra lines. The main parameters of the orca_mapspc program are described in section 7.18.1. A typical usage of the orca_mapspc program for MCD spectra calculation for the current example may look as the following:

orca_mapspc input.1.mcd MCD -x020000 -x150000 -w2000

Here the interval for the spectra generation is set from 20000 cm\(^{-1}\) to 50000 cm\(^{-1}\), and the line shape parameter is set to 2000 cm\(^{-1}\).

Very often, it is desirable to assign different line width parameters to different peaks of the spectra to obtain a better fitting to experiment. orca_mapspc can read the line shape parameters from a simple text file named as input.1.mcd.inp. This file should contain the energy intervals (in cm\(^{-1})\) and the line shape parameters for this energy interval in the form of:

20000 35000 1000
35000 40000 2000
40000 50000 1000

This file should not be specified in the executing command; orca_mapspc checks for its presence automatically:

orca_mapspc input.1.mcd MCD -x020000 -x150000
Mode is MCD
Number of peaks            ... 276001
Start wavenumber [cm-1]    ...  20000.0
Stop wavenumber [cm-1]     ...  50000.0
Line width parameters are taken from the file:input.1.mcd.inp
Number of points           ... 1024

Finally, the orca_mapspc program generates the output text file input.1.mcd.dat which contains seven columns of numbers: transition energies, intensities of MCD transitions (the MCD spectrum), intensities of absorption transitions (the absorption spectrum), the ratio between the MCD and absorption intensities, and the last three columns represent the “sticks” of the corresponding transitions.

Energy      C         D         C/D       C        D       E/D
24310.8   0.6673    980.2678   0.0006   0.0000   0.0000   0.0000
24340.1   0.8471   1174.3637   0.0007  -0.0001   0.0129  -0.0112
24369.5   1.0664   1408.5788   0.0007   0.0001   0.0281   0.0033
24398.8   1.3325   1690.5275   0.0007   0.0000   0.0000   0.0000
24428.1   1.6542   2029.0152   0.0008   0.0000   0.0000   0.0000
24457.4   2.0416   2434.1699   0.0008   0.0000   0.0332   0.0003

Now the MCD and the absorption spectra can be plotted with a suitable graphical program, for instance with the Origin program.

../../_images/713.svg

Fig. 7.41 Calculated MCD and absorption spectra of [Fe(CN) \(_6\)]\(^{3-}\) (dash lines) compared to experimental spectra (solid lines).

7.39.2.7. Addition of Magnetic Fields

The inclusion of the Zeeman contribution into the QDPT procedure allows to obtain the splittings of the magnetic levels in an external magnetic field. The switch for this calculation and the magnetic field strength are defined in the soc subblock of the mrci block. Optionally the wave function decomposition can be printed for MagneticField_PrintLevel larger 0. The latter employs the thresh TPrint to omit small contributions from the printing:

%mrci
  soc
    DoSOC true          #
    DoSSC true          #
    MagneticField true  # default false
    B 1,10,100,1000     # Strengh of the magnetic field in Gauss.
                        # 4000 is the default value
    
    # Optional printing of the wave function for each
    # magnetic field settings
    MagneticField_PrintLevel 0 # default (disabled)
    TPrint                    1e-3 
  end
end

Then, the output contains three sets of data of splittings of the magnetic levels with the magnetic field applied parallel to x, y, and z directions:

End B (Gauss)  Energy levels (cm-1) and populations for B || x

     1.0   -0.030   0.333   0.012   0.333   0.018   0.333
    10.0   -0.030   0.333   0.012   0.333   0.018   0.333
   100.0   -0.031   0.333   0.012   0.333   0.020   0.333
  1000.0   -0.102   0.333   0.012   0.333   0.091   0.333

B (Gauss)  Energy levels (cm-1) and populations for B || y

     1.0   -0.030   0.333   0.012   0.333   0.018   0.333
    10.0   -0.030   0.333   0.012   0.333   0.018   0.333
   100.0   -0.032   0.333   0.014   0.333   0.018   0.333
  1000.0   -0.105   0.334   0.018   0.333   0.087   0.333

B (Gauss)  Energy levels (cm-1) and populations for B || z

     1.0   -0.030   0.333   0.012   0.333   0.018   0.333
    10.0   -0.030   0.333   0.011   0.333   0.018   0.333
   100.0   -0.030   0.333   0.005   0.333   0.025   0.333
  1000.0   -0.079   0.333  -0.030   0.333   0.108   0.333

Here the number in a row represents the strength of the magnetic field (in Gauss), and the following pairs of numbers denote the energy of the magnetic level (in cm\(^{-1})\) with its occupation number. This table can be readily plotted with any suitable graphical program.

7.39.2.8. Relativistic Picture Change in Douglas-Kroll-Hess SOC and Zeeman Operators

The DKH correction to the SOC operator is implemented in ORCA as a correction to the one-electron part of the SOMF operator. The DKH transformation is performed up to the second order, and the two-electron part in our implementation is left untransformed. However, the electronic density employed for evaluating the SOMF matrix elements is obtained from a scalar relativistic calculation. The inclusion of the DKH correction is controlled by the picturechange key in the rel block:

%rel method DKH              # relativistic method
     picturechange 2         # include the DKH correction to SOC
     end

The “picturechange” key can be set to 0, 1, and 2 for no picture change, the first order, and the second order DKH transformations of the SOC operator.

With “picturechange” set to 1 or 2 the DKH correction are applied in the first order to the Zeeman operator. This correction has a visible effect on calculated g-tensors for molecules containing third-row and heavier atoms.

7.39.2.9. X-ray Spectroscopy

Likewise to the CASCI/NEVPT2 computational protocol presented in section Core excited states with (C/R)ASCI/NEVPT2 starting from ORCA 4.2 the MRCI module can be used to compute core excited spectra, namely X-ray absorption (XAS) and resonant inelastic scattering (RIXS) spectra.

As discussed in the case of CASCI/NEVPT2 protocol Core excited states with (C/R)ASCI/NEVPT2 a similar strategy is followed to compute XAS/RIXS spectra within the MRCI module. In principle the XAS/RIXS spectra calculations require two steps:

  • In a first step one needs to optimize the valence active space orbitals in the framework of SA-CASSCF calculations, e.g. including valence excited states in the range between 6 to 15 eV.

  • In a second step the relevant core orbitals are rotated into the active space and the MRCI problem is solved by saturating the excitation space with singly core-excited electronic configurations using the previously optimized sets of orbitals.

  • The core orbitals are also included in the XASMOs definition. The use of this keyword is two fold. At first it effecteively reduces the number of the generated configuration state functions (CSFs) to those that exclusively contain contributions from the defined core orbitals. In the case of RIXS also XES (see below) the specified XASMOs are used to define intermediate or core ionized states.

A representative input for the case of Fe(Cl) \(_4\) is provided bellow:

  • In the first step one performs a SA-CASSCF calculation for the 5 and 15 quintet and triplet states (FeIICl4.casscf.inp).

!CC-PWCVTZ-DK cc-pVTZ/C  RIJCOSX SARC/J TightSCF DKH2

%rel
  FiniteNuc true
end

%basis
  newgto Cl "cc-pVTZ-DK" end
  newauxgto Cl "cc-pVTZ/C" end
end

%method FrozenCore FC_NONE
end

%casscf nel   6
  norb   5
  mult   5,3
  nroots 5,15
  switchstep nr
end

* xyz -2 5
Fe  -17.84299991694815     -0.53096694321123      6.09104775508499
Cl  -19.84288422845700      0.31089495619796      7.04101319789001
Cl  -17.84298666758073      0.11868125024595      3.81067954087770
Cl  -17.84301352218429     -2.87052442818457      6.45826391412877
Cl  -15.84311566482982      0.31091516495189      7.04099559201853
*
  • In a second step the core orbitals are rotated in the active space and the MRCI problem is solved by saturating the excitation space with all the quintet and triplet states that involve single excitations from the core orbitals (FeIICl4-mrci.inp)

!MORead CC-PWCVTZ-DK cc-pVTZ/C  RIJCOSX SARC/J TightSCF DKH2 

%moinp "FeIICl4-casscf.gbw"

%rel
  FiniteNuc true
end

%method FrozenCore FC_NONE
end

%scf
  rotate { 6,42,90} { 7,43,90} { 8,44,90}  end
end

%basis
  newgto Cl "cc-pVTZ-DK" end
  newauxgto Cl "cc-pVTZ/C" end
end


%casscf
  nel    12
  norb   8
  mult   5,3
  nroots 34,195
  maxiter 1
  switchstep nr
end


%mrci
  CIType MRCI
  intmode fulltrafo
  XASMOs 42, 43, 44
  newblock 5 *
  nroots 34
  excitations cisd
  refs CAS(12,8)
  end
  end
  newblock 3 *
  nroots 195
  excitations cisd
  refs CAS(12,8)
  end
  end
  maxiter 100
  soc
  printlevel 3
  dosoc true
  end
end


* xyz -2 5
Fe  -17.84299991694815     -0.53096694321123      6.09104775508499
Cl  -19.84288422845700      0.31089495619796      7.04101319789001
Cl  -17.84298666758073      0.11868125024595      3.81067954087770
Cl  -17.84301352218429     -2.87052442818457      6.45826391412877
Cl  -15.84311566482982      0.31091516495189      7.04099559201853
*

In a similar fashion Multireference Equation of Motopn Couple Cluster MR-EOM-CC (see next section) can also be used to compute X-ray spectra. Further information can be found in reference[549]

As it is explicitly described in the respective ROCIS section RIXS spectra can be requested by the following keywords:

RIXS    true    # Request RIXS calculation (NoSOC)   
RIXSSOC true    # Request RIXS calculation (with SOC)   
Elastic true    # Request RIXS calculation (Elastic)

Please consult section Resonant Inelastic Scattering Spectroscopy for processing and analyzing the generated spectra

Likewise to TDDFT (Use of TD-DFT for the Calculation of X-ray Absorption Spectra) ROCIS (General Use) and CASSCF (Core excited states with (C/R)ASCI/NEVPT2) the computed transition densities also in the presence of SOC can be taken beyond the dipole approximation by using the OPS tool for details.

  1. by performing a multiple expantion up to second order

  2. by computing the exact transition moments

The whole set of spectroscopy tables can be requested with the following commands:

%mrci
 DoDipoleLength        true
 DoDipoleVelocity      true
 DoHigherMoments       true
 DecomposeFoscLength   true
 DecomposeFoscVelocity true
 DoFullSemiclassical   true
end

More details can be found in TDDFT (Use of TD-DFT for the Calculation of X-ray Absorption Spectra) ROCIS (General Use) and CASSCF (Core excited states with (C/R)ASCI/NEVPT2) sections.

Starting from ORCA 4.2 the previously reported RASCI-XES protocol reference[690], which can compute K\(_\beta\) Mainline XES spectra, can be processed entirely within the ORCA modules. In ORCA 5.0 a similar protocol (CASCI-XES) exist in the CASSCF module (Core excited states with (C/R)ASCI/NEVPT2)

  • Like above or in the CASCI/NEVPT2 case in a first step one needs to optimize the valence active space orbitals in the framework of SA-CASSCF calculations, e.g. including valence excited states in the range between 6 to 15 eV for the N electron system.

  • In a second step the metal 1s and 3p orbitals are rotated in the active space and the 1s MO is defined in the XASMOs list

  • Computes the XES spectrum in the RASCI framework for the N-1 electron system in the presence of SOC if the XESSOC keyword for all the states that are dominated by 3p-1s electron decays.

A representative input sequence for the case of Fe(Cl) \(_6\) is provided bellow:

As described in reference[690] at first for a CAS(5,5) the excitation space is saturated by the sextet as well as the 24 quartet and the 75 doublet states which are optimized in the SA-CASSCF fashion.

!ZORA  def2-TZVP  def2-TZVP/C

%cpcm
  epsilon 80
  refrac 1.33
  surfacetype gepol_ses
end

%scf
  MaxDisk 40000
end

%casscf
  nel 5
  norb 5
  mult 6,4,2
  nroots 1,24,75
  shiftup 0.5
  shiftdn 0.5
  trafostep RI
  maxiter 150
end

*xyzfile -3  6 
Fe      0.0000       0.0000     0.000000
Cl      2.478        0.0000     0.000
Cl     -2.478        0.0000     0.000
Cl      0.000005     2.478      0.00000
Cl      0.000005    -2.478     -0.0000
Cl     -0.000       -0.000      2.478
Cl      0.000       -0.0000    -2.478
*

In following the 1s and 3p Fe based MOs are rotated in the active space and the XES spectra are computed for the [Fe(Cl) \(_6\)]\(^+\) system for the 4 septet and 81 quintet states.

! ZORA  def2-TZVP  def2-TZVP/C noiter moread AllowRHF

%moinp "fecl6_casscf.gbw"

%cpcm
  epsilon 80
  refrac 1.33
  surfacetype gepol_ses
end

%scf
  MaxDisk 40000
end

%scf
  rotate {0,59,90} {36, 60, 90} {37,61,90} {38,62,90}  end
end

%mrci citype mrci
  UseIVOs false
  Etol 1e-5
  newblock 5 *
  excitations none
  nroots 81
  refs ras(12:4 1/5/ 0 0) end
  end
  newblock 7 *
  excitations none
  nroots 4
  refs ras(12:4 1/5/ 0 0) end
  end
  XASMOs 59
  soc
  dosoc true
  XESSOC true
  end
end

*xyzfile -2  7 
Fe      0.0000       0.0000     0.000000
Cl      2.478        0.0000     0.000
Cl     -2.478        0.0000     0.000
Cl      0.000005     2.478      0.00000
Cl      0.000005    -2.478     -0.0000
Cl     -0.000       -0.000     2.478
Cl      0.000       -0.0000    -2.478
*

As a result the X-ray emission spectrum is calculated and the intensities are computed on the basis of the transition electric dipole moments

Printing the XES spectrum ...
	
	-------------------------------------------------------------------------------------
	SPIN-ORBIT X-RAY EMISSION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS
	-------------------------------------------------------------------------------------
	Transition          Energy           INT             TX        TY        TZ 
	1  421 ->    5      7228.632     0.000000000000     0.00000   0.00000   0.00000
	2  422 ->    5      7228.632     0.000000000000     0.00000   0.00000   0.00000
	3  423 ->    5      7228.632     0.000000000000     0.00000   0.00000   0.00000
	4  424 ->    5      7228.632     0.000000000000     0.00000   0.00000   0.00000
	5  425 ->    5      7228.632     0.000000000000     0.00000   0.00000   0.00000
	...
	242  422 ->   25      7177.286     0.000917305388     0.00025   0.00171   0.00149
	243  423 ->   25      7177.286     0.002043577370     0.00197   0.00211   0.00181
	244  424 ->   25      7177.286     0.000789769987     0.00114   0.00133   0.00119
	245  425 ->   25      7177.286     0.000026130790     0.00018   0.00034   0.00002
	246  426 ->   25      7177.286     0.000035191741     0.00034   0.00028   0.00003
	247  427 ->   25      7177.286     0.005143175830     0.00294   0.00345   0.00296
	248  428 ->   25      7177.341     0.000000000000     0.00000   0.00000   0.00000
	249  429 ->   25      7177.341     0.000000000001     0.00000   0.00000   0.00000
	250  430 ->   25      7177.341     0.000000000001     0.00000   0.00000   0.00000
	251  431 ->   25      7177.341     0.000000000000     0.00000   0.00000   0.00000
	252  432 ->   25      7177.341     0.000000000000     0.00000   0.00000   0.00000
	...
	4991  431 ->  420      7153.111     0.000195885011     0.00106   0.00000   0.00000
	4992  432 ->  420      7153.111     0.002719228427     0.00256   0.00299   0.00002
	
	All Done
	-------------------------------------------------------------------------------------

The resulted XES spectrum can be visualized by processing the above output file with the orca_mapspc

orca_mapspc fecl6_xes.out XESSOC -x07140 -x17190 -w4.0 -eV -n10000

This will result in Fig. 7.42.

../../_images/RASCI_XESSOC.svg

Fig. 7.42 Calculated RASCI K\(_\beta\) XES spectrum of [Fe(Cl) \(_6\)]\(^{+}\) .