``getVxc'' and ``getDepsxc''

getVxc(double Vxc[], double Rho[], double r[], double dr[], int N)

getDepsxc(double Depsxc[], double Rho[], double r[], double dr[], int N)

Output:

• Vxc[k]: $V_{xc}(r[k])$ for all grid points r[k]
• Depsxc[k]: $\Delta\epsilon_{xc}(n(r[k])) \equiv -\epsilon_{xc}'(n(\vec x)) n(\vec x)$ for all grid points r[k]

Input:

• Rho[k]: $4 \pi r^2 n(r[k])$ for all grid points r[k]
• r, dr, N: standard grid information

Write routines of the above prototypes which (a) loop over all grid points k, (b) convert our variable Rho[] to the true value of the density, (c) evaluate rs, (d) call upon excPZ and excpPZ, and (e) evaluate $V_{xc}$ and $\Delta\epsilon_{xc}$ according to Eqs. (1,2). There are two special cases to consider to avoid NaNs corrupting your data. At the origin (k=0), n=Rho[0]/r[0] gives NaN, and, at infinity, (k=N), Rho=0. Thus, if either k=0 or Rho$<$SMALL, your routine should set Vxc=0 and Depsxc=0. (Recall that we #define SMALL 1e-300.)

Hint: In Eqs. (1,2) $\epsilon_{xc}'(n)$ refers to $d\epsilon_{xc}/dn \equiv (d\epsilon_{xc}/dr_s)(dr_s/dn)$, whereas ``excpPZ'' returns $\epsilon_{xc}'(r_s)\equiv d\epsilon_{xc}/d r_s$. You therefore must multiply the return value of ``excpPZ'' by $(dr_s/dn)=(-r_s/(3 n))$ to evaluate $\epsilon_{xc}'(n)$ in Eqs. (1,2).