Convergence issues for DFT

1

I am trying to compare PSI4 with couple other DFT codes and found a strange problem with SCF convergence.
I was able to simplify the non-working example to the LiF: If I take LiF and try to compute DFT energy with PBE, when the distance is sufficiently large the SCF fails. And it looks exactly like it does look for more complex cases. Sometimes similar behavior I see when try to make the geometry optimization, even near the optimal geometry.

Example (I am using python API, this is fresh 1.3 version):

import psi4
mol = """
   Li
   F 1 4.0
"""
geom = psi4.geometry(mol)
e = psi4.energy("PBE/aug-cc-pvqz",molecule=geom)

The cycle ends with

@DF-RKS iter 96: -105.93407344771342 3.05420e-10 1.26750e-02 DIIS
@DF-RKS iter 97: -105.93407344818426 -4.70848e-10 1.26750e-02 DIIS
@DF-RKS iter 98: -105.93407344775187 4.32394e-10 1.26750e-02 DIIS
@DF-RKS iter 99: -105.93407344819431 -4.42441e-10 1.26750e-02 DIIS
@DF-RKS iter 100: -105.93407344820115 -6.83542e-12 1.26750e-02 DIIS

which is obviously oscillation around non-ground state because somewhere soon after beginning it was around correct value

@DF-RKS iter 1: -107.00420082761590 -1.40277e-01 8.58256e-03 DIIS
@DF-RKS iter 2: -94.40386218333836 1.26003e+01 2.49009e-02 DIIS
@DF-RKS iter 3: -106.99605273796766 -1.25922e+01 8.58719e-03 DIIS
Occupation by irrep:
A1 A2 B1 B2
DOCC [ 5, 0, 1, 0 ]
@DF-RKS iter 4: -106.99127569765074 4.77704e-03 8.70503e-03 DIIS
Occupation by irrep:
A1 A2 B1 B2
DOCC [ 4, 0, 1, 1 ]
@DF-RKS iter 5: -105.54067118242492 1.45060e+00 1.68138e-02 DIIS
Occupation by irrep:
A1 A2 B1 B2
DOCC [ 5, 0, 0, 1 ]
@DF-RKS iter 6: -107.20369889781466 -1.66303e+00 8.57220e-04 DIIS
Occupation by irrep:
A1 A2 B1 B2
DOCC [ 5, 0, 1, 0 ]
@DF-RKS iter 7: -105.73636336480513 1.46734e+00 1.47178e-02 DIIS

Surprisingly, triplet state does converge. However, according to the result, the ground state should be singlet.

Could you explain, where this problem comes from and is it possible to tune some parameters to avoid this? I tried to play myself, but with no success.

2

Triplett instabilities occasional appear for bond breaking PES.

You have two degenerate orbitals and the occupation is jumping in-between them. You can try the SOSCF procedure (soscf true) to help convergence of the density. Worked for me.

3

It worked. I should double check more complex case when I did try the same and it did not work: maybe I messed with some other options.

Is there any reference where this kind of instabilities is explained in more details?

4

triplet instabilities for RHF are in many textbooks and of course papers. I don’t have references at hand.
“Molecular Electronic-Structure Theory” from Helgaker/Jorgensen/Olsen discusses it somewhere.

5

Oh, I realized why I had the problem even with SOSCF. I did iterations over the distance. And the density at the next distance seem to be taken from previous one. Then the SCF is stuck. If I start from fresh SAD guess, the algorithm does converge. Still has problems when sometimes it falls back to DIIS, but definitely convergence is much better.

Forgive me my naive question, but if this is a RHF problem, why setting UHF does not help?

6

It helps a lot with the SCF convergence to make PES scans starting from the long distance and getting smaller if you reuse the orbitals.

UHF removes the spatial restrictions of RHF, but the SCF may not always be able to escape them, e.g. when you have a very stable RHF guess.

In any case it is instructive to plot both the RHF and UHF curves and compare. Also with different methods. One can learn a lot :slight_smile:

7

Thank you for the suggestions.

Indeed, I am trying to plot PES, but the main idea was to use LiF as a minimal example to understand how to proceed when there are convergence problems.

Correct me if I am wrong, but this particular problem with LiF does not seem to be a real problem on a chemistry/physics side: the equations are complicated and non-linear, but they do have a solution within the given method/basis. Which SCF is unable to find without a lot of extra work to force it into correct solution.