MultiSource/Benchmarks/DOE-ProxyApps-C/CoMD/eam.c

/// \file
/// Compute forces for the Embedded Atom Model (EAM).
///
/// The Embedded Atom Model (EAM) is a widely used model of atomic
/// interactions in simple metals.
/// 
/// http://en.wikipedia.org/wiki/Embedded_atom_model
///
/// In the EAM, the total potential energy is written as a sum of a pair
/// potential and the embedding energy, F:
///
/// \f[
///   U = \sum_{ij} \varphi(r_{ij}) + \sum_i F({\bar\rho_i})
/// \f]
///
/// The pair potential \f$\varphi_{ij}\f$ is a two-body inter-atomic
/// potential, similar to the Lennard-Jones potential, and
/// \f$F(\bar\rho)\f$ is interpreted as the energy required to embed an
/// atom in an electron field with density \f$\bar\rho\f$.  The local
/// electon density at site i is calulated by summing the "effective
/// electron density" due to all neighbors of atom i:
///
/// \f[
/// \bar\rho_i = \sum_j \rho_j(r_{ij})
/// \f]
///
/// The force on atom i, \f${\bf F}_i\f$ is given by
///
/// \f{eqnarray*}{
///   {\bf F}_i & = & -\nabla_i \sum_{jk} U(r_{jk})\\
///       & = & - \sum_j\left\{
///                  \varphi'(r_{ij}) +
///                  [F'(\bar\rho_i) + F'(\bar\rho_j)]\rho'(r_{ij})
///                \right\} \hat{r}_{ij}
/// \f}
///
/// where primes indicate the derivative of a function with respect to
/// its argument and \f$\hat{r}_{ij}\f$ is a unit vector in the
/// direction from atom i to atom j.
///
/// The form of this force expression has two significant consequences.
/// First, unlike with a simple pair potential, it is not possible to
/// compute the potential energy and the forces on the atoms in a single
/// loop over the pairs.  The terms involving \f$ F'(\bar\rho) \f$
/// cannot be calculated until \f$ \bar\rho \f$ is known, but
/// calculating \f$ \bar\rho \f$ requires a loop over the pairs.  Hence
/// the EAM force routine contains three loops.
///
///   -# Loop over all pairs, compute the two-body
///   interaction and the electron density at each atom
///   -# Loop over all atoms, compute the embedding energy and its
///   derivative for each atom
///   -# Loop over all pairs, compute the embedding
///   energy contribution to the force and add to the two-body force
///
/// The second loop over pairs doubles the data motion requirement
/// relative to a simple pair potential.
///
/// The second consequence of the force expression is that computing the
/// forces on all atoms requires additional communication beyond the
/// coordinates of all remote atoms within the cutoff distance.  This is
/// again because of the terms involving \f$ F'(\bar\rho_j) \f$.  If
/// atom j is a remote atom, the local task cannot compute \f$
/// \bar\rho_j \f$.  (Such a calculation would require all the neighbors
/// of atom j, some of which can be up to 2 times the cutoff distance
/// away from a local atom---outside the typical halo exchange range.)
///
/// To obtain the needed remote density we introduce a second halo
/// exchange after loop number 2 to communicate \f$ F'(\bar\rho) \f$ for
/// remote atoms.  This provides the data we need to complete the third
/// loop, but at the cost of introducing a communication operation in
/// the middle of the force routine.
///
/// At least two alternate methods can be used to deal with the remote
/// density problem.  One possibility is to extend the halo exchange
/// radius for the atom exchange to twice the potential cutoff distance.
/// This is likely undesirable due to large increase in communication
/// volume.  The other possibility is to accumulate partial force terms
/// on the tasks where they can be computed.  In this method, tasks will
/// compute force contributions for remote atoms, then communicate the
/// partial forces at the end of the halo exchange.  This method has the
/// advantage that the communication is deffered until after the force
/// loops, but the disadvantage that three times as much data needs to
/// be set (three components of the force vector instead of a single
/// scalar \f$ F'(\bar\rho) \f$.


#include "eam.h"

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <assert.h>

#include "constants.h"
#include "memUtils.h"
#include "parallel.h"
#include "linkCells.h"
#include "CoMDTypes.h"
#include "performanceTimers.h"
#include "haloExchange.h"

#define MAX(A,B) ((A) > (B) ? (A) : (B))

/// Handles interpolation of tabular data.
///
/// \see initInterpolationObject
/// \see interpolate
typedef struct InterpolationObjectSt 
{
   int n;          //!< the number of values in the table
   real_t x0;      //!< the starting ordinate range
   real_t invDx;   //!< the inverse of the table spacing
   real_t* values; //!< the abscissa values
} InterpolationObject;

/// Derived struct for an EAM potential.
/// Uses table lookups for function evaluation.
/// Polymorphic with BasePotential.
/// \see BasePotential
typedef struct EamPotentialSt 
{
   real_t cutoff;          //!< potential cutoff distance in Angstroms
   real_t mass;            //!< mass of atoms in intenal units
   real_t lat;             //!< lattice spacing (angs) of unit cell
   char latticeType[8];    //!< lattice type, e.g. FCC, BCC, etc.
   char  name[3];	   //!< element name
   int	 atomicNo;	   //!< atomic number  
   int  (*force)(SimFlat* s); //!< function pointer to force routine
   void (*print)(FILE* file, BasePotential* pot);
   void (*destroy)(BasePotential** pot); //!< destruction of the potential
   InterpolationObject* phi;  //!< Pair energy
   InterpolationObject* rho;  //!< Electron Density
   InterpolationObject* f;    //!< Embedding Energy

   real_t* rhobar;        //!< per atom storage for rhobar
   real_t* dfEmbed;       //!< per atom storage for derivative of Embedding
   HaloExchange* forceExchange;
   ForceExchangeData* forceExchangeData;
} EamPotential;

// EAM functionality
static int eamForce(SimFlat* s);
static void eamPrint(FILE* file, BasePotential* pot);
static void eamDestroy(BasePotential** pot); 
static void eamBcastPotential(EamPotential* pot);


// Table interpolation functionality
static InterpolationObject* initInterpolationObject(
   int n, real_t x0, real_t dx, real_t* data);
static void destroyInterpolationObject(InterpolationObject** table);
static void interpolate(InterpolationObject* table, real_t r, real_t* f, real_t* df);
static void bcastInterpolationObject(InterpolationObject** table);
static void printTableData(InterpolationObject* table, const char* fileName);


// Read potential tables from files.
static void eamReadSetfl(EamPotential* pot, const char* dir, const char* potName);
static void eamReadFuncfl(EamPotential* pot, const char* dir, const char* potName);
static void fileNotFound(const char* callSite, const char* filename);
static void notAlloyReady(const char* callSite);
static void typeNotSupported(const char* callSite, const char* type);


/// Allocate and initialize the EAM potential data structure.
///
/// \param [in] dir   The directory in which potential table files are found.
/// \param [in] file  The name of the potential table file.
/// \param [in] type  The file format of the potential file (setfl or funcfl).
BasePotential* initEamPot(const char* dir, const char* file, const char* type)
{
   EamPotential* pot = comdMalloc(sizeof(EamPotential));
   assert(pot);
   pot->force = eamForce;
   pot->print = eamPrint;
   pot->destroy = eamDestroy;
   pot->phi = NULL;
   pot->rho = NULL;
   pot->f   = NULL;

   // Initialization of the next three items requires information about
   // the parallel decomposition and link cells that isn't available
   // with the potential is initialized.  Hence, we defer their
   // initialization until the first time we call the force routine.
   pot->dfEmbed = NULL;
   pot->rhobar  = NULL;
   pot->forceExchange = NULL;

   if (getMyRank() == 0)
   {
      if (strcmp(type, "setfl" ) == 0)
         eamReadSetfl(pot, dir, file);
      else if (strcmp(type,"funcfl") == 0)
         eamReadFuncfl(pot, dir, file);
      else
         typeNotSupported("initEamPot", type);
   }
   eamBcastPotential(pot);
   
   return (BasePotential*) pot;
}

/// Calculate potential energy and forces for the EAM potential.
///
/// Three steps are required:
///
///   -# Loop over all atoms and their neighbors, compute the two-body
///   interaction and the electron density at each atom
///   -# Loop over all atoms, compute the embedding energy and its
///   derivative for each atom
///   -# Loop over all atoms and their neighbors, compute the embedding
///   energy contribution to the force and add to the two-body force
/// 
int eamForce(SimFlat* s)
{
   EamPotential* pot = (EamPotential*) s->pot;
   assert(pot);

   // set up halo exchange and internal storage on first call to forces.
   if (pot->forceExchange == NULL)
   {
      int maxTotalAtoms = MAXATOMS*s->boxes->nTotalBoxes;
      pot->dfEmbed = comdMalloc(maxTotalAtoms*sizeof(real_t));
      pot->rhobar  = comdMalloc(maxTotalAtoms*sizeof(real_t));
      pot->forceExchange = initForceHaloExchange(s->domain, s->boxes);
      pot->forceExchangeData = comdMalloc(sizeof(ForceExchangeData));
      pot->forceExchangeData->dfEmbed = pot->dfEmbed;
      pot->forceExchangeData->boxes = s->boxes;
   }
   
   real_t rCut2 = pot->cutoff*pot->cutoff;

   // zero forces / energy / rho /rhoprime
   real_t etot = 0.0;
   memset(s->atoms->f,  0, s->boxes->nTotalBoxes*MAXATOMS*sizeof(real3));
   memset(s->atoms->U,  0, s->boxes->nTotalBoxes*MAXATOMS*sizeof(real_t));
   memset(pot->dfEmbed, 0, s->boxes->nTotalBoxes*MAXATOMS*sizeof(real_t));
   memset(pot->rhobar,  0, s->boxes->nTotalBoxes*MAXATOMS*sizeof(real_t));

   int nbrBoxes[27];
   // loop over local boxes
   for (int iBox=0; iBox<s->boxes->nLocalBoxes; iBox++)
   {
      int nIBox = s->boxes->nAtoms[iBox];
      int nNbrBoxes = getNeighborBoxes(s->boxes, iBox, nbrBoxes);
      // loop over neighbor boxes of iBox (some may be halo boxes)
      for (int jTmp=0; jTmp<nNbrBoxes; jTmp++)
      {
         int jBox = nbrBoxes[jTmp];
         if (jBox < iBox ) continue;

         int nJBox = s->boxes->nAtoms[jBox];
         // loop over atoms in iBox
         for (int iOff=MAXATOMS*iBox,ii=0; ii<nIBox; ii++,iOff++)
         {
            // loop over atoms in jBox
            for (int jOff=MAXATOMS*jBox,ij=0; ij<nJBox; ij++,jOff++)
            {
               if ( (iBox==jBox) &&(ij <= ii) ) continue;

               double r2 = 0.0;
               real3 dr;
               for (int k=0; k<3; k++)
               {
                  dr[k]=s->atoms->r[iOff][k]-s->atoms->r[jOff][k];
                  r2+=dr[k]*dr[k];
               }
               if(r2>rCut2) continue;

               double r = sqrt(r2);

               real_t phiTmp, dPhi, rhoTmp, dRho;
               interpolate(pot->phi, r, &phiTmp, &dPhi);
               interpolate(pot->rho, r, &rhoTmp, &dRho);

               for (int k=0; k<3; k++)
               {
                  s->atoms->f[iOff][k] -= dPhi*dr[k]/r;
                  s->atoms->f[jOff][k] += dPhi*dr[k]/r;
               }

               // update energy terms
               // calculate energy contribution based on whether
               // the neighbor box is local or remote
               if (jBox < s->boxes->nLocalBoxes)
                  etot += phiTmp;
               else
                  etot += 0.5*phiTmp;

               s->atoms->U[iOff] += 0.5*phiTmp;
               s->atoms->U[jOff] += 0.5*phiTmp;

               // accumulate rhobar for each atom
               pot->rhobar[iOff] += rhoTmp;
               pot->rhobar[jOff] += rhoTmp;

            } // loop over atoms in jBox
         } // loop over atoms in iBox
      } // loop over neighbor boxes
   } // loop over local boxes

   // Compute Embedding Energy
   // loop over all local boxes
   for (int iBox=0; iBox<s->boxes->nLocalBoxes; iBox++)
   {
      int iOff;
      int nIBox =  s->boxes->nAtoms[iBox];

      // loop over atoms in iBox
      for (int iOff=MAXATOMS*iBox,ii=0; ii<nIBox; ii++,iOff++)
      {
         real_t fEmbed, dfEmbed;
         interpolate(pot->f, pot->rhobar[iOff], &fEmbed, &dfEmbed);
         pot->dfEmbed[iOff] = dfEmbed; // save derivative for halo exchange
         etot += fEmbed; 
         s->atoms->U[iOff] += fEmbed;
      }
   }

   // exchange derivative of the embedding energy with repsect to rhobar
   startTimer(eamHaloTimer);
   haloExchange(pot->forceExchange, pot->forceExchangeData);
   stopTimer(eamHaloTimer);

   // third pass
   // loop over local boxes
   for (int iBox=0; iBox<s->boxes->nLocalBoxes; iBox++)
   {
      int nIBox =  s->boxes->nAtoms[iBox];
      int nNbrBoxes = getNeighborBoxes(s->boxes, iBox, nbrBoxes);
      // loop over neighbor boxes of iBox (some may be halo boxes)
      for (int jTmp=0; jTmp<nNbrBoxes; jTmp++)
      {
         int jBox = nbrBoxes[jTmp];
         if(jBox < iBox) continue;

         int nJBox = s->boxes->nAtoms[jBox];
         // loop over atoms in iBox
         for (int iOff=MAXATOMS*iBox,ii=0; ii<nIBox; ii++,iOff++)
         {
            // loop over atoms in jBox
            for (int jOff=MAXATOMS*jBox,ij=0; ij<nJBox; ij++,jOff++)
            { 
               if ((iBox==jBox) && (ij <= ii))  continue;

               double r2 = 0.0;
               real3 dr;
               for (int k=0; k<3; k++)
               {
                  dr[k]=s->atoms->r[iOff][k]-s->atoms->r[jOff][k];
                  r2+=dr[k]*dr[k];
               }
               if(r2>=rCut2) continue;

               real_t r = sqrt(r2);

               real_t rhoTmp, dRho;
               interpolate(pot->rho, r, &rhoTmp, &dRho);

               for (int k=0; k<3; k++)
               {
                  s->atoms->f[iOff][k] -= (pot->dfEmbed[iOff]+pot->dfEmbed[jOff])*dRho*dr[k]/r;
                  s->atoms->f[jOff][k] += (pot->dfEmbed[iOff]+pot->dfEmbed[jOff])*dRho*dr[k]/r;
               }

            } // loop over atoms in jBox
         } // loop over atoms in iBox
      } // loop over neighbor boxes
   } // loop over local boxes

   s->ePotential = (real_t) etot;

   return 0;
}

void eamPrint(FILE* file, BasePotential* pot)
{
   EamPotential *eamPot = (EamPotential*) pot;
   fprintf(file, "  Potential type  : EAM\n");
   fprintf(file, "  Species name    : %s\n", eamPot->name);
   fprintf(file, "  Atomic number   : %d\n", eamPot->atomicNo);
   fprintf(file, "  Mass            : "FMT1" amu\n", eamPot->mass/amuToInternalMass); // print in amu
   fprintf(file, "  Lattice type    : %s\n", eamPot->latticeType);
   fprintf(file, "  Lattice spacing : "FMT1" Angstroms\n", eamPot->lat);
   fprintf(file, "  Cutoff          : "FMT1" Angstroms\n", eamPot->cutoff);
}

void eamDestroy(BasePotential** pPot)
{
   if ( ! pPot ) return;
   EamPotential* pot = *(EamPotential**)pPot;
   if ( ! pot ) return;
   destroyInterpolationObject(&(pot->phi));
   destroyInterpolationObject(&(pot->rho));
   destroyInterpolationObject(&(pot->f));
   destroyHaloExchange(&(pot->forceExchange));
   comdFree(pot);
   *pPot = NULL;

   return;
}

/// Broadcasts an EamPotential from rank 0 to all other ranks.
/// If the table coefficients are read from a file only rank 0 does the
/// read.  Hence we need to broadcast the potential to all other ranks.
void eamBcastPotential(EamPotential* pot)
{
   assert(pot);
   
   struct 
   {
      real_t cutoff, mass, lat;
      char latticeType[8];
      char name[3];
      int atomicNo;
   } buf;
   if (getMyRank() == 0)
   {
      buf.cutoff   = pot->cutoff;
      buf.mass     = pot->mass;
      buf.lat      = pot->lat;
      buf.atomicNo = pot->atomicNo;
      strcpy(buf.latticeType, pot->latticeType);
      strcpy(buf.name, pot->name);
   }
   bcastParallel(&buf, sizeof(buf), 0);
   pot->cutoff   = buf.cutoff;
   pot->mass     = buf.mass;
   pot->lat      = buf.lat;
   pot->atomicNo = buf.atomicNo;
   strcpy(pot->latticeType, buf.latticeType);
   strcpy(pot->name, buf.name);

   bcastInterpolationObject(&pot->phi);
   bcastInterpolationObject(&pot->rho);
   bcastInterpolationObject(&pot->f);
}

/// Builds a structure to store interpolation data for a tabular
/// function.  Interpolation must be supported on the range
/// \f$[x_0, x_n]\f$, where \f$x_n = n*dx\f$.
///
/// \see interpolate
/// \see bcastInterpolationObject
/// \see destroyInterpolationObject
///
/// \param [in] n    number of values in the table.
/// \param [in] x0   minimum ordinate value of the table.
/// \param [in] dx   spacing of the ordinate values.
/// \param [in] data abscissa values.  An array of size n. 
InterpolationObject* initInterpolationObject(
   int n, real_t x0, real_t dx, real_t* data)
{
   InterpolationObject* table =
      (InterpolationObject *)comdMalloc(sizeof(InterpolationObject)) ;
   assert(table);

   table->values = (real_t*)comdCalloc(1, (n+3)*sizeof(real_t));
   assert(table->values);

   table->values++; 
   table->n = n;
   table->invDx = 1.0/dx;
   table->x0 = x0;

   for (int ii=0; ii<n; ++ii)
      table->values[ii] = data[ii];
   
   table->values[-1] = table->values[0];
   table->values[n+1] = table->values[n] = table->values[n-1];

   return table;
}

void destroyInterpolationObject(InterpolationObject** a)
{
   if ( ! a ) return;
   if ( ! *a ) return;
   if ( (*a)->values)
   {
      (*a)->values--;
      comdFree((*a)->values);
   }
   comdFree(*a);
   *a = NULL;

   return;
}

/// Interpolate a table to determine f(r) and its derivative f'(r).
///
/// The forces on the particle are much more sensitive to the derivative
/// of the potential than on the potential itself.  It is therefore
/// absolutely essential that the interpolated derivatives are smooth
/// and continuous.  This function uses simple quadratic interpolation
/// to find f(r).  Since quadric interpolants don't have smooth
/// derivatives, f'(r) is computed using a 4 point finite difference
/// stencil.
///
/// Interpolation is used heavily by the EAM force routine so this
/// function is a potential performance hot spot.  Feel free to
/// reimplement this function (and initInterpolationObject if necessay)
/// with any higher performing implementation of interpolation, as long
/// as the alternate implmentation that has the required smoothness
/// properties.  Cubic splines are one common alternate choice.
///
/// \param [in] table Interpolation table.
/// \param [in] r Point where function value is needed.
/// \param [out] f The interpolated value of f(r).
/// \param [out] df The interpolated value of df(r)/dr.
void interpolate(InterpolationObject* table, real_t r, real_t* f, real_t* df)
{
   const real_t* tt = table->values; // alias
   
   if ( r < table->x0 ) r = table->x0;

   r = (r-table->x0)*(table->invDx) ;
   int ii = (int)floor(r);
   if (ii > table->n)
   {
      ii = table->n;
      r = table->n / table->invDx;
   }
   // reset r to fractional distance
   r = r - floor(r);

   real_t g1 = tt[ii+1] - tt[ii-1];
   real_t g2 = tt[ii+2] - tt[ii];

   *f = tt[ii] + 0.5*r*(g1 + r*(tt[ii+1] + tt[ii-1] - 2.0*tt[ii]) );

   *df = 0.5*(g1 + r*(g2-g1))*table->invDx;
}

/// Broadcasts an InterpolationObject from rank 0 to all other ranks.
///
/// It is commonly the case that the data needed to create the
/// interpolation table is available on only one task (for example, only
/// one task has read the data from a file).  Broadcasting the table
/// eliminates the need to put broadcast code in multiple table readers.
///
/// \see eamBcastPotential
void bcastInterpolationObject(InterpolationObject** table)
{
   struct
   {
      int n;
      real_t x0, invDx;
   } buf;

   if (getMyRank() == 0)
   {
      buf.n     = (*table)->n;
      buf.x0    = (*table)->x0;
      buf.invDx = (*table)->invDx;
   }
   bcastParallel(&buf, sizeof(buf), 0);

   if (getMyRank() != 0)
   {
      assert(*table == NULL);
      *table = comdMalloc(sizeof(InterpolationObject));
      (*table)->n      = buf.n;
      (*table)->x0     = buf.x0;
      (*table)->invDx  = buf.invDx;
      (*table)->values = comdMalloc(sizeof(real_t) * (buf.n+3) );
      (*table)->values++;
   }
   
   int valuesSize = sizeof(real_t) * ((*table)->n+3);
   bcastParallel((*table)->values-1, valuesSize, 0);
}

void printTableData(InterpolationObject* table, const char* fileName)
{
   if (!printRank()) return;

   FILE* potData;
   potData = fopen(fileName,"w");
   real_t dR = 1.0/table->invDx;
   for (int i = 0; i<table->n; i++)
   {
      real_t r = table->x0+i*dR;
      fprintf(potData, "%d %e %e\n", i, r, table->values[i]);
   }
   fclose(potData);
}
   
/// Reads potential data from a setfl file and populates
/// corresponding members and InterpolationObjects in an EamPotential.
///
/// setfl is a file format for tabulated potential functions used by
/// the original EAM code DYNAMO.  A setfl file contains EAM
/// potentials for multiple elements.
///
/// The contents of a setfl file are:
///
/// | Line Num | Description
/// | :------: | :----------
/// | 1 - 3    | comments
/// | 4        | ntypes type1 type2 ... typen
/// | 5        | nrho     drho     nr   dr   rcutoff
/// | F, rho   | Following line 5 there is a block for each atom type with F, and rho.
/// | b1       | ielem(i)   amass(i)     latConst(i)    latType(i)
/// | b2       | embedding function values F(rhobar) starting at rhobar=0
/// |   ...    | (nrho values. Multiple values per line allowed.)
/// | bn       | electron density, starting at r=0
/// |   ...    | (nr values. Multiple values per line allowed.)
/// | repeat   | Return to b1 for each atom type.
/// | phi      | phi_ij for (1,1), (2,1), (2,2), (3,1), (3,2), (3,3), (4,1), ..., 
/// | p1       | pair potential between type i and type j, starting at r=0
/// |   ...    | (nr values. Multiple values per line allowed.)
/// | repeat   | Return to p1 for each phi_ij
///
/// Where:
///    -  ntypes        :      number of element types in the potential  
///    -  nrho          :      number of points the embedding energy F(rhobar)
///    -  drho          :      table spacing for rhobar 
///    -  nr            :      number of points for rho(r) and phi(r)
///    -  dr            :      table spacing for r in Angstroms
///    -  rcutoff       :      cut-off distance in Angstroms
///    -  ielem(i)      :      atomic number for element(i)
///    -  amass(i)      :      atomic mass for element(i) in AMU
///    -  latConst(i)   :      lattice constant for element(i) in Angstroms
///    -  latType(i)    :      lattice type for element(i)  
///
/// setfl format stores r*phi(r), so we need to converted to the pair
/// potential phi(r).  In the file, phi(r)*r is in eV*Angstroms.
/// NB: phi is not defined for r = 0
///
/// F(rhobar) is in eV.
///
void eamReadSetfl(EamPotential* pot, const char* dir, const char* potName)
{
   char tmp[4096];
   sprintf(tmp, "%s/%s", dir, potName);

   FILE* potFile = fopen(tmp, "r");
   if (potFile == NULL)
      fileNotFound("eamReadSetfl", tmp);
   
   // read the first 3 lines (comments)
   fgets(tmp, sizeof(tmp), potFile);
   fgets(tmp, sizeof(tmp), potFile);
   fgets(tmp, sizeof(tmp), potFile);

   // line 4
   fgets(tmp, sizeof(tmp), potFile);
   int nElems;
   sscanf(tmp, "%d", &nElems);
   if( nElems != 1 )
      notAlloyReady("eamReadSetfl");

   //line 5
   int nRho, nR;
   double dRho, dR, cutoff;
   //  The same cutoff is used by all alloys, NB: cutoff = nR * dR is redundant
   fgets(tmp, sizeof(tmp), potFile);
   sscanf(tmp, "%d %le %d %le %le", &nRho, &dRho, &nR, &dR, &cutoff);
   pot->cutoff = cutoff;

   // **** THIS CODE IS RESTRICTED TO ONE ELEMENT
   // Per-atom header 
   fgets(tmp, sizeof(tmp), potFile);
   int nAtomic;
   double mass, lat;
   char latticeType[8];
   sscanf(tmp, "%d %le %le %s", &nAtomic, &mass, &lat, latticeType);
   pot->atomicNo = nAtomic;
   pot->lat = lat;
   pot->mass = mass * amuToInternalMass;  // file has mass in AMU.
   strcpy(pot->latticeType, latticeType);
   
   // allocate read buffer
   int bufSize = MAX(nRho, nR);
   real_t* buf = comdMalloc(bufSize * sizeof(real_t));
   real_t x0 = 0.0;

   // Read embedding energy F(rhobar)
   for (int ii=0; ii<nRho; ++ii)
      fscanf(potFile, FMT1, buf+ii);
   pot->f = initInterpolationObject(nRho, x0, dRho, buf);

   // Read electron density rho(r)
   for (int ii=0; ii<nR; ++ii)
      fscanf(potFile, FMT1, buf+ii);
   pot->rho = initInterpolationObject(nR, x0, dR, buf);

   // Read phi(r)*r and convert to phi(r)
   for (int ii=0; ii<nR; ++ii)
      fscanf(potFile, FMT1, buf+ii);
   for (int ii=1; ii<nR; ++ii)
   {
      real_t r = x0 + ii*dR;
      buf[ii] /= r;
   }
   buf[0] = buf[1] + (buf[1] - buf[2]); // Linear interpolation to get phi[0].
   pot->phi = initInterpolationObject(nR, x0, dR, buf);

   comdFree(buf);

   // write to text file for comparison, currently commented out
/*    printPot(pot->f, "SetflDataF.txt"); */
/*    printPot(pot->rho, "SetflDataRho.txt"); */
/*    printPot(pot->phi, "SetflDataPhi.txt");  */
}

/// Reads potential data from a funcfl file and populates
/// corresponding members and InterpolationObjects in an EamPotential.
/// 
/// funcfl is a file format for tabulated potential functions used by
/// the original EAM code DYNAMO.  A funcfl file contains an EAM
/// potential for a single element.
/// 
/// The contents of a funcfl file are:
///
/// | Line Num | Description
/// | :------: | :----------
/// | 1        | comments
/// | 2        | elem amass latConstant latType
/// | 3        | nrho   drho   nr   dr    rcutoff
/// | 4        | embedding function values F(rhobar) starting at rhobar=0
/// |    ...   | (nrho values. Multiple values per line allowed.)
/// | x'       | electrostatic interation Z(r) starting at r=0
/// |    ...   | (nr values. Multiple values per line allowed.)
/// | y'       | electron density values rho(r) starting at r=0
/// |    ...   | (nr values. Multiple values per line allowed.)
///
/// Where:
///    -  elem          :   atomic number for this element
///    -  amass         :   atomic mass for this element in AMU
///    -  latConstant   :   lattice constant for this elemnent in Angstroms
///    -  lattticeType  :   lattice type for this element (e.g. FCC) 
///    -  nrho          :   number of values for the embedding function, F(rhobar)
///    -  drho          :   table spacing for rhobar
///    -  nr            :   number of values for Z(r) and rho(r)
///    -  dr            :   table spacing for r in Angstroms
///    -  rcutoff       :   potential cut-off distance in Angstroms
///
/// funcfl format stores the "electrostatic interation" Z(r).  This needs to
/// be converted to the pair potential phi(r).
/// using the formula 
/// \f[phi = Z(r) * Z(r) / r\f]
/// NB: phi is not defined for r = 0
///
/// Z(r) is in atomic units (i.e., sqrt[Hartree * bohr]) so it is
/// necesary to convert to eV.
///
/// F(rhobar) is in eV.
///
void eamReadFuncfl(EamPotential* pot, const char* dir, const char* potName)
{
   char tmp[4096];

   sprintf(tmp, "%s/%s", dir, potName);
   FILE* potFile = fopen(tmp, "r");
   if (potFile == NULL)
      fileNotFound("eamReadFuncfl", tmp);

   // line 1
   fgets(tmp, sizeof(tmp), potFile);
   char name[3];
   sscanf(tmp, "%s", name);
   strcpy(pot->name, name);

   // line 2
   int nAtomic;
   double mass, lat;
   char latticeType[8];
   fgets(tmp,sizeof(tmp),potFile);
   sscanf(tmp, "%d %le %le %s", &nAtomic, &mass, &lat, latticeType);
   pot->atomicNo = nAtomic;
   pot->lat = lat;
   pot->mass = mass*amuToInternalMass; // file has mass in AMU.
   strcpy(pot->latticeType, latticeType);

   // line 3
   int nRho, nR;
   double dRho, dR, cutoff;
   fgets(tmp,sizeof(tmp),potFile);
   sscanf(tmp, "%d %le %d %le %le", &nRho, &dRho, &nR, &dR, &cutoff);
   pot->cutoff = cutoff;
   real_t x0 = 0.0; // tables start at zero.

   // allocate read buffer
   int bufSize = MAX(nRho, nR);
   real_t* buf = comdMalloc(bufSize * sizeof(real_t));

   // read embedding energy
   for (int ii=0; ii<nRho; ++ii)
      fscanf(potFile, FMT1, buf+ii);
   pot->f = initInterpolationObject(nRho, x0, dRho, buf);

   // read Z(r) and convert to phi(r)
   for (int ii=0; ii<nR; ++ii)
      fscanf(potFile, FMT1, buf+ii);
   for (int ii=1; ii<nR; ++ii)
   {
      real_t r = x0 + ii*dR;
      buf[ii] *= buf[ii] / r;
      buf[ii] *= hartreeToEv * bohrToAngs; // convert to eV
   }
   buf[0] = buf[1] + (buf[1] - buf[2]); // linear interpolation to get phi[0].
   pot->phi = initInterpolationObject(nR, x0, dR, buf);

   // read electron density rho
   for (int ii=0; ii<nR; ++ii)
      fscanf(potFile, FMT1, buf+ii);
   pot->rho = initInterpolationObject(nR, x0, dR, buf);

   comdFree(buf);
   
/*    printPot(pot->f,   "funcflDataF.txt"); */
/*    printPot(pot->rho, "funcflDataRho.txt"); */
/*    printPot(pot->phi, "funcflDataPhi.txt"); */
}

void fileNotFound(const char* callSite, const char* filename)
{
   fprintf(screenOut,
           "%s: Can't open file %s.  Fatal Error.\n", callSite, filename);
   exit(-1);
}

void notAlloyReady(const char* callSite)
{
   fprintf(screenOut,
          "%s: CoMD 1.1 does not support alloys and cannot\n"
           "   read setfl files with multiple species.  Fatal Error.\n", callSite);
   exit(-1);
}

void typeNotSupported(const char* callSite, const char* type)
{
   fprintf(screenOut,
          "%s: Potential type %s not supported. Fatal Error.\n", callSite, type);
   exit(-1);
}