``Geatoms.m''

To calculate the electronic structure of the eight atom cell of germanium, make a copy of ``Hatoms.m'' called ``Geatoms.m''. The eight atom cell of germanium has 32 electrons distributed among 16 states; thus, modify ``Geatoms.m'' so that gbl_f=2 and Ns=16.

Next, we need a Fourier-transformed ionic potential for Vps. The Starkloff-Joannopoulos local pseudopotential for germanium is

$$V_{ps}(r)=-\frac{Z}{r} \frac{1-e^{-\lambda r}}{1+e^{-\lambda(r-r_c)}},$$

where $Z \equiv 4$, $r_c\equiv 1.052$ bohr, and $\lambda\equiv 18.5$ bohr$^{-1}$. Arias's (unfamous) analytic result for the Fourier transform of this potential is the 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}{G^2}+\frac{4\pi Z (1+e^{-\lambda r_c})}{G^2}\left(... ...^{-\pi G/\lambda} \cos(G r_c) \frac{G}{\lambda}}{1-e^{-2\pi G/\lambda}} \right.$$ (1)

$$\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 to use for the case $\vec G=0$ (the first element in the FFT box),

$$\lim_{G \rightarrow 0} \left( \frac{4 \pi Z}{G^2}+ \int e^{-... ...}^\infty{(-1)^n \frac{e^{-\lambda r_c n}}{n^2}}\right)\right)$$ (2)

To avoid needless heartache, feel free to replace the original computation of Vps for the hydrogen atom,

Vps=-4*pi*Z./G2; Vps(1)=0.;

by cutting and pasting the octave form for the resulting formula below directly into ``Geatoms.m'',

%# Ge 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)));

Finally, to allow for viewing of the total charge density, add the code block below to the very end of ``Geatoms.m'' (deleting any ``pause'' statements that may remain in your code or statements for viewing individual wave functions).

%# Compute density for final viewing
Y=W*inv(sqrtm(W'*O(W))); %# Orthonormal wave functions
n=getn(Y,gbl_f); %# Charge density

%# Basic slice from ``100'' edge of cell
sl=reshape(n(1:S(1)*S(2)),S(1),S(2));
%# Make and view image
ppm("100.ppm",sl*0.3,sl,sl); system("display 100.ppm &");

%# 110 slice cutting bonds (assumes cube)
sl=reshape(n(find(M(:,2)==M(:,3))),S(1),S(2));
%# Expand by 2, drop data to restore (approximate) aspect ratio
li=find(rem([1:size(sl,1)],3)!=0); sl=sl(li,:);
%# Make and view image
ppm("110.ppm",sl,sl,sl*0.3); system("display 110.ppm &");