diff --git a/src/MANYBODY/pair_bop.cpp b/src/MANYBODY/pair_bop.cpp
index 4b5f4c7d17..c98e87600e 100644
--- a/src/MANYBODY/pair_bop.cpp
+++ b/src/MANYBODY/pair_bop.cpp
@@ -46,10 +46,11 @@
 #include "memory.h"
 #include "error.h"
 #include "utils.h"
+#include "tokenizer.h"
+#include "potential_file_reader.h"
 
 using namespace LAMMPS_NS;
 
-#define MAXLINE 1024
 #define EPSILON 1.0e-6
 
 /* ---------------------------------------------------------------------- */
@@ -294,7 +295,6 @@ void PairBOP::compute(int eflag, int vflag)
   int nlocal = atom->nlocal;
   int nall = nlocal + atom->nghost;
   int cnt1;
-  MPI_Comm_rank(world,&me);
 
   inum = list->inum;
   ilist = list->ilist;
@@ -616,7 +616,6 @@ void PairBOP::coeff(int narg, char **arg)
 {
   int i,j;
   int n = atom->ntypes;
-  MPI_Comm_rank(world,&me);
 
   delete[] map;
   map = new int[n+1];
@@ -645,7 +644,7 @@ void PairBOP::coeff(int narg, char **arg)
   read_table(arg[2]);
 
   // match element names to BOP word types
-  if (me == 0) {
+  if (comm->me == 0) {
     for (i = 3; i < narg; i++) {
       if (strcmp(arg[i],"NULL") == 0) {
         map[i-2] = -1;
@@ -659,7 +658,7 @@ void PairBOP::coeff(int narg, char **arg)
   }
 
   MPI_Bcast(&map[1],atom->ntypes,MPI_INT,0,world);
-  if (me == 0) {
+  if (comm->me == 0) {
     if (elements) {
       for (i = 0; i < bop_types; i++) delete [] elements[i];
       delete [] elements;
@@ -1137,7 +1136,6 @@ double PairBOP::sigmaBo(int itmp, int jtmp)
   nlocal = atom->nlocal;
   ilist = list->ilist;
   firstneigh = list->firstneigh;
-  MPI_Comm_rank(world,&me);
 
   if(nb_sg==0) {
     nb_sg=(maxneigh)*(maxneigh/2);
@@ -4873,101 +4871,11 @@ double PairBOP::PiBo(int itmp, int jtmp)
 }
 
 /* ----------------------------------------------------------------------
-   read BOP potential file
+   allocate BOP tables
 ------------------------------------------------------------------------- */
 
-/* ---------------------------------------------------------------------- */
-
-void _noopt PairBOP::read_table(char *filename)
+void PairBOP::allocate_tables()
 {
-  int i,j,k,n,m;
-  int buf1,pass;
-  int nws,ws;
-  double buf2;
-  char s[MAXLINE],buf[2];
-
-  MPI_Comm_rank(world,&me);
-  if (me == 0) {
-    FILE *fp = force->open_potential(filename);
-    if (fp == NULL) {
-      char str[128];
-      snprintf(str,128,"Cannot open BOP potential file %s",filename);
-      error->one(FLERR,str);
-    }
-    utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);  // skip first comment line
-    utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-    sscanf(s,"%d",&bop_types);
-    elements = new char*[bop_types];
-    for(i=0;i<bop_types;i++) elements[i]=NULL;
-    for(i=0;i<bop_types;i++) {
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      nws=0;
-      ws=1;
-      for(j=0;j<(int)strlen(s);j++) {
-        if(ws==1) {
-          if(isspace(s[j])) {
-            ws=1;
-          } else {
-            ws=0;
-          }
-        } else {
-          if(isspace(s[j])) {
-            ws=1;
-            nws++;
-          } else {
-            ws=0;
-          }
-        }
-      }
-      if(nws!=3){
-         error->all(FLERR,"Incorrect table format check for element types");
-      }
-      sscanf(s,"%d %lf %s",&buf1,&buf2,buf);
-      n= strlen(buf)+1;
-      elements[i] = new char[n];
-      strcpy(elements[i],buf);
-    }
-    nws=0;
-    ws=1;
-    utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-    for(j=0;j<(int)strlen(s);j++) {
-      if(ws==1) {
-        if(isspace(s[j])) {
-          ws=1;
-        } else {
-          ws=0;
-        }
-      } else {
-        if(isspace(s[j])) {
-          ws=1;
-          nws++;
-        } else {
-          ws=0;
-        }
-      }
-    }
-    if (nws==3) {
-      sscanf(s,"%d %d %d",&nr,&ntheta,&nBOt);
-      npower=2;
-      if(ntheta<=10) npower=ntheta;
-    } else if (nws==2) {
-      sscanf(s,"%d %d",&nr,&nBOt);
-      ntheta=0;
-      npower=3;
-    } else {
-      error->one(FLERR,"Unsupported BOP potential file format");
-    }
-    fclose(fp);
-    npairs=bop_types*(bop_types+1)/2;
-  }
-
-  MPI_Bcast(&nr,1,MPI_INT,0,world);
-  MPI_Bcast(&nBOt,1,MPI_INT,0,world);
-  MPI_Bcast(&ntheta,1,MPI_INT,0,world);
-  MPI_Bcast(&bop_types,1,MPI_INT,0,world);
-  MPI_Bcast(&npairs,1,MPI_INT,0,world);
-  MPI_Bcast(&npower,1,MPI_INT,0,world);
-
   memory->destroy(pi_a);
   memory->destroy(pro_delta);
   memory->destroy(pi_delta);
@@ -5011,330 +4919,350 @@ void _noopt PairBOP::read_table(char *filename)
   memory->create(gfunc5,bop_types,bop_types,bop_types,ntheta,"BOP:gfunc5");
   memory->create(gfunc6,bop_types,bop_types,bop_types,ntheta,"BOP:gfunc6");
   memory->create(gpara,bop_types,bop_types,bop_types,npower+1,"BOP:gpara");
+}
 
-  allocate();
-  if (me == 0) {
-    FILE *fp = force->open_potential(filename);
-    if (fp == NULL) {
-      char str[128];
-      snprintf(str,128,"Cannot open BOP potential file %s",filename);
-      error->one(FLERR,str);
-    }
-    utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);  // skip first comment line
-    for(i=0;i<bop_types+2;i++) {
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-    }
-    utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-    sscanf(s,"%lf%lf%lf%lf%lf%lf%lf",&small1,&small2,&small3g
-        ,&small4,&small5,&small6,&small7);
-    for(i=0;i<bop_types;i++) {
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf",&pi_p[i]);
-    }
-    cutmax=0;
-    for(i=0;i<npairs;i++) {
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf",&rcut[i]);
-      if(rcut[i]>cutmax)
-        cutmax=rcut[i];
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf%lf%lf%lf",&sigma_c[i],&sigma_a[i],&pi_c[i],&pi_a[i]);
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf%lf",&sigma_delta[i],&pi_delta[i]);
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf%lf%lf",&sigma_f[i],&sigma_k[i],&small3[i]);
-    }
-    if(nws==3) {
-      for(i=0;i<bop_types;i++)
-        for(j=0;j<bop_types;j++)
-          for(k=j;k<bop_types;k++) {
-            if(npower<=2) {
-              for(m=0;m<ntheta;m++) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gfunc[j][i][k][n],&gfunc[j][i][k][n+1]
-                    ,&gfunc[j][i][k][n+2],&gfunc[j][i][k][n+3],&gfunc[j][i][k][n+4]);
-                n+=4;
+/* ----------------------------------------------------------------------
+   read BOP potential file
+------------------------------------------------------------------------- */
+
+/* ---------------------------------------------------------------------- */
+
+void _noopt PairBOP::read_table(char *filename)
+{
+  if (comm->me == 0) {
+    try {
+      PotentialFileReader reader(lmp, filename, "BOP");
+      char * line = nullptr;
+
+      bop_types = reader.next_int();
+      elements = new char*[bop_types];
+      for(int i=0; i < bop_types; i++) {
+        line = reader.next_line(3);
+        ValueTokenizer values(line);
+        values.next_int();
+        values.next_double();
+        std::string name = values.next_string();
+
+        elements[i] = new char[name.length()+1];
+        strcpy(elements[i], name.c_str());
+      }
+
+      line = reader.next_line();
+      ValueTokenizer values(line);
+      int format = values.count();
+
+      switch(format) {
+        case 3:
+          nr = values.next_int();
+          ntheta = values.next_int();
+          nBOt = values.next_int();
+          if(ntheta <= 10)
+            npower = ntheta;
+          else
+            npower = 2;
+          break;
+
+        case 2:
+          nr = values.next_int();
+          nBOt = values.next_int();
+          ntheta = 0;
+          npower = 3;
+          break;
+
+        default:
+          error->one(FLERR,"Unsupported BOP potential file format");
+      }
+
+      npairs = bop_types*(bop_types+1)/2;
+
+      allocate_tables();
+      allocate();
+
+      line = reader.next_line(7);
+      values = ValueTokenizer(line);
+
+      small1  = values.next_double();
+      small2  = values.next_double();
+      small3g = values.next_double();
+      small4  = values.next_double();
+      small5  = values.next_double();
+      small6  = values.next_double();
+      small7  = values.next_double();
+
+      for(int i = 0; i < bop_types; i++) {
+        pi_p[i] = reader.next_double();
+      }
+
+      cutmax = 0.0;
+
+      for(int i = 0; i < npairs; i++) {
+        rcut[i] = reader.next_double();
+        if (rcut[i] > cutmax)
+          cutmax = rcut[i];
+
+        line = reader.next_line(4);
+        values = ValueTokenizer(line);
+        sigma_c[i] = values.next_double();
+        sigma_a[i] = values.next_double();
+        pi_c[i]    = values.next_double();
+        pi_a[i]    = values.next_double();
+
+        line = reader.next_line(2);
+        values = ValueTokenizer(line);
+        sigma_delta[i] = values.next_double();
+        pi_delta[i]    = values.next_double();
+
+        line = reader.next_line(3);
+        values = ValueTokenizer(line);
+        sigma_f[i] = values.next_double();
+        sigma_k[i] = values.next_double();
+        small3[i]  = values.next_double();
+      }
+
+      if(format == 3) {
+        for(int i = 0; i < bop_types; i++)
+          for(int j = 0; j < bop_types; j++)
+            for(int k = j; k < bop_types; k++) {
+              if(npower <= 2) {
+                reader.next_dvector(ntheta, &gfunc[j][i][k][0]);
+              } else {
+                reader.next_dvector(npower+1, &gpara[j][i][k][0]);
+              }
+            }
+      } else {
+        for(int i = 0; i < bop_types; i++)
+          for(int j = 0; j < bop_types; j++)
+            for(int k = 0; k < bop_types; k++) {
+              reader.next_dvector(3, &gpara[i][j][k][0]);
+              gpara[j][i][k][3] = 0;
+            }
+      }
+
+      for(int i = 0; i < npairs; i++) {
+        reader.next_dvector(nr, &pRepul[i][0]);
+      }
+
+      for(int i = 0; i < npairs; i++) {
+        reader.next_dvector(nr, &pBetaS[i][0]);
+      }
+
+      for(int i = 0; i < npairs; i++) {
+        reader.next_dvector(nr, &pBetaP[i][0]);
+      }
+
+      for(int i = 0; i < npairs; i++) {
+        reader.next_dvector(nBOt, &FsigBO[i][0]);
+      }
+
+      for(int i = 0; i < bop_types; i++) {
+        pro_delta[i] = reader.next_double();
+      }
+
+      for(int i = 0; i < bop_types; i++) {
+        pro[i] = reader.next_double();
+      }
+
+      for(int i=0;i<npairs;i++) {
+        rcut3[i]=0.0;
+      }
+
+      if(format == 3) {
+        for(int i = 0; i < npairs; i++) {
+          rcut3[i] = reader.next_double();
+        }
+        for(int i = 0; i < npairs; i++) {
+          reader.next_dvector(nr, &pLong[i][0]);
+        }
+      }
+
+      rcutall=0.0;
+      for(int i=0; i < npairs; i++) {
+        if(rcut[i]>rcutall)
+          rcutall=rcut[i];
+        if(rcut3[i]>rcutall)
+          rcutall=rcut3[i];
+        rcutsq[i]=rcut[i]*rcut[i];
+        rcutsq3[i]=rcut3[i]*rcut3[i];
+        dr[i]=rcut[i]/((double)nr-1.0);
+        rdr[i]=1.0/dr[i];
+        dr3[i]=rcut3[i]/((double)nr-1.0);
+        rdr3[i]=1.0/dr3[i];
+      }
+
+      rctroot=rcutall;
+      dtheta=2.0/((double)ntheta-1.0);
+      rdtheta=1.0/dtheta;
+      dBO=1.0/((double)nBOt-1.0);
+      rdBO=1.0/(double)dBO;
+
+      for (int i = 0; i < npairs; i++) {
+        pBetaS1[i][0]=pBetaS[i][1]-pBetaS[i][0];
+        pBetaS1[i][1]=0.5*(pBetaS[i][2]-pBetaS[i][0]);
+        pBetaS1[i][nr-2]=0.5*(pBetaS[i][nr-1]-pBetaS[i][nr-3]);
+        pBetaS1[i][nr-1]=pBetaS[i][nr-1]-pBetaS[i][nr-2];
+        pBetaP1[i][0]=pBetaP[i][1]-pBetaP[i][0];
+        pBetaP1[i][1]=0.5*(pBetaP[i][2]-pBetaP[i][0]);
+        pBetaP1[i][nr-2]=0.5*(pBetaP[i][nr-1]-pBetaP[i][nr-3]);
+        pBetaP1[i][nr-1]=pBetaP[i][nr-1]-pBetaP[i][nr-2];
+        pRepul1[i][0]=pRepul[i][1]-pRepul[i][0];
+        pRepul1[i][1]=0.5*(pRepul[i][2]-pRepul[i][0]);
+        pRepul1[i][nr-2]=0.5*(pRepul[i][nr-1]-pRepul[i][nr-3]);
+        pRepul1[i][nr-1]=pRepul[i][nr-1]-pRepul[i][nr-2];
+        FsigBO1[i][0]=FsigBO[i][1]-FsigBO[i][0];
+        FsigBO1[i][1]=0.5*(FsigBO[i][2]-FsigBO[i][0]);
+        FsigBO1[i][nBOt-2]=0.5*(FsigBO[i][nBOt-1]-FsigBO[i][nBOt-3]);
+        FsigBO1[i][nBOt-1]=FsigBO[i][nBOt-1]-FsigBO[i][nBOt-2];
+        pLong1[i][0]=pLong[i][1]-pLong[i][0];
+        pLong1[i][1]=0.5*(pLong[i][2]-pLong[i][0]);
+        pLong1[i][nBOt-2]=0.5*(pLong[i][nr-1]-pLong[i][nr-3]);
+        pLong1[i][nBOt-1]=pLong[i][nr-1]-pLong[i][nr-2];
+        for (int k = 2; k < nr-2; k++) {
+          pBetaS1[i][k]=((pBetaS[i][k-2]-pBetaS[i][k+2])
+              +8.0*(pBetaS[i][k+1]-pBetaS[i][k-1]))/12.0;
+          pBetaP1[i][k]=((pBetaP[i][k-2]-pBetaP[i][k+2])
+              +8.0*(pBetaP[i][k+1]-pBetaP[i][k-1]))/12.0;
+          pRepul1[i][k]=((pRepul[i][k-2]-pRepul[i][k+2])
+              +8.0*(pRepul[i][k+1]-pRepul[i][k-1]))/12.0;
+          pLong1[i][k]=((pLong[i][k-2]-pLong[i][k+2])
+              +8.0*(pLong[i][k+1]-pLong[i][k-1]))/12.0;
+        }
+        for (int k=2; k < nr-2; k++) {
+          FsigBO1[i][k]=((FsigBO[i][k-2]-FsigBO[i][k+2])
+              +8.0*(FsigBO[i][k+1]-FsigBO[i][k-1]))/12.0;
+        }
+        for (int k = 0; k < nr-1; k++) {
+          pBetaS2[i][k]=3.0*(pBetaS[i][k+1]-pBetaS[i][k])
+              -2.0*pBetaS1[i][k]-pBetaS1[i][k+1];
+          pBetaS3[i][k]=pBetaS1[i][k]+pBetaS1[i][k+1]
+              -2.0*(pBetaS[i][k+1]-pBetaS[i][k]);
+          pBetaP2[i][k]=3.0*(pBetaP[i][k+1]-pBetaP[i][k])
+              -2.0*pBetaP1[i][k]-pBetaP1[i][k+1];
+          pBetaP3[i][k]=pBetaP1[i][k]+pBetaP1[i][k+1]
+              -2.0*(pBetaP[i][k+1]-pBetaP[i][k]);
+          pRepul2[i][k]=3.0*(pRepul[i][k+1]-pRepul[i][k])
+              -2.0*pRepul1[i][k]-pRepul1[i][k+1];
+          pRepul3[i][k]=pRepul1[i][k]+pRepul1[i][k+1]
+              -2.0*(pRepul[i][k+1]-pRepul[i][k]);
+          pLong2[i][k]=3.0*(pLong[i][k+1]-pLong[i][k])
+              -2.0*pLong1[i][k]-pLong1[i][k+1];
+          pLong3[i][k]=pLong1[i][k]+pLong1[i][k+1]
+              -2.0*(pLong[i][k+1]-pLong[i][k]);
+        }
+
+        for (int k = 0; k < nBOt-1; k++) {
+          FsigBO2[i][k]=3.0*(FsigBO[i][k+1]-FsigBO[i][k])
+              -2.0*FsigBO1[i][k]-FsigBO1[i][k+1];
+          FsigBO3[i][k]=FsigBO1[i][k]+FsigBO1[i][k+1]
+              -2.0*(FsigBO[i][k+1]-FsigBO[i][k]);
+        }
+
+        pBetaS2[i][nr-1]=0.0;
+        pBetaS3[i][nr-1]=0.0;
+        pBetaP2[i][nr-1]=0.0;
+        pBetaP3[i][nr-1]=0.0;
+        pRepul2[i][nr-1]=0.0;
+        pRepul3[i][nr-1]=0.0;
+        pLong2[i][nr-1]=0.0;
+        pLong3[i][nr-1]=0.0;
+        FsigBO2[i][nBOt-1]=0.0;
+        FsigBO3[i][nBOt-1]=0.0;
+
+        for (int k=0; k < nr; k++) {
+          pBetaS4[i][k]=pBetaS1[i][k]/dr[i];
+          pBetaS5[i][k]=2.0*pBetaS2[i][k]/dr[i];
+          pBetaS6[i][k]=3.0*pBetaS3[i][k]/dr[i];
+          pBetaP4[i][k]=pBetaP1[i][k]/dr[i];
+          pBetaP5[i][k]=2.0*pBetaP2[i][k]/dr[i];
+          pBetaP6[i][k]=3.0*pBetaP3[i][k]/dr[i];
+          pRepul4[i][k]=pRepul1[i][k]/dr[i];
+          pRepul5[i][k]=2.0*pRepul2[i][k]/dr[i];
+          pRepul6[i][k]=3.0*pRepul3[i][k]/dr[i];
+          pLong4[i][k]=pLong1[i][k]/dr3[i];
+          pLong5[i][k]=2.0*pLong2[i][k]/dr3[i];
+          pLong6[i][k]=3.0*pLong3[i][k]/dr3[i];
+        }
+        for (int k=0; k < nBOt; k++) {
+          FsigBO4[i][k]=FsigBO1[i][k]/dBO;
+          FsigBO5[i][k]=2.0*FsigBO2[i][k]/dBO;
+          FsigBO6[i][k]=3.0*FsigBO3[i][k]/dBO;
+        }
+      }
+
+      if (npower <= 2) {
+        for (int i = 0; i < bop_types; i++) {
+          for (int j = 0; j < bop_types; j++) {
+            for (int k = j; k < bop_types; k++) {
+              gfunc1[j][i][k][0] = gfunc[j][i][k][1] - gfunc[j][i][k][0];
+              gfunc1[j][i][k][1] = 0.5 * (gfunc[j][i][k][2] - gfunc[j][i][k][0]);
+              gfunc1[j][i][k][ntheta - 2] = 0.5 * (gfunc[j][i][k][ntheta - 1] - gfunc[j][i][k][ntheta - 3]);
+              gfunc1[j][i][k][ntheta - 1] = 0.5 * (gfunc[j][i][k][ntheta - 1] - gfunc[j][i][k][ntheta - 2]);
+
+              for (int m = 2; m < ntheta - 2; m++) {
+                gfunc1[j][i][k][m] = ((gfunc[j][i][k][m - 2] - gfunc[j][i][k][m + 2]) +
+                                      8.0 * (gfunc[j][i][k][m + 1] - gfunc[j][i][k][m + 1] - gfunc[j][i][k][m - 1])) /
+                                    12.0;
+              }
+
+              for (int m = 0; m < ntheta - 1; m++) {
+                gfunc2[j][i][k][m] = 3.0 * (gfunc[j][i][k][m + 1] - gfunc[j][i][k][m]) -
+                                    2.0 * gfunc1[j][i][k][m] - gfunc1[j][i][k][m + 1];
+                gfunc3[j][i][k][m] = gfunc1[j][i][k][m] + gfunc1[j][i][k][m + 1] -
+                                    2.0 * (gfunc[j][i][k][m + 1] - gfunc[j][i][k][m]);
+              }
+
+              gfunc2[j][i][k][ntheta - 1] = 0.0;
+              gfunc3[j][i][k][ntheta - 1] = 0.0;
+
+              for (int m = 0; m < ntheta; m++) {
+                gfunc4[j][i][k][ntheta - 1] = gfunc1[j][i][k][m] / dtheta;
+                gfunc5[j][i][k][ntheta - 1] = 2.0 * gfunc2[j][i][k][m] / dtheta;
+                gfunc6[j][i][k][ntheta - 1] = 3.0 * gfunc3[j][i][k][m] / dtheta;
+              }
+            }
+          }
+        }
+      }
+
+      for (int i = 0; i < bop_types; i++) {
+        for (int j = 0; j < bop_types; j++) {
+          for (int k = 0; k < j; k++) {
+            if (npower <= 2) {
+              for (int n = 0; n < ntheta; n++) {
+                gfunc[j][i][k][n] = gfunc[k][i][j][n];
+                gfunc1[j][i][k][n] = gfunc1[k][i][j][n];
+                gfunc2[j][i][k][n] = gfunc2[k][i][j][n];
+                gfunc3[j][i][k][n] = gfunc3[k][i][j][n];
+                gfunc4[j][i][k][n] = gfunc4[k][i][j][n];
+                gfunc5[j][i][k][n] = gfunc5[k][i][j][n];
+                gfunc6[j][i][k][n] = gfunc6[k][i][j][n];
               }
             } else {
-              if(npower==3) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3]);
+              for (int n = 0; n < npower + 1; n++) {
+                gpara[j][i][k][n] = gpara[k][i][j][n];
               }
-              else if(npower==4) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-              }
-              else if(npower==5) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf",&gpara[j][i][k][5]);
-              }
-              else if(npower==6) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6]);
-              }
-              else if(npower==7) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7]);
-              }
-              else if(npower==8) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7],&gpara[j][i][k][8]);
-              }
-              else if(npower==9) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7],&gpara[j][i][k][8],&gpara[j][i][k][9]);
-              }
-              else if(npower==10) {
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][0],&gpara[j][i][k][1],&gpara[j][i][k][2],&gpara[j][i][k][3],&gpara[j][i][k][4]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf%lf%lf%lf%lf",&gpara[j][i][k][5],&gpara[j][i][k][6],&gpara[j][i][k][7],&gpara[j][i][k][8],&gpara[j][i][k][9]);
-                utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-                sscanf(s,"%lf",&gpara[j][i][k][10]);
-              }
-            }
-          }
-    } else {
-      for(i=0;i<bop_types;i++)
-        for(j=0;j<bop_types;j++)
-          for(k=0;k<bop_types;k++) {
-            utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-            sscanf(s,"%lf%lf%lf",&gpara[i][j][k][0],&gpara[i][j][k][1],&gpara[i][j][k][2]);
-            gpara[j][i][k][3]=0;
-          }
-    }
-    for(i=0;i<npairs;i++) {
-      for(j=0;j<nr;j++) {
-        utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-        sscanf(s,"%lf%lf%lf%lf%lf",&pRepul[i][j],&pRepul[i][j+1]
-            ,&pRepul[i][j+2],&pRepul[i][j+3],&pRepul[i][j+4]);
-        j+=4;
-      }
-    }
-    for(i=0;i<npairs;i++) {
-      for(j=0;j<nr;j++) {
-        utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-        sscanf(s,"%lf%lf%lf%lf%lf",&pBetaS[i][j],&pBetaS[i][j+1]
-            ,&pBetaS[i][j+2],&pBetaS[i][j+3],&pBetaS[i][j+4]);
-        j+=4;
-      }
-    }
-    for(i=0;i<npairs;i++) {
-      for(j=0;j<nr;j++) {
-        utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-        sscanf(s,"%lf%lf%lf%lf%lf",&pBetaP[i][j],&pBetaP[i][j+1]
-            ,&pBetaP[i][j+2],&pBetaP[i][j+3],&pBetaP[i][j+4]);
-        j+=4;
-      }
-    }
-    for(i=0;i<npairs;i++) {
-      for(j=0;j<nBOt;j++) {
-        utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-        sscanf(s,"%lf%lf%lf%lf%lf",&FsigBO[i][j],&FsigBO[i][j+1]
-            ,&FsigBO[i][j+2],&FsigBO[i][j+3],&FsigBO[i][j+4]);
-        j+=4;
-      }
-    }
-    for(i=0;i<bop_types;i++) {
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf",&pro_delta[i]);
-    }
-    for(i=0;i<bop_types;i++) {
-      utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-      sscanf(s,"%lf",&pro[i]);
-    }
-    for(i=0;i<npairs;i++) {
-      rcut3[i]=0.0;
-    }
-    pass=0;
-    i=0;
-    if(nws==3) {
-      for(i=0;i<npairs;i++) {
-        utils::sfgets(FLERR,s,MAXLINE,fp,filename,error);
-        sscanf(s,"%lf",&rcut3[i]);
-      }
-      for(i=0;i<npairs;i++) {
-        for(j=0;j<nr;j++) {
-          pass=0;
-          while(fgets(s,MAXLINE,fp)!=NULL&&pass==0) {
-            sscanf(s,"%lf%lf%lf%lf%lf",&pLong[i][j],&pLong[i][j+1]
-              ,&pLong[i][j+2],&pLong[i][j+3],&pLong[i][j+4]);
-            j+=4;
-            pass=1;
-          }
-        }
-      }
-    }
-
-    rcutall=0.0;
-    for(i=0;i<npairs;i++) {
-      if(rcut[i]>rcutall)
-        rcutall=rcut[i];
-      if(rcut3[i]>rcutall)
-        rcutall=rcut3[i];
-      rcutsq[i]=rcut[i]*rcut[i];
-      rcutsq3[i]=rcut3[i]*rcut3[i];
-      dr[i]=rcut[i]/((double)nr-1.0);
-      rdr[i]=1.0/dr[i];
-      dr3[i]=rcut3[i]/((double)nr-1.0);
-      rdr3[i]=1.0/dr3[i];
-    }
-    rctroot=rcutall;
-    dtheta=2.0/((double)ntheta-1.0);
-    rdtheta=1.0/dtheta;
-    dBO=1.0/((double)nBOt-1.0);
-    rdBO=1.0/(double)dBO;
-    for(i=0;i<npairs;i++) {
-      pBetaS1[i][0]=pBetaS[i][1]-pBetaS[i][0];
-      pBetaS1[i][1]=0.5*(pBetaS[i][2]-pBetaS[i][0]);
-      pBetaS1[i][nr-2]=0.5*(pBetaS[i][nr-1]-pBetaS[i][nr-3]);
-      pBetaS1[i][nr-1]=pBetaS[i][nr-1]-pBetaS[i][nr-2];
-      pBetaP1[i][0]=pBetaP[i][1]-pBetaP[i][0];
-      pBetaP1[i][1]=0.5*(pBetaP[i][2]-pBetaP[i][0]);
-      pBetaP1[i][nr-2]=0.5*(pBetaP[i][nr-1]-pBetaP[i][nr-3]);
-      pBetaP1[i][nr-1]=pBetaP[i][nr-1]-pBetaP[i][nr-2];
-      pRepul1[i][0]=pRepul[i][1]-pRepul[i][0];
-      pRepul1[i][1]=0.5*(pRepul[i][2]-pRepul[i][0]);
-      pRepul1[i][nr-2]=0.5*(pRepul[i][nr-1]-pRepul[i][nr-3]);
-      pRepul1[i][nr-1]=pRepul[i][nr-1]-pRepul[i][nr-2];
-      FsigBO1[i][0]=FsigBO[i][1]-FsigBO[i][0];
-      FsigBO1[i][1]=0.5*(FsigBO[i][2]-FsigBO[i][0]);
-      FsigBO1[i][nBOt-2]=0.5*(FsigBO[i][nBOt-1]-FsigBO[i][nBOt-3]);
-      FsigBO1[i][nBOt-1]=FsigBO[i][nBOt-1]-FsigBO[i][nBOt-2];
-      pLong1[i][0]=pLong[i][1]-pLong[i][0];
-      pLong1[i][1]=0.5*(pLong[i][2]-pLong[i][0]);
-      pLong1[i][nBOt-2]=0.5*(pLong[i][nr-1]-pLong[i][nr-3]);
-      pLong1[i][nBOt-1]=pLong[i][nr-1]-pLong[i][nr-2];
-      for(k=2;k<nr-2;k++) {
-        pBetaS1[i][k]=((pBetaS[i][k-2]-pBetaS[i][k+2])
-            +8.0*(pBetaS[i][k+1]-pBetaS[i][k-1]))/12.0;
-        pBetaP1[i][k]=((pBetaP[i][k-2]-pBetaP[i][k+2])
-            +8.0*(pBetaP[i][k+1]-pBetaP[i][k-1]))/12.0;
-        pRepul1[i][k]=((pRepul[i][k-2]-pRepul[i][k+2])
-            +8.0*(pRepul[i][k+1]-pRepul[i][k-1]))/12.0;
-        pLong1[i][k]=((pLong[i][k-2]-pLong[i][k+2])
-            +8.0*(pLong[i][k+1]-pLong[i][k-1]))/12.0;
-      }
-      for(k=2;k<nr-2;k++) {
-        FsigBO1[i][k]=((FsigBO[i][k-2]-FsigBO[i][k+2])
-            +8.0*(FsigBO[i][k+1]-FsigBO[i][k-1]))/12.0;
-      }
-      for(k=0;k<nr-1;k++) {
-        pBetaS2[i][k]=3.0*(pBetaS[i][k+1]-pBetaS[i][k])
-            -2.0*pBetaS1[i][k]-pBetaS1[i][k+1];
-        pBetaS3[i][k]=pBetaS1[i][k]+pBetaS1[i][k+1]
-            -2.0*(pBetaS[i][k+1]-pBetaS[i][k]);
-        pBetaP2[i][k]=3.0*(pBetaP[i][k+1]-pBetaP[i][k])
-            -2.0*pBetaP1[i][k]-pBetaP1[i][k+1];
-        pBetaP3[i][k]=pBetaP1[i][k]+pBetaP1[i][k+1]
-            -2.0*(pBetaP[i][k+1]-pBetaP[i][k]);
-        pRepul2[i][k]=3.0*(pRepul[i][k+1]-pRepul[i][k])
-            -2.0*pRepul1[i][k]-pRepul1[i][k+1];
-        pRepul3[i][k]=pRepul1[i][k]+pRepul1[i][k+1]
-            -2.0*(pRepul[i][k+1]-pRepul[i][k]);
-        pLong2[i][k]=3.0*(pLong[i][k+1]-pLong[i][k])
-            -2.0*pLong1[i][k]-pLong1[i][k+1];
-        pLong3[i][k]=pLong1[i][k]+pLong1[i][k+1]
-            -2.0*(pLong[i][k+1]-pLong[i][k]);
-      }
-      for(k=0;k<nBOt-1;k++) {
-        FsigBO2[i][k]=3.0*(FsigBO[i][k+1]-FsigBO[i][k])
-            -2.0*FsigBO1[i][k]-FsigBO1[i][k+1];
-        FsigBO3[i][k]=FsigBO1[i][k]+FsigBO1[i][k+1]
-            -2.0*(FsigBO[i][k+1]-FsigBO[i][k]);
-      }
-      pBetaS2[i][nr-1]=0.0;
-      pBetaS3[i][nr-1]=0.0;
-      pBetaP2[i][nr-1]=0.0;
-      pBetaP3[i][nr-1]=0.0;
-      pRepul2[i][nr-1]=0.0;
-      pRepul3[i][nr-1]=0.0;
-      pLong2[i][nr-1]=0.0;
-      pLong3[i][nr-1]=0.0;
-      FsigBO2[i][nBOt-1]=0.0;
-      FsigBO3[i][nBOt-1]=0.0;
-      for(k=0;k<nr;k++) {
-        pBetaS4[i][k]=pBetaS1[i][k]/dr[i];
-        pBetaS5[i][k]=2.0*pBetaS2[i][k]/dr[i];
-        pBetaS6[i][k]=3.0*pBetaS3[i][k]/dr[i];
-        pBetaP4[i][k]=pBetaP1[i][k]/dr[i];
-        pBetaP5[i][k]=2.0*pBetaP2[i][k]/dr[i];
-        pBetaP6[i][k]=3.0*pBetaP3[i][k]/dr[i];
-        pRepul4[i][k]=pRepul1[i][k]/dr[i];
-        pRepul5[i][k]=2.0*pRepul2[i][k]/dr[i];
-        pRepul6[i][k]=3.0*pRepul3[i][k]/dr[i];
-        pLong4[i][k]=pLong1[i][k]/dr3[i];
-        pLong5[i][k]=2.0*pLong2[i][k]/dr3[i];
-        pLong6[i][k]=3.0*pLong3[i][k]/dr3[i];
-      }
-      for(k=0;k<nBOt;k++) {
-        FsigBO4[i][k]=FsigBO1[i][k]/dBO;
-        FsigBO5[i][k]=2.0*FsigBO2[i][k]/dBO;
-        FsigBO6[i][k]=3.0*FsigBO3[i][k]/dBO;
-      }
-    }
-    if(npower<=2) {
-      for(i=0;i<bop_types;i++) {
-        for(j=0;j<bop_types;j++) {
-          for(k=j;k<bop_types;k++) {
-            gfunc1[j][i][k][0]=gfunc[j][i][k][1]-gfunc[j][i][k][0];
-            gfunc1[j][i][k][1]=0.5*(gfunc[j][i][k][2]-gfunc[j][i][k][0]);
-            gfunc1[j][i][k][ntheta-2]=0.5*(gfunc[j][i][k][ntheta-1]-gfunc[j][i][k][ntheta-3]);
-            gfunc1[j][i][k][ntheta-1]=0.5*(gfunc[j][i][k][ntheta-1]-gfunc[j][i][k][ntheta-2]);
-            for(m=2;m<ntheta-2;m++) {
-              gfunc1[j][i][k][m]=((gfunc[j][i][k][m-2]-gfunc[j][i][k][m+2])+
-                  8.0*(gfunc[j][i][k][m+1]-gfunc[j][i][k][m+1]-gfunc[j][i][k][m-1]))/12.0;
-            }
-            for(m=0;m<ntheta-1;m++) {
-              gfunc2[j][i][k][m]=3.0*(gfunc[j][i][k][m+1]-gfunc[j][i][k][m])-
-                  2.0*gfunc1[j][i][k][m]-gfunc1[j][i][k][m+1];
-              gfunc3[j][i][k][m]=gfunc1[j][i][k][m]+gfunc1[j][i][k][m+1]-
-                  2.0*(gfunc[j][i][k][m+1]-gfunc[j][i][k][m]);
-            }
-            gfunc2[j][i][k][ntheta-1]=0.0;
-            gfunc3[j][i][k][ntheta-1]=0.0;
-            for(m=0;m<ntheta;m++) {
-              gfunc4[j][i][k][ntheta-1]=gfunc1[j][i][k][m]/dtheta;
-              gfunc5[j][i][k][ntheta-1]=2.0*gfunc2[j][i][k][m]/dtheta;
-              gfunc6[j][i][k][ntheta-1]=3.0*gfunc3[j][i][k][m]/dtheta;
             }
           }
         }
       }
+    } catch (TokenizerException & e) {
+      error->one(FLERR, e.what());
+    } catch (FileReaderException & fre) {
+      error->one(FLERR, fre.what());
     }
-    for(i=0;i<bop_types;i++) {
-      for(j=0;j<bop_types;j++) {
-        for(k=0;k<j;k++) {
-          if(npower<=2) {
-            for(n=0;n<ntheta;n++) {
-              gfunc[j][i][k][n]=gfunc[k][i][j][n];
-              gfunc1[j][i][k][n]=gfunc1[k][i][j][n];
-              gfunc2[j][i][k][n]=gfunc2[k][i][j][n];
-              gfunc3[j][i][k][n]=gfunc3[k][i][j][n];
-              gfunc4[j][i][k][n]=gfunc4[k][i][j][n];
-              gfunc5[j][i][k][n]=gfunc5[k][i][j][n];
-              gfunc6[j][i][k][n]=gfunc6[k][i][j][n];
-            }
-          } else {
-            for(n=0;n<npower+1;n++) {
-              gpara[j][i][k][n]=gpara[k][i][j][n];
-            }
-          }
-        }
-      }
-    }
-    fclose(fp);
   }
+
+  MPI_Bcast(&nr,1,MPI_INT,0,world);
+  MPI_Bcast(&nBOt,1,MPI_INT,0,world);
+  MPI_Bcast(&ntheta,1,MPI_INT,0,world);
+  MPI_Bcast(&bop_types,1,MPI_INT,0,world);
+  MPI_Bcast(&npairs,1,MPI_INT,0,world);
+  MPI_Bcast(&npower,1,MPI_INT,0,world);
+
+  if (comm->me != 0){
+    allocate_tables();
+    allocate();
+  }
+
   MPI_Bcast(&rdBO,1,MPI_DOUBLE,0,world);
   MPI_Bcast(&dBO,1,MPI_DOUBLE,0,world);
   MPI_Bcast(&rdtheta,1,MPI_DOUBLE,0,world);
diff --git a/src/MANYBODY/pair_bop.h b/src/MANYBODY/pair_bop.h
index 6c4b255d25..2d6962e27c 100644
--- a/src/MANYBODY/pair_bop.h
+++ b/src/MANYBODY/pair_bop.h
@@ -41,7 +41,6 @@ class PairBOP : public Pair {
   double memory_usage();
 
  private:
-  int me;
   int maxneigh;                 // maximum size of neighbor list on this processor
   int maxneigh3;                // maximum size of neighbor list on this processor
   int update_list;              // check for changing maximum size of neighbor list
@@ -207,6 +206,7 @@ class PairBOP : public Pair {
 
   void read_table(char *);
   void allocate();
+  void allocate_tables();
   void create_pi(int);
   void create_sigma(int);
   void destroy_pi();