optimization blow-up for a sugar molecule

[jhrmnn]

Jordan recently found a geometry optimization problem for sugar molecule (see the attached script). The total energy decreases at the beginning but then it is increased and moves to a strange geometry which causes HF convergence issue. I have tried different basis and convergence tolerance but basically the optimization ends up with the same error. Do you have any idea what the problem is?

from pyscf import gto, scf
from pyscf.geomopt import berny_solver

atomic_coords = '''
C     0.000000000000     0.000000000000     0.000000000000
C     1.452298373280    -0.401861167948    -0.173251347546
C     2.344872240953     0.029889263625     0.983826513995
C     3.252971665673    -0.081516775774    -1.741828636347
C     3.785092679180    -0.364002548175     0.704925081565
C     4.261045334110     0.182591362187    -0.629685589792
H    -0.535619939370    -0.159724413449    -0.936650807056
H    -0.489037279760    -0.591086739000     0.774965311405
H    -0.083341235870     1.058092475028     0.261895660275
H     1.192849869220    -0.212083373080     2.542522206215
H     1.514208266430    -1.493211501840    -0.273323426729
H     2.302321816343     1.125838759480     1.071873981885
H     2.589264335408    -1.663816507960    -2.582233475835
H     3.474172941150     0.600176993441    -2.571962584455
H     3.837063092580    -1.460264413730     0.704205972685
H     4.371981733230    -0.168137321331     2.536015514695
H     4.828078993830     1.694125015840     0.348005496195
H     5.197981129350    -0.331518876442    -0.889920929506
O     1.923676623248     0.212550105645    -1.369362313082
O     2.005075417661    -0.591027362977     2.206647133165
O     3.399618207531    -1.410549203160    -2.138110953605
O     4.470095013330     1.563849680340    -0.536986586354
O     4.633997693990     0.188396169539     1.684719646075
'''

#basis_set = "631+g*"
basis_set = "sto3g"
#spin = 0
#charge = 0

mol = gto.M(atom=atomic_coords, basis=basis_set, spin=spin, charge=charge, verbose=4)

mf  = scf.RHF(mol)
mf.conv_tol = 1e-9
mf.max_cycle = 100

new_mol = berny_solver.optimize(mf, assert_convergence=True)

[mjw99]

Seeing a similar fail on cyclopentene with PySCF (B3LYP/6-31G*) (pyberny ea3038ef), however it is fine with ethene at the same theory/basis.

[jhrmnn]

I've been looking into this problem briefly, and it seems that there is a problem in general with cyclic structures. Once the cycle becomes almost perfectly planar, the structure of the eigenvalues of the Hessian changes, and the optimization blows up. I'm not sure why it happens though. The algorithm is supposed to handle redundant coordinates. I feel that this may be solved by adding additional internal coordinates to the points of planarity, but I'm not sure how yet.

[mjw99]

Thanks for the overview and also thanks for coding this project all up. You have a nice logging system in place; I'll try and look into this to see if I can see anything as well.

[mjw99]

Is there any reason for the difference between eq. 1.10 and https://github.com/azag0/pyberny/blob/6ba5a4cc113b146a3d67217d8b809f7c81fb59ed/berny/berny.py#L140

[jhrmnn]

Ah, that seems to be just a typo in my thesis. The algorithm is taken from Fletcher's book and there it is as (eq 5.1.6):

[mjw99]

Many thanks for the clarification.

[mjw99]

OK, I think I have localised this. There is an issue in the general angle Dihedral differentiating code starting here: https://github.com/azag0/pyberny/blob/6ba5a4cc113b146a3d67217d8b809f7c81fb59ed/berny/coords.py#L166 If one removes this and the conditioning statements and the just uses the code from here https://github.com/azag0/pyberny/blob/6ba5a4cc113b146a3d67217d8b809f7c81fb59ed/berny/coords.py#L145 for differentiating, I can optimise the example sugar above, as well as other examples I have locally that were not working.

[jhrmnn]

Amazing! This makes sense given that planar molecules are problematic. I guess the general-angle code itself is ok, so are the thresholds of 1e-6 incorrect?

[mjw99]

I did have a little play with that, but it did not make much difference. I think it's something specific to the general angle case. But, by all means, do have a play with this aspect.

[mjw99]

I also suspect there may be a problem here: https://github.com/azag0/pyberny/blob/6ba5a4cc113b146a3d67217d8b809f7c81fb59ed/berny/geomlib.py#L109

Should it be -1?

[jhrmnn]

I don't think so, this comes from eq 6 in Swart06.

[jhrmnn]

If one removes this and the conditioning statements and the just uses the code from here

I've looked into this too, and I don't think it's that simple. By using only the phi = pi limit for calculating the derivatives, I suspect that the dihedral angles simply do not change at all, and so the molecule never gets to the particular geometry where it blows up. So rather than pointing to the problem, it rather sidesteps it.

[mjw99]

So rather than pointing to the problem, it rather sidesteps it.

Agreed.

[mjw99]

OK, I think I have it. The following: https://github.com/azag0/pyberny/blob/6ba5a4cc113b146a3d67217d8b809f7c81fb59ed/berny/coords.py#L129

should be: v2 = (coords[self.l]-coords[self.k])/bohr

With the above sugar example, and only this change, the following is seen:

(pyscf) mw529@rabbit:~/code/pyscf/examples/opt$ python sugar.py | grep "converged SCF energy" 
converged SCF energy = -600.63274719127
converged SCF energy = -600.648966473777
converged SCF energy = -600.650860701134
converged SCF energy = -600.651526777412
converged SCF energy = -600.651688496354
converged SCF energy = -600.651714384974
converged SCF energy = -600.651721501987
converged SCF energy = -600.65172298699
converged SCF energy = -600.651723317869

[jhrmnn]

Yep! Also fixes #6. Great catch, I'm surprised it causes problems only for planar molecules.

Feel free to create a pull request, I'll be happy to merge it.

[jhrmnn]

Ok, I think I understand. The difference is basically between measuring the dihedral angle clockwise or anticlockwise, so it doesn't matter as long as the formulas for general and ~0, ~pi angles are consistent, and they weren't, they differed in this convention. So once any dihedral swapped from the general to the limiting regime, the angle jumped by 2*pi (or perhaps pi, not sure), which blowed up the algorithm.

At least I think that's what was happening.