Extracting orbitals for a state averaged calculation

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

Can you use triple ``` blocks so that your code is formatted? I cannot quite see what is happening with the current formatting.

Is this okay??

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()

I am sorry. Somehow, it is not showing the spaces and thus, it the format is messed up. Is there any way I can attach my python script here? That would be really helpful.

surround code with triple backticks (like on GitHub) for formatting. https://discourse.stonehearth.net/t/discourse-guide-code-formatting/30587

Thanks for the link. Here is the code

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()

I played around with this some. My hypothesis is that when you supply ref_wfn to psi4.energy, the code blindly takes the overlap matrix from the ref_wfn. You pass in a ref_wfn from a different geometry, so it uses a bad overlap matrix. The proper test for this would be to pass in no reference wavefunction or a reference wavefunction from the correct geometry and see if that cures your convergence problem.

Unfortunately, the CAS won’t even converge when I try either of those, so I’ll leave further investigation to somebody who has used Psi’s CAS code before.

This might help: https://github.com/psi4/psi4/issues/758