January 16, 2018, 10:46am
how do I obtain manually the energy from the matrices H, J, K of a HF calculation. I use
J = wfn.jk().J().to_array()
K = wfn.jk().K().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,
January 16, 2018, 1:28pm
You have to use the expression for spatial orbitals. This would do the job (excuse the poor python):
import numpy as np
h2o = psi4.geometry("""
H 1 0.96
H 1 0.96 2 104.5
eng, wfn = psi4.energy("scf/def2-svp", molecule=h2o, return_wfn=True)
Enuc = h2o.nuclear_repulsion_energy()
j = wfn.jk().J().to_array()
k = wfn.jk().K().to_array()
h = wfn.H().to_array()
c = wfn.Ca_subset("AO", "OCC").to_array()
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:
January 16, 2018, 7:29pm
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?
January 16, 2018, 8:10pm
Not sure. Maybe your r_ij is in Angstrom and not in bohr?
Just guessing, obviously.
January 16, 2018, 8:23pm
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?
January 16, 2018, 8:40pm
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).
np.trace(np.dot(D , H)) you can call
np.vdot(D, H) for an inner product without trace requirements.
July 24, 2019, 5:10pm
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.
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.