5.24. Electron Paramagnetic Resonance (EPR) Parameters

Electron paramagnetic resonance probes the magnetic sublevels of the electronic ground state. In the %eprnmr input block, calculations of the parameters of a spin Hamiltonian of the form

(5.102)HS=βeBgS+NSANIN+SDS

can be requested. The properties are evaluated via analytical derivative/response theory methods. For the analytic response approach based on CC, please see AUTOCI Response Properties via Analytic Derivatives. For evaluation of spin Hamiltonian parameters using quasi-degenerate perturbation theory, you may be interested in Magnetic Properties Through Quasi Degenerate Perturbation Theory A complete list of possible keywords for the eprnmr module can be found in EPRNMR - keywords for magnetic properties. The cartesian index convention for EPR can be found in Cartesian Index Conventions for EPR and NMR Tensors.

5.24.1. The g-Tensor

The g-tensor determines the position of the signal in an EPR spectrum. It can be calculated as a second derivative of the energy and it is implemented as such in ORCA for the SCF methods (HF and DFT), CASSCF, as well as all-electron MP2 (or RI-MP2) and double-hybrid DFT, plus Coupled Cluster (CCSD and CCSD(T)).

5.24.1.1. Theory

The g-tensor has four contributions which may be evaluated in ORCA,

(5.103)g=ge1+gRMC+gDSO+gPSO.

These contributions are:

  • The spin Zeeman contribution, which is isotropic and characterized by a simple constant, the free-electron g-value ge=2.002319

    (5.104)gμνSZ=δμνge

    meaning that it doesn’t require any computation. The remaining terms contribute to the deviation from the free-electron value, which is sometimes called the g-shift.

  • The first-order relativistic mass correction, which is a usually rather small scalar-relativistic correction and can be evaluated as an expectation value over the spin density

    (5.105)gRMC=α2ge2Sk,lPklαβϕk|T^|ϕl.
  • The diamagnetic spin-orbit (also called gauge correction) term, another first-order correction which is often small

    (5.106)gμνDSO=α2ge4Sk,lPklαβϕk|Aξ(rA)[rArOrA,μrO,ν]|ϕl.
  • The second-order orbital Zeeman/SOC contribution, usually the main source of deviation from the free-electron g-value in a molecule,

    (5.107)gμνPSO=ge2Sk,lPklαβBμϕk|hνSOC|ϕl.

    Its calculation requires the solution of response equations to obtain the perturbed spin density, this is done for the magnetic field perturbation in ORCA. Precise details depend on the level of theory. At the SCF level, it is achieved by solving coupled-perturbed SCF equations (see CP-SCF Options). At the CASSCF level, the CP-CASSCF equations are solved (CASSCF Linear Response).

In these equations, S is the total spin, α the fine structure constant, Pαβ is the spin density matrix, ϕ is the basis set, T^ is the kinetic energy operator, ξ(rA) an approximate radial operator, hSOC the spatial part of an effective one-electron spin-orbit operator and Bμ is a component of the magnetic field.

5.24.1.2. Basic usage

As an example, consider the following simple g-tensor job:

! BP86 Def2-SVP TightSCF g-tensor SOMF(1X)
* int 1 2
   O  0  0  0  0      0  0
   H  1  0  0  1.1056 0  0
   H  1  2  0  1.1056 109.62 0
*

The simplest way is to call the g-tensor property in the simple input line as shown above. It can also be specified in the %eprnmr block with gtensor true. SOMF(1X) defines the chosen spin-orbit coupling (SOC) operator.

The output looks like the following. It contains information on the contributions to the g-tensor (relativistic mass correction, diamagnetic spin-orbit term (= gauge-correction), paramagnetic spin-orbit term (= OZ/SOC)), the isotropic g-value and the orientation of the total tensor.

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

The g-matrix:
              2.0104321   -0.0031354   -0.0000000
             -0.0031354    2.0081968   -0.0000000
             -0.0000000   -0.0000000    2.0021275

Breakdown of the contributions
 gel          2.0023193    2.0023193    2.0023193
 gRMC        -0.0003174   -0.0003174   -0.0003174
 gDSO(tot)    0.0000808    0.0001539    0.0001515
 gPSO(tot)    0.0000449    0.0038301    0.0104898
             ----------   ----------   ----------
 g(tot)       2.0021275    2.0059858    2.0126431 iso=  2.0069188
 Delta-g     -0.0001917    0.0036665    0.0103238 iso=  0.0045995

 Orientation:
  X           0.0000000    0.5762906   -0.8172448
  Y           0.0000000    0.8172448    0.5762906
  Z           1.0000000   -0.0000000    0.0000000

g-tensor calculations at the SCF level are not highly demanding in terms of basis set size. Basis sets that give reliable SCF results (at least valence double-zeta plus polarization) usually also give reliable g-tensor results. For many molecules the Hartree-Fock approximation will give reasonable predictions. In a number of cases, however, it breaks down completely. DFT is more robust and the number of molecules where it fails is much smaller. Among the density functionals, the hybrid functionals seem to be the most accurate.

Note

There are different options for treatment of the spin-orbit coupling operator available in the program. The defaults should be quite reliable. Check out The Spin-Orbit Coupling Operator for more information.

5.24.1.3. Gauge origin treatment

The g-tensor is not gauge-independent. Without taking special precautions, the results will have an unphysical dependence on where the gauge origin is placed with respect to the molecule. Unless a fully invariant procedure (such as GIAOs) is used, this undesirable aspect is always present in the calculations.

Starting from ORCA 5.0, the default treatment for g-tensor calculations gauge is the GIAO approach (GIAO stands for “gauge including atomic orbitals”). GIAOs are available at the SCF level, RI-MP2 as well as Coupled Cluster up to CCSD(T) level. GIAOs are not currently available with CASSCF linear response and a common gauge origin must be provided in the %eprnmr block.

The GIAO one-electron integrals are done analytically by default whereas the treatment of the GIAO two-electron integrals is chosen to be same as for the SCF. The available options which can be set with giao_1el / giao_2el in the %eprnmr block can be found in section EPRNMR - keywords for magnetic properties. Concerning the computational time, for small systems, e.g. phenyl radical (41 electrons), the rijk-approximation is good to use for the SCF-procedures as well as the GIAO two-electron integrals. Going to larger systems, e.g. chlorophyll radical (473 electrons), the rijcosx-approximation reduces the computational time enormously. While the new default grid settings in ORCA 5.0 (defgrid2) should be sufficient in most cases, certain cases might need the use of defgrid3. Note that for the current implementation of CC-NMR/EPR, both giao_1el and giao_2el need to be evaluated fully analytically, which is also the default when CC-NMR/EPR is requested.

If the choice of the gauge origin is not outrageously poor, a common gauge origin often gives reasonable results for g-tensors (much more reasonable than in NMR, where it shouldn’t be used). This is especially true with a large basis set. The origin may be modified with the keyword ori inside the %eprnmr input block. If you are using a common gauge origin, it is wise to check the sensitivity of the results with respect to origin location, especially when small g-shifts on the order of only a few hundred ppm are calculated.

5.24.1.4. Keywords

Table 5.13 %eprnmr block input keywords relevant for g-tensor calculations. For more GIAO-related options, see EPRNMR - keywords for magnetic properties

Keyword

Options

Description

gtensor

true/false

Calculate the g-tensor

gtensor_1el2el

true/false

Print the 1- and 2-electron contributions to the g-tensor

ori

GIAO

use the GIAO formalism (default)

CenterOfElCharge

Common gauge origin at center of electronic charge

CenterOfNucCharge

Common gauge origin at center of nuclear charge

CenterOfSpinDens

Common gauge origin at center of spin density

CenterOfMass

Common gauge origin at center of mass

N

Common gauge origin at atom N

x, y, z

Common gauge origin at position x, y, z (in coordinate input units, default Angstrom)

do_giao_soc2el

true/false

whether to include the 2-el contribution to GIAO term from the SOMF operator (usually small but expensive, disabled by default)

5.24.2. Hyperfine Coupling

Hyperfine couplings characterize the interaction between the electronic spin and the spin of a nucleus in the molecule. In typical EPR spectra, they lead to a splitting of the signal.

5.24.2.1. Theory

The hyperfine coupling has four contributions which may be evaluated in ORCA,

(5.108)AN=aiso1+Adip+Aorb+AGC.

These contributions are:

  • The isotropic Fermi contact term that arises from finite spin density on the nucleus under investigation. It is calculated for nucleus N from:

    (5.109)aiso(N)=(43πSz1)gegNβeβNρ(RN)

    Here, Sz is the expectation value of the z-component of the total spin, ge and gN are the electron and nuclear g-factors and βe and βN are the electron and nuclear magnetons respectively. ρ(RN) is the spin density at the nucleus. The proportionality factor PN=gegNβeβN is commonly used and has the dimensions MHz bohr3 in ORCA.

  • The spin dipole term that arises from the through-space dipole-dipole interaction of the magnetic nucleus with the magnetic moment of the electron. It is also calculated as an expectation value over the spin density as:

    (5.110)Aμνdip(N)=PNklρklϕk|rN5(3rNμrNνδμνrN2)|ϕl

    where ρ is the spin-density matrix and rN is a vector of magnitude rN that points from the nucleus in question to the electron ({ϕ} is the set of basis functions).

  • The second-order spin-orbit contribution, which arises as a cross-term between the spin-orbit and nucleus-orbit coupling operators. It requires the solution of coupled-perturbed SCF equations and is consequently more computationally demanding. The contribution can be written as:

    (5.111)Aμνorb(N)=12SPNklρklIνϕk|hμSOC|ϕl

    The derivative of the spin density is computed by solving the coupled-perturbed SCF equations with respect to the nucleus-orbit coupling perturbation, which is represented by the operator

    (5.112)hνNOC(N)=iriA3li,ν(N)

    where the sum is over electrons and N is the nucleus in question.

  • The gauge-correction contribution (sometimes also called diamagnetic). This term is often small. However, it is needed in order to get exactly gauge-invariant results. We recently showed that the gauge correction can become crucial in the long-distance limit between the nuclear spin and the electron spin. This is relevant for pseudocontact NMR chemical shifts (PCS).[702]

Note

Second-order HFCs can be quite significant for heavier nuclei and are certainly good to include for transition metal complexes. Available treatments of the spin-orbit coupling operator are described under The Spin-Orbit Coupling Operator.

5.24.2.2. Basic usage

Hyperfine and quadrupole couplings can be requested in the %eprnmr input block. Since there may be several nuclei that you are interested in, the input is relatively sophisticated.

An example how to calculate the hyperfine and field gradient tensors for the CN radical is given below:

! PBE0 def2-MSVP TightSCF
* int 0 2
   C  0  0  0  0     0  0
   N  1  0  0  1.170 0  0
*
%eprnmr  Nuclei = all C { aiso, adip }
         Nuclei = 2 { aiso, adip, fgrad }
         end

In this example, the hyperfine tensor is calculated for all carbon atoms and atom 2, which is nitrogen in this case. (Additionally, the electric field gradient is calculated for nitrogen.)

Note

  • Counting of atom numbers starts from 1.

  • Beware - all nuclei mentioned in one line will be assigned the same isotopic mass! If several different nuclei are desired (such as C and H), there has to be a new line for each of them.

  • You have to specify the Nuclei statement after the definition of the atomic coordinates or the program will not figure out what is meant by “all”.

The output looks like the following. It contains detailed information about the individual contributions to the hyperfine couplings, its orientation, its eigenvalues, the isotropic part and (if requested) also the quadrupole coupling tensor.

-----------------------------------------
ELECTRIC AND MAGNETIC HYPERFINE STRUCTURE
-----------------------------------------

 -----------------------------------------------------------
 Nucleus   0C : A:ISTP=   13 I=  0.5 P=134.1903 MHz/au**3
                Q:ISTP=   13 I=  0.5 Q=  0.0000 barn
 -----------------------------------------------------------
 Raw HFC matrix (all values in MHz):
 -----------------------------------
               695.8952               0.0000              -0.0000
                 0.0000             543.0617              -0.0000
                -0.0000              -0.0000             543.0617

 A(FC)         594.0062             594.0062             594.0062
 A(SD)         -50.9445             -50.9445             101.8890
             ----------           ----------           ----------
 A(Tot)        543.0617             543.0617             695.8952    A(iso)=  594.0062
 Orientation:
  X           0.0000000            0.0000000           -1.0000000
  Y          -0.8111216           -0.5848776           -0.0000000
  Z          -0.5848776            0.8111216            0.0000000

Notes:  (1) The A matrix conforms to the "SAI" spin Hamiltonian convention.
        (2) Tensor is right-handed.

 -----------------------------------------------------------
 Nucleus   1N : A:ISTP=   14 I=  1.0 P= 38.5677 MHz/au**3
                Q:ISTP=   14 I=  1.0 Q=  0.0204 barn
 -----------------------------------------------------------
 Raw HFC matrix (all values in MHz):
 -----------------------------------
                13.2095              -0.0000               0.0000
                -0.0000             -45.6036              -0.0000
                 0.0000              -0.0000             -45.6036

 A(FC)         -25.9993             -25.9993             -25.9993
 A(SD)          39.2088             -19.6044             -19.6044
             ----------           ----------           ----------
 A(Tot)         13.2095             -45.6036             -45.6036    A(iso)=  -25.9993
 Orientation:
  X           1.0000000            0.0000000           -0.0000000
  Y          -0.0000000            0.9996462            0.0265986
  Z           0.0000000           -0.0265986            0.9996462

Notes:  (1) The A matrix conforms to the "SAI" spin Hamiltonian convention.
        (2) Tensor is right-handed.

 Raw EFG matrix (all values in a.u.**-3):
 -----------------------------------
            -0.1832      -0.0000       0.0000
            -0.0000       0.0916       0.0000
             0.0000       0.0000       0.0916

 V(El)       0.6468       0.6468      -1.2935
 V(Nuc)     -0.5551      -0.5551       1.1103
         ----------   ----------   ----------
 V(Tot)      0.0916       0.0916      -0.1832
 Orientation:
  X      -0.0000003    0.0000002    1.0000000
  Y       0.9878165    0.1556229    0.0000003
  Z      -0.1556229    0.9878165   -0.0000002

Note: Tensor is right-handed


 Quadrupole tensor eigenvalues (in MHz;Q= 0.0204 I=  1.0)
  e**2qQ            =    -0.880 MHz
  e**2qQ/(4I*(2I-1))=    -0.220 MHz
  eta               =     0.000
  NOTE: the diagonal representation of the SH term I*Q*I = e**2qQ/(4I(2I-1))*[-(1-eta),-(1+eta),2]

If also EPR g-tensor or D-tensor calculations (see next section) are carried out in the same job, ORCA automatically prints the orientation between the hyperfine/quadrupole couplings and the molecular g- or D-tensor. For more information on this see section orca_euler.

Note

For heavy nuclei you may want to consider the possibility of relativistic effects. Scalar relativistic effects can be handled with several quasi-relativistic Hamiltonians in ORCA. An overview of the possibilities and some recommendations can be found in Relativistic Calculations. Note that relativistic calculations may have special requirements on basis sets. In relativistic property calculations, you should be aware of the importance of picture change corrections (see Picture-Change Effects). In quasi-relativistic calculations with DFT, one should also be very cautious about accuracy of the numerical integration, especially for heavier (transition metal) nuclei.

Note

For the calculation of HFCCs using DLPNO-CCSD, it is recommended to use the tailored truncation settings !DLPNO-HFC1 or !DLPNO-HFC2 in the simple keyword line which correspond to the “Default1” and “Default2” setting in Ref. [703].

If you wish to extract the A tensor for an oligonuclear transition metal complex, the A(iso) value in the output can be processed according to the method described in ref. [704].

5.24.2.3. A note on basis sets

For hyperfine (and quadrupole) couplings, standard basis sets designed for energies and geometry optimizations and are often not satisfactory (especially for atoms heavier than Ne). You should look into a basis set capable of describing the electronic structure near the nucleus. One option is to use a basis set tailored towards “core-property” calculations. The following are good options:

  • EPR-II basis of Barone and co-workers: It is only available for a few light atoms (H, B, C, N, O, F) and is essentially of double-zeta plus polarization quality with added flexibility in the core region, which should give reasonable results.

  • IGLO-II and IGLO-III bases of Kutzelnigg and co-workers: They are fairly accurate but also only available for some first and second row elements.

  • CP basis: They are accurate for first row transition metals as well.

  • uncontracted Partridge basis: They are general purpose HF-limit basis sets and will probably be too expensive for routine use, but are very useful for calibration purposes.

  • the pcH basis sets by Jensen tailored towards hyperfine coupling calculations, also with limited availability (see Jensen Basis Sets).

If ORCA does not yet have a dedicated basis set for your element, you will likely have to tailor the basis set to your needs. You can do this by manually adding s-type primitives with large exponents, or by decontracting core parts of the basis set (see Basis Sets for description of the decontraction feature). You can start by examining the basis set you already have - if you add the statement Print[p_basis] 2 in the %output block (or PrintBasis in the simple input line) the program will print the basis set in input format (for the basis block). You can then add or remove primitives, uncontract core functions, etc. For example, here is a printout of the carbon basis DZP in input format:

# Basis set for element : C
  NewGTO 6
  s 5
   1  3623.8613000000      0.0022633312
   2   544.0462100000      0.0173452633
   3   123.7433800000      0.0860412011
   4    34.7632090000      0.3022227208
   5    10.9333330000      0.6898436475
  s 1
   1     3.5744765000      1.0000000000
  s 1
    1    0.5748324500      1.0000000000
  s 1
    1    0.1730364000      1.0000000000
  p 3
    1    9.4432819000      0.0570590790
    2    2.0017986000      0.3134587330
    3    0.5462971800      0.7599881644
  p 1
    1    0.1520268400      1.0000000000
  d 1
    1    0.8000000000      1.0000000000
  end;

Reading this information is relatively straightforward. For example: the “s 5” stands for the angular momentum and the number of primitives in the first basis function. The following five lines then contain the number of the primitive, the exponent and the contraction coefficient (unnormalized) for each primitive in the function.

Remember that if you add very steep functions, you must increase the size of the integration grid if your calculation uses DFT! Otherwise, your results will be inaccurate. To some extent, the program now automatically adapts the grids when it detects steep functions. If you further need to increase integration accuracy, you could globally increase the radial grid size by increasing IntAcc in the Method block, or increase integration accuracy for individual atoms only, with the latter option being less expensive. More detail can be found in section Other Details and Options. In the present example the changes caused by larger basis sets in the core region and more accurate integration are relatively modest – on the order of 3%. This is, however, still significant if you are a little puristic.

5.24.2.4. Hyperfine Coupling Keywords

Warning

All nuclei specified on one line will be assigned with the same isotopic mass (and other parameters)!

Table 5.14 %eprnmr block input keywords relevant for HFC calculations

Keyword

Selection

What to do

Description

nuclei

all <type>

Selects all nuclei of element <type>

i1,i2,i3...

Selects atoms i1, i2, i3… (starts at 1!)

{ aiso  }

Calculates the isotropic Fermi-contact contribution

{ adip  }

Calculates the spin-dipole contribution

{ aorb  }

Calculates the spin-orbit contribution

{ ist=<N> }

Selects isotope <N> for selection

{ PPP=<N> }

Sets HFC proportionality factor gegNβEβN to <N> for selection

{ III=<N> }

Sets nuclear spin to <N> for selection

{ aiso, adip, aorb }

example requesting multiple actions for selection

Table 5.15 %eprnmr block: possible tweaks for the diamagnetic gauge correction. Unless you have a special reason, it is unlikely that you need to use these settings.

Keyword

Options

Description

hfcgaugecorrection_zeff

true/false

use effective nuclear charges for the gauge correction

hfcgaugecorrection_numeric

true/false

calculate DSO integrals numerically (faster)

hfcgaugecorrection_angulargrid

-1

< 0 means to use the DFT grid setting

hfcgaugecorrection_intacc

-1

< 0 means to use the DFT grid setting

hfcgaugecorrection_prunegrid

-1

< 0 means to use the DFT grid setting

hfcgaugecorrection_bfcutoff

-1

< 0 means to use the DFT grid setting

hfcgaugecorrection_wcutoff

-1

< 0 means to use the DFT grid setting

5.24.3. Zero-Field Splitting

The zero-field splitting (ZFS) is typically the leading term in the Spin-Hamiltonian (SH) for transition metal complexes with a total ground state spin S>1/2 (for reviews and references see chapter Publications Related to ORCA). Its effect is to introduce a splitting of the 2S+1 MS sublevels, which are exactly degenerate at the level of the Born-Oppenheimer Hamiltonian, even in the absence of an external magnetic field.

5.24.3.1. Theory

There are two important contributions to zero-field splitting [705]:

  • A first order term arising from the direct spin-spin interaction

    (5.113)DKLSS=12α2S(2S1)0SS|ijirij2δKL3(rij)K(rij)Lrij5{2s^zis^zjs^xis^xjs^yis^yj}|0SS

    where K,L=x,y,z, α is the fine structure constant ( 1/137 in atomic units), rij is the electronic distance vector with magnitude rij and s^i is the spin-vector operator for the i’th electron. |0SS is the exact ground state eigenfunction of the Born-Oppenheimer Hamiltonian with total spin S and projection quantum number MS=S. Since the spin-spin interaction is of first order, its evaluation presents no particular difficulties.

  • A second-order spin-orbit contribution, which is substantially more complicated. Under the assumption that the spin-orbit coupling (SOC) operator can to a good approximation be represented by an effective one-electron operator (H^SOC=ih^iSOCs^i), ref [656] has derived the following sum-over-states (SOS) equations for the SOC contribution to the ZFS tensor:

    (5.114)DKLSOC(0)=1S2b(Sb=S)Δb10SS|ih^iK;SOCs^i,0|bSSbSS|ih^iL;SOCs^i,0|0SS
    (5.115)DKLSOC(1)=1S(2S1)b(Sb=S1)Δb10SS|ih^iK;SOCs^i,+1|bS1S1bS1S1|ih^iL;SOCs^i,1|0SS
    (5.116)DKLSOC(+1)=1(S+1)(2S+1)b(Sb=S+1)Δb10SS|ih^iK;SOCs^i,1|bS+1S+1bS+1S+1|ih^iL;SOCs^i,+1|0SS

    Here the one-electron spin-operator for electron i has been written in terms of spherical vector operator components si,m with m=0,±1 and Δb=EbE0 is the excitation energy to the excited state multiplet |bSS (all MS components are degenerate at the level of the BO Hamiltonian).

Note

The second-order spin-orbit term may be evaluated using the above sum-over-states formulation in the context of wavefunction theory methods. The following section focuses on a formulation based on SCF analytical derivative/response theory. You may also be interested in QDPT, which is more reliable particularly in systems where zero-field splitting is large.

Note

Available treatments of the spin-orbit coupling operator are described under The Spin-Orbit Coupling Operator.

5.24.3.2. Basic usage

For example, consider the following job on a hypothetical Mn(III)-complex.

! BP86 def2-SVP SOMF(1X)

%eprnmr DTensor  ssandso  
        DSOC     cp  # qro, pk, cvw
        DSS      uno # direct
        end

* int  1 5
Mn 0 0 0  0 0 0
 O 1 0 0  2.05   0   0
 O 1 2 0  2.05  90   0
 O 1 2 3  2.05  90 180
 O 1 2 3  2.05 180   0
 F 1 2 3  1.90  90  90
 F 1 2 3  1.90  90 270
 H 2 1 6  1.00 127   0
 H 2 1 6  1.00 127 180
 H 3 1 6  1.00 127   0
 H 3 1 6  1.00 127 180
 H 4 1 6  1.00 127   0
 H 4 1 6  1.00 127 180
 H 5 1 6  1.00 127   0
 H 5 1 6  1.00 127 180
*

The output documents the individual contributions to the D-tensor, which also contains (unlike the g-tensor) contributions from spin-flip terms.

Note

There are four different variants of the SOC-contribution, which alone should demonstrate that this is a difficult property. The options are:

  • The QRO method is fully documented[360] and is based on a theory developed earlier.[656] The QRO method is reasonable but somewhat simplistic and is superseded by the CP method described below.

  • The Pederson-Khanna model was brought forward in 1999 from qualitative reasoning.[706] It also contains incorrect prefactors for the spin-flip terms. We have nevertheless implemented the method for comparison. In the original form it is only valid for local functionals. In ORCA it is extended to hybrid functionals and HF.

  • The default coupled-perturbed method is a generalization of the DFT method for ZFSs; it uses revised prefactors for the spin-flip terms and solves a set of coupled-perturbed equations for the SOC perturbation. Therefore it is valid for hybrid functionals. It has been described in detail.[707]

The present implementation in ORCA is valid for HF, DFT and hybrid DFT.

The DSS part is an expectation value that involves the spin density of the system. In detailed calibration work[708] it was found that the spin-unrestricted DFT methods behave somewhat erratically and that much more accurate values were obtained from open-shell spin-restricted DFT. Therefore the “UNO” option allows the calculation of the SS term with a restricted spin density obtained from the singly occupied unrestricted natural orbitals.

The DSS part contains an erratic self-interaction term for UKS/UHF wavefunction and canonical orbitals. Thus, UNO is recommended for these types of calculations.[709] If the option DIRECT is used nevertheless, ORCA will print a warning in the respective part of the output.

Note

In case that D-tensor is calculated using the correlated wave function methods such as (DLPNO-/LPNO-)CCSD, one should not use DSS=UNO option.

5.24.3.3. Detailed theory of the coupled-perturbed formulation

In 2007, we have developed a procedure that makes the ZFS calculation compatible with the language of analytic derivatives.[707] Perhaps the most transparent route is to start from the exact solutions of the Born-Oppenheimer Hamiltonian. To this end, we look at the second derivative of the ground state energy (E=0SS|H^|0SS) with respect to a spin-dependent one-electron operator of the general form:

(5.117)h^K;(m)=xK(m)pqhpqKS^pq(m)

Where hpqK is the matrix of the K’th component of the spatial part of the operator (assumed to be imaginary Hermitian as is the case for the spatial components of the SOC operator) and S^pq(m) is the second quantized form of the spin vector operator (m=0,±1). The quantity xK(m) is a formal perturbation parameter. Using the exact eigenfunctions of the BO operator, the first derivative is:

(5.118)ExK(m)|xK(m)=0=pqhpqKPpq(m)

With the components of the spin density:

(5.119)Ppq(m)=0SS|S^pq(m)|0SS

The second derivative with respect to a second component for m=m is:

(5.120)2ExK(m)xL(m)|xK(m)=xL(m)=0=pqhpqKPpq(m)xL(m)

The derivative of the spin density may be written as:

(5.121)Ppq(m)xL(m)=0LSS(m)|S^pq(m)|0SS+0SS|S^pq(m)|0LSS(m)

Expanding the perturbed wavefunction in terms of the unperturbed states gives to first order:

(5.122)|0LSS(m)=n0Δn1|nn|h^L;(m)|0SS

Where |n is any of the |bSM. Thus, one gets:

(5.123)2ExK(m)xL(m)=pqhpqKPpq(m)xL(m)
(5.124)=n0Δn1[0SS|h^L;(m)|nn|h^K;(m)|0SS+0SS|h^K;(m)|nn|h^L;(m)|0SS]

The equality holds for exact states. For approximate electronic structure treatments, the analytic derivative approach is more attractive since an infinite sum over states can never be performed in practice and the calculation of analytic derivative is computationally less demanding than the calculation of excited many electron states.

Using eq. (5.123), the components of the SOC-contribution to the D-tensor are reformulated as

(5.125)DKLSOC(0)=12S2pqhpqK;SOCPpq(0)xL(0)
(5.126)DKLSOC(1)=1S(2S1)pqhpqK;SOCPpq(+1)xL(1)
(5.127)DKLSOC(+1)=1(S+1)(2S+1)pqhpqK;SOCPpq(1)xL(+1)

These are general equations that can be applied together with any non-relativistic or scalar relativistic electronic structure method that can be cast in second quantized form. Below, the formalism is applied to the case of a self-consistent field (HF, DFT) reference state.

For DFT or HF ground states, the equations are further developed as follows:

The SCF energy is:

(5.128)ESCF=VNN+Ph++12ρ(r1)ρ(r2)|r1r2|dr1dr212aXμνκτσPμκσPντσ(μν|κτ)+cDFEXC[ρα,ρβ]

Here VNN is the nuclear repulsion energy and hμν is a matrix element of the one-electron operator which contains the kinetic energy and electron-nuclear attraction terms (ab denotes the trace of the matrix product ab). As usual, the molecular spin-orbitals ψpσ are expanded in atom centered basis functions (σ=α,β):

(5.129)ψpσ(r)=μcμpσϕμ(r)

with MO coefficients cμpσ. The two-electron integrals are defined as:

(5.130)(μν|κτ)=ϕμ(r1)ϕν(r1)r121ϕκ(r2)ϕτ(r2)dr1dr2

The mixing parameter aX controls the fraction of Hartree-Fock exchange and is of a semi-empirical nature. EXC[ρα,ρβ] represent the exchange-correlation energy. The parameter cDF is an overall scaling factor that allows one to proceed from Hartree-Fock theory (aX=1,cDF=0) to pure DFT (aX=0,cDF=1) to hybrid DFT (0<aX<1,cDF=1). The orbitals satisfy the spin-unrestricted SCF equations:

(5.131)Fμνσ=hμν+κτPκτ(μν|κτ)aXPκτσ(μκ|ντ)+cDF(μ|VXCα|ν)

With VXCσ=δEXCδρσ(r) and Pμν=Pμνα+Pμνβ being the total electron density. For the SOC perturbation it is customary to regard the basis set as perturbation independent. In a spin-unrestricted treatment, the first derivative is:

(5.132)ESCFxK(m)=iα(iα|hKsm|iα)+iβ(iβ|hKsm|iβ)=0

For the second derivative, the perturbed orbitals are required. However, in the presence of a spin-dependent perturbation they can no longer be taken as pure spin-up or spin-down orbitals. With respect to the L’th component of the perturbation for spin-component m, the orbitals are expanded as:

(5.133)ψiα;(m)L(r)=aαUaαiα(m);Lψaα(r)+aβUaβiα(m);Lψaβ(r)
(5.134)ψiβ;(m)L(r)=aαUaαiβ(m);Lψaα(r)+aβUaβiβ(m);Lψaβ(r)

Since the matrix elements of the spin-vector operator components are purely real and the spatial part of the SOC operator has purely imaginary matrix elements, it follows that the first order coefficients are purely imaginary. The second derivative of the total SCF energy becomes:

(5.135)2ESCFxK(m)xL(m)=iαaα{Uaαiα(m);L(aα|hKsm|iα)+Uaαiα(m);L(iα|hKsm|aα)}+iαaβ{Uaβiα(m);L(aβ|hKsm|iα)+Uaβiα(m);L(iα|hKsm|aβ)}+iβaα{Uaαiβ(m);L(aα|hKsm|iβ)+Uaαiβ(m);L(iβ|hKsm|aα)}+iβaβ{Uaβiβ(m);L(aβ|hKsm|iβ)+Uaβiβ(m);L(iβ|hKsm|aβ)}

Examination of the three cases m=0,±1 leads to the following equations for the D-tensor components:

(5.136)DKL(0)=14S2μνPμν(0)xL(0)(μ|hK;SOC|ν)
(5.137)DKL(+1)=12(S+1)(2S+1)μν(μ|hK;SOC|ν)Pμν(1)xL(+1)
(5.138)DKL(1)=12S(2S1)μν(μ|hK;SOC|ν)Pμν(+1)xL(1)

Where a special form of the perturbed densities has been chosen. They are given in the atomic orbital basis as:

(5.139)Pμν(0)xL(0)=iαaαUaαiα(0);Lcμiαcνaα+iβaβUaβiβ(0);Lcμiβcνaβ
(5.140)Pμν(+1)xL(1)=iαaβUaβiα(1);LcμiαcνaβiβaαUaαiβ(1);Lcμaαcνiβ
(5.141)Pμν(1)xL(+1)=iαaβUaβiα(+1);Lcμaβcνiα+iβaαUaαiβ(+1);Lcμiβcνaα

The special form of the coupled perturbed equations are implemented in ORCArun as follows: The perturbation is of the general form hKs^m. The equations (5.136)-(5.141) and (5.142)-(5.147) below have been written in such a way that the spin integration has been performed but that the spin-dependent factors have been dropped from the right-hand sides and included in the prefactors of eqs. (5.136)-(5.138). The explicit forms of the linear equations to be solved are:

m=0:

(5.142)(εaα(0)εiα(0))UaαiαK(0)+aXjαbαUbαjαK(m){(bαiα|aαjα)(jαiα|aαbα)}=(aα|hK|iα)
(5.143)(εaβ(0)εiβ(0))UaβiβK(0)+aXjβbβUbβjβK(m){(bβiβ|aβjβ)(jβiβ|aβbβ)}=(aβ|hK|iβ)

m=+1:

(5.144)(εaα(0)εiβ(0))UaαiβK(+1)+aXjαbαUbβjαK(+1)(bβiβ|aαjα)aXbαjβUbβjαK(+1)(jβiβ|aαbα)=(aα|hK|iβ)
(5.145)(εaβ(0)εiα(0))UaβiαK(+1)+aXjβbαUbαjβK(+1)(bαiα|abjβ)aXbβjαUbβjαK(+1)(jαiα|aβbβ)=0

m=1:

(5.146)(εaβ(0)εiα(0))UaβiαK(1)+aXjβbαUbαjβK(1)(bαiα|abjβ)aXbβjαUbβjαK(1)(jαiα|aβbβ)=(aβ|hK|iα)
(5.147)(εaα(0)εiβ(0))UaαiβK(1)+aXjαbαUbβjαK(1)(bβiβ|aαjα)aXbαjβUbβjαK(1)(jβiβ|aαbα)=0

Note that these coupled-perturbed (CP) equations contain no contribution from the Coulomb potential or any other local potential such as the exchange-correlation potential in DFT. Hence, in the absence of HF exchange, the equations are trivially solved:

(5.148)UaαiαK(0)=(aα|hK|iα)εaα(0)εiα(0)
(5.149)UaβiβK(0)=(aβ|hK|iβ)εaβ(0)εiβ(0)
(5.150)UaαiβK(+1)=(aα|hK|iβ)εaα(0)εiβ(0)
(5.151)UaβiαK(+1)=0
(5.152)UaβiαK(1)=(aβ|hK|iα)εaβ(0)εiα(0)
(5.153)UaαiβK(1)=0

It is interesting that the “reverse spin flip coefficients” UaβiαK(+1) and UaαiβK(1) are only nonzero in the presence of HF exchange. In a perturbation expansion of the CP equations they arise at second order (V2/Δε2) while the other coefficients are of first order (V/Δε; V represents the matrix elements of the perturbation). Hence, these contributions are of the order of α4 and one could conceive dropping them from the treatment in order to stay consistently at the level of α2. These terms were nevertheless kept in the present treatment.

Equations (5.142)-(5.147) are referred to as CP-SOC (coupled-perturbed spin-orbit coupling) equations. They can be solved by standard techniques and represent the desired analogue of the CP-SCF magnetic response equations solved for the determination of the g-tensor and discussed in detail earlier [710]. It is readily confirmed that in the absence of HF exchange, eqs. (5.148)-(5.153) inserted into eqs. (5.136)-(5.141) lead back to a modified Pederson-Khanna type treatment of the SOC contributions to the D-tensor [206]. In the framework of the formalism developed above, the Pederson-Khanna formula can be re-written in the form:

(5.154)DKL(SOC;PK)=14S2iβ,aβ(ψiβ|hK;SOC|ψaβ)UaβiβL;(0)+14S2iα,aα(ψiα|hK;SOC|ψaα)UaαiαL;(0)14S2iα,aβ(ψiα|hK;SOC|ψaβ)UaβiαL;(1)14S2iβ,aα(ψiα|hK;SOC|ψaα)UaαiβL;(+1)

This equation was derived from second-order non-self-consistent perturbation theory without recourse to spin-coupling. For the special case of no Hartree-Fock exchange, the main difference to the treatment presented here is that the correct prefactors from eqs. (5.125)-(5.127) occur in front of the spin-flip contributions rather than ± 1/(4S2) in eq. (5.154). In the presence of HF exchange it is suggested that the consistent generalization of the PK method are eqs. (5.136)-(5.138) with the ± 1/(4S2) prefactors and this way the method has been implemented as an option into the ORCA program.

For completeness, the evaluation of the spin-spin term in the SCF case proceeds conveniently through:

(5.155)DKL(SS)=ge24α2S(2S1)μνκτ{PμναβPκταβPμκαβPνταβ}μν|r125{3r12,Kr12,LδKLr122}|κτ

as derived by McWeeny and Mizuno and discussed in some detail by Sinnecker and Neese.[708] In this reference it was found that DFT methods tend to overestimate the spin-spin contribution if the calculations are based on a spin-unrestricted SCF treatment. A much better correlation with experiment was found for open-shell spin restricted calculations. The origin of this effect proved to be difficult to understand but it was shown in ref [207] that in the case of small spin-contamination, the results of ROKS calculations and of those that are obtained on the basis of the spin-unrestricted natural orbital (UNO) determinant are virtually indistinguishable. It is therefore optionally possible in the ORCA program to calculate the spin-spin term on the basis of the UNO determinant.

5.24.3.4. ZFS Keywords

Table 5.16 %eprnmr block input keywords relevant for ZFS calculations

Keyword

Options

Description

dtensor

ss

Calculate spin-spin part

so

Calculate spin-orbit part

ssandso

Calculate both parts

DSOC

qro

quasi-restricted method for SOC, must be combined with !UNO

pk

Pederson-Khanna method for SOC

cp

coupled-perturbed method for SOC (default)

cvw

van Wüllen method for SOC

DSS

direct

directly use the canonical orbitals for the spin density

uno

use spin density from UNOs

printEuler

true/false

Whether to calculate Euler angles via orca_euler