SCF Bug

[jacharrym]

I found that the SCF is failing to get the correct minima for PsH_2. However, Lowdin is able to get a lower energy after restarting the calculation multiple times.

This is the input PsH2.test.lowdin

GEOMETRY
	e-(H)   SHARON-E-6S2P 	0.00	0.00 	0.00 addParticles=2
	e-(H)   SHARON-E-6S2P 	0.00	0.00 	3.20 
	H	dirac		0.00	0.00 	0.00
	H	dirac		0.00	0.00 	3.20
	e+A	SHARON-E+6S2P	0.00	0.00	0.00 addParticles=-1
	e+A	SHARON-E+6S2P	0.00	0.00	3.20 
	e+B	SHARON-E+6S2P	0.00	0.00	0.00 addParticles=-1
	e+B	SHARON-E+6S2P	0.00	0.00	3.20
END GEOMETRY

TASKS
	method = "UHF"
END TASKS

CONTROL
	readCoefficients=VALUE
END CONTROL

I ran that input 100 times restarting from scratch and then 100 times more restarting from a previous calculation.

These are the results for Total HF energy vs step ID using openLowdin with branch TCP, commit eaaface

However, this is not new... with lowdin2 I got the same but at a different number of restarts

I attached the input and scripts here test-scf.zip

I found the problem because I was creating a test for dipole moment under an external electric field. The effect is the same but the dipole changes are even larger

[fsmoncadaa]

I don't know if this is a bug or a feature.

I've run the test and obtained a similar result.
There are two SCF solutions, one with symmetric e+ orbitals (higher in energy), one with asymmetric e+ orbitals (lower in energy).

Occupied Eigenvectors for: E+A
Symmetric Asymmetric
1 E+A S -0.002872 -0.001070
2 E+A S -0.016874 -0.005798
3 E+A S -0.037588 -0.014229
4 E+A S 0.055513 0.017982
5 E+A S 0.252795 0.037020
6 E+A S 0.301577 0.025365
7 E+A Px 0.000000 0.000000
8 E+A Py 0.000000 0.000000
9 E+A Pz -0.003540 -0.000769
10 E+A Px 0.000000 0.000000
11 E+A Py 0.000000 0.000000
12 E+A Pz 0.075057 0.022701
13 E+A S -0.002872 -0.004386
14 E+A S -0.016874 -0.026765
15 E+A S -0.037588 -0.055910
16 E+A S 0.055513 0.107124
17 E+A S 0.252795 0.504433
18 E+A S 0.301577 0.436373
19 E+A Px 0.000000 0.000000
20 E+A Py 0.000000 0.000000
21 E+A Pz 0.003540 0.000885
22 E+A Px 0.000000 0.000000
23 E+A Py 0.000000 0.000000
24 E+A Pz -0.075057 -0.008407

Depending on the guess, it can converge to either solution. The guess from the HCORE procedure is symmetric, and it seems that breaking this symmetry is not easy. When it reads the converged orbital, it does a two-step SCF, in which it updates the orbitals a little bit. Around the trial 30, it introduces enough asymmetry, probably from round-off errors, to slightly shift the energy, although it is still within the convergence threshold. At trial 40, it converges to the lower energy asymmetric solution, as observed in this sample from the read=true log:

index energy
25 -1.323403924781
26 -1.323403924782
27 -1.323403924782
28 -1.323403924783
29 -1.323403924787
30 -1.323403924794
31 -1.323403924808
32 -1.323403924837
33 -1.323403925041
34 -1.323403925358
35 -1.323403927554
36 -1.323403930965
37 -1.323403954673
38 -1.323403994258
39 -1.323404305040
40 -1.330414232372
41 -1.330414234894
42 -1.330414235600
43 -1.330414235789
44 -1.330414235861
45 -1.330414235881
46 -1.330414235887
47 -1.330414235888
48 -1.330414235889
49 -1.330414235889
50 -1.330414235889

At this point, I think a quick solution is to provide ,vec files with the asymmetric guess.

[jacharrym]

It's probably an "unexpected feature"

The asymmetric guess should get the correct solution. By the way, I fixed the density mixed guess in a previous commit. However, I'm a bit surprised that an asymmetric guess was necessary for those distances.

This is what I got when I was doing the dipole test for the same system under an external electric field

Here the dipole changes even after the first restart, although the behaviour in the energy is the same as the system without the field.

Which is a bit weird because it's getting the same energy with different MO

[jacharrym]

Note to myself @jacharrym Reading the HF orbitals is changing the LUMO coefficients (maybe some kind of rotation?)