Hello everyone,
I am working on bimolecular reaction dynamics and am interested in studying non-adiabatic intersystem crossing (ISC) between singlet and triplet states.
I would like to know “Is on-the-fly non-adiabatic intersystem crossing (singlet–triplet) dynamics for bimolecular reactions possible in Q-Chem?”
What is available is described here:
https://manual.q-chem.com/latest/Ch7.S3.SS9.html
If you don’t have access to Q-Chem currently, you can always try something out using IQmol.
@jherbert Thanks for the suggestion. We have a Q-Chem 6.1 license.
I want to perform quasi-classical trajectory (QCT) dynamics (NAMD) starting from a triplet ground state and accessing excited singlet states. Since the system has multiplicity 3 (open-shell), I used SPIN_FLIP = TRUE in my input to find the singlet states.
This is my input:
$molecule
0 3
H -0.71303900 -1.20333300 0.00000000
O -0.96077500 -2.14993200 0.00000000
N 0.00000000 0.98242700 0.00000000
O 1.04990500 1.44072500 0.00000000
$end
$rem
JOBTYPE OPT
METHOD B3LYP
BASIS 6-311++G**
$end
@@@
$molecule
read
$end
$rem
JOBTYPE Freq
METHOD B3LYP
BASIS 6-311++G**
$end
@@@
$molecule
read
$end
$rem
JOBTYPE AIMD
METHOD B3LYP
BASIS 6-311++G**
MAX_SCF_CYCLES 100
MAX_CIS_CYCLES 300
MEM_TOTAL 10000
MEM_STATIC 2000
CIS_N_ROOTS 5
AIMD_METHOD BOMD
TIME_STEP 10
AIMD_STEPS 1000
AIMD_INIT_VELOC quasiclassical
FSSH_LOWESTSURFACE 1
FSSH_NSURFACES 2
FSSH_INITIALSURFACE 1
SPIN_FLIP True
CIS_SINGLETS True
CIS_TRIPLETS True
CALC_SOC 2
CALC_NAC True
AIMD_TEMP 298
UNRESTRICTED True
SCF_CONVERGENCE 5
AIMD_PRINT 2
$end
SF-DFT Excitation Energies
Excited state 1: excitation energy (eV) = 0.5815
Total energy for state 1: -205.67492925 au
<S**2> : 0.9712
S( 2) → S( 1) amplitude = 0.2633 alpha
S( 2) → S( 2) amplitude = 0.9625 alpha
Excited state 2: excitation energy (eV) = 0.9989
Total energy for state 2: -205.65958977 au
<S**2> : 0.9208
S( 1) → S( 1) amplitude = 0.9429 alpha
S( 2) → S( 1) amplitude = -0.3018 alpha
Excited state 3: excitation energy (eV) = 1.1716
Total energy for state 3: -205.65324128 au
<S**2> : 0.1599
S( 1) → S( 1) amplitude = 0.3168 alpha
S( 2) → S( 1) amplitude = 0.9145 alpha
S( 2) → S( 2) amplitude = -0.2431 alpha
Excited state 4: excitation energy (eV) = 1.4082
Total energy for state 4: -205.64454564 au
<S**2> : 1.0114
S( 2) → V( 1) amplitude = 0.9951 alpha
Excited state 5: excitation energy (eV) = 2.4873
Total energy for state 5: -205.60489000 au
<S**2> : 1.0027
D( 11) → S( 1) amplitude = 0.9935
Further
To reduce spin contamination, I tried adding:
SASF_CIS = 1
and keeping UNRESTRICTED = False
However, this results in the following error:
“This JobType is not supported by SA-SF-DFT.”
Could you please suggest how to reduce spin contamination in this type of QCT NAMD simulation? Also, is there anything else I should modify to make this calculation more reliable?
SA-SF-TDDFT is not supported for nonadiabatic dynamics and the spin contamination is a real issue. There are state-tracking/state-following procedures that you can try, see the manual. Starting with Q-Chem 6.3 there is MRSF-TDDFT which suffers less from spin contamination.
@jherbert We only have a license for Q-Chem 6.1.
state-following is not work.
I tried EOM-SF-CCSD but still contamination
EOMSF-CCSD transition 1/A
S^2 calculation will be performed in double precision
EOMSF-CCSD transition 2/A
S^2 calculation will be performed in double precision
Excited state properties for EOMSF-CCSD transition 1/A
Dipole moment (a.u.): 0.891142 (X 0.285140, Y 0.844262, Z 0.007176)
R-squared (a.u.): 346.669757 (XX 69.306046, YY 265.003674, ZZ 12.360038)
Gauge origin (a.u.): (0.000000, 0.000000, 0.000000)
Angular momentum (a.u.) against gauge origin:
(X -0.000567i, Y 0.000242i, Z 0.025965i)
Traces of the OPDMs: Tr(AA) 12.000000, Tr(BB) 12.000000
<S^2> = 1.012614
Excited state properties for EOMSF-CCSD transition 2/A
Dipole moment (a.u.): 0.890287 (X 0.299213, Y 0.838472, Z 0.006935)
R-squared (a.u.): 346.439406 (XX 69.011589, YY 264.138891, ZZ 13.288925)
Gauge origin (a.u.): (0.000000, 0.000000, 0.000000)
Angular momentum (a.u.) against gauge origin:
(X -0.000547i, Y 0.000198i, Z 0.015503i)
Traces of the OPDMs: Tr(AA) 12.000000, Tr(BB) 12.000000
<S^2> = 1.002846
I then tried CASSCF with the Wigner–Eckart SOC
NOONs
0.0097
0.0266
0.0736
1.0000
1.0176
1.9269
1.9717
1.9739
kinetic energy: 204.7581705
num elec: 24.0000000
*********SPIN-ORBIT COUPLING JOB USING WIGNER–ECKART THEOREM BEGINS HERE*********
Calculating one electron integrals: Completed calculation of one electron integrals!
Calculating MF integrals: Completed calculation of MF integrals!
Calculating S^2 values of the spin-states: S^2 values computed!
Computing transition densities and evaluating SOCME with Libwe:
*********SPIN-ORBIT COUPLING JOB ENDS HERE*********
But no values print. I could not find any relevant keyword in the manual. Is casscf soc for namd ok?
Thanks
Spin-flip methods are well known for their spin contamination,
https://chemrxiv.org/doi/full/10.26434/chemrxiv.15000030/v1
I’m not sure that the SOCs are available for CASSCF and even more sure that Q-Chem’s CASSCF code is not set up for nonadiabatic MD simulations.
You could try a trivial example with state tracking in IQmol (time-limited to 10 min) to see whether it’s worth upgrading your version of Q-Chem.