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!