Spin-Orbit coupling with SF-TDDFT

Dear all,

Is it possible to compute the spin-orbit coupling constant (SOCC) between the S0 and T1 states ? In my calculations, I set “0 3” at the beginning so my “ground state” is the T1 and my first excited state is the singlet ground state (S0).

When I tried to compute the SOCC, I obtained values of SOCC but only between the excited states, not with the ground state.

Thank you in advance !

If spin contamination is small, one of the spin-flipped states corresponds to the same triplet with zero spin projection. It is printed in one the “excited states”. The excitation energy is not strictly zero since the method is not spin adapted; in most cases it can be identified cleanly. So most likely the couplings that you seek are already in the output. If you post your output here, we could help identifying this state.

Hello Pavel,

Thank you for your answer ! Here is my input :

$molecule
0 3
C 2.5966427852 -1.2393492541 -1.0784209758
C 1.5220381266 -0.3772216297 -0.9065356038
C 1.6105825520 0.6433461108 0.0475429007
C 2.7763544628 0.8026336343 0.8082420550
C 3.8398801228 -0.0575595086 0.6165646981
C 3.7531765490 -1.0820557447 -0.3265728871
H 2.5267508339 -2.0453491755 -1.8072421342
H 0.6027320544 -0.4911963779 -1.4764328109
H 2.8091906124 1.6028763341 1.5450295564
H 4.7481966034 0.0641157018 1.2037892847
H 4.5910026797 -1.7630090528 -0.4659678738
N 0.5045242746 1.4285371024 0.3585204696
N -0.5006068750 1.4189124433 -0.4068150802
C -1.6092046713 0.6472707078 -0.0720281171
C -2.7753566269 0.7882678347 -0.8358264729
C -1.5231998581 -0.3450503298 0.9116910256
C -3.8419561037 -0.0616172545 -0.6172558572
H -2.8060802388 1.5657124367 -1.5967316505
C -2.6009099028 -1.1974252193 1.1105170490
H -0.6034930745 -0.4458318677 1.4833771069
C -3.7579596902 -1.0579259464 0.3559171740
H -4.7505667016 0.0460316975 -1.2067298330
H -2.5331432601 -1.9812425107 1.8633562809
H -4.5984424602 -1.7307891268 0.5168515430
$end

$rem
BASIS = 6-311G**
CIS_N_ROOTS = 4
CIS_MAX_CYCLES = 200
CIS_CONVERGENCE = 6
EXCHANGE = M11
GUI = 0
JOB_TYPE = SP
SCF_CONVERGENCE = 6
SCF_MAX_CYCLES = 200
SPIN_FLIP = 1
CALC_SOC = 2
$end

Here the output of the spin-flip :

Excited state 1: excitation energy (eV) = -0.0012
Total energy for state 1: -572.48746231 au
<S**2> : 0.1099
D( 44) → S( 1) amplitude = 0.2560
S( 1) → S( 2) amplitude = 0.2513 alpha
S( 2) → S( 1) amplitude = 0.8831 alpha
S( 2) → V( 3) amplitude = 0.1516 alpha

Excited state 2: excitation energy (eV) = 0.0888
Total energy for state 2: -572.48415504 au
<S**2> : 2.1657
D( 45) → S( 1) amplitude = 0.2420
S( 1) → S( 1) amplitude = 0.6443 alpha
S( 2) → S( 2) amplitude = 0.6187 alpha

Excited state 3: excitation energy (eV) = 0.7423
Total energy for state 3: -572.46013920 au
<S**2> : 0.1452
D( 44) → S( 2) amplitude = 0.2318
D( 45) → S( 1) amplitude = -0.3635
S( 1) → S( 1) amplitude = -0.5350 alpha
S( 2) → S( 2) amplitude = 0.6323 alpha
S( 2) → V( 6) amplitude = -0.1605 alpha

Excited state 4: excitation energy (eV) = 1.5175
Total energy for state 4: -572.43165418 au
<S**2> : 0.1423
D( 44) → S( 1) amplitude = -0.1844
D( 45) → S( 2) amplitude = 0.3374
S( 1) → S( 2) amplitude = 0.8263 alpha
S( 1) → V( 6) amplitude = -0.1819 alpha
S( 2) → S( 1) amplitude = -0.1968 alpha

So if I understood well, I should look at the SOC between excited state 1 and 2 for the SOCC between T1 and S0 ?

Here the ouput for SOC :

State A: Root 1
State B: Root 2

Analysing Sz ans S^2 of the pair of states…
Ket state: Computed S^2 = 0.109897241 will be treated as 0.000000000 Sz = 0.000000000
Bra state: Computed S^2 = 2.165674754 will be treated as 2.000000000 Sz = 0.000000000
Clebsh-Gordan coefficient: <0.000,0.000;1.000,0.000|1.000,0.000> = 1.000


One-electron SO (cm-1)
Reduced matrix elements:
<S|| Hso(L-) ||S’> = (-7.026738,29.142119)
<S|| Hso(L0) ||S’> = (0.000000,0.000772)
<S|| Hso(L+) ||S’> = (-7.026738,-29.142119)

SOCC = 42.394296

Actual matrix elements:
|Sz=0.00>
<Sz=-1.00|(-7.026738,-29.142119)
<Sz=0.00|(0.000000,0.000772)
<Sz=1.00|(-7.026738,29.142119)


Mean-field SO (cm-1)
Reduced matrix elements:
<S|| Hso(L-) ||S’> = (-4.024831,16.978996)
<S|| Hso(L0) ||S’> = (0.000000,0.000448)
<S|| Hso(L+) ||S’> = (-4.024043,-16.971654)

Singlet part of <S|| Hso(L0) ||S’> = (-0.000000,-0.000001) (excluded from all matrix elements)
L-/L+ Averaged reduced matrix elements:
<S|| Hso(L-) ||S’> = (-4.024437,16.975325)
<S|| Hso(L+) ||S’> = (-4.024437,-16.975325)

SOCC = 24.672160

Actual matrix elements:
|Sz=0.00>
<Sz=-1.00|(-4.024437,-16.975325)
<Sz=0.00|(0.000000,0.000448)
<Sz=1.00|(-4.024437,16.975325)

Should I take the 24.67 (Mean-field SO) or the 42.39 (One-electron SO) ? I’m not sure which one to use.

Thank you in advance !

Yes, it looks like the “Excited state 2” is the same triplet as the reference state, but with Sz=0. Spin contamination is not too big to distort the picture. I usually confirm the character by looking at NOs.

The Breit–Pauli spin-orbit operator has one and two-electron terms. The one-electron term overestimates the coupling, especially for light elements. The spin-orbit mean-field approximation captures the missing treatment of the two-electron term reasonably well. If your interest is in the overall spin-orbit coupling constant between the singlet and triplet states, the mean-field estimate is more appropriate.

Thank you for your answer ! It helps me a lot.

It might sounds like a naive question but, what is “NOs” ? How do you look at it ?

Thank you again !

NOs = natural orbitals, which are the eigenvectors of the one-particle density matrix expressed through atomic orbitals. They provide a compact representation of the wave function. That’s why they are very convenient to for analysis of electronic states. Q-Chem can write them either to fchk or to molden format, where they can be visualized by IQmol, Jmol, Gabedit, etc.

Thank you, all of your answers were very helpful !