SCF convergence at pathological geometries

Hello everyone!
I have been trying some things in PSI4 recently, and I am interested if you have any advice on making pathological geometries converge. I am trying some new ways of fitting global PESs, and ran into the issue that high energy regions need to be “anchored” by having a few very high energy points with more or less pathological geometries, like this:

memory 8 GB

set basis aug-cc-pVDZ
set scf_type pk
set maxiter 200
set DAMPING_PERCENTAGE 20

set guess sad
set ints_tolerance 1.0E-10
set freeze_core True

molecule {
-1 1
H -3.7800000000 3.1400000000 4.4000000000
H -3.1900000000 -3.9900000000 3.5500000000
H 3.4200000000 2.2500000000 3.3500000000
C -0.2900000000 0.6500000000 1.8700000000
F 0.6100000000 -2.4900000000 -7.2100000000
Cl -0.2400000000 0.8500000000 3.4800000000
}
energy(‘OMP2’)

The thing I noticed, was that this geometry produced a rather distinct failure of the SCF convergence:

@RHF iter 194: -596.76208061116688 -1.13687e-13 1.01356e-02 DIIS/DAMP
@RHF iter 195: -596.76208061116665 2.27374e-13 1.01356e-02 DIIS/DAMP
@RHF iter 196: -596.76208061116677 -1.13687e-13 1.01356e-02 DIIS/DAMP
@RHF iter 197: -596.76208061116665 1.13687e-13 1.01356e-02 DIIS/DAMP
@RHF iter 198: -596.76208061116643 2.27374e-13 1.01356e-02 DIIS/DAMP
@RHF iter 199: -596.76208061116688 -4.54747e-13 1.01356e-02 DIIS/DAMP
@RHF iter 200: -596.76208061116688 0.00000e+00 1.01356e-02 DIIS/DAMP

As you can see the energy seems to have converged to a rather good accuracy, but the density is nowhere near convergence. I have tried to play a bit with damping, MOM and DIIS options, but I could not find anything that can make this converge.
I am wondering if any of you are familiar with this type of SCF non-convergence, and have some tips on solving it. I am well aware that SCF convergence can be a very hard problem, and this geometry is likely multireference in character, and I can work around this issue by some extra work in my input generator, but I am still interested what would it take to converge this.
PS: I also tried this with Molpro, and that implementation also fails at this geometry.

Ive seen this before, but it usually happens in ROHF where some rotations are not optimized due to how the Fock matrix is formalized.

Does this converge in a smaller basis?

If yes, then we can probably project up. If no, then its likely something weird in the gradient and will need some investigating.

I have just tried with STO-3G, 6-31G, 6-31G* and 6-311G, none of them converges.

Using Molpro I could converge a calculation with their small basis set called MINAO, but something went horribly wrong with that, the results were clearly unphysical, for example the HOMO-LUMO gap was negative, and after casting up, the second SCF did not converge. Probably the occupations were completely wrong.

This may be a terrible idea, but the cation (with 2 fewer electrons) should be easier to converge. If so, that could be used as a guess. Another possible guess for small HOMO-LUMO gaps is to use the triplet ROHF anion.

The cation does not converge either, nor does the triplet ROHF anion. (tried both with AVDZ and 6-31G)

Also, on this system SOSCF is completely useless. I tried using it, but whenever the density RMS got low enough for it to start doing anything, it never managed to converge the microiterations. When I gave it 100 microiterations, the SOSCF just diverged to infinity.

Yea, its never close enough to convergence for full SOSCF to kick in and approximate updates are not really helping. Looking at the values of the gradient it appears that the pathological errors are in the coupling of the carbon and chlorine atoms. Not entirely sure why they are not mixing in the diagonalization procedure.

Out of curiosity, what is this for? Plotting it out shows a pretty crazy geometry.

@andysim converging an cation is a valid strategy many people forget about.

Using the GWH guess helps with the problematic density.
With SOSCF and no damping it seems to go somewhere, but even converging the density to 1e-6 is tough.
I also turned off the integral screening, it often helps (“set ints_tolerance 1.0E-16”).

But this seems more of a spin-state problem with all those atoms apart.

The septett converges without any problems. If this physical I am very unsure :wink:
E.g. if all isolated fragments are up spin and the C-Cl is a doublet, it could be a septett… maybe…??

I am working on an automated PES builder tool, that can start from a small number of energies at different geometries, fit a PES, test it for unphysical or inaccurate regions and generate input files with geometries that have a good chance of improving the fitted PES. Unphysical regions on the fitted PES usually arise from not having enough energy points in the right places, and of course early in the PES building there are a lot of regions that are completely non-anchored, and exhibit nasty things like barrierless dissociation to free atoms.

The input I have shown here was generated by my tool, as a candidate for anchoring such an unphysical region. I could try implementing some constraints on what sort of geometry can be generated, but that can get a bit messy, and would also ruin the advantage of being able to use this tool with zero a priori assumptions about the reactions that can or cannot take place.

@hokru Thanks for the tips about guess and the integral screening.
About the spins, since this should be a singulet PES I should use singulet points. However this does give me an idea, I could flag these points as “dirty” and remove them later. After I have enough other points these very high energy anchors will no longer be needed.

Thats what I would do in your situation. The points are probably pretty bad in absolute terms, but the energy is likely high enough that in relative terms its fine. You can always check the deviation between your interpolation and actual result to confirm this.