Computing the frozen core energy

1

Is there a simple way of getting the frozen core energy?
Looking into the documentation I found that psi4.core.Wavefunction has an efzc() method but calling it after setting freeze_core to true always returns 0.
For example:

import psi4
mol = psi4.geometry("""
O                   
H 1 0.958           
H 1 0.958 2 104.4776
                 
symmetry c1        
""")                   
psi4.set_options(                
    {                        
        'basis': 'cc-pvdz',   
        'scf_type': 'direct',
        'reference': 'rhf',  
        'freeze_core': 'true'
    }                        
)                            
_, wfn = psi4.energy('mp2', return_wfn=True)
print(wfn.efzc())

Returns 0.0.

2

efzc() is not currently used by Psi.

What do you mean by “the frozen core energy,” exactly?

3

I am not sure I fully understand what it is that I want, I was hoping that psi4 would give me a well-defined quantity and then go from there.

I am trying to understand how one goes about freezing some core electrons in a Full CI approach. From what I understand, you take out some of the core electrons from the CI Hamiltonian and then diagonalize it. In the end you have to add to the CI energy the frozen core energy, which I was hoping Psi4 would give me without having to think about it. I guess you could compute it using
E_{frozen-core} = sum_a^Nf <a|h|a> + 1/2 sum_{a,b}^N <ab||ab>,
but I’m not sure.

4

Your understanding of the frozen-core energy appears to be correct. The integral-transformation code computes this quantity automatically when it sets up the two-electron integral transformation so that, if one needs to compute the energy associated with the reference determinant, the contribution of the frozen MOs can be included correctly. (See libtrans/integraltransform_sort_so_tei.cc). However, I don’t think we currently return this quantity to the Python layer. It wouldn’t be hard to do this, though.

5

In theory, it shouldn’t be. Just put the libtrans member variable frozen_core_energy_ on the wavefunction object. However, I have the memory that the handling of frozen core in libtrans is wonky, and it equivocates between “the frozen core space” and “the space I don’t want to transform,” which makes the logic messy. However adds it needs to be very careful that libtrans does what they think it does.

For reference, libtrans is still on my list to refactor, but I’m spending more of my Psi time reviewing than coding, nowadays.

6

Thank you @crawdad and @jmisiewicz!

For future reference here’s a small program that computes the frozen core energy. I modified the Psi4 code to print the frozen_core_energy_ variable and it does seem to give the same result.

import psi4
import numpy as np

mol = psi4.geometry('''
        H 0 0 0 
        Cl 0 0 2

        symmetry c1
        ''')

psi4.set_options(
        {   
            'basis': 'sto-3g',
            'scf_type': 'direct', 
            'reference': 'rhf',
            'freeze_core': 'true'
        }
    )

# Do SCF to get the molecular orbitals (MOs)
scf_e, scf_wfn = psi4.energy('scf', return_wfn=True)

# These are the MOs
Ca = scf_wfn.Ca()

# Helper for integrals
basis = scf_wfn.basisset()
mints = psi4.core.MintsHelper(basis)

# Get the atomic orbital 1-el integrals <i|H|j>
Hao = mints.ao_kinetic()
Hao.add(mints.ao_potential())

# Transform them into molecular orbital representation
Hmo = psi4.core.triplet(Ca, Hao, Ca, True, False, False)

# Number of frozen orbitals
Nf = scf_wfn.nfrzc()

# ERIs in the MO basis 
Rm = mints.mo_eri(Ca, Ca, Ca, Ca).np

efc = 2*np.trace(Hmo.np[:Nf,:Nf]) +\
        2*np.einsum('iijj', Rm[:Nf,:Nf,:Nf,:Nf]) -\
        np.einsum('ijij', Rm[:Nf,:Nf,:Nf,:Nf])

print('Frozen core energy: ', efc)

It doesn’t return the exact same number, for this particular example, the frozen_core_energy_ variable is -442.70246464960132471, whereas the code above returns -442.70246462996215.
This is a good reference for the matter https://doi.org/10.1016/S0065-3276(08)60532-8.

7

My guess is that the disagreement is because the SCF equations need to be converged more tightly.

And yes, the Psi developers are very familiar with Prof. Sherrill’s work.