Electronic structure of Ge

The electronic configuration of the Ge atom ($Z=32$) is [1s$^{2}$2s$^{2}$2p$^{6}$3s$^{2}$3p$^{6}$3d$^{10}$]4s$^{2}$4p$^{2}$. Verify that your atomic solver functions properly by running it with Itmx=20, N=4000, and the following specification for the above electronic configuration.

  /* Specs for Ge */
 Z=32.;
 lmax=2;

 nmax=ivector(0,lmax);
 nmax[0]=3;
 nmax[1]=2;
 nmax[2]=0;

 nmaxmax=0;
 for (l=0; l<=lmax; l++)
   if (nmax[l]>nmaxmax) nmaxmax=nmax[l];

 F=dmatrix(0,lmax,0,nmaxmax);
 F[0][0]=2.;
 F[0][1]=2.;
 F[0][2]=2.;
 F[0][3]=2.;
 F[1][0]=6.;
 F[1][1]=6.;
 F[1][2]=2.;
 F[2][0]=10.;

Confirm your results by comparing with the table at the NIST web site.

To ease the change from the true potential to the pseudopotential, add the following function to the top of the file containing your function getg(),

double Vnuc(double r, double Z)
{
  if (r==0.)
    return 0.;
  else
    return -Z/r;
}

and replace all occurrences along the lines of -Z/r[k] in your code with Vnuc(r[k],Z).
Hint: There is likely only one occurrence, the computation of the total self-consistent potential V[] in getg(). Rerun ``broyden.c'' and verify that the results are unchanged.

More detailed information than simply the total energy is needed to verify the pseudopotential. Modify your code to provide this information as follows:

1. Change the function prototype of getg() to include an additional argument ``flag'' that we will use to get an optional detailed report of the electronic structure. The new prototype will be

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

2. Add to getg() the following code block that produces a detailed report of the electronic structure only if flag$>$0. You should add this block after all of the computations but before the freeing of the dvectors and various other objects.

     /* Get detailed report if flag>0 */
     if (flag>0) {
       char str[6]="spdfg"; /* Strings for angular momentum states */
       FILE *fp;
   
       /* Output various total energy contributions to stdout */
       printf("Etot : %12.6f\n",Etot);
       for (k=0; k<=N; k++)
         integrand[k]=(phi[k]/2+Vnuc(r[k],Z)+Depsxc[k]+Vxc[k])*Rho[k];
       printf("Ekin : %12.6f\n",Etot-simpint(integrand,r,dr,N));
       for (k=0; k<=N; k++)
         integrand[k]=phi[k]*Rho[k]/2;
       printf("Ecoul: %12.6f\n",simpint(integrand,r,dr,N));
       for (k=0; k<=N; k++)
         integrand[k]=Vnuc(r[k],Z)*Rho[k];
       printf("Eenuc: %12.6f\n",simpint(integrand,r,dr,N));
       for (k=0; k<=N; k++)
         integrand[k]=(Depsxc[k]+Vxc[k])*Rho[k];
       printf("Exc  : %12.6f\n",simpint(integrand,r,dr,N));
   
       /* Output state energies to stdout */
       for (l=0; l<=lmax; l++)
         for (n=0; n<=nmax[l]; n++)
   	printf("%d%c[%.0f]   : %12.6f\n",n+l+1,str[l],F[l][n],E[l][n]);
   
       /* Output wave functions to file "psiout" */
       /*   First number is r, the rest are the wave functions in the same */
       /*   order as the states above come out. */
       if ( (fp=fopen("psis","w"))==NULL ) {
         printf("\nError: Cannot open file \"psis\" for output!\n\n");
         exit(1);
       }
       for (k=0; k<N; k++) {
         fprintf(fp,"%le  ",r[k]);
         for (l=0; l<=lmax; l++)
   	for (n=0; n<=nmax[l]; n++)
   	  fprintf(fp," %le",Psi[l][n][k]);
         fprintf(fp,"\n");
       }
       fclose(fp);
     }
   

3. To every current call to getg() in main(), add an additional argument with value ``0'', so that your current code is unaffected by these changes.

4. Add a final call to getg() to get the detailed report. This call will be along the lines of

   getg(g[Itmx+1],Rho[Itmx+1],Z,lmax,nmax,nmaxmax,F,r,dr,N,1);
   

   It should use the most updated version of Rho[], and be placed in main() after all of the computations, outside of all of the loops, and just before the freeing of the various data objects.

Rerun ``broyden.c'' with these modifications, and rename the file containing the output wave functions (``psis'') to ``psis_real'' with a statement along the lines of cp psis psis_real. To view a plot of these functions, run the octave file ``pseudo.m'' from Appendix B. This program will crash looking for ``psis_pseudo'', but will first provide a plot of all of the electronic wave functions of Ge.