How to Verify if Optimized SF-TDDFT MECP (Penalty Function) is a True MECP or Artificial?

Aarzoo · September 9, 2025, 6:26am

Hi everyone,

I recently optimized a Minimum Energy Crossing Point (MECP) using the penalty-constrained search algorithm within Spin-Flip TDDFT (SF-TDDFT) in Q-Chem. The optimization successfully converged. However, I’m not sure whether the result represents a true conical intersection (MECP) or an artificial crossing.

My Questions:

How can I verify whether the optimized structure is a true MECP and not an artificial or incorrect topology?
Is there a way to compute or extract the CI branching plane vectors at the optimized MECP geometry using SF-TDDFT?

Input file:

$molecule
0 3
C       1.8197791310     2.8255337158     7.8326295405
H       1.3082850562     3.7809679608     8.0220619840
H       2.8964033392     2.9896932802     8.0227395371
C       1.6246780268     2.4517998650     6.3760003309
H       0.5421836536     2.3910199357     6.1680562282
C       2.2574584145     1.1636444922     5.9038684250
C       2.0799985494     0.8907811670     4.5368381892
H       1.5523909735     1.6357268541     3.9323850471
C       2.5341063392    -0.2718892425     3.9347125546
H       2.3771203077    -0.4423729755     2.8668438042
C       3.1633727738    -1.2272192076     4.7323128559
H       3.4955923538    -2.1792076388     4.3101507256
C       3.3457771097    -0.9734940794     6.0820610310
H       3.8018802929    -1.7501613104     6.6976003887
C       2.9546701944     0.2348481381     6.7073388890
C       3.2217298023     0.2964730245     8.1809272101
C       2.1611236279     0.6311165487     9.0507041847
C       1.9341653422    -0.1223758624    10.2812442587
H       2.1599800090    -1.1909456738    10.2303555654
C       0.8397684316     0.2215175570    11.1336703792
H       0.6593976690    -0.4056758674    12.0115494118
C       0.0638370528     1.3185942117    10.8902469329
H      -0.7421181643     1.6186398673    11.5610123050
C       0.3664030642     2.1091485844     9.7420694044
H      -0.1864760908     3.0449737119     9.6048556697
C       1.3526728674     1.8151753690     8.8368934330
C       6.0224464662    -2.7367194708     8.8007450354
H       6.6017920805    -3.6403705134     8.5611815134
H       4.9548648224    -3.0148805971     8.7900004574
C       6.4180426466    -2.2277081369    10.1861268755
H       7.3980101978    -1.7255838630    10.0948707431
C       5.4473986713    -1.2832675684    10.8676799831
C       5.3953956234    -1.3676949809    12.2265695041
H       6.0331250841    -2.0868578728    12.7451118259
C       4.5066975470    -0.5424647460    12.9934608548
H       4.5722934498    -0.5680038585    14.0863282662
C       3.6030703543     0.2729362496    12.4092752004
H       2.9445505530     0.9084469722    13.0051355349
C       3.5416006149     0.3819934669    10.9473334960
H       3.5187439571     1.4415953111    10.6309245718
C       4.5481534416    -0.3909751832    10.1506443612
C       4.4194477477    -0.2695558678     8.7640758794
C       5.5946254194    -0.4640131874     7.8573653192
C       5.9874700151     0.6080831779     7.0437393983
H       5.4364035325     1.5495923895     7.1002039292
C       7.0586974227     0.4872875208     6.1644573995
H       7.3486012064     1.3377466559     5.5423530173
C       7.7464824500    -0.7212576170     6.0736068955
H       8.5780101677    -0.8340416891     5.3733042417
C       7.3773152044    -1.7817313702     6.8970651155
H       7.9289849873    -2.7253549747     6.8523078940
C       6.3180583270    -1.6620828914     7.8004381889
H       1.9994524229     3.2774759736     5.7480521440
H       6.5765479200    -3.0788690635    10.8654099490
$end

$rem
BASIS  =  cc-pvdz
METHOD  =  LRC-wPBEh
MECP_OPT = true
mecp_methods = penalty_function
JOB_TYPE  =  optimization
MAX_CIS_CYCLES = 400
SCF_CONVERGENCE = 8
MAX_SCF_CYCLES = 400
GEOM_OPT_MAX_CYCLES = 300
THRESH = 14
UNRESTRICTED = true
SPIN_FLIP = true
CIS_N_ROOTS = 5
mecp_state1 = [0,1]
mecp_state2 = [0,2]
SYMMETRY_IGNORE = true
SYMMETRY = false
SOLVENT_METHOD  =  PCM
$end

$solvent
DIELECTRIC  7.5
OPTICALDIELECTRIC  1.40
$end

Here, some part of output:

 ---------------------------------------------------
            SF-DFT Excitation Energies              
 (The first "excited" state might be the ground state) 
 ---------------------------------------------------

 Excited state   1: excitation energy (eV) =    0.5667
 Total energy for state  1:                 -1156.42475579 au
    <S**2>     :  0.2854
    S(    1) --> S(    1) amplitude = -0.4912 alpha
    S(    2) --> S(    2) amplitude =  0.8304 alpha

 Excited state   2: excitation energy (eV) =    0.6331
 Total energy for state  2:                 -1156.42231420 au
    <S**2>     :  1.3597
    S(    1) --> S(    1) amplitude =  0.5492 alpha
    S(    2) --> S(    1) amplitude =  0.7371 alpha
    S(    2) --> S(    2) amplitude =  0.2983 alpha

 Excited state   3: excitation energy (eV) =    0.6986
 Total energy for state  3:                 -1156.41990756 au
    <S**2>     :  0.8402
    S(    1) --> S(    1) amplitude = -0.6390 alpha
    S(    2) --> S(    1) amplitude =  0.6478 alpha
    S(    2) --> S(    2) amplitude = -0.3572 alpha

 Excited state   4: excitation energy (eV) =    1.3104
 Total energy for state  4:                 -1156.39742343 au
    <S**2>     :  0.2680
    S(    1) --> S(    2) amplitude =  0.9528 alpha
    S(    2) --> S(    2) amplitude =  0.1825 alpha

 Excited state   5: excitation energy (eV) =    2.3889
 Total energy for state  5:                 -1156.35779144 au
    <S**2>     :  1.0991
    D(   99) --> S(    1) amplitude = -0.1686
    D(   99) --> S(    2) amplitude =  0.2368
    D(  101) --> S(    1) amplitude =  0.3510
    D(  101) --> S(    2) amplitude = -0.1993
    S(    1) --> V(    1) amplitude =  0.2802 alpha
    S(    2) --> V(    1) amplitude =  0.6895 alpha
    S(    2) --> V(    3) amplitude = -0.1936 alpha
 
A threshold of <S**2> = 1.20 is used to identify singlet states
	Excited state 1, <S^2> = 0.2854: is it a singlet? YES.
	Excited state 2, <S^2> = 1.3597: is it a singlet? No.
	Excited state 3, <S^2> = 0.8402: is it a singlet? YES.

 CI_Sigma   (Lagrange multiplier) = 3.500 
 CI_Alpha   (stability parameter) = 0.020 a.u.
 Lower state energy (state 1) = -1156.4247557881 a.u.
 Upper state energy (state 2) = -1156.4199075612 a.u.
 Energy gap = 0.004848 a.u. (0.132 eV)
 Final energy is -1156.422331674677

	******************************
	**  OPTIMIZATION CONVERGED  **
	******************************
 ----------------------------------------------------------------
             Standard Nuclear Orientation (Angstroms)
    I     Atom           X                Y                Z
 ----------------------------------------------------------------
    1      C       1.8198907272     2.8254866621     7.8326933111
    2      H       1.3085415319     3.7810049855     8.0220710095
    3      H       2.8965572864     2.9894396958     8.0227736407
    4      C       1.6247167265     2.4517310271     6.3760735816
    5      H       0.5422183472     2.3911178539     6.1681126429
    6      C       2.2573771060     1.1635127440     5.9039154968
    7      C       2.0799993292     0.8907616070     4.5368530425
    8      H       1.5524029762     1.6357346210     3.9324219442
    9      C       2.5342136927    -0.2718277212     3.9346459204
   10      H       2.3773616600    -0.4421873708     2.8667373620
   11      C       3.1634877298    -1.2271951500     4.7321976622
   12      H       3.4958123786    -2.1791143884     4.3099597043
   13      C       3.3457743093    -0.9735835869     6.0819837702
   14      H       3.8017693994    -1.7503250763     6.6975122349
   15      C       2.9545648214     0.2346800021     6.7073538127
   16      C       3.2216644190     0.2963553737     8.1809309749
   17      C       2.1610357939     0.6311168461     9.0507386037
   18      C       1.9340631387    -0.1223933563    10.2812492763
   19      H       2.1597912167    -1.1909876001    10.2303616677
   20      C       0.8398326418     0.2216271819    11.1337648835
   21      H       0.6594852279    -0.4054699981    12.0117143337
   22      C       0.0638994599     1.3187130759    10.8903227126
   23      H      -0.7420206473     1.6187726719    11.5611222685
   24      C       0.3664113421     2.1091710298     9.7421221283
   25      H      -0.1864693924     3.0449895471     9.6048663594
   26      C       1.3527285812     1.8151981557     8.8369915038
   27      C       6.0225531018    -2.7366973461     8.8007779246
   28      H       6.6020029403    -3.6402851375     8.5611864827
   29      H       4.9550175254    -3.0150120919     8.7900381130
   30      C       6.4180879027    -2.2276781240    10.1861699645
   31      H       7.3980365477    -1.7255178622    10.0949627235
   32      C       5.4473999482    -1.2832907481    10.8677220427
   33      C       5.3953808742    -1.3677006204    12.2266557781
   34      H       6.0330520444    -2.0868813105    12.7452474334
   35      C       4.5066223137    -0.5424580317    12.9934727323
   36      H       4.5721402107    -0.5680203006    14.0863439482
   37      C       3.6030311451     0.2729727478    12.4092241039
   38      H       2.9443985283     0.9084542086    13.0049848697
   39      C       3.5415728020     0.3819776302    10.9472279168
   40      H       3.5187725175     1.4415777078    10.6308543123
   41      C       4.5481335929    -0.3910036680    10.1506430112
   42      C       4.4194008401    -0.2696709587     8.7640861567
   43      C       5.5945374352    -0.4640455458     7.8573473820
   44      C       5.9872855024     0.6080458680     7.0436627324
   45      H       5.4361022159     1.5494896594     7.1000673077
   46      C       7.0584979961     0.4873019805     6.1643569454
   47      H       7.3482930538     1.3377446543     5.5421787102
   48      C       7.7463926459    -0.7211817433     6.0735495245
   49      H       8.5779308631    -0.8339128242     5.3732503275
   50      C       7.3773252106    -1.7816565470     6.8970593273
   51      H       7.9290490282    -2.7252491038     6.8522999007
   52      C       6.3180845631    -1.6620604474     7.8004481490
   53      H       1.9996985300     3.2773313867     5.7481499617
   54      H       6.5765527778    -3.0788655735    10.8654382505

jherbert · September 9, 2025, 10:37am

I’m not sure what you mean by an “artificial crossing”. By construction, spin-flip methods have correct topology at a conical intersection because both target states are in the response space. See, e.g.:

The g and h vectors can be extracted from the output although they aren’t labeled as such. The h vector is the nonadiabatic coupling vector, which can be computed and printed as indicated in the manual,
https://manual.q-chem.com/5.4/subsec_deriv-couple-jobcontrol.html
The g vector is the difference of the gradients on the two states in question, either one of which can be printed using JOBTYPE = FORCE in conjunction with CIS_STATE_DERIV to tell Q-Chem which adiabatic electronic state you want.

Aarzoo · September 10, 2025, 6:20am

Thank you for your response, @jherbert

By artificial crossing, I meant an accidental degeneracy, rather than a true MECI.

I was wondering whether the penalty-constrained optimization algorithm ensures that the optimized point is a true conical intersection with proper topology, or whether it could converge to such an accidental degeneracy that doesn’t represent a physically meaningful MECI.

jherbert · September 10, 2025, 12:37pm

Any degeneracy between two states that are in the response space must have the topology of a conical intersection, that’s implicit in the material in the book chapter cited above. (If you don’t have access to it, you can contact me privately for a copy.)