Reading in ROHF orbitals as SCF guess

In the past, I have successfully used UHF orbitals as an initial guess for a second SCF procedure. For example, with this input:

set basis cc-pvdz
set reference uhf

molecule h2cn {
0 2
N          0.0000000000        0.0000000000       -0.6661452159
C          0.0000000000        0.0000000000        0.5871908626
H          0.0000000000       -0.9366014853        1.1547227488
H          0.0000000000        0.9366014853        1.1547227488
}

psi4_io.set_specific_path(PSIF_SCF_MOS, '.')
psi4_io.set_specific_retention(PSIF_SCF_MOS, True)
energy('scf')

molecule h2cn {
0 2
N          0.0000000000        0.0000000000       -0.6661452159
C          0.0000000000        0.0000000000        0.5871908626
H          0.0000000000       -0.9366014853        1.1547227488
H          0.0000000000        0.9366014853        1.1547227488
symmetry c1
}

set scf guess read
energy('scf')

The second SCF job takes one iteration to converge, as I would expect it to.

SCF Guess: Guess MOs from previously saved Fock matrix.

  ==> Iterations <==

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

   @DF-UHF iter   0:   -93.44076220039035   -9.34408e+01   4.09110e-07 
   @DF-UHF iter   1:   -93.44076220094324   -5.52888e-10   3.53635e-07 

But when I submit the same input file with an ROHF reference:

set basis cc-pvdz
set reference rohf

molecule h2cn {
0 2
N          0.0000000000        0.0000000000       -0.6661452159
C          0.0000000000        0.0000000000        0.5871908626
H          0.0000000000       -0.9366014853        1.1547227488
H          0.0000000000        0.9366014853        1.1547227488
}

psi4_io.set_specific_path(PSIF_SCF_MOS, '.')
psi4_io.set_specific_retention(PSIF_SCF_MOS, True)
energy('scf')

molecule h2cn {
0 2
N          0.0000000000        0.0000000000       -0.6661452159
C          0.0000000000        0.0000000000        0.5871908626
H          0.0000000000       -0.9366014853        1.1547227488
H          0.0000000000        0.9366014853        1.1547227488
symmetry c1
}

set scf guess read
energy('scf')

the energy at the start of the second SCF procedure does not match the output energy of the first, and it takes fourteen steps to converge.

 SCF Guess: Guess MOs from previously saved Fock matrix.

  ==> Iterations <==

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

   @DF-ROHF iter   0:   -42.02176950510786   -4.20218e+01   6.55809e-02
   @DF-ROHF iter   1:   -77.79769171277506   -3.57759e+01   7.05907e-02
   @DF-ROHF iter   2:   -78.40032309173856   -6.02631e-01   6.50899e-02 DIIS
   @DF-ROHF iter   3:   -92.86416319126418   -1.44638e+01   1.92976e-02 DIIS
   @DF-ROHF iter   4:   -93.34112108956496   -4.76958e-01   8.03414e-03 DIIS
   @DF-ROHF iter   5:   -93.42147297934090   -8.03519e-02   2.06201e-03 DIIS
   @DF-ROHF iter   6:   -93.42983321175151   -8.36023e-03   4.12600e-04 DIIS
   @DF-ROHF iter   7:   -93.43015790673918   -3.24695e-04   8.55039e-05 DIIS
   @DF-ROHF iter   8:   -93.43017186325338   -1.39565e-05   1.72777e-05 DIIS
   @DF-ROHF iter   9:   -93.43017246776162   -6.04508e-07   6.86456e-06 DIIS
   @DF-ROHF iter  10:   -93.43017254558441   -7.78228e-08   3.65034e-06 DIIS
   @DF-ROHF iter  11:   -93.43017256035043   -1.47660e-08   1.82749e-06 DIIS
   @DF-ROHF iter  12:   -93.43017256448108   -4.13065e-09   2.16025e-07 DIIS

I have had this problem on two separate machines with a version of psi that is up-to-date as of this morning. I am wondering if this is a bug, or if there is something else I need to do to successfully use ROHF orbitals as an initial guess? Maybe @jgonthier or @dgasmith would have some ideas about this??

Kind of odd, on the first iteration I get:

@DF-ROHF iter   0:   -93.43017256459265   -9.34302e+01   0.00000e+00

Which is correct, but the diagonalization step takes it away from the correct result. Could be that the Wavefunction isn’t stable, could be a phase issue, etc. So nothing specific comes to mind. So two things:

  • Why do you specify the molecule twice?
  • Is there a specific reason you want this to work?

Ultimately, I am interested in computing frequencies for the doublet A_1 excited state. For symmetric geometries, I can access this state by specifying orbital occupations, but for non-symmetric displacements I need to converge the correct SCF solution using orbitals from a nearby symmetric structure as an initial guess. I have used this procedure successfully for both UHF and ROHF references in Molpro, but I would like to use Psi4’s interface to MRCC.

Can you give me an example input? The symmetric ROHF orbitals should be fine and easily cast to C1 symmetry, we can hack together a new wavefunction from there without running a second SCF computation and then pass it on as normal.

Here is a minimal input for what I am trying to do with the excited state, only using UHF reference instead of ROHF. The second geometry has a small displacement that breaks symmetry, but the second SCF solution still converges to the excited state surface.

set basis cc-pvdz
set reference uhf 
set SOCC [1,0,0,0]
set DOCC [4,0,1,2]

molecule h2cn {
0 2
N            0.000000000000     0.000000000000    -0.678990573068
C           -0.000000000000     0.000000000000     0.594584322830
H            0.000000000000    -0.953046245022     1.177259590584
H           -0.000000000000     0.953046245022     1.177259590584
}

psi4_io.set_specific_path(PSIF_SCF_MOS, '.')
psi4_io.set_specific_retention(PSIF_SCF_MOS, True)                                                                                          
energy('scf')

set SOCC [1] 
set DOCC [7] 

molecule h2cn {
0 2
symmetry c1
N            0.000000000000     0.000000000000    -0.678990573068
C           -0.000000000000     0.000000000000     0.594584322830
H            0.001000000000    -0.953046245022     1.177259590584
H           -0.000000000000     0.953046245022     1.177259590584
}

set scf guess read
energy('ccsd')

Converging the second SCF wavefunction with a perturbed geometry is a different story than just requiring C1 orbitals. I don’t think we have anything to complete the above. What we can do in the future however is to use updates through rotations instead of diagonalization starting at the first iteration. Rotations should keep you within the same well on the SCF potential energy surface.

If you are adventurous you can change line src/lib/lib_scfsolver/hf.cc:1858 from:

if (soscf_enabled_ && (Drms_ < soscf_r_start_) && (iteration_ > 3)){

to:

if (soscf_enabled_ && (Drms_ < soscf_r_start_)){

and then place set soscf true in your keywords for the second SCF run. This may trick the SOSCF module into take rotation steps on the first iteration.