``physics2.c''

 

To complete this assignment we will need to evaluate the interaction energy of nuclei in a periodic crystal. 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 Homework Assignment # 9.

To do this, copy your program physics.c into a new program physics2.c and make the following modifications.

• Of the two Gaussians, make one centered at the origin (as we had done before) and the other centered at a point r=(x,0,0) along the x-axis. You can accomplish this with the following code fragment in the main program along with an new function GaussR2() modeled appropriately on your original GaussR(). For this exercise, you may ``hardwire'' the value of the parameter x into the GaussR2() routine.

                 .
                 .
    rho1.evalr(GaussR,sigma1);
    rho2.evalr(GaussR2,sigma2);
                 .
                 .

• 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, rho.
• 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 U-Uself where Uself .
• You may remove the 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.02 millihartree for (n=32, dx=0.5) as you do for (n=64, dx=0.25). More importantly, your result should change by only about 1 millihartree when you change to sigma1=sigma2=0.25 and evaluate at (n=64, dx=0.25). All of these results should be about -0.094 hartree.

Finally, using the parameters sigma1=sigma2=0.25, and (n=64, dx=0.25), calculate and save in a file for later the values of U-Uself for nuclei at separation x=0.50, 1.00, 1.25, 1.50, 1.75, 2.00, 4.00 and 6.00 bohr.