Minimum-energy crossing point between the S0 and S1

Hello everyone,

I am trying to locate the minimum-energy crossing point (MECP) between the first excited state and the ground state of a molecular photoswitch (a dithienylethene) but so far all the attempts have not achieved the minimum. I am using spin-flip TDDFT (BHHLYP/cc-pVDZ). I have used the branching plane method and penalty function method. My initial structures have been the ground state and the highest energy barrier geometry obtained during a potential energy surface scan. When the calculations have achieved the maximum time allowed, I have taken the final structure and have submitted again the calculation but neither of them has converged to the minimum. Do you have any suggestions to locate this MECP? I am attaching the rem section for both methods. ​

Any suggestions are welcome.
Thanks!!

$rem
JOBTYPE = opt
MECP_OPT = true
MECP_METHODS = branching_plane
MECP_PROJ_HESS = true
GEOM_OPT_COORDS = 0
MECP_STATE1 = [0,1]
MECP_STATE2 = [0,2]
UNRESTRICTED = true
SPIN_FLIP = true
CIS_N_ROOTS = 5
CALC_NAC = true
CIS_DER_NUMSTATE = 2
SET_ITER = 150
CIS_S2_THRESH = 120
EXCHANGE = bhhlyp
BASIS = cc-PVDZ
MAX_SCF_CYCLES = 500
GEOM_OPT_MAX_CYCLES = 500
MEM_STATIC = 2000
MEM_TOTAL = 10000
SYMMETRY_IGNORE = true
DFT_D = EMPIRICAL_GRIMME
$end

$rem
JOBTYPE = opt
MECP_OPT = true
MECP_METHODS = penalty_function
MECP_STATE1 = [0,1]
MECP_STATE2 = [0,2]
UNRESTRICTED = true
SPIN_FLIP = true
CIS_N_ROOTS = 5
CALC_NAC = true
CIS_DER_NUMSTATE = 2
SET_ITER = 150
CIS_S2_THRESH = 120
EXCHANGE = bhhlyp
BASIS = cc-PVDZ
MAX_SCF_CYCLES = 500
GEOM_OPT_MAX_CYCLES = 500
MEM_STATIC = 2000
MEM_TOTAL = 10000
SYMMETRY_IGNORE = true
DFT_D = EMPIRICAL_GRIMME
$end

Hard to guess without a structure. Can you provide coordinates for your molecule?

Sure. I am attaching the coordinates of the ground state.

$molecule
0 3
C -1.9426968527 0.7739896311 -0.5676212039
C -0.7286599020 1.1930192979 -0.1541059399
S -2.1004724089 -0.9318748547 -0.9788487037
C -0.4080368434 -1.2599707797 -0.2742974746
C 0.2393270962 0.1258141043 -0.1476845937
C 1.5811208151 0.1663584945 -0.0085448696
C 2.4721040166 1.3245814749 0.3128446033
C 3.6646396227 0.6453981881 1.0413827374
C 3.7878156069 -0.7480157990 0.3604907498
C 2.3936678970 -1.0370919619 -0.1033768826
C 1.8860472638 -2.1926004737 -0.5882486315
C 0.4760839974 -2.1412614539 -1.1821545907
S -0.0639027212 -3.9144001523 -1.2618696393
C 1.5901454742 -4.4839571722 -0.9862545926
C 2.4688712143 -3.5011961919 -0.6603696417
Cl -3.3329815999 1.7774673524 -0.7324244159
C -0.5866153562 -1.8228466595 1.1500086791
F 1.8982474617 2.2791659528 1.0776144182
F 2.9370008629 1.9357256762 -0.8099760166
F 3.3265676295 0.4674632723 2.3366981280
F 4.7943032953 1.3594079819 0.9838926086
F 4.2849667590 -1.6721315646 1.2173487586
F 4.6664006947 -0.6630656934 -0.6732064551
H -0.4958512324 2.2246523341 0.0991039090
C 0.5981253719 -1.6151805973 -2.6266433962
C 1.8901992093 -5.9194528356 -1.1073479856
H 3.5065332518 -3.6928111113 -0.3929529097
C 0.8995068031 -6.8864687108 -0.8858953000
C 3.1859739244 -6.3421779355 -1.4417059797
C 3.4833394991 -7.6968858031 -1.5392592179
C 2.4923048506 -8.6512037044 -1.3087299814
C 1.2005140212 -8.2417327858 -0.9829315930
H -0.1127354397 -6.5765101599 -0.6175887085
H 3.9598749162 -5.6012659789 -1.6476704743
H 4.4940303596 -8.0092443899 -1.8066601609
H 0.4209062158 -8.9828044461 -0.8002366480
H 2.7260454913 -9.7141053723 -1.3889790732
H -1.1730051622 -1.1142237826 1.7503349239
H 0.3940305888 -1.9578626618 1.6297181278
H -1.1074114305 -2.7883384688 1.1325796885
H 0.9987601365 -0.5907026693 -2.6262892300
H -0.3747459098 -1.6145113407 -3.1339088329
H 1.2886405120 -2.2590782497 -3.1877341896
$end

I started running your input except that I made the basis set 6-31G* rather than cc-pvdz in order to get better throughput. (This is usually a good idea, at least at first, in order to get the gross “lay of the landscape”.) I’ve only run a few cycles but one thing I noticed right away is that for the first few optimization cycles the <S^2> threshold identifies excited states 2 and 3 as singlets whereas the first excited state is identified as triplet. However, this changes in opt cycle #6:

A threshold of <S**2> = 1.20 is used to identify singlet states
Excited state 1, <S^2> = 1.4059: is it a singlet? No.
Excited state 2, <S^2> = 0.2236: is it a singlet? YES.
Excited state 3, <S^2> = 0.9678: is it a singlet? YES

Now the first state is identified as triplet, although given the intermediate value of <S^2> that’s not altogether clear-cut and it is possible that the optimization algorithm is root-flipping. You would need to take a careful look at the orbitals involved in the excited state at this geometry in order to make that determination. The bottom line is that the spin contamination in spin-flip can mean that these MECP searches are not black-box. Another way to check is to use the spin-adapted spin-flip DFT (SA-SF-DFT) that is described in the manual, which lacks gradients (so can only be used to spot-check) but is free of spin contamination so you can get proper singlets as described here: https://aip.scitation.org/doi/full/10.1063/1.4937571

Thanks for the information!! I will have a look at the paper on spin-adapted.