``Ge.m''

Begin by making a new copy of ``dft.m'' called ``Ge.m'' and confirming that ``setup.m'' is currently set to run with a problem size of S=[32; 32; 32]; R=diag([8 8 8]). Because ``dft.m'' was already running with fillings f=[2; 2/3; 2/3; 2/3], the only change needed to calculate Ge is to replace the simple harmonic oscillator potential with the pseudopotential.

The final form used for the potential in the code is the vector of dual coefficients, Vdual=cJdag(O(cJ(V))), where V is the vector of sample values of the potential. For the pseudopotential, our strategy will be to compute the expansion coefficients $\hat V_{ps}$=cJ(V) directly by analytically Fourier transforming the pseudopotential in Eq. (1).

Keeping in mind that we have a periodic copy of the atom at all positions ${\bf R}\vec m$, where $\vec m$ is any vector of three integers, the actual total pseudopotential is $\sum_{\vec m} V_{ps}(\vec r-{\bf R}\vec m)$. Thus, the relation between the expansion coefficients Vps and the pseudopotential is

$$\sum_{\vec G} (\hat V_{ps})_{\vec G}\, e^{i\vec G \cdot \vec r} = \sum_{\vec m} V_{ps}(\vec r-{\bf R} \vec m),$$

which, upon multiplication of both sides by $e^{-i\vec G \cdot \vec r}$ and integration over the supercell, gives

$$(\hat V_{ps})_{\vec G} = \frac{1}{\det{\bf R}} \sum_{\vec m} \int_{cell} e^{-i\vec G \cdot (\vec r-{\bf R} \vec m)} V_{ps}(\vec r-{\bf R}\vec m) \, d^3r$$

$$= \frac{1}{\det{\bf R}} \int_{\mbox{all space}} e^{-i\vec G \cdot \vec x} V_{ps}(\vec x),$$ (2)

where on the first line we use the fact that always $e^{i \vec G \cdot {\bf R} \vec m}=1$ and, in going to the second line, we note that the range of points covered by all $\vec r-{\bf R}\vec m$, where $\vec r$ ranges over one supercell, is in fact all of space. Finally, the integral over all of space can be done almost analytically. Arias's (unfamous) result for this integral is in terms of a rapidly convergent asymptotic series (four terms suffice to give 26 digits!),

$$\int e^{-i\vec G \cdot \vec x} V_{ps}(\vec x) = -\frac{4\pi Z_v}{G^2}+\frac{4\pi Z_v (1+e^{-\lambda r_c})}{G^2}\l... ...^{-\pi G/\lambda} \cos(G r_c) \frac{G}{\lambda}}{1-e^{-2\pi G/\lambda}} \right.$$ (3)

$$\mbox{\ \ } \left. + \sum_{n=0}^\infty (-1)^n \frac{e^{-\lambda r_c n}}{1+\left(\frac{n \lambda}{G}\right)^2} \right)$$

Along with this formula, we need the following limit for the case $\vec G=0$ (the first element in the FFT box),

$$\lim_{G \rightarrow 0} \int e^{-i\vec G \cdot \vec x} V_{ps}... ...}^\infty{(-1)^n \frac{e^{-\lambda r_c n}}{n^2}}\right)\right)$$ (4)

Combining Eqs. (2-4) gives $\hat V_{ps}$, from which we then get the final result needed by the code Vdual=cJdag(O($\hat V_{ps}$)). To avoid needless heartache, feel free to cut and past the octave form for the resulting formula,

%# Pseudopotential
Z=4;
lambda=18.5;
rc=1.052;
Gm=sqrt(G2);
Vps=-2*pi*exp(-pi*Gm/lambda).*cos(Gm*rc).*(Gm/lambda)./(1-exp(-2*pi*Gm/lambda));
for n=0:4
  Vps=Vps+(-1)^n*exp(-lambda*rc*n)./(1+(n*lambda./Gm).^2);
end
Vps=Vps.*4*pi*Z./Gm.^2*(1+exp(-lambda*rc))-4*pi*Z./Gm.^2;

n=[1:4];
Vps(1)=4*pi*Z*(1+exp(-lambda*rc))*(rc^2/2+1/lambda^2* ...
                     (pi^2/6+sum((-1).^n.*exp(-lambda*rc*n)./n.^2)));
Vps=Vps/det(R);
Vdual=cJdag(O(Vps)); %# Dual coefficients
gbl_Vdual=Vdual; %# Export value of global variable