# 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)

## 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))))

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
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))))

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