diff --git a/lib/gpu/lal_amoeba.cu b/lib/gpu/lal_amoeba.cu
index 66926721cb..da5c6f0c3c 100644
--- a/lib/gpu/lal_amoeba.cu
+++ b/lib/gpu/lal_amoeba.cu
@@ -1922,6 +1922,185 @@ __kernel void k_amoeba_fphi_uind(const __global numtyp4 *restrict thetai1,
   }
 }
 
+
+/* ----------------------------------------------------------------------
+   fphi_mpole = multipole potential from grid
+   fphi_mpole extracts the permanent multipole potential from
+   the particle mesh Ewald grid
+------------------------------------------------------------------------- */
+
+__kernel void k_amoeba_fphi_mpole(const __global numtyp4 *restrict thetai1,
+                          const __global numtyp4 *restrict thetai2,
+                          const __global numtyp4 *restrict thetai3,
+                          const __global int *restrict igrid,
+                          const __global numtyp *restrict grid,
+                          __global numtyp *restrict fphi,
+                          const int bsorder, const int inum,
+                          const int nzlo_out, const int nylo_out,
+                          const int nxlo_out, const int ngridxy,
+                          const int ngridx)
+{
+  int tid=THREAD_ID_X;
+  int ii=tid+BLOCK_ID_X*BLOCK_SIZE_X;
+
+  if (ii<inum) {
+
+    int nlpts = (bsorder-1) / 2;
+    
+    int istart = fast_mul(ii,4);
+    int igridx = igrid[istart];
+    int igridy = igrid[istart+1];
+    int igridz = igrid[istart+2];
+    
+    // now istart is used to index thetai1, thetai2 and thetai3
+    istart = fast_mul(ii,bsorder);
+
+    // extract the permanent multipole field at each site
+
+    numtyp tuv000 = (numtyp)0.0;
+    numtyp tuv001 = (numtyp)0.0;
+    numtyp tuv010 = (numtyp)0.0;
+    numtyp tuv100 = (numtyp)0.0;
+    numtyp tuv200 = (numtyp)0.0;
+    numtyp tuv020 = (numtyp)0.0;
+    numtyp tuv002 = (numtyp)0.0;
+    numtyp tuv110 = (numtyp)0.0;
+    numtyp tuv101 = (numtyp)0.0;
+    numtyp tuv011 = (numtyp)0.0;
+    numtyp tuv300 = (numtyp)0.0;
+    numtyp tuv030 = (numtyp)0.0;
+    numtyp tuv003 = (numtyp)0.0;
+    numtyp tuv210 = (numtyp)0.0;
+    numtyp tuv201 = (numtyp)0.0;
+    numtyp tuv120 = (numtyp)0.0;
+    numtyp tuv021 = (numtyp)0.0;
+    numtyp tuv102 = (numtyp)0.0;
+    numtyp tuv012 = (numtyp)0.0;
+    numtyp tuv111 = (numtyp)0.0;
+
+    int k = (igridz - nzlo_out) - nlpts;
+    for (int kb = 0; kb < bsorder; kb++) {
+      /*
+      v0 = thetai3[m][kb][0];
+      v1 = thetai3[m][kb][1];
+      v2 = thetai3[m][kb][2];
+      v3 = thetai3[m][kb][3];
+      */
+      int i3 = istart + kb;
+      numtyp4 tha3 = thetai3[i3];
+      numtyp v0 = tha3.x;
+      numtyp v1 = tha3.y;
+      numtyp v2 = tha3.z;
+      numtyp v3 = tha3.w;
+      numtyp tu00 = (numtyp)0.0;
+      numtyp tu10 = (numtyp)0.0;
+      numtyp tu01 = (numtyp)0.0;
+      numtyp tu20 = (numtyp)0.0;
+      numtyp tu11 = (numtyp)0.0;
+      numtyp tu02 = (numtyp)0.0;
+      numtyp tu30 = (numtyp)0.0;
+      numtyp tu21 = (numtyp)0.0;
+      numtyp tu12 = (numtyp)0.0;
+      numtyp tu03 = (numtyp)0.0;
+
+      int j = (igridy - nylo_out) - nlpts;
+      for (int jb = 0; jb < bsorder; jb++) {
+        /*
+        u0 = thetai2[m][jb][0];
+        u1 = thetai2[m][jb][1];
+        u2 = thetai2[m][jb][2];
+        u3 = thetai2[m][jb][3];
+        */
+        int i2 = istart + jb;
+        numtyp4 tha2 = thetai2[i2];
+        numtyp u0 = tha2.x;
+        numtyp u1 = tha2.y;
+        numtyp u2 = tha2.z;
+        numtyp u3 = tha2.w;
+        numtyp t0 = (numtyp)0.0;
+        numtyp t1 = (numtyp)0.0;
+        numtyp t2 = (numtyp)0.0;
+        numtyp t3 = (numtyp)0.0;
+
+        int i = (igridx - nxlo_out) - nlpts;
+        for (int ib = 0; ib < bsorder; ib++) {
+          int i1 = istart + ib;
+          numtyp4 tha1 = thetai1[i1];
+          int gidx = 2*(k*ngridxy + j*ngridx + i);
+          numtyp tq = grid[gidx];
+          t0 += tq*tha1.x;
+          t1 += tq*tha1.y;
+          t2 += tq*tha1.z;
+          t3 += tq*tha1.w;
+          i++;
+        }
+
+        tu00 += t0*u0;
+        tu10 += t1*u0;
+        tu01 += t0*u1;
+        tu20 += t2*u0;
+        tu11 += t1*u1;
+        tu02 += t0*u2;
+        tu30 += t3*u0;
+        tu21 += t2*u1;
+        tu12 += t1*u2;
+        tu03 += t0*u3;
+        j++;
+      }
+
+      tuv000 += tu00*v0;
+      tuv100 += tu10*v0;
+      tuv010 += tu01*v0;
+      tuv001 += tu00*v1;
+      tuv200 += tu20*v0;
+      tuv020 += tu02*v0;
+      tuv002 += tu00*v2;
+      tuv110 += tu11*v0;
+      tuv101 += tu10*v1;
+      tuv011 += tu01*v1;
+      tuv300 += tu30*v0;
+      tuv030 += tu03*v0;
+      tuv003 += tu00*v3;
+      tuv210 += tu21*v0;
+      tuv201 += tu20*v1;
+      tuv120 += tu12*v0;
+      tuv021 += tu02*v1;
+      tuv102 += tu10*v2;
+      tuv012 += tu01*v2;
+      tuv111 += tu11*v1;
+      k++;
+    }
+
+    numtyp buf[20];
+    buf[0] = tuv000;
+    buf[1] = tuv100;
+    buf[2] = tuv010;
+    buf[3] = tuv001;
+    buf[4] = tuv200;
+    buf[5] = tuv020;
+    buf[6] = tuv002;
+    buf[7] = tuv110;
+    buf[8] = tuv101;
+    buf[9] = tuv011;
+    buf[10] = tuv300;
+    buf[11] = tuv030;
+    buf[12] = tuv003;
+    buf[13] = tuv210;
+    buf[14] = tuv201;
+    buf[15] = tuv120;
+    buf[16] = tuv021;
+    buf[17] = tuv102;
+    buf[18] = tuv012;
+    buf[19] = tuv111;
+
+    int idx = ii;    
+    for (int m = 0; m < 20; m++) {
+      fphi[idx] = buf[m];
+      idx += inum;
+    }
+  }
+}
+
 /* ----------------------------------------------------------------------
    scan standard neighbor list and make it compatible with 1-5 neighbors
    if IJ entry is a 1-2,1-3,1-4 neighbor then adjust offset to SBBITS15
diff --git a/src/AMOEBA/amoeba_kspace.cpp b/src/AMOEBA/amoeba_kspace.cpp
index 6f76c7ec51..0f74f19797 100644
--- a/src/AMOEBA/amoeba_kspace.cpp
+++ b/src/AMOEBA/amoeba_kspace.cpp
@@ -62,25 +62,31 @@ void PairAmoeba::moduli()
   int maxfft = MAX(nfft1,nfft2);
   maxfft = MAX(maxfft,nfft3);
 
-  double *array = new double[bsorder];
-  double *bsarray = new double[maxfft];
-
+  //double *array = new double[bsorder];
+  //double *bsarray = new double[maxfft];
+  if (maxfft > _nfft_max) {
+    memory->destroy(_moduli_bsarray);
+    _nfft_max = maxfft;
+    memory->create(_moduli_bsarray,_nfft_max,"amoeba:_moduli_bsarray");
+  }
+  
   // compute and load the moduli values
 
   double x = 0.0;
-  bspline(x,bsorder,array);
+  //bspline(x,bsorder,array);
+  bspline(x,bsorder,_moduli_array);
 
-  for (i = 0; i < maxfft; i++) bsarray[i] = 0.0;
-  for (i = 0; i < bsorder; i++) bsarray[i+1] = array[i];
+  for (i = 0; i < maxfft; i++) _moduli_bsarray[i] = 0.0;
+  for (i = 0; i < bsorder; i++) _moduli_bsarray[i+1] = _moduli_array[i];
 
-  dftmod(bsmod1,bsarray,nfft1,bsorder);
-  dftmod(bsmod2,bsarray,nfft2,bsorder);
-  dftmod(bsmod3,bsarray,nfft3,bsorder);
+  dftmod(bsmod1,_moduli_bsarray,nfft1,bsorder);
+  dftmod(bsmod2,_moduli_bsarray,nfft2,bsorder);
+  dftmod(bsmod3,_moduli_bsarray,nfft3,bsorder);
 
   // perform deallocation of local arrays
 
-  delete[] array;
-  delete[] bsarray;
+  //delete[] array;
+  //delete[] bsarray;
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp
index 9890904e42..d301a86cdb 100644
--- a/src/AMOEBA/pair_amoeba.cpp
+++ b/src/AMOEBA/pair_amoeba.cpp
@@ -68,67 +68,71 @@ PairAmoeba::PairAmoeba(LAMMPS *lmp) : Pair(lmp)
   // force field settings
 
   nmax = 0;
-  xaxis2local = yaxis2local = zaxis2local = NULL;
-  rpole = NULL;
-  tq = NULL;
+  xaxis2local = yaxis2local = zaxis2local = nullptr;
+  rpole = nullptr;
+  tq = nullptr;
 
-  red2local = NULL;
-  xred = NULL;
+  red2local = nullptr;
+  xred = nullptr;
 
-  uind = uinp = udirp = NULL;
-  uopt = uoptp = NULL;
-  fopt = foptp = NULL;
-  field = fieldp = NULL;
-  ufld = dufld = NULL;
-  rsd = rsdp = NULL;
-  zrsd = zrsdp = NULL;
+  uind = uinp = udirp = nullptr;
+  uopt = uoptp = nullptr;
+  fopt = foptp = nullptr;
+  field = fieldp = nullptr;
+  ufld = dufld = nullptr;
+  rsd = rsdp = nullptr;
+  zrsd = zrsdp = nullptr;
 
-  cmp = fmp = NULL;
-  cphi = fphi = NULL;
+  cmp = fmp = nullptr;
+  cphi = fphi = nullptr;
 
-  poli = NULL;
-  conj = conjp = NULL;
-  vec = vecp = NULL;
-  udir = usum = usump = NULL;
+  _moduli_array = nullptr;
+  _moduli_bsarray = nullptr;
+  _nfft_max = 0;
 
-  fuind = fuinp = NULL;
-  fdip_phi1 = fdip_phi2 = fdip_sum_phi = NULL;
-  dipfield1 = dipfield2 = NULL;
+  poli = nullptr;
+  conj = conjp = nullptr;
+  vec = vecp = nullptr;
+  udir = usum = usump = nullptr;
 
-  fphid = fphip = NULL;
-  fphidp = cphidp = NULL;
+  fuind = fuinp = nullptr;
+  fdip_phi1 = fdip_phi2 = fdip_sum_phi = nullptr;
+  dipfield1 = dipfield2 = nullptr;
+
+  fphid = fphip = nullptr;
+  fphidp = cphidp = nullptr;
 
   bsordermax = 0;
-  thetai1 = thetai2 = thetai3 = NULL;
-  bsmod1 = bsmod2 = bsmod3 = NULL;
-  bsbuild = NULL;
-  igrid = NULL;
-  m_kspace = p_kspace = pc_kspace = d_kspace = NULL;
-  i_kspace = ic_kspace = NULL;
+  thetai1 = thetai2 = thetai3 = nullptr;
+  bsmod1 = bsmod2 = bsmod3 = nullptr;
+  bsbuild = nullptr;
+  igrid = nullptr;
+  m_kspace = p_kspace = pc_kspace = d_kspace = nullptr;
+  i_kspace = ic_kspace = nullptr;
 
-  numneigh_dipole = NULL;
-  firstneigh_dipole = NULL;
-  firstneigh_dipdip = NULL;
-  ipage_dipole = NULL;
-  dpage_dipdip = NULL;
+  numneigh_dipole = nullptr;
+  firstneigh_dipole = nullptr;
+  firstneigh_dipdip = nullptr;
+  ipage_dipole = nullptr;
+  dpage_dipdip = nullptr;
 
-  numneigh_precond = NULL;
-  firstneigh_precond = NULL;
-  ipage_precond = NULL;
+  numneigh_precond = nullptr;
+  firstneigh_precond = nullptr;
+  ipage_precond = nullptr;
 
-  firstneigh_pcpc = NULL;
-  dpage_pcpc = NULL;
+  firstneigh_pcpc = nullptr;
+  dpage_pcpc = nullptr;
 
-  qfac = NULL;
-  gridfft1 = NULL;
+  qfac = nullptr;
+  gridfft1 = nullptr;
 
   initialize_type_class();
   initialize_vdwl();
   initialize_smallsize();
 
-  forcefield = NULL;
+  forcefield = nullptr;
 
-  id_pole = id_udalt = id_upalt = NULL;
+  id_pole = id_udalt = id_upalt = nullptr;
 
   nualt = 0;
   first_flag = 1;
@@ -220,6 +224,9 @@ PairAmoeba::~PairAmoeba()
   memory->destroy(fphidp);
   memory->destroy(cphidp);
 
+  memory->destroy(_moduli_array);
+  memory->destroy(_moduli_bsarray);
+
   memory->destroy(thetai1);
   memory->destroy(thetai2);
   memory->destroy(thetai3);
@@ -2312,6 +2319,8 @@ void PairAmoeba::grow_local()
     firstneigh_pcpc = (double **)
       memory->smalloc(nmax*sizeof(double *),"induce:firstneigh_pcpc");
   }
+
+  memory->create(_moduli_array,bsordermax,"amoeba:_moduli_array");
 }
 
 /* ----------------------------------------------------------------------
diff --git a/src/AMOEBA/pair_amoeba.h b/src/AMOEBA/pair_amoeba.h
index 24ce6fcfbc..91ec8faf0c 100644
--- a/src/AMOEBA/pair_amoeba.h
+++ b/src/AMOEBA/pair_amoeba.h
@@ -337,7 +337,11 @@ class PairAmoeba : public Pair {
   double *gridfft1;                // copy of p_kspace FFT grid
 
   double **cmp,**fmp;              // Cartesian and fractional multipoles
-  double **cphi,**fphi;     
+  double **cphi,**fphi;
+
+  double *_moduli_array;           // buffers for moduli
+  double *_moduli_bsarray;
+  int _nfft_max;
 
   // params for current KSpace solve and FFT being worked on
 
@@ -347,8 +351,12 @@ class PairAmoeba : public Pair {
   double ctf[10][10];      // indices NOT flipped vs Fortran
   double ftc[10][10];      // indices NOT flipped vs Fortran
 
-  class AmoebaConvolution *m_kspace,*p_kspace,*pc_kspace,*d_kspace;
-  class AmoebaConvolution *i_kspace,*ic_kspace;
+  class AmoebaConvolution *m_kspace;   // multipole KSpace
+  class AmoebaConvolution *p_kspace;   // polar KSpace
+  class AmoebaConvolution *pc_kspace;  
+  class AmoebaConvolution *d_kspace;   // dispersion KSpace 
+  class AmoebaConvolution *i_kspace;   // induce KSpace
+  class AmoebaConvolution *ic_kspace;
 
   // FFT grid size factors