diff --git a/src/MANYBODY/pair_eam.cpp b/src/MANYBODY/pair_eam.cpp
index c28268173d..de2c10501e 100644
--- a/src/MANYBODY/pair_eam.cpp
+++ b/src/MANYBODY/pair_eam.cpp
@@ -225,7 +225,7 @@ void PairEAM::compute(int eflag, int vflag)
 
   // communicate derivative of embedding function
 
-  comm->comm_pair(this);
+  comm->forward_comm_pair(this);
 
   // compute forces on each atom
   // loop over neighbors of my atoms
diff --git a/src/MEAM/pair_meam.cpp b/src/MEAM/pair_meam.cpp
index b5587a0960..0d8cc586cd 100644
--- a/src/MEAM/pair_meam.cpp
+++ b/src/MEAM/pair_meam.cpp
@@ -1,919 +1,919 @@
-/* ----------------------------------------------------------------------
-   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: Greg Wagner (SNL)
-------------------------------------------------------------------------- */
-
-#include "math.h"
-#include "stdio.h"
-#include "stdlib.h"
-#include "string.h"
-#include "pair_meam.h"
-#include "atom.h"
-#include "force.h"
-#include "comm.h"
-#include "memory.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.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
-
-enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11,L12};
-int nkeywords = 16;
-char *keywords[] = {"Ec","alpha","rho0","delta","lattce",
-		    "attrac","repuls","nn2","Cmin","Cmax","rc","delr",
-		    "augt1","gsmooth_factor","re","ialloy"};
-
-/* ---------------------------------------------------------------------- */
-
-PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp)
-{
-  single_enable = 0;
-  one_coeff = 1;
-
-  nmax = 0;
-  rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
-  gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
-  arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
-
-  maxneigh = 0;
-  scrfcn = dscrfcn = fcpair = NULL;
-
-  nelements = 0;
-  elements = NULL;
-  mass = NULL;
-
-  // set comm size needed by this Pair
-
-  comm_forward = 35;
-  comm_reverse = 27;
-}
-
-/* ----------------------------------------------------------------------
-   free all arrays
-   check if allocated, since class can be destructed when incomplete
-------------------------------------------------------------------------- */
-
-PairMEAM::~PairMEAM()
-{
-  meam_cleanup_();
-
-  memory->sfree(rho);
-  memory->sfree(rho0);
-  memory->sfree(rho1);
-  memory->sfree(rho2);
-  memory->sfree(rho3);
-  memory->sfree(frhop);
-  memory->sfree(gamma);
-  memory->sfree(dgamma1);
-  memory->sfree(dgamma2);
-  memory->sfree(dgamma3);
-  memory->sfree(arho2b);
-
-  memory->destroy_2d_double_array(arho1);
-  memory->destroy_2d_double_array(arho2);
-  memory->destroy_2d_double_array(arho3);
-  memory->destroy_2d_double_array(arho3b);
-  memory->destroy_2d_double_array(t_ave);
-  memory->destroy_2d_double_array(tsq_ave);
-
-  memory->sfree(scrfcn);
-  memory->sfree(dscrfcn);
-  memory->sfree(fcpair);
-  
-  for (int i = 0; i < nelements; i++) delete [] elements[i];
-  delete [] elements;
-  delete [] mass;
-
-  if (allocated) {
-    memory->destroy_2d_int_array(setflag);
-    memory->destroy_2d_double_array(cutsq);
-    delete [] map;
-    delete [] fmap;
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairMEAM::compute(int eflag, int vflag)
-{
-  int i,j,ii,n,inum_half,itype,jtype,errorflag;
-  double evdwl;
-  int *ilist_half,*jlist_half,*numneigh_half,**firstneigh_half;
-  int *numneigh_full,**firstneigh_full;
-
-  evdwl = 0.0;
-  if (eflag || vflag) ev_setup(eflag,vflag);
-  else evflag = vflag_fdotr = 0;
-
-  int newton_pair = force->newton_pair;
-
-  // grow local arrays if necessary
-
-  if (atom->nmax > nmax) {
-    memory->sfree(rho);
-    memory->sfree(rho0);
-    memory->sfree(rho1);
-    memory->sfree(rho2);
-    memory->sfree(rho3);
-    memory->sfree(frhop);
-    memory->sfree(gamma);
-    memory->sfree(dgamma1);
-    memory->sfree(dgamma2);
-    memory->sfree(dgamma3);
-    memory->sfree(arho2b);
-    memory->destroy_2d_double_array(arho1);
-    memory->destroy_2d_double_array(arho2);
-    memory->destroy_2d_double_array(arho3);
-    memory->destroy_2d_double_array(arho3b);
-    memory->destroy_2d_double_array(t_ave);
-    memory->destroy_2d_double_array(tsq_ave);
-
-    nmax = atom->nmax;
-
-    rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho");
-    rho0 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho0");
-    rho1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho1");
-    rho2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho2");
-    rho3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho3");
-    frhop = (double *) memory->smalloc(nmax*sizeof(double),"pair:frhop");
-    gamma = (double *) memory->smalloc(nmax*sizeof(double),"pair:gamma");
-    dgamma1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma1");
-    dgamma2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma2");
-    dgamma3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma3");
-    arho2b = (double *) memory->smalloc(nmax*sizeof(double),"pair:arho2b");
-    arho1 = memory->create_2d_double_array(nmax,3,"pair:arho1");
-    arho2 = memory->create_2d_double_array(nmax,6,"pair:arho2");
-    arho3 = memory->create_2d_double_array(nmax,10,"pair:arho3");
-    arho3b = memory->create_2d_double_array(nmax,3,"pair:arho3b");
-    t_ave = memory->create_2d_double_array(nmax,3,"pair:t_ave");
-    tsq_ave = memory->create_2d_double_array(nmax,3,"pair:tsq_ave");
-  }
-
-  // neighbor list info
-
-  inum_half = listhalf->inum;
-  ilist_half = listhalf->ilist;
-  numneigh_half = listhalf->numneigh;
-  firstneigh_half = listhalf->firstneigh;
-  numneigh_full = listfull->numneigh;
-  firstneigh_full = listfull->firstneigh;
-
-  // check size of scrfcn based on half neighbor list
-
-  int nlocal = atom->nlocal;
-  int nall = nlocal + atom->nghost;
-
-  n = 0;
-  for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
-
-  if (n > maxneigh) {
-    memory->sfree(scrfcn);
-    memory->sfree(dscrfcn);
-    memory->sfree(fcpair);
-    maxneigh = n;
-    scrfcn =
-      (double *) memory->smalloc(maxneigh*sizeof(double),"pair:scrfcn");
-    dscrfcn = 
-      (double *) memory->smalloc(maxneigh*sizeof(double),"pair:dscrfcn");
-    fcpair = 
-      (double *) memory->smalloc(maxneigh*sizeof(double),"pair:fcpair");
-  }
-
-  // zero out local arrays
-
-  for (i = 0; i < nall; i++) {
-    rho0[i] = 0.0;
-    arho2b[i] = 0.0;
-    arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
-    for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
-    for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
-    arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
-    t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
-    tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
-  }
-
-  double **x = atom->x;
-  double **f = atom->f;
-  int *type = atom->type;
-  int ntype = atom->ntypes;
-
-  // change neighbor list indices to Fortran indexing
-
-  neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half);
-  neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full);
-
-  // 3 stages of MEAM calculation
-  // loop over my atoms followed by communication
-
-  int ifort;
-  int offset = 0;
-  errorflag = 0;
-
-  for (ii = 0; ii < inum_half; ii++) {
-    i = ilist_half[ii];
-    ifort = i+1;
-    meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0],
-		    &numneigh_half[i],firstneigh_half[i],
-		    &numneigh_full[i],firstneigh_full[i],
-		    &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
-		    rho0,&arho1[0][0],&arho2[0][0],arho2b,
-		    &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],
-                    &errorflag);
-    if (errorflag) {
-      char str[128];
-      sprintf(str,"MEAM library error %d",errorflag);
-      error->one(str);
-    }
-    offset += numneigh_half[i];
-  }
-
-  comm->reverse_comm_pair(this);
-
-  meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
-		   &eng_vdwl,eatom,&ntype,type,fmap,
-		   &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
-		   &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
-		   dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
-  if (errorflag) {
-    char str[128];
-    sprintf(str,"MEAM library error %d",errorflag);
-    error->one(str);
-  }
-
-  comm->comm_pair(this);
-
-  offset = 0;
-
-  // vptr is first value in vatom if it will be used by meam_force()
-  // else vatom may not exist, so pass dummy ptr
-
-  double *vptr;
-  if (vflag_atom) vptr = &vatom[0][0];
-  else vptr = &cutmax;
-
-  for (ii = 0; ii < inum_half; ii++) {
-    i = ilist_half[ii];
-    ifort = i+1;
-    meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom,
-		&vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
-		&numneigh_half[i],firstneigh_half[i],
-		&numneigh_full[i],firstneigh_full[i],
-		&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
-		dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
-		&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
-		&t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag);
-    if (errorflag) {
-      char str[128];
-      sprintf(str,"MEAM library error %d",errorflag);
-      error->one(str);
-    }
-    offset += numneigh_half[i];
-  }
-
-  // change neighbor list indices back to C indexing
-
-  neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half);
-  neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full);
-
-  if (vflag_fdotr) virial_compute();
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairMEAM::allocate()
-{
-  allocated = 1;
-  int n = atom->ntypes;
-
-  setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
-  cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
-
-  map = new int[n+1];
-  fmap = new int[n];
-}
-
-/* ----------------------------------------------------------------------
-   global settings 
-------------------------------------------------------------------------- */
-
-void PairMEAM::settings(int narg, char **arg)
-{
-  if (narg != 0) error->all("Illegal pair_style command");
-}
-
-/* ----------------------------------------------------------------------
-   set coeffs for one or more type pairs
-------------------------------------------------------------------------- */
-
-void PairMEAM::coeff(int narg, char **arg)
-{
-  int i,j,m,n;
-
-  if (!allocated) allocate();
-
-  if (narg < 6) 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 MEAM element names between 2 filenames
-  // nelements = # of MEAM elements
-  // elements = list of unique element names
-
-  if (nelements) {
-    for (i = 0; i < nelements; i++) delete [] elements[i];
-    delete [] elements;
-    delete [] mass;
-  }
-  nelements = narg - 4 - atom->ntypes;
-  if (nelements < 1) error->all("Incorrect args for pair coefficients");
-  elements = new char*[nelements];
-  mass = new double[nelements];
-  
-  for (i = 0; i < nelements; i++) {
-    n = strlen(arg[i+3]) + 1;
-    elements[i] = new char[n];
-    strcpy(elements[i],arg[i+3]);
-  }
-
-  // read MEAM library and parameter files
-  // pass all parameters to MEAM package
-  // tell MEAM package that setup is done
-
-  read_files(arg[2],arg[2+nelements+1]);
-  meam_setup_done_(&cutmax);
-  
-  // read args that map atom types to MEAM elements
-  // map[i] = which element the Ith atom type is, -1 if not mapped
-
-  for (i = 4 + nelements; i < narg; i++) {
-    m = i - (4+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 (int i = 1; i <= n; i++)
-    for (int j = i; j <= n; j++)
-      setflag[i][j] = 0;
-
-  // set setflag i,j for type pairs where both are mapped to elements
-  // set mass for i,i in atom class
-
-  int count = 0;
-  for (int i = 1; i <= n; i++)
-    for (int j = i; j <= n; j++)
-      if (map[i] >= 0 && map[j] >= 0) {
-	setflag[i][j] = 1;
-	if (i == j) atom->set_mass(i,mass[map[i]]);
-	count++;
-      }
-
-  if (count == 0) error->all("Incorrect args for pair coefficients");
-}
-
-/* ----------------------------------------------------------------------
-   init specific to this pair style
-------------------------------------------------------------------------- */
-
-void PairMEAM::init_style()
-{
-  if (force->newton_pair == 0)
-    error->all("Pair style MEAM requires newton pair on");
-
-  // need full and half neighbor list
-
-  int irequest_full = neighbor->request(this);
-  neighbor->requests[irequest_full]->id = 1;
-  neighbor->requests[irequest_full]->half = 0;
-  neighbor->requests[irequest_full]->full = 1;
-  int irequest_half = neighbor->request(this);
-  neighbor->requests[irequest_half]->id = 2;
-  neighbor->requests[irequest_half]->half = 0;
-  neighbor->requests[irequest_half]->half_from_full = 1;
-  neighbor->requests[irequest_half]->otherlist = irequest_full;
-
-  // setup Fortran-style mapping array needed by MEAM package
-  // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index
-  // if type I is not a MEAM atom, fmap stores a 0
-
-  for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1;
-}
-
-/* ----------------------------------------------------------------------
-   neighbor callback to inform pair style of neighbor list to use
-   half or full
-------------------------------------------------------------------------- */
-
-void PairMEAM::init_list(int id, NeighList *ptr)
-{
-  if (id == 1) listfull = ptr;
-  else if (id == 2) listhalf = ptr;
-}
-
-/* ----------------------------------------------------------------------
-   init for one type pair i,j and corresponding j,i
-------------------------------------------------------------------------- */
-
-double PairMEAM::init_one(int i, int j)
-{
-  return cutmax;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairMEAM::read_files(char *globalfile, char *userfile)
-{
-  // open global meamf file on proc 0
-
-  FILE *fp;
-  if (comm->me == 0) {
-    fp = fopen(globalfile,"r");
-    if (fp == NULL) {
-      char str[128];
-      sprintf(str,"Cannot open MEAM potential file %s",globalfile);
-      error->one(str);
-    }
-  }
-
-  // allocate parameter arrays
-
-  int params_per_line = 19;
-
-  int *lat = new int[nelements];
-  int *ielement = new int[nelements];
-  int *ibar = new int[nelements];
-  double *z = new double[nelements];
-  double *atwt = new double[nelements];
-  double *alpha = new double[nelements];
-  double *b0 = new double[nelements];
-  double *b1 = new double[nelements];
-  double *b2 = new double[nelements];
-  double *b3 = new double[nelements];
-  double *alat = new double[nelements];
-  double *esub = new double[nelements];
-  double *asub = new double[nelements];
-  double *t0 = new double[nelements];
-  double *t1 = new double[nelements];
-  double *t2 = new double[nelements];
-  double *t3 = new double[nelements];
-  double *rozero = new double[nelements];
-
-  bool *found = new bool[nelements];
-  for (int i = 0; i < nelements; i++) found[i] = false;
-
-  // read each set of params from global MEAM file
-  // one set of params can span multiple lines
-  // store params if element name is in element list
-  // if element name appears multiple times, only store 1st entry
-
-  int i,n,nwords;
-  char **words = new char*[params_per_line+1];
-  char line[MAXLINE],*ptr;
-  int eof = 0;
-
-  int nset = 0;
-  while (1) {
-    if (comm->me == 0) {
-      ptr = fgets(line,MAXLINE,fp);
-      if (ptr == NULL) {
-	eof = 1;
-	fclose(fp);
-      } else n = strlen(line) + 1;
-    }
-    MPI_Bcast(&eof,1,MPI_INT,0,world);
-    if (eof) break;
-    MPI_Bcast(&n,1,MPI_INT,0,world);
-    MPI_Bcast(line,n,MPI_CHAR,0,world);
-
-    // strip comment, skip line if blank
-
-    if (ptr = strchr(line,'#')) *ptr = '\0';
-    nwords = atom->count_words(line);
-    if (nwords == 0) continue;
-
-    // concatenate additional lines until have params_per_line words
-
-    while (nwords < params_per_line) {
-      n = strlen(line);
-      if (comm->me == 0) {
-        ptr = fgets(&line[n],MAXLINE-n,fp);
-        if (ptr == NULL) {
-	  eof = 1;
-	  fclose(fp);
-        } else n = strlen(line) + 1;
-      }
-      MPI_Bcast(&eof,1,MPI_INT,0,world);
-      if (eof) break;
-      MPI_Bcast(&n,1,MPI_INT,0,world);
-      MPI_Bcast(line,n,MPI_CHAR,0,world);
-      if (ptr = strchr(line,'#')) *ptr = '\0';
-      nwords = atom->count_words(line);
-    }
-
-    if (nwords != params_per_line)
-      error->all("Incorrect format in MEAM potential file");
-
-    // words = ptrs to all words in line
-    // strip single and double quotes from words
-
-    nwords = 0;
-    words[nwords++] = strtok(line,"' \t\n\r\f");
-    while (words[nwords++] = strtok(NULL,"' \t\n\r\f")) continue;
-
-    // skip if element name isn't in element list
-
-    for (i = 0; i < nelements; i++)
-      if (strcmp(words[0],elements[i]) == 0) break;
-    if (i == nelements) continue;
-
-    // skip if element already appeared
-
-    if (found[i] == true) continue;
-    found[i] = true;
-
-    // map lat string to an integer
-
-    if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
-    else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
-    else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
-    else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
-    else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND;
-    else error->all("Unrecognized lattice type in MEAM file 1");
-
-    // store parameters
-
-    z[i] = atof(words[2]);
-    ielement[i] = atoi(words[3]);
-    atwt[i] = atof(words[4]);
-    alpha[i] = atof(words[5]);
-    b0[i] = atof(words[6]);
-    b1[i] = atof(words[7]);
-    b2[i] = atof(words[8]);
-    b3[i] = atof(words[9]);
-    alat[i] = atof(words[10]);
-    esub[i] = atof(words[11]);
-    asub[i] = atof(words[12]);
-    t0[i] = atof(words[13]);
-    t1[i] = atof(words[14]);
-    t2[i] = atof(words[15]);
-    t3[i] = atof(words[16]);
-    rozero[i] = atof(words[17]);
-    ibar[i] = atoi(words[18]);
-
-    nset++;
-  }
-
-  // error if didn't find all elements in file
-
-  if (nset != nelements)
-    error->all("Did not find all elements in MEAM library file");
-
-  // pass element parameters to MEAM package
-
-  meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
-  		     alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
-
-  // set element masses
-
-  for (i = 0; i < nelements; i++) mass[i] = atwt[i];
-
-  // clean-up memory
-
-  delete [] words;
-
-  delete [] lat;
-  delete [] ielement;
-  delete [] ibar;
-  delete [] z;
-  delete [] atwt;
-  delete [] alpha;
-  delete [] b0;
-  delete [] b1;
-  delete [] b2;
-  delete [] b3;
-  delete [] alat;
-  delete [] esub;
-  delete [] asub;
-  delete [] t0;
-  delete [] t1;
-  delete [] t2;
-  delete [] t3;
-  delete [] rozero;
-  delete [] found;
-
-  // done if user param file is NULL
-
-  if (strcmp(userfile,"NULL") == 0) return;
-
-  // open user param file on proc 0
-
-  if (comm->me == 0) {
-    fp = fopen(userfile,"r");
-    if (fp == NULL) {
-      char str[128];
-      sprintf(str,"Cannot open MEAM potential file %s",userfile);
-      error->one(str);
-    }
-  }
-
-  // read settings
-  // pass them one at a time to MEAM package
-  // match strings to list of corresponding ints
-
-  int which;
-  double value;
-  int nindex,index[3];
-  int maxparams = 6;
-  char **params = new char*[maxparams];
-  int nparams;
-
-  eof = 0;
-  while (1) {
-    if (comm->me == 0) {
-      ptr = fgets(line,MAXLINE,fp);
-      if (ptr == NULL) {
-	eof = 1;
-	fclose(fp);
-      } else n = strlen(line) + 1;
-    }
-    MPI_Bcast(&eof,1,MPI_INT,0,world);
-    if (eof) break;
-    MPI_Bcast(&n,1,MPI_INT,0,world);
-    MPI_Bcast(line,n,MPI_CHAR,0,world);
-
-    // strip comment, skip line if blank
-
-    if (ptr = strchr(line,'#')) *ptr = '\0';
-    nparams = atom->count_words(line);
-    if (nparams == 0) continue;
-
-    // words = ptrs to all words in line
-
-    nparams = 0;
-    params[nparams++] = strtok(line,"=(), '\t\n\r\f");
-    while (nparams < maxparams && 
-	   (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f")))
-      continue;
-    nparams--;
-
-    for (which = 0; which < nkeywords; which++)
-      if (strcmp(params[0],keywords[which]) == 0) break;
-    if (which == nkeywords) {
-      char str[128];
-      sprintf(str,"Keyword %s in MEAM parameter file not recognized",
-	      params[0]);
-      error->all(str);
-    }
-    nindex = nparams - 2;
-    for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]);
-
-    // map lattce_meam value to an integer
-
-    if (which == 4) {
-      if (strcmp(params[nparams-1],"fcc") == 0) value = FCC;
-      else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC;
-      else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP;
-      else if (strcmp(params[nparams-1],"dim") == 0) value = DIM;
-      else if (strcmp(params[nparams-1],"dia") == 0) value = DIAMOND;
-      else if (strcmp(params[nparams-1],"b1")  == 0) value = B1;
-      else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
-      else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
-      else error->all("Unrecognized lattice type in MEAM file 2");
-    }
-    else value = atof(params[nparams-1]);
-
-    // pass single setting to MEAM package
-
-    int errorflag = 0;
-    meam_setup_param_(&which,&value,&nindex,index,&errorflag);
-    if (errorflag) {
-      char str[128];
-      sprintf(str,"MEAM library error %d",errorflag);
-      error->all(str);
-    }
-  }
-
-  delete [] params;
-}
-
-/* ---------------------------------------------------------------------- */
-
-int PairMEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
-{
-  int i,j,k,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    buf[m++] = rho0[j];
-    buf[m++] = rho1[j];
-    buf[m++] = rho2[j];
-    buf[m++] = rho3[j];
-    buf[m++] = frhop[j];
-    buf[m++] = gamma[j];
-    buf[m++] = dgamma1[j];
-    buf[m++] = dgamma2[j];
-    buf[m++] = dgamma3[j];
-    buf[m++] = arho2b[j];
-    buf[m++] = arho1[j][0];
-    buf[m++] = arho1[j][1];
-    buf[m++] = arho1[j][2];
-    buf[m++] = arho2[j][0];
-    buf[m++] = arho2[j][1];
-    buf[m++] = arho2[j][2];
-    buf[m++] = arho2[j][3];
-    buf[m++] = arho2[j][4];
-    buf[m++] = arho2[j][5];
-    for (k = 0; k < 10; k++) buf[m++] = arho3[j][k];
-    buf[m++] = arho3b[j][0];
-    buf[m++] = arho3b[j][1];
-    buf[m++] = arho3b[j][2];
-    buf[m++] = t_ave[j][0];
-    buf[m++] = t_ave[j][1];
-    buf[m++] = t_ave[j][2];
-    buf[m++] = tsq_ave[j][0];
-    buf[m++] = tsq_ave[j][1];
-    buf[m++] = tsq_ave[j][2];
-  }
-  return 38;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairMEAM::unpack_comm(int n, int first, double *buf)
-{
-  int i,k,m,last;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    rho0[i] = buf[m++];
-    rho1[i] = buf[m++];
-    rho2[i] = buf[m++];
-    rho3[i] = buf[m++];
-    frhop[i] = buf[m++];
-    gamma[i] = buf[m++];
-    dgamma1[i] = buf[m++];
-    dgamma2[i] = buf[m++];
-    dgamma3[i] = buf[m++];
-    arho2b[i] = buf[m++];
-    arho1[i][0] = buf[m++];
-    arho1[i][1] = buf[m++];
-    arho1[i][2] = buf[m++];
-    arho2[i][0] = buf[m++];
-    arho2[i][1] = buf[m++];
-    arho2[i][2] = buf[m++];
-    arho2[i][3] = buf[m++];
-    arho2[i][4] = buf[m++];
-    arho2[i][5] = buf[m++];
-    for (k = 0; k < 10; k++) arho3[i][k] = buf[m++];
-    arho3b[i][0] = buf[m++];
-    arho3b[i][1] = buf[m++];
-    arho3b[i][2] = buf[m++];
-    t_ave[i][0] = buf[m++];
-    t_ave[i][1] = buf[m++];
-    t_ave[i][2] = buf[m++];
-    tsq_ave[i][0] = buf[m++];
-    tsq_ave[i][1] = buf[m++];
-    tsq_ave[i][2] = buf[m++];
-  }
-}
-
-/* ---------------------------------------------------------------------- */
-
-int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
-{
-  int i,k,m,last,size;
-
-  m = 0;
-  last = first + n;
-  for (i = first; i < last; i++) {
-    buf[m++] = rho0[i];
-    buf[m++] = arho2b[i];
-    buf[m++] = arho1[i][0];
-    buf[m++] = arho1[i][1];
-    buf[m++] = arho1[i][2];
-    buf[m++] = arho2[i][0];
-    buf[m++] = arho2[i][1];
-    buf[m++] = arho2[i][2];
-    buf[m++] = arho2[i][3];
-    buf[m++] = arho2[i][4];
-    buf[m++] = arho2[i][5];
-    for (k = 0; k < 10; k++) buf[m++] = arho3[i][k];
-    buf[m++] = arho3b[i][0];
-    buf[m++] = arho3b[i][1];
-    buf[m++] = arho3b[i][2];
-    buf[m++] = t_ave[i][0];
-    buf[m++] = t_ave[i][1];
-    buf[m++] = t_ave[i][2];
-    buf[m++] = tsq_ave[i][0];
-    buf[m++] = tsq_ave[i][1];
-    buf[m++] = tsq_ave[i][2];
-  }
-
-  return 30;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
-{
-  int i,j,k,m;
-
-  m = 0;
-  for (i = 0; i < n; i++) {
-    j = list[i];
-    rho0[j] += buf[m++];
-    arho2b[j] += buf[m++];
-    arho1[j][0] += buf[m++];
-    arho1[j][1] += buf[m++];
-    arho1[j][2] += buf[m++];
-    arho2[j][0] += buf[m++];
-    arho2[j][1] += buf[m++];
-    arho2[j][2] += buf[m++];
-    arho2[j][3] += buf[m++];
-    arho2[j][4] += buf[m++];
-    arho2[j][5] += buf[m++];
-    for (k = 0; k < 10; k++) arho3[j][k] += buf[m++];
-    arho3b[j][0] += buf[m++];
-    arho3b[j][1] += buf[m++];
-    arho3b[j][2] += buf[m++];
-    t_ave[j][0] += buf[m++];
-    t_ave[j][1] += buf[m++];
-    t_ave[j][2] += buf[m++];
-    tsq_ave[j][0] += buf[m++];
-    tsq_ave[j][1] += buf[m++];
-    tsq_ave[j][2] += buf[m++];
-  }
-}
-
-/* ----------------------------------------------------------------------
-   memory usage of local atom-based arrays 
-------------------------------------------------------------------------- */
-
-double PairMEAM::memory_usage()
-{
-  double bytes = 11 * nmax * sizeof(double);
-  bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double);
-  bytes += 3 * maxneigh * sizeof(double);
-  return bytes;
-}
-
-/* ----------------------------------------------------------------------
-   toggle neighbor list indices between zero- and one-based values
-   needed for access by MEAM Fortran library
-------------------------------------------------------------------------- */
-
-void PairMEAM::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh)
-{
-  int i,j,ii,jnum;
-  int *jlist;
-
-  for (ii = 0; ii < inum; ii++) {
-    i = ilist[ii];
-    jlist = firstneigh[i];
-    jnum = numneigh[i];
-    for (j = 0; j < jnum; j++) jlist[j]--;
-  }
-}
-
-void PairMEAM::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh)
-{
-  int i,j,ii,jnum;
-  int *jlist;
-
-  for (ii = 0; ii < inum; ii++) {
-    i = ilist[ii];
-    jlist = firstneigh[i];
-    jnum = numneigh[i];
-    for (j = 0; j < jnum; j++) jlist[j]++;
-  }
-}
+/* ----------------------------------------------------------------------
+   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: Greg Wagner (SNL)
+------------------------------------------------------------------------- */
+
+#include "math.h"
+#include "stdio.h"
+#include "stdlib.h"
+#include "string.h"
+#include "pair_meam.h"
+#include "atom.h"
+#include "force.h"
+#include "comm.h"
+#include "memory.h"
+#include "neighbor.h"
+#include "neigh_list.h"
+#include "neigh_request.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
+
+enum{FCC,BCC,HCP,DIM,DIAMOND,B1,C11,L12};
+int nkeywords = 16;
+char *keywords[] = {"Ec","alpha","rho0","delta","lattce",
+		    "attrac","repuls","nn2","Cmin","Cmax","rc","delr",
+		    "augt1","gsmooth_factor","re","ialloy"};
+
+/* ---------------------------------------------------------------------- */
+
+PairMEAM::PairMEAM(LAMMPS *lmp) : Pair(lmp)
+{
+  single_enable = 0;
+  one_coeff = 1;
+
+  nmax = 0;
+  rho = rho0 = rho1 = rho2 = rho3 = frhop = NULL;
+  gamma = dgamma1 = dgamma2 = dgamma3 = arho2b = NULL;
+  arho1 = arho2 = arho3 = arho3b = t_ave = tsq_ave = NULL;
+
+  maxneigh = 0;
+  scrfcn = dscrfcn = fcpair = NULL;
+
+  nelements = 0;
+  elements = NULL;
+  mass = NULL;
+
+  // set comm size needed by this Pair
+
+  comm_forward = 35;
+  comm_reverse = 27;
+}
+
+/* ----------------------------------------------------------------------
+   free all arrays
+   check if allocated, since class can be destructed when incomplete
+------------------------------------------------------------------------- */
+
+PairMEAM::~PairMEAM()
+{
+  meam_cleanup_();
+
+  memory->sfree(rho);
+  memory->sfree(rho0);
+  memory->sfree(rho1);
+  memory->sfree(rho2);
+  memory->sfree(rho3);
+  memory->sfree(frhop);
+  memory->sfree(gamma);
+  memory->sfree(dgamma1);
+  memory->sfree(dgamma2);
+  memory->sfree(dgamma3);
+  memory->sfree(arho2b);
+
+  memory->destroy_2d_double_array(arho1);
+  memory->destroy_2d_double_array(arho2);
+  memory->destroy_2d_double_array(arho3);
+  memory->destroy_2d_double_array(arho3b);
+  memory->destroy_2d_double_array(t_ave);
+  memory->destroy_2d_double_array(tsq_ave);
+
+  memory->sfree(scrfcn);
+  memory->sfree(dscrfcn);
+  memory->sfree(fcpair);
+  
+  for (int i = 0; i < nelements; i++) delete [] elements[i];
+  delete [] elements;
+  delete [] mass;
+
+  if (allocated) {
+    memory->destroy_2d_int_array(setflag);
+    memory->destroy_2d_double_array(cutsq);
+    delete [] map;
+    delete [] fmap;
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAM::compute(int eflag, int vflag)
+{
+  int i,j,ii,n,inum_half,itype,jtype,errorflag;
+  double evdwl;
+  int *ilist_half,*jlist_half,*numneigh_half,**firstneigh_half;
+  int *numneigh_full,**firstneigh_full;
+
+  evdwl = 0.0;
+  if (eflag || vflag) ev_setup(eflag,vflag);
+  else evflag = vflag_fdotr = 0;
+
+  int newton_pair = force->newton_pair;
+
+  // grow local arrays if necessary
+
+  if (atom->nmax > nmax) {
+    memory->sfree(rho);
+    memory->sfree(rho0);
+    memory->sfree(rho1);
+    memory->sfree(rho2);
+    memory->sfree(rho3);
+    memory->sfree(frhop);
+    memory->sfree(gamma);
+    memory->sfree(dgamma1);
+    memory->sfree(dgamma2);
+    memory->sfree(dgamma3);
+    memory->sfree(arho2b);
+    memory->destroy_2d_double_array(arho1);
+    memory->destroy_2d_double_array(arho2);
+    memory->destroy_2d_double_array(arho3);
+    memory->destroy_2d_double_array(arho3b);
+    memory->destroy_2d_double_array(t_ave);
+    memory->destroy_2d_double_array(tsq_ave);
+
+    nmax = atom->nmax;
+
+    rho = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho");
+    rho0 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho0");
+    rho1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho1");
+    rho2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho2");
+    rho3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:rho3");
+    frhop = (double *) memory->smalloc(nmax*sizeof(double),"pair:frhop");
+    gamma = (double *) memory->smalloc(nmax*sizeof(double),"pair:gamma");
+    dgamma1 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma1");
+    dgamma2 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma2");
+    dgamma3 = (double *) memory->smalloc(nmax*sizeof(double),"pair:dgamma3");
+    arho2b = (double *) memory->smalloc(nmax*sizeof(double),"pair:arho2b");
+    arho1 = memory->create_2d_double_array(nmax,3,"pair:arho1");
+    arho2 = memory->create_2d_double_array(nmax,6,"pair:arho2");
+    arho3 = memory->create_2d_double_array(nmax,10,"pair:arho3");
+    arho3b = memory->create_2d_double_array(nmax,3,"pair:arho3b");
+    t_ave = memory->create_2d_double_array(nmax,3,"pair:t_ave");
+    tsq_ave = memory->create_2d_double_array(nmax,3,"pair:tsq_ave");
+  }
+
+  // neighbor list info
+
+  inum_half = listhalf->inum;
+  ilist_half = listhalf->ilist;
+  numneigh_half = listhalf->numneigh;
+  firstneigh_half = listhalf->firstneigh;
+  numneigh_full = listfull->numneigh;
+  firstneigh_full = listfull->firstneigh;
+
+  // check size of scrfcn based on half neighbor list
+
+  int nlocal = atom->nlocal;
+  int nall = nlocal + atom->nghost;
+
+  n = 0;
+  for (ii = 0; ii < inum_half; ii++) n += numneigh_half[ilist_half[ii]];
+
+  if (n > maxneigh) {
+    memory->sfree(scrfcn);
+    memory->sfree(dscrfcn);
+    memory->sfree(fcpair);
+    maxneigh = n;
+    scrfcn =
+      (double *) memory->smalloc(maxneigh*sizeof(double),"pair:scrfcn");
+    dscrfcn = 
+      (double *) memory->smalloc(maxneigh*sizeof(double),"pair:dscrfcn");
+    fcpair = 
+      (double *) memory->smalloc(maxneigh*sizeof(double),"pair:fcpair");
+  }
+
+  // zero out local arrays
+
+  for (i = 0; i < nall; i++) {
+    rho0[i] = 0.0;
+    arho2b[i] = 0.0;
+    arho1[i][0] = arho1[i][1] = arho1[i][2] = 0.0;
+    for (j = 0; j < 6; j++) arho2[i][j] = 0.0;
+    for (j = 0; j < 10; j++) arho3[i][j] = 0.0;
+    arho3b[i][0] = arho3b[i][1] = arho3b[i][2] = 0.0;
+    t_ave[i][0] = t_ave[i][1] = t_ave[i][2] = 0.0;
+    tsq_ave[i][0] = tsq_ave[i][1] = tsq_ave[i][2] = 0.0;
+  }
+
+  double **x = atom->x;
+  double **f = atom->f;
+  int *type = atom->type;
+  int ntype = atom->ntypes;
+
+  // change neighbor list indices to Fortran indexing
+
+  neigh_c2f(inum_half,ilist_half,numneigh_half,firstneigh_half);
+  neigh_c2f(inum_half,ilist_half,numneigh_full,firstneigh_full);
+
+  // 3 stages of MEAM calculation
+  // loop over my atoms followed by communication
+
+  int ifort;
+  int offset = 0;
+  errorflag = 0;
+
+  for (ii = 0; ii < inum_half; ii++) {
+    i = ilist_half[ii];
+    ifort = i+1;
+    meam_dens_init_(&ifort,&nmax,&ntype,type,fmap,&x[0][0],
+		    &numneigh_half[i],firstneigh_half[i],
+		    &numneigh_full[i],firstneigh_full[i],
+		    &scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
+		    rho0,&arho1[0][0],&arho2[0][0],arho2b,
+		    &arho3[0][0],&arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],
+                    &errorflag);
+    if (errorflag) {
+      char str[128];
+      sprintf(str,"MEAM library error %d",errorflag);
+      error->one(str);
+    }
+    offset += numneigh_half[i];
+  }
+
+  comm->reverse_comm_pair(this);
+
+  meam_dens_final_(&nlocal,&nmax,&eflag_either,&eflag_global,&eflag_atom,
+		   &eng_vdwl,eatom,&ntype,type,fmap,
+		   &arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],
+		   &arho3b[0][0],&t_ave[0][0],&tsq_ave[0][0],gamma,dgamma1,
+		   dgamma2,dgamma3,rho,rho0,rho1,rho2,rho3,frhop,&errorflag);
+  if (errorflag) {
+    char str[128];
+    sprintf(str,"MEAM library error %d",errorflag);
+    error->one(str);
+  }
+
+  comm->forward_comm_pair(this);
+
+  offset = 0;
+
+  // vptr is first value in vatom if it will be used by meam_force()
+  // else vatom may not exist, so pass dummy ptr
+
+  double *vptr;
+  if (vflag_atom) vptr = &vatom[0][0];
+  else vptr = &cutmax;
+
+  for (ii = 0; ii < inum_half; ii++) {
+    i = ilist_half[ii];
+    ifort = i+1;
+    meam_force_(&ifort,&nmax,&eflag_either,&eflag_global,&eflag_atom,
+		&vflag_atom,&eng_vdwl,eatom,&ntype,type,fmap,&x[0][0],
+		&numneigh_half[i],firstneigh_half[i],
+		&numneigh_full[i],firstneigh_full[i],
+		&scrfcn[offset],&dscrfcn[offset],&fcpair[offset],
+		dgamma1,dgamma2,dgamma3,rho0,rho1,rho2,rho3,frhop,
+		&arho1[0][0],&arho2[0][0],arho2b,&arho3[0][0],&arho3b[0][0],
+		&t_ave[0][0],&tsq_ave[0][0],&f[0][0],vptr,&errorflag);
+    if (errorflag) {
+      char str[128];
+      sprintf(str,"MEAM library error %d",errorflag);
+      error->one(str);
+    }
+    offset += numneigh_half[i];
+  }
+
+  // change neighbor list indices back to C indexing
+
+  neigh_f2c(inum_half,ilist_half,numneigh_half,firstneigh_half);
+  neigh_f2c(inum_half,ilist_half,numneigh_full,firstneigh_full);
+
+  if (vflag_fdotr) virial_compute();
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAM::allocate()
+{
+  allocated = 1;
+  int n = atom->ntypes;
+
+  setflag = memory->create_2d_int_array(n+1,n+1,"pair:setflag");
+  cutsq = memory->create_2d_double_array(n+1,n+1,"pair:cutsq");
+
+  map = new int[n+1];
+  fmap = new int[n];
+}
+
+/* ----------------------------------------------------------------------
+   global settings 
+------------------------------------------------------------------------- */
+
+void PairMEAM::settings(int narg, char **arg)
+{
+  if (narg != 0) error->all("Illegal pair_style command");
+}
+
+/* ----------------------------------------------------------------------
+   set coeffs for one or more type pairs
+------------------------------------------------------------------------- */
+
+void PairMEAM::coeff(int narg, char **arg)
+{
+  int i,j,m,n;
+
+  if (!allocated) allocate();
+
+  if (narg < 6) 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 MEAM element names between 2 filenames
+  // nelements = # of MEAM elements
+  // elements = list of unique element names
+
+  if (nelements) {
+    for (i = 0; i < nelements; i++) delete [] elements[i];
+    delete [] elements;
+    delete [] mass;
+  }
+  nelements = narg - 4 - atom->ntypes;
+  if (nelements < 1) error->all("Incorrect args for pair coefficients");
+  elements = new char*[nelements];
+  mass = new double[nelements];
+  
+  for (i = 0; i < nelements; i++) {
+    n = strlen(arg[i+3]) + 1;
+    elements[i] = new char[n];
+    strcpy(elements[i],arg[i+3]);
+  }
+
+  // read MEAM library and parameter files
+  // pass all parameters to MEAM package
+  // tell MEAM package that setup is done
+
+  read_files(arg[2],arg[2+nelements+1]);
+  meam_setup_done_(&cutmax);
+  
+  // read args that map atom types to MEAM elements
+  // map[i] = which element the Ith atom type is, -1 if not mapped
+
+  for (i = 4 + nelements; i < narg; i++) {
+    m = i - (4+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 (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      setflag[i][j] = 0;
+
+  // set setflag i,j for type pairs where both are mapped to elements
+  // set mass for i,i in atom class
+
+  int count = 0;
+  for (int i = 1; i <= n; i++)
+    for (int j = i; j <= n; j++)
+      if (map[i] >= 0 && map[j] >= 0) {
+	setflag[i][j] = 1;
+	if (i == j) atom->set_mass(i,mass[map[i]]);
+	count++;
+      }
+
+  if (count == 0) error->all("Incorrect args for pair coefficients");
+}
+
+/* ----------------------------------------------------------------------
+   init specific to this pair style
+------------------------------------------------------------------------- */
+
+void PairMEAM::init_style()
+{
+  if (force->newton_pair == 0)
+    error->all("Pair style MEAM requires newton pair on");
+
+  // need full and half neighbor list
+
+  int irequest_full = neighbor->request(this);
+  neighbor->requests[irequest_full]->id = 1;
+  neighbor->requests[irequest_full]->half = 0;
+  neighbor->requests[irequest_full]->full = 1;
+  int irequest_half = neighbor->request(this);
+  neighbor->requests[irequest_half]->id = 2;
+  neighbor->requests[irequest_half]->half = 0;
+  neighbor->requests[irequest_half]->half_from_full = 1;
+  neighbor->requests[irequest_half]->otherlist = irequest_full;
+
+  // setup Fortran-style mapping array needed by MEAM package
+  // fmap is indexed from 1:ntypes by Fortran and stores a Fortran index
+  // if type I is not a MEAM atom, fmap stores a 0
+
+  for (int i = 1; i <= atom->ntypes; i++) fmap[i-1] = map[i] + 1;
+}
+
+/* ----------------------------------------------------------------------
+   neighbor callback to inform pair style of neighbor list to use
+   half or full
+------------------------------------------------------------------------- */
+
+void PairMEAM::init_list(int id, NeighList *ptr)
+{
+  if (id == 1) listfull = ptr;
+  else if (id == 2) listhalf = ptr;
+}
+
+/* ----------------------------------------------------------------------
+   init for one type pair i,j and corresponding j,i
+------------------------------------------------------------------------- */
+
+double PairMEAM::init_one(int i, int j)
+{
+  return cutmax;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAM::read_files(char *globalfile, char *userfile)
+{
+  // open global meamf file on proc 0
+
+  FILE *fp;
+  if (comm->me == 0) {
+    fp = fopen(globalfile,"r");
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open MEAM potential file %s",globalfile);
+      error->one(str);
+    }
+  }
+
+  // allocate parameter arrays
+
+  int params_per_line = 19;
+
+  int *lat = new int[nelements];
+  int *ielement = new int[nelements];
+  int *ibar = new int[nelements];
+  double *z = new double[nelements];
+  double *atwt = new double[nelements];
+  double *alpha = new double[nelements];
+  double *b0 = new double[nelements];
+  double *b1 = new double[nelements];
+  double *b2 = new double[nelements];
+  double *b3 = new double[nelements];
+  double *alat = new double[nelements];
+  double *esub = new double[nelements];
+  double *asub = new double[nelements];
+  double *t0 = new double[nelements];
+  double *t1 = new double[nelements];
+  double *t2 = new double[nelements];
+  double *t3 = new double[nelements];
+  double *rozero = new double[nelements];
+
+  bool *found = new bool[nelements];
+  for (int i = 0; i < nelements; i++) found[i] = false;
+
+  // read each set of params from global MEAM file
+  // one set of params can span multiple lines
+  // store params if element name is in element list
+  // if element name appears multiple times, only store 1st entry
+
+  int i,n,nwords;
+  char **words = new char*[params_per_line+1];
+  char line[MAXLINE],*ptr;
+  int eof = 0;
+
+  int nset = 0;
+  while (1) {
+    if (comm->me == 0) {
+      ptr = fgets(line,MAXLINE,fp);
+      if (ptr == NULL) {
+	eof = 1;
+	fclose(fp);
+      } else n = strlen(line) + 1;
+    }
+    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    if (eof) break;
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+    // strip comment, skip line if blank
+
+    if (ptr = strchr(line,'#')) *ptr = '\0';
+    nwords = atom->count_words(line);
+    if (nwords == 0) continue;
+
+    // concatenate additional lines until have params_per_line words
+
+    while (nwords < params_per_line) {
+      n = strlen(line);
+      if (comm->me == 0) {
+        ptr = fgets(&line[n],MAXLINE-n,fp);
+        if (ptr == NULL) {
+	  eof = 1;
+	  fclose(fp);
+        } else n = strlen(line) + 1;
+      }
+      MPI_Bcast(&eof,1,MPI_INT,0,world);
+      if (eof) break;
+      MPI_Bcast(&n,1,MPI_INT,0,world);
+      MPI_Bcast(line,n,MPI_CHAR,0,world);
+      if (ptr = strchr(line,'#')) *ptr = '\0';
+      nwords = atom->count_words(line);
+    }
+
+    if (nwords != params_per_line)
+      error->all("Incorrect format in MEAM potential file");
+
+    // words = ptrs to all words in line
+    // strip single and double quotes from words
+
+    nwords = 0;
+    words[nwords++] = strtok(line,"' \t\n\r\f");
+    while (words[nwords++] = strtok(NULL,"' \t\n\r\f")) continue;
+
+    // skip if element name isn't in element list
+
+    for (i = 0; i < nelements; i++)
+      if (strcmp(words[0],elements[i]) == 0) break;
+    if (i == nelements) continue;
+
+    // skip if element already appeared
+
+    if (found[i] == true) continue;
+    found[i] = true;
+
+    // map lat string to an integer
+
+    if (strcmp(words[1],"fcc") == 0) lat[i] = FCC;
+    else if (strcmp(words[1],"bcc") == 0) lat[i] = BCC;
+    else if (strcmp(words[1],"hcp") == 0) lat[i] = HCP;
+    else if (strcmp(words[1],"dim") == 0) lat[i] = DIM;
+    else if (strcmp(words[1],"dia") == 0) lat[i] = DIAMOND;
+    else error->all("Unrecognized lattice type in MEAM file 1");
+
+    // store parameters
+
+    z[i] = atof(words[2]);
+    ielement[i] = atoi(words[3]);
+    atwt[i] = atof(words[4]);
+    alpha[i] = atof(words[5]);
+    b0[i] = atof(words[6]);
+    b1[i] = atof(words[7]);
+    b2[i] = atof(words[8]);
+    b3[i] = atof(words[9]);
+    alat[i] = atof(words[10]);
+    esub[i] = atof(words[11]);
+    asub[i] = atof(words[12]);
+    t0[i] = atof(words[13]);
+    t1[i] = atof(words[14]);
+    t2[i] = atof(words[15]);
+    t3[i] = atof(words[16]);
+    rozero[i] = atof(words[17]);
+    ibar[i] = atoi(words[18]);
+
+    nset++;
+  }
+
+  // error if didn't find all elements in file
+
+  if (nset != nelements)
+    error->all("Did not find all elements in MEAM library file");
+
+  // pass element parameters to MEAM package
+
+  meam_setup_global_(&nelements,lat,z,ielement,atwt,alpha,b0,b1,b2,b3,
+  		     alat,esub,asub,t0,t1,t2,t3,rozero,ibar);
+
+  // set element masses
+
+  for (i = 0; i < nelements; i++) mass[i] = atwt[i];
+
+  // clean-up memory
+
+  delete [] words;
+
+  delete [] lat;
+  delete [] ielement;
+  delete [] ibar;
+  delete [] z;
+  delete [] atwt;
+  delete [] alpha;
+  delete [] b0;
+  delete [] b1;
+  delete [] b2;
+  delete [] b3;
+  delete [] alat;
+  delete [] esub;
+  delete [] asub;
+  delete [] t0;
+  delete [] t1;
+  delete [] t2;
+  delete [] t3;
+  delete [] rozero;
+  delete [] found;
+
+  // done if user param file is NULL
+
+  if (strcmp(userfile,"NULL") == 0) return;
+
+  // open user param file on proc 0
+
+  if (comm->me == 0) {
+    fp = fopen(userfile,"r");
+    if (fp == NULL) {
+      char str[128];
+      sprintf(str,"Cannot open MEAM potential file %s",userfile);
+      error->one(str);
+    }
+  }
+
+  // read settings
+  // pass them one at a time to MEAM package
+  // match strings to list of corresponding ints
+
+  int which;
+  double value;
+  int nindex,index[3];
+  int maxparams = 6;
+  char **params = new char*[maxparams];
+  int nparams;
+
+  eof = 0;
+  while (1) {
+    if (comm->me == 0) {
+      ptr = fgets(line,MAXLINE,fp);
+      if (ptr == NULL) {
+	eof = 1;
+	fclose(fp);
+      } else n = strlen(line) + 1;
+    }
+    MPI_Bcast(&eof,1,MPI_INT,0,world);
+    if (eof) break;
+    MPI_Bcast(&n,1,MPI_INT,0,world);
+    MPI_Bcast(line,n,MPI_CHAR,0,world);
+
+    // strip comment, skip line if blank
+
+    if (ptr = strchr(line,'#')) *ptr = '\0';
+    nparams = atom->count_words(line);
+    if (nparams == 0) continue;
+
+    // words = ptrs to all words in line
+
+    nparams = 0;
+    params[nparams++] = strtok(line,"=(), '\t\n\r\f");
+    while (nparams < maxparams && 
+	   (params[nparams++] = strtok(NULL,"=(), '\t\n\r\f")))
+      continue;
+    nparams--;
+
+    for (which = 0; which < nkeywords; which++)
+      if (strcmp(params[0],keywords[which]) == 0) break;
+    if (which == nkeywords) {
+      char str[128];
+      sprintf(str,"Keyword %s in MEAM parameter file not recognized",
+	      params[0]);
+      error->all(str);
+    }
+    nindex = nparams - 2;
+    for (i = 0; i < nindex; i++) index[i] = atoi(params[i+1]);
+
+    // map lattce_meam value to an integer
+
+    if (which == 4) {
+      if (strcmp(params[nparams-1],"fcc") == 0) value = FCC;
+      else if (strcmp(params[nparams-1],"bcc") == 0) value = BCC;
+      else if (strcmp(params[nparams-1],"hcp") == 0) value = HCP;
+      else if (strcmp(params[nparams-1],"dim") == 0) value = DIM;
+      else if (strcmp(params[nparams-1],"dia") == 0) value = DIAMOND;
+      else if (strcmp(params[nparams-1],"b1")  == 0) value = B1;
+      else if (strcmp(params[nparams-1],"c11") == 0) value = C11;
+      else if (strcmp(params[nparams-1],"l12") == 0) value = L12;
+      else error->all("Unrecognized lattice type in MEAM file 2");
+    }
+    else value = atof(params[nparams-1]);
+
+    // pass single setting to MEAM package
+
+    int errorflag = 0;
+    meam_setup_param_(&which,&value,&nindex,index,&errorflag);
+    if (errorflag) {
+      char str[128];
+      sprintf(str,"MEAM library error %d",errorflag);
+      error->all(str);
+    }
+  }
+
+  delete [] params;
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairMEAM::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc)
+{
+  int i,j,k,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    buf[m++] = rho0[j];
+    buf[m++] = rho1[j];
+    buf[m++] = rho2[j];
+    buf[m++] = rho3[j];
+    buf[m++] = frhop[j];
+    buf[m++] = gamma[j];
+    buf[m++] = dgamma1[j];
+    buf[m++] = dgamma2[j];
+    buf[m++] = dgamma3[j];
+    buf[m++] = arho2b[j];
+    buf[m++] = arho1[j][0];
+    buf[m++] = arho1[j][1];
+    buf[m++] = arho1[j][2];
+    buf[m++] = arho2[j][0];
+    buf[m++] = arho2[j][1];
+    buf[m++] = arho2[j][2];
+    buf[m++] = arho2[j][3];
+    buf[m++] = arho2[j][4];
+    buf[m++] = arho2[j][5];
+    for (k = 0; k < 10; k++) buf[m++] = arho3[j][k];
+    buf[m++] = arho3b[j][0];
+    buf[m++] = arho3b[j][1];
+    buf[m++] = arho3b[j][2];
+    buf[m++] = t_ave[j][0];
+    buf[m++] = t_ave[j][1];
+    buf[m++] = t_ave[j][2];
+    buf[m++] = tsq_ave[j][0];
+    buf[m++] = tsq_ave[j][1];
+    buf[m++] = tsq_ave[j][2];
+  }
+  return 38;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAM::unpack_comm(int n, int first, double *buf)
+{
+  int i,k,m,last;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    rho0[i] = buf[m++];
+    rho1[i] = buf[m++];
+    rho2[i] = buf[m++];
+    rho3[i] = buf[m++];
+    frhop[i] = buf[m++];
+    gamma[i] = buf[m++];
+    dgamma1[i] = buf[m++];
+    dgamma2[i] = buf[m++];
+    dgamma3[i] = buf[m++];
+    arho2b[i] = buf[m++];
+    arho1[i][0] = buf[m++];
+    arho1[i][1] = buf[m++];
+    arho1[i][2] = buf[m++];
+    arho2[i][0] = buf[m++];
+    arho2[i][1] = buf[m++];
+    arho2[i][2] = buf[m++];
+    arho2[i][3] = buf[m++];
+    arho2[i][4] = buf[m++];
+    arho2[i][5] = buf[m++];
+    for (k = 0; k < 10; k++) arho3[i][k] = buf[m++];
+    arho3b[i][0] = buf[m++];
+    arho3b[i][1] = buf[m++];
+    arho3b[i][2] = buf[m++];
+    t_ave[i][0] = buf[m++];
+    t_ave[i][1] = buf[m++];
+    t_ave[i][2] = buf[m++];
+    tsq_ave[i][0] = buf[m++];
+    tsq_ave[i][1] = buf[m++];
+    tsq_ave[i][2] = buf[m++];
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
+int PairMEAM::pack_reverse_comm(int n, int first, double *buf)
+{
+  int i,k,m,last,size;
+
+  m = 0;
+  last = first + n;
+  for (i = first; i < last; i++) {
+    buf[m++] = rho0[i];
+    buf[m++] = arho2b[i];
+    buf[m++] = arho1[i][0];
+    buf[m++] = arho1[i][1];
+    buf[m++] = arho1[i][2];
+    buf[m++] = arho2[i][0];
+    buf[m++] = arho2[i][1];
+    buf[m++] = arho2[i][2];
+    buf[m++] = arho2[i][3];
+    buf[m++] = arho2[i][4];
+    buf[m++] = arho2[i][5];
+    for (k = 0; k < 10; k++) buf[m++] = arho3[i][k];
+    buf[m++] = arho3b[i][0];
+    buf[m++] = arho3b[i][1];
+    buf[m++] = arho3b[i][2];
+    buf[m++] = t_ave[i][0];
+    buf[m++] = t_ave[i][1];
+    buf[m++] = t_ave[i][2];
+    buf[m++] = tsq_ave[i][0];
+    buf[m++] = tsq_ave[i][1];
+    buf[m++] = tsq_ave[i][2];
+  }
+
+  return 30;
+}
+
+/* ---------------------------------------------------------------------- */
+
+void PairMEAM::unpack_reverse_comm(int n, int *list, double *buf)
+{
+  int i,j,k,m;
+
+  m = 0;
+  for (i = 0; i < n; i++) {
+    j = list[i];
+    rho0[j] += buf[m++];
+    arho2b[j] += buf[m++];
+    arho1[j][0] += buf[m++];
+    arho1[j][1] += buf[m++];
+    arho1[j][2] += buf[m++];
+    arho2[j][0] += buf[m++];
+    arho2[j][1] += buf[m++];
+    arho2[j][2] += buf[m++];
+    arho2[j][3] += buf[m++];
+    arho2[j][4] += buf[m++];
+    arho2[j][5] += buf[m++];
+    for (k = 0; k < 10; k++) arho3[j][k] += buf[m++];
+    arho3b[j][0] += buf[m++];
+    arho3b[j][1] += buf[m++];
+    arho3b[j][2] += buf[m++];
+    t_ave[j][0] += buf[m++];
+    t_ave[j][1] += buf[m++];
+    t_ave[j][2] += buf[m++];
+    tsq_ave[j][0] += buf[m++];
+    tsq_ave[j][1] += buf[m++];
+    tsq_ave[j][2] += buf[m++];
+  }
+}
+
+/* ----------------------------------------------------------------------
+   memory usage of local atom-based arrays 
+------------------------------------------------------------------------- */
+
+double PairMEAM::memory_usage()
+{
+  double bytes = 11 * nmax * sizeof(double);
+  bytes += (3 + 6 + 10 + 3 + 3 + 3) * nmax * sizeof(double);
+  bytes += 3 * maxneigh * sizeof(double);
+  return bytes;
+}
+
+/* ----------------------------------------------------------------------
+   toggle neighbor list indices between zero- and one-based values
+   needed for access by MEAM Fortran library
+------------------------------------------------------------------------- */
+
+void PairMEAM::neigh_f2c(int inum, int *ilist, int *numneigh, int **firstneigh)
+{
+  int i,j,ii,jnum;
+  int *jlist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    for (j = 0; j < jnum; j++) jlist[j]--;
+  }
+}
+
+void PairMEAM::neigh_c2f(int inum, int *ilist, int *numneigh, int **firstneigh)
+{
+  int i,j,ii,jnum;
+  int *jlist;
+
+  for (ii = 0; ii < inum; ii++) {
+    i = ilist[ii];
+    jlist = firstneigh[i];
+    jnum = numneigh[i];
+    for (j = 0; j < jnum; j++) jlist[j]++;
+  }
+}
diff --git a/src/MEAM/pair_meam.h b/src/MEAM/pair_meam.h
index ba484f9fe3..c8a5f43648 100644
--- a/src/MEAM/pair_meam.h
+++ b/src/MEAM/pair_meam.h
@@ -1,97 +1,97 @@
-/* ----------------------------------------------------------------------
-   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.
-------------------------------------------------------------------------- */
-
-#ifndef PAIR_MEAM_H
-#define PAIR_MEAM_H
-
-extern "C" {
-  void meam_setup_global_(int *, int *, double *, int *, double *, double *,
-			 double *, double *, double *, double *, double *,
-			 double *, double *, double *, double *, double *,
-			 double *, double *, int *);
-  void meam_setup_param_(int *, double *, int *, int *, int *);
-  void meam_setup_done_(double *);
-
-  void meam_dens_init_(int *, int *, int *, int *, int *,
-		       double *, int *, int *, int *, int *,
-		       double *, double *, double *, double *,
-		       double *, double *,
-		       double *, double *, double *, double *, double *,
-                       int *);
-  
-  void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
-			int *, int *, int *,
-			double *, double *, double *, double *, 
-			double *, double *, double *,
-			double *, double *, double *, double *, 
-			double *, double *, 
-			double *, double *, double *, double *, int *);
- 
-  void meam_force_(int *, int *, int *, int *, int *, int *,
-		   double *, double *, int *, int *, int *, 
-		   double *, int *, int *, int *, int *, double *, double *,
-		   double *, double *, double *, double *, double *, double *,
-		   double *, double *, double *, double *, double *, double *,
-		   double *, double *, double *, double *, double *, double *, int *);
-
-  void meam_cleanup_();
-}
-
-
-#include "pair.h"
-
-namespace LAMMPS_NS {
-
-class PairMEAM : public Pair {
- public:
-  PairMEAM(class LAMMPS *);
-  ~PairMEAM();
-  void compute(int, int);
-  void settings(int, char **);
-  void coeff(int, char **);
-  void init_style();
-  void init_list(int, class NeighList *);
-  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 cutmax;                // max cutoff for all elements
-  int nelements;                // # of unique elements
-  char **elements;              // names of unique elements
-  double *mass;                 // mass of each element
-
-  int *map;                     // mapping from atom types to elements
-  int *fmap;                    // Fortran version of map array for MEAM lib
-
-  int maxneigh;
-  double *scrfcn,*dscrfcn,*fcpair;
-
-  int nmax;
-  double *rho,*rho0,*rho1,*rho2,*rho3,*frhop;
-  double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b;
-  double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave;
-
-  void allocate();
-  void read_files(char *, char *);
-  void neigh_f2c(int, int *, int *, int **);
-  void neigh_c2f(int, int *, int *, int **);
-};
-
-}
-
-#endif
+/* ----------------------------------------------------------------------
+   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.
+------------------------------------------------------------------------- */
+
+#ifndef PAIR_MEAM_H
+#define PAIR_MEAM_H
+
+extern "C" {
+  void meam_setup_global_(int *, int *, double *, int *, double *, double *,
+			 double *, double *, double *, double *, double *,
+			 double *, double *, double *, double *, double *,
+			 double *, double *, int *);
+  void meam_setup_param_(int *, double *, int *, int *, int *);
+  void meam_setup_done_(double *);
+
+  void meam_dens_init_(int *, int *, int *, int *, int *,
+		       double *, int *, int *, int *, int *,
+		       double *, double *, double *, double *,
+		       double *, double *,
+		       double *, double *, double *, double *, double *,
+                       int *);
+  
+  void meam_dens_final_(int *, int *, int *, int *, int *, double *, double *,
+			int *, int *, int *,
+			double *, double *, double *, double *, 
+			double *, double *, double *,
+			double *, double *, double *, double *, 
+			double *, double *, 
+			double *, double *, double *, double *, int *);
+ 
+  void meam_force_(int *, int *, int *, int *, int *, int *,
+		   double *, double *, int *, int *, int *, 
+		   double *, int *, int *, int *, int *, double *, double *,
+		   double *, double *, double *, double *, double *, double *,
+		   double *, double *, double *, double *, double *, double *,
+		   double *, double *, double *, double *, double *, double *, int *);
+
+  void meam_cleanup_();
+}
+
+
+#include "pair.h"
+
+namespace LAMMPS_NS {
+
+class PairMEAM : public Pair {
+ public:
+  PairMEAM(class LAMMPS *);
+  ~PairMEAM();
+  void compute(int, int);
+  void settings(int, char **);
+  void coeff(int, char **);
+  void init_style();
+  void init_list(int, class NeighList *);
+  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 cutmax;                // max cutoff for all elements
+  int nelements;                // # of unique elements
+  char **elements;              // names of unique elements
+  double *mass;                 // mass of each element
+
+  int *map;                     // mapping from atom types to elements
+  int *fmap;                    // Fortran version of map array for MEAM lib
+
+  int maxneigh;
+  double *scrfcn,*dscrfcn,*fcpair;
+
+  int nmax;
+  double *rho,*rho0,*rho1,*rho2,*rho3,*frhop;
+  double *gamma,*dgamma1,*dgamma2,*dgamma3,*arho2b;
+  double **arho1,**arho2,**arho3,**arho3b,**t_ave,**tsq_ave;
+
+  void allocate();
+  void read_files(char *, char *);
+  void neigh_f2c(int, int *, int *, int **);
+  void neigh_c2f(int, int *, int *, int **);
+};
+
+}
+
+#endif
diff --git a/src/MOLECULE/fix_bond_break.cpp b/src/MOLECULE/fix_bond_break.cpp
index 62f93e1ed9..0a656c80bc 100755
--- a/src/MOLECULE/fix_bond_break.cpp
+++ b/src/MOLECULE/fix_bond_break.cpp
@@ -163,7 +163,7 @@ void FixBondBreak::post_integrate()
 
   // need updated ghost atom positions
 
-  comm->communicate();
+  comm->forward_comm();
 
   // resize bond partner list and initialize it
   // probability array overlays distsq array
@@ -235,7 +235,7 @@ void FixBondBreak::post_integrate()
       if (partner[i]) probability[i] = random->uniform();
   }
 
-  comm->comm_fix(this);
+  comm->forward_comm_fix(this);
 
   // break bonds
   // if both atoms list each other as winning bond partner
diff --git a/src/MOLECULE/fix_bond_create.cpp b/src/MOLECULE/fix_bond_create.cpp
index f759b68523..ff0f2df8ee 100755
--- a/src/MOLECULE/fix_bond_create.cpp
+++ b/src/MOLECULE/fix_bond_create.cpp
@@ -273,12 +273,12 @@ void FixBondCreate::post_integrate()
 
   // need updated ghost atom positions
 
-  comm->communicate();
+  comm->forward_comm();
 
   // forward comm of bondcount, so ghosts have it
 
   commflag = 0;
-  comm->comm_fix(this);
+  comm->forward_comm_fix(this);
 
   // resize bond partner list and initialize it
   // probability array overlays distsq array
@@ -376,7 +376,7 @@ void FixBondCreate::post_integrate()
   }
  
   commflag = 1;
-  comm->comm_fix(this);
+  comm->forward_comm_fix(this);
 
   // create bonds for atoms I own
   // if other atom is owned by another proc, it should create same bond
diff --git a/src/Makefile.package b/src/Makefile.package
index 0c8aa9d850..e8b9bb6bd2 100644
--- a/src/Makefile.package
+++ b/src/Makefile.package
@@ -1,9 +1,9 @@
 # Settings for libraries used by specific LAMMPS packages
 # this file is auto-edited when those packages are included/excluded
 
-PKG_INC =  -I../../lib/reax -I../../lib/poems -I../../lib/meam 
-PKG_PATH = -L../../lib/reax -L../../lib/poems -L../../lib/meam 
-PKG_LIB =  -lreax -lpoems -lmeam 
+PKG_INC =  -I../../lib/atc -I../../lib/reax -I../../lib/poems -I../../lib/meam 
+PKG_PATH = -L../../lib/atc -L../../lib/reax -L../../lib/poems -L../../lib/meam -L../../lib/gpu 
+PKG_LIB =  -latc -lreax -lpoems -lmeam -lgpu 
 
-PKG_SYSPATH = $(reax_SYSPATH) $(meam_SYSPATH) 
-PKG_SYSLIB =  $(reax_SYSLIB) $(meam_SYSLIB) 
+PKG_SYSPATH = $(user-atc_SYSPATH) $(reax_SYSPATH) $(meam_SYSPATH) $(gpu_SYSPATH) 
+PKG_SYSLIB =  $(user-atc_SYSLIB) $(reax_SYSLIB) $(meam_SYSLIB) $(gpu_SYSLIB) 
diff --git a/src/OPT/pair_eam_opt.h b/src/OPT/pair_eam_opt.h
index b98c18de03..875cd5c885 100644
--- a/src/OPT/pair_eam_opt.h
+++ b/src/OPT/pair_eam_opt.h
@@ -224,7 +224,7 @@ void PairEAMOpt::eval()
   
   // communicate derivative of embedding function
   
-  comm->comm_pair(this);
+  comm->forward_comm_pair(this);
   
   // compute forces on each atom
   // loop over neighbors of my atoms
diff --git a/src/REAX/pair_reax.cpp b/src/REAX/pair_reax.cpp
index c14291f988..4d4990f234 100644
--- a/src/REAX/pair_reax.cpp
+++ b/src/REAX/pair_reax.cpp
@@ -168,7 +168,7 @@ void PairREAX::compute(int eflag, int vflag)
   // communicate local atomic bond order to ghost atomic bond order
 
   packflag = 0;
-  comm->comm_pair(this);
+  comm->forward_comm_pair(this);
 
   FORTRAN(molec, MOLEC)();
   FORTRAN(encalc, ENCALC)();
@@ -941,7 +941,7 @@ void PairREAX::cg_solve(const int & nlocal, const int & nghost,
   
   packflag = 1;
   comm->reverse_comm_pair(this);
-  comm->comm_pair(this);
+  comm->forward_comm_pair(this);
   
   MPI_Allreduce(&w[n-1], &sumtmp, 1, MPI_DOUBLE, MPI_SUM, world);
   w[n-1] = sumtmp;
@@ -978,7 +978,7 @@ void PairREAX::cg_solve(const int & nlocal, const int & nghost,
     }
       
     comm->reverse_comm_pair(this);
-    comm->comm_pair(this);
+    comm->forward_comm_pair(this);
 
     MPI_Allreduce(&w[n-1], &sumtmp, 1, MPI_DOUBLE, MPI_SUM, world);
     w[n-1] = sumtmp;
diff --git a/src/USER-ATC/fix_atc.cpp b/src/USER-ATC/fix_atc.cpp
index b3539e71da..215d78ab93 100644
--- a/src/USER-ATC/fix_atc.cpp
+++ b/src/USER-ATC/fix_atc.cpp
@@ -515,7 +515,7 @@ void FixATC::initial_integrate(int vflag)
 void FixATC::final_integrate()
 {
   // need updated ghost atom positions
-  comm->comm_fix(this);
+  comm->forward_comm_fix(this);
 
   try {
     atcTransfer_->pre_final_integrate();
diff --git a/src/USER-CD-EAM/pair_cdeam.cpp b/src/USER-CD-EAM/pair_cdeam.cpp
index 05028392af..9a4b9c9ee1 100644
--- a/src/USER-CD-EAM/pair_cdeam.cpp
+++ b/src/USER-CD-EAM/pair_cdeam.cpp
@@ -193,7 +193,7 @@ void PairCDEAM::compute(int eflag, int vflag)
 	// Communicate derivative of embedding function and densities
 	// and D_values (this for one-site formulation only).
 	communicationStage = 2;
-	comm->comm_pair(this);
+	comm->forward_comm_pair(this);
 
 	// The electron densities may not drop to zero because then the concentration would no longer be defined.
 	// But the concentration is not needed anyway if there is no interaction with another atom, which is the case
@@ -262,7 +262,7 @@ void PairCDEAM::compute(int eflag, int vflag)
 			comm->reverse_comm_pair(this);
 		}
 		communicationStage = 4;
-		comm->comm_pair(this);
+		comm->forward_comm_pair(this);
 	}
 
 	// Stage III
diff --git a/src/comm.cpp b/src/comm.cpp
index 8d85cd125c..ba676c8627 100644
--- a/src/comm.cpp
+++ b/src/comm.cpp
@@ -403,11 +403,11 @@ void Comm::setup()
 }
 
 /* ----------------------------------------------------------------------
-   communication of atom coords every timestep
+   forward communication of atom coords every timestep
    other per-atom attributes may also be sent via pack/unpack routines
 ------------------------------------------------------------------------- */
 
-void Comm::communicate()
+void Comm::forward_comm()
 {
   int n;
   MPI_Request request;
@@ -473,7 +473,7 @@ void Comm::communicate()
    other per-atom attributes may also be sent via pack/unpack routines
 ------------------------------------------------------------------------- */
       
-void Comm::reverse_communicate()
+void Comm::reverse_comm()
 {
   int n;
   MPI_Request request;
@@ -808,7 +808,7 @@ void Comm::borders()
    forward communication invoked by a Pair
 ------------------------------------------------------------------------- */
 
-void Comm::comm_pair(Pair *pair)
+void Comm::forward_comm_pair(Pair *pair)
 {
   int iswap,n;
   double *buf;
@@ -877,7 +877,7 @@ void Comm::reverse_comm_pair(Pair *pair)
    forward communication invoked by a Fix
 ------------------------------------------------------------------------- */
 
-void Comm::comm_fix(Fix *fix)
+void Comm::forward_comm_fix(Fix *fix)
 {
   int iswap,n;
   double *buf;
@@ -946,7 +946,7 @@ void Comm::reverse_comm_fix(Fix *fix)
    forward communication invoked by a Compute
 ------------------------------------------------------------------------- */
 
-void Comm::comm_compute(Compute *compute)
+void Comm::forward_comm_compute(Compute *compute)
 {
   int iswap,n;
   double *buf;
diff --git a/src/comm.h b/src/comm.h
index 79d60f5b4a..e51d14870f 100644
--- a/src/comm.h
+++ b/src/comm.h
@@ -43,16 +43,16 @@ class Comm : protected Pointers {
   void init();
   void set_procs();                 // setup 3d grid of procs
   void setup();                     // setup 3d communication pattern
-  void communicate();               // communication of atom coords
-  void reverse_communicate();       // reverse communication of forces
+  void forward_comm();              // forward communication of atom coords
+  void reverse_comm();              // reverse communication of forces
   void exchange();                  // move atoms to new procs
   void borders();                   // setup list of atoms to communicate
 
-  void comm_pair(class Pair *);                // forward comm from a Pair
+  void forward_comm_pair(class Pair *);        // forward comm from a Pair
   void reverse_comm_pair(class Pair *);        // reverse comm from a Pair
-  void comm_fix(class Fix *);                  // forward comm from a Fix
+  void forward_comm_fix(class Fix *);          // forward comm from a Fix
   void reverse_comm_fix(class Fix *);          // reverse comm from a Fix
-  void comm_compute(class Compute *);          // forward comm from a Compute
+  void forward_comm_compute(class Compute *);  // forward comm from a Compute
   void reverse_comm_compute(class Compute *);  // reverse comm from a Compute
 
   void irregular();                 // irregular communication across all procs
diff --git a/src/fix_orient_fcc.cpp b/src/fix_orient_fcc.cpp
index b87edcc578..5a71ac65c1 100644
--- a/src/fix_orient_fcc.cpp
+++ b/src/fix_orient_fcc.cpp
@@ -335,7 +335,7 @@ void FixOrientFCC::post_force(int vflag)
   
   // communicate to acquire nbr data for ghost atoms
 
-  comm->comm_fix(this);
+  comm->forward_comm_fix(this);
 
   // compute grain boundary force on each owned atom
   // skip atoms not in group
diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp
index b400c3eaa7..a7280f44dc 100644
--- a/src/fix_shake.cpp
+++ b/src/fix_shake.cpp
@@ -503,7 +503,7 @@ void FixShake::post_force(int vflag)
 
   // communicate results if necessary
 
-  if (nprocs > 1) comm->comm_fix(this);
+  if (nprocs > 1) comm->forward_comm_fix(this);
 
   // virial setup
 
@@ -2331,7 +2331,7 @@ void FixShake::post_force_respa(int vflag, int ilevel, int iloop)
 
   // communicate results if necessary
 
-  if (nprocs > 1) comm->comm_fix(this);
+  if (nprocs > 1) comm->forward_comm_fix(this);
 
   // virial setup
 
diff --git a/src/min.cpp b/src/min.cpp
index adf8fbeb95..2df35433f1 100644
--- a/src/min.cpp
+++ b/src/min.cpp
@@ -233,7 +233,7 @@ void Min::setup()
     force->kspace->compute(eflag,vflag);
   }
 
-  if (force->newton) comm->reverse_communicate();
+  if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
   output->setup(1);
@@ -287,7 +287,7 @@ void Min::setup_minimal(int flag)
     force->kspace->compute(eflag,vflag);
   }
 
-  if (force->newton) comm->reverse_communicate();
+  if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
 
@@ -394,7 +394,7 @@ double Min::energy_force(int resetflag)
 
   if (nflag == 0) {
     timer->stamp();
-    comm->communicate();
+    comm->forward_comm();
     timer->stamp(TIME_COMM);
   } else {
     if (modify->n_min_pre_exchange) modify->min_pre_exchange();
@@ -438,7 +438,7 @@ double Min::energy_force(int resetflag)
   }
 
   if (force->newton) {
-    comm->reverse_communicate();
+    comm->reverse_comm();
     timer->stamp(TIME_COMM);
   }
 
diff --git a/src/respa.cpp b/src/respa.cpp
index 143c5372ee..01bfdb35d9 100644
--- a/src/respa.cpp
+++ b/src/respa.cpp
@@ -353,7 +353,7 @@ void Respa::setup()
       force->kspace->setup();
       force->kspace->compute(eflag,vflag);
     }
-    if (newton[ilevel]) comm->reverse_communicate();
+    if (newton[ilevel]) comm->reverse_comm();
     copy_f_flevel(ilevel);
   }
   
@@ -413,7 +413,7 @@ void Respa::setup_minimal(int flag)
       force->kspace->setup();
       force->kspace->compute(eflag,vflag);
     }
-    if (newton[ilevel]) comm->reverse_communicate();
+    if (newton[ilevel]) comm->reverse_comm();
     copy_f_flevel(ilevel);
   }
   
@@ -513,7 +513,7 @@ void Respa::recurse(int ilevel)
 
     } else if (ilevel == 0) {
       timer->stamp();
-      comm->communicate();
+      comm->forward_comm();
       timer->stamp(TIME_COMM);
     }
 
@@ -560,7 +560,7 @@ void Respa::recurse(int ilevel)
     }
 
     if (newton[ilevel]) {
-      comm->reverse_communicate();
+      comm->reverse_comm();
       timer->stamp(TIME_COMM);
     }
   
diff --git a/src/style_user_ackland.h b/src/style_user_ackland.h
index e69de29bb2..6e7483a9f7 100644
--- a/src/style_user_ackland.h
+++ b/src/style_user_ackland.h
@@ -0,0 +1,20 @@
+/* ----------------------------------------------------------------------
+   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 ComputeInclude
+#include "compute_ackland_atom.h"
+#endif
+
+#ifdef ComputeClass
+ComputeStyle(ackland/atom,ComputeAcklandAtom)
+#endif
diff --git a/src/style_user_ewaldn.h b/src/style_user_ewaldn.h
index e69de29bb2..3eafa50744 100644
--- a/src/style_user_ewaldn.h
+++ b/src/style_user_ewaldn.h
@@ -0,0 +1,30 @@
+/* ----------------------------------------------------------------------
+   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 KSpaceInclude
+#include "ewald_n.h"
+#endif
+
+#ifdef KSpaceClass
+KSpaceStyle(ewald/n,EwaldN)
+#endif
+
+#ifdef PairInclude
+#include "pair_buck_coul.h"
+#include "pair_lj_coul.h"
+#endif
+
+#ifdef PairClass
+PairStyle(buck/coul,PairBuckCoul)
+PairStyle(lj/coul,PairLJCoul)
+#endif
diff --git a/src/verlet.cpp b/src/verlet.cpp
index e6670057f0..2f7ecaf8af 100644
--- a/src/verlet.cpp
+++ b/src/verlet.cpp
@@ -115,7 +115,7 @@ void Verlet::setup()
     force->kspace->compute(eflag,vflag);
   }
 
-  if (force->newton) comm->reverse_communicate();
+  if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
   output->setup(1);
@@ -165,7 +165,7 @@ void Verlet::setup_minimal(int flag)
     force->kspace->compute(eflag,vflag);
   }
 
-  if (force->newton) comm->reverse_communicate();
+  if (force->newton) comm->reverse_comm();
 
   modify->setup(vflag);
 }
@@ -201,7 +201,7 @@ void Verlet::run(int n)
 
     if (nflag == 0) {
       timer->stamp();
-      comm->communicate();
+      comm->forward_comm();
       timer->stamp(TIME_COMM);
     } else {
       if (n_pre_exchange) modify->pre_exchange();
@@ -250,7 +250,7 @@ void Verlet::run(int n)
     // reverse communication of forces
 
     if (force->newton) {
-      comm->reverse_communicate();
+      comm->reverse_comm();
       timer->stamp(TIME_COMM);
     }