About ROHF orbital rotation

Dear developers,

I have felt the pressing need of ROHF orbital rotation, but given the PSI4 user manual I am a bit confused about how to apply the functionality.

I have found two ways of rotating molecular orbitals in the manual:

  1. MCSCF_ROTATE (DETCI) option, this option is referring to the DETCI module, which does not contain information about the feasibility for ROHF MOs. If I am going to perform a calculation like the following:

molecule Cr {
0 7
Cr
}

set {
reference rohf
basis diatomic_awcvqzdk
basis_relativistic diatomic_awcvqzdk_decon
DAMPING_PERCENTAGE 95
die_if_not_converged False
maxiter 500
docc [3,0,0,0,0,2,2,2]
socc [3,1,1,1,0,0,0,0]
relativistic x2c
frozen_docc [2,0,0,0,0,1,1,1]
units angstrom
MCSCF_ROTATE [1,14,19,90]
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}

energy(‘scf’)

set {
reference rohf
basis diatomic_awcvqzdk
basis_relativistic diatomic_awcvqzdk_decon
DAMPING_PERCENTAGE 95
die_if_not_converged False
maxiter 500
docc [3,0,0,0,0,2,2,2]
socc [3,1,1,1,0,0,0,0]
relativistic x2c
frozen_docc [2,0,0,0,0,1,1,1]
units angstrom
guess read
print_MOs true
MOM_START 10
print 2
SCF_TYPE PK
}
energy(‘scf’)

Is the second calculation going to read the rotated orbital of the first one?


2.A combination of options: ROTATE_MO_ANGLE (MCSCF), ROTATE_MO_IRREP (MCSCF), ROTATE_MO_P (MCSCF), ROTATE_MO_Q (MCSCF). This one is specifically for MCSCF(Including ROHF according to the manual), but I am not quite sure how to use this approach, do I just assign the values to the above options and leave it there? Should there be some kind of option=true to active the orbital rotation?


This is all I have, any help would be deeply appreciated.

Thank you!
-Rulin

Those keywords are specific to the DETCI and MCSCF modules and will not rotate the SCF orbitals. Unfortunately, I do not believe we have a conventional way of rotating SCF orbitals to satisfy your requirements.

Is your goal the global minima or a local state? I can likely show you how to hack something into the driver if you wish.

Thanks for the reply, my purpose is to calculate the global minima, but given the basis sets and level of theory it sometimes inevitably falls into a local trap(which if you looked at the electronic occupation by print MOs, is completely wrong), this often happens when lots of openshell electrons are present or the different angular momentum functions are required in the same irrep(For example, both d0 and d+2 and s are in Ag symmetry, it is sometimes easy to get s^1d^1 occupation while what you really need is d^2).

Yes, show me the trick if it is convenient for you.

-Rulin

If you are after the global minima, have you tried the stability_analysis follow command yet?

Yes, I have tried. The problem is scf itself is extremely difficult to converge, so I am hoping to use a better initial guess which involves orbital rotation to get the correct occupation to start with, this is usually how a problem of this kind is tackled with MOLPRO. The “stability_analysis follow” does not really help as it is still giving me a similar occupation which is not what I want, this could be partially due to the unconverged scf orbitals.

What we can do is replace the SCF read data:

import numpy as np
molecule Cr {
0 1
Cr
}

set {
reference rohf
basis diatomic_awcvqzdk
basis_relativistic diatomic_awcvqzdk_decon
DAMPING_PERCENTAGE 95
die_if_not_converged False
maxiter 500
docc [3,0,0,0,0,2,2,2]
socc [3,1,1,1,0,0,0,0]
relativistic x2c
frozen_docc [2,0,0,0,0,1,1,1]
units angstrom
print_MOs true
MOM_START 10
print 2
SCF_TYPE pk
}

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

Then we can write a quick python-based orbital rotation function:

# Quick orbital rotation code
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

Rotate any desired columns.

# Rotate columns here
# Switches the 1Ag irrep, third occupied orbital with fourth singly occupied orbital
orb_rotate(scf_wfn.Ca(), 0, 2, 3, 90.0)

Write out the SCF read data:

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 read and recompute!

set {
guess read
}

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

No need to specify the options twice, they are carried over. Combine all of the above pieces into a single input file and you should be good to go. This seemed slightly easier than hacking away at the driver.

1 Like

This is great!
Just one little question, I am not so sure about how wfn.Ca() stores MOs, does orb_rotate(scf_wfn.Ca(), 0, 2, 3, 90.0) mean to rotate the third and fourth orbitals within Ag irrep or within the entire MOs, namely is the number 2 and 3 here referring to the orbital number in each irrep or the orbital number in all the irreps in the order of energies?

Only orbitals are counted within the 1Ag irrep. So if I wanted to rotate the last irrep:

orb_rotate(scf_wfn.Ca(), 7, 1, 2, 90.0)

This would switch the B3u HOMO/LUMO.

Splendid!
Thanks a lot!

Hey all.

I am attempting orbital rotations for my work on the X and A states of the ethylperoxy radical, but I’ve been unsuccessful thus far in replicating your success. I have a feeling that the answer may lie in the “guess” option, because when I pass the converged ground state as my next guess, I obtain the same pattern of no convergence - whether or not I’ve attempted to rotate the orbitals using this code, or my own. Any suggestions?

Thanks in advance!

For reference, here is my input sans orbital rotation. The first SCF procedure convergences, but the second does not.

set { 
    basis 6-31G
    reference rohf
    die_if_not_converged False
    maxiter 500
    e_convergence 10
    d_convergence 10
    memory 60000
    scf_type pk
    freeze_core True
    docc [16]
    socc [1]
}

molecule A_G {
0 2
C  0.1464604158 -0.5367025185 -1.5160893616
C -0.3244778801  0.5832145568 -0.6104005954
H  1.2299209247 -0.5061553528 -1.6267460923
H -0.3096873161 -0.4204070726 -2.5008003884
H -0.1368517924 -1.5098957044 -1.1180395184
O  0.3164337434  0.5619703322  0.6808322459
H -1.4064302258  0.5620822257 -0.4605196798
H -0.0393765366  1.5620908218 -0.9933453042
O -0.1411298020 -0.5771854991  1.3366659546
}

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

set scf {
    guess read
}

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

Andrew

Also, I’ve noticed that cc-pVDZ is able to converge on the second run, but similar large jumps in energy occur within the first few iterations (as with 6-31G). I would imagine that, were this functioning as I think it should, converged guesses should work just fine - independent of basis set choice.

As far as I know, the problem has not been fixed, I encountered the same issue:

And this closed/fixed request on Git may help, but it did not seem to work either:

-Rulin

The ROHF Wavefunctions are not stable so diagonalization pops them out of any minimum that they found and is not guaranteed to return there. Im still not quite sure what we should do about this. It may be that we need to implement a new SCF solver that only does rotations and avoids diagonalization. However, this would suggest that the initial ROHF Wavefunction is not a global minima to being with, and you may not want that solution.

Unfortunately, I dont know how to solve this for ROHF. If anyone has suggestions and its fairly straightforward I would be happy to code it up.

Maybe this is a silly example, but if specify an rohf reference for a closed shell molecule (CO), and then use these orbitals as an initial guess, the 1st calculation looks fine (just like the rhf case) but in the guess=read case at the 2nd iteration things go crazy and the calculation actually never converges (the 1st iteration yields the converged RHF/ROHF/UHF result). Certainly this particular ROHF wavefunction is stable. Both RHF (of course) and UHF references work as expected.

#! test
memory 450 mb

molecule CO {
0 1
C
O 1 r
r=1.6
}

set {
reference rohf
basis cc-pvdz
maxiter 150
docc [5,0,1,1]
units angstrom
SCF_TYPE PK
}

energy(‘scf’)

set scf {
guess read
}

energy(‘scf’)

@Kirk Thanks for pushing on this a bit, you inspired me to go bug hunting and sure enough: https://github.com/psi4/psi4/pull/570

Fantastic, thanks!
Have a great weekend

Using Daniel’s fix, I got both the user-defined SCF guess and his orbital rotation code to work for my system. Thanks a bunch!!!

I compiled the latest changes this afternoon, and my second set of SCF iterations won’t converge (just as before). I worked my way back to Daniel’s commit ( https://github.com/psi4/psi4/commit/5b0aba2a3e933796eef398c1bf9120e13379b6a3 ), but no luck. Any thoughts?

With the latest version of Psi4 I get the following for the second ROHF call:

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

@ROHF iter   0:  -228.14502747694382   -2.28145e+02   1.42814e-11
@ROHF iter   1:  -228.14502747694371    1.13687e-13   2.24159e-11

Using your A_G molecule input script supplied above. Do you have a different input file to test out?

I’ve isolated the problem. The path to the Psi4 executable was recently changed, and I hadn’t realized that. Thank you for taking the time to double check yours, though - it helped me figure this out!