``ewald.m''

Make a copy of ``poisson.m'' called ``ewald.m'' and modify ``ewald.m'' as follows.

  1. Modify the calculation of the first Gaussian g1 so that (a) its width is 0.25 bohr instead of 0.75 bohr and (b) its norm is Z instead of unity. (To avoid real-space truncation effects at the cell edges, g1 should remain at the center of the cell.)
  2. Remove all references to g2.
  3. The total charge should be evaluated as
    %# Use structure factor to create all atoms
    n=cI(cJ(g1).*Sf); n=real(n);
    
    rather than g2-g1. (The statement n=real(n); simply removes tiny imaginary parts due to rounding.)
  4. Keep the printouts which check the normalization of the Gaussian g1 and the total charge n in the cell.
  5. Replace the computation of Uanal with the following statement that computes the self-energies of the Gaussians,
    Uself=Z^2/(2*sqrt(pi))*(1/sigma1)*size(X,1);
    
  6. The final printout should give Unum-Uself as the final result for the Ewald energy.

Debugging: The visualizations should show a ``dimer'' oriented along the x direction, with one ``bump'' in the n1-slice and two in the n2- and n3- slices. Finally, the Ewald energy should evaluate to approximately -0.333 hartree. Once you have completed this quick verification of ``ewald.m'', you should comment out the two visualization code blocks as they dramatically slow the run time for larger calculations.

As final confirmation, you should find that at a separation of 4.00 bohr along the x-direction and Gaussian width sigma1=0.5, you get results which are the same to within about 0.025 millihartree for setup parameters S=[32; 32; 32] and R=diag([16; 16; 16]) as you do for setup parameters S=[64; 64; 64] and R=diag([16; 16; 16]). (Be sure to modify and rerun ``setup.m'' and ``ewald.m'' to change these parameters and to make sure that the changes take effect!) More importantly, your result should change by only about 1 millihartree when you change to sigma1=0.25 and evaluate at S=[64; 64; 64], R=diag([16; 16; 16]). All of these results should be about -0.094 hartree.

To ensure that everyone works with the same parameters for future problems, be sure, after running these tests, to restore S and R in ``setup.m'' to their original values of S=[20; 25; 30], R=diag([6 6 6]), and sigma1 in ``ewald.m'' to its original value of sigma1=0.25.

Tomas Arias 2006-04-05