Spin-contamination

Questions
#1

Hi
I am optimizing a dimer system with SF-TDDFT. I ran a job for that keeping one monomer fix:

$molecule
0 3
coordinates
$end

$rem
BASIS = cc-pVDZ
EXCHANGE = omegaB97XD
JOB_TYPE = optimization
sts_mom = true
SCF_CONVERGENCE = 8
MAX_SCF_CYCLES = 700
MAX_CIS_CYCLES = 700
UNRESTRICTED = true
SPIN_FLIP = 1
CIS_N_ROOTS = 4
THRESH = 14
SYMMETRY_IGNORE = true
SYMMETRY = false
CIS_STATE_DERIVATIVE = 1
SOLVENT_METHOD = PCM
geom_opt_max_cycles = 100
$end

$solvent
DIELECTRIC 32.613000
OPTICALDIELECTRIC 1.765709
$end

$opt
FIXED
1:36 XYZ
ENDFIXED
$end

In 100 optimization cycles criteria were not converged so I took the last cycle coordinates and put a job with the same keywords again but I faced the issue with the values.
In the last optimization cycle (100), the situation was this:

        SF-DFT Excitation Energies

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

Excited state 1: excitation energy (eV) = -2.1398
Total energy for state 1: -1797.97035658 au
<S**2> : 0.0535
S( 2) → S( 1) amplitude = 0.9887 alpha

Excited state 2: excitation energy (eV) = 0.4743
Total energy for state 2: -1797.87429054 au
<S**2> : 2.0103
D( 139) → S( 1) amplitude = -0.6612
S( 2) → V( 1) amplitude = 0.7321 alpha

Excited state 3: excitation energy (eV) = 1.5510
Total energy for state 3: -1797.83472418 au
<S**2> : 0.1472
D( 139) → S( 1) amplitude = 0.6878
S( 2) → V( 1) amplitude = 0.6398 alpha

Excited state 4: excitation energy (eV) = 2.0505
Total energy for state 4: -1797.81636797 au
<S**2> : 1.0668
D( 135) → S( 1) amplitude = -0.4719
S( 2) → V( 2) amplitude = -0.2218 alpha
S( 2) → V( 3) amplitude = 0.8262 alpha

But in the second time optimization (I took the last cycle coordinates i.e., 100th optimized cycle coordinates) the situation of excited states is this:

        SF-DFT Excitation Energies

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

Excited state 1: excitation energy (eV) = -1.3635
Total energy for state 1: -1797.90126864 au
<S**2> : 0.1736
S( 2) → S( 1) amplitude = 0.9917 alpha

Excited state 2: excitation energy (eV) = 0.4626
Total energy for state 2: -1797.83416326 au
<S**2> : 1.1583
S( 2) → V( 2) amplitude = 0.9561 alpha
S( 2) → V( 3) amplitude = 0.1820 alpha

Excited state 3: excitation energy (eV) = 0.7540
Total energy for state 3: -1797.82345211 au
<S**2> : 1.2327
D( 136) → S( 1) amplitude = -0.9457

Excited state 4: excitation energy (eV) = 1.1579
Total energy for state 4: -1797.80860917 au
<S**2> : 1.0686
S( 2) → S( 2) amplitude = 0.9805 alpha

What should I have to do to resolve this issue?

#2

I’m not sure what your question is. SF-TDDFT often suffers from spin contamination, as described here:

The spin-adapted spin-flip (SASF) version ensures proper spin eigenstates but does not have gradients.

#3

Hi @jherbert
In starting I was not facing the spin contamination issue with my geometry. In the input, I gave only geom_opt_max_cycles = 100 but my geometry was not optimized within 100 cycles and the job got terminated. I resubmitted the job by taking the last cycle coordinates from the terminated job, I faced this spin contamination issue. Why is it so? How can I remove this?

#4

It is often the case that spin contamination is less severe at the Franck-Condon point (ground state geometry) but increases as one optimizes on a potential surface and encounters incipient bond break. For example, pulling apart a sigma bond is a good way to induce a singlet/triplet near degeneracy and that often exacerbates spin contamination. There’s no easy way to avoid it in methods that suffer from spin contamination.

#5

One can not simply “remove” spin contamination in the SF-TDDFT method, since it is an inherent problem/feature of SF-TDDFT.

To reduce the possible spin contamination, you can use ROKS (aka RODFT) wave function as the reference of SF-TDDFT, instead of using UKS (aka UDFT). The ROKS reference wave function is spin pure, but the SF-TDDFT will still have spin contamination (but smaller than that of UKS-based SF-TDDFT). As far as I know, Q-Chem does not support ROKS-based SF-TDDFT currently.

For your situation, my suggestion is just optimizing the geometry. If the <S**2> is 0.1 or 0.2, it can be viewed as a singlet state. After the geometry optimization is accomplished, you can use the SA-SF-DFT method in Q-Chem to perform a single-point calculation. Every electronic state is spin pure in this method.