HF two electron energy based on two electron integral

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

You need twice the coulomb energy for the J (Coulomb) type integrals:

en2b += 2.0 * mo_I[i,i,j,j] - mo_I[i,j,j,i]

Also, so that your energy and the Psi energy match exactly you should set:

set scf_type pk

You might want to check out psi4numpy for more examples.

1 Like