Debugging

Change the initialization parts of your sch2.C code to now initialize the potential to that of a hydrogen atom at the origin:

  // Initialize for H atom
  pvector Jrhonuc(i3vec(n,n,n),d3vec(d,d,d));
  Jrhonuc.evalG(Nuc1,1./Vol);
  Vtilde=Jdag(O(Linv((4*pi)*O(Jrhonuc))));

Have the program run Itmx=25 iterations of preconditioned conjugate gradients on a grid of size n=32 and spacing dx=0.5. At the end, have the program normalize u and output the result in real space:

  u=u*(1./u.dat[0][0][0]);
  u=u*(1./sqrt(u^O(u)));
  pvector psi=I(u);
  psi.printr("u");

(The first of these statements makes picks the complex phase of the wave function to be real).

Make a comparison plot of the wave function in file ``u'' (first column is r, second column is psi, third column is the imaginary part of psi which should be very close to zero) with the analytic solution for the 1s state, . You should find a quite close match, except for the point right on top of the nucleus (r=0). Show the plot to the grader to get credit for this problem!

Note: I find an energy of about -0.321. You should not worry that this isn't -0.5 because Poisson's equation only gives us the potential (and the the electron's energy) up to an arbitrary constant.