Hi, I’m trying to use mintshelper to obtain two electron integrals from hartree fock orbitals and confirm that I can reconstruct the HF energy. I was able to successfully do so for one electron integrals (thanks to online resources made available from psi4!)
I’ve obtained the transformed two electron integrals through using psi4’s implementation from here
I’m using the equation 15 from here
For my trivial water molecule in ground state, first 5 orbitals are occupied so I tried to do so by
en2b = 0.0
for i in range(5):
for j in range(5):
en2b += mo_I[i,i,j,j] - mo_I[i,j,j,i]
However, this does not match the number reported by the SCF procedure.
This is what I did for one electron energy:
T = np.array(mints.ao_kinetic())
V = np.array(mints.ao_potential())
H = T+V
mo_H = np.dot(np.dot(C.T, H), C)
en2a = 0.0
for i in range(5):
en2a += mo_H[i,i]
en2a *= 2.0 # pairs of electrons
This matches the number reported by the SCF procedure.
What am I doing wrong?
This is the full script:
#! 6-31G H2O Test FCI Energy Point
memory 250 mb
molecule h2o {
O .0000000000 .0000000000 -.0742719254
H .0000000000 -1.4949589982 -1.0728640373
H .0000000000 1.4949589982 -1.0728640373
units bohr
symmetry c1
}
set {
basis 6-31G
}
en, wfn = energy('scf', return_wfn=True)
Matrix_C = wfn.Ca()
C = np.asarray(Matrix_C)
nbf = wfn.nso()
print("Number as basis functions: %d" % nbf)
Enuc = h2o.nuclear_repulsion_energy()
print "nuclear repulsion energy:", Enuc
mints = psi4.core.MintsHelper(wfn.basisset())
Matrix_I = mints.ao_eri()
I = np.asarray(Matrix_I)
mo_I = np.asarray(mints.mo_transform(Matrix_I, Matrix_C, Matrix_C, Matrix_C, Matrix_C))
T = np.array(mints.ao_kinetic())
V = np.array(mints.ao_potential())
H = T+V
mo_H = np.dot(np.dot(C.T, H), C)
en2a = 0.0
for i in range(5):
en2a += mo_H[i,i]
en2a *= 2.0 # pairs of electrons
print "one-electron energy:", en2a
en2b = 0.0
for i in range(5):
for j in range(5):
en2b += mo_I[i,i,j,j] - mo_I[i,j,j,i]
print "two-electron energy:", en2b
print "total energy:", Enuc+en2a+en2b