Energy from integrals and density matrices does not match output

1

Hi all,

I have been using Psi4 to extract some integrals and density matrices, which I always use to calculate the total energy of the system I’m working with as a check to make sure I’m using the right quantities. It has always worked, until I came across the calculation for which I include the input below.

molecule mol {
-1 1
N       -0.3622910716      4.2591226882     -0.6798687874
C       -1.8101472287      4.0936094603     -0.7410374944
H       -0.1337808573      5.2168840640     -1.0319642161
H       -2.1439853850      4.2344989128     -1.7932206134
H       -2.2888756047      4.8838755020     -0.1208070642
C       -3.7763324888      2.5830122505     -0.3248365808
C       -1.6067696806      1.5922925607     -1.1080654733
C       -1.8096049953      2.5013519297      1.2331396929
C       -2.2415581448      2.6944048337     -0.2351204659
H       -1.9280514078      0.5856778314     -0.7630076235
H       -1.9121443771      1.7097815626     -2.1702744377
H       -0.4976790064      1.6331714566     -1.0576834118
H       -2.2631015007      3.2826832786      1.8806020633
H       -2.1335733473      1.5067800256      1.6092133372
H       -0.7049774398      2.5622338559      1.3350360325
H       -4.1192290822      2.7129906865     -1.3741914757
H       -4.1204579774      1.5879498596      0.0317189524
H       -4.2642154008      3.3627845355      0.2992976988

symmetry c1
}

set {
  scf_type pk
  save_jk true
  soscf true
}

eng, wfn = psi4.energy('pbe/aug-cc-pvdz', return_wfn=True)
nre = mol.nuclear_repulsion_energy()
ks = VBase.quadrature_values(wfn.V_potential())
E_xc = ks["FUNCTIONAL"]

print(2.0*wfn.Da().vector_dot(wfn.jk().J()[0]))
print(2.0*wfn.Da().vector_dot(wfn.H()))
print(E_xc)
print(nre)

print(2.0*(wfn.Da().vector_dot(wfn.H())+wfn.Da().vector_dot(wfn.jk().J()[0]))+E_xc+nre)

The corresponding energies as found in the output are:

    Nuclear Repulsion Energy =            255.1348092728447909
    One-Electron Energy =                -853.3919961572901229
    Two-Electron Energy =                 384.3724083783367860
    DFT Exchange-Correlation Energy =     -38.2622811292476115
    Empirical Dispersion Energy =           0.0000000000000000
    VV10 Nonlocal Energy =                  0.0000000000000000
    Total Energy =                       -252.1470596353561575

The HF energy converges just fine, and the energies match, as expected. If I run the input above, the KS orbitals do not converge unless I resort to second-order convergence, but the energies I compute using the density matrices and integrals do not match, as can be seen below:

print(nre): 255.134809273 
print(2.0*wfn.Da().vector_dot(wfn.H())): -841.0806467281856
print(2.0*wfn.Da().vector_dot(wfn.jk().J()[0])): 371.99678840941533
print(E_xc): -38.2622550156

With the total energy being -252.2113301751721.

Since the nuclear repulsion energy and the E_xc term do not depend on the density matrices, and seem to be correct to some degree, I guess it may be something wrong with the KS orbitals/density matrices, and since the HF energy does not exhibit this problem, I’m tempted to think this is a consequence of having to turn to second-order convergence. I would like to know if anyone has run into the same problem or knows what the problem is.

Any help is appreciated!

Thanks!

2

Sure, you are in a UKS wavefunction while your density accessors are only considering RKS wavefunction data:

Dt = wfn.Da().clone()
Dt.axpy(1.0, wfn.Db())

J = wfn.jk().J()[0].clone()
J.axpy(1.0, wfn.jk().J()[1])

print("one_el:  %12.10f" % Dt.vector_dot(wfn.H()))
print("two_el:  %12.10f" % (0.5 * Dt.vector_dot(J)))

This should get you patched up. You should notice that the energies will be slightly different because the iterative process isn’t possible to end “cleanly.” After the last Fock build where the energies are printed, the orbitals and densities are updated. The tighter the convergence, the smaller the energies will deviate.

3

Thanks for your reply.

I tried what you suggested, and it did not work. I am having some trouble in understanding why it is a UKS wavefunction, as I am dealing with an anion and have set the multiplicity to one. The reason why it did not work is because there is no wfn.jk().J()[1] stored, which I assume makes sense, because it is a restricted reference.

In trying to find the source of the problem, I tried running the same calculation using the SOSCF converged orbitals as guess, and they not only do not converge in the first iteration, but they take more iterations to converge than the first calculation. If I turn SOSCF off in this second calculation but still use the SOSCF form the first calculation orbitals, it does not even converge. So, my guess is that is not storing the optimal orbitals in wfn.Ca() after the SCF is converged when SOSCF is on.

I have tried the same with water and ethanol, and SOSCF works as expected. At this point, I do not know if this is case-specific or something more general.

Thanks!