Calculating the first order electrostatic energy in the monomer basis using Psi4Numpy

SAPT
1

Hi all,

I would like to calculate the first order electrostatic energy in the monomer-centered basis (MCBS) using Psi4Numpy. I am aware that it is possible to obtain this in the main Psi4 code but I would like to obtain this energy directly without calculating the other components of the SAPT energy.

I have been using the tutorial notebook 7b (AO-SAPT) to do this. In helper_SAPT.py the fragments are obtained in the dimer basis using dimer.extract_subsets(1, 2) i.e. with ghost functions. In the main Psi4 code, the MCBS calculation grabs the fragments without ghost functions using dimer.extract_subsets(1).

Are SAPT-AO and MCBS SAPT the same methods? This difference in how the fragments are obtained and the AO to MO integral transformation in tutorial 7a leads me to think that these methods are not the same. The SAPT-AO method gives essentially the same electrostatic energy as the energy obtained with the main Psi4 code using energy('sapt0', sapt_basis='dimer'). This is mentioned in the tutorial, where the author mentions “we expect exactly the same numerical answers” as SAPT-MO.

Assuming that performing the calculation in the MCBS is different to SAPT-AO, I have been trying to calculate the electrostatic energy in the MCBS using Psi4Numpy, where I obtain the dimer fragments in the monomer basis without ghost functions.

I need the nuclear potential integrals between the nuclei of one system and the basis functions of the other, as well as the ERIs between the basis functions of each system, in the monomer basis. Is it possible to get these “intermolecular” integrals using MintsHelper i.e. ao_potential() or ao.eri() or some other alternative method in Psi4/Psi4Numpy?

I would assume that it is possible to get these integrals as they would be used in the MCBS SAPT calculation, but I am unsure of how the main Psi4 code obtains them and whether I could also obtain them in Psi4Numpy.

Thank you for your time.

2

MCBS and DCBS are the primary differentiators where M (monomer) computations are only performed in the respective monomer basis (only basis functions of A or B) while D (dimer) computations are performed on the entire space and account for BSSE. AO representations have no effect on monomer or dimer basis sets as the AO space can either span monomer or dimer space.

Taking A = monomer A’s basis without ghost and B = monomer B’s basis without ghosts we can compute arbitrary AO integrals in various forms AO_aabb = mints.ao_eri(A, A, B, B) would for example compute electostatic-like integrals between monomers A and B in MCBS form. You can compute mints.ao_potential(A, B), but these are not typically use since one-electron operators spanning two different spaces yields somewhat odd quantities (though they do appear is some higher order theories).

I’m not sure I fully answered your question, but hopefully this helps.

1 Like
3

Thanks a lot for your answer, this was very helpful. I can now calculate the first order electrostatic energy and it matches the one obtained with Psi4.

I have been obtaining the potential integrals in the dimer basis and taking a submatrix to obtain the interaction between the nuclei of A and the electrons of B.

Is there a way to directly obtain mintsA.ao_potential(B, B) to get this interaction? I have tried this but got the same nuclear potentials as mintsA.ao_potential(A, A), where both of these MintsHelper objects were obtained with the basis of A.

4

If mintsA was initialized with Monomer A’s basis mintsA.ao_potential(B, B) should work. Perhaps try mintsA.ao_potential(B) to mimic what is in Psi4Numpy. Though if this does work the (B, B) form is likely a bug and should be made into an issue on GitHub.

5

I have enclosed my code for the He dimer to see if I am obtaining the nuclear integrals correctly. The relevant section is under # potential integrals.

# He dimer 
dimer = psi4.geometry("""
He 0.0 0.0 0.0
--
He 3.0 0.0 0.0
symmetry c1
no_reorient
no_com
""")

psi4.set_options({'basis': '6-31g',
                  'e_convergence': 1e-8,
                  'd_convergence': 1e-8})

dimer.reset_point_group('c1')
dimer.fix_orientation(True)
dimer.fix_com(True)
dimer.update_geometry()

monomerA_DCBS = dimer.extract_subsets(1,2)

monomerA = dimer.extract_subsets(1)
monomerB = dimer.extract_subsets(2)

# SCF
psi4.core.clean()
rhf_dimA,wfn_dimA = psi4.energy('SCF', return_wfn=True, molecule=monomerA_DCBS)
psi4.core.clean()
rhfA, wfnA = psi4.energy('SCF', return_wfn=True, molecule=monomerA)
psi4.core.clean()
rhfB, wfnB = psi4.energy('SCF', return_wfn=True, molecule=monomerB)

# initialise mintshelper

mints_dimA = psi4.core.MintsHelper(wfn_dimA.basisset())

mints_A = psi4.core.MintsHelper(wfnA.basisset())
mints_B = psi4.core.MintsHelper(wfnB.basisset())

# potential integrals

# potential integrals in dimer basis
V_dimA = np.asarray(mints_dimA.ao_potential())
print("V_dimA","\n",V_dimA)

# potential integrals in monomer basis (using mints_A)
V_AA = np.asarray(mints_A.ao_potential(wfnA.basisset(), wfnA.basisset()))
V_AB = np.asarray(mints_A.ao_potential(wfnA.basisset(), wfnB.basisset()))
V_BA = np.asarray(mints_A.ao_potential(wfnB.basisset(), wfnA.basisset()))
V_BB = np.asarray(mints_A.ao_potential(wfnB.basisset(), wfnB.basisset()))

print("V_AA","\n",V_AA,"\n","V_AB","\n",V_AB,"\n","V_BA","\n",V_BA,"\n","V_BB","\n",V_BB)

V_B = np.asarray(mints_A.ao_potential(wfnB.basisset()))
print(V_B)

This gives the following output:

V_dimA 
 [[-4.80528110e+00 -2.01363618e+00 -9.57020756e-10 -4.67437262e-04]
 [-2.01363618e+00 -1.74213558e+00 -1.10937910e-04 -5.86388978e-03]
 [-9.57020756e-10 -1.10937910e-04 -3.52784807e-01 -2.23717688e-01]
 [-4.67437262e-04 -5.86388978e-03 -2.23717688e-01 -3.52784807e-01]]
V_AA 
 [[-4.8052811  -2.01363618]
 [-2.01363618 -1.74213558]] 
 V_AB 
 [[-9.57020756e-10 -4.67437262e-04]
 [-1.10937910e-04 -5.86388978e-03]] 
 V_BA 
 [[-9.57020756e-10 -4.67437262e-04]
 [-1.10937910e-04 -5.86388978e-03]] 
 V_BB 
 [[-4.8052811  -2.01363618]
 [-2.01363618 -1.74213558]]
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-97666792cad0> in <module>
     54 print("V_AA","\n",V_AA,"\n","V_AB","\n",V_AB,"\n","V_BA","\n",V_BA,"\n","V_BB","\n",V_BB)
     55 
---> 56 V_B = np.asarray(mints_A.ao_potential(wfnB.basisset()))

TypeError: ao_potential(): incompatible function arguments. The following argument types are supported:
    1. (self: psi4.core.MintsHelper) -> psi4.core.Matrix
    2. (self: psi4.core.MintsHelper, arg0: psi::BasisSet, arg1: psi::BasisSet) -> psi4.core.Matrix

Invoked with: <psi4.core.MintsHelper object at 0x7f7529e68cf0>, <psi4.core.BasisSet object at 0x7f752a0820b0>

V_AA, V_AB and V_BA give the expected matrices when referring to V_dimA. Is there an error in how I am initialising MintsHelper?
It seems like the ao_potential() method accepts either 0 or 2 arguments.

6

Old thread, but this question does not seem to have been answered.
Is there a way to directly return the ‘BB’ block of the ao_potential integrals (i.e. nuclear attraction integrals for basis functions on MonomerB with nuclei of MonomerA)? I too need a simple way to do this.