Charge distribution: ``poisson.m''

The charge density for our example solution to Poisson's equation will be

$$n = g_1(r)-g_2(r)$$ (4)

$$\equiv \frac{e^{-r^2/(2\sigma_2^2)}}{(2\pi\sigma_2^2)^{3/2}} -\frac{e^{-r^2/(2\sigma_1^2)}}{(2\pi\sigma_1^2)^{3/2}}$$ (5)

where $r$ is the distance from the center of the cell, $g_1(r)$ and $g_2(r)$ are normalized three-dimensional Gaussian distributions (each containing a net charge of unity), and $\sigma_1=0.75$ bohr and $\sigma_2=0.50$ bohr, respectively.

Copy the code ``poisson.m'' from Appendix C to ``$\sim$/A1'', and complete the program by replacing ``...'' on the line labelled ``CODE INSERTION # 1'' with a code fragment to evaluate in the column vector dr the distance between each sampling point in r and the center of the cell,

sum(R,2)/2

Hints: Consider the matrix

ones(prod(S),1)*sum(R,2)'

and use a similar procedure to what you used to compute G2. Feel free to type ones(prod(S),1)*sum(R,2)' at the command prompt to see the value of this expression. Finally, note that in octave, ``sqrt(A)'' creates an object of the same dimensions as ``A'' and with elements equal to the square root of each corresponding element of ``A''.

To verify your code fragment, run ``poisson.m''. This octave program will use the distances which you compute in dr to evaluate $g_1$, $g_2$ and $n$ according to Eq. (6). The program then checks your distributions by integrating each of these functions by a simple Riemann sum $\int f(\vec r)\,d^3r \approx \sum_i f(r_i) \Delta V$, where the volume per sample point is $\Delta V=\mbox{det\,}R/\prod_{k=1}^3 S_k$. For these tests, your results should be good to at least 3 decimal places. Finally, to aid in debugging, the program will plot the charge density in three perpendicular planes passing through the center of your cell. Each of these should appear as a positive Gaussian with a noticeable negative ripple as you move away from the main peak.