I am trying to compute integrals of the form $(pq|rs)$ where the orbitals $p$ and $q$ are orthogonalized atomic orbitals for a fragment $A$ and $r,s$ are orthogonalized on fragment $B$. I was thinking of a setup similar to this
from pyscf import gto, scf, dft, ao2mo, lo
import numpy as np
#BASIS SET: (4s) -> [2s]
mol1 = gto.M(atom='He 0 0 0', basis='6-31g')
#Lowdin orthogonalization of the ao basis
C1=lo.orth.orth_ao(mol1,method='lowdin')
mol2 = gto.M(atom='He 2 0 0', basis='6-31g')
#Lowdin orthogonalization of the ao basis
C2=lo.orth.orth_ao(mol2,method='lowdin')
This snippet will create two separate molecule instances with their corresponding $C$ orthogonalization matrix (the same in this situation).
Is there a way to either take these two molecule instances and calculate the two-electron integrals between them? Or do I create a "supramolecular" system instead, compute all two-electron integrals, and somehow sort them with respect to the centers of the AOs? If the latter, how to do this?
