SF-TDDFT optimization job

Questions
#1

I tried to use SF-TDDFT to locate a TS of C=C bond rotation.
IQmol shows that there includes over 200 geometries (exactly the same). But the first job in input file is freq job which should only include one geometry. It also shows that total energy keeps oscillating.
Output file only shows “Running Job 1” and SCF is converged.
Thank you for your attention!

image
image1440×857 110 KB

Below is my input:

$molecule
    1 3
    C      -2.34670090      2.18774046     -0.13639317
    C      -3.84941067      2.15887725      0.16964645
    H      -4.03756893      1.84100530      1.20723421
    H      -4.27265805      3.16567783      0.04291284
    H      -4.39789125      1.48616254     -0.50611988
    C      -1.68447025      3.28315842      0.70527030
    H      -1.82535991      3.09694204      1.78109040
    H      -0.60570422      3.33969536      0.50352715
    H      -2.12719006      4.26085500      0.46543670
    C      -2.12986013      2.50266517     -1.62107132
    H      -1.06091177      2.59242920     -1.86738468
    H      -2.58068201      1.74967350     -2.28404096
    H      -2.59735229      3.46883652     -1.85992157
    C      -1.71542906      0.82363564      0.26387076
    H      -1.81977149      0.74631333      1.35759011
    C      -2.45531376     -0.36329497     -0.37300413
    H      -2.55924602     -0.22307705     -1.45757520
    H      -3.45127513     -0.49717068      0.06331145
    N      -0.29730476      0.70628576     -0.02363026
    H      -0.71741177     -1.37137932     -0.41306818
    N      -1.69134842     -1.63343169     -0.18231473
    C      -2.11950534     -2.70675255     -1.14113227
    H      -2.02022912     -2.27131413     -2.14286594
    H      -1.38564327     -3.51630335     -1.05704129
    C      -1.67474706     -2.06127449      1.25835964
    H      -1.16128728     -1.25798899      1.80141888
    H      -2.71913138     -2.08880338      1.59549774
    C      -0.96798011     -3.38232561      1.48514994
    H      -1.52870845     -4.23330746      1.07869969
    H       0.04189379     -3.38060100      1.05133462
    H      -0.86755416     -3.53397920      2.56692360
    C      -3.52934872     -3.20190745     -0.89184506
    H      -4.27347339     -2.40170363     -0.99795954
    H      -3.75881317     -3.96972941     -1.64121851
    H      -3.63849643     -3.66249220      0.09914565
    H      -0.01738734      0.96931859     -0.96675210
    C       0.67084141      0.91414632      0.96605985
    C       2.01161356      1.36167801      0.55504807
    C       3.01052758      0.40369586      0.16332932
    C       4.23177347      0.79451007     -0.44650337
    C       2.80304342     -0.98718365      0.36192721
    C       5.17669869     -0.14874454     -0.83104243
    H       4.43152531      1.84964041     -0.63435545
    C       3.75508274     -1.92207433     -0.01872871
    H       1.87875341     -1.31187033      0.84196618
    C       4.95025022     -1.51265093     -0.62090132
    H       6.10295846      0.18221470     -1.30404047
    H       3.56746497     -2.98321709      0.15515709
    H       5.69602612     -2.24858312     -0.92402646
    H       0.29337490      1.07975312      1.97949923
    C       2.25999072      2.83983395      0.43581221
    H       1.62569566      3.39105271      1.14325432
    H       3.30663687      3.10691387      0.63785551
    H       2.01765477      3.21266191     -0.57698196
$end

$rem
    jobtype freq
    method bhhlyp
    basis def2-svp
    spin_flip true
    unrestricted true
    cis_n_roots 5
    cis_triplets false
    cis_state_deriv 1
	solvent_method pcm
$end

$pcm
	theory iefpcm
$end

$solvent
	dielectric 35.9
	opticaldielectric 1.803649
$end

@@@

$molecule
    read
$end

$rem
    jobtype ts
    method bhhlyp
    basis def2-svp
    geom_opt_hessian read
    geom_opt_max_cycles 128
    spin_flip true
    unrestricted true
    cis_n_roots 5
    cis_triplets false
    cis_state_deriv 1
	solvent_method pcm
$end

$pcm
	theory iefpcm
$end

$solvent
	dielectric 35.9
	opticaldielectric 1.803649
$end

@@@

$molecule
    read
$end

$rem
    jobtype freq
    scf_guess read
    method bhhlyp
    basis def2-svp
    spin_flip true
    unrestricted true
    cis_n_roots 5
    cis_triplets false
    cis_state_deriv 1
	solvent_method pcm
$end

$pcm
	theory iefpcm
$end

$solvent
	dielectric 35.9
	opticaldielectric 1.803649
$end

@@@

$molecule
    read
$end

$rem
    jobtype rpath
    scf_guess read
    method bhhlyp
    basis def2-svp
    spin_flip true
    unrestricted true
    cis_n_roots 5
    cis_triplets false
    cis_state_deriv 1
    rpath_max_cycles 50
	solvent_method pcm
$end

$pcm
	theory iefpcm
$end

$solvent
	dielectric 35.9
	opticaldielectric 1.803649
$end
#2

It’s always a good idea to look at the output file itself, not just IQmol, as there’s a lot more information in the output file. In this case, what you would discover is that analytic frequencies are not available for spin-flip so your first job is being evaluated by finite difference. There’s a warning to this effect at the top of the output file:

!!********************************************.
!!******************Warning*******************.
!!Desired Analytical derivatives not available.
!!Finite difference job might take a long time.
!!********************************************.
!!********************************************.

Also note that analytic frequencies are not available for IEF-PCM solvation model, though they are available for C-PCM. I suggest that you try this job again with the the frequency calculation (i.e., start with the TS search). Please note that SF-TDDFT can be subject to state-switching due to spin contamination so you will need to be on the lookout for that.

#3

Thank you for your advice!
I started the job from optimization for TS and the problem was solved.
This time the output shows that:

 ---------------------------------------------------
 Calculating Relaxed Density
 ---------------------------------------------------
 Iter    Rts Conv    Rts Left    Ttl Dev     Max Dev
 ---------------------------------------------------
   1         0           5      0.000538    0.000163
   2         0           5      0.000577    0.000179
   3         0           5      0.000681    0.000252
   4         0           5      0.000687    0.000233
   5         0           5      0.000463    0.000135
   6         0           5      0.000307    0.000080
   7         0           5      0.000255    0.000063
   8         0           5      0.000258    0.000068
   9         0           5      0.000256    0.000077
  10         0           5      0.000177    0.000056
  11         0           5      0.000129    0.000040
  12         0           5      0.000128    0.000042
  13         0           5      0.000115    0.000037
  14         0           5      0.000102    0.000033
  15         0           5      0.000070    0.000023
  16         0           5      0.000053    0.000017
  17         0           5      0.000051    0.000017
  18         0           5      0.000046    0.000015
  19         0           5      0.000035    0.000011
  20         0           5      0.000025    0.000008
  21         0           5      0.000018    0.000005
  22         0           5      0.000017    0.000005
  23         0           5      0.000015    0.000005
  24         0           5      0.000012    0.000003
  25         0           5      0.000010    0.000003
  26         0           5      0.000007    0.000002
  27         3           2      0.000003    0.000002
  28         3           2      0.000003    0.000001
  29         3           2      0.000002    0.000001
  30         3           2      0.000002    0.000001

 Q-Chem fatal error occurred in module 0, line  194:

 Out of Iterations- IterZ

It seems that the max iterations for relaxed density calculation is 30. However, the iteration detail really shows a trend of convergence. I guess that convergence criterion could be met if larger number of iterations given. Is there any keyword can be used for this?

#4

It looks like the convergence threshold is hard-coded to 1e-6, however the max number of cycles can be changed (MAX_CIS_CYCLES or CIS_MAX_CYCLES in $rem – both are aliased to the same internal variable, which defaults to 30).

#5

Thank you jherbert! I tried it and succeeded. Also as you mentioned before, spin contamination really disturbs the calculation so I chose to use spin-adapted SF-TDDFT. For the same molecule, after 120 optimization steps, the structure only changed a little. Total energies of S0 and S1 did not have notable change as well. (also energy gap between S0 and S1, as shown in the screenshot below)
Does this mean that the initial structure is close to the CI point?

image
image845×474 23 KB

#6

The SA-SF-TDDFT method does not have analytic gradients. Usually there is a block of text at the top of the Q-Chem output to tell you this (i.e., a warning that Q-Chem is going to try to evaluate the gradients by finite difference). However, I have just checked and for some reason this doesn’t get printed for SA-SF though it seems like it’s going through the finite-diff procedure anyway. I think this is a bug, and should have at least printed the warning. It also leads me to suspect that the very small energy changes you are seeing are not optimization steps but rather the small atomic displacements that are being used to evaluate the gradient.

UPDATE: I’m now certain that the gradients are evaluate finite-diff, it’s just that the warning at the top is missing. (There is a warning after the first SCF.) Geometry optimization does seem to work, in a simple test job, but you do need to be sure that you’re grabbing energies that correspond to optimization steps not to finite-diff displacement steps, and the output isn’t printed in a way that makes this easy. What you can do is to set MOLDEN_FORMAT=TRUE and then you will get a MolDen-formatted output file at the end that allows you to step through the geometry optimization (or just read out the energies at each step). The formatted checkpoint file used by IQmol doesn’t seem to read this properly.

#7

For the Q-Chem 6.0 release, I added the usual warning at the top of the output about finite-diff gradients. Q-Chem prints this whenever the user requests JOBTYPE=OPT with a method that does not have analytic gradients, but this message did does not get printed for SA-SF-TDDFT. That will be fixed in 6.0