Spin state changes during calculation

[susilehtola]

Testing some calculations, I ran into very odd behavior in Psi4 1.3.2. The input

molecule {
0 2
Y
}

set basis seg-cc-pv5z-pp
set scf_type direct
set df_scf_guess false
set reference uhf
energy('scf')

leads to the output

  ==> Iterations <==

                        Total Energy        Delta E     RMS |[F,P]|

    Occupation by irrep:
             Ag   B1g   B2g   B3g    Au   B1u   B2u   B3u 
    DOCC [     2,    0,    0,    0,    0,    1,    1,    1 ]
    SOCC [     0,    0,    0,    0,    0,    1,    0,    0 ]

   @UHF iter   1:   -33.78761858648352   -3.37876e+01   4.53671e-02 DIIS
    Occupation by irrep:
             Ag   B1g   B2g   B3g    Au   B1u   B2u   B3u 
    DOCC [     1,    0,    0,    0,    0,    1,    1,    1 ]
    SOCC [     1,    1,    0,    1,    0,    0,    0,    0 ]

   @UHF iter   2:   -36.39804695906182   -2.61043e+00   2.16710e-02 DIIS

For some reason, the occupation update changes the spin state from a doublet (one unpaired electron) to hextet (3 unpaired electrons). Basis set is attached

seg-cc-pv5z-pp.gbs.txt

[susilehtola]

Reproduced with current master with set scf_initial_accelerator none.

[susilehtola]

On further thought, I am not sure this is an error. In UHF there's no guarantee that the orbitals are spin-paired; the beta spatial orbitals are allowed to differ from the alpha ones. This naturally leads to spin contamination, but that's a feature not a bug!

[JonathonMisiewicz]

You said it's a doublet in the molecule input. Psi should honor that. This is a bug.

[susilehtola]

You're right. The problem is the use of docc and socc: https://github.com/psi4/psi4/blob/1ffb67bf52dcca70e357d23714c7d975f9d1356e/psi4/src/psi4/libscf_solver/hf.cc#L474

After the recent addition of the sanity check in #2488, the code crashes in the iteration where some symmetries have more beta electrons and alpha electrons.

Occupations
sym 0: 1 alpha 2 beta
sym 1: 0 alpha 0 beta
sym 2: 1 alpha 0 beta
sym 3: 1 alpha 0 beta
sym 4: 0 alpha 0 beta
sym 5: 1 alpha 1 beta
sym 6: 1 alpha 1 beta
sym 7: 1 alpha 1 beta

The docc and socc syntax is unable to handle this, which leads to spin flips and the wrong spin state.

[susilehtola]

The alpha and beta occupations need to be specified separately for spin-unrestricted theories.

[JonathonMisiewicz]

When I attempt this, L2 complains that the angular momentum limit exceeded.

Is this expected? I should be able to fix this with a higher AM L2, but I'd rather avoid building myself, and we certainly aren't having the test suite build L2 itself.