In Psi4/Numpy it is possible to easily compute a potential (I assume xc remainder potential) for an arbitrary density matrix using VBase, per the example here.

Is there a similarly easy way to calculate the xc energy?

EDIT: If there is a similar way to obtain the total energy of a HF, DFT or hybrid-DFT approximation at arbitrary density matrix that would be as good or better!

Obviously I can do this using blocks, a grid, and direct calls to XC routines. But the internal C++ code will be faster and (hopefully) easier.

import psi4
# Obtain a SCF Wavefunction
mol = psi4.geometry("He\nsymmetry c1")
scf_e, scf_wfn = psi4.energy("PBE/sto-3g", return_wfn=True)
# Grab a few matrices
D = scf_wfn.Da().clone()
V_out = scf_wfn.Da().clone()
V_out.np[:] = 0
Vpotential = scf_wfn.V_potential()
# Set arbitrary density
D.np[:] = 1
# Compute V which evaluates the energy as well
Vpotential.set_D([D])
Vpotential.compute_V([V_out])
print(V_out.np)
# Actual values exist on this dictionary
print(Vpotential.quadrature_values()["FUNCTIONAL"])

Building out a Vpotential can be done without the wavefunction, but it’s slightly tricky to do so. In principle the signature is here, but (and I am recalling) there is some sort of issue:

Perfect. Can indeed confirm that seems to work for obtaining the energy with varying density matrix (which was my main goal).

In future I’ll see if I can work out a sneaky way to use it for changing the functional itself (e.g. for density-driven / functional-driven analysis and density corrected DFT). But that’s lower priority for me right now.