Debugging

Run your test_atoms.c for hydrogen (which has 1 1s electron) with (a) the VWN form for the exchange-correlation functional, (b) N=4000 intervals, (c) the change of variables $r(x)=1/(1-x)-1-x-x^2-x^3$, and (d) $\alpha=0.25$. After about about 50 iterations, you should converge to a value of about -0.445671. Compare your result with the ``official'' result at

http://math.nist.gov/DFTdata/atomdata/tables/01H.html

(Note that we are using the local density approximation (LDA) to density functional theory.)

Hint: The following code fragment gives an example of a suitable initialization of our arrays for hydrogen:

  Z=1.; /* Charge on nucleus */
  lmax=0; /* Only s states are filled in H */

  nmax=ivector(0,lmax); /* Make space for nmax[] */
  nmax[0]=0; /* For s states, we only need the first, 1s, with zero nodes */

  /* Find nmaxmax=max{nmax} */
  nmaxmax=0;
  for (l=0; l<=lmax; l++)
    if (nmax[l]>nmaxmax) nmaxmax=nmax[l];

  F=dmatrix(0,lmax,0,nmaxmax);
  F[0][0]=1.; /* There is only one electron in the 1s for H */