Previously, you solved the Kohn-Sham equations in the very simple
potential coming from a simple harmonic oscillator. For most
electronic structure calculations, however, the potential comes from
the interaction between electrons and nuclei with the form
Our approach to funding the potential Vdual from the nuclei for our periodic calculations is based on our standard approach to solving Poisson's equation,
phinuc=cI(Linv((-4*pi)*O(cJ(rhonuc))));and our the definition of Vdual from Problem Set #8,
Vdual=cJdag(O(cJ(-phinuc)));where the extra `-' in front of phinuc corresponds to the minus sign in Eq. (1). Combining the above two results gives
Vdual=cJdag(O(Linv((4*pi)*O(cJ(rhonuc)))));where we have used the fact that
The only difficulty with this approach is that rhonuc represents a
point charge, whose density is zero everywhere but at a single point
where the density is infinite. We can, however, find the expansion
coefficients
cJ(rhonuc)
mathematically. To do this, we note that

In the case of a single nucleus,
is a point charge
centered at location
with total integral
. Thus, we
have
rhohat=ZI/det(R)*exp(-i*G*rI);where ZI is the charge of the nucleus and rI is a column vector containing its location. Finally, for the case of multiple nuclei, one should sum a series of such contributions, one for each nucleus.