X

Search Results

Searching....

7.3 Time-Dependent Density Functional Theory (TDDFT)

7.3.7 Spin-Orbit Coupling

(September 1, 2024)

7.3.7.1 SOCs Between TDDFT States

Several options for computing spin-orbit couplings (SOCs) between TDDFT states are available:

  1. (i)

    one-electron part of the Breit-Pauli Hamiltonian,

  2. (ii)

    one-electron SOC using scaled nuclear charges, and

  3. (iii)

    full SOC using the mean-field treatment of the two-electron part.

Options (i) and (ii) are are available for both TDA and RPA variants (including TDHF and CIS states), for restricted Kohn-Sham references only. Option (iii) is available for both restricted and unrestricted variants, but only within TDA. Calculations of SOC for SF-TDDFT are also possible within TDA using option (iii).

The implementation of one-electron SOC, options (i) and (ii), is based on the following. The SOCs are computed by evaluating matrix elements of the one-electron part of the Breit-Pauli Hamiltonian:

H^SO=-α022i,AZAriA3(𝐫iA×𝐩i)𝐬i (7.20)

where i denotes electrons, A denotes nuclei, α0=1/137.037 is the fine structure constant, and ZA is the bare positive charge on nucleus A. In the second quantization representation, the spin-orbit Hamiltonian in different directions can be expressed as

H^SO,x =-α022pqL~x,pq2(a^pa^q¯+a^p¯a^q) (7.21a)
H^SO,y =-α022pqL~y,pq2i(a^pa^q¯-a^p¯a^q) (7.21b)
H^SO,z =-α022pqL~z,pq2(a^pa^q-a^p¯a^q¯) (7.21c)

where L~^α=L^αr-3 for α{x,y,z} and L~α,pq are the matrix elements of this operator. The single-reference ab initio excited states (within the TDA) are given by

|ΦsingletI =i,asiIa(a^aa^i+a^a¯a^i¯)|ΦHF (7.22a)
|ΦtripletI,Ms=0 =i,atiIa(a^aa^i-a^a¯a^i¯)|ΦHF (7.22b)
|ΦtripletI,Ms=1 =2i,atiIaa^aa^i¯|ΦHF (7.22c)
|ΦtripletI,Ms=-1 =2i,atiIaa^a¯a^i|ΦHF (7.22d)

where siIa and tiIa are singlet and triplet excitation coefficients of the Ith singlet or triplet state respectively, with the normalization

ia(siIa)2=ia(tiIa)2=12. (7.23)

The quantity |ΦHF refers to the Hartree-Fock ground state. Thus the SOC constant from the singlet state to different triplet manifolds are

ΦsingletI|H^SO|ΦtripletJ,Ms=0=α022(i,a,bL~z,absiIatiJb-i,j,aL~z,ijsiIatjJa) (7.24)

or

ΦsingletI|H^SO|ΦtripletJ,Ms=±1=α0222(i,a,bL~x,absiIatiJb-i,j,aL~x,ijsiIatjJa)+α0222i(i,a,bL~y,absiIatiJb-i,j,aL~y,ijsiIatjJa). (7.25)

The SOC constant between different triplet manifolds can be obtained as

ΦtripletI,Ms=0|H^SO|ΦtripletJ,Ms=±1=α0222(i,a,bL~x,abtiIatiJb+i,j,aL~x,ijtiIatjJa)α0222i(i,a,bL~y,abtiIatiJb+i,j,aL~y,ijtiIatjJa) (7.26)

or

ΦtripletI,Ms=±1|H^SO|ΦtripletJ,Ms=±1=±α022(i,a,bL~z,abtiIatiJb+i,j,aL~z,ijtiIatjJa). (7.27)

Note that

ΦtripletI,Ms=0|H^SO|ΦtripletJ,Ms=0=0=ΦtripletI,Ms=±1|H^SO|ΦtripletJ,Ms=1. (7.28)

The total (root-mean-square) spin-orbit coupling is

ΦsingletI|H^SO|ΦtripletJ =(Ms=0,±1ΦsingletI|H^SO|ΦtripletJ,Ms2)1/2 (7.29a)
ΦtripletI|H^SO|ΦtripletJ =(Ms=0,±1ΦtripletI,Ms|H^SO|ΦtripletJ,Ms2)1/2. (7.29b)

For RPA states, the SOC constant can simply be obtained by replacing siIatjJb with Xi,tripletIaXj,tripletJb+Yi,singletIaYj,tripletJb and tiIatjJb with Xi,tripletIaXj,tripletJb+Yi,tripletIaYj,tripletJb.

The calculation of SOCs using effective nuclear charges, option (ii), is described in Section 7.10.20.4. The calculations of SOCs using option (iii), with a mean-field treatment of the two-electron part, is implemented following the algorithm described in Refs.  1028 Pokhilko P., Epifanovsky E., Krylov A. I.
J. Chem. Phys.
(2019), 151, pp. 034106.
Link
, 675 Kotaru S., Pokhilko P., Krylov A. I.
J. Chem. Phys.
(2022), 157, pp. 224110.
Link
and outlined in Section 7.10.20.4.

The SOC calculation is activated by $rem variable CALC_SOC: CALC_SOC = 1 activates option (i), CALC_SOC = 4 activates option (ii), and CALC_SOC = 2 activates option (iii).

Note:  Setting CALC_SOC = TRUE activates a one-electron calculation using old algorithm, i.e., option (i).

CALC_SOC

CALC_SOC
       Controls whether to calculate the SOC constants for EOM-CC, RAS-CI, CIS, TDDFT/TDA and TDDFT/RPA.
TYPE:
       INTEGER/LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not perform the SOC calculation. TRUE Perform the SOC calculation.
RECOMMENDATION:
       Although TRUE and FALSE values will work, EOM-CC code has more variants of SOC evaluations. For details, consult with the EOM section. For TDDFT/CIS, one can use values 1, 2, and 4, as explained above.

Examples 7.3.7.1, 7.3.7.1, and 7.3.7.1 illustrate calculations of SOCs for (SF)-TDDFT states using the above features. These calculations can also be carried out for CIS states by modifying METHOD appropriately.

The libwfa analysis of the spinless one-particle transition density matrices (as described in Ref. 1030 Pokhilko P., Krylov A. I.
J. Phys. Chem. Lett.
(2019), 10, pp. 4857–4862.
Link
) is implemented for the TDDFT and SF-TDDFT calculations. To activate this analysis, use the following: CALC_SOC = 2, STATE_ANALYSIS = TRUE, MOLDEN_FORMAT = TRUE, and NTO_PAIRS = N (see Section 10.2 for details of the libwfa package). Note:  This analysis differs from the NTO analysis of the regular transition density matrix between the singlet and triplet states.

Example 7.16  Calculation of one-electron SOCs for using TDDFT/TDA.

$comment
   This sample input calculates the spin-orbit coupling constants for water
   between its ground state and its TDDFT/TDA excited triplets as well as the
   coupling between its TDDFT/TDA singlets and triplets.  Results are given in
   cm-1.
$end

$molecule
   0 1
   H       0.000000    -0.115747     1.133769
   H       0.000000     1.109931    -0.113383
   O       0.000000     0.005817    -0.020386
$end

$rem
   EXCHANGE             b3lyp
   BASIS                6-31G
   CIS_N_ROOTS          4
   CIS_CONVERGENCE      8
   MAX_SCF_CYCLES       600
   MAX_CIS_CYCLES       50
   SCF_ALGORITHM        diis
   MEM_STATIC           300
   MEM_TOTAL            2000
   CIS_SINGLETS         true
   CIS_TRIPLETS         true
   CALC_SOC             true
   MAX_CIS_CYCLES       300
   INTEGRAL_SYMMETRY    false
   POINT_GROUP_SYMMETRY false
$end

Example 7.17  Calculation of full SOCs for water molecule including mean-field treatment of the two-electron part of the Breit-Pauli Hamilton and UHF/TDDFT/TDA.

$comment
Calculation of full SOCs for water molecule inlcuding mean-field treatment
of the two-electron part of the Breit-Pauli Hamiltonian and Wigner-Eckart theorem.
UHF/TDDFT/B3LYP/6-31G within the TDA.
$end

$molecule
0 1
H       0.000000    -0.115747     1.133769
H       0.000000     1.109931    -0.113383
O       0.000000     0.005817    -0.020386
$end

$rem
jobtype                 sp
unrestricted            true
method                  b3lyp
basis                   6-31G
cis_n_roots             4
cis_convergence         8
cis_singlets            true
cis_triplets            true
calc_soc                2
$end

Example 7.18  Calculation of SOCs for methylene using non-collinear SF-TDDFT/PBE0.

$comment
Calculation of SOCs for methylene using non-collinear SF-TDDFT/PBE0,
with tight convergence.
$end

$molecule
0 3
H1
C  H1 1.0775
H2 C  1.0775 H1 133.29
$end

$rem
method = pbe0
basis = cc-pvtz
scf_convergence = 12
cis_convergence = 12
THRESH = 14
cis_n_roots = 2
calc_soc = 2                    Compute full SOC with mean-field treatment of 2el part
WANG_ZIEGLER_KERNEL =  TRUE     Important for 1,1 diradicals
spin_flip = TRUE
$end

7.3.7.2 TDDFT+SOC for Spin-Adiabats

For a system exhibiting strong SOC, it can be extremely helpful to calculate the electronic eigenfunctions (called “spin-adiabats”) of the entire (Coulomb + SOC) Hamiltonian. Computation of such spin-adiabatic states has been implemented within the CIS or TDDFT/TDA framework, with the SOC treated by the one-electron Breit-Pauli Hamiltonian at the post-SCF level.

Within a CIS or TDDFT/TDA-SOC, the wavefunction is represented as

|Ψ=ϵiatia(ϵ)|Φia(ϵ) (7.30)

where ϵ is a label of spin, and

|Φia(ϵ){12(|Ψiαaα+|Ψiβaβ)if ϵ=singlet12(|Ψiαaα-|Ψiβaβ)if ϵ=triplet,ms=0|Ψiβaαif ϵ=triplet,ms=+1|Ψiαaβif ϵ=triplet,ms=-1 (7.31)

The spin-adiabatic wavefunctions in Eq. (7.30) are obtained as the eigenstates of following Hamiltonian having (4NONV)2 total matrix elements:

H=A+VSO+Vm (7.32)

In Eq. (7.32), A is the usual CIS-like TDDFT/TDA linear response matrix (including the one-electron and two-electron terms):

Φ𝒊𝒂|A|Φ𝒋𝒃 =F𝒂𝒃δ𝒊𝒋-F𝒋𝒊δ𝒂𝒃+𝒂𝒋|𝒊𝒃-cHF𝒂𝒋|𝒃𝒊
  +Ω𝒂𝒊𝒃𝒋+Egsδ𝒂𝒃δ𝒊𝒋 (7.33)

where the indices 𝒂,𝒃,𝒊,𝒋 are the spin-orbitals, F is the Fock matrix, Egs is the SCF energy, and Ω is the DFT exchange-correlation functional. More details can be found in Sec. 7.3.1 and Ref.  100 Bellonzi N. et al.
J. Chem. Phys.
(2020), 152, pp. 044112.
Link
.

The VSO in Eq. (7.32) is the one-electron SOC Hamiltonian appearing in Eq. (7.20), and the Vm=gμB𝐁𝐒 is the spin-Zeeman term. Here g=2 is assumed. (The coupling between orbital motion and magnetic field is not yet included.) In matrix form, VSO+Vm reads

Ψ𝒊𝒂|VSO+Vm|Ψ𝒋𝒃=V𝒂𝒃δ𝒊𝒋-V𝒋𝒊δ𝒂𝒃 (7.34)

The matrix elements of V are given by

Vpαqα =μBBzδpq-α024L~pqz (7.35)
Vpβqβ =-μBBzδpq+α024L~pqz (7.36)
Vpαqβ =(μBBxδpq-α024L~pqx)-i(μBByδpq-α024L~pqy) (7.37)
Vpβqα =(μBBxδpq-α024L~pqx)+i(μBByδpq-α024L~pqy) (7.38)

where α0 is the fine-structure constant and 𝐋~=𝐋/r3.

Note:  Currently, this CIS or TDDFT/TDA-SOC implementation supports only restricted SCF calculations (which would produce excited state singlets and triplets in the absence of SOC). Also, the ground state is not coupled to the excited states through SOC (or otherwise).

The code is triggered by setting $rem variable CIS_SOC = N, where N is the number of roots desired:

CIS_SOC

CIS_SOC
       Controls the roots of performing TDDFT/TDA-SOC calculation.
TYPE:
       INTEGER/LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not perform the calculation. N Solve the N lowest spin-adiabatic states of TDDFT/TDA-SOC.
RECOMMENDATION:
       Less than or equal to 4×CIS_N_ROOTS. TDDFT/TDA-SOC first performs a standard TDDFT/TDA calculation so as to generate an initial guess before rerunning the diagonalization to generate the spin adiabats. Therefore, it is a good idea to perform a stand-alone normal TDDFT/TDA calculation and generate the excited spin-diabats and check the desired range of energies, before generating the excited spin-adiabats.

To include a post-SCF magnetic field, the $bfield_postscf section must be initialized:

    $bfield_postscf
    B_x B_y B_z
    $end

Here, the Bx,By,Bz are the spatial components of the external magnetic field in Tesla. Their format must be float (i.e., “0” should be “0.0”).

Derivative coupling calculations are possible between spin-adiabats, but note that the value of the real and imaginary components may vary between different calculations – because the wavefunctions are complex-valued and may acquire a different global phase (gauge) for each calculation. The norm of derivative coupling should not change.

Examples 7.3.7.2, 7.3.7.2 and 7.3.7.2 illustrate the energy, gradient and derivative coupling calculations using TDDFT/TDA-SOC. In Example 7.3.7.2, a post-SCF external magnetic field is also included.

Example 7.19  TDDFT/TDA-SOC calculation of formaldehyde using the ωB97X functional for the lowest 7 spin-adiabats.

$rem
jobtype                 sp
exchange                wb97x
basis                   6-31g
cis_convergence         8
cis_n_roots             2
symmetry                false
sym_ignore              true
cis_soc                 7
$end

$molecule
0   1
C 0 0 0
H 0 0 1.101224
H 1.00171657 0 -0.45744749
O -1.00461935 0 -0.645642
$end

Example 7.20  Nuclear gradient of the 4th state in the TDDFT/TDA-SOC calculation of formaldehyde with the ωB97X functional.

$rem
jobtype                 force
exchange                wb97x
basis                   6-31g
cis_convergence         8
cis_n_roots             2
symmetry                false
sym_ignore              true
cis_soc                 4
cis_state_deriv         4
new_dft                 1 ! This is necessary
$end

$molecule
0   1
C 0 0 0
H 0 0 1.101224
H 1.00171657 0 -0.45744749
O -1.00461935 0 -0.645642
$end

Example 7.21  Derivative coupling between the 4th and 7th TDDFT/TDA-SOC states for a calculation of benzaldehyde with the ωB97X functional, in the presence of 10 Tesla external magnetic field along the z-direction.

$rem
jobtype                 sp
exchange                wb97x
basis                   6-31g
scf_convergence         10
cis_convergence         10
cis_n_roots             2
symmetry                false
sym_ignore              true
cis_soc                 7
calc_nac                true
cis_der_numstate        2
$end

$derivative_coupling
...
4 7
$end

$bfield_postscf
0.0 0.0 10.0
$end

$molecule
0   1
C -0.0827555811 -0.0021927111 0.0742168160
C 0.0163712008 -0.0019362341 1.3503524350
H 0.9426483579 -0.0021873514 1.9554892885
C -1.2462538108 -0.0012178895 2.1035447415
H -1.1928615664 -0.0009432042 3.1942205420
C -2.4029093145 -0.0008942470 1.4089112839
H -3.3571329386 -0.0003915665 1.9718859042
C -2.4418611557 -0.0012099958 0.0672877891
H -3.3674530317 -0.0009676703 -0.5174127708
C -1.1678601986 -0.0019316068 -0.6798528411
H -1.2448466539 -0.0022668669 -1.7738119236
C 1.1607172056 -0.0028394281 -0.8113481262
H 1.2563061018 -0.0031882755 -1.8250763954
O 2.4950317353 -0.0032860998 -0.1735593327
$end