Dihedral scan optimisation fails partway through due to too many bad steps

I’ve recently been fiddling about with potential energy scans, namely scanning through torsional angles, e.g. the inter-ring torsion of biphenyl.

Using this manual entry (the last example) as a starring point, and eventually working out from this post that I needed to increase
FIXED_COORD_FORCE_CONSTANT or it immediately fell over horribly, I’ve discovered what I would describe as an “unexpected feature”.

The number of optimisation steps is not reset on subsequent calls to opt(). Hence, my scans died fairly quickly after (the standard) 50 optimisation steps, which happened over, maybe, 10 or so steps along my torsional coordinate. I.e., GEOM_MAXITER is the maximum total number of opt steps and not the maximum per call to opt() (as I was expecting).

I’m not sure whether this is expected behaviour or not. In driver.py, it looks like the previous wfn is used as a guess for all steps where opt_iter > 1 (i.e., all but the first step). Presumably, this is also useful in a PES scan, assuming that subsequent geometries are not too different!

Is there any way of resetting the counter between opt() calls (or setting it to 2, so that the guess comes from the previous step), or is the only real option to set GEOM_MAXITER to something massive (that may or may not be obvious at the start of the run!)?

I was also wandering whether the gradients from the previous step are reused. Essentially the same PES calculation on PSI4 took, maybe, 5 opt steps at each new geometry, whereas another code (that will remain nameless) took only a couple of steps at each point.

Could you show me your input file? The last time I ran a scan, the counter reset as expected.

I don’t understand what you mean by “the gradients from the previous step are reused”. What would they be reused for?

Hmmmm… I stand corrected!

It actually fell over due to “too many bad steps”, even though I’m totally convinced that it failed much much earlier before I increased GEOM_MAXITER to several thousand (and I overwrote the original output file!). (Post title changed to reflect a more general problem.)

It could be an issue I’ve come across before where I’ve started a job over ssh and then logged out. Any(?) python-side print() functions cause an exception because the terminal has closed. (I keep on meaning to look into that properly, rather than just commenting out a few print() lines in driver.py, etc.!)

Anyway, this is an input file I was using (based on the example mentioned above):

memory 34GB

import numpy as np

molecule mol { 
 C     0.000000     0.000000     0.000000
 C     0.000000     0.000000     1.400000
 C     1.212436     0.000000     2.100000
 C     2.424871     0.000000     1.400000
 C     2.424871     0.000000     0.000000
 C     1.212436     0.000000    -0.700000
 H    -0.943102     0.000000     1.944500
 H     1.212436     0.000000     3.189000
 H     3.367973     0.000000     1.944500
 H     3.367973     0.000000    -0.544500
 H     1.212436     0.000000    -1.789000
 C    -1.299038     0.000000    -0.750000
 C    -2.511474     0.000000    -0.050000
 C    -3.723909     0.000000    -0.750000
 C    -3.723909     0.000000    -2.150000
 C    -2.511474     0.000000    -2.850000
 C    -1.299038     0.000000    -2.150000
 H    -2.511474     0.000000     1.039000
 H    -4.667011     0.000000    -0.205500
 H    -4.667011     0.000000    -2.694500
 H    -2.511474     0.000000    -3.939000
 H    -0.355936     0.000000    -2.694500
symmetry c1

set {
  basis cc-pvdz
  intrafrag_step_limit 0.1
# Needed or it falls over after a couple of steps!
  fixed_coord_force_constant = 160.

dihedrals = np.arange(0, 360)

PES = []

set optking {
  dynamic_level = 1

for phi in dihedrals:
  my_string = "2 1 12 13 " + str(phi)
  set optking fixed_dihedral = $my_string
  E = opt('mp2')
  PES.append((phi, E))

print("\n\tcc-pVDZ SCF energy as a function of phi\n")
for point in PES:
  print("\t%5.1f%20.10f" % (point[0], point[1]))

It eventually exits with:

	Energy change for the previous step:
		Projected    :        -0.7803579448
		Actual       :         0.0001661772
	The BAD_STEP_EXCEPTion handler:
	Energy has increased in a minimization.

	Dynamic level is 7.
	Consecutive backsteps is 1.
	The INTCO_EXCEPTion handler:
	Too many bad steps.
	Dynamic level is 7.
	exc.g_really_quit() is 0.

  **** Optimization has failed! (in 5 steps) ****
			 OPTKING Finished Execution 
	Removing binary optimization data file.
	Cleaning optimization helper files.

PsiException: Could not converge geometry optimization in 56 iterations.

Traceback (most recent call last):
  File "/home/laz/miniconda3/bin/psi4", line 287, in <module>
  File "<string>", line 59, in <module>
  File "/home/laz/miniconda3/lib//python3.6/site-packages/psi4/driver/driver.py", line 1155, in optimize
    raise OptimizationConvergenceError("""geometry optimization""", n - 1, wfn)

psi4.driver.p4util.exceptions.OptimizationConvergenceError: Could not converge geometry optimization in 56 iterations.

This is after optimising 181 points along the set of dihedrals, with each point taking 5 or fewer steps to converge.

I’ll dig into it more deeply to see if I can work out what’s failing at that point. (I suspect wrapping the opt() call in a try/except pair could be prudent anyway!)

I’ve set it going again using scf instead of mp2 (and not adjusting GEOM_MAXITER) to see if it fails at the same point…

Thanks for the update. At this point, this is definitely an issue for @Rollin_King.

Hmmmm… Using try / except catches the exception but it then fails for all subsequent steps.

I think the calculation proceeds as follows:

  • the starting geometry is taken and SCF solved and gradients calulated
  • the geometry is modified accordingly, including any fixed parameters being incorporated at this point
  • repeat until opt criteria are met
  • for each subsequent step along my dihedral scan, the initial geometry is the final geometry from the previous step
  • SCF and gradient calculated
  • then the fixed dihedral is applied
    Hence, once my calculation has reached something that opt() doesn’t like, the failing geometry is used in subsequent steps which also fail.

An alternative would be to take the geometry from the previous step and then modify the dihedral before SCF and gradient calculations for the new step. I thought this would be how it would be working but it looks like it isn’t. I.e., take previous geometry, modify fixed coordinate to new value, and then optimise from there.

I see that psi4.core.Molecule has a set_variable() function that looks like I would be able to use to alter the geometry before SCF and gradient calculations. However, it doesn’t look the format of the variable strings are documented anywhere, and anything I’ve tried has been rejected.

What form do these variables take (or where would I find that out?)?

Below is a very ill-conceived but minimal example. The trick is to insert variables into the geometry definition and then modify those with the .variablename syntax later. Yes, the syntax is confusing, but we’re in the process of rewriting the optimizer, and I’m not sure what the syntax will be like in the rewritten version. It’s also possible to use both Cartesian and ZMAT coordinates; see here for an example.

molecule water {
H 1 r1
H 1 r2 2 104.5

r1 = 1.0 
r2 = 1.01

dist = [1.0, 1.1, 1.2]

set basis def2-TZVP

for r in dist:
    r_string = "1 2 " + str(r)
    water.r1 = r # Resets the value of r1 for the current geometry
    set optking fixed_distance = $r_string # Fix the 1 2 distance to the specified value
    optimize("b3lyp") # Now optimize

Ahhh…fair enough. I knew I could use variables like that within a z-matrix but I was thinking it was going to let me target / modify internal coordinates created as part of the optimisation algorithm (or maybe define and modify my own internal coordinates).

Using a z-matrix variable as you suggest would probably do the job for me (will give it a go a bit later on if I get the chance…) but for systems with more than a few atoms, it can be tricky to define a z-matrix where modifying one variable doesn’t have unexpected effects on adjacent bits of the molecule!


No, the variables that set_variable affects are separate from the internal coordinates created by the optimizer.

Yes, I agree this is not an effective method for large systems, and what I proposed is absolutely a hack. I’ll defer to our optimizer expert on less hacky solutions.

I’ve bodged a solution! It’s fairly straightforward to use the openbabel python module to modify the relevant dihedral angle in between subsequent constrained opt() calls. (I do like the way psi4 makes it so easy to string things together in this way!)

Running happily so far…hopefully, it shouldn’t fall over due to a single failing structure.