Problems with geometry convergence in dihedral scan

I am a first time undergraduate user of Psi4 and am trying to collect information on the potential energy vs dihedral plot of phenol. However, I keep getting problems with convergence at some dihedral angles, even after expanding the geom_maxiter to 250 steps. Additionally, when I try to add damping to my scf calculations, I still am unable to converge. I have copied my input file as well as my error messages and would appreciate your help!

molecule phenol {
0 1
C -1.40249 -0.06649 0.04385
C -0.76579 1.16517 0.00350
C 0.62963 1.22075 -0.05185
C 1.45228 0.08210 -0.06982
C 0.75327 -1.13575 -0.02720
C -0.64037 -1.22534 0.02850
O 2.72273 0.14823 -0.12041
H -2.43229 -0.12009 0.08486
H -1.31956 2.03594 0.01430
H 1.08021 2.14873 -0.08111
H 1.29850 -2.01190 -0.03759
H -1.10000 -2.14886 0.05806
H 3.14772 0.98287 -0.14750
}
set {
basis cc-pvdz
intrafrag_step_limit 0.1
}

set scf {
soscf_max_iter 35
basis_guess true
damping_percentage 30
}

dihedrals = [i for i in range(0,185, 5)]
PES = []

for phi in dihedrals:
my_string = "5 4 7 13 " + str(phi)
set optking geom_maxiter = 250
set optking fixed_dihedral = $my_string
print(phi)
E = optimize(‘scf’)
print(E)
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]))

And my error message is:

psi4.driver.p4util.exceptions.SCFConvergenceError: Could not converge SCF iterations in 100 iterations.

Please let me know if you have any questions or would like to see more of my output to better address this problem!

Could you upload the output file, or make it otherwise available?

In addition, please read our guidelines on asking questions. Your input file needs backticks, and we need Psi4 version numbers.

Yes, I have saved my output file (phenol_damp.out) here on google drive: phenol_damp.out.

And I am using the Psi4 1.3.2 release! Thank you for your help

First, I find it helpful to set ensure_bt_convergence true when I see “Could not converge backtransformation.” appearing in the geometry optimization section. The optimizer needs to perform nonlinear transformation between internal coordinates and Cartesian coordinates, and sometimes it needs to spend some extra effort to perform the transformation correctly. It adds negligible computational cost within an optimization cycle, and with that keyword, you should require fewer cycles overall.

That said, the real issue is with the SCF settings. I’ll need to do a more involved look at the SCF logic to figure out exactly what is happening, but as a solution for your specific case, try removing the keyword basis_guess true. There are some unexpected interactions between the various SCF keywords that are crippling the SCF.

Let me know if that fixes the issue.

I did some more digging. This is related to a known problem in Psi4. I updated our notes on this with the additional information from your report - thanks for that.

Unfortunately, I have no estimate of when this issue will be fixed. Again, let us know if removing basis_guess works for you. I think that will be enough, here.

Thank you for looking into this! Removing basis_guess helped me get the energies associated with several dihedral angles in phenol, but when I adapted my input file to perform a similar calculation with methylaniline, I get another Optimization Convergence Error.

`molecule methylaniline {

0 1
C -0.68557 -0.03323 0.55530
N 0.36336 0.48508 -0.31632
C 0.98155 -0.42438 -1.10319
C 1.93535 0.03826 -2.02051
C 2.62358 -0.82129 -2.87750
C 2.37184 -2.18664 -2.83507
C 1.43786 -2.68530 -1.93673
C 0.75728 -1.81140 -1.08794
H -1.10209 0.76888 1.12805
H -0.27075 -0.76545 1.21611
H -1.45316 -0.48324 -0.03901
H 0.60632 1.43903 -0.33070
H 2.13544 1.04977 -2.06345
H 3.31782 -0.44359 -3.54113
H 2.87630 -2.82704 -3.46786
H 1.24764 -3.69883 -1.89690
H 0.06550 -2.19706 -0.42633
}

set {
basis cc-pvdz
intrafrag_step_limit 0.1
}

set scf {
soscf_max_iter 35
damping_percentage 30.0
maxiter 300
fail_on_maxiter false
}

dihedrals = [i for i in range(0,360, 5)]
PES = []

for phi in dihedrals:
my_string = "1 2 3 4 " + str(phi)
set optking geom_maxiter = 250
set optking ensure_bt_convergence true
set opking opt_coordinates = cartesian
set optking fixed_dihedral = $my_string
print(phi)
try:
E = optimize(‘scf’)
print(E)
PES.append((phi, E))
except SCFConvergenceError:
continue

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]))
`

Here is my output file for the calculation methylaniline.dat and I hope this is a simple fix! Thank you again for your help!

The SCF seems fine this time. There is just something going very, very wrong in the optimizer.

My first thought is that you should fix the typo in set opking opt_coordinates = cartesian. I believe you mean optking. If you look at the output, you’re still optimizing in interal coordinates. I can believe those would struggle for systems with rings.

This makes me glad that we are going to be moving the optimizer logic soon. No guarantees on when, but I am hoping for Psi4 1.5.

After changing this typo to optking opt_coordinates = cartesian, I continue to get an OptimizationConvergenceError and here is my new output file methylaniline_update.out.

If this is a calculation Psi4 cannot perform, can you recommend other programs that could calculate the energies associated with manipulating the dihedrals of ring containing compounds?

The latest input file has some SCF convergence issues, due to the fact that you have SCF damping on (I strongly you recommend you remove that), but the optimizer is acting badly even before that. I’m going to ping @Rollin_King at this point, who can hopefully say something more useful. Even using an initial hessian, this optimization is getting nowhere. Possibly, this is due to the fact that the starting geometry is nowhere close to the starting constraint?

I don’t work with ring-containing compounds, so I can’t give meaningful recommendations.

Setting a fixed dihedral does not construct a geometry with this dihedral. If I see it right, you start out with 180deg and constraint the angle to 0deg! The optimizer then starts to enforce this constraint with an external force that does not care about your molecule. In this case it leads to an out-of-plane bending of a carbon atom that rightfully leads to all sort of problems.(the last geometry in your output).

Either you do this scan in a Z-matrix that allows you define the starting geometry, or your starting geometry for the optimization is close to the constrained value. In the latter case you also want to restart the next optimization from the previous geometry, not from the initial.

Hope that clears some things up.

I’m happy to help, but I can’t look in detail until tomorrow. It may help to clarify the terminology/keywords. A “FROZEN” coordinate remains unchanged. A “FIXED” coordinate is on that the user specifies (fixed by the user input) that may not be satisfied by the input geometry. If you are able to generate for your molecule reasonable guesses (e.g. in the case of a simple dihedral scan) then FROZEN is probably what you want. The capability of fixing a coordinate is primitive (the code as I recall imposes an increasing harmonic potential) but may allow you to push a system toward a more complicated endpoint.

One problem I found was being caused by the discontinuity in the dihedral definition, and the force being applied in the wrong direction - but even after I changed the input to avoid that I still had problems. In particular, converging tightly at each point. I don’t want to further investigate making this work in the C++ optimizer, but I’ll add it to a collection of problems for future testing.

But for what you are doing this is a waste of computer resources anyway. You should construct geometries with the precise desired dihedral angles (which is not part of the ring system, thankfully), and then freeze them in a series of constrained optimizations. You could build your starting geometries in a program (such as Spartan) and export them as xyz. Or, as mentioned previously, you could put your geometry in z-matrix form. It’s ugly (and not good for optimizing rings), but there are online tools to generate a z-matrix, and you only need the z-matrix to specify your input geometry. Then your loop should be over optimizations with frozen coordinates.