How to obtain correct HF energy from matrices

Hi all,

how do I obtain manually the energy from the matrices H, J, K of a HF calculation. I use

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()

From theory we know the energy without nucleus Coulomb part is

E =  trace(C.T * H * C) + 0.5 * (trace(C.T * J * C) - trace(C.T * K * C))

Somehow, this result is not identical to the HF energy (without nucleus Coulomb).
How do I get the resulting energy of a HF calculation?

Thanks in advance,
Wiktor

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

Thanks. This was one source of error.

The second source of error is that I calculate the nuclear repulsion energy manually and wrong.
Is the nuclear repulsion energy not just

.5 * sum(Zi * Zj / r_ij) forall i != j

What am I missing?

Not sure. Maybe your r_ij is in Angstrom and not in bohr?
Just guessing, obviously.

That was exactly it, thanks.

Just to be clear, the default units of Psi4 are Angstrom, but the nuclear repulsion energy has to be computed using Bohr?

electronic structure methods normally always use atomic units (bohr, Hartree, etc.) internally and you need to match the units of all parts.

The “default” unit for lengths internally will be bohr, but Angstrom is more convenient for the user interface.

wfn.molecule().nuclear_repulsion_energy() will give you back the nuclear energy in Hartree.

The other thing to note is that the orbitals are one iteration in advance of the current J/K matrices so your energy will differ slightly (usually within the current delta E).

Instead of np.trace(np.dot(D , H)) you can call np.vdot(D, H) for an inner product without trace requirements.

I ran the input provided by hokru above and the energies did not match exactly.
The printout is:
reference : -75.9607980647
HF energy : -75.9607971762
The output file matches the reference energy:
Total Energy = -75.9607980647016205

The discrepancy goes down when I crank up the convergence, but it’s still about four orders of magnitude greater than the current delta E.

I won’t bother tracking the error down further.

As noted above the J and K matrices are for the second to last iteration. We always diagonalize the Fock matrices and update the orbital during the last iteration, but do not recompute the J/K matrices. So there will always be a different in the example Hokru gave.

The only way to compute the exact same energy is to recompute the J/K matrices from the current Density matrix.

@jmisiewicz 4-orders of magnitude is somewhat concerning. The example by @bge6 is about 1.e-7 Hartree which is expected with default convergence metrics.

@dgasmith

On adding set e_convergence 12 to Holger’s script, I see

reference :  -76.02665366188839
HF energy :  -76.02665362881304

The error is in the eighth decimal place.