diff --git a/src/MANYBODY/Install.sh b/src/MANYBODY/Install.sh
index efa85e3aae..c0e2f4d956 100644
--- a/src/MANYBODY/Install.sh
+++ b/src/MANYBODY/Install.sh
@@ -11,6 +11,7 @@ if (test $1 = 1) then
   cp pair_eam.cpp ..
   cp pair_eam_alloy.cpp ..
   cp pair_eam_fs.cpp ..
+  cp pair_eim.cpp ..
   cp pair_sw.cpp ..
   cp pair_tersoff.cpp ..
   cp pair_tersoff_zbl.cpp ..
@@ -21,6 +22,7 @@ if (test $1 = 1) then
   cp pair_eam.h ..
   cp pair_eam_alloy.h ..
   cp pair_eam_fs.h ..
+  cp pair_eim.h ..
   cp pair_sw.h ..
   cp pair_tersoff.h ..
   cp pair_tersoff_zbl.h ..
@@ -37,6 +39,7 @@ elif (test $1 = 0) then
   rm ../pair_eam.cpp
   rm ../pair_eam_alloy.cpp
   rm ../pair_eam_fs.cpp
+  rm ../pair_eim.cpp
   rm ../pair_sw.cpp
   rm ../pair_tersoff.cpp
   rm ../pair_tersoff_zbl.cpp
@@ -47,6 +50,7 @@ elif (test $1 = 0) then
   rm ../pair_eam.h
   rm ../pair_eam_alloy.h
   rm ../pair_eam_fs.h
+  rm ../pair_eim.h
   rm ../pair_sw.h
   rm ../pair_tersoff.h
   rm ../pair_tersoff_zbl.h
diff --git a/src/MANYBODY/pair_eim.cpp b/src/MANYBODY/pair_eim.cpp
new file mode 100644
index 0000000000..68a83c1f52
--- /dev/null
+++ b/src/MANYBODY/pair_eim.cpp
@@ -0,0 +1,1172 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+/* ----------------------------------------------------------------------
+   Contributing author: Xiaowang Zhou (SNL)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_eim.h"
+#include "atom.h"
+#include "force.h"
+#include "comm.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "memory.h"
+#include "error.h"
+
+using namespace LAMMPS_NS;
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+#define MAXLINE 1024
+
+/* ---------------------------------------------------------------------- */
+
+PairEIM::PairEIM(LAMMPS *lmp) : Pair(lmp)
+{
+  single_enable = 0;
+  one_coeff = 1;
+
+  setfl = NULL;
+  nmax = 0;
+  rho = NULL;
+  fp = NULL;
+
+  nelements = 0;
+  elements = NULL;
+
+  negativity = NULL;
+  q0 = NULL;
+  cutforcesq = NULL;
+  Fij = NULL;
+  Gij = NULL;
+  phiij = NULL;
+
+  Fij_spline = NULL;
+  Gij_spline = NULL;
+  phiij_spline = NULL;
+
+  // set comm size needed by this Pair
+
+  comm_forward = 1;
+  comm_reverse = 1;
+}
+
+/* ----------------------------------------------------------------------
+   check if allocated, since class can be destructed when incomplete
+------------------------------------------------------------------------- */
+
+PairEIM::~PairEIM()
+{
+  memory->sfree(rho);
+  memory->sfree(fp);
+
+  if (allocated) {
+    memory->destroy_2d_int_array(setflag);
+    memory->destroy_2d_double_array(cutsq);
+    delete [] map;
+    memory->destroy_2d_int_array(type2Fij);
+    memory->destroy_2d_int_array(type2Gij);
+    memory->destroy_2d_int_array(type2phiij);
+  }
+
+  for (int i = 0; i < nelements; i++) delete [] elements[i];
+  delete [] elements;
+  
+  deallocate_setfl();
+
+  delete [] negativity;
+  delete [] q0;
+  memory->destroy_2d_double_array(cutforcesq);
+  memory->destroy_2d_double_array(Fij);
+  memory->destroy_2d_double_array(Gij);
+  memory->destroy_2d_double_array(phiij);
+
+  memory->destroy_3d_double_array(Fij_spline);
+  memory->destroy_3d_double_array(Gij_spline);
+  memory->destroy_3d_double_array(phiij_spline);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEIM::compute(int eflag, int vflag)
+{
+  int i,j,ii,jj,m,inum,jnum,itype,jtype;
+  double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
+  double rsq,r,p,rhoip,rhojp,phip,phi,coul,coulp,recip,psip;
+  double *coeff;
+  int *ilist,*jlist,*numneigh,**firstneigh;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = eflag_global = eflag_atom = 0;
+
+  // grow energy array if necessary
+
+  if (atom->nmax > nmax) {
+    memory->sfree(rho);
+    memory->sfree(fp);
+    nmax = atom->nmax;
+    rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho");
+    fp = (double *) memory->smalloc(nmax*sizeof(double),"pair:fp");
+  }
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int nlocal = atom->nlocal;
+  int newton_pair = force->newton_pair;
+
+  inum = list->inum;
+  ilist = list->ilist;
+  numneigh = list->numneigh;
+  firstneigh = list->firstneigh;
+  
+  // zero out density
+
+  if (newton_pair) {
+    m = nlocal + atom->nghost;
+    for (i = 0; i < m; i++) {
+      rho[i] = 0.0;
+      fp[i] = 0.0;
+    }
+  } else {
+    for (i = 0; i < nlocal; i++) {
+      rho[i] = 0.0;
+      fp[i] = 0.0;
+    }
+  }
+
+  // rho = density at each atom
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      jtype = type[j];
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      if (rsq < cutforcesq[itype][jtype]) {
+	p = sqrt(rsq)*rdr + 1.0;
+	m = static_cast<int> (p);
+	m = MIN(m,nr-1);
+	p -= m;
+	p = MIN(p,1.0);
+	coeff = Fij_spline[type2Fij[itype][jtype]][m];
+	rho[i] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+	if (newton_pair || j < nlocal) {
+	  coeff = Fij_spline[type2Fij[jtype][itype]][m];
+	  rho[j] += ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+	}
+      }
+    }
+  }
+
+  // communicate and sum densities
+
+  rhofp = 1;
+  if (newton_pair) comm->reverse_comm_pair(this);
+  comm->forward_comm_pair(this);
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+ 
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      jtype = type[j];
+ 
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+ 
+      if (rsq < cutforcesq[itype][jtype]) {
+        p = sqrt(rsq)*rdr + 1.0;
+        m = static_cast<int> (p);
+        m = MIN(m,nr-1);
+        p -= m;
+        p = MIN(p,1.0);
+        coeff = Gij_spline[type2Gij[itype][jtype]][m];
+        fp[i] += rho[j]*(((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6]);
+        if (newton_pair || j < nlocal) {
+          fp[j] += rho[i]*(((coeff[3]*p + coeff[4])*p + coeff[5])*p + 
+			   coeff[6]);
+        }
+      }
+    }
+  }
+
+  // communicate and sum modified densities
+
+  rhofp = 2;
+  if (newton_pair) comm->reverse_comm_pair(this);
+  comm->forward_comm_pair(this);
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    itype = type[i];
+    if (eflag) {
+      phi = 0.5*rho[i]*fp[i];
+      if (eflag_global) eng_vdwl += phi;
+      if (eflag_atom) eatom[i] += phi;
+    }
+  }
+
+  // compute forces on each atom
+  // loop over neighbors of my atoms
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    xtmp = x[i][0];
+    ytmp = x[i][1];
+    ztmp = x[i][2];
+    itype = type[i];
+
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+
+    for (jj = 0; jj < jnum; jj++) {
+      j = jlist[jj];
+      jtype = type[j];
+      delx = xtmp - x[j][0];
+      dely = ytmp - x[j][1];
+      delz = ztmp - x[j][2];
+      rsq = delx*delx + dely*dely + delz*delz;
+
+      if (rsq < cutforcesq[itype][jtype]) {
+        r = sqrt(rsq);
+        p = r*rdr + 1.0;
+        m = static_cast<int> (p);
+        m = MIN(m,nr-1);
+        p -= m;
+        p = MIN(p,1.0);
+
+        // rhoip = derivative of (density at atom j due to atom i)
+        // rhojp = derivative of (density at atom i due to atom j)
+        // phi = pair potential energy
+        // phip = phi'
+
+        coeff = Fij_spline[type2Fij[jtype][itype]][m];
+        rhoip = (coeff[0]*p + coeff[1])*p + coeff[2];
+        coeff = Fij_spline[type2Fij[itype][jtype]][m];
+        rhojp = (coeff[0]*p + coeff[1])*p + coeff[2];
+        coeff = phiij_spline[type2phiij[itype][jtype]][m];
+        phip = (coeff[0]*p + coeff[1])*p + coeff[2];
+        phi = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+        coeff = Gij_spline[type2Gij[itype][jtype]][m];
+        coul = ((coeff[3]*p + coeff[4])*p + coeff[5])*p + coeff[6];
+        coulp = (coeff[0]*p + coeff[1])*p + coeff[2];
+        psip = phip + (rho[i]*rho[j]-q0[itype]*q0[jtype])*coulp +
+               fp[i]*rhojp + fp[j]*rhoip;
+        recip = 1.0/r;
+        fpair = -psip*recip;
+        f[i][0] += delx*fpair;
+        f[i][1] += dely*fpair;
+        f[i][2] += delz*fpair;
+        if (newton_pair || j < nlocal) {
+          f[j][0] -= delx*fpair;
+          f[j][1] -= dely*fpair;
+          f[j][2] -= delz*fpair;
+        }
+
+        if (eflag) evdwl = phi-q0[itype]*q0[jtype]*coul;
+        if (evflag) ev_tally(i,j,nlocal,newton_pair,
+                             evdwl,0.0,fpair,delx,dely,delz);
+      }
+    }
+  }
+
+  if (vflag_fdotr) virial_compute();
+}
+
+/* ----------------------------------------------------------------------
+   allocate all arrays
+------------------------------------------------------------------------- */
+
+void PairEIM::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+
+  map = new int[n+1];
+  for (int i = 1; i <= n; i++) map[i] = -1;
+
+  type2Fij = memory->create_2d_int_array(n+1,n+1,"pair:type2Fij");
+  type2Gij = memory->create_2d_int_array(n+1,n+1,"pair:type2Gij");
+  type2phiij = memory->create_2d_int_array(n+1,n+1,"pair:type2phiij");
+}
+
+/* ----------------------------------------------------------------------
+   global settings
+------------------------------------------------------------------------- */
+
+void PairEIM::settings(int narg, char **arg)
+{
+  if (narg > 0) error->all("Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs from set file
+------------------------------------------------------------------------- */
+
+void PairEIM::coeff(int narg, char **arg)
+{
+  int i,j,m,n;
+
+  if (!allocated) allocate();
+
+  if (narg < 5) error->all("Incorrect args for pair coefficients");
+
+  // insure I,J args are * *
+
+  if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
+    error->all("Incorrect args for pair coefficients");
+
+  // read EIM element names before filename
+  // nelements = # of EIM elements to read from file
+  // elements = list of unique element names
+
+  if (nelements) {
+    for (i = 0; i < nelements; i++) delete [] elements[i];
+    delete [] elements;
+  }
+  nelements = narg - 3 - atom->ntypes;
+  if (nelements < 1) error->all("Incorrect args for pair coefficients");
+  elements = new char*[nelements];
+  
+  for (i = 0; i < nelements; i++) {
+    n = strlen(arg[i+2]) + 1;
+    elements[i] = new char[n];
+    strcpy(elements[i],arg[i+2]);
+  }
+
+  // read EIM file
+
+  deallocate_setfl();
+  setfl = new Setfl();
+  read_file(arg[2+nelements]);
+
+  // read args that map atom types to elements in potential file
+  // map[i] = which element the Ith atom type is, -1 if NULL
+
+  for (i = 3 + nelements; i < narg; i++) {
+    m = i - (3+nelements) + 1;
+    for (j = 0; j < nelements; j++)
+      if (strcmp(arg[i],elements[j]) == 0) break;
+    if (j < nelements) map[m] = j;
+    else if (strcmp(arg[i],"NULL") == 0) map[m] = -1;
+    else error->all("Incorrect args for pair coefficients");
+  }
+
+  // clear setflag since coeff() called once with I,J = * *
+
+  n = atom->ntypes;
+  for (i = 1; i <= n; i++)
+    for (j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  // set setflag i,j for type pairs where both are mapped to elements
+  // set mass of atom type if i = j
+
+  int count = 0;
+  for (i = 1; i <= n; i++)
+    for (j = i; j <= n; j++)
+      if (map[i] >= 0 && map[j] >= 0) {
+	setflag[i][j] = 1;
+	if (i == j) atom->set_mass(i,setfl->mass[map[i]]);
+	count++;
+      }
+
+  if (count == 0) error->all("Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairEIM::init_style()
+{
+  // convert read-in file(s) to arrays and spline them
+
+  file2array();
+  array2spline();
+
+  int irequest = neighbor->request(this);
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairEIM::init_one(int i, int j)
+{
+  cutmax = sqrt(cutforcesq[i][j]);
+  return cutmax;
+}
+
+/* ----------------------------------------------------------------------
+   read potential values from a set file
+------------------------------------------------------------------------- */
+
+void PairEIM::read_file(char *filename)
+{
+  // open potential file
+
+  int me = comm->me;
+  FILE *fptr;
+
+  if (me == 0) {
+    fptr = fopen(filename,"r");
+    if (fptr == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open EIM potential file %s",filename);
+      error->one(str);
+    }
+  }
+
+  int npair = nelements*(nelements+1)/2;
+  setfl->ielement = new int[nelements];
+  setfl->mass = new double[nelements];
+  setfl->negativity = new double[nelements];
+  setfl->ra = new double[nelements];
+  setfl->ri = new double[nelements];
+  setfl->Ec = new double[nelements];
+  setfl->q0 = new double[nelements];
+  setfl->rcutphiA = new double[npair];
+  setfl->rcutphiR = new double[npair];
+  setfl->Eb = new double[npair];
+  setfl->r0 = new double[npair];
+  setfl->alpha = new double[npair];
+  setfl->beta = new double[npair];
+  setfl->rcutq = new double[npair];
+  setfl->Asigma = new double[npair];
+  setfl->rq = new double[npair];
+  setfl->rcutsigma = new double[npair];
+  setfl->Ac = new double[npair];
+  setfl->zeta = new double[npair];
+  setfl->rs = new double[npair];
+  setfl->tp = new int[npair];
+
+  if (me == 0)
+    if (!grabglobal(fptr))
+      error->one("Could not grab global entry from EIM potential file");
+  MPI_Bcast(&setfl->division,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&setfl->rbig,1,MPI_DOUBLE,0,world);
+  MPI_Bcast(&setfl->rsmall,1,MPI_DOUBLE,0,world);
+
+  for (int i = 0; i < nelements; i++) {
+    if (me == 0)
+      if (!grabsingle(fptr,i))
+	error->one("Could not grab element entry from EIM potential file");
+    MPI_Bcast(&setfl->ielement[i],1,MPI_INT,0,world);
+    MPI_Bcast(&setfl->mass[i],1,MPI_DOUBLE,0,world);
+    MPI_Bcast(&setfl->negativity[i],1,MPI_DOUBLE,0,world);
+    MPI_Bcast(&setfl->ra[i],1,MPI_DOUBLE,0,world);
+    MPI_Bcast(&setfl->ri[i],1,MPI_DOUBLE,0,world);
+    MPI_Bcast(&setfl->Ec[i],1,MPI_DOUBLE,0,world);
+    MPI_Bcast(&setfl->q0[i],1,MPI_DOUBLE,0,world);
+  }
+
+  for (int i = 0; i < nelements; i++) {
+    for (int j = i; j < nelements; j++) {
+      int ij;
+      if (i == j) ij = i;
+      else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+      else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
+      if (me == 0)
+        if (grabpair(fptr,i,j) == 0)
+	  error->one("Could not grab pair entry from EIM potential file");
+      MPI_Bcast(&setfl->rcutphiA[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->rcutphiR[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->Eb[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->r0[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->alpha[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->beta[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->rcutq[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->Asigma[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->rq[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->rcutsigma[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->Ac[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->zeta[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->rs[ij],1,MPI_DOUBLE,0,world);
+      MPI_Bcast(&setfl->tp[ij],1,MPI_INT,0,world);
+    }
+  }
+
+  setfl->nr = 5000;
+  setfl->cut = 0.0;
+  for (int i = 0; i < npair; i++) {
+    if (setfl->cut < setfl->rcutphiA[i]) setfl->cut = setfl->rcutphiA[i];
+    if (setfl->cut < setfl->rcutphiR[i]) setfl->cut = setfl->rcutphiR[i];
+    if (setfl->cut < setfl->rcutq[i]) setfl->cut = setfl->rcutq[i];
+    if (setfl->cut < setfl->rcutsigma[i]) setfl->cut = setfl->rcutsigma[i];
+  }
+  setfl->dr = setfl->cut/(setfl->nr-1.0);
+
+  setfl->cuts = memory->create_2d_double_array(nelements,
+					      nelements,"pair:cuts");
+  for (int i = 0; i < nelements; i++) {
+    for (int j = 0; j < nelements; j++) {
+      if (i > j) {
+        setfl->cuts[i][j] = setfl->cuts[j][i];
+      } else {
+        int ij;
+        if (i == j) {
+          ij = i;
+        } else {
+          ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+        }
+        setfl->cuts[i][j] = setfl->rcutphiA[ij];
+        if (setfl->cuts[i][j] < setfl->rcutphiR[ij]) 
+	  setfl->cuts[i][j] = setfl->rcutphiR[ij];
+        if (setfl->cuts[i][j] < setfl->rcutq[ij]) 
+	  setfl->cuts[i][j] = setfl->rcutq[ij];
+        if (setfl->cuts[i][j] < setfl->rcutsigma[ij]) 
+	  setfl->cuts[i][j] = setfl->rcutsigma[ij];
+      }
+    }
+  }
+  
+  setfl->Fij = memory->create_3d_double_array(nelements,nelements,
+                                             setfl->nr+1,"pair:Fij");
+  setfl->Gij = memory->create_3d_double_array(nelements,nelements,
+                                             setfl->nr+1,"pair:Gij");
+  setfl->phiij = memory->create_3d_double_array(nelements,nelements,
+					     setfl->nr+1,"pair:phiij");
+  for (int i = 0; i < nelements; i++)  
+    for (int j = 0; j < nelements; j++) {
+      for (int k = 0; k < setfl->nr; k++) {
+        if (i > j) {
+          setfl->phiij[i][j][k+1] = setfl->phiij[j][i][k+1];
+        } else {
+          double r = k*setfl->dr;
+          setfl->phiij[i][j][k+1] = funcphi(i,j,r);
+        }
+      }
+    }
+
+  for (int i = 0; i < nelements; i++)
+    for (int j = 0; j < nelements; j++) {
+      for (int k = 0; k < setfl->nr; k++) {
+        double r = k*setfl->dr;
+        setfl->Fij[i][j][k+1] = funcsigma(i,j,r);
+      }
+    }
+
+  for (int i = 0; i < nelements; i++)
+    for (int j = 0; j < nelements; j++) {
+      for (int k = 0; k < setfl->nr; k++) {
+        if (i > j) {
+          setfl->Gij[i][j][k+1] = setfl->Gij[j][i][k+1];
+        } else {
+          double r = k*setfl->dr;
+          setfl->Gij[i][j][k+1] = funccoul(i,j,r);
+        }
+      }
+    }
+
+  // close the potential file
+
+  if (me == 0) fclose(fptr);
+}
+
+/* ----------------------------------------------------------------------
+   deallocate data associated with setfl file
+------------------------------------------------------------------------- */
+
+void PairEIM::deallocate_setfl()
+{
+  if (!setfl) return;
+  delete [] setfl->ielement;
+  delete [] setfl->mass;
+  delete [] setfl->negativity;
+  delete [] setfl->ra;
+  delete [] setfl->ri;
+  delete [] setfl->Ec;
+  delete [] setfl->q0;
+  delete [] setfl->rcutphiA;
+  delete [] setfl->rcutphiR;
+  delete [] setfl->Eb;
+  delete [] setfl->r0;
+  delete [] setfl->alpha;
+  delete [] setfl->beta;
+  delete [] setfl->rcutq;
+  delete [] setfl->Asigma;
+  delete [] setfl->rq;
+  delete [] setfl->rcutsigma;
+  delete [] setfl->Ac;
+  delete [] setfl->zeta;
+  delete [] setfl->rs;
+  delete [] setfl->tp;
+  memory->destroy_2d_double_array(setfl->cuts);
+  memory->destroy_3d_double_array(setfl->Fij);
+  memory->destroy_3d_double_array(setfl->Gij);
+  memory->destroy_3d_double_array(setfl->phiij);
+  delete setfl;
+}
+
+/* ----------------------------------------------------------------------
+   convert read-in potentials to standard array format
+   interpolate all file values to a single grid and cutoff
+------------------------------------------------------------------------- */
+
+void PairEIM::file2array()
+{
+  int i,j,m,n;
+  int irow,icol;
+  int ntypes = atom->ntypes;
+
+  delete [] negativity;
+  delete [] q0;
+  delete [] cutforcesq;
+  negativity = new double[ntypes+1];
+  q0 = new double[ntypes+1];
+  cutforcesq = memory->create_2d_double_array(ntypes+1,ntypes+1,
+					      "pair:cutforcesq");
+  for (i = 1; i <= ntypes; i++) {
+    if (map[i] == -1) {
+      negativity[i]=0.0;
+      q0[i]=0.0;
+    } else {
+      negativity[i]=setfl->negativity[map[i]];
+      q0[i]=setfl->q0[map[i]];
+    }
+  }
+
+  for (i = 1; i <= ntypes; i++)
+    for (j = 1; j <= ntypes; j++) {
+      if (map[i] == -1 || map[j] == -1) {
+        cutforcesq[i][j] = setfl->cut;
+        cutforcesq[i][j] =  cutforcesq[i][j]*cutforcesq[i][j];
+      } else {
+        cutforcesq[i][j] = setfl->cuts[map[i]][map[j]];
+        cutforcesq[i][j] =  cutforcesq[i][j]*cutforcesq[i][j];
+      }
+    }
+
+  nr = setfl->nr;
+  dr = setfl->dr;
+
+  // ------------------------------------------------------------------
+  // setup Fij arrays
+  // ------------------------------------------------------------------
+
+  nFij = nelements*nelements + 1;
+  memory->destroy_2d_double_array(Fij);
+  Fij = (double **) memory->create_2d_double_array(nFij,nr+1,"pair:Fij");
+
+  // copy each element's Fij to global Fij
+
+  n=0;
+  for (i = 0; i < nelements; i++)
+    for (j = 0; j < nelements; j++) {
+      for (m = 1; m <= nr; m++) Fij[n][m] = setfl->Fij[i][j][m];
+      n++;
+    }
+
+  // add extra Fij of zeroes for non-EIM types to point to (pair hybrid)
+ 
+  for (m = 1; m <= nr; m++) Fij[nFij-1][m] = 0.0;
+ 
+  // type2Fij[i][j] = which Fij array (0 to nFij-1) each type pair maps to
+  // setfl of Fij arrays 
+  // value = n = sum over rows of matrix until reach irow,icol
+  // if atom type doesn't point to element (non-EIM atom in pair hybrid)
+  // then map it to last Fij array of zeroes
+
+  for (i = 1; i <= ntypes; i++) {
+    for (j = 1; j <= ntypes; j++) {
+      irow = map[i];
+      icol = map[j];
+      if (irow == -1 || icol == -1) {
+        type2Fij[i][j] = nFij-1;
+      } else {
+        n = 0;
+        for (m = 0; m < irow; m++) n += nelements;
+        n += icol;
+        type2Fij[i][j] = n;
+      }
+    }
+  }
+
+  // ------------------------------------------------------------------
+  // setup Gij arrays
+  // ------------------------------------------------------------------
+
+  nGij = nelements * (nelements+1) / 2 + 1;
+  memory->destroy_2d_double_array(Gij);
+  Gij = (double **) memory->create_2d_double_array(nGij,nr+1,"pair:Gij");
+
+  // copy each element's Gij to global Gij, only for I >= J
+
+  n=0;
+  for (i = 0; i < nelements; i++)
+    for (j = 0; j <= i; j++) {
+      for (m = 1; m <= nr; m++) Gij[n][m] = setfl->Gij[i][j][m];
+      n++;
+    }
+
+  // add extra Gij of zeroes for non-EIM types to point to (pair hybrid)
+
+  for (m = 1; m <= nr; m++) Gij[nGij-1][m] = 0.0;
+
+  // type2Gij[i][j] = which Gij array (0 to nGij-1) each type pair maps to
+  // setfl of Gij arrays only fill lower triangular Nelement matrix
+  // value = n = sum over rows of lower-triangular matrix until reach irow,icol
+  // swap indices when irow < icol to stay lower triangular
+  // if atom type doesn't point to element (non-EIM atom in pair hybrid)
+  // then map it to last Gij array of zeroes
+
+  for (i = 1; i <= ntypes; i++) {
+    for (j = 1; j <= ntypes; j++) {
+      irow = map[i];
+      icol = map[j];
+      if (irow == -1 || icol == -1) {
+        type2Gij[i][j] = nGij-1;
+      } else {
+        if (irow < icol) {
+          irow = map[j];
+          icol = map[i];
+        }
+        n = 0;
+        for (m = 0; m < irow; m++) n += m + 1;
+        n += icol;
+        type2Gij[i][j] = n;
+      }
+    }
+  }
+
+  // ------------------------------------------------------------------
+  // setup phiij arrays
+  // ------------------------------------------------------------------
+
+  nphiij = nelements * (nelements+1) / 2 + 1;
+  memory->destroy_2d_double_array(phiij);
+  phiij = (double **) memory->create_2d_double_array(nphiij,nr+1,"pair:phiij");
+
+  // copy each element pair phiij to global phiij, only for I >= J
+
+  n = 0;
+  for (i = 0; i < nelements; i++)
+    for (j = 0; j <= i; j++) {
+      for (m = 1; m <= nr; m++) phiij[n][m] = setfl->phiij[i][j][m];
+      n++;
+    }
+
+  // add extra phiij of zeroes for non-EIM types to point to (pair hybrid)
+
+  for (m = 1; m <= nr; m++) phiij[nphiij-1][m] = 0.0;
+
+  // type2phiij[i][j] = which phiij array (0 to nphiij-1) 
+  //                    each type pair maps to
+  // setfl of phiij arrays only fill lower triangular Nelement matrix
+  // value = n = sum over rows of lower-triangular matrix until reach irow,icol
+  // swap indices when irow < icol to stay lower triangular
+  // if atom type doesn't point to element (non-EIM atom in pair hybrid)
+  // then map it to last phiij array of zeroes
+
+  for (i = 1; i <= ntypes; i++) {
+    for (j = 1; j <= ntypes; j++) {
+      irow = map[i];
+      icol = map[j];
+      if (irow == -1 || icol == -1) {
+        type2phiij[i][j] = nphiij-1;
+      } else {
+        if (irow < icol) {
+          irow = map[j];
+          icol = map[i];
+        }
+        n = 0;
+        for (m = 0; m < irow; m++) n += m + 1;
+        n += icol;
+        type2phiij[i][j] = n;
+      }
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEIM::array2spline()
+{
+  rdr = 1.0/dr;
+
+  memory->destroy_3d_double_array(Fij_spline);
+  memory->destroy_3d_double_array(Gij_spline);
+  memory->destroy_3d_double_array(phiij_spline);
+
+  Fij_spline = memory->create_3d_double_array(nFij,nr+1,7,"pair:Fij");
+  Gij_spline = memory->create_3d_double_array(nGij,nr+1,7,"pair:Gij");
+  phiij_spline = memory->create_3d_double_array(nphiij,nr+1,7,"pair:phiij");
+
+  for (int i = 0; i < nFij; i++)
+    interpolate(nr,dr,Fij[i],Fij_spline[i],0.0);
+
+  for (int i = 0; i < nGij; i++)
+    interpolate(nr,dr,Gij[i],Gij_spline[i],0.0);
+
+  for (int i = 0; i < nphiij; i++)
+    interpolate(nr,dr,phiij[i],phiij_spline[i],0.0);
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEIM::interpolate(int n, double delta, double *f, 
+			  double **spline, double origin)
+{
+  for (int m = 1; m <= n; m++) spline[m][6] = f[m];
+
+  spline[1][5] = spline[2][6] - spline[1][6];
+  spline[2][5] = 0.5 * (spline[3][6]-spline[1][6]);
+  spline[n-1][5] = 0.5 * (spline[n][6]-spline[n-2][6]);
+  spline[n][5] = 0.0;
+  
+  for (int m = 3; m <= n-2; m++)
+    spline[m][5] = ((spline[m-2][6]-spline[m+2][6]) + 
+		    8.0*(spline[m+1][6]-spline[m-1][6])) / 12.0;
+  
+  for (int m = 1; m <= n-1; m++) {
+    spline[m][4] = 3.0*(spline[m+1][6]-spline[m][6]) - 
+      2.0*spline[m][5] - spline[m+1][5];
+    spline[m][3] = spline[m][5] + spline[m+1][5] - 
+      2.0*(spline[m+1][6]-spline[m][6]);
+  }
+  
+  spline[n][4] = 0.0;
+  spline[n][3] = 0.0;
+  
+  for (int m = 1; m <= n; m++) {
+    spline[m][2] = spline[m][5]/delta;
+    spline[m][1] = 2.0*spline[m][4]/delta;
+    spline[m][0] = 3.0*spline[m][3]/delta;
+  }
+}
+
+/* ----------------------------------------------------------------------
+   grab global line from file and store info in setfl
+   return 0 if error
+------------------------------------------------------------------------- */
+
+int PairEIM::grabglobal(FILE *fptr)
+{
+  char line[MAXLINE];
+
+  char *pch = NULL, *data = NULL;
+  while (pch == NULL) {
+    if (fgets(line,MAXLINE,fptr) == NULL) break;
+    pch = strstr(line,"global");
+    if (pch != NULL) {
+      data = strtok (line," \t\n\r\f");
+      data = strtok (NULL,"?");
+      sscanf(data,"%lg %lg %lg",&setfl->division,&setfl->rbig,&setfl->rsmall);
+    }
+  }
+  if (pch == NULL) return 0;
+  return 1;
+}
+
+/* ----------------------------------------------------------------------
+   grab elemental line from file and store info in setfl
+   return 0 if error
+------------------------------------------------------------------------- */
+
+int PairEIM::grabsingle(FILE *fptr, int i)
+{
+  char line[MAXLINE];
+
+  rewind(fptr);
+
+  char *pch1 = NULL, *pch2 = NULL, *data = NULL;
+  while (pch1 == NULL || pch2 == NULL) {
+    if (fgets(line,MAXLINE,fptr) == NULL) break;
+    pch1 = strtok (line," \t\n\r\f");
+    pch1 = strstr(pch1,"element:");
+    if (pch1 != NULL) {
+      pch2 = strtok(NULL, " \t\n\r\f");
+      if (pch2 != NULL) data = strtok (NULL, "?");
+      if (strcmp(pch2,elements[i]) == 0) {
+        sscanf(data,"%d %lg %lg %lg %lg %lg %lg",&setfl->ielement[i],
+          &setfl->mass[i],&setfl->negativity[i],&setfl->ra[i],
+          &setfl->ri[i],&setfl->Ec[i],&setfl->q0[i]);
+      } else {
+        pch2 = NULL;
+      }
+    }
+  }
+  if (pch1 == NULL || pch2 == NULL) return 0;
+  return 1;
+}
+
+/* ----------------------------------------------------------------------
+   grab pair line from file and store info in setfl
+   return 0 if error
+------------------------------------------------------------------------- */
+
+int PairEIM::grabpair(FILE *fptr, int i, int j)
+{
+  char line[MAXLINE];
+
+  rewind(fptr);
+
+  int ij;
+  if (i == j) ij = i;
+  else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+  else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
+
+  char *pch1 = NULL, *pch2 = NULL, *pch3 = NULL, *data = NULL;
+  while (pch1 == NULL || pch2 == NULL || pch3 == NULL) {
+    if (fgets(line,MAXLINE,fptr) == NULL) break;
+    pch1 = strtok (line," \t\n\r\f");
+    pch1 = strstr(pch1,"pair:");
+    if (pch1 != NULL) {
+      pch2 = strtok (NULL, " \t\n\r\f");
+      if (pch2 != NULL) pch3 = strtok (NULL, " \t\n\r\f");
+      if (pch3 != NULL) data = strtok (NULL, "?");
+      if ((strcmp(pch2,elements[i]) == 0 &&
+        strcmp(pch3,elements[j]) == 0) ||
+        (strcmp(pch2,elements[j]) == 0 &&
+        strcmp(pch3,elements[i]) == 0)) {
+        sscanf(data,"%lg %lg %lg %lg %lg",
+          &setfl->rcutphiA[ij],&setfl->rcutphiR[ij],
+          &setfl->Eb[ij],&setfl->r0[ij],&setfl->alpha[ij]);
+        fgets(line,MAXLINE,fptr);
+        sscanf(line,"%lg %lg %lg %lg %lg",
+          &setfl->beta[ij],&setfl->rcutq[ij],&setfl->Asigma[ij],
+          &setfl->rq[ij],&setfl->rcutsigma[ij]);
+        fgets(line,MAXLINE,fptr);
+        sscanf(line,"%lg %lg %lg %d",
+          &setfl->Ac[ij],&setfl->zeta[ij],&setfl->rs[ij],
+          &setfl->tp[ij]);
+      } else {
+         pch1 = NULL;
+         pch2 = NULL;
+         pch3 = NULL;
+      }
+    }
+  }
+  if (pch1 == NULL || pch2 == NULL || pch3 == NULL) return 0;
+  return 1;
+}
+
+/* ----------------------------------------------------------------------
+   cutoff function
+------------------------------------------------------------------------- */
+
+double PairEIM::funccutoff(double rp, double rc, double r)
+{
+  double rbig = setfl->rbig;
+  double rsmall = setfl->rsmall;
+  
+  double a = (rsmall-rbig)/(rc-rp)*(r-rp)+rbig;
+  a = erfc(a);
+  double b = erfc(rbig);
+  double c = erfc(rsmall);
+  return (a-c)/(b-c);
+}
+
+/* ----------------------------------------------------------------------
+   pair interaction function phi
+------------------------------------------------------------------------- */
+
+double PairEIM::funcphi(int i, int j, double r)
+{
+  int ij;
+  double value = 0.0;
+  if (i == j) ij = i;
+  else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+  else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
+  if (r < 0.2) r = 0.2;
+  if (setfl->tp[ij] == 1) {
+    double a = setfl->Eb[ij]*setfl->alpha[ij] / 
+      (setfl->beta[ij]-setfl->alpha[ij]);
+    double b = setfl->Eb[ij]*setfl->beta[ij] / 
+      (setfl->beta[ij]-setfl->alpha[ij]);
+    if (r < setfl->rcutphiA[ij]) {
+      value -= a*exp(-setfl->beta[ij]*(r/setfl->r0[ij]-1.0))*
+        funccutoff(setfl->r0[ij],setfl->rcutphiA[ij],r);
+    }
+    if (r < setfl-> rcutphiR[ij]) {
+      value += b*exp(-setfl->alpha[ij]*(r/setfl->r0[ij]-1.0))*
+        funccutoff(setfl->r0[ij],setfl->rcutphiR[ij],r);
+    }
+  } else if (setfl->tp[ij] == 2) {
+    double a=setfl->Eb[ij]*setfl->alpha[ij]*pow(setfl->r0[ij],setfl->beta[ij])/
+      (setfl->beta[ij]-setfl->alpha[ij]);
+    double b=a*setfl->beta[ij]/setfl->alpha[ij]*
+      pow(setfl->r0[ij],setfl->alpha[ij]-setfl->beta[ij]);
+    if (r < setfl->rcutphiA[ij]) {
+      value -= a/pow(r,setfl->beta[ij])*
+        funccutoff(setfl->r0[ij],setfl->rcutphiA[ij],r);
+    }
+    if (r < setfl-> rcutphiR[ij]) {
+      value += b/pow(r,setfl->alpha[ij])*
+        funccutoff(setfl->r0[ij],setfl->rcutphiR[ij],r);
+    }
+  }
+  return value;
+}
+
+/* ----------------------------------------------------------------------
+   ion propensity function sigma
+------------------------------------------------------------------------- */
+
+double PairEIM::funcsigma(int i, int j, double r)
+{
+  int ij;
+  double value = 0.0;
+  if (i == j) ij = i;
+  else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+  else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
+  if (r < 0.2) r = 0.2;
+  if (r < setfl->rcutq[ij]) {
+    value = setfl->Asigma[ij]*(setfl->negativity[j]-setfl->negativity[i]) *
+      funccutoff(setfl->rq[ij],setfl->rcutq[ij],r);
+  }
+  return value;
+}
+
+/* ----------------------------------------------------------------------
+   charge-charge interaction function sigma
+------------------------------------------------------------------------- */
+
+double PairEIM::funccoul(int i, int j, double r)
+{
+  int ij;
+  double value = 0.0;
+  if (i == j) ij = i;
+  else if (i < j) ij = nelements*(i+1) - (i+1)*(i+2)/2 + j;
+  else ij = nelements*(j+1) - (j+1)*(j+2)/2 + i;
+  if (r < 0.2) r = 0.2;
+  if (r < setfl->rcutsigma[ij]) {
+    value = setfl->Ac[ij]*exp(-setfl->zeta[ij]*r)*
+      funccutoff(setfl->rs[ij],setfl->rcutsigma[ij],r);
+  }
+  return value;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairEIM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+{
+  int i,j,m;
+ 
+  m = 0;
+  if (rhofp == 1) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = rho[j];
+    }
+  }
+  if (rhofp == 2) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      buf[m++] = fp[j];
+    }
+  }
+  return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEIM::unpack_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+ 
+  m = 0;
+  last = first + n;
+  if (rhofp == 1) {
+    for (i = first; i < last; i++) rho[i] = buf[m++];
+  }
+  if (rhofp == 2) {
+    for (i = first; i < last; i++) fp[i] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairEIM::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,m,last;
+
+  m = 0;
+  last = first + n;
+  if (rhofp == 1) {
+    for (i = first; i < last; i++) buf[m++] = rho[i];
+  }
+  if (rhofp == 2) {
+    for (i = first; i < last; i++) buf[m++] = fp[i];
+  }
+  return 1;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairEIM::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,m;
+
+  m = 0;
+  if (rhofp == 1) {
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      rho[j] += buf[m++];
+    }
+  }
+  if (rhofp == 2) { 
+    for (i = 0; i < n; i++) {
+      j = list[i];
+      fp[j] += buf[m++];
+    }
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays 
+------------------------------------------------------------------------- */
+
+double PairEIM::memory_usage()
+{
+  double bytes = maxeatom * sizeof(double);
+  bytes += maxvatom*6 * sizeof(double);
+  bytes += 2 * nmax * sizeof(double);
+  return bytes;
+}
diff --git a/src/MANYBODY/pair_eim.h b/src/MANYBODY/pair_eim.h
new file mode 100644
index 0000000000..22b89f2b98
--- /dev/null
+++ b/src/MANYBODY/pair_eim.h
@@ -0,0 +1,101 @@
+/* ----------------------------------------------------------------------
+   LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
+   http://lammps.sandia.gov, Sandia National Laboratories
+   Steve Plimpton, sjplimp@sandia.gov
+
+   Copyright (2003) Sandia Corporation.  Under the terms of Contract
+   DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
+   certain rights in this software.  This software is distributed under 
+   the GNU General Public License.
+
+   See the README file in the top-level LAMMPS directory.
+------------------------------------------------------------------------- */
+
+#ifdef PAIR_CLASS
+
+PairStyle(eim,PairEIM)
+
+#else
+
+#ifndef LMP_PAIR_EIM_H
+#define LMP_PAIR_EIM_H
+
+#include "stdio.h"
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairEIM : public Pair {
+ public:
+  PairEIM(class LAMMPS *);
+  ~PairEIM();
+  void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  double init_one(int, int);
+
+  int pack_comm(int, int *, double *, int, int *);
+  void unpack_comm(int, int, double *);
+  int pack_reverse_comm(int, int, double *);
+  void unpack_reverse_comm(int, int *, double *);
+  double memory_usage();
+
+ private:
+  double **cutforcesq,cutmax;
+  int nmax;
+  double *rho,*fp;
+  int rhofp;
+  int *map;                   // which element each atom type maps to
+
+  int nelements;              // # of elements to read from potential file
+  char **elements;            // element names
+
+  struct Setfl {
+    double division,rbig,rsmall;
+    int nr;
+    int *ielement,*tp;
+    double *mass,*negativity,*ra,*ri,*Ec,*q0;
+    double *rcutphiA,*rcutphiR,*Eb,*r0,*alpha,*beta,
+           *rcutq,*Asigma,*rq,*rcutsigma,*Ac,*zeta,
+           *rs;
+    double dr,cut;
+    double ***Fij,***Gij,***phiij;
+    double **cuts;
+  };
+  Setfl *setfl;
+
+  // potentials as array data
+
+  int nr;
+  int nFij,nGij,nphiij;
+  double **Fij,**Gij,**phiij;
+  int **type2Fij,**type2Gij,**type2phiij;
+  
+  // potentials in spline form used for force computation
+
+  double dr,rdr;
+  double *negativity,*q0;
+  double ***Fij_spline,***Gij_spline,***phiij_spline;
+
+  void allocate();
+  void array2spline();
+  void interpolate(int, double, double *, double **, double);
+  int grabglobal(FILE *); 
+  int grabsingle(FILE *, int);
+  int grabpair(FILE *, int, int);
+
+  double funccutoff(double, double, double);
+  double funcphi(int, int, double);
+  double funcsigma(int, int, double);
+  double funccoul(int, int, double);
+
+  void read_file(char *);
+  void deallocate_setfl();
+  void file2array();
+};
+
+}
+
+#endif
+#endif