``test_atoms.c''

Write a main program which attains self-consistency with the algorithm,

• Initialize:
  1. Rho$_0$: Rho[k]=0 for all k
• Iterate:
  1. get phi and Vxc from Rho$_n$
  2. compute V[k]=-1/r[k]+phi[k]+Vxc[k] for all k
  3. get E's from V
  4. get Psi's from E's and V
  5. get Rhonew from Psi's
  6. Rho$_{n+1}$[k]=$(1-\alpha)$ Rho$_n$[k]+ $\alpha$ Rhonew[k] for all k

where $\alpha$ is any value between 0 and 1.

Note that if $\alpha=1$, Rho$_{n+1}$=Rhonew, and we are just putting the output charge density back in as the input, and iterating. This process, however, often is unstable because the electron density Rho changes too quickly. Taking $\alpha < 1$ helps calm this instability by changing Rho more slowly. For instance, for $\alpha=0.25$ (which seems to work in most cases), we are ``mixing'' 75% of the previous estimate of Rho with 25% of the new estimate.

Finally, to monitor the progress of your program, you should have it print out the total energy at each iteration. A suitable code fragment for computing this is

                        .
                        .
      /* Compute and output total energy */

      /* Get correction to sum of electron energies */
      getDepsxc(Depsxc,Rho,r,dr,N);
      for (k=0; k<=N; k++)
        integrand[k]=(-0.5*phi[k]+Depsxc[k])*Rho[k];
      Etot=simpint(integrand,r,dr,N);

      /* Add on the sum of the electron energies times the occupancies */
      for (l=0; l<=lmax; l++)
        for (n=0; n<=nmax[l]; n++)
          Etot+=F[l][n]*E[l][n];
      
      printf("Etot: %20.15f\n",Etot);
                        .
                        .

Note that, because, your solution is converging (hopefully!) and all the quantities eventually stop changing, it doesn't really matter if you use Rho[] or Rhonew[] or if you compute phi[] from Rho[] or Rhonew[] in the above code fragment.