Optimisation failure with flat potential energy surface near minimum

[mjw99]

The following example, with PySCF 1.5.3 and Pyberny 0.4.0:

# c1(F)c(F)c(F)c(F)c(F)c1C#CCN(C)(C(=O)c1ccccc1)
# ZCUWMRDRNFJWEK-UHFFFAOYSA

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

atomic_coords = '''
C      7.403569010821  -0.674065308632   1.565827066815
C      4.823525932952  -0.157981104014   1.717950019842
C      3.134110777591  -1.439026443856   0.153256788702
C      4.032675549822  -3.206109342937  -1.598708301382
C      2.217593607177  -4.724126338800  -3.094237556363
O      2.354787723821  -7.035828306981  -2.910745149668
N      0.369252484740  -3.497694083957  -4.511343177174
C      0.568618590882  -0.907446485016  -5.440143567398
C     -1.109836152957   0.826188261660  -4.080485620773
C     -2.431510604478   2.256144020118  -2.912634875792
C     -3.953496025203   3.916646365774  -1.422963771797
C     -4.364700429908   3.393570174494   1.130434167715
C     -5.838686807069   4.995679982900   2.603475681813
C     -6.911295355372   7.139763243832   1.540504736745
C     -6.508783690839   7.684760258156  -0.994184914134
C     -5.036120121966   6.081138668850  -2.471572798319
C     -1.902765234825  -4.911587170357  -5.125882112883
C      6.614041435978  -3.751673275099  -1.724753033891
C      8.298921248640  -2.471950743544  -0.146453774654
H      4.120736787227   1.197141499912   3.092536802851
H      1.117584030068  -1.117584030068   0.367362758615
H      0.095431169291  -0.929556280674  -7.457426205371
H      2.528075609443  -0.258325561228  -5.311642190927
F     -3.327807705359   1.337359178355   2.193027167558
F     -6.216632031982   4.468635366759   5.051237930962
F     -8.329345839245   8.683858460214   2.960822891969
F     -7.536416757378   9.760057488154  -2.019739281935
F     -4.663088184977   6.643143218296  -4.915744567831
H     -2.876730079425  -4.076895141137  -6.744810483798
H     -1.450175827991  -6.883327408728  -5.548424874335
H      7.303791471444  -5.202793966153  -3.005798373733
H     10.305621420316  -2.908666450931  -0.215806723425
H     -3.148094750913  -4.864344017243  -3.473505589563
H      8.714660996044   0.291017823183   2.820038295688

'''

basis_set = "sto3g"
mol = gto.M(atom=atomic_coords, unit="Bohr", basis=basis_set, verbose=8)

method = dft.RKS(mol)
method.grids.prune = dft.gen_grid.treutler_prune
method.xc = "b3lypg"

opt_mol = berny_solver.optimize(method, verbose=5, assert_convergence=True)

fails with the following error:

  File "/home/mw529/.conda/envs/rdkit/lib/python2.7/site-packages/berny/berny.py", line 144, in send
    raise RuntimeError('The trust radius got too small, check forces?')
RuntimeError: The trust radius got too small, check forces?

Do we need to add another criteria for whether and optimisation has succeeded or not?

[jhrmnn]

Will check. This may or may not be related to the issues that Ask ran into when he was testing pyberny before including it into ASE: https://gitlab.com/ase/ase/merge_requests/889

[mjw99]

Hi Jan, Thanks for the reply and the link; it makes interesting reading. I think it is similar to what I am seeing. I will investigate more.

[mjw99]

I can reproduce the fail with the smaller molecule mentioned in the above link.

# Cyanogen
# N#CC#N
# JMANVNJQNLATNU-UHFFFAOYSA-N

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

atomic_coords = '''
N      3.545830    3.669192    7.228181
C      3.601888    3.624940    6.062501
C      3.671915    3.575700    4.697549
N      3.727670    3.537778    3.532496
'''

basis_set = "sto3g"
#basis_set = "6-31g*"
mol = gto.M(atom=atomic_coords, basis=basis_set, verbose=8)

method = dft.RKS(mol)
method.grids.prune = dft.gen_grid.treutler_prune
method.xc = "b3lypg"

opt_mol = berny_solver.optimize(method, verbose=5, assert_convergence=True)

What is interesting here, is that if I disable dihedral in Berny's defaults{}, the molecule optimises within 5 steps. It seems that selection of the internal coordinate terms is important here. Also, perhaps variants of the angle terms, i.e. linear bend, may be needed in place of linear dihedrals.

[mjw99]

Discussion of linear dihedrals within this context and possible solution on page 19 of this thesis.

[mjw99]

The exact algorithm for the "Linear Bend" is described on p10 of 12 of http://dx.doi.org/10.1007/s00214-016-1847-3 by Schlegel et al. Essentially, it could be implemented during the angle assignment phase, given the criteria on page p9 of 12. However, it will need the concept of a dummy atom.

[jhrmnn]

A useful reference, thanks! Pity they haven't published the referenced Mathematica program. I have created a Github project to track all the potential enhancement from that paper, the linear bends being one of them.

[mjw99]

Whilst this is not a fix for #30, changing the following in: https://github.com/azag0/pyberny/blob/8c07a463e82515953f7959dcc3e33c53a2041d93/berny/coords.py#L387 and: https://github.com/azag0/pyberny/blob/8c07a463e82515953f7959dcc3e33c53a2041d93/berny/coords.py#L390 to:

    nonlinear_l = [
        n for n, ang in zip(neigh_l, angles_l) if ang < pi-1e-1 and ang >= 2e-1
    ]
    nonlinear_r = [
        n for n, ang in zip(neigh_r, angles_r) if ang < pi-1e-1 and ang >= 2e-1
    ]

fixes both the examples in this ticket. It does not break the example in #10 ; hence it might be worth applying and seeing if it helps with the ASE issues.

Edit; pasted incorrect "fix"

[jhrmnn]

Sorry for the late response. So I can confirm that it fixes the issue with cyanogen, with MOPAC as well. But with MOPAC the large molecule still doesn't converge. I've tried also changing the thresholds for linear_l and linear_r, but that didn't help. Looking at the optimization, it seems that the linear dihedral is still an issue.

I see that you picked 0.1 rad for the ~0° angles and 0.2 rad for the ~180° angles. Is the fix sensitive to the exact value and that they are different? I'd try myself, but I have PySCF only on my laptop, and this molecule is way too large.

[mjw99]

Hi Jan, This was quite an ad-hoc attempt for this and these values are really quite empirical; I was also testing this with another larger system as well.

I've just been trying to reproduce the ZCUWMRDRNFJWEK-UHFFFAOYSA example above with my changes on another machine, but it seems not to be working. I'm just going to recheck this on the original machine I ran it on.

[jhrmnn]

I've just pushed a0ca9d909785c8fae338b6a9058c6d47d4e54905 that increases the linear threshold to 5°, and this definitely fixes the cyanogen, but I'm afraid the issue with the large molecule runs deeper than that. Hopefully the linear bends will fix it. If you're not on it, I might find some time to add them over Christmas.

[mjw99]

I think I've found the issue I described in https://github.com/azag0/pyberny/issues/23#issuecomment-445235620 ; I had used a larger basis set as well (i.e. 6-31g*). Apologies for the misreporting.

I also agree that this is a deeper issue and that the linear bends approach will help this. I am hoping to find some more time for this soon.