Obtain the energy from DFT at arbitrary density matrix

Psi4Education
3

Similarly you can use the built in C++ routines:

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:

Maybe I’m wrong, that looks like enough to stitch it together from the Python side. I don’t have time to try it out, but let me know if it works!