/*
    DFT++ is a density functional package developed by the research group
    of Professor Tomas Arias

    Copyright 1996-2003 Sohrab Ismail-Beigi

    This file is part of DFT++.

    DFT++ is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    DFT++ is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with DFT++; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

    Please see the file CREDITS for a list of authors.

    For academic users, we request that publications using results obtained with
    this software reference

    "New algebraic formulation of density functional calculation," by Sohrab Ismail-Beigi
    and T.A. Arias, Computer Physics Communications 128:1-2, 1-45 (June 2000).

    and, if using the wavelet basis, further reference

    "Multiresolution analysis of electronic structure: semicardinal and wavelet bases,"
    T.A. Arias, Reviews of Modern Physics 71:1, 267-311 (January 1999).

    and 

    "Robust ab initio calculation of condensed matter: transparent convergence through
    semicardinal multiresolution analysis,'' I.P. Daykov, T.A. Arias, and
    Torkel D. Engeness, Physical Review Letters, 90:21, 216402 (May 2003).

    For your convenience, preprints of the above articles may be obtained from
    http://arXiv.org/abs/cond-mat/9909130, 9805262, and 0204411, respectively.
*/

#include "header.h"

PW_IonicPotential::PW_IonicPotential()
{
  potential_type = DFT_MIT_PSP;

  dEperNatoms_dNGperVol = 0.0;

  strcpy(potfilename,"");
  strcpy(pulayfilename,"");

  
  ngrid_loc = 0;
  dq_loc = 0.0;
  V_loc = 0;

  nlm = 0;
  l = m = ngamma = 0;
  M = 0;
  ngrid_nl = 0;
  dq_nl = 0;
  flq = 0;

  p_nl = 0;
  p_l = p_np = 0;
  p_ngrid_nl = 0;
  p_dq_nl = 0;
  p_flq = 0;
  Z = 0;

  // for FHI-type PSP
  fhi_max_l = fhi_loc_l = 'd';  // default max-l is d, local component is d
  fhi_ngrid_loc = 2004;  // default Vloc total number of grid points.
  fhi_dq_loc    = 0.02;  // default Vloc grid seperation.
  fhi_ngrid_nl  = 501;   // default flq total number of grid points.
  fhi_dq_nl     = 0.02;  // default flq grid seperation.
  
}


PW_IonicPotential::~PW_IonicPotential()
{
  myfree(V_loc);
  for (int lm=0; lm < nlm; lm++)
    {
      myfree(ngrid_nl[lm]);
      myfree(dq_nl[lm]);
      for (int g=0; g < ngamma[lm]; g++)
	myfree(flq[lm][g]);
      myfree(flq[lm]);
      M[lm].freemem();
    }
  if (nlm>0)
    {
      myfree(l);
      myfree(m);
      myfree(ngamma);
      myfree(M);
      myfree(ngrid_nl);
      myfree(dq_nl);
      myfree(flq);
    }
  for (int ll=0; ll < p_nl; ll++)
    {
      myfree(p_ngrid_nl[ll]);
      myfree(p_dq_nl[ll]);
      for (int p=0; p < p_np[ll]; p++)
	myfree(p_flq[ll][p]);
	  myfree(p_flq[ll]);
    }
  if (p_nl>0)
    {
      myfree(p_l);
      myfree(p_np);
      myfree(p_ngrid_nl);
      myfree(p_dq_nl);
      myfree(p_flq);
    }
}


//
// Read the pseudopotential for this species from file
//
/*
 * Read the file <potfilename> to get all the specifications for the
 * ionic pseudopotential
 *
 ***************************************************************************
 * (3) The format of <potfilename> is:
 *        [#comment]
 *        ngrid_loc dq_loc <locvqfile>
 *        nlm
 *        # for i = 1 to nlm {
 *            l_i m_i ngamma_i     [# comment]
 *            M_i
 *            # for j = 1 to ngamma_i {
 *                ngrid_ij dq_ij <flqfile_ij> [# comment]
 *            # }
 *        # }
 *  if we have projectors, then after this,
 *        p_nl
 *        # for i = 1 to p_nl {
 *            p_l_i p_np_i
 *            # for j = 1 to p_np_i {
 *                p_ngrid_ij p_dq_ij <p_flqfile_ij>
 *            # }
 *        # }
 *  
 * -> ngrid_loc is an integer specifying the number of 
 *    grid points specifying the local pseudopot.
 * -> dq_loc is a real number specifying the spacing in q-space of the grid.
 * -> <locvqfile> contains the actual values of the local pseudopotential.
 * -> nlm is an integer specifying the number of (l,m) states for the 
 *    non-local pseudopotentials.
 * -> l and m are integers and are the values of l and m.
 * -> ngamma is an integer specifying the number of multiple-projector
 *    states for the current (l,m).
 * -> M is a ngamma x ngamma matrix of real numbers specifying the 
 *    details of the multiple-projector states at the current (l,m).
 *    M is read a as a string of real numbers separated by spaces, with
 *    the column-index of M changing fastest.  The last number specifying
 *    M MUST BE FOLLOWED BY A CARRIAGE RETURN (i.e must end the line).
 * -> ngrid (int) is an integer specifying the number of grid-points defining
 *    the non-local pseudopotential for the current (l,m,j) state.
 * -> dq (real) is the regular q-space spacing for the potential in the
 *    (l,m,j) states.
 * -> <flqfile> contains the actual values of the pseudopotential for (l,m,j)
 *
 * Note:
 * -> [# comment] means any text from here to the end-of-line is ignored
 * -> Lines starting with '#' are commenst and are ignored.
 *
 * -> IMPORTANT NOTE REGARDING M:
 *    THERE CAN BE NO COMMENTS AT ALL SEPARATING THE VALUES THAT SPECIFY
 *    THE M-MATRIX.  FURTHREMORE, THE VERY LAST VALUE OF M SPECIFIED MUST
 *    BE FOLLOWED BY A CARRIAGE RETURN (I.E. MUST END THE LINE IT IS ON).
 *
 ***************************************************************************
 * (4) The format of the <locvqfile> and <flqfile> files:
 *
 *            0     v(0)
 *            dq    v(dq)
 *            2*dq  v(2*dq)
 *            3*dq  v(3*dq)
 *               .....
 *
 * where dq is the spacing specified for the locvq/flq files, and v(q) is
 * value of the pseudopotential at q.
 * There must be ngrid pairs of numbers (q,v(q)).
 * 
 * NO COMMENTS ALLOWED.
 *
 */
void PW_IonicPotential::read_dft_psp()
{
  dft_text_FILE *potfile,*vqfile,*flqfile;
  char line[DFT_LINE_LEN],
       vqfilename[DFT_FILENAME_LEN],
       flqfilename[DFT_FILENAME_LEN];

  /* Open the potential file for current species */
  if ( (potfile = dft_text_fopen(potfilename,"r")) == NULL)
    die("Can't open potential file for reading.\n");

  /*
   * Read gridding information for local pseudopotential and file name
   * containing atomic potential in q-space
   */
  do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
  sscanf(line,"%d %lg %s",&ngrid_loc,&dq_loc,vqfilename);
  dft_log("Reading potential file %s\n",potfilename);
  dft_log("Local pseudopotential: ngrid = %d   dq = %g  file = %s\n",
	    ngrid_loc,dq_loc,vqfilename);
  dft_log_flush();

  /* Allocate space to store local pseudopotential on grid */
  V_loc = (real *)mymalloc(sizeof(real)*ngrid_loc,
				    "V_loc","Speciesinfo::read_dft_psp()");
      
  /* Read the file containg local pseudopotential in q-space: the
   * format is ngrid triplets of real values (q,V(q),V'(q)) one after the
   * other in ASCII text format. */
  if ( (vqfile = dft_text_fopen(vqfilename,"r")) == (dft_text_FILE *)0)
    die("Can't open vqfile for reading.\n");
  real r; // not used... just to read the file below
  for (int n=0; n < ngrid_loc; n++)
    dft_text_fscanf(vqfile,"%lg %lg",&r,&(V_loc[n]));
  dft_text_fclose(vqfile);

  /*
   * Read non-local pseudopotential information:
   *
   * first read nlm on a line, the number of (l,m) states described
   */
  do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
  sscanf(line,"%d",&nlm);
  dft_log("Non-local pseudopotential(s): nlm = %d\n", nlm);
  dft_log_flush();

  /* If there are non-local pseudopots to read, read them! */
  if (nlm != 0 ){

    /* Allocate space for l,m,ng,M,ngrid_nl,dq_nl,
     * and flq/flqprime at each (l,m) */
    l = (int *)mymalloc(sizeof(int)*nlm,"l","Speciesinfo::read_dft_psp()");
    m = (int *)mymalloc(sizeof(int)*nlm,"m","Speciesinfo::read_dft_psp()");
    ngamma = (int *)mymalloc(sizeof(int)*nlm,
			     "ngamma","Speciesinfo::read_dft_psp()");
    M = (Matrix *)mymalloc(sizeof(Matrix)*nlm,
			   "M","Speciesinfo::read_dft_psp()");
    ngrid_nl = (int **)mymalloc(sizeof(int *)*nlm,
			      "ngrid_nl","Speciesinfo::read_dft_psp()");
    dq_nl = (real **)mymalloc(sizeof(real *)*nlm,
			      "dq_nl","Speciesinfo::read_dft_psp()");
    flq = (real ***)mymalloc(sizeof(real **)*nlm,
			     "flq","Speciesinfo::read_dft_psp()");
    
    /* Loop over the (l,m) states */
    for (int lm = 0; lm < nlm; lm++)
      {
	/* Read l, m , and ngamma for this (l,m) state */
	do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
	sscanf(line,"%d %d %d",&(l[lm]),&(m[lm]),&(ngamma[lm]));
	int ng = ngamma[lm];
	dft_log("lm = %d:\n  l = %d   m = %d   ngamma = %d\n",
		lm,l[lm],m[lm],ng);
	dft_log_flush();
	
	/* Allocate space for ngrid_nl[lm], dq[lm], and flq/flqprime[lm]
	 * at each gamma */ 
	ngrid_nl[lm] = (int *)mymalloc(sizeof(int)*ng,
				       "ngrid_nl[lm]",
				       "Speciesinfo::read_dft_psp()");
	dq_nl[lm] = (real *)mymalloc(sizeof(real)*ng,
				     "dq_nl[lm]",
				     "Speciesinfo::read_dft_psp()");
	flq[lm] = (real **)mymalloc(sizeof(real *)*ng,
				    "flq[lm]","Speciesinfo::read_dft_psp()");
	
	/* Allocate and read the ng x ng M[lm] matrix */
	M[lm].init(ng,ng);
	int gamma;
	for (gamma = 0; gamma < ng; gamma++)
	  for (int gammap = 0; gammap < ng; gammap++)
	    {
	    dft_text_fscanf(potfile,"%lg",&r);
	    M[lm](gamma,gammap) = (scalar)r;
	    }
	dft_log("M = \n");
	for (gamma = 0; gamma < ng; gamma++)
	  {
	    dft_log("[ ");
	    for (int gammap = 0; gammap < ng; gammap++)
	      dft_log("%g ",REAL(M[lm](gamma,gammap)));
	    dft_log(" ]\n");
	  }
	dft_log_flush();
	
	/* Skip to next line after last M-entry */
	dft_text_fgets(line,DFT_LINE_LEN,potfile);
	
	/* For each gamma, read the gridding info and flq filename */
	for (gamma = 0; gamma < ng; gamma++)
	  {
	    /* Read gridding information and flq file name for 
	     * the current [lm][gamma] pseudopotential */
	    do
	      { dft_text_fgets(line,DFT_LINE_LEN,potfile); }
	    while(line[0] == '#');
	    sscanf(line,"%d %lg %s",
		   &(ngrid_nl[lm][gamma]),
		   &(dq_nl[lm][gamma]),
		   flqfilename);
	    int ng_nl = ngrid_nl[lm][gamma];
	    dft_log("  gamma = %d:  ngrid = %d  dq = %g  file = %s\n",
		    gamma,ng_nl,dq_nl[lm][gamma],flqfilename);
	    dft_log_flush();
	    
	    /* Allocate space to store the gridded flq[lm][gamma] */
	    flq[lm][gamma] = (real *)mymalloc(sizeof(real)*ng_nl,
					      "flq[lm][gamma]",
					      "Speciesinfo::read_dft_psp()");
	    
	    /* Open flq file for flq[lm][gamma] and read it:
	     * the format is ngrid doublets (q,flq(q)) */
	    if ( (flqfile = dft_text_fopen(flqfilename,"r")) ==
		 (dft_text_FILE *)0)
	      die("Can't open flq file for reading.\n");
	    for (int n=0; n < ng_nl; n++)
	      dft_text_fscanf(flqfile,"%lg %lg",&r,&(flq[lm][gamma][n]));
	    dft_text_fclose(flqfile);
	  }
      } /* of lm: loop on (l,m) states */
  }

  /* if there are projectors to read, read them! */
  if( projectorp == 1){
    /*
     * Read projector information:
     *
     * first read nl on a line, the number of l states described
     */
    do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
    sscanf(line,"%d",&p_nl);
    dft_log("Projector(s): p_nl = %d\n", p_nl);
    dft_log_flush();
    
    /* If there are no projectors to read, we're done! */
    if (p_nl==0)
      return;
    
    /* Allocate space for p_l,p_np,p_ngrid_nl,p_dq_nl,
     * and flq/flqprime at each l */
    p_l = (int *)mymalloc(sizeof(int)*p_nl,"p_l","Speciesinfo::read_dft_psp()");
    p_np = (int *)mymalloc(sizeof(int)*p_nl,
			   "p_np","Speciesinfo::read_dft_psp()");
    p_ngrid_nl = (int **)mymalloc(sizeof(int *)*p_nl,
				  "p_ngrid_nl","Speciesinfo::read_dft_psp()");
    p_dq_nl = (real **)mymalloc(sizeof(real *)*p_nl,
				"p_dq_nl","Speciesinfo::read_dft_psp()");
    p_flq = (real ***)mymalloc(sizeof(real **)*p_nl,
			       "p_flq","Speciesinfo::read_dft_psp()");
    
    /* Loop over the l states */
    for (int ll = 0; ll < p_nl; ll++){
      /* Read ll , and np for this ll state */
      do { dft_text_fgets(line,DFT_LINE_LEN,potfile); } while(line[0] == '#');
      sscanf(line,"%d %d",&(p_l[ll]),&(p_np[ll]));
      int np = p_np[ll];
      dft_log("\n l = %d   p_np = %d\n", p_l[ll],np);
      dft_log_flush();
      
      /* Allocate space for p_ngrid_nl[ll], p_dq[ll], and p_flq[ll]
       * at each np */ 
      p_ngrid_nl[ll] = (int *)mymalloc(sizeof(int)*np,
				       "p_ngrid_nl[ll]",
				       "Speciesinfo::read_dft_psp()");
      p_dq_nl[ll] = (real *)mymalloc(sizeof(real)*np,
				     "p_dq_nl[ll]",
				     "Speciesinfo::read_dft_psp()");
      p_flq[ll] = (real **)mymalloc(sizeof(real *)*np,
				    "p_flq[ll]","Speciesinfo::read_dft_psp()");
      
      
      /* For each p, read the gridding info and flq filename */
      for (int p = 0; p < np; p++){
	/* Read gridding information and p_flq file name for 
	 * the current [ll][p] projector */
	do
	  { dft_text_fgets(line,DFT_LINE_LEN,potfile); }
	while(line[0] == '#');
	sscanf(line,"%d %lg %s",
	       &(p_ngrid_nl[ll][p]),
	       &(p_dq_nl[ll][p]),
	       flqfilename);
	int p_ng_nl = p_ngrid_nl[ll][p];
	dft_log("  p = %d:  p_ngrid = %d  p_dq = %g  p_file = %s\n",
		p,p_ng_nl,p_dq_nl[ll][p],flqfilename);
	dft_log_flush();
	
	/* Allocate space to store the gridded flq[lm][gamma] */
	p_flq[ll][p] = (real *)mymalloc(sizeof(real)*p_ng_nl,
					"p_flq[ll][p]",
					"Speciesinfo::read_dft_psp()");
	
	/* Open p_flq file for p_flq[ll][p] and read it:
	 * the format is ngrid doublets (q,flq(q)) */
	if ( (flqfile = dft_text_fopen(flqfilename,"r")) ==
	     (dft_text_FILE *)0)
	  die("Can't open p_flq file for reading.\n");
	for (int n=0; n < p_ng_nl; n++)
	  dft_text_fscanf(flqfile,"%lg %lg",&r,&(p_flq[ll][p][n]));
	dft_text_fclose(flqfile);
      }
    } /* of ll: loop on l states */
  }
}


//
// Read pulay stuff from file (if needed)
//
/***************************************************************************
 *
 * (5) The format of the <pulayfile> is:
 *
 *         [ # comment ]
 *         NEcuts  [ # comment ]
 *         Ecut1      dEperNatoms_dNGperVol1  [ # comment]
 *         Ecut2      dEperNatoms_dNGperVol2  [ # comment]
 *          ...              ...
 *
 * -> NEcuts: the number of pairs of (Ecut,dE...perVol) provided
 * -> Ecut: planwave cutoff energy in Hartrees
 * -> dEperNatoms_dNGperVol: exactly what is says (NG=# of planewaves)!
 *
 ***************************************************************************
 */

void PW_IonicPotential::setup_pulay(Everything &e)
{
  real einfoEcut = e.basis_spec.Ecut;

  /* If Pulay correction filename is "none", then set the
   * Pulay correction to zero and skip reading the file 
   */
  if ( strcmp(pulayfilename,"none") == 0 ) {
    dft_log("Pulay file is specified as 'none':\n");
    dft_log("setting Pulay correction to zero.\n");
    dEperNatoms_dNGperVol = 0.0;
  } else {
    dft_text_FILE *pulay_file;

    /* Read file of Pulay correction factors for the species 
     */

    if ( (pulay_file = dft_text_fopen(pulayfilename,"r")) == 0)
      die("Can't open pulay file for reading.\n");
    dft_log("Reading pulay file %s\n",pulayfilename);

    /*
     * Read number of energy cutoff values listed in file
     */

    int NEcuts,n,flag;
    real Ecut;
    char line[DFT_LINE_LEN];

    do
      { dft_text_fgets(line,DFT_LINE_LEN,pulay_file); }
    while (line[0] == '#');
    sscanf(line,"%d",&NEcuts);
    dft_log("NEcuts=%d\n",NEcuts);
    flag = 0;
    for (n=0; n < NEcuts; n++) {
      do
	{ dft_text_fgets(line,DFT_LINE_LEN,pulay_file); }
      while (line[0] == '#');
      sscanf(line,"%lg %lg",
	     &Ecut,&dEperNatoms_dNGperVol);
      dft_log("Ecut = %g     dEperNatoms_dNGperVol = %g\n",
		Ecut,dEperNatoms_dNGperVol);
      if (Ecut == einfoEcut) {
	dft_log("Using Ecut=%g as it agrees with Elecinfo Ecut\n",Ecut);
	flag = 1;
	break;
      }
    }
    if (flag==0)
      die("\nCan't find appropriate energy cutoff in Pulay file!!!\n\n");
  }
}

/*
 * simps  integration for  f  with step length  h  for  n points
 * n  should be odd.
 *
 * Used by   read_fhi_psp()  below.
 */
real simps(int n, real* f, real h)
{
  real result;
  int nm12 = (n-1)/2, i;

  if ( (nm12*2) != (n-1) )
    die("Function simps: n must be odd, but is not.");

  result = f[0];
  for (i = 1; i < (n-1); i += 2)
    result += 4 * f[i] + 2 * f[i+1];
  result -= f[n-1];

  result *= (h/3);
  
  return(result);
}

/*
 * Read Pseudo potential information from file
 *
 * 2. The berlin format from package  FHI98PPMD
 *
 *    URL:   http://www.fhi-berlin.mpg.de/th/fhi98md/fhi98PP/index.html
 *
 *    max_l indicates which angular momentum components are used
 *          max_l = 0,  s
 *                  1,  s, p
 *                  2,  s, p, d
 *    loc_l indicates which one is being used as the 
 *          local component.
 *
 *    The file format is:
 *     >>>>>>>( Start )>>>>>>>>>>
 *        Zv  L_max
 *          (10 lines of number for fitted psp.
 *           not used here.)
 *        ngrid  clog       #  (for each angular momentum)
 *          i  r(i)  Psi(i)  Vion(i)
 *          ...
 *          ...
 *        (ngrid of core electron density for 
 *         nonlinear core correction
 *           --not yet supported)
 *     <<<<<<<(  End  )<<<<<<<<<<
 */
void PW_IonicPotential::read_fhi_psp()
{
  dft_text_FILE *potfile;
  char line[DFT_LINE_LEN];

  int max_l = 2, loc_l = 2;
  real zv, clog; 
  int n, i, nl, lm, gamma, gammap; 
  int ll, mm, nngamma;

  // Temporary variables.
  // maximum angular momenta is 3.
  real *r[3], *rPsi[3], *Vion[3], fhi_clog[3];
  int fhi_ngrid[3];
  real *r1, *rPsi1, *V1, *V0, *fun, q, q1, arg;
  int ngrid_r = 0;

  const int N_Lines_to_Skip = 10;
  const real fourpi = 4 * M_PI;

  // only support for  -1 < loc_l <= max_l < 3
  switch (fhi_max_l) {
  case 's' : max_l = 0; break;
  case 'p' : max_l = 1; break;
  case 'd' : max_l = 2; break;
  default:
    die("Angular momenta other than s, p, and d (0,1,2) are not supported: %d",
	max_l);
  }
  switch (fhi_loc_l) {
  case 's' : loc_l = 0; break;
  case 'p' : loc_l = 1; break;
  case 'd' : loc_l = 2; break;
  default:
    die("Angular momenta other than s, p, and d (0,1,2) are not supported: %d",
	loc_l);
  }
  if (loc_l > max_l)
    die("Illegal local component: max_l = %d, loc_l = %d",max_l,loc_l);


  /* Open the potential file for current species */
  if ( (potfile = dft_text_fopen(potfilename,"r")) == (dft_text_FILE *)0)
    die("Can't open potential file for reading.\n");

  ngrid_loc = fhi_ngrid_loc;
  dq_loc    = fhi_dq_loc;
  dft_log("Reading potential file %s (FHI CPI format)\n", potfilename);
  dft_log("Local pseudopotential: ngrid = %d   dq = %g  file = %s",
	  ngrid_loc, dq_loc, potfilename);

  switch (loc_l) {
  case 0 : dft_log(" (s) \n"); break;
  case 1 : dft_log(" (p) \n"); break;
  case 2 : dft_log(" (d) \n"); break;
  default: die("Illegal loc_l value!"); break;
  }
  
  /* Allocate space to store local pseudopotential on grid */
  V_loc = (real *)mymalloc(sizeof(real)*ngrid_loc,
			   "V_loc","setup_ioninfo");


  /*
   * The first line has two numbers:
   *   Zv   number_of_angular_momentum_components
   */
  dft_text_fgets(line,DFT_LINE_LEN,potfile);
  sscanf(line,"%lg %d", &zv, &nl);
  // verify number of valence charges, and number of angular momentum
  if (zv != Z) {
    dft_log("Valence charge unmatched in PSP file %s : %5.1f != %5.1f !!",
	    potfilename,zv,Z);
    die("Fatal error in read_fhi_psp()!");
  }
  if (nl != (max_l+1)) {
    /*
     * In the FHI pseudopotentials, nl actually takes values: 1, 2, 3
     * (perhaps due to the fact that FORTRAN index starts from 1 instead of 0)
     * corresponding to our  max_l = 0, 1, 2. 
     */
    dft_log("Number of angular momentum mismatch: %d != %d !!", nl, max_l+1 );
    die("Fatal error in read_fhi_psp()!");
  }


  // there are 10 lines of data for pseudopotentials that 
  // are fitted to expressions. Skip.
  for (n = 0; n < N_Lines_to_Skip; n++) 
    dft_text_fgets(line,DFT_LINE_LEN,potfile);

  // only allow up to 3 angular momentum components.
  // for (i = 0; i < 3; i++)
  //  sp_info->cpi_ngrid[i] = 0;

  /*
   * Read in all wavefunctions and potentials.
   * All are in real space with logrithmic grid spacing.
   */
  for (ll = 0; ll <= max_l; ll++) {
    dft_text_fgets(line,DFT_LINE_LEN,potfile);
    sscanf(line, "%d %lg", &fhi_ngrid[ll], &clog);
    // clog[l] == R[l][1]/R[l][0]  is the scale of the log grid spacing.
    fhi_clog[ll] = log(clog);
    // now clog[l] = D log(R), used for integration later.
    dft_log(" l = %d:  ngrid = %d, clog = %f\n",
	    ll, fhi_ngrid[ll], fhi_clog[ll]);

    r[ll]    = (real*)mymalloc(sizeof(real)*fhi_ngrid[ll],"r1","read_fhi_psp");
    rPsi[ll] = (real*)mymalloc(sizeof(real)*fhi_ngrid[ll],"rPsi","read_fhi_psp");
    Vion[ll] = (real*)mymalloc(sizeof(real)*fhi_ngrid[ll],"Vion","read_fhi_psp");
    if (fhi_ngrid[ll] > ngrid_r) 
      ngrid_r = fhi_ngrid[ll];

    for (i = 0; i < fhi_ngrid[ll]; i++) {
      dft_text_fgets(line,DFT_LINE_LEN,potfile);
      sscanf(line,"%*d %lg %lg %lg",
	     &(r[ll][i]), 
	     &(rPsi[ll][i]), 
	     &(Vion[ll][i]));
    }
  }

  /*
   * Read core electron density, for nonlinear core correction.
   * (not yet supported)
   */


  /*
   * Convert the local potential from log-grid real space to equal-grid k-space
   *
   */
  r1    = r[loc_l]; 
  rPsi1 = rPsi[loc_l];
  V1    = Vion[loc_l];
  fun   = (real*)mymalloc(sizeof(real)* ngrid_r,"fun","read_fhi_psp");
  for (i = 0; i < fhi_ngrid[loc_l]; i++)
    fun[i] = ( V1[i]*r1[i]*r1[i] + Z*r1[i] ) * r1[i];
  V_loc[0] = simps(fhi_ngrid[loc_l],fun,fhi_clog[loc_l]) * fourpi;
  for (n = 1; n < ngrid_loc; n++) {
    q = n * dq_loc;
    for (i = 0; i < fhi_ngrid[loc_l]; i++) {
      arg = r1[i] * q;
      fun[i] = ( V1[i]*r1[i]*r1[i] + Z*r1[i])*r1[i]*sin(arg)/arg;
    }
    // do the logarithmic-grid integral.
    V_loc[n] = simps(fhi_ngrid[loc_l],fun,fhi_clog[loc_l]) * fourpi;
  }


  /*
   * Initialize nlm, l, m, ngamma, and M
   */
  nlm = 0;
  for (ll=0; ll<=max_l; ll++) {
    if (ll != loc_l)
      nlm += (2*ll + 1);
  }
  dft_log("\nNon-local pseudopotential(s): nlm = %d\n", nlm);
  
  /* If there are no non-local pseudopots to read, we're done! */
  if (nlm==0)
    return;

  l = (int *)mymalloc(sizeof(int)*nlm,"l","setup_ioninfo");
  m = (int *)mymalloc(sizeof(int)*nlm,"m","setup_ioninfo");
  ngamma = (int *)mymalloc(sizeof(int)*nlm,"ngamma","setup_ioninfo");
  M = (Matrix *)mymalloc(sizeof(Matrix)*nlm,"M","setup_ioninfo");
  ngrid_nl = (int **)mymalloc(sizeof(int *)*nlm,
				       "ngrid_nl","setup_ioninfo");
  dq_nl = (real **)mymalloc(sizeof(real *)*nlm,
				     "dq_nl","setup_ioninfo");
  flq = (real ***)mymalloc(sizeof(real **)*nlm,
				    "flq","setup_ioninfo");
  lm = 0;  // set lm counter to 0.
  for (ll=0; ll<= max_l; ll++) { // go through all l values.
    if (ll != loc_l) { // for nonlocal cases
      // Calculate Mnl = ( <phi_l| dV_l |phi_l> )^(-1)
      ngrid_r = fhi_ngrid[ll];
      r1    = &(r[ll][0]);
      V1    = &(Vion[ll][0]);
      V0    = &(Vion[loc_l][0]);
      rPsi1 = &(rPsi[ll][0]);

      for (i = 0; i < ngrid_r; i++)
	fun[i] = (V1[i] - V0[i])*r1[i]*(rPsi1[i]*rPsi1[i]);
      q = simps(ngrid_r,fun,fhi_clog[ll]);

      // for all m values of the current l value
      for (mm= ll; mm > (-ll-1); mm--) {
	l[lm] = ll;
	m[lm] = mm;
	nngamma = ngamma[lm] = 1;  // single projector for this type of pseudopotential
	M[lm].init(nngamma,nngamma);
        for (gamma = 0; gamma < nngamma; gamma++)
          for (gammap = 0; gammap < nngamma; gammap++)
	    M[lm](gamma,gammap) = 1.0/q;

	dft_log("lm = %d:\n  l = %d   m = %d   ngamma = %d\n",
		lm,ll,mm,nngamma);
	dft_log("M = \n");
	for (gamma = 0; gamma < nngamma; gamma++) {
	  dft_log("[ ");
	  for (gammap = 0; gammap < nngamma; gammap++)
	    dft_log("%g ", REAL(M[lm](gamma,gammap)));
	  dft_log(" ]\n");
	}
	
	/* Allocate space for ngrid_nl[lm], dq[lm], and flq/flqprime[lm]
	 * at each gamma */ 
	ngrid_nl[lm] = (int *)mymalloc(sizeof(int)*nngamma,
				       "ngrid_nl[lm]","setup_ioninfo");
	dq_nl[lm] = (real *)mymalloc(sizeof(real)*nngamma,
				     "dq_nl[lm]","setup_ioninfo");
	flq[lm] = (real **)mymalloc(sizeof(real *)*nngamma,
				    "flq[lm]","setup_ioninfo");
	// ngamma must be 1 in this case.
	ngrid_nl[lm][0] = fhi_ngrid_nl;
	dq_nl[lm][0]    = fhi_dq_nl;
	flq[lm][0] = (real *)mymalloc(sizeof(real)*fhi_ngrid_nl,
				      "flq[lm][gamma]","setup_ioninfo");

	lm++; // increment counter.
      }

      for (n = 0; n < fhi_ngrid_nl; n++) {
	q = n * fhi_dq_nl;
	// construct flq from FHI pseudopotential data.
	// integrate from pseudoatom wavefunctions.
	// int_d(log(r)) r^2 * dV_l(r) * u(r) * j_l(qr)
	for (i = 0; i < ngrid_r; i++) {
	  fun[i] = (V1[i] - V0[i])*(r1[i]*r1[i])*rPsi1[i];
	  arg = r1[i] * q;
	  
	  /* The special cases for evaluating  j_l(x) near x = 0
	   * need to be taken care of. 
	   * What cut off value to use and how exactly to do this 
	   * need further thoughts.
	   */
	  q1 = 0;
	  switch (ll) {
	  case 0 : // j_0(qr) = sin(qr)/qr
	    if (arg < 1e-10) // if (qr < 1e-10), j_0(qr) == 1
	      q1 = 1;
	    else
	      q1 = (sin(arg)/arg);
	    break;
	  case 1 : // j_1(qr) = sin(qr)/(qr)^2 - cos(qr)/qr
	    if (arg < 1e-10) // j_1(x) = x/3 - x^3/30 + ...
	      q1 = arg/3;
	    else
	      q1 = (sin(arg)/arg - cos(arg))/arg;
	    break;
	  case 2 : // j_2(qr) = (3/(qr)^3-1/qr)sin(qr) - 3cos(qr)/(qr)^2
	    if (arg < 1e-10) // j_2(x) = x^2/15 - x^4/210 + ...
	      q1 = arg*arg /15;
	    else {
	      q1 = (3/(arg*arg)-1)*sin(arg)-3*cos(arg)/arg;
	      q1 = q1/arg;
	    }
	    break;
	  default :
	    die("Illegal l value!");
	    break;
	  }
	  fun[i] *= q1;
	}

	q1 = simps(ngrid_r,fun,fhi_clog[ll]);
	for (mm = 0; mm < (2*ll+1); mm++)
	  flq[lm-mm-1][0][n] = q1;
      }
    }
  }

  // now free up temporary memory.
  myfree(fun);
  for (ll = 0; ll <= max_l; ll++) {
    myfree(r[ll]);
    myfree(rPsi[ll]);
    myfree(Vion[ll]);
  }
}