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!
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
}
- 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.
- Your input file isn’t complete. You never even call your orbital rotation function or for an SCF.
- 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.
- Replace your orbital rotation function with
scf_wfn.Ca().rotate_columns(0, 4, 5, np.deg2rad(90)). The arguments torotate_columnsare h, i, j, rotation-in-radians, so this is a drop-in replacement for your existing function. - 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.