Ab initio calculation of H bond length!

 

Use the following code fragment to initialize Vtilde for a molecule (two nuclei) instead of just one:

  // Initialize for H2 molecule
  pvector Jrhonuc1(i3vec(n,n,n),d3vec(d,d,d));
  Jrhonuc1.evalG(Nuc1,1./Vol);
  pvector Jrhonuc2(i3vec(n,n,n),d3vec(d,d,d));
  Jrhonuc2.evalG(Nuc2,1./Vol);
  Vtilde=Jdag(O(Linv((4*pi)*O(Jrhonuc1+Jrhonuc2))));

Because of your coding for Nuc2, this represents a hydrogen molecule with the nuclei separated by 1.5 bohr. Running at n=32 and spacing dx=0.5 you should find (after about 25 iterations) a total energy for the electrons of about -1.443237 hartree. Adding this to the result for the nuclear energy U in Part 1.1 (0.313589 hartree) gives a total energy for the molecule of E=-1.129648 hartree.

Using this approach, compute the total energy of the H molecule for the separations 0.50, 1.00, 1.25, 1.50, 1.75, 2.00, 4.00 and 6.00 bohr.

Make a plot of your results. From the plot, determine the bond-length of H and compare with the experimental result from a standard reference such as the CRC handbook. Show the plot to the grader to get credit for this problem!