DKH does not play nicely with restart

[loriab]

http://forum.psicode.org/t/scf-problems-please-help/765/14

First ROHF (DZ):
Computing 2-order Douglas-Kroll-Hess integrals.
Adding Douglas-Kroll-Hess corrections to the potential integrals.
…
@ROHF Final Energy: -1271.85284181863290

Second ROHF (pwCVTZ):
Computing 2-order Douglas-Kroll-Hess integrals.
Adding Douglas-Kroll-Hess corrections to the potential integrals.
…
SCF Guess: Orbitals guess was supplied from a previous computation.
@ROHF Final Energy: -1267.93014988281766

Even though converged docc and socc are the same, we get a difference of 4 Hartrees in the initial versus final ROHF.

Below is a simpler input that shows the problem.

molecule {
0 4
Fe
H 1 1.566665
}

set {
reference rohf
docc [7, 0, 3, 2]
socc [1, 1, 0, 1]
scf_type pk
relativistic dkh
dkh_order 2
print_mos true
maxiter=500
freeze_core false
}

basis mine2{
spherical
****
Fe 0
S 20 1.00
4.316265E+06 8.048803E-06
6.463424E+05 6.258306E-05
1.470897E+05 3.290239E-04
4.166152E+04 1.387355E-03
1.359077E+04 5.023256E-03
4.905750E+03 1.610140E-02
1.912746E+03 4.590034E-02
7.926043E+02 1.136154E-01
3.448065E+02 2.283869E-01
1.558999E+02 3.221159E-01
7.223091E+01 2.383661E-01
3.272506E+01 7.404667E-02
1.566762E+01 9.214197E-02
7.503483E+00 9.339790E-02
3.312223E+00 1.573965E-02
1.558471E+00 -4.186682E-04
6.839140E-01 5.376318E-05
1.467570E-01 -3.816654E-05
7.058300E-02 4.319603E-05
3.144900E-02 -3.401019E-06
S 20 1.00
4.316265E+06 -4.155954E-06
6.463424E+05 -3.231401E-05
1.470897E+05 -1.699525E-04
4.166152E+04 -7.171369E-04
1.359077E+04 -2.603625E-03
4.905750E+03 -8.399109E-03
1.912746E+03 -2.434109E-02
7.926043E+02 -6.251948E-02
3.448065E+02 -1.365929E-01
1.558999E+02 -2.312707E-01
7.223091E+01 -2.383734E-01
3.272506E+01 3.123837E-02
1.566762E+01 5.086818E-01
7.503483E+00 4.987695E-01
3.312223E+00 9.033552E-02
1.558471E+00 -6.005337E-03
6.839140E-01 2.312454E-04
1.467570E-01 -5.643680E-04
7.058300E-02 4.992260E-04
3.144900E-02 -1.015293E-04
S 20 1.00
4.316265E+06 9.532178E-07
6.463424E+05 7.414605E-06
1.470897E+05 3.898393E-05
4.166152E+04 1.647152E-04
1.359077E+04 5.985980E-04
4.905750E+03 1.942390E-03
1.912746E+03 5.687237E-03
7.926043E+02 1.501329E-02
3.448065E+02 3.452455E-02
1.558999E+02 6.495820E-02
7.223091E+01 7.716194E-02
3.272506E+01 -1.873411E-02
1.566762E+01 -3.009185E-01
7.503483E+00 -4.554661E-01
3.312223E+00 1.286463E-01
1.558471E+00 7.183316E-01
6.839140E-01 4.051743E-01
1.467570E-01 2.168227E-02
7.058300E-02 -8.343566E-03
3.144900E-02 3.658979E-03
S 20 1.00
4.316265E+06 -2.063008E-07
6.463424E+05 -1.604169E-06
1.470897E+05 -8.438437E-06
4.166152E+04 -3.563151E-05
1.359077E+04 -1.295998E-04
4.905750E+03 -4.201534E-04
1.912746E+03 -1.231954E-03
7.926043E+02 -3.248922E-03
3.448065E+02 -7.493717E-03
1.558999E+02 -1.410102E-02
7.223091E+01 -1.691600E-02
3.272506E+01 4.218996E-03
1.566762E+01 6.833810E-02
7.503483E+00 1.098201E-01
3.312223E+00 -4.009005E-02
1.558471E+00 -2.174739E-01
6.839140E-01 -2.465135E-01
1.467570E-01 2.731435E-01
7.058300E-02 5.748321E-01
3.144900E-02 3.012713E-01
S 20 1.00
4.316265E+06 -4.009367E-07
6.463424E+05 -3.189255E-06
1.470897E+05 -1.623079E-05
4.166152E+04 -7.157920E-05
1.359077E+04 -2.463958E-04
4.905750E+03 -8.544907E-04
1.912746E+03 -2.307593E-03
7.926043E+02 -6.728292E-03
3.448065E+02 -1.366165E-02
1.558999E+02 -3.062240E-02
7.223091E+01 -2.631137E-02
3.272506E+01 -9.760183E-03
1.566762E+01 1.801906E-01
7.503483E+00 1.529634E-01
3.312223E+00 5.505413E-02
1.558471E+00 -9.551364E-01
6.839140E-01 2.586813E-01
1.467570E-01 1.834049E+00
7.058300E-02 -9.333240E-01
3.144900E-02 -6.981605E-01
S 1 1.00
3.144900E-02 1.000000E+00
P 16 1.00
1.774569E+04 4.100000E-05
4.200721E+03 3.690000E-04
1.364429E+03 2.129000E-03
5.220806E+02 9.369000E-03
2.214595E+02 3.309700E-02
1.009096E+02 9.443100E-02
4.840115E+01 2.080770E-01
2.398536E+01 3.323330E-01
1.218250E+01 3.329870E-01
6.242298E+00 1.568430E-01
3.110944E+00 2.154900E-02
1.509958E+00 -2.095000E-03
7.108450E-01 -1.739000E-03
2.731900E-01 -3.000000E-04
1.042330E-01 2.900000E-05
3.829100E-02 -1.100000E-05
P 16 1.00
1.774569E+04 -1.500000E-05
4.200721E+03 -1.300000E-04
1.364429E+03 -7.510000E-04
5.220806E+02 -3.329000E-03
2.214595E+02 -1.191200E-02
1.009096E+02 -3.493300E-02
4.840115E+01 -7.998900E-02
2.398536E+01 -1.346360E-01
1.218250E+01 -1.385980E-01
6.242298E+00 3.027800E-02
3.110944E+00 3.332160E-01
1.509958E+00 4.561530E-01
7.108450E-01 2.850510E-01
2.731900E-01 4.614400E-02
1.042330E-01 -3.249000E-03
3.829100E-02 1.357000E-03
P 16 1.00
1.774569E+04 3.000000E-06
4.200721E+03 2.900000E-05
1.364429E+03 1.650000E-04
5.220806E+02 7.340000E-04
2.214595E+02 2.626000E-03
1.009096E+02 7.725000E-03
4.840115E+01 1.773300E-02
2.398536E+01 3.005500E-02
1.218250E+01 3.109400E-02
6.242298E+00 -1.004800E-02
3.110944E+00 -8.830600E-02
1.509958E+00 -1.298240E-01
7.108450E-01 -7.693700E-02
2.731900E-01 2.126610E-01
1.042330E-01 5.730610E-01
3.829100E-02 3.696510E-01
P 16 1.00
1.774569E+04 5.000000E-06
4.200721E+03 4.200000E-05
1.364429E+03 2.410000E-04
5.220806E+02 1.085000E-03
2.214595E+02 3.831000E-03
1.009096E+02 1.142300E-02
4.840115E+01 2.579200E-02
2.398536E+01 4.481800E-02
1.218250E+01 4.459800E-02
6.242298E+00 -1.117700E-02
3.110944E+00 -1.381340E-01
1.509958E+00 -1.882850E-01
7.108450E-01 -1.073990E-01
2.731900E-01 4.448630E-01
1.042330E-01 6.402390E-01
3.829100E-02 6.445700E-02
P 1 1.00
3.829100E-02 1.000000E+00
D 8 1.00
1.133440E+02 3.530000E-03
3.364140E+01 2.578400E-02
1.233100E+01 9.911900E-02
4.994780E+00 2.390730E-01
2.072800E+00 3.571990E-01
8.307530E-01 3.621880E-01
3.091780E-01 2.364610E-01
1.001300E-01 6.011800E-02
D 8 1.00
1.133440E+02 -3.890000E-03
3.364140E+01 -2.844200E-02
1.233100E+01 -1.124290E-01
4.994780E+00 -2.742570E-01
2.072800E+00 -3.155460E-01
8.307530E-01 5.710900E-02
3.091780E-01 5.636040E-01
1.001300E-01 3.846370E-01
D 1 1.00
1.001300E-01 1.000000E+00
F 2 1.00
3.224300E+00 4.222490E-01
7.758000E-01 7.714680E-01
****
H 0
S 3 1.00
13.0100000 0.0196850
1.9620000 0.1379770
0.4446000 0.4781480
S 1 1.00
0.1220000 1.0000000
P 1 1.00
0.7270000 1.0000000
****
}

set basis mine2
set basis_relativistic mine2

set scf df_scf_guess false
set scf guess gwh
set ccenergy r_convergence 10
energy('scf')

clean()

basis{
assign Fe cc-pwcvtz-dk
assign H cc-pvtz-dk
}

set scf guess read
energy('cc2')

[jturney]

I looked into this yesterday. On my laptop when I disable DKH, I still obtain a large total energy difference:

  @ROHF Final Energy: -1262.91530042362160
  @ROHF Final Energy: -1259.00365041946634

I emailed one of the issue reporters directly and waiting to hear a response. At least with my compilation (master/debug) it doesn't appear to be a DKH issue.

[andysim]

I'm guessing the second basis is larger than the first. If so, isn't this just a case of converging on the wrong state? Stability analysis could confirm.

[jturney]

Yes, that's it. Negative eigenvalue in the stability analysis.

[natedey]

hey all! I don't think at all think this is fundamentally an issue with the SCF converging to the wrong state, even if the stability analysis is showing a negative eigenvalue. IMO, there could be an issue with the basis set projection.

For this electronic state in this system, it is indeed difficult to converge to the appropriate state. However, the multireference character results in low-lying LUMOs. So if PSI4 was populating the wrong orbitals, we would see the energy of the larger basis set ROHF increasing by (this is all back-of-the-envelope, btw) 5 mEh. We know that the smaller basis set (cc-pwCVDZ-DK, which is explicitly defined in the input file because that basis set is not in the PSI4 library for iron) is converging to the right state.

If DZ->TZ projection is going wrong by 3 Hartrees, then there is something totally haywire with PSI4. Molpro's SCF guess from atomic densities sucks, but I have never seen it crap the bed with incorrect states on the order of Hartrees!

Try this test: run the ROHF energy calculation using just the TZ set, with no DZ step and no "set scf guess read"
When I do this, I get
@rohf Final Energy: -1271.65162910275285

HERE is where PSI4 has the incorrect state. I haven't run an orbital stability analysis, but I can tell by looking at the MO coefficients (thanks Yukio!) --- 10A1 (doubly occupied d2+) needs to be swapped with 14A1 (singly occupied d0) to change the excited Phi state to the correct Pi state.

Another clue that leads me to believe there's an actual bug somewhere: look at the SCF orbital energies for the initial DZ guess versus the screwed up TZ projection. The only significant difference is
DZ ROHF:
1A1 -263.972182
and
TZ ROHF:
1A1 -261.966897

Since 1A1 is doubly occupied, there's your 4 Hartrees worth of error.
Looking at the MO coefficients
DZ ROHF:
FE1 s0 0.9825555 (duh, first orbital should be iron 1s)

TZ ROHF:
1 FE1 s0 0.8531471 0.3343274 0.0000000 0.0000000 0.0003960
2 FE1 s0 -0.5169916 1.4108806 0.0000000 0.0000000 0.0017387
All sorts of crazy mixing of 1s/2s iron basis functions in the first two MOs.

To me, something messed up with core s orbitals in the SCF indicates a problem with the basis set projection, a problem with the DKH, or perhaps even an index start error if only the first MO is wrong.
Fun!

[natedey]

Justin, can you send me the input where you disabled DKH? An SCF read of converged orbitals from the same molecule should never provide such a bad guess that you get an energy error on the order of Hartrees. Even if it's a naughty molecule like FeH. This points to an error in the way that PSI4 is reading orbitals for the set scf guess read procedure.

[jturney]

It was in the email that I sent you early last week. I've attached here, too.
input.txt

[susilehtola]

I had a look at the basis projection code in psi4/src/psi4/libmints/wavefunction.cc, and I don't really understand why it is so complicated. Namely, the whole point about the basis set projection is that you use the resolution of the identity

 \sum_{UV} |U> <U|V>^{-1} <V|

to project between basis sets. That is, given orbitals in an original basis

  c_{ui} |u>

you get the new orbitals as

   c_{ui} |U> <U|V>^{-1} <V|u>
= [<U|V>^{-1} <V|u> c_{ui} ]  |U>

that is, a simple matrix multiply

   CA = [X_A^2 S_{AB} C_B ]

where X_A is the orthogonalizing matrix in the new basis, S_{AB} is the overlap between the old and new bases, and C_B are the coefficients in the old basis.

Now, the real question is just what you have to do when the dimension of the basis changes. When the new basis is smaller than the old one, there's no problem. If it's bigger than the old one, you just need to run a SVD on CSX and pad these to the end of the MO coefficient matrix.