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.
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.