diff --git a/src/nwpw/nwpwlib/pseudopotential/atom.c b/src/nwpw/nwpwlib/pseudopotential/atom.c
index 0f7420848c..163064349a 100644
--- a/src/nwpw/nwpwlib/pseudopotential/atom.c
+++ b/src/nwpw/nwpwlib/pseudopotential/atom.c
@@ -1,6 +1,6 @@
 /* atom.c -
    author - Eric Bylaska
-   $Id: atom.c,v 1.11 2007-04-10 18:43:01 d3p708 Exp $
+   $Id: atom.c,v 1.12 2007-04-10 23:30:33 d3p708 Exp $
 */
 
 #include	<stdio.h>
@@ -778,6 +778,21 @@ state_Atom (int nt, int lt)
   return i;
 }
 
+int
+state_RelAtom (int nt, int lt, int st)
+{
+  int i;
+  i = 0;
+  while (((nt != n[i]) || (lt != l[i]) || (st!=s2[i])) && (i <= Ncv))
+    ++i;
+
+  /* Error */
+  if (i > Ncv)
+    printf ("Error: state_RelAtom\n");
+
+  return i;
+}
+
 char *
 name_Atom ()
 {
diff --git a/src/nwpw/nwpwlib/pseudopotential/atom.h b/src/nwpw/nwpwlib/pseudopotential/atom.h
index 2b79917d55..7bed583b35 100644
--- a/src/nwpw/nwpwlib/pseudopotential/atom.h
+++ b/src/nwpw/nwpwlib/pseudopotential/atom.h
@@ -35,6 +35,7 @@ extern int turning_point_Atom (int);
 extern double Zion_Atom ();
 extern double Amass_Atom ();
 extern int state_Atom (int, int);
+extern int state_RelAtom (int, int, int);
 extern char *name_Atom ();
 extern int solver_Type_Atom();
 /* used for setting solver parameters */
diff --git a/src/nwpw/nwpwlib/pseudopotential/pspsolve.c b/src/nwpw/nwpwlib/pseudopotential/pspsolve.c
index d347648406..945499a3b8 100644
--- a/src/nwpw/nwpwlib/pseudopotential/pspsolve.c
+++ b/src/nwpw/nwpwlib/pseudopotential/pspsolve.c
@@ -1,5 +1,5 @@
 /*
- $Id: pspsolve.c,v 1.16 2007-04-10 18:43:01 d3p708 Exp $
+ $Id: pspsolve.c,v 1.17 2007-04-10 23:30:34 d3p708 Exp $
 */
 #include <stdlib.h>
 #include <stdio.h>
@@ -124,7 +124,8 @@ void FATR pspsolve_
 
   if (isRelativistic_Atom ())
     {
-      rpspsolve_ (print_ptr, debug_ptr, lmax_ptr, locp_ptr, rlocal_ptr,
+      printf("Relativstic Atom detected!\n");
+      rpspsolve_(print_ptr, debug_ptr, lmax_ptr, locp_ptr, rlocal_ptr,
 		     sdir_name, n9, dir_name, n0, in_filename, n1,
 		     out_filename, n2);
       return;
@@ -493,7 +494,7 @@ void FATR pspsolve_
       printf ("Outputing semicore density gradient: %s\n", full_filename);
       fp = fopen (full_filename, "w+");
       for (k = 0; k < Ngrid; ++k)
-	fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], drho_semicore_Psp ()[k]);
+	fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], (drho_semicore_Psp ())[k]);
       fclose (fp);
 
       /* output all-electron potential infile.pot */
diff --git a/src/nwpw/nwpwlib/pseudopotential/rpsp.c b/src/nwpw/nwpwlib/pseudopotential/rpsp.c
index 5e924d34b7..d6a2ab90b5 100644
--- a/src/nwpw/nwpwlib/pseudopotential/rpsp.c
+++ b/src/nwpw/nwpwlib/pseudopotential/rpsp.c
@@ -1,6 +1,6 @@
 /* psp.c -
-   author - Patrick Nichols
-   Spin-Orbit Pseudopotential Generation
+   author - Eric Bylaska
+   $Id: rpsp.c,v 1.2 2007-04-10 23:30:34 d3p708 Exp $
 */
 
 #include	<stdio.h>
@@ -13,14 +13,12 @@
 
 #include	"dft.h"
 #include	"atom.h"
-#include	"hamann.h"
+#include	"rhamann.h"
 #include	"troullier.h"
 #include        "vanderbilt.h"
 #include	"generate_rho_semicore.h"
 #include	"rpsp.h"
 #include        "psp.h"
-
-
 #define	False	0
 #define	True	1
 #define	Max(x,y)	((x>y) ? x : y)
@@ -28,37 +26,42 @@
 
 /* atom structure variables */
 
-static int Nvalence;
-static int npsp_states;
-static int *n;
-static int *l;
-static int *spin;
-static int lmax;
-static double *fill;
-static double *rcut;
-static double *peak;
-static double Zion;
-static double Total_E,
-  E_Hartree, P_Hartree, E_exchange, P_exchange, E_correlation, P_correlation;
-static double *eigenvalue;
-static double **r_psi;
-static double **r_psi_prime;
-static double *rho;
-static double *rho_semicore;
-static double r_semicore;
-static double **V_psp;
-static const char *comment="Hamann Pseudopotential";
-static int aSolver;
+static	int	Nvalence;
+static  int     npsp_states;
+static	int	*n;
+static	int	*l;
+static  int     *spin; 
+static	int	lmax;
+static 	double	*fill;
+static 	double	*rcut;
+static  double	*peak;
+static	double	Zion;
+static	double	Total_E,
+		E_Hartree, P_Hartree,
+		E_exchange,
+		P_exchange,
+		E_correlation,
+		P_correlation;
+static	double  *eigenvalue;
+static	double	**r_psi;
+static	double	**r_psi_prime;
+static  double  *rho;
+static  double  *rho_semicore;
+static  double  *drho_semicore;
+static  double  r_semicore;
+static	double	**V_psp;
+static  char    comment[80];
+
 /* extra Vanderbilt parameters */
-static double rlocal, clocal;
-static int ns[10], indx_il[4][10], indx_ijl[4][4][10];
+static double rlocal,clocal;
+static int    ns[10],indx_il[4][10],indx_ijl[4][4][10];
 static double *Vlocal;
 static double **r_hard_psi;
 static double *D0;
 static double *q;
 
 /* solver parameters: this is the Kawai-Weare default */
-static int Solver_Type = Hamann;
+static	int	Solver_Type      = Hamann;
 
 
 /********************************
@@ -67,207 +70,243 @@ static int Solver_Type = Hamann;
  *				*
  ********************************/
 
-void
-init_RelPsp (char *filename)
+void	init_RelPsp(char *filename)
 {
-  int p, p1, p2;
-  int ltmp;
-  double rctmp;
-  char *w;
-  FILE *fp;
+   int	  p,p1,p2;
+   int	  ltmp,semicore_type;
+   double rctmp;
+   char   *w,*tc;
+   FILE	  *fp;
+
+   /* find the comment */
+   strcpy(comment,"Hamann pseudopotential");
+   fp = fopen(filename,"r+");
+   w = get_word(fp);
+   while ((w!=NIL) && (strcmp("<comment>",w)!=0))
+      w = get_word(fp);
+
+   if(w!=NIL)
+   {
+      w = get_word(fp);
+      p  = 0;
+      tc = comment;
+      while ((w!=NIL)&&(strcmp("<end>",w) != 0))
+      {
+        p = (strlen(w));
+        strcpy(tc, w);
+        for(p1=0;p1<p; ++p1) ++tc;
+        strcpy(tc, " ");
+        ++tc;
+
+        w = get_word(fp);
+      }
+   }
+   fclose(fp);
 
 
 
-  Solver_Type = Hamann;
-
-  /* set lmax  */
-  lmax = lmax_Atom () + 1;
-
-  /* set the number psp projectors */
-  npsp_states = 2 * lmax + 1;
-  fp = fopen (filename, "r+");
-  w = get_word (fp);
-  while ((w != NIL) && (strcmp ("<npsp-states>", w) != 0))
-    w = get_word (fp);
-
-  if (w != NIL)
-    {
-      w = get_word (fp);
-      while ((w != NIL) && (strcmp ("<end>", w) != 0))
-	{
-	  sscanf (w, "%d", &ltmp);
-	  w = get_word (fp);
-	  npsp_states = ltmp;
-	}
-    }
-  fclose (fp);
 
 
-  /* allocate memory for n,l,fill,rcut,peak, and eigenvalue */
-  n = (int *) malloc ((npsp_states) * sizeof (int));
-  l = (int *) malloc ((npsp_states) * sizeof (int));
-  spin = (int *) malloc ((npsp_states) * sizeof (int));
-  fill = (double *) malloc ((npsp_states) * sizeof (double));
-  rcut = (double *) malloc ((npsp_states) * sizeof (double));
-  peak = (double *) malloc ((npsp_states) * sizeof (double));
-  eigenvalue = (double *) malloc ((npsp_states) * sizeof (double));
+   /* set lmax  */
+   lmax = lmax_Atom() + 1;
+
+   /* set the number psp projectors */
+   npsp_states = 2*lmax + 2;
+   /* allocate memory for n,l,fill,rcut,peak, and eigenvalue */
+   n    	 = (int *)    malloc((npsp_states)*sizeof(int));
+   l    	 = (int *)    malloc((npsp_states)*sizeof(int));
+   spin          = (int *)    malloc((npsp_states)*sizeof(int));
+   fill 	 = (double *) malloc((npsp_states)*sizeof(double));
+   rcut 	 = (double *) malloc((npsp_states)*sizeof(double));
+   peak 	 = (double *) malloc((npsp_states)*sizeof(double));
+   eigenvalue = (double *) malloc((npsp_states)*sizeof(double));
+
+   
+   /* allocate memory for r_psi, V_psp, and rho */ 
+   r_psi       = (double **) malloc((npsp_states)*sizeof(double*));
+   r_psi_prime = (double **) malloc((npsp_states)*sizeof(double*));
+   V_psp       = (double **) malloc((npsp_states)*sizeof(double*));
+   for (p=0; p<npsp_states; ++p)
+   {
+      r_psi[p]       = alloc_LogGrid();
+      r_psi_prime[p] = alloc_LogGrid();
+      V_psp[p]       = alloc_LogGrid();
+   }
+   rho           = alloc_LogGrid();
+   rho_semicore  = alloc_LogGrid();
+   drho_semicore = alloc_LogGrid();
 
 
-  /* allocate memory for r_psi, V_psp, and rho */
-  r_psi = (double **) malloc ((npsp_states) * sizeof (double *));
-  r_psi_prime = (double **) malloc ((npsp_states) * sizeof (double *));
-  V_psp = (double **) malloc ((npsp_states) * sizeof (double *));
-  for (p = 0; p < npsp_states; ++p)
-    {
-      r_psi[p] = alloc_LogGrid ();
-      r_psi_prime[p] = alloc_LogGrid ();
-      V_psp[p] = alloc_LogGrid ();
-    }
-  rho = alloc_LogGrid ();
-  rho_semicore = alloc_LogGrid ();
+ 
+   /* get the psp info */ 
+   Suggested_Param_RelHamann(&Nvalence,n,l,spin,eigenvalue,fill,rcut);
+   /* set the number psp projectors */
+   fp = fopen(filename,"r+");
+   w = get_word(fp);
+   while ((w!=NIL) && (strcmp("<npsp-states>",w)!=0))
+      w = get_word(fp);
 
-  /* get the psp info */
-  Suggested_Param_RelHamann (&Nvalence, n, l, spin, eigenvalue, fill, rcut);
-  /* set the number psp projectors */
-  fp = fopen (filename, "r+");
-  w = get_word (fp);
-  while ((w != NIL) && (strcmp ("<npsp-states>", w) != 0))
-    w = get_word (fp);
-
-  if (w != NIL)
-    {
-      w = get_word (fp);
-      while ((w != NIL) && (strcmp ("<end>", w) != 0))
-	{
-	  sscanf (w, "%d", &ltmp);
-	  w = get_word (fp);
-	  Nvalence = ltmp;
-	}
-    }
-  fclose (fp);
+   if (w!=NIL)
+   {
+      w = get_word(fp);
+      while ((w!=NIL) && (strcmp("<end>",w) != 0))
+      {
+         sscanf(w,"%d", &ltmp);
+         w = get_word(fp);
+         Nvalence = 2*ltmp+2;
+      }
+   }
+   fclose(fp);
 
 
-  /* get rcut */
-  fp = fopen (filename, "r+");
-  w = get_word (fp);
-  while ((w != NIL) && (strcmp ("<rcut>", w) != 0))
-    w = get_word (fp);
+      /* get rcut */
+      fp = fopen(filename,"r+");
+      w = get_word(fp);
+      while ((w!=NIL) && (strcmp("<rcut>",w)!=0))
+         w = get_word(fp);
 
-  if (w != NIL)
-    {
-      w = get_word (fp);
-      while ((w != NIL) && (strcmp ("<end>", w) != 0))
-	{
-	  sscanf (w, "%d", &ltmp);
-	  w = get_word (fp);
-	  sscanf (w, "%lf", &rctmp);
-	  w = get_word (fp);
-	  rcut[ltmp] = rctmp;
-	}
-    }
-  fclose (fp);
+      if (w!=NIL)
+      { 
+         w = get_word(fp);
+         while ((w!=NIL) && (strcmp("<end>",w) != 0))
+         {
+            sscanf(w,"%d", &ltmp);
+            w = get_word(fp);
+            sscanf(w,"%lf", &rctmp);
+            w = get_word(fp);
+            rcut[2*ltmp] = rctmp;
+            rcut[2*ltmp+1] = rctmp;
+         }
+      }
+      fclose(fp);
 
-  /* get ecut */
-  fp = fopen (filename, "r+");
-  w = get_word (fp);
-  while ((w != NIL) && (strcmp ("<ecut>", w) != 0))
-    w = get_word (fp);
+      /* get ecut */
+      fp = fopen(filename,"r+");
+      w = get_word(fp);
+      while ((w!=NIL) && (strcmp("<ecut>",w)!=0))
+         w = get_word(fp);
 
-  if (w != NIL)
-    {
-      w = get_word (fp);
-      while ((w != NIL) && (strcmp ("<end>", w) != 0))
-	{
-	  sscanf (w, "%d", &ltmp);
-	  w = get_word (fp);
-	  sscanf (w, "%lf", &rctmp);
-	  w = get_word (fp);
-	  eigenvalue[ltmp] = rctmp;
-	}
-    }
-  fclose (fp);
+      if (w!=NIL)
+      {
+         w = get_word(fp);
+         while ((w!=NIL) && (strcmp("<end>",w) != 0))
+         {
+            sscanf(w,"%d", &ltmp);
+            w = get_word(fp);
+            sscanf(w,"%lf", &rctmp);
+            w = get_word(fp);
+            eigenvalue[2*ltmp] = rctmp;
+            eigenvalue[2*ltmp+1] = rctmp;
+         }
+      }
+      fclose(fp);
+      
+   /* get r_semicore - if zero then no core corrections added */
+   r_semicore = 0.0;
+   fp = fopen(filename,"r+");
+   w = get_word(fp);
+   while ((w!=NIL) && (strcmp("<semicore>",w)!=0))
+      w = get_word(fp);
+
+   if (w!=NIL)
+   {
+      w = get_word(fp);
+      while ((w!=NIL) && (strcmp("<end>",w) != 0))
+      {
+         sscanf(w,"%lf", &rctmp);
+         w = get_word(fp);
+         r_semicore = rctmp;
+      }
+   }
+   fclose(fp);
+
+   /* find the semicore_type */
+   semicore_type = 2;
+   fp = fopen(filename,"r+");
+   w = get_word(fp);
+   while ((w!=NIL) && (strcmp("<semicore_type>",w)!=0))
+      w = get_word(fp);
+   if (w!=NIL)
+   {
+      w = get_word(fp);
+      if (strcmp("quadratic",w)==0)   semicore_type = 0;
+      if (strcmp("louie",w)==0)       semicore_type = 1;
+      if (strcmp("fuchs",w)==0)       semicore_type = 2;
+   }
+   fclose(fp);
 
 
 
-  /* get r_semicore - if zero then no core corrections added */
-  r_semicore = 0.0;
-  fp = fopen (filename, "r+");
-  w = get_word (fp);
-  while ((w != NIL) && (strcmp ("<semicore>", w) != 0))
-    w = get_word (fp);
 
-  if (w != NIL)
-    {
-      w = get_word (fp);
-      while ((w != NIL) && (strcmp ("<end>", w) != 0))
-	{
-	  sscanf (w, "%lf", &rctmp);
-	  w = get_word (fp);
-	  r_semicore = rctmp;
-	}
-    }
-  fclose (fp);
+   /* generate non-zero rho_semicore */
+   if (r_semicore > 0.0)
+   {
+      /*
+      printf("\n\n");
+      printf("Generating non-zero semicore density\n");
+      */
+      generate_rho_semicore(semicore_type,rho_core_Atom(),r_semicore,rho_semicore);
+      Derivative_LogGrid(rho_semicore,drho_semicore);
+   }
 
-  /* generate non-zero rho_semicore */
-  if (r_semicore > 0.0)
-    {
-      printf ("\n\n");
-      printf ("Generating non-zero semicore density\n");
-      generate_rho_semicore (rho_core_Atom (), r_semicore, rho_semicore);
-    }
+   /* define the ion charge */
+   /*
+     Zion=0.0;
+     for (p=Ncore_Atom(); p<(Ncore_Atom()+Nvalence_Atom()); ++p)
+        Zion += fill_Atom(p);
+   */
 
-  /* define the ion charge */
-  Zion = 0.0;
-  for (p = Ncore_Atom (); p < (Ncore_Atom () + Nvalence_Atom ()); ++p)
-    Zion += fill_Atom (p);
+   Zion = Zion_Atom();
+   for (p=0; p<(Ncore_Atom()); ++p)
+      Zion -= fill_Atom(p);
 
 
-}				/* init_RelPsp */
+
+   
+} /* init_RelPsp */
 
 
 
 /********************************
  *				*
- *	   solve_RelPsp 	*
+ *	   solve_RelPsp 		*
  *				*
  ********************************/
 
-void
-solve_RelPsp ()
+void	solve_RelPsp()
 {
 
-  int p, k, Ngrid;
-  double *rgrid;
+   int 		p,k,Ngrid;
+   double	*rgrid;
 
-  /* get loggrid variables */
-  Ngrid = N_LogGrid ();
-  rgrid = r_LogGrid ();
+   /* get loggrid variables */
+   Ngrid = N_LogGrid();
+   rgrid = r_LogGrid();
 
-  solve_RelHamann (Nvalence, n, l, spin, eigenvalue, fill, rcut,
-		r_psi, r_psi_prime, rho, rho_semicore, V_psp,
+      solve_RelHamann(Nvalence,n,l,spin,eigenvalue,fill,rcut,
+		r_psi,r_psi_prime,rho,rho_semicore,V_psp,
 		&Total_E,
-		&E_Hartree, &P_Hartree,
-		&E_exchange, &P_exchange, &E_correlation, &P_correlation);
-
+		&E_Hartree,&P_Hartree,
+		&E_exchange,&P_exchange,
+		&E_correlation,&P_correlation);
    /******************************************************/
-  /* find the outermost peak of the pseudowavefunctions */
+   /* find the outermost peak of the pseudowavefunctions */
    /******************************************************/
-  for (p = 0; p < Nvalence; ++p)
-    {
+   for (p=0; p<Nvalence; ++p)
+   {
       if (fill[p] != 0.0)
-	{
-	  k = Ngrid - 2;
-	  while ((r_psi_prime[p][k] * r_psi_prime[p][k + 1] >= 0.0)
-		 && (k >= 0))
-	    --k;
-	  peak[p] = rgrid[k];
-	}
+      {
+         k = Ngrid-2;
+         while ((r_psi_prime[p][k]*r_psi_prime[p][k+1] >= 0.0) && (k>=0))
+            --k;
+         peak[p] = rgrid[k];
+      }
       else
-	peak[p] = rgrid[Ngrid - 1];
-    }
+         peak[p] = rgrid[Ngrid-1];
+   }
+      
 
-
-}				/* solve_RelPsp */
+} /* solve_RelPsp */
 
 
 
@@ -277,81 +316,79 @@ solve_RelPsp ()
  *				*
  ********************************/
 
-char *
-solver_Name_RelPsp ()
+char	*solver_Name_RelPsp()
 {
-  char *s;
-  if (Solver_Type == Hamann)
-    s = "Hamann";
-  else if (Solver_Type == Troullier)
-    s = "Troullier-Martins";
-  else if (Solver_Type == Vanderbilt)
-    s = "Vanderbilt";
-  else
-    s = "unknown?";
-
-  return s;
+   char *s;
+      s = "Hamann";
+   return s;
 }
 
 
-void
-print_RelPsp (fp)
-     FILE *fp;
+void	print_RelPsp(fp)
+
+FILE 	*fp;
 {
-  int i;
+   int i;
+   
+   fprintf(fp,"\n\n");
+   fprintf(fp,"PSP solver information\n\n");
+   fprintf(fp,"Atom name: %s\n",name_Atom());
+   fprintf(fp,"Zcharge  : %le\n",Zion);
+   fprintf(fp,"Nvalence : %d\n",Nvalence);
+   fprintf(fp,"         : restricted calculation\n");
 
-  fprintf (fp, "\n\n");
-  fprintf (fp, "PSP solver information\n\n");
-  fprintf (fp, "Atom name: %s\n", name_Atom ());
-  fprintf (fp, "Zcharge  : %le\n", Zion);
-  fprintf (fp, "Nvalence : %d\n", Nvalence);
-  fprintf (fp, "         : restricted calculation\n");
+   fprintf(fp,"\n------ Solver information ------\n");
+   fprintf(fp,"solver type      : %s\n",solver_Name_RelPsp());
+   fprintf(fp,"hartree type     : %s\n",hartree_Name_DFT());
+   fprintf(fp,"exchange type    : %s\n",exchange_Name_DFT());
+   if (strcmp(exchange_Name_DFT(),"Dirac") == 0)
+     fprintf(fp,"           alpha : %lf\n",Dirac_alpha());
+   fprintf(fp,"correlation type : %s\n",correlation_Name_DFT());
 
-  fprintf (fp, "\n------ Solver information ------\n");
-  fprintf (fp, "solver type      : %s\n", solver_Name_RelPsp ());
-  fprintf (fp, "hartree type     : %s\n", hartree_Name_DFT ());
-  fprintf (fp, "exchange type    : %s\n", exchange_Name_DFT ());
-  if (strcmp (exchange_Name_DFT (), "Dirac") == 0)
-    fprintf (fp, "           alpha : %lf\n", Dirac_alpha ());
-  fprintf (fp, "correlation type : %s\n", correlation_Name_DFT ());
+   fprintf(fp,"----------------------------------------------------------------------------\n");
+   fprintf(fp,"n\tl\tpopulation\tEcut\t\tRcut\t\tOuter Peak\n");
+   
+   for (i=0; i<Nvalence; ++i)
+   {
+        fprintf(fp,"%d\t%s\t%.1lf\t%.2lf\t\t%le\t%le\t%le\n",l[i]+1,spd_Name(l[i]),
+					0.5*spin[i], 
+				        fill[i],eigenvalue[i],rcut[i],peak[i]);
+   }
+   fprintf(fp,"----------------------------------------------------------------------------\n");
+   if (r_semicore > 0.0)
+   {
+     fprintf(fp,"SemiCore Corrections Added\n");
+     fprintf(fp,"    rcore                      : %lf\n",r_semicore);
+     fprintf(fp,"    Semicore Charge            : %lf\n",
+                 Integrate_LogGrid(rho_semicore));
+     fprintf(fp,"    Semicore Charge gradient   : %lf\n",
+                 Integrate_LogGrid(drho_semicore));
+   }
 
-  fprintf (fp,
-	   "----------------------------------------------------------------------------\n");
-  fprintf (fp, "n\tl\tpopulation\tEcut\t\tRcut\t\tOuter Peak\n");
+   if (Solver_Type==Vanderbilt) 
+      fprintf(fp,"Using AE valence density for descreening\n");
 
-  for (i = 0; i < Nvalence; ++i)
-    {
-      fprintf (fp, "%d\t%s\t%3.1lf\t%s\t\t%le\t%le\t%le\n", l[i] + 1,
-	       spd_Name (l[i]), 0.5 * spin[i], "  -", eigenvalue[i], rcut[i],
-	       peak[i]);
-    }
-  fprintf (fp,
-	   "----------------------------------------------------------------------------\n");
-  if (r_semicore > 0.0)
-    {
-      fprintf (fp, "SemiCore Corrections Added\n");
-      fprintf (fp, "    rcore             : %lf\n", r_semicore);
-      fprintf (fp, "    Semicore Charge   : %lf\n",
-	       Integrate_LogGrid (rho_semicore));
-    }
-  fprintf (fp, "Pseudopotential Charge: %lf\n", Integrate_LogGrid (rho));
-  fprintf (fp, "\nTotal E       = %le\n", Total_E);
-  fprintf (fp, "\n");
+   fprintf(fp,  "Pseudopotential ion charge       = %lf\n", Zion);
+   fprintf(fp,  "Pseudopotential electronic charge= %lf\n", -Integrate_LogGrid(rho));
+   fprintf(fp,  "Pseudopotential atom charge      = %lf\n", Zion-Integrate_LogGrid(rho));
 
+   fprintf(fp,"\nTotal E       = %le\n",Total_E);
+   fprintf(fp,"\n");
+   
 
-  fprintf (fp, "E_Hartree     = %le\n", E_Hartree);
-  fprintf (fp, "<Vh>          = %le\n", P_Hartree);
-  fprintf (fp, "\n");
+   fprintf(fp,"E_Hartree     = %le\n",E_Hartree);
+   fprintf(fp,"<Vh>          = %le\n",P_Hartree);
+   fprintf(fp,"\n");
 
-  fprintf (fp, "E_exchange    = %le\n", E_exchange);
-  fprintf (fp, "<Vx>          = %le\n", P_exchange);
-  fprintf (fp, "\n");
+   fprintf(fp,"E_exchange    = %le\n",E_exchange);
+   fprintf(fp,"<Vx>          = %le\n",P_exchange);
+   fprintf(fp,"\n");
 
-  fprintf (fp, "E_correlation = %le\n", E_correlation);
-  fprintf (fp, "<Vc>          = %le\n\n", P_correlation);
-
-
-}				/* print_Atom */
+   fprintf(fp,"E_correlation = %le\n",E_correlation);
+   fprintf(fp,"<Vc>          = %le\n\n",P_correlation);
+   
+   
+} /* print_Atom */
 
 /********************************
  *				*
@@ -359,23 +396,30 @@ print_RelPsp (fp)
  *				*
  ********************************/
 
-void
-set_Solver_RelPsp (solver)
-     int solver;
+void set_Solver_RelPsp(solver)
+int solver;
 {
-  Solver_Type = solver;
+   fprintf(stdout,"RelRelPsp::Cannot set the solver yet!\n");
 }
 
 
-int
-Vanderbilt_RelPsp ()
+int Vanderbilt_RelPsp()
 {
-  int value;
+   int value;
 
-  value = 0;
-  if (Solver_Type == Vanderbilt)
-    value = 1;
-  return value;
+   value = 0;
+   if (Solver_Type==Vanderbilt) value = 1;
+   return value;
+}
+
+int NormConserving_RelPsp()
+{
+   int value;
+
+   value = 0;
+   if (Solver_Type==Hamann)    value = 1;
+   if (Solver_Type==Troullier) value = 1;
+   return value;
 }
 
 /********************************
@@ -384,169 +428,149 @@ Vanderbilt_RelPsp ()
  *				*
  ********************************/
 
-double
-E_RelPsp ()
+double	E_RelPsp()
 {
-  return Total_E;
-}				/*E_Atom */
+   return Total_E;
+} /*E_Atom*/
 
 
-double
-eigenvalue_RelPsp (int i)
+double	eigenvalue_RelPsp(int i)
 {
-  return eigenvalue[i];
+   return eigenvalue[i];
 }
 
-double *
-rho_RelPsp ()
+double	*rho_RelPsp()
 {
-  return rho;
+   return rho;
 }
 
-double *
-rho_semicore_RelPsp ()
+double	*rho_semicore_RelPsp()
 {
-  return rho_semicore;
+   return rho_semicore;
 }
 
-double
-r_semicore_RelPsp ()
+double * drho_semicore_RelPsp()
 {
-  return r_semicore;
+   return drho_semicore;
+}
+
+double	r_semicore_RelPsp()
+{
+   return r_semicore;
 }
 
 
-double *
-Beta_RelPsp (int i, int l)
+double *Beta_RelPsp(int i,int l)
 {
-  return V_psp[indx_il[i][l]];
+   return V_psp[indx_il[i][l]];
 }
 
-double *
-r_psi_il_RelPsp (int i, int l)
+double *r_psi_il_RelPsp(int i, int l)
 {
-  return r_psi[indx_il[i][l]];
+   return r_psi[indx_il[i][l]];
 }
 
-double *
-r_hard_psi_il_RelPsp (int i, int l)
+double *r_hard_psi_il_RelPsp(int i, int l)
 {
-  return r_hard_psi[indx_il[i][l]];
+   return r_hard_psi[indx_il[i][l]];
 }
 
-int
-ns_RelPsp (int l)
+int ns_RelPsp(int l)
 {
-  return ns[l];
+   return ns[l];
 }
 
-double
-D0_RelPsp (int i, int j, int l)
+double D0_RelPsp(int i, int j, int l)
 {
-  return D0[indx_ijl[i][j][l]];
+   return D0[indx_ijl[i][j][l]];
 }
 
-double
-q_RelPsp (int i, int j, int l)
+double q_RelPsp(int i, int j, int l)
 {
-  return q[indx_ijl[i][j][l]];
+   return q[indx_ijl[i][j][l]];
 }
 
-double *
-Vlocal_RelPsp ()
+double *Vlocal_RelPsp()
 {
-  return Vlocal;
+   return Vlocal;
 }
 
 
 
-double *
-V_RelPsp (int i)
+double	*V_RelPsp(int i)
 {
-  return V_psp[i];
+   return V_psp[i];
 }
 
 
-double *
-r_psi_RelPsp (int i)
+double	*r_psi_RelPsp(int i)
 {
-  return r_psi[i];
+   return r_psi[i];
 }
 
-int
-n_RelPsp (int i)
+int	n_RelPsp(int i)
 {
-  return n[i];
+   return n[i];
 }
 
-int
-l_RelPsp (int i)
+int	l_RelPsp(int i)
 {
-  return l[i];
+   return l[i];
 }
 
-int
-lmax_RelPsp ()
+int	lmax_RelPsp()
 {
-  return lmax;
+   return lmax;
 }
 
-double
-fill_RelPsp (int i)
+double	fill_RelPsp(int i)
 {
-  return fill[i];
+   return fill[i];
 }
 
 
-int
-Nvalence_RelPsp ()
+int	Nvalence_RelPsp()
 {
-  return Nvalence;
+   return Nvalence;
 }
 
-double
-peak_RelPsp (int i)
+double	peak_RelPsp(int i)
 {
-  return peak[i];
+   return peak[i];
 }
 
-double
-rcut_RelPsp (int i)
+double	rcut_RelPsp(int i)
 {
-  return rcut[i];
+   return rcut[i];
 }
 
-double
-rcut_il_RelPsp (int i, int l)
+double	rcut_il_RelPsp(int i, int l)
 {
-  return rcut[indx_il[i][l]];
+   return rcut[ indx_il[i][l] ];
 }
 
-double
-Zion_RelPsp ()
+double	Zion_RelPsp()
 {
-  return Zion;
+   return Zion;
 }
 
-int
-state_RelPsp (int nt, int lt)
+int	state_RelPsp(int nt, int lt)
 {
-  int i;
-  i = 0;
-  while (((nt != n[i]) || (lt != l[i])) && (i < Nvalence))
-    ++i;
+   int i;
+   i = 0;
+   while ( ((nt != n[i]) || (lt != l[i])) && (i<Nvalence) )
+     ++i;
 
-  /* Error */
-  if (i >= Nvalence)
-    printf ("Error: state_RelPsp\n");
+   /* Error */
+   if (i>=Nvalence)
+      printf("Error: state_RelPsp\n");
 
-  return i;
+   return i;
 }
 
 
-
-const char *
-comment_RelPsp ()
+char    *comment_RelPsp()
 {
-  return comment;
+   return comment;
 }
+
diff --git a/src/nwpw/nwpwlib/pseudopotential/rpsp.h b/src/nwpw/nwpwlib/pseudopotential/rpsp.h
index ab0a0f607e..191d0105cc 100644
--- a/src/nwpw/nwpwlib/pseudopotential/rpsp.h
+++ b/src/nwpw/nwpwlib/pseudopotential/rpsp.h
@@ -1,39 +1,39 @@
 #ifndef _RPSP_H_
 #define _RPSP_H_
 #include <stdio.h>
-/* 
- * rpsp.c 
- *   author: Patrick Nichols
- */
-extern void init_RelPsp (char *filename);
-extern void solve_RelPsp (void);
-extern char *solver_Name_RelPsp (void);
-extern void print_RelPsp (FILE * fp);
-extern void set_Solver_RelPsp (int solver);
-extern intVanderbilt_RelPsp (void);
-extern double E_RelPsp (void);
-extern double eigenvalue_RelPsp (int i);
-extern double *rho_RelPsp (void);
-extern double *rho_semicore_RelPsp (void);
-extern double r_semicore_RelPsp (void);
-extern double *Beta_RelPsp (int i, int l);
-extern double *r_psi_il_RelPsp (int i, int l);
-extern double *r_hard_psi_il_RelPsp (int i, int l);
-extern int ns_RelPsp (int l);
-extern double D0_RelPsp (int i, int j, int l);
-extern double q_RelPsp (int i, int j, int l);
-extern double *Vlocal_RelPsp (void);
-extern double *V_RelPsp (int i);
-extern double *r_psi_RelPsp (int i);
-extern int n_RelPsp (int i);
-extern int l_RelPsp (int i);
-extern int lmax_RelPsp (void);
-extern double fill_RelPsp (int i);
-extern int Nvalence_RelPsp (void);
-extern double peak_RelPsp (int i);
-extern double rcut_RelPsp (int i);
-extern double rcut_il_RelPsp (int i, int l);
-extern double Zion_RelPsp (void);
-extern int state_RelPsp (int nt, int lt);
-extern const char *comment_RelPsp (void);
+
+/* rpsp.c */
+extern void init_RelPsp(char *filename);
+extern void solve_RelPsp(void);
+extern char *solver_Name_RelPsp(void);
+extern void print_RelPsp(FILE *fp);
+extern void set_Solver_RelPsp(int solver);
+extern int Vanderbilt_RelPsp(void);
+extern int NormConserving_RelPsp(void);
+extern double E_RelPsp(void);
+extern double eigenvalue_RelPsp(int i);
+extern double *rho_RelPsp(void);
+extern double *rho_semicore_RelPsp(void);
+extern double *drho_semicore_RelPsp(void);
+extern double r_semicore_RelPsp(void);
+extern double *Beta_RelPsp(int i, int l);
+extern double *r_psi_il_RelPsp(int i, int l);
+extern double *r_hard_psi_il_RelPsp(int i, int l);
+extern int ns_RelPsp(int l);
+extern double D0_RelPsp(int i, int j, int l);
+extern double q_RelPsp(int i, int j, int l);
+extern double *Vlocal_RelPsp(void);
+extern double *V_RelPsp(int i);
+extern double *r_psi_RelPsp(int i);
+extern int n_RelPsp(int i);
+extern int l_RelPsp(int i);
+extern int lmax_RelPsp(void);
+extern double fill_RelPsp(int i);
+extern int Nvalence_RelPsp(void);
+extern double peak_RelPsp(int i);
+extern double rcut_RelPsp(int i);
+extern double rcut_il_RelPsp(int i, int l);
+extern double Zion_RelPsp(void);
+extern int state_RelPsp(int nt, int lt);
+extern char *comment_RelPsp(void);
 #endif
diff --git a/src/nwpw/nwpwlib/pseudopotential/rpspsolve.c b/src/nwpw/nwpwlib/pseudopotential/rpspsolve.c
index 7d0819df1b..f3282d1028 100644
--- a/src/nwpw/nwpwlib/pseudopotential/rpspsolve.c
+++ b/src/nwpw/nwpwlib/pseudopotential/rpspsolve.c
@@ -106,15 +106,15 @@ void FATR rpspsolve_
 
 /* we have already solved for the Relativsitic AE wavefunctions using atom */
   set_debug_print (debug);
-  init_Psp (infile);
-  solve_Psp ();
+  init_RelPsp (infile);
+  solve_RelPsp ();
   if (debug)
-    print_Psp (stdout);
+    print_RelPsp (stdout);
 
   init_Linear (infile);
 
   /* allocate linear meshes */
-  Nvalence = Nvalence_Psp ();
+  Nvalence = Nvalence_RelPsp ();
   Nlinear = nrl_Linear ();
   psil = (double **) malloc (Nvalence * sizeof (double *));
   pspl = (double **) malloc (Nvalence * sizeof (double *));
@@ -129,16 +129,16 @@ void FATR rpspsolve_
   /* Norm-conserving output */
   for (p = 0; p < Nvalence; ++p)
     {
-      Log_to_Linear (r_psi_Psp (p), rl, psil[p]);
-      Log_to_Linear_zero (V_Psp (p), rl, pspl[p]);
+      Log_to_Linear (r_psi_RelPsp (p), rl, psil[p]);
+      Log_to_Linear_zero (V_RelPsp (p), rl, pspl[p]);
 
       /* normalize scattering state */
-      if (fill_Psp (p) == 0.0)
+      if (fill_RelPsp (p) == 0.0)
 	{
 	  normalize_Linear (psil[p]);
 	}
     }
-  Log_to_Linear (rho_Psp (), rl, rhol);
+  Log_to_Linear (rho_RelPsp (), rl, rhol);
 
   if (debug)
     {
@@ -211,23 +211,23 @@ void FATR rpspsolve_
     }
   fp = fopen (outfile, "w+");
   fprintf (fp, "%s\n", name_Atom ());
-  fprintf (fp, "%lf %lf %d   %d %d %lf\n", Zion_Psp (), Amass_Atom (),
-	   lmax_Psp (), lmax_out, locp_out, rlocal_out);
-  for (p = 0; p <= lmax_Psp (); ++p)
-    fprintf (fp, "%lf ", rcut_Psp (p));
+  fprintf (fp, "%lf %lf %d   %d %d %lf\n", Zion_RelPsp (), Amass_Atom (),
+	   lmax_RelPsp (), lmax_out, locp_out, rlocal_out);
+  for (p = 0; p <= lmax_RelPsp (); ++p)
+    fprintf (fp, "%lf ", rcut_RelPsp (p));
   fprintf (fp, "\n");
   fprintf (fp, "%d %lf\n", nrl_Linear (), drl_Linear ());
-  fprintf (fp, "%s\n", comment_Psp ());
+  fprintf (fp, "%s\n", comment_RelPsp ());
 
   if (print)
     {
       printf ("  + Appending pseudopotentials:    %s thru %s\n",
-	      spd_Name (0), spd_Name (lmax_Psp ()));
+	      spd_Name (0), spd_Name (lmax_RelPsp ()));
     }
   for (k = 0; k < Nlinear; ++k)
     {
       fprintf (fp, "%12.8lf", rl[k]);
-      for (p = 0; p <= lmax_Psp (); ++p)
+      for (p = 0; p <= lmax_RelPsp (); ++p)
 	{
 	  vx = p * pspl[2 * p][k] + (p + 1.) * pspl[2 * p + 1][k];
 	  vx *= 1. / (2. * p + 1.);
@@ -238,25 +238,25 @@ void FATR rpspsolve_
   if (print)
     {
       printf ("  + Appending pseudowavefunctions: %s thru %s\n",
-	      spd_Name (0), spd_Name (lmax_Psp ()));
+	      spd_Name (0), spd_Name (lmax_RelPsp ()));
     }
   for (k = 0; k < Nlinear; ++k)
     {
       fprintf (fp, "%12.8lf", rl[k]);
-      for (p = 0; p <= lmax_Psp (); ++p)
+      for (p = 0; p <= lmax_RelPsp (); ++p)
 	fprintf (fp, " %12.8lf", psil[p][k]);
       fprintf (fp, "\n");
     }
 
   /* append semicore corrections */
-  if (r_semicore_Psp () != 0.0)
+  if (r_semicore_RelPsp () != 0.0)
     {
       if (print)
 	{
 	  printf ("  + Appending semicore density\n");
 	}
-      Log_to_Linear (rho_semicore_Psp (), rl, rhol);
-      fprintf (fp, "%lf\n", r_semicore_Psp ());
+      Log_to_Linear (rho_semicore_RelPsp (), rl, rhol);
+      fprintf (fp, "%lf\n", r_semicore_RelPsp ());
       for (k = 0; k < Nlinear; ++k)
 	fprintf (fp, "%12.8lf %12.8lf\n", rl[k],
 		 fabs (rhol[k] * over_fourpi));
@@ -265,7 +265,7 @@ void FATR rpspsolve_
 	{
 	  printf ("  + Appending semicore density gradient\n");
 	}
-      Log_to_Linear (drho_semicore_Psp (), rl, rhol);
+      Log_to_Linear (drho_semicore_RelPsp (), rl, rhol);
       for (k = 0; k < Nlinear; ++k)
 	fprintf (fp, "%12.8lf %12.8lf\n", rl[k], (rhol[k] * over_fourpi));
     }
@@ -283,23 +283,23 @@ void FATR rpspsolve_
   strncat (soutfile, ".so", 3);
   fp = fopen (soutfile, "w+");
   fprintf (fp, "%s\n", name_Atom ());
-  fprintf (fp, "%lf %lf %d   %d %d %lf\n", Zion_Psp (), Amass_Atom (),
-	   lmax_Psp (), lmax_out, locp_out, rlocal_out);
-  for (p = 0; p <= lmax_Psp (); ++p)
-    fprintf (fp, "%lf ", rcut_Psp (p));
+  fprintf (fp, "%lf %lf %d   %d %d %lf\n", Zion_RelPsp (), Amass_Atom (),
+	   lmax_RelPsp (), lmax_out, locp_out, rlocal_out);
+  for (p = 0; p <= lmax_RelPsp (); ++p)
+    fprintf (fp, "%lf ", rcut_RelPsp (p));
   fprintf (fp, "\n");
   fprintf (fp, "%d %lf\n", nrl_Linear (), drl_Linear ());
-  fprintf (fp, "%s\n", comment_Psp ());
+  fprintf (fp, "%s\n", comment_RelPsp ());
 
   if (print)
     {
       printf ("  + Appending Spin-Orbit pseudopotentials:    %s thru %s\n",
-	      spd_Name (0), spd_Name (lmax_Psp ()));
+	      spd_Name (0), spd_Name (lmax_RelPsp ()));
     }
   for (k = 0; k < Nlinear; ++k)
     {
       fprintf (fp, "%12.8lf", rl[k]);
-      for (p = 1; p <= lmax_Psp (); ++p)
+      for (p = 1; p <= lmax_RelPsp (); ++p)
 	{
 	  vx = pspl[2 * p][k] - pspl[2 * p + 1][k];
 	  vx *= 2. / (2. * p + 1.);
@@ -413,7 +413,7 @@ void FATR rpspsolve_
       printf ("Outputing semicore density: %s\n", full_filename);
       fp = fopen (full_filename, "w+");
       for (k = 0; k < Ngrid; ++k)
-	fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_semicore_Psp ()[k]);
+	fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], rho_semicore_RelPsp ()[k]);
       fclose (fp);
 
       /* output semicore density gradient infile.sddns */
@@ -431,7 +431,10 @@ void FATR rpspsolve_
       printf ("Outputing semicore density gradient: %s\n", full_filename);
       fp = fopen (full_filename, "w+");
       for (k = 0; k < Ngrid; ++k)
-	fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], drho_semicore_Psp ()[k]);
+      {
+        vx=drho_semicore_RelPsp()[k];
+	fprintf (fp, "%12.8lf %12.8lf\n", rgrid[k], vx);
+      }
       fclose (fp);
 
       /* output all-electron potential infile.pot */
@@ -454,7 +457,7 @@ void FATR rpspsolve_
       c = 0.0;
       for (k = 0; k < Ngrid; ++k)
 	{
-	  x = rgrid[k] / rcut_Psp (0);
+	  x = rgrid[k] / rcut_RelPsp (0);
 	  x = pow (x, 3.5);
 	  x = exp (-x);
 	  y = 1.0 - x;