Nuclear-nuclear potential energy

An important part of the energy of a molecule or solid is the potential energy between the nuclei

$$U_{nn} = \frac{1}{2} \sum_{I,J} \frac{Z_I Z_J}{\vert\vec r_I-\vec r_J\vert}.$$

For us, this sum is somewhat tricky because we work in periodic boundary conditions and must sum over all of the the repeated copies of the nuclei. However, outside of a small Gaussian distribution of charge, the potential is identical to that of a point charge. Thus, we can calculate the energy of a periodic array of charges by solving for the electrostatic energy of Gaussian distributions as we did in Problem Set # 7.

To do this, make a new copy your program ``poisson.m'' called ``poisson2.m'' and make the following modifications.

• In addition to the current set of distances dr from each grid point to the point (diag(R)'/2), compute a new set of distances dr2 from each grid point to the point (diag(R)'/2+[x,0,0]).
• Because these Gaussian charges are to represent the two nuclei in the hydrogen molecule, be sure to add them together (not subtract) to form the total density n.

• Because we are interested only in the energy of interaction of the two Gaussians and not their self-energy, you should printout the value of Unum-Uself where Uself $\equiv 1/(2\sqrt{\pi})(1/\sigma_1+1/\sigma_2)$.

• You may remove any plotting and analytic tests of phi.

Debugging: You should find that at a separation of x=4.00 and widths sigma1=sigma2=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'' 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=sigma2=0.25 and evaluate at S=[64; 64; 64], R=diag([16; 16; 16]). All of these results should be about -0.094 hartree.

Finally, using the last of the above sets of parameters, (sigma1=sigma2=0.25, S=[64;64;64], R=diag([16; 16; 16])), calculate and save in a file for later use the values of Unum-Uself for nuclei at separations x=0.50, 1.00, 1.25, 1.50, 1.75, 2.00, 4.00 and 6.00 bohr.