Clarification on SF-DFT MECP Optimization and ⟨S²⟩ Values

Questions
#1

Dear Support Team,

I am new to SF-DFT and currently working on MECP optimization for S₀/S₁ geometries. I have a question regarding the CIS_S2_THRESH parameter. Based on my understanding, when CIS_S2_THRESH is set to ~1.20, if the value is less than 1.2, the state is considered a singlet.

However, in my output, I noticed that for excited states 2 and 3 (Ex:2 & Ex:3), the ⟨S²⟩ values are 0.5 and 0.6, respectively. According to theory, ⟨S²⟩ = 0 corresponds to a pure singlet and ⟨S²⟩ = 1 to a pure triplet. Given this, should the MECP geometries in my cases still be considered singlets?

Could you please clarify how the ⟨S²⟩ values in my output should be interpreted in the context of SF-DFT?

$molecule
read S0.txt
$end

$rem
JOBTYPE opt
MECP_OPT true
MECP_METHODS penalty_function
METHOD bhhlyp
SPIN_FLIP true
UNRESTRICTED true
MAX_CIS_CYCLES 300
GEOM_OPT_MAX_CYCLES 300
BASIS 6-31G(d,p)
CIS_N_ROOTS 5
STATE_FOLLOW true
MECP_STATE1 [0,1]
MECP_STATE2 [0,2]
CIS_S2_THRESH 120
CIS_MAX_CYCLES 500
$end

output

        SF-DFT Excitation Energies

(The first “excited” state might be the ground state)

Excited state 1: excitation energy (eV) = 0.8844
Total energy for state 1: -3358.05186503 au
<S**2> : 1.8961
D( 79) → S( 2) amplitude = 0.1674
S( 2) → S( 1) amplitude = 0.8718 alpha
S( 2) → S( 2) amplitude = -0.2787 alpha

Excited state 2: excitation energy (eV) = 1.0452
Total energy for state 2: -3358.04599168 au
<S**2> : 0.5551
D( 78) → S( 2) amplitude = 0.1642
D( 79) → S( 2) amplitude = -0.3898
S( 1) → S( 2) amplitude = 0.4238 alpha
S( 2) → S( 1) amplitude = 0.3471 alpha
S( 2) → S( 2) amplitude = 0.6002 alpha

Excited state 3: excitation energy (eV) = 1.0635
Total energy for state 3: -3358.04524463 au
<S**2> : 0.6952
D( 76) → S( 2) amplitude = 0.3636
D( 78) → S( 2) amplitude = -0.1604
D( 79) → S( 2) amplitude = -0.4458
D( 80) → S( 2) amplitude = 0.2994
D( 81) → S( 2) amplitude = -0.1634
D( 82) → S( 2) amplitude = -0.2033
S( 2) → S( 2) amplitude = -0.5261 alpha

Thank you for your help in advance.

#2

Hi Suhani,
Your understanding is largely correct except that <S^2> = S(S+1) (in a.u.), so for a triplet that’s S = 1 so <S^2> = 2. The CIS_S2_THRESH value is set between the singlet and triplet values.

The issue is that SF-TDDFT excited states can often be significantly spin-contaminated, especially at bond-breaking geometries (such as near a MECP). This is definitely reflected in your output. These are not pure spin states, they are mixtures, so there’s a question of how to classify them and that’s where the threshold (CIS_S2_THRESH) comes in. If you have questions about which state is which, you can run spin-adapted SF-TDDFT for single points (no gradients), which affords spin-pure states:
https://manual.q-chem.com/latest/Ch7.S2.SS3.html

#3

Thank you so much @jherbert