How to obtain correct HF energy from matrices

2

You have to use the expression for spatial orbitals. This would do the job (excuse the poor python):

import psi4
import numpy as np

h2o = psi4.geometry("""
symmetry c1
O
H 1 0.96
H 1 0.96 2 104.5
""")

set {
  save_jk true
}

eng, wfn = psi4.energy("scf/def2-svp", molecule=h2o, return_wfn=True)
Enuc = h2o.nuclear_repulsion_energy()

j = wfn.jk().J()[0].to_array()
k = wfn.jk().K()[0].to_array()
h = wfn.H().to_array()
c = wfn.Ca_subset("AO", "OCC").to_array()

J=np.array(j)
K=np.array(k)
C=np.array(c)
H=np.array(h)


D=np.dot(C,C.transpose())

Eel =  2*np.trace(np.dot(D , H)) +  2*np.trace(np.dot(D , J)) - np.trace(np.dot(D , K))

print 'reference : ',eng
print 'HF energy : ', Eel+Enuc

I think you would enjoy this: https://github.com/psi4/psi4numpy