Several options for computing spin-orbit couplings (SOCs) between TDDFT states are available:
one-electron part of the Breit-Pauli Hamiltonian,
one-electron SOC using scaled nuclear charges, and
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:
ˆHSO=-α202∑i,AZAr3iA(𝐫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
ˆHSO,x | =-α202∑pq˜Lx,pqℏ2(ˆa†pˆaˉq+ˆa†ˉpˆaq) | (7.21a) | ||
ˆHSO,y | =-α202∑pq˜Ly,pqℏ2i(ˆa†pˆaˉq-ˆa†ˉpˆaq) | (7.21b) | ||
ˆHSO,z | =-α202∑pq˜Lz,pqℏ2(ˆa†pˆaq-ˆ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
|ΦIsinglet⟩ | =∑i,asIai(ˆa†aˆai+ˆa†ˉaˆaˉi)|ΦHF⟩ | (7.22a) | ||
|ΦI,Ms=0triplet⟩ | =∑i,atIai(ˆa†aˆai-ˆa†ˉaˆaˉi)|ΦHF⟩ | (7.22b) | ||
|ΦI,Ms=1triplet⟩ | =√2∑i,atIaiˆa†aˆaˉi|ΦHF⟩ | (7.22c) | ||
|ΦI,Ms=-1triplet⟩ | =√2∑i,atIaiˆa†ˉaˆai|ΦHF⟩ | (7.22d) |
where sIai and tIai are singlet and triplet excitation coefficients of the Ith singlet or triplet state respectively, with the normalization
∑ia(sIai)2=∑ia(tIai)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
⟨ΦIsinglet|ˆHSO|ΦJ,Ms=0triplet⟩=α20ℏ2(∑i,a,b˜Lz,absIaitJbi-∑i,j,a˜Lz,ijsIaitJaj) | (7.24) |
or
⟨ΦIsinglet|ˆHSO|ΦJ,Ms=±1triplet⟩=∓α20ℏ2√2(∑i,a,b˜Lx,absIaitJbi-∑i,j,a˜Lx,ijsIaitJaj)+α20ℏ2√2i(∑i,a,b˜Ly,absIaitJbi-∑i,j,a˜Ly,ijsIaitJaj). | (7.25) |
The SOC constant between different triplet manifolds can be obtained as
⟨ΦI,Ms=0triplet|ˆHSO|ΦJ,Ms=±1triplet⟩=∓α20ℏ2√2(∑i,a,b˜Lx,abtIaitJbi+∑i,j,a˜Lx,ijtIaitJaj)α20ℏ2√2i(∑i,a,b˜Ly,abtIaitJbi+∑i,j,a˜Ly,ijtIaitJaj) | (7.26) |
or
⟨ΦI,Ms=±1triplet|ˆHSO|ΦJ,Ms=±1triplet⟩=±α20ℏ2(∑i,a,b˜Lz,abtIaitJbi+∑i,j,a˜Lz,ijtIaitJaj). | (7.27) |
Note that
⟨ΦI,Ms=0triplet|ˆHSO|ΦJ,Ms=0triplet⟩=0=⟨ΦI,Ms=±1triplet|ˆHSO|ΦJ,Ms=∓1triplet⟩. | (7.28) |
The total (root-mean-square) spin-orbit coupling is
⟨ΦIsinglet|ˆHSO|ΦJtriplet⟩ | =(∑Ms=0,±1∥⟨ΦIsinglet|ˆHSO|ΦJ,Mstriplet⟩∥2)1/2 | (7.29a) | ||
⟨ΦItriplet|ˆHSO|ΦJtriplet⟩ | =(∑Ms=0,±1∥⟨ΦI,Mstriplet|ˆHSO|ΦJ,Mstriplet⟩∥2)1/2. | (7.29b) |
For RPA states, the SOC constant can simply be obtained by replacing sIaitJbj with XIai,tripletXJbj,triplet+YIai,singletYJbj,triplet and tIaitJbj with XIai,tripletXJbj,triplet+YIai,tripletYJbj,triplet.
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
J. Chem. Phys.
(2019),
151,
pp. 034106.
Link
,
675
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
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
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
|Ψ⟩=∑ϵ∑iat(ϵ)ia|Φa(ϵ)i⟩ | (7.30) |
where ϵ is a label of spin, and
|Φa(ϵ)i⟩≡{1√2(|Ψaαiα⟩+|Ψaβiβ⟩)if ϵ=singlet1√2(|Ψaαiα⟩-|Ψaβiβ⟩)if ϵ=triplet,ms=0|Ψaαiβ⟩if ϵ=triplet,ms=+1|Ψaβiα⟩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
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-ℏα204˜Lzpq | (7.35) | ||
Vpβqβ | =-μBBzδpq+ℏα204˜Lzpq | (7.36) | ||
Vpαqβ | =(μBBxδpq-ℏα204˜Lxpq)-i(μBByδpq-ℏα204˜Lypq) | (7.37) | ||
Vpβqα | =(μBBxδpq-ℏα204˜Lxpq)+i(μBByδpq-ℏα204˜Lypq) | (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