Looking for global minima X2C module (SCF convergence issue)

Dear developers,

Through days of learning, I have been trying my best to get the global minima of different atoms&molecules that have subtle openshell structures, which involves rotation of orbitals, parsing density guesses from different calculations, decontracting basis sets and so on, to get there. For many times, when I thought I catch sight of the light of the other side of the tunnel, the SCF convergence always beats me.
For example, the 3d transition metal Ti, the groundstate occupation is triplet d^2 s^2, if you do a calculation as follows:

molecule Ti {
0 3
Ti
}

set {
reference rohf
basis cc-pvdz
maxiter 300
die_if_not_converged False
docc [4,0,0,0,0,2,2,2]
socc [0,1,0,1,0,0,0,0]
units angstrom
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}

scf_e, scf_wfn = energy(‘scf’, return_wfn=True)


The SCF just goes:
@ROHF iter 280: -848.16630505675914 -1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 281: -848.16630505675892 2.27374e-13 1.04937e-03 DIIS/MOM
@ROHF iter 282: -848.16630505675903 -1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 283: -848.16630505675903 0.00000e+00 1.04937e-03 DIIS/MOM
@ROHF iter 284: -848.16630505675903 0.00000e+00 1.04937e-03 DIIS/MOM
@ROHF iter 285: -848.16630505675926 -2.27374e-13 1.04937e-03 DIIS/MOM
@ROHF iter 286: -848.16630505675880 4.54747e-13 1.04937e-03 DIIS/MOM
@ROHF iter 287: -848.16630505675903 -2.27374e-13 1.04937e-03 DIIS/MOM
@ROHF iter 288: -848.16630505675892 1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 289: -848.16630505675926 -3.41061e-13 1.04937e-03 DIIS/MOM
@ROHF iter 290: -848.16630505675892 3.41061e-13 1.04937e-03 DIIS/MOM
@ROHF iter 291: -848.16630505675903 -1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 292: -848.16630505675892 1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 293: -848.16630505675926 -3.41061e-13 1.04937e-03 DIIS/MOM
@ROHF iter 294: -848.16630505675892 3.41061e-13 1.04937e-03 DIIS/MOM
@ROHF iter 295: -848.16630505675857 3.41061e-13 1.04937e-03 DIIS/MOM
@ROHF iter 296: -848.16630505675869 -1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 297: -848.16630505675937 -6.82121e-13 1.04937e-03 DIIS/MOM
@ROHF iter 298: -848.16630505675903 3.41061e-13 1.04937e-03 DIIS/MOM
@ROHF iter 299: -848.16630505675914 -1.13687e-13 1.04937e-03 DIIS/MOM
@ROHF iter 300: -848.16630505675926 -1.13687e-13 1.04937e-03 DIIS/MOM

And if you look at the MOs of the density-not-converged state, you will actually see it is stuck in a triplet d^3 s^1 occupation(because two of the d orbitals are in the same symmetry as s, so instead of doubly occupying s, it chooses to occupy d).
This is fine, @dgasmith Daniel cooked something for me to rotate the orbitals, then I proceed to switch the doubly occupied d and single occupied s and did this:

set {
reference rohf
basis cc-pvdz
maxiter 300
die_if_not_converged False
docc [4,0,0,0,0,2,2,2]
socc [0,1,0,1,0,0,0,0]
guess read
units angstrom
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}
scf_e, scf_wfn = energy(‘scf’, return_wfn=True)

The result of this calculation is perfect, the convergence meets the thresholds, the density is optimized, the occupation is the expected triplet d^2 s^2. This is great!
And the next step is to try to feed this MOs as an perfect initial guess for the following calculation:

set {
reference rohf
basis diatomic_awcvqzdk
basis_relativistic diatomic_awcvqzdk_decon
die_if_not_converged False
maxiter 300
docc [4,0,0,0,0,2,2,2]
socc [0,1,0,1,0,0,0,0]
relativistic x2c
frozen_docc [2,0,0,0,0,1,1,1]
units angstrom
print_MOs true
guess read
MOM_START 10
print 2
SCF_TYPE PK
}

This is the first 30 iterations I got:
@ROHF iter 0: -852.71248186439072 -8.52712e+02 0.00000e+00
@ROHF iter 1: -506.88722039301888 3.45825e+02 1.47207e+00
@ROHF iter 2: -849.47743169553962 -3.42590e+02 6.57636e-02 DIIS
@ROHF iter 3: -850.18921825984944 -7.11787e-01 3.50998e-02 DIIS
@ROHF iter 4: -850.91687102350818 -7.27653e-01 1.52076e-02 DIIS
@ROHF iter 5: -851.88208308166327 -9.65212e-01 1.35810e-02 DIIS
@ROHF iter 6: -852.39043556728734 -5.08352e-01 4.90455e-03 DIIS
@ROHF iter 7: -851.85424643806550 5.36189e-01 7.17867e-03 DIIS
@ROHF iter 8: -852.21747060536939 -3.63224e-01 4.60652e-03 DIIS
@ROHF iter 9: -852.33400181298271 -1.16531e-01 3.32486e-03 DIIS
@ROHF iter 10: -852.53026112962118 -1.96259e-01 1.02129e-03 DIIS/MOM
@ROHF iter 11: -852.55830448938229 -2.80434e-02 1.20901e-03 DIIS/MOM
@ROHF iter 12: -852.56398493340032 -5.68044e-03 1.14473e-03 DIIS/MOM
@ROHF iter 13: -852.54352860453218 2.04563e-02 1.75843e-03 DIIS/MOM
@ROHF iter 14: -852.45683197908886 8.66966e-02 2.43253e-03 DIIS/MOM
@ROHF iter 15: -852.43842850666647 1.84035e-02 2.35650e-03 DIIS/MOM
@ROHF iter 16: -852.51899267604995 -8.05642e-02 1.27321e-03 DIIS/MOM
@ROHF iter 17: -852.54265933866191 -2.36667e-02 1.17389e-03 DIIS/MOM
@ROHF iter 18: -852.55935493554387 -1.66956e-02 1.08803e-03 DIIS/MOM
@ROHF iter 19: -852.41305256959436 1.46302e-01 2.88337e-03 DIIS/MOM
@ROHF iter 20: -852.39514794386821 1.79046e-02 3.26967e-03 DIIS/MOM
@ROHF iter 21: -852.43916678618598 -4.40188e-02 2.38863e-03 DIIS/MOM
@ROHF iter 22: -852.46722872063992 -2.80619e-02 1.89799e-03 DIIS/MOM
@ROHF iter 23: -852.48020779595902 -1.29791e-02 2.54558e-03 DIIS/MOM
@ROHF iter 24: -852.37194530347324 1.08262e-01 4.32037e-03 DIIS/MOM
@ROHF iter 25: -852.23982354355758 1.32122e-01 5.27231e-03 DIIS/MOM
@ROHF iter 26: -852.35991488315858 -1.20091e-01 3.75578e-03 DIIS/MOM
@ROHF iter 27: -852.32448253019447 3.54324e-02 4.34466e-03 DIIS/MOM
@ROHF iter 28: -852.35397532458012 -2.94928e-02 3.86708e-03 DIIS/MOM
@ROHF iter 29: -852.35773139595483 -3.75607e-03 3.79735e-03 DIIS/MOM
@ROHF iter 30: -852.32079243450130 3.69390e-02 4.39597e-03 DIIS/MOM

It’s stuck again! And the occupation is going back to triplet d3s1 too!
Here allow me to quote the ROHF-SCF energy I got from MOLPRO:
-852.734292838881
(For most of the closed-shell and trivial open-shell cases, the two values, SCF energies of PSI4 and MOLPRO, agree within 10-E8, which is perfect.)
If you look at the very first iteration, you can see it actually reads in the correct occupation and did one iteration with it:
@ROHF iter 0: -852.71248186439072 -8.52712e+02 0.00000e+00
@ROHF iter 1: -506.88722039301888 3.45825e+02 1.47207e+00
The energy went down to -852.71248186439072, it almost seems like after 10 iterations it is going to converge.
But the next iteration destroys it by bring back the energy to nowhere and start to re-optimize from a wrong occupation again, I finished the this calculation with die_if_not_converged False and maxiter 300, and the MOs are just as wrong as the very first calculation I did.

So in a word, I really want to figure out what PSI4 is doing when it hits second iteration above
@ROHF iter 0: -852.71248186439072 -8.52712e+02 0.00000e+00
@ROHF iter 1: -506.88722039301888 3.45825e+02 1.47207e+00
Why it could not proceed as the energy is already much lower than even in the end of the iteration
@ROHF iter 300: -852.32084708579191 3.69937e-02 4.39524e-03 DIIS/MOM
Is there a hidden option somewhere I don’t know to really decrease the step size to a point that it does not completely mess the “perfect” initial guess? Or is it an algorithmic thing that can’t be tackled within this week?

Thanks for reading through such a long post!
Any help would be appreciated!
-Rulin

I think I can post a calculation which is more generic, that might help:


memory 1050 mb

def orb_rotate(C, h, i, j, theta):

if h > C.nirrep():
raise p4util.PsiException(“ROTATE: Expression %s irrep number is larger than the number of irreps” %
(str(h)))

if max(i, j) > C.coldim()[h]:
raise p4util.PsiException(“ROTATE: Expression %s orbital number exceeds number of orbitals in irrep” %
(str(max(i, j))))

theta = np.deg2rad(theta)

x = C.nph[h][:, i].copy()
y = C.nph[h][:, j].copy()

xp = np.cos(theta) * x - np.sin(theta) * y
yp = np.sin(theta) * x + np.cos(theta) * y

C.nph[h][:, i] = xp
C.nph[h][:, j] = yp

molecule Mn {
0 6
Mn
}

set {
reference rohf
basis cc-pvdz
maxiter 250
docc [4,0,0,0,0,2,2,2]
socc [2,1,1,1,0,0,0,0]
units angstrom
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}

scf_e, scf_wfn = energy(‘scf’, return_wfn=True)
clean()

orb_rotate(scf_wfn.Ca(), 0, 3, 4, 90.0)

filename = core.get_writer_file_prefix(scf_wfn.molecule().name()) + “.180.npz”
data = {}
data.update(scf_wfn.Ca().np_write(None, prefix=“Ca”))
data.update(scf_wfn.Cb().np_write(None, prefix=“Cb”))

Ca_occ = scf_wfn.Ca_subset(“SO”, “OCC”)
data.update(Ca_occ.np_write(None, prefix=“Ca_occ”))

Cb_occ = scf_wfn.Cb_subset(“SO”, “OCC”)
data.update(Cb_occ.np_write(None, prefix=“Cb_occ”))

data[“reference”] = core.get_option(‘SCF’, ‘REFERENCE’)
data[“nsoccpi”] = scf_wfn.soccpi().to_tuple()
data[“ndoccpi”] = scf_wfn.doccpi().to_tuple()
data[“nalphapi”] = scf_wfn.nalphapi().to_tuple()
data[“nbetapi”] = scf_wfn.nbetapi().to_tuple()
data[“symmetry”] = scf_wfn.molecule().schoenflies_symbol()
data[“BasisSet”] = scf_wfn.basisset().name()
data[“BasisSet PUREAM”] = scf_wfn.basisset().has_puream()
np.savez(filename, **data)

set {
reference rohf
basis cc-pvdz
maxiter 250
docc [4,0,0,0,0,2,2,2]
socc [2,1,1,1,0,0,0,0]
units angstrom
guess read
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}

scf_e, scf_wfn = energy(‘scf’, return_wfn=True)
clean()

This example can be carried out by a generic user, and again, the second calculation reads in the correctly rotated orbital of the first calculation then proceed:


==> Iterations <==

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

@ROHF iter 0: -1149.76511992142537 -1.14977e+03 0.00000e+00
@ROHF iter 1: -1146.83260510739615 2.93251e+00 1.25067e-01
@ROHF iter 2: -1148.37716908495418 -1.54456e+00 8.80257e-02 DIIS
@ROHF iter 3: -1148.73807692682840 -3.60908e-01 6.59311e-02 DIIS
@ROHF iter 4: -1149.62092330349151 -8.82846e-01 2.52979e-02 DIIS
@ROHF iter 5: -1149.72238827031811 -1.01465e-01 1.06204e-02 DIIS
@ROHF iter 6: -1149.74272482668516 -2.03366e-02 3.29663e-03 DIIS
@ROHF iter 7: -1149.74602651547298 -3.30169e-03 1.14398e-03 DIIS
@ROHF iter 8: -1149.74643330082722 -4.06785e-04 4.82710e-04 DIIS
@ROHF iter 9: -1149.74659284793074 -1.59547e-04 2.15934e-04 DIIS
@ROHF iter 10: -1149.74661965576229 -2.68078e-05 1.05742e-04 DIIS/MOM
@ROHF iter 11: -1149.74663288684405 -1.32311e-05 6.65788e-05 DIIS/MOM
@ROHF iter 12: -1149.74663747837712 -4.59153e-06 6.76452e-05 DIIS/MOM
@ROHF iter 13: -1149.74658085518308 5.66232e-05 8.06763e-05 DIIS/MOM
@ROHF iter 14: -1149.74657387083630 6.98435e-06 2.32184e-05 DIIS/MOM
@ROHF iter 15: -1149.74657393380426 -6.29680e-08 5.65871e-06 DIIS/MOM
@ROHF iter 16: -1149.74657393364441 1.59844e-10 2.52556e-06 DIIS/MOM
@ROHF iter 17: -1149.74657393575762 -2.11321e-09 8.85975e-07 DIIS/MOM

Again, allow me to quote MOLPRO result. A value of -1149.86469181 was obtained from MOLPRO ROHF-SCF.
The first iteration was approaching the correct state first, then suddenly the second iteration destroys everything and brought the occupation back to the wrong one.

-Rulin

Dear developers,

I have been following GitHub, too.
There is a newly updated fix #541, which talks about the occupation in ROHF. Is that related to this problem at all?
That fix seems like to me fixes the occupation so it would not read the occupation from the READ in orbitals.

Thanks!
-Rulin

Formatting the provided input file to make it easier for others to assist.

memory 1050 mb

def orb_rotate(C, h, i, j, theta):

    if h > C.nirrep():
        raise p4util.PsiException("ROTATE: Expression %s irrep number is larger than the number of irreps" % (str(h)))

    if max(i, j) > C.coldim()[h]:
        raise p4util.PsiException("ROTATE: Expression %s orbital number exceeds number of orbitals in irrep" % (str(max(i, j))))

    theta = np.deg2rad(theta)

    x = C.nph[h][:, i].copy()
    y = C.nph[h][:, j].copy()

    xp = np.cos(theta) * x - np.sin(theta) * y
    yp = np.sin(theta) * x + np.cos(theta) * y

    C.nph[h][:, i] = xp
    C.nph[h][:, j] = yp

molecule Mn {
0 6
Mn
}

set {
reference rohf
basis cc-pvdz
maxiter 250
docc [4,0,0,0,0,2,2,2]
socc [2,1,1,1,0,0,0,0]
units angstrom
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}

scf_e, scf_wfn = energy('scf', return_wfn=True)
clean()

orb_rotate(scf_wfn.Ca(), 0, 3, 4, 90.0)

filename = core.get_writer_file_prefix(scf_wfn.molecule().name()) + ".180.npz"
data = {}
data.update(scf_wfn.Ca().np_write(None, prefix="Ca"))
data.update(scf_wfn.Cb().np_write(None, prefix="Cb"))

Ca_occ = scf_wfn.Ca_subset("SO", "OCC")
data.update(Ca_occ.np_write(None, prefix="Ca_occ"))

Cb_occ = scf_wfn.Cb_subset("SO", "OCC")
data.update(Cb_occ.np_write(None, prefix="Cb_occ"))

data["reference"] = core.get_option('SCF', 'REFERENCE')
data["nsoccpi"] = scf_wfn.soccpi().to_tuple()
data["ndoccpi"] = scf_wfn.doccpi().to_tuple()
data["nalphapi"] = scf_wfn.nalphapi().to_tuple()
data["nbetapi"] = scf_wfn.nbetapi().to_tuple()
data["symmetry"] = scf_wfn.molecule().schoenflies_symbol()
data["BasisSet"] = scf_wfn.basisset().name()
data["BasisSet PUREAM"] = scf_wfn.basisset().has_puream()
np.savez(filename, **data)

set {
reference rohf
basis cc-pvdz
maxiter 250
docc [4,0,0,0,0,2,2,2]
socc [2,1,1,1,0,0,0,0]
units angstrom
guess read
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}

scf_e, scf_wfn = energy('scf', return_wfn=True)
clean()