I have been playing around with the Psi4 python interface (in particular, MintsHelper), assessing whether I can integrate it more smoothly into a new fragment-based method I have been working on.
It seems like what I need is to be able to get the ao_potential integrals broken out by nucleus. That is to say that, if the usual matrix contains elements
U_pq = \sum_N <p| Z_N(r_N)^-1 |q>
where r_N is the distance from the electron to nucleus N (and Z_N is its charge), then what I want are matrices with elements
U_Npq = <p| Z_N(r_N)^-1 |q>
such that these matrices would sum up to the usual nuclear potential integrals.
Is this possible?
Perhaps another way of phrasing the question in broader strokes, which might open additional paths to a solution would be: Is it possible to get the integrals corresponding to the attraction of electrons on one subsystem to the nuclei on a different subsystem (and excluding the attraction to its own nuclei)? If the capability described above exists, then clearly I can build this higher-level capability, but maybe something like this is already built in.
I have so far been using horrible hacks involving multiple calls with different atoms “ghosted.”
Have you found out how to obtain these nuclear potential integrals with Psi4NumPy?
If anyone else knows if this is possible, I would also be interested in knowing how to directly get these integrals between the nuclei of one system and the electrons of the other, to calculate the electrostatic energy in monomer-centered basis SAPT.
Sorry, I still don’t know how to do it cleanly. I have a kludge where I do it by subtraction (using 1e integrals matrices calculated for individual atoms and then pairs of atoms). If you want more detail on that to get you by for now I can elaborate . . . might take some time since I am teaching over the summer (though it is not too involved theoretically).
I just had a quick look at my code. It would not make any sense to send it to you because you would need to read 1000 other lines to understand the 20 that are relevant, but here is roughly what I did. I am simplifying and specializing to the case where individual atoms are fragments, but you can abstract the idea and apply it more generally.
loop over atoms a in system:
geometry_a <- build psi4 input (i.e., object for python interface) where all atoms in system except a are ghost atoms
integrals[a] <- invoke psi4 to build all 1e- integrals for geometry_a
When done you will have built an array integrals, where each element contains both the kinetic energy integrals and the nuclear attraction integrals where only one nucleus is “on”. You can do a couple of sanity checks. For example, the kinetic energy integrals should be the same each time (this is not meant to be efficient, just get the job done for experimental code). Secondly, if you add up all of the nuclear attraction integrals, they should give you the same answer as if you had just done one psi4 invocation with no ghosted atoms.
Sorry I misled you about doing things by subtraction. I think I tried something like that but deprecated it in favor of this approach.
Hope that helps.
BTW, do I get the pleasure of knowing who I am helping out?
I missed this thread when it was originally posted, but just saw @Lenora’s followup. Here’s a simple example that will hopefully help anybody needing this functionality in the future (assuming I correctly interpreted the question):
import psi4
import numpy as np
mol = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
""")
e, wfn = psi4.energy('hf/dz', return_wfn=True)
mints = psi4.core.MintsHelper(wfn.basisset())
# Get the density matrix in the AO basis
D = wfn.Da_subset('AO')
nuc_energy = D.vector_dot(mints.ao_potential())
individual_terms = []
for atom in range(mol.natom()):
origin = [mol.x(atom), mol.y(atom), mol.z(atom)]
potential_ints = mints.ao_multipole_potential(origin)[0]
charge = mol.Z(atom)
individual_terms.append(-charge*D.vector_dot(potential_ints))
print(f'individual terms {individual_terms}')
print(f'nuclear energy from full integrals: {nuc_energy}')
print(f'nuclear energy from individual terms: {np.sum(individual_terms)}')
Happy to see the activity on this thread. Glad my hack worked for you Lenora.
Thanks for the posts andysim and jmisiewicz. I look forward to trying this when I get back to this project later this summer.
I will be wanting the full matrix of integrals in the AO basis, but it looks like the adaptation I would need to the code above just involves setting up MintsHelper with the AOs (and then not building and contracting with a density matrix).
. I am also interested in printing out individual Ven.