Implementation

First, make a copy of ``$\sim$/A3/dft.m'' called ``$\sim$/A3/Hatoms.m''. Then, remove from ``Hatoms.m'' the code block which displays the wave functions, and replace the code block that relaxes the wave functions with

format long

%# Converge
W=sd(W,20); %# 20 iterations of simple sd() to get nearer to the minimum
W=W*inv(sqrtm(W'*O(W))); %# Restart as orthonormal functions
[W, Elist]=pccg(W,50,1); %# 50 iterations of pclm from same W

At this point, you may wish to rerun ``setup.m'' and ``Hatoms.m'' to verify that the code functions properly as before.

Next, delete the computation of V and Vdual in ``Hatoms.m'', and replace it with the following. First, compute the ionic potential in Fourier space as

%# H atoms
Vps=-4*pi*Z./G2; Vps(1)=0.;

Note that this code is careful to set the infinite $\vec G=0$ component of the ionic potential to zero, a term which exactly balances with the $\vec G=0$ components of the Ewald and Hartree energies. Also, we are calling this Fourier-transformed potential from a single ion ``Vps'' because it will become the Fourier transform of something called a ``pseudopotential'' when we do solid germanium. After computing Vps, next evaluate $\tilde V$ according to Eqs. (4-5) as

Vdual=cJ(Vps.*Sf);

Finally, be sure to store the result in gbl_Vdual for communication to the various functions.