AFSSH-NAMD calculation between ground state and S1 - Bad Overlap

Questions
#1

I am trying to run NAMD in Q-Chem to study the surface hopping from the S1 state of retinal to its ground state. To do so I have defined the $rem section as follows:

$rem
JOBTYPE aimd
EXCHANGE B3LYP
BASIS 6-31g
MEM_TOTAL 10000
CIS_N_ROOTS 1
SYMMETRY off
SYM_IGNORE true
CIS_SINGLETS true
CIS_TRIPLETS false
PROJ_TRANSROT true
FSSH_LOWESTSURFACE 0
FSSH_NSURFACES 2 ! hop between S0 and S1
FSSH_INITIALSURFACE 1 ! start on S1
AFSSH 1 ! decoherence
CALC_NAC true
AIMD_STEPS 1000 ! Typically more would be used
TIME_STEP 14
AIMD_SHORT_TIME_STEP 2
AIMD_TIME_STEP_CONVERSION 1 ! Do not alter time_step duration
AIMD_PRINT 1
AIMD_INIT_VELOC thermal
AIMD_TEMP 300
AIMD_INTEGRATION vverlet
FOCK_EXTRAP_ORDER 6
FOCK_EXTRAP_POINTS 12
$end

but I end up with a Bad Overlap error after 1 or 2 steps:
overlap matrix:
1.8957e-05 1.4558e-02 -6.8652e-03 -1.5981e-04
-4.0988e-04 4.8832e-04 -7.0997e-04 6.3056e-04
3.2727e-03 -6.1944e-04 -3.6709e-03 -5.0990e-03
-1.6554e-03 -7.6313e-03 -2.0576e-02 7.0119e-05

Q-Chem fatal error occurred in module aimdman/finite_difference_dc.C, line 754:

Ending Q-Chem run because of bad overlap.

I have tried already to increase the number of roots/fssh_nsurfaces (up to 4) and to reduce the time step ranging from 40 to 6. Despite this I always get the same error.

Interestingly, I didn’t get any errors when not including the ground state (FSSH_LOWESTSURFACE = 1) in the calculations, i.e. NAMD with surface hopping from S1 to higher states.

Is there any other suggestion about how to resolve this?

#2

Not sure. It is possible this is a “feature” (of your system), not a bug. Specifically the fact that it works when you exclude the ground state makes me suspicious. Conventional LR-TDDFT has topology problems in the vicinity of conical intersections involving ground state, leading to warped potential surfaces in a way that doesn’t affect intersections between excited states. (See here: https://doi.org/10.1063/5.0062757).

Beyond that, it would be useful if you could strip down your input file to the simplest example that manifests the problem, rather than something with 25 different REMs. What happens if you turn off Fock matrix extrapolation, multiple time steps, etc?

#3

I was running into this problem too. When I looked at the coordinates of the trajectory, I noticed that “Standard Nuclear Orientation” was very different between the initial SCF+TDDFT and the remaining trajectory. On a hunch I changed my coordinates in the job file to the trajectory orientation from the 1st step, and dynamics proceeded fine. Might be just luck, but I suspect Q-Chem has a bug (otherwise either the coordinates would stay in about the same frame, or I would see the jump regardless of starting center of mass).

Maybe give it a try and see.

#4

In that case maybe try SYM_IGNORE = TRUE? That prevents Q-Chem from rotating the coordinates into the Standard Nuclear Orientation and instead uses whatever coordinate frame is specified in $molecule section. If that calculation works okay then there may be an issue with this change of frame.

#5

Thanks John. It’s looking like luck though–the job that failed before is now running fine. Maybe some velocity assignments in this system are worse than others…

#6

Thank you so much John for the suggestion. But it seems I already had SYM_IGNORE = true in the input.

#7

Thank you so much for your suggestion. I tried starting the job with the coordinates from the trajectory orientation of the first step as you suggested, still the job failed with the bad overlap error after the first step. I will keep trying tuning different input parameters.

#8

Jumping back in on this, since I can’t yet rationalize my problem with bad overlap.

For context, I’m using Zhang & Herbert DOI:10.1063/5.0062757 as a starting point to get experience with the Q-Chem FSSH capabilities. I am trying out the spin-flip code, as open-shell low-spin systems might be of interest to us. Some other possible modifications as exemplified in the $REM section:

  jobtype aimd
  method b5050lyp
  basis 6-31g
  spin_flip true
  rpa false
  cis_n_roots 8
  max_cis_cycles 100
  symmetry off
  sym_ignore true
  proj_transrot true
  fssh_initialsurface 3
  fssh_lowestsurface 1
  fssh_nsurfaces 5
  afssh 0 ! decoherence correction
  calc_nac true
  aimd_steps 100 ! Typically more would be used
  time_step 10
  aimd_time_step_conversion 1 ! 1 means do not alter time_step duration
  aimd_print 1
  aimd_init_veloc thermal
  aimd_temp 300 # k
  aimd_integration vverlet
  fock_extrap_order 6
  fock_extrap_points 12

The job stopped at time step 34 (8 fs in) with the message

ERROR: very little overlap between states at the current
       geometry and states at the previous geometry. Consider
       decreasing the nuclear time step, or increasing the
       number of states included in this calculation.
              overlap matrix:
   0.9715   0.0034   0.0009  -0.0006  -0.0029
  -0.0021   0.0765  -0.0080  -0.0033   0.0015
  -0.0001   0.0069   0.4929  -0.0046  -0.0021
   0.0008   0.0024   0.0150   0.3185   0.0019
  -0.0034  -0.0021   0.0025  -0.0055   0.5515

I’ve already tried minimal nuclear timesteps in other runs and expanded the state space to 8 (evolving on 5 surfaces) here, and always see the same termination. So, thought I’d dive in.

I’m assuming the diagonal elements reflect the overlap between the current state and the corresponding state (i.e., ordered by energy) from the last time step. If so, I’d kind of expect a change in the SCF reference, a change in the TDDFT excitation makeup, or both. But what I see from the last few steps are

SCF energies:

 SCF   energy in the final basis set =     -249.4841937865
 SCF   energy in the final basis set =     -249.4832061774
 SCF   energy in the final basis set =     -249.4822000937

TDDFT (taking state 3 as an arbitrary data point):

--
 Excited state   3: excitation energy (eV) =    2.2342
 Total energy for state  3:                  -249.40209025 au
    <S**2>     :  0.2125
    D(   21) --> S(    1) amplitude = -0.3789
    S(    1) --> S(    1) amplitude =  0.5659 alpha
    S(    2) --> S(    1) amplitude = -0.2500 alpha
    S(    2) --> S(    2) amplitude =  0.6011 alpha
    S(    2) --> V(    1) amplitude = -0.1835 alpha

--
 Excited state   3: excitation energy (eV) =    2.2145
 Total energy for state  3:                  -249.40182298 au
    <S**2>     :  0.2149
    D(   21) --> S(    1) amplitude = -0.3910
    S(    1) --> S(    1) amplitude =  0.5565 alpha
    S(    2) --> S(    1) amplitude = -0.2588 alpha
    S(    2) --> S(    2) amplitude =  0.5924 alpha
    S(    2) --> V(    1) amplitude = -0.1941 alpha

--
 Excited state   3: excitation energy (eV) =    2.1941
 Total energy for state  3:                  -249.40156902 au
    <S**2>     :  0.2165
    D(   21) --> S(    1) amplitude =  0.4004
    S(    1) --> S(    1) amplitude = -0.5480 alpha
    S(    2) --> S(    1) amplitude = -0.2671 alpha
    S(    2) --> S(    2) amplitude =  0.5848 alpha
    S(    2) --> V(    1) amplitude = -0.2028 alpha

No big discontinuities; the only thing that sticks out from prior time steps (beyond these 3) is that the D(21) → S(1) and S(1) → S(1) excitations swapped sign. That seems like a phase change that shouldn’t affect the state overlaps though.

Is there some other signature that I should be looking for? Thanks for any help.

#9

I also found similar errors. It seemed that if S0 is not included in, similar job can ended normally. However when perform this job for S1 and S0, “bad overlap” occurs.