``getg''

double getg(double g[], double Rho[],
	 double Z, int lmax, int nmax[], int nmaxmax, double **F,
	 double r[], double dr[], int N)

Return value:

• Etot: total energy of atom

Output:

• g[k]: at the grid points r[k]

Input:

• Rho[k]: for all grid points r[k]
• Z: nuclear charge
• lmax: maximum angular momentum to consider (l=0..lmax)
• nmax[l]: maximum number of nodes sought for angular momentum l
• nmaxmax: maximum of all nmax[l]
• F[l][n]: shell occupancies ( )
• r[], dr[], N: standard grid information

Using the pieces of your main routine from the Homework #6, construct a subroutine of the above prototype, which takes Rho[] as input, and then computes from it the total potential V=Vnuc+phi+Vxc, the wave functions Psi[][][], and the density corresponding to those wave functions, Rhonew[]. Finally, the routine should put the difference Rhonew[]-Rho[] in g[] and return the value Etot.

Hint: Note that getg() will have to allocate space for many variables such as Vxc[], phi[], Psi[][][]. You will find the Numerical Recipes utilities dvector, dmatrix, d3tensor very useful for this. Now that you are allocating this space in a subroutine, it is essential that you deallocate everything before you exit the routine.

• Debugging