<jupyter_start><jupyter_text>DFT: GGA & Meta GGA Kernels  TheoryIn density functional theory, we are interested in the energy expression:$$E_{xc} = D_{\mu \nu}^{T}(T_{\mu \nu}  + V_{\mu \nu}) + \frac{1}{2} D_{\mu \nu }^{T} D_{\lambda \sigma}^T (\mu \nu|\lambda \sigma) + E_{xc}[\rho_{\alpha}({\vec{r})}, \rho_{\beta}({\vec{r})}]$$Although thet exchange correlation energy $E_{xc}$ is a functional of the density alone, the dependance on  $\rho_{\sigma}(\vec{r})$ is highly non-local. Because of this, small variations of the densities may cause large variations of the exchange correlation potential $v_{xc}$. Additionally, $v_{xc}$ at a given point $\vec{r}_i$ may be sensitive to changes at very distant points of $\vec{r}_j$. In order to overcome this, both semilocal and nonlocal ingredients must be added to the energy density. As expressed on the grid, these include:The gradient of the density:$$g = |\nabla \rho(\vec(r)) |$$The laplacian of the density:$$l = \nabla^2 \rho (\vec{r})$$And the non-interacting kinetic energy density:$$\tau = \frac{1}{2} \sum_k^{occ.} | \nabla \psi_k (\vec{r}) |^2$$Different functional approximations are defined by what ingredients are required to be created. Here, we concentrate on GGAs and meta-GGAs:$$E_{xc}^{GGA}[\rho] = \int f(\rho, g) \cdot d\vec{r}$$$$E_{xc}^{MGGA}[\rho] = \int f(\rho, g, l, \tau) \cdot d\vec{r}$$Where the integrands are known as the kernel, or the exchange-correlation energy density. It is clear here that, the more sophisticated density functional approximations, the more components are added to the exchange-correlation energy and potential. In practice we need to build a Kohn-Sham matrix:$$F_{\mu \nu}^{\alpha} = H_{\mu \nu} + J_{\mu \nu} + V_{\mu \nu}^{xc}$$Where the last therm is the exchange-correlation contribution $V_{\mu \nu}^{xc}$ that is defined as the functional derivative of the energy with respect to the density. $$V_{xc} = \frac{\partial E_{xc}}{\partial D_{ab}} $$Once we have the $V_{xc}$ we can build the Kohn-Sham matrix and solve self consistently.<jupyter_code>import psi4
import numpy as np

mol = psi4.geometry("""
He

symmetry c1
""")
psi4.set_options({'BASIS':               'cc-pVDZ',
                  'DFT_SPHERICAL_POINTS': 6,
                  'DFT_RADIAL_POINTS':    5})

gga_energy, gga_wfn = psi4.energy("PBE", return_wfn=True)
meta_energy, meta_wfn = psi4.energy("TPSS", return_wfn=True)

Vpot = gga_wfn.V_potential()
Vpot_meta = meta_wfn.V_potential()<jupyter_output><empty_output><jupyter_text>Building the GGA kernelFrom the LDA tutorial, we have seen how to obtain the density from the basis functions:$$\rho_{\sigma}(\vec{r}) = D_{\mu \nu}^{\sigma} \phi_{\mu}(\vec{r}) \phi_{\nu}(\vec{r}) $$The GGA depends on the gradient of the density. Using a basis set it is calculated as:$$\nabla_{\sigma} \rho{} (\vec{r}) = 2  D_{\mu \nu}^{\sigma}  \phi_{\mu}(\vec{r}) \nabla \phi_{\nu}(\vec{r})$$ So that we can produce $\gamma$:$$\gamma_{\alpha \alpha}(\vec{r}) = \nabla \rho_{\alpha}(\vec{r}) \cdot \nabla \rho_{\alpha}(\vec{r})$$We then need to get the energy from the kernel by doing a numerical integration:$$e_{\rm{xc}} = w^a f^a_{\rm{xc}}$$Where the $w_{\alpha}$ correspond to the combined Truetler and Lebedev weights at each point needen for the numerical quadrature.Finally, the potential on the grid will be given by the derivative of the kernel with respect to gamma:$$V^{\gamma} = 2 \frac{\partial f}{\partial \gamma_{\alpha\alpha}} \nabla \rho_{\alpha} + \frac{\partial f}{\partial \gamma_{\alpha\beta}} \nabla \rho_{\beta}$$And it can be added and then it needs to get reintegrated back to atomic orbital space:$$V_{ab}^{\gamma} = \int_{\mathbb{R}^3} V^{\gamma} \nabla (\phi_{\mu} \phi_{\nu}) d\vec{r}$$The next calculation assumes that the density matrix $D$ is symmetric. This means that $ \nabla \phi(\vec{r}) D_{\mu \nu} \phi(\vec{r})= \phi(\vec{r}) D_{\mu \nu} \nabla \phi(\vec{r})$.One then ought to be careful with systems where this condition is not met, for example CPHF. Specifically,  here we this matrix contribution as:$$V_{ab}^{\gamma} = 4 \cdot \nabla \phi_{\mu}(\vec{r}) \cdot V_{\alpha \alpha}^{\gamma} \cdot \nabla n_{p}(\vec{r}) \cdot w_{\alpha} \cdot \phi_{\nu} (\vec{r})$$<jupyter_code>#GGA Kernel

D = np.array(gga_wfn.Da())

V = np.zeros_like(D)
xc_e = 0.0

rho = []
points_func = Vpot.properties()[0]
superfunc = Vpot.functional()

# Loop over the blocks
for b in range(Vpot.nblocks()):
    
    # Obtain block information
    block = Vpot.get_block(b)
    points_func.compute_points(block)
    npoints = block.npoints()
    lpos = np.array(block.functions_local_to_global())
    
    # Obtain the grid weight
    w = np.array(block.w())

    # Compute phi!
    phi = np.array(points_func.basis_values()["PHI"])[:npoints, :lpos.shape[0]]
    
    
    phi_x = np.array(points_func.basis_values()["PHI_X"])[:npoints, :lpos.shape[0]]
    phi_y = np.array(points_func.basis_values()["PHI_Y"])[:npoints, :lpos.shape[0]]
    phi_z = np.array(points_func.basis_values()["PHI_Z"])[:npoints, :lpos.shape[0]]
    
    # Build a local slice of D
    lD = D[(lpos[:, None], lpos)]
    
    # Compute rho
    rho = 2.0 * np.einsum('pm,mn,pn->p', phi, lD, phi, optimize=True)
    
    # 2.0 for Px D P + P D Px, 2.0 for non-spin Density
    rho_x = 4.0 * np.einsum('pm,mn,pn->p', phi, lD, phi_x, optimize=True)
    rho_y = 4.0 * np.einsum('pm,mn,pn->p', phi, lD, phi_y, optimize=True)
    rho_z = 4.0 * np.einsum('pm,mn,pn->p', phi, lD, phi_z, optimize=True)
    gamma = rho_x ** 2 + rho_y ** 2 + rho_z ** 2
    
    inp = {}
    inp["RHO_A"] = psi4.core.Vector.from_array(rho)
    inp["GAMMA_AA"] = psi4.core.Vector.from_array(gamma)
    
    # Compute the kernel
    ret = superfunc.compute_functional(inp, -1)
    
    # Compute the XC energy
    vk = np.array(ret["V"])[:npoints]
    xc_e += np.einsum('a,a->', w, vk, optimize=True)
        
    # Compute the XC derivative.
    v_rho_a = np.array(ret["V_RHO_A"])[:npoints]
    Vtmp = 0.5 * np.einsum('pb,p,p,pa->ab', phi, v_rho_a, w, phi, optimize=True)

    #Comute gamma and its associated potential
    v_gamma_aa = np.array(ret["V_GAMMA_AA"])[:npoints]
    Vtmp += 2.0 *np.einsum('pb,p,p,p,pa->ab', phi_x, v_gamma_aa, rho_x, w, phi, optimize=True)
    Vtmp += 2.0 *np.einsum('pb,p,p,p,pa->ab', phi_y, v_gamma_aa, rho_y, w, phi, optimize=True)
    Vtmp += 2.0 *np.einsum('pb,p,p,p,pa->ab', phi_z, v_gamma_aa, rho_z, w, phi, optimize=True)


    # Add the temporary back to the larger array by indexing, ensure it is symmetric
    V[(lpos[:, None], lpos)] += (Vtmp + Vtmp.T)


print("XC Energy %16.8f" % xc_e)
print("V matrix:")
print(V)

print("\nMatches Psi4 V: %s"% np.allclose(V, gga_wfn.Va()))<jupyter_output>XC Energy      -1.03630199
V matrix:
[[-8.73472973e-01 -4.17611601e-01 -7.59172011e-18  0.00000000e+00
  -4.38186681e-18]
 [-4.17611601e-01 -4.02149673e-01 -1.66608209e-18  0.00000000e+00
   2.51135617e-18]
 [-7.59172011e-18 -1.66608209e-18 -5.28713482e-01  9.62964972e-35
   1.46806936e-35]
 [ 0.00000000e+00  0.00000000e+00  9.62964972e-35 -5.28713482e-01
   4.11732573e-36]
 [-4.38186681e-18  2.51135617e-18  1.46806936e-35  4.11732573e-36
  -5.28713482e-01]]

Matches Psi4 V: True<jupyter_text>Building the meta-GGA kernelJust like we did with GGA, meta-GGA requires an extra component to be added. In this case is the kinetic energy density. $$\tau_{\sigma} (\vec{r}) = D_{\mu \nu}^{\sigma} \nabla \phi_{\mu}(\vec{r}) \nabla \phi_{\nu}(\vec{r})$$We calculate the $E_{xc}$ again with $ w^a f^a_{\rm{xc}}$.And finally, the $\tau$ potential contribution can be calculated as:$$V^{\tau} = \frac{\partial f}{\partial \tau}$$Which is expressed as a matrix in the basis sest as:$$V_{\mu \nu}^{\tau} = \int_{\mathbb{R}^3} V^{\tau} \nabla \phi_{\mu} \nabla \phi_{\nu} d\vec{r}$$In the code we calculate this contribution like so:$$V_{\mu \nu}^{\tau} = \frac{1}{2}  \cdot \nabla \phi_{\mu} \cdot V^{\tau}_{a} \cdot w_{a} \cdot \nabla \phi_{\nu} $$<jupyter_code>#meta-GGA Kernel

D = np.array(meta_wfn.Da())

V = np.zeros_like(D)
xc_e = 0.0


rho = []
points_func = Vpot_meta.properties()[0]
superfunc = Vpot_meta.functional()

# Loop over the blocks
for b in range(Vpot.nblocks()):
    

    
    # Obtain block information
    block = Vpot.get_block(b)
    points_func.compute_points(block)
    npoints = block.npoints()
    lpos = np.array(block.functions_local_to_global())
    
    tau = np.zeros(npoints)
    
    
    # Obtain the grid weight
    w = np.array(block.w())

    # Compute phi!
    phi = np.array(points_func.basis_values()["PHI"])[:npoints, :lpos.shape[0]]
    
    
    phi_x = np.array(points_func.basis_values()["PHI_X"])[:npoints, :lpos.shape[0]]
    phi_y = np.array(points_func.basis_values()["PHI_Y"])[:npoints, :lpos.shape[0]]
    phi_z = np.array(points_func.basis_values()["PHI_Z"])[:npoints, :lpos.shape[0]]
    
    # Build a local slice of D
    lD = D[(lpos[:, None], lpos)]
    
    # Compute rho
    rho = 2.0 * np.einsum('pm,mn,pn->p', phi, lD, phi, optimize=True)
    
    # 2.0 for Px D P + P D Px, 2.0 for non-spin Density
    rho_x = 4.0 * np.einsum('pm,mn,pn->p', phi, lD, phi_x, optimize=True)
    rho_y = 4.0 * np.einsum('pm,mn,pn->p', phi, lD, phi_y, optimize=True)
    rho_z = 4.0 * np.einsum('pm,mn,pn->p', phi, lD, phi_z, optimize=True)
    gamma = rho_x ** 2 + rho_y ** 2 + rho_z ** 2
    
    #Compute Tau
    tau = np.einsum('pm, mn, pn->p', phi_x,lD, phi_x, optimize=True)
    tau += np.einsum('pm, mn, pn->p', phi_y,lD, phi_y, optimize=True)
    tau += np.einsum('pm, mn, pn->p', phi_z,lD, phi_z, optimize=True)
        
    
    inp = {}
    inp["RHO_A"] = psi4.core.Vector.from_array(rho)
    inp["GAMMA_AA"] = psi4.core.Vector.from_array(gamma)
    inp["TAU_A"]= psi4.core.Vector.from_array(tau)
    
    # Compute the kernel
    ret = superfunc.compute_functional(inp, -1)
    
    # Compute the XC energy
    vk = np.array(ret["V"])[:npoints]
    xc_e += np.einsum('a,a->', w, vk, optimize=True)
        
    # Compute the XC derivative.
    v_rho_a = np.array(ret["V_RHO_A"])[:npoints]
    Vtmp = 0.5 * np.einsum('pb,p,p,pa->ab', phi, v_rho_a, w, phi, optimize=True)

    
    #Compute gamma and its potential matrix
    v_gamma_aa = np.array(ret["V_GAMMA_AA"])[:npoints]
    Vtmp += 2.0 *np.einsum('pb,p,p,p,pa->ab', phi_x, v_gamma_aa, rho_x, w, phi, optimize=True)
    Vtmp += 2.0 *np.einsum('pb,p,p,p,pa->ab', phi_y, v_gamma_aa, rho_y, w, phi, optimize=True)
    Vtmp += 2.0 *np.einsum('pb,p,p,p,pa->ab', phi_z, v_gamma_aa, rho_z, w, phi, optimize=True)
    
    #Compute V_Tau
    v_tau_a = np.array(ret["V_TAU_A"])[:npoints]
    Vtmp += 0.5 * np.einsum( 'pb, p, p, pa -> ab' , phi_x, v_tau_a, w, phi_x, optimize=True)
    Vtmp += 0.5 * np.einsum( 'pb, p, p, pa -> ab' , phi_y, v_tau_a, w, phi_y, optimize=True)
    Vtmp += 0.5 * np.einsum( 'pb, p, p, pa -> ab' , phi_z, v_tau_a, w, phi_z, optimize=True)

    # Add the temporary back to the larger array by indexing, ensure it is symmetric
    V[(lpos[:, None], lpos)] += (Vtmp + Vtmp.T)


print("XC Energy %16.8f" % xc_e)
print("V matrix:")
print(V)

print("\nMatches Psi4 V: %s"% np.allclose(V, meta_wfn.Va()))<jupyter_output>XC Energy      -1.05818563
V matrix:
[[-8.96867572e-01 -4.49476052e-01 -5.47236918e-17 -5.55111512e-17
  -1.12362033e-18]
 [-4.49476052e-01 -4.09315374e-01 -2.79211769e-17 -1.38777878e-17
  -2.19764937e-18]
 [-5.47236918e-17 -2.79211769e-17 -5.52174114e-01 -1.63704045e-33
   2.69345334e-35]
 [-5.55111512e-17 -1.38777878e-17 -1.63704045e-33 -5.52174114e-01
   3.40730391e-35]
 [-1.12362033e-18 -2.19764937e-18  2.69345334e-35  3.40730391e-35
  -5.52174114e-01]]

Matches Psi4 V: True
{'fraction_non_alphanumeric': 0.17842395044735032, 'fraction_numerical': 0.057811424638678596, 'mean_word_length': 3.613095238095238, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '591089', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 5110, 'n_tokens_neox': 4639, 'n_words': 1143, 'path': 'Example/Psi4Numpy/04-DensityFunctionalTheory/4c_GGA_and_Meta_GGA.ipynb', 'repo_name': 'yychuang/109-2-compchem-lite', 'size': '11624'}