About ROHF orbital rotation

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!

Daniel - will this python code be expected to work when one uses an ECP? We’ve recently tried and it doesn’t seem to do anything (the 2nd SCF just returns the original energy at the 2nd iteration).

Daniel isn’t active among the Psi developers anymore.

Incompatibilities with ECPs aren’t expected. I can take a look at this, but I’ll need a complete example of what you’re doing.

So I’m trying to run RuB, which I have as a doublet delta. When I run this input file, the open-shell electron, before and after rotating, is in a d0 orbital. In my output I’m also getting a warning about using ECP:

Blockquote
!!! WARNING: ECP capability is in beta. Please check occupations closely. !!!

Input file:

#avdz rub
memory 500 mb
basis {
assign vdz-pp
assign o avdz+d
assign b avdz+d
assign s avdz+d
assign se avdz+d
assign cl avdz+d
assign f avdz+d
assign si avdz+d
assign c avdz+d
}
#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

molecule = {
units bohr
0 2
Ru
B 1 R2
R2 = 3.0290000000000004
}

  1. In future, please enclose code in backticks. The forum messes with the formatting otherwise. I had to edit your code to get it to run, and I’d rather not have the niggling doubt about whether my edits changed something.
  2. Your input file isn’t complete. You never even call your orbital rotation function or for an SCF.
  3. When I add an SCF energy call to the input file, I get an error
psi4.driver.qcdb.exceptions.BasisSetNotFound: BasisSet::construct: Unable to find a basis set for atom 1 for key BASIS among:
  Shell Entries: ['RU']
  Basis Sets: [('vdz-pp', 'vdz-pp', None)]

Are you passing extra basis set information?

Very sorry, I thought the entire file had copied over.

#avdz rub

memory 500 mb


basis {
  assign vdz-pp
  assign o avdz+d
  assign b avdz+d
  assign s avdz+d
  assign se avdz+d
  assign cl avdz+d
  assign f avdz+d
  assign si avdz+d
  assign c avdz+d
}

# 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





molecule = {
    units bohr
    0 2
    Ru
    B 1 R2
    R2 = 3.0290000000000004
    }
set {
    docc [5,1,2,2]
    socc [1,0,0,0]
    reference rohf
    print_MOs true
    print 2
    scf_type pk
    guess sad
    freeze_core true
    maxiter 200
}
R2 = 3.0290000000000004
scf_e, scf_wfn = energy('scf', return_wfn=True)


# Rotate columns here
# Switches the A1 irrep, 5th occupied orbital with first singly occupied orbital
orb_rotate(scf_wfn.Ca(), 0, 4, 5, 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 {
guess read
}
ehf_2, scf_wfn = energy('scf', return_wfn=True)
ecc_2 = energy('bccd(t)')


psi4.core.print_out("    R [bohr]       HF                   CCSD(T)\n")
psi4.core.print_out("    {:1.4f}        {:7.6f}              {:7.6f}\n".format(R2, ehf_2, ecc_2))

Yes, I am passing my own basis sets for Ru.

vdz-pp.gbs for Ru:

!

Ru     0
S    7   1.00
    159.3720000              0.0006730
     20.9016000             -0.0371610
     13.0741000              0.1665170
      4.7147400             -0.6159540
      1.1361200              0.8428300
      0.5154640              0.4477630
      0.1067670              0.0184450
S    7   1.00
    159.3720000              0.0008950
     20.9016000             -0.0396350
     13.0741000              0.1442430
      4.7147400             -0.5192870
      1.1361200              1.5854040
      0.5154640             -0.9082030
      0.1067670             -1.4720020
S    7   1.00
    159.3720000             -0.0001690
     20.9016000              0.0093920
     13.0741000             -0.0477130
      4.7147400              0.2011130
      1.1361200             -0.3734210
      0.5154640             -0.2738590
      0.1067670              0.6093840
S    1   1.00
      0.0392180              1.0000000
P    6   1.00
     12.8738000              0.0299680
      6.1475000             -0.1879300
      1.5652400              0.4999900
      0.7361940              0.4957320
      0.3263170              0.1475150
      0.1124070              0.0060880
P    6   1.00
     12.8738000             -0.0135020
      6.1475000              0.0937330
      1.5652400             -0.3159100
      0.7361940             -0.3618000
      0.3263170              0.4008770
      0.1124070              0.7537660
P    6   1.00
     12.8738000             -0.0078010
      6.1475000              0.0536560
      1.5652400             -0.1719400
      0.7361940             -0.1991700
      0.3263170              0.0744110
      0.1124070              0.5971370
P    1   1.00
      0.0391720              1.0000000
D    5   1.00
      6.2126500             -0.0245970
      2.7853900              0.1273340
      1.2685100              0.3570810
      0.5517900              0.4082070
      0.2247310              0.2812400
D    5   1.00
      6.2126500              0.0261160
      2.7853900             -0.1459400
      1.2685100             -0.4425210
      0.5517900             -0.1690200
      0.2247310              0.5052370
D    1   1.00
      0.0824110              1.0000000
F    1   1.00
      0.8500000              1.0000000
****

RU     0
RU-ECP     4     28
g potential
  1
2      1.0000000              0.0000000
s-g potential
  2
2     11.50059000           209.78649300
2      5.06857500            30.21430700
p-g potential
  4
2     10.53263400            48.75124400
2     10.19201000            97.49652900
2      4.73489200             7.86018800
2      4.50906500            15.32975100
d-g potential
  4
2      8.87797700            26.96750600
2      8.76612200            40.43230300
2      3.17019600             3.34075800
2      3.22885100             5.25635200
f-g potential
  2
2      7.82024900            -8.84752500
2      7.83964700           -11.83551800

and avdzpd.gbs for B:

!

B     0
S    9   1.00
      4.570000D+03           6.960000D-04
      6.859000D+02           5.353000D-03
      1.565000D+02           2.713400D-02
      4.447000D+01           1.013800D-01
      1.448000D+01           2.720550D-01
      5.131000D+00           4.484030D-01
      1.898000D+00           2.901230D-01
      3.329000D-01           1.432200D-02
      1.043000D-01          -3.486000D-03
S    9   1.00
      4.570000D+03          -1.390000D-04
      6.859000D+02          -1.097000D-03
      1.565000D+02          -5.444000D-03
      4.447000D+01          -2.191600D-02
      1.448000D+01          -5.975100D-02
      5.131000D+00          -1.387320D-01
      1.898000D+00          -1.314820D-01
      3.329000D-01           5.395260D-01
      1.043000D-01           5.807740D-01
S    1   1.00
      1.043000D-01           1.000000D+00
S    1   1.00
      0.0310500              1.0000000
P    4   1.00
      6.001000D+00           3.548100D-02
      1.241000D+00           1.980720D-01
      3.364000D-01           5.052300D-01
      9.538000D-02           4.794990D-01
P    1   1.00
      9.538000D-02           1.000000D+00
P    1   1.00
      0.0237800              1.0000000
D    1   1.00
      3.430000D-01           1.0000000
D    1   1.00
      0.0904000              1.0000000

I’ve been able to reproduce the problem even without ECPs. My diagnosis is that this script is broken due to Psi changing its mechanism for wavefunction reading and has been broken for multiple Psi versions.

The old script is dead, long live two lines of code.

  1. Replace your orbital rotation function with scf_wfn.Ca().rotate_columns(0, 4, 5, np.deg2rad(90)). The arguments to rotate_columns are h, i, j, rotation-in-radians, so this is a drop-in replacement for your existing function.
  2. Replace your wavefunction write with scf_wfn.to_file(scf_wfn.get_scratch_filename(180)).

When I make the second change, I get expected behavior on my ECP-free test system. I strongly encourage both changes: you’re down from 40 lines of code to 2, it’ll be easier to get support from Psi developers in future, and you’re less likely to need support from Psi developers in future.

Test it on your actual system and please report back.

This looks great, thank you!

I implemented both changes and I can rotate the orbitals! Thank you

1 Like

This topic is a monolith, so here’s a summary for people-from-the-future:

You may want to swap the HOMO and LUMO of a converged SCF to target a different SCF solution, especially in ROHF.

Let’s now assume that you have the first SCF wavefunction in-memory as scf_wfn.

You can accomplish a swap between orbitals I and J of irrep h via
scf_wfn.Ca().rotate_columns(h, i, j, np.deg2rad(90))
You can change 90 to any number of degrees, but 90 will accomplish a swap. For RHF and ROHF solutions, this is enough, but for UHF, you need to use scf_wfn.Cb() to swap the beta orbitals.

It remains to pass this off as a guess to the next set attempt at SCF. To do that, write scf_wfn.to_file(scf_wfn.get_scratch_filename(180)). Set the option guess read, run another SCF computation, and your altered orbitals should be used as the initial guess.

By the time you’re reading this, this topic is probably locked. If so, feel free to make a new topic with follow-up questions.

1 Like

This topic was automatically closed 60 days after the last reply. New replies are no longer allowed.