Hi,
I am looking at a way to use the previously determined MCSCF orbitals (like SA-CASSCF for instance, or state-specific CASSCF) to improve speed and convergence in my calculations.
Typically when I plot a 1D-PES, I would like to use the MCSCF orbitals of point n-1 as a starting guess for the MCSCF of point n.
For SCF calculations, this is easily done by setting set guess read. Is there any analog for mcscf ? I couldn’t find it in the doc.
When you say “previously determined,” do you mean from (potentially) a previous geometry? This is an important point.
Yes exactly. Let’s take a simple example: H2. You compute the MCSCF orbitals for R=1.2 and then you want to use this guess for the next point, i.e. R = 1.25 for instance. I expect that in some particular cases, it might help to use this guess as a starting point (?).
To put things in a context, I’m looking at a conical intersection for formaldimine with a CAS(4,4). SA-CASSCF is required here, but I’d like to show that state-specific CASSCF will fail. It indeed fails and have a very erratic behavior. I think I can improve a bit this behavior by starting with previously optimized orbital from the last geometry. And sometimes, the MCSCF (being state-average or not) also struggle to converge, depending on the number of virtual that I let free to be optimized. Maybe setting this will help as well.
Oh, it should help as a guess, no question.
The only mechanism for orbital passing into a correlated method is to get them from the wavefunction… Which takes us right back to the fact that serialized wavefunctions with CASSCF aren’t quite working. I can’t think of a way to make this work offhand. @york0822 may also have thoughts on this.
Yes ok, I was hopping that we could bypass this problem but indeed that’s actually the same issue… 
Are you running energy('casscf') or energy('mcscf')? Two different modules.
The latter does not accept a reference wavefunction.
If you can generate the PES within a single python instance you do not need to serialize the wavefunction which should avoid the linked issue, i believe.
Typically I’m using such kind of code (here I give an example for the ethylene twisted pyramidalization)
import psi4
import numpy as np
def get_energies(alpha,ref_wfn_casscf=False):
R_CC = 1.334
R_CH = 1.092
a_HCH = 116.8*np.pi/180.0
# hydrogens of the first carbon atom
y_H1 = - R_CC - R_CH*np.cos(a_HCH/2.0) # = y_H2
z_H1 = R_CH*np.sin(a_HCH/2.0) # = - z_H2
# x_H1 = x_H2 = 0
# hydrogens of the second carbon atom evolve on the zy plane centered on (x_H3,0,0) and radius R_CH*cos(phi)
# so one has to fulfil y^2 + z^2 = (R_CH*cos(phi))^2 for a given pyramidal angle alpha... phi is fixed
# say alpha = 0 corresponds to the twisted and non pyramidal ethylene
# increasing alpha towards 90° put the C-H parallele to z axis, H on top
phi = a_HCH/2.0
y_H3 = R_CH*np.cos(phi)*np.cos(alpha) # = y_H4
x_H3 = R_CH*np.sin(phi) # = -x_H4
z_H3 = R_CH*np.cos(phi)*np.sin(alpha)
# all variables:
variables = [-R_CC,y_H1,z_H1,x_H3,y_H3,z_H3]
mol = psi4.geometry(
"""0 1
C 0. {0} 0.
H 0. {1} {2}
H 0. {1} -{2}
C 0. 0. 0.
H {3} {4} {5}
H -{3} {4} {5}
symmetry c1""".format(*variables))
# Computation options
psi4.core.clean_options()
psi4.set_options({'basis': 'cc-pVDZ',
'reference': 'rhf',
'scf_type': 'pk', # set e_convergence and d_convergence to 1e-8 instead of 1e-6
'num_roots': 2,
'frozen_docc':[0],
'restricted_docc': [7],
'active':[2],
'restricted_uocc':[39],
'frozen_uocc':[0],
'mcscf_diis_start': 10, # default is 3
'mcscf_maxiter': 300})
casscf,casscf_wfn = psi4.energy('casscf',return_wfn=True)
# WHAT I TRIED TO DO BUT I DOESN'T ALLOW NEW OPTIMISATION OF ORBITALS
#if ref_wfn_casscf is not False:
# casscf,casscf_wfn = psi4.energy('casscf',molecule=mol,ref_wfn = ref_wfn_casscf,return_wfn=True)
E0_casscf = psi4.variable('CI ROOT 0 TOTAL ENERGY')
E1_casscf = psi4.variable('CI ROOT 1 TOTAL ENERGY')
return E0_casscf, E1_casscf, casscf_wfn
angle = []
E0_casscf = []
E1_casscf = []
for alpha in np.linspace(0,10,5):
angle.append(alpha)
if alpha == 0:
en0_casscf,en1_casscf,casscf_wfn = get_energies(alpha*np.pi/180.0)
else:
en0_casscf,en1_casscf,casscf_wfn = get_energies(alpha*np.pi/180.0,ref_wfn_casscf=casscf_wfn)
E0_casscf.append(en0_casscf)
E1_casscf.append(en1_casscf)
print(E0_casscf)
print(E1_casscf)
If you uncomment the attempt I did with the restart wave function, we will normally get the exact same CASSCF energy whatever the geometry. So I guess it just do not reoptimize the orbitals again such that you always end up with the same result.
Indeally one can pass the previous wavefunction (say R = 1.2) to the current computation (R = 1.25) by doing
# old geometry
molecule H2 {
H
H 1 r
}
H2.r = 1.2
Ecas, wfn = energy("casscf", return_wfn=True)
# new geometry
H2.r = 1.25
energy("casscf", ref_wfn=wfn)
However, I am not sure the extra orthogonalization of the basis is now implemented correctly in Psi4 @jmisiewicz. If not, please refer to the script in this issue https://github.com/psi4/psi4/issues/758.
So as you can see in my code, I tried to do that already here:
# WHAT I TRIED TO DO BUT I DOESN'T ALLOW NEW OPTIMISATION OF ORBITALS
if ref_wfn_casscf is not False:
casscf,casscf_wfn = psi4.energy('casscf',molecule=mol,ref_wfn = ref_wfn_casscf,return_wfn=True)
but unfortunately, I just obtain the exact same CASSCF energy as for the previous geometry. So I’m assuming that passing the reference wavefunction works, but the orbitals are then kept as is and not reoptimized by a new casscf with the new geometry.
I will take a close look at the script you mention.
Yes, I missed that part. However, you should not get the same energy but a wrong energy (presumably lower in energy).
molecule H2{
H
H 1 r
}
set {
reference rhf
basis 3-21g
active [1,0,0,0,0,1,0,0]
restricted_docc [0,0,0,0,0,0,0,0]
}
H2.r = 1.0
Ecas, wfn = energy('casscf', return_wfn=True)
H2.r = 1.1
energy('casscf', ref_wfn=wfn)
The energies I got are -1.120277767817416 for 1.0 and -1.168384786969245 for 1.1, which seems incorrect. I think the correct energy for 1.1 should be -1.105415016659182.
The reason is that the orbitals at r = 1.0 are orthogonal C_1^+ S_1 C_1 = 1 and when passing to r = 1.1 we get C_1^+ S_2 C_1 != 1. That script is just trying to recover this property.
Indeed I get the same as you with your example. I don’t know why mine doesn’t work (I really get the same energies), I’ll try to make it simpler and see what is going wrong.
So I made the same script as you did but with import psi4 in a python script and I indeed miss something:
r = 1.0
psi4.geometry="""0 1
H
H 1 {0}""".format(r)
psi4.set_options({'basis': '3-21g',
'reference':'rhf',
'active': [1,0,0,0,0,1,0,0],
'restricted_docc':[0,0,0,0,0,0,0,0]})
Ecas, wfn = psi4.energy('casscf', return_wfn=True)
r = 1.1
psi4.geometry="""0 1
H
H 1 {0}""".format(r)
psi4.energy('casscf', ref_wfn=wfn)
gives the following error:
Traceback (most recent call last):
File "test_H2.py", line 14, in <module>
Ecas, wfn = psi4.energy('casscf', return_wfn=True)
File "/Users/bsenjean/Git_repositories/psi4/objdir/stage/lib/psi4/driver/driver.py", line 530, in energy
molecule.update_geometry()
AttributeError: 'NoneType' object has no attribute 'update_geometry'
So that’s why I don’t get different result: the geometry is certainly not updated, that’s all, when I use a restart wavefunction. I tried to get this update here but without any success, how should I proceed?
Also, I tried your script to orthogonalize the orbitals for the second geometry of H2 here, and it didn’t work straightforwardly. From this input:
molecule H2{
H
H 1 r
}
set {
reference rhf
basis 3-21g
active [1,0,0,0,0,1,0,0]
restricted_docc [0,0,0,0,0,0,0,0]
}
H2.r = 1.0
Ecas, wfn = energy('casscf', return_wfn=True)
print("CASSCF ORBITALS GEOMETRY 1")
print(wfn.Ca().nph)
H2.r = 1.1
Escf, wfnSCF = energy('scf', return_wfn=True)
print("SCF ORBITALS GEOMETRY 2")
print(wfnSCF.Ca().nph)
# import the script <== Change HERE
sys.path.insert(0, '/Users/bsenjean/Git_repositories/qchem_simulation/Notes/Embedding_project/AVAS')
from orthogonalize_orbitals import ortho_orbs
wfnSCF.Ca().copy(ortho_orbs(wfn,wfnSCF))
print("CASSCF ORBITALS FROM GEOMETRY 1 BUT ORTHOGONALIZED FOR GEOMETRY 2")
print(wfnSCF.Ca().nph)
Ecas2, wfn = energy('casscf', ref_wfn=wfnSCF, return_wfn=True)
print(Ecas)
print(Ecas2)
I get this output:
CASSCF ORBITALS GEOMETRY 1
(array([[ 0.37092414, 1.24482581],
[ 0.51019803, -0.9571024 ]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([[ 0.71712113, 1.15700169],
[ 0.91643852, -2.11514574]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64))
SCF ORBITALS GEOMETRY 2
(array([[ 0.33795149, 1.28231161],
[ 0.54922924, -0.94958934]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([[ 0.22654256, 1.31537109],
[ 1.50900811, -1.5438462 ]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64))
CASSCF ORBITALS FROM GEOMETRY 1 BUT ORTHOGONALIZED FOR GEOMETRY 2
(array([[ 0.37764796, 1.27118697],
[ 0.5194465 , -0.96620288]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64), array([[ 0.67716336, 1.15020539],
[ 0.86537484, -1.97780002]]), array([], shape=(0, 0), dtype=float64), array([], shape=(0, 0), dtype=float64))
invalid value encountered in less_equal
Traceback (most recent call last):
File "/Users/bsenjean/Git_repositories/psi4/objdir/stage/bin/psi4", line 331, in <module>
exec(content)
File "<string>", line 40, in <module>
File "/Users/bsenjean/Git_repositories/psi4/objdir/stage/lib/psi4/driver/driver.py", line 570, in energy
wfn = procedures['energy'][lowername](lowername, molecule=molecule, **kwargs)
File "/Users/bsenjean/Git_repositories/psi4/objdir/stage/lib/psi4/driver/procrouting/proc.py", line 4804, in run_detcas
ciwfn = mcscf.mcscf_solver(ref_wfn)
File "/Users/bsenjean/Git_repositories/psi4/objdir/stage/lib/psi4/driver/procrouting/mcscf/mcscf_solver.py", line 231, in mcscf_solver
total_step = diis_obj.extrapolate()
File "/Users/bsenjean/Git_repositories/psi4/objdir/stage/lib/psi4/driver/p4util/solvers.py", line 296, in extrapolate
invB.power(-1.0, 1.e-12)
RuntimeError:
Fatal Error: Matrix::power: C_DSYEV failed
Error occurred in file: /Users/bsenjean/Git_repositories/psi4/psi4/src/psi4/libmints/matrix.cc on line: 2280
The most recent 1 function calls were:
Printing out the relevant lines from the Psithon --> Python processed input file:
sys.path.insert(0, '/Users/bsenjean/Git_repositories/qchem_simulation/Notes/Embedding_project/AVAS')
from orthogonalize_orbitals import ortho_orbs
wfnSCF.Ca().copy(ortho_orbs(wfn,wfnSCF))
print("CASSCF ORBITALS FROM GEOMETRY 1 BUT ORTHOGONALIZED FOR GEOMETRY 2")
print(wfnSCF.Ca().nph)
--> Ecas2, wfn = energy('casscf', ref_wfn=wfnSCF, return_wfn=True)
print(Ecas)
print(Ecas2)
!----------------------------------------------------------------------------------!
! !
! Fatal Error: Matrix::power: C_DSYEV failed !
! Error occurred in file: !
! /Users/bsenjean/Git_repositories/psi4/psi4/src/psi4/libmints/matrix.cc on line: !
! 2280 !
! The most recent 1 function calls were: !
! !
!----------------------------------------------------------------------------------!
I’m going to leave this topic to York, since he has far more relevant experience. Ping me if I’m needed.
As for basis set orthogonalization, York is correct. Psi currently does not modify the orbitals on an incoming wavefunction for orthogonality. I need to change this because an upcoming paper of mine needs this feature, but that paper needs a lot of features in Psi, so it may be some time before I get to it. 
When it’s time for me to focus on that paper again, I’ll put this at the top of my list for you, but I need to finish another paper first (for York, this is the one I needed the UCC code for). I’m not expecting that to happen for at least a month. Sorry! Until then - or if somebody else fixes it - you’re stuck with this kludgier way of doing things.
I have encountered this error before but I did not figure out why.
This problem is likely not coming from the script but the initial CI becomes nan, which appears somewhat random. For example, if you change the first geometry to 0.9 or change themcscf_type to AO, everything works fine.
I will let you know if I figure out how to fix it.
Edits: Sorry, I cannot reproduce the error using the latest Psi4 and I attached my output. Please let me know if you find a problematic case.
output.dat (26.6 KB)
Erf, that’s weird. I have the exact same psi4 version but I still get the error.
Do you use the exact same orthogonalization script that described in https://github.com/psi4/psi4/issues/758 or do you have some local modifications ?
I enclose my output: test_H2.txt (27.7 KB)
I could also try on my formaldimine case, but do you know how I should fix the update_geometry issue that I mentioned in previous replies ?
I think there is a typo in your python code.
r = 1.0
psi4.geometry("""0 1
H
H 1 {0}""".format(r)) # <= here is the error, where you set psi4.geometry as a string
psi4.set_options({'basis': '3-21g',
'reference':'rhf',
'active': [1,0,0,0,0,1,0,0],
'restricted_docc':[0,0,0,0,0,0,0,0]})
Ecas, wfn = psi4.energy('casscf', return_wfn=True)
r = 1.1
psi4.geometry="""0 1
H
H 1 {0}""".format(r) # I left the typo unchanged here.
psi4.energy('casscf', ref_wfn=wfn)
Oh right, stupid mistake !!! Now that it’s fixed, it works when I use your orthogonalization script.
Without this script, the geometry is not updated by just setting psi4.geometry again.
See the enclosed output: test_H2_withimport.txt (43.2 KB)
Generated with the following input:
import psi4
import sys
import numpy as np
psi4.core.set_output_file("test_H2_withimport.out", False)
r = 1.0
psi4.geometry("""0 1
H
H 1 {0}""".format(r))
psi4.set_options({'basis': '3-21g',
'reference':'rhf',
'active': [1,0,0,0,0,1,0,0],
'restricted_docc':[0,0,0,0,0,0,0,0]})
Ecas, wfn = psi4.energy('casscf', return_wfn=True)
# import the script <== Change HERE
sys.path.insert(0, '/Users/bsenjean/Git_repositories/qchem_simulation/Notes/Embedding_project/AVAS')
from orthogonalize_orbitals import ortho_orbs
r = 1.1
psi4.geometry("""0 1
H
H 1 {0}""".format(r))
Escf, wfnSCF = psi4.energy('scf', return_wfn=True)
wfnSCF.Ca().copy(ortho_orbs(wfn, wfnSCF, True))
Ecas2 = psi4.energy('casscf', ref_wfn=wfn)
Ecas3 = psi4.energy('casscf', ref_wfn=wfnSCF)
Ecas4 = psi4.energy('casscf')
print(Ecas)
print(Ecas2)
print(Ecas3)
print(Ecas4)
which prints:
-1.1202777678174165
-1.120277767817428
-1.1054150166591654
-1.1054150166591818
So it works 
. And we can see that the first iteration of MCSCF with the previous MCSCF orbitals is much closer to the converged MCSCF energy than with HF orbitals.
However I still don’t know why the psi4 script (called like psi4 test_H2.in) does not work the same as the python script. But it does not matter so much to me as I’m using the python scripts anyway.
That is good to know it works! I am treaking the script to handle frozen orbitals. I will upload it shortly.
Here is the updated script that should work when there are frozen orbitals. I tested on N2, which seems to work.
@bsenjean @jmisiewicz
ortho_orb.py (5.3 KB)
Thanks ! It seems to work fine !
Also, I noticed that one should set the no_reorient keyword (as well as no_com maybe) otherwise the convergence is even worst with the restarting orbitals.
Here is the script I use to test it with the pyramidalization of ethylene, if you want to try out:
def geometry(angle):
alpha = angle*np.pi/180
R_CC = 1.334
R_CH = 1.092
a_HCH = 116.8*np.pi/180.0
# hydrogens of the first carbon atom
y_H1 = - R_CC - R_CH*np.cos(a_HCH/2.0) # = y_H2
z_H1 = R_CH*np.sin(a_HCH/2.0) # = - z_H2
# x_H1 = x_H2 = 0
phi = a_HCH/2.0
y_H3 = R_CH*np.cos(phi)*np.cos(alpha) # = y_H4
x_H3 = R_CH*np.sin(phi) # = -x_H4
z_H3 = R_CH*np.cos(phi)*np.sin(alpha)
# all variables:
variables = [-R_CC,y_H1,z_H1,x_H3,y_H3,z_H3]
return variables
# cc-pVDZ : 48 orbs
psi4.core.set_output_file("ethylene_pyramidalized_withfrozen_SA.txt", False)
angle = 0
mol = psi4.geometry(
"""0 1
C 0. {0} 0.
H 0. {1} {2}
H 0. {1} -{2}
C 0. 0. 0.
H {3} {4} {5}
H -{3} {4} {5}
symmetry c1
no_reorient
no_com""".format(*geometry(angle)))
psi4.set_options({'basis': 'cc-pVDZ',
'reference': 'rhf',
'scf_type': 'pk', # set e_convergence and d_convergence to 1e-8 instead of 1e-6
'num_roots': 2,
'frozen_docc':[2],
'restricted_docc': [5],
'active':[2],
'restricted_uocc':[39],
'frozen_uocc':[0],
'mcscf_maxiter': 300,
'mcscf_diis_start': 10, # default is 3
'avg_states': [0,1],
'avg_weights': [0.5,0.5]})
casscf,casscf_wfn = psi4.energy('casscf',return_wfn=True)
angle = 1
psi4.core.set_output_file("ethylene_pyramidalized_withfrozen_restart_SA.txt", False)
mol = psi4.geometry(
"""0 1
C 0. {0} 0.
H 0. {1} {2}
H 0. {1} -{2}
C 0. 0. 0.
H {3} {4} {5}
H -{3} {4} {5}
symmetry c1
no_reorient
no_com""".format(*geometry(angle)))
Escf, wfnSCF = psi4.energy('scf', return_wfn=True)
wfnSCF.Ca().copy(ortho_orbs(casscf_wfn, wfnSCF, True))
casscf,casscf_wfn = psi4.energy('casscf',ref_wfn=wfnSCF,return_wfn=True)
I also enclose the different outputs:
without reorientation it works fine:
ethylene_pyramidalized_withfrozen_SA.txt (14.3 KB)
ethylene_pyramidalized_withfrozen_restart_SA.txt (13.2 KB)
And with reorientation it doesn’t:
ethylene_pyramidalized_withfrozen_SA_reorient.txt (14.6 KB)
ethylene_pyramidalized_withfrozen_restart_SA_reorient.txt (14.8 KB)
Note also that the threshold for the SCF calculation change from 1e-8 to 1e-6 with the restart, despite the pk option is set. But I guess it does not matter for the MCSCF calculation as the HF orbitals are not used anyway. But still, I think it might be important to have the 1e-8 for consistency.