# How to evaluate Hartree--Fock exchange energy density on grids in pyscf?

I am trying to evaluate the Hartree--Fock exchange energy density on grids in pyscf with an instance of class "UKS", in order to test some functionals depending on the exchange energy density such as Becke05. I have noticed that there is an example in 31-xc_value_on_grid.py, but it seems that the function dft.libxc.eval_xc does not contain the Hartree--Fock exchange energy, that

ao_value = numint.eval_ao(mol, coords, deriv=1)
rho = numint.eval_rho(mol, ao_value, dm, xctype='GGA')
exc, vxc = dft.libxc.eval_xc('.2*HF + .08*SLATER + .72*B88, .81*LYP + .19*VWN', rho)[:2]
print('Exc = %.12f  ref = -7.520014202688' % numpy.einsum('i,i,i->', exc, rho[0], weights))
exc, vxc = dft.libxc.eval_xc('.08*SLATER + .72*B88, .81*LYP + .19*VWN', rho)[:2]
print('Exc = %.12f  ref = -7.520014202688' % numpy.einsum('i,i,i->', exc, rho[0], weights))


Exc = -7.520014211805 ref = -7.520014202688

Exc = -7.520014211805 ref = -7.520014202688

So is there some functions in pyscf to evalue the HF exchange energy density on grids? Thanks in advance for the help!

Exact exchange is not included in the DFT grid since exact exchange is not DFT. PySCF does implement seminumerical schemes for evaluating exchange energies; this is documented in the user guide in https://pyscf.org/user/sgx.html

• Thanks for the help. Actually I was wondering if there are any functions to evaluate the exchange energy on grids, with or without the exchange matrix. SGX seems to be an efficient way to calculate the matrix, but is there any way I can use the matrix to (or directly) calculate the exchange energy density (or exchange potential) on some specific grids with implementations in pyscf? Commented Oct 26, 2022 at 12:07
• @LexLee this is exactly what is done in the SGX code, just look into it for details. Commented Oct 28, 2022 at 14:34

I find this can be calculated by a similar procedure as in pyscf.tools.cubegen.mep, where a delta function at a grid is approximated by a Gaussian with efficiently large exponent. Here is my code:

def get_ExxDen(ks,coords=None):
if(type(coords)==type(None)): coords = ks.grids.coords
mol = ks.mol
dm = ks.make_rdm1()
fakemol = gto.fakemol_for_charges(coords)
ints = df.incore.aux_e2(mol, fakemol)
if dm.ndim == 2:  # RHF DM
dm = numpy.asarray((dm*.5,dm*.5))
ni = ks._numint
ao = ni.eval_ao(mol,coords)
ExxDen = [0,0] # ExxDen for alpha and beta spin
for idm in range(2):
ExxDen[idm] = einsum('pi,ij->pj', ao, dm[idm])
ExxDen[idm] = einsum('pj,jkp->pk', ExxDen[idm], ints)
ExxDen[idm] = einsum('pk,lk->pl', ExxDen[idm], dm[idm])
ExxDen[idm] = einsum('pl,pl->p', ExxDen[idm], ao)
return ExxDen[0],ExxDen[1]

• That is exactly the same approach as used in SGX: one can approximate nuclear attraction integrals <u|V_A|v> used in seminumerical exchange with a two-electron integral <uv|A> with a dummy S-type function on the nucleus A with a very large exponent. Commented Oct 28, 2022 at 14:34