Dear PSI4 Developers,
I am quite new to the PSI4 package. Here I am trying to extract the orbitals (mainly the AOMO coefficient matrix for a state averaged calculation). I am using PsiAPI interface. I am doing it the following way
First, I setup the PSI4 input for a diatomic molecule
########################################
########################################
psi4_option = {
‘basis’: ‘cc-pvtz’,
‘scf_type’: ‘pk’,
‘reference’: ‘rhf’,
‘frozen_docc’: [3, ],
‘active’: [7, ],
‘num_roots’: 2,
‘avg_states’: [0,1],
‘avg_weights’: [1,1],
‘MCSCF_MAXITER’: 100,
}
mol = psi4.geometry("""
0 1
Li
F 1 R
symmetry c1
units bohr
“”")
########################################
########################################
Then I do a SCF calculation for a reference geometry
########################################
ref_point = 2.0
psi4out = ‘psi4out_ref.out’
psi4.core.set_output_file(psi4out, False)
mol.R = ref_point
mol.update_geometry()
psi4.set_options(psi4_option)
try:
scf_e, scf_wfn = psi4.energy(‘scf’, mol=mol, return_wfn=True)
except Exception as inst:
print("SCF Failed for bond length ==> "+ str(ref_point))
exit()
########################################
########################################
########################################
Then I do CASSCF calculations using ref_wfn .
########################################
########################################
ref_wfn=scf_wfn
bond_info = [2.0, 4.0, 1.0]
npoint = abs(int(round((bond_info[1]-bond_info[0])/bond_info[2])))+1
for i in range(npoint):
########################################
bondlen = bond_info[0]+i*bond_info[2]
ext = format(bondlen, ‘.2f’)
########################################
psi4out = ‘psi4out_’+ext+’.out’
psi4.core.set_output_file(psi4out, False)
mol.R = bondlen
mol.update_geometry()
psi4.set_options(psi4_option)
########################################
print("calculation running for bond length ==> ", bondlen)
casscf_e, casscf_wfn = psi4.energy(‘casscf’, mol=mol, return_wfn=True, ref_wfn=ref_wfn)
ref_wfn = casscf_wfn
c_casscf = np.asarray(casscf_wfn.Ca())
########################################
########################################
However, I am afraid that the C matrix I am extracting (c_casscf) is not the the C matrix corresponding to each geometry! To check that I have calculated the overlap matrix (C_transposeS_overlapC), but i am not getting the Identity matrix which corresponds to the MO overlap matrix.
Can anyone please point out what I am doing wrong or how to extract the C matrix of a state CASSCF calculation correctly.
As I am a new user, I cannot upload my python script, So, I am copy pasting the same here!
################################################################################
import psi4
import numpy as np
bondlen_init = 2.0
bondlen_fin = 4.0
bond_gap = 2.0
bond_info = [None]*3
bond_info[0] = bondlen_init
bond_info[1] = bondlen_fin
bond_info[2] = bond_gap
psi4_option = {
‘basis’: ‘cc-pvtz’,
‘scf_type’: ‘pk’,
‘reference’: ‘rhf’,
‘frozen_docc’: [3, ],
‘active’: [7, ],
‘num_roots’: 2,
‘avg_states’: [0,1],
‘avg_weights’: [1,1],
‘MCSCF_MAXITER’: 100,
}
mol = psi4.geometry("""
0 1
Li
F 1 R
symmetry c1
units bohr
“”")
########################################
Reference point calculation
########################################
ref_point = 2.0
########################################
psi4out = ‘psi4out_ref.out’
psi4.core.set_output_file(psi4out, False)
mol.R = ref_point
mol.update_geometry()
psi4.set_options(psi4_option)
########################################
########################################
try:
scf_e, scf_wfn = psi4.energy(‘scf’, mol=mol, return_wfn=True)
except Exception as inst:
print("SCF Failed for bond length ==> "+ str(ref_point))
exit()
########################################
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option(‘BASIS’))
mints = psi4.core.MintsHelper(wfn.basisset())
s_ref = np.asarray(mints.ao_overlap())
c_ref = np.asarray(scf_wfn.Ca())
ovl_ref = np.matmul(np.transpose(c_ref), np.matmul(s_ref, c_ref))
fw = open(‘ovl_ref’, ‘w’)
for irow in range(ovl_ref.shape[0]):
for icol in range(ovl_ref.shape[1]):
oneline = str(irow)+’ ‘+str(icol)+’ ‘+’{:{width}.{prec}f}’.format(ovl_ref[irow,icol], width=12, prec=9)+’\n’
fw.write(oneline)
fw.close()
########################################
ref_wfn=scf_wfn
########################################
########################################
Calculation at different points
########################################
npoint = abs(int(round((bond_info[1]-bond_info[0])/bond_info[2])))+1
for i in range(npoint):
########################################
bondlen = bond_info[0]+i*bond_info[2]
ext = format(bondlen, ‘.2f’)
########################################
psi4out = ‘psi4out_’+ext+’.out’
psi4.core.set_output_file(psi4out, False)
mol.R = bondlen
mol.update_geometry()
psi4.set_options(psi4_option)
########################################
print("calculation running for bond length ==> ", bondlen)
casscf_e, casscf_wfn = psi4.energy(‘casscf’, mol=mol, return_wfn=True, ref_wfn=ref_wfn)
ref_wfn = casscf_wfn
########################################
print(“CASSCF ENERGY=”,casscf_e)
########################################
cfile = ‘C_’+ext
fc = open(cfile,‘w’)
c_casscf = np.asarray(casscf_wfn.Ca())
for irow in range(c_casscf.shape[0]):
for icol in range(c_casscf.shape[1]):
oneline = str(irow)+’ ‘+str(icol)+’ ‘+’{:{width}.{prec}f}’.format(c_casscf[irow,icol], width=12, prec=9)+’\n’
fc.write(oneline)
fc.close()
########################################
wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option(‘BASIS’))
mints = psi4.core.MintsHelper(wfn.basisset())
s = np.asarray(mints.ao_overlap())
########################################
# overlap calcuation at each points
ovl = np.matmul(np.transpose(c_casscf), np.matmul(s, c_casscf))
fx = 'ovl_'+ext
fw = open(fx, 'w')
for irow in range(ovl.shape[0]):
for icol in range(ovl.shape[1]):
oneline = str(irow)+' '+str(icol)+' '+'{:{width}.{prec}f}'.format(ovl[irow,icol], width=12, prec=9)+'\n'
fw.write(oneline)
fw.close()
########################################
# Overlap calculation using reference overlap matrix
ovl = np.matmul(np.transpose(c_casscf), np.matmul(s_ref, c_casscf))
fx = 'ovl_s_ref_'+ext
fw = open(fx, 'w')
for irow in range(ovl.shape[0]):
for icol in range(ovl.shape[1]):
oneline = str(irow)+' '+str(icol)+' '+'{:{width}.{prec}f}'.format(ovl[irow,icol], width=12, prec=9)+'\n'
fw.write(oneline)
fw.close()
########################################
########################################
It would be great help if you could shed some light in this regards.
Best regards,
Sudip