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