diff --git a/src/AMOEBA/amoeba_hal.cpp b/src/AMOEBA/amoeba_hal.cpp
index 21bb7ad099..896595945e 100644
--- a/src/AMOEBA/amoeba_hal.cpp
+++ b/src/AMOEBA/amoeba_hal.cpp
@@ -16,12 +16,17 @@
 
 #include "atom.h"
 #include "error.h"
+#include "math_special.h"
 #include "neigh_list.h"
 
 #include <cmath>
 
 using namespace LAMMPS_NS;
 
+using MathSpecial::square;
+using MathSpecial::cube;
+using MathSpecial::powint;
+
 enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
 
 /* ----------------------------------------------------------------------
@@ -114,14 +119,14 @@ void PairAmoeba::hal()
       }
       eps *= factor_hal;
 
-      rv7 = pow(rv,7.0);
-      rik6 = pow(rik2,3.0);
+      rv7 = powint(rv,7);
+      rik6 = cube(rik2);
       rik7 = rik6 * rik;
       rho = rik7 + ghal*rv7;
       tau = (dhal+1.0) / (rik + dhal*rv);
-      tau7 = pow(tau,7.0);
+      tau7 = powint(tau,7);
       dtau = tau / (dhal+1.0);
-      gtau = eps*tau7*rik6*(ghal+1.0)*pow(rv7/rho,2.0);
+      gtau = eps*tau7*rik6*(ghal+1.0)*square(rv7/rho);
       e = eps*tau7*rv7*((ghal+1.0)*rv7/rho-2.0);
       de = -7.0 * (dtau*e+gtau);
 
diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp
index 3d9d7809cc..43688cef94 100644
--- a/src/AMOEBA/amoeba_induce.cpp
+++ b/src/AMOEBA/amoeba_induce.cpp
@@ -22,6 +22,7 @@
 #include "fft3d_wrap.h"
 #include "fix_store.h"
 #include "math_const.h"
+#include "math_special.h"
 #include "memory.h"
 #include "my_page.h"
 #include "neigh_list.h"
@@ -32,6 +33,8 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
+using MathSpecial::cube;
+
 enum{INDUCE,RSD,SETUP_AMOEBA,SETUP_HIPPO,KMPOLE,AMGROUP};   // forward comm
 enum{FIELD,ZRSD,TORQUE,UFLD};                               // reverse comm
 enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
@@ -732,7 +735,7 @@ void PairAmoeba::uscale0b(int mode, double **rsd, double **rsdp,
         damp = pdi * pdamp[jtype];
         if (damp != 0.0) {
           pgamma = MIN(pti,thole[jtype]);
-          damp = -pgamma * pow((r/damp),3.0);
+          damp = -pgamma * cube(r/damp);
           if (damp > -50.0) {
             expdamp = exp(damp);
             scale3 *= 1.0 - expdamp;
@@ -1332,7 +1335,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
             }
           } else {
             pgamma = MIN(pti,thole[jtype]);
-            damp = pgamma * pow(r/damp,3.0);
+            damp = pgamma * cube(r/damp);
             if (damp < 50.0) {
               expdamp = exp(-damp);
               scale3 = 1.0 - expdamp;
@@ -1384,7 +1387,7 @@ void PairAmoeba::udirect2b(double **field, double **fieldp)
           damp = pdi * pdamp[jtype];
           if (damp != 0.0) {
             pgamma = MIN(pti,thole[jtype]);
-            damp = pgamma * pow(r/damp,3.0);
+            damp = pgamma * cube(r/damp);
             if (damp < 50.0) {
               expdamp = exp(-damp);
               scale3 = 1.0 - expdamp;
diff --git a/src/AMOEBA/amoeba_kspace.cpp b/src/AMOEBA/amoeba_kspace.cpp
index 5e12deae1c..51b206b6d8 100644
--- a/src/AMOEBA/amoeba_kspace.cpp
+++ b/src/AMOEBA/amoeba_kspace.cpp
@@ -17,6 +17,7 @@
 #include "atom.h"
 #include "domain.h"
 #include "math_const.h"
+#include "math_special.h"
 #include "memory.h"
 
 #include <cmath>
@@ -24,6 +25,8 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
+using MathSpecial::powint;
+
 #define ANINT(x) ((x)>0 ? floor((x)+0.5) : ceil((x)-0.5))
 
 /* ----------------------------------------------------------------------
@@ -173,13 +176,13 @@ void PairAmoeba::dftmod(double *bsmod, double *bsarray, int nfft, int order)
       factor = MY_PI * k / nfft;
       for (j = 1; j <= jcut; j++) {
         arg = factor / (factor + MY_PI*j);
-        sum1 += pow(arg,order);
-        sum2 += pow(arg,order2);
+        sum1 += powint(arg,order);
+        sum2 += powint(arg,order2);
       }
       for (j = 1; j <= jcut; j++) {
         arg = factor / (factor - MY_PI*j);
-        sum1 += pow(arg,order);
-        sum2 += pow(arg,order2);
+        sum1 += powint(arg,order);
+        sum2 += powint(arg,order2);
       }
       zeta = sum2 / sum1;
     }
diff --git a/src/AMOEBA/amoeba_multipole.cpp b/src/AMOEBA/amoeba_multipole.cpp
index 5d11bde1ab..8466a8fe1d 100644
--- a/src/AMOEBA/amoeba_multipole.cpp
+++ b/src/AMOEBA/amoeba_multipole.cpp
@@ -20,6 +20,7 @@
 #include "domain.h"
 #include "fft3d_wrap.h"
 #include "math_const.h"
+#include "math_special.h"
 #include "memory.h"
 #include "neigh_list.h"
 
@@ -29,6 +30,8 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
+using MathSpecial::square;
+
 enum{FIELD,ZRSD,TORQUE,UFLD};                          // reverse comm
 enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
 
@@ -670,7 +673,7 @@ void PairAmoeba::multipole_kspace()
   nzlo = m_kspace->nzlo_fft;
   nzhi = m_kspace->nzhi_fft;
 
-  pterm = pow((MY_PI/aewald),2.0);
+  pterm = square(MY_PI/aewald);
   volterm = MY_PI * volbox;
 
   n = 0;
diff --git a/src/AMOEBA/amoeba_polar.cpp b/src/AMOEBA/amoeba_polar.cpp
index b8020cb513..ad6b585f25 100644
--- a/src/AMOEBA/amoeba_polar.cpp
+++ b/src/AMOEBA/amoeba_polar.cpp
@@ -20,6 +20,7 @@
 #include "domain.h"
 #include "fft3d_wrap.h"
 #include "math_const.h"
+#include "math_special.h"
 #include "memory.h"
 #include "neigh_list.h"
 
@@ -29,6 +30,9 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 
+using MathSpecial::square;
+using MathSpecial::cube;
+
 enum{FIELD,ZRSD,TORQUE,UFLD};                          // reverse comm
 enum{MUTUAL,OPT,TCG,DIRECT};
 enum{VDWL,REPULSE,QFER,DISP,MPOLE,POLAR,USOLV,DISP_LONG,MPOLE_LONG,POLAR_LONG};
@@ -82,7 +86,7 @@ void PairAmoeba::polar()
 
   // compute the Ewald self-energy torque and virial terms
 
-  term = (4.0/3.0) * felec * pow(aewald,3.0) / MY_PIS;
+  term = (4.0/3.0) * felec * cube(aewald) / MY_PIS;
 
   for (i = 0; i < nlocal; i++) {
     dix = rpole[i][1];
@@ -454,7 +458,7 @@ void PairAmoeba::polar_real()
         damp = pdi * pdamp[jtype];
         if (damp != 0.0) {
           pgamma = MIN(pti,thole[jtype]);
-          damp = pgamma * pow(r/damp,3.0);
+          damp = pgamma * cube(r/damp);
           if (damp < 50.0) {
             expdamp = exp(-damp);
             sc3 = 1.0 - expdamp;
@@ -1272,7 +1276,7 @@ void PairAmoeba::polar_kspace()
   int nlocal = atom->nlocal;
 
   double volbox = domain->prd[0] * domain->prd[1] * domain->prd[2];
-  pterm = pow((MY_PI/aewald),2.0);
+  pterm = square(MY_PI/aewald);
   volterm = MY_PI * volbox;
 
   // initialize variables required for the scalar summation
diff --git a/src/AMOEBA/amoeba_repulsion.cpp b/src/AMOEBA/amoeba_repulsion.cpp
index 0784a32d0b..041d74e54d 100644
--- a/src/AMOEBA/amoeba_repulsion.cpp
+++ b/src/AMOEBA/amoeba_repulsion.cpp
@@ -504,7 +504,7 @@ void PairAmoeba::damprep(double r, double r2, double rr1, double rr3,
     dmpk24 = dmpk23 * dmpk2;
     dmpk25 = dmpk24 * dmpk2;
     term = dmpi22 - dmpk22;
-    pre = 8192.0 * dmpi23 * dmpk23 / pow(term,4.0);
+    pre = 8192.0 * dmpi23 * dmpk23 / (term*term*term*term);
     tmp = 4.0 * dmpi2 * dmpk2 / term;
     s = (dampi-tmp)*expk + (dampk+tmp)*expi;
 
diff --git a/src/AMOEBA/pair_amoeba.cpp b/src/AMOEBA/pair_amoeba.cpp
index be5f9c73df..e8b7b71753 100644
--- a/src/AMOEBA/pair_amoeba.cpp
+++ b/src/AMOEBA/pair_amoeba.cpp
@@ -25,6 +25,7 @@
 #include "force.h"
 #include "gridcomm.h"
 #include "group.h"
+#include "math_special.h"
 #include "memory.h"
 #include "modify.h"
 #include "my_page.h"
@@ -40,6 +41,8 @@
 
 using namespace LAMMPS_NS;
 
+using MathSpecial::powint;
+
 enum{INDUCE,RSD,SETUP_AMOEBA,SETUP_HIPPO,KMPOLE,AMGROUP,PVAL};  // forward comm
 enum{FIELD,ZRSD,TORQUE,UFLD};                                   // reverse comm
 enum{ARITHMETIC,GEOMETRIC,CUBIC_MEAN,R_MIN,SIGMA,DIAMETER,HARMONIC,HHG,W_H};
@@ -1956,7 +1959,7 @@ void PairAmoeba::choose(int which)
 
   // taper coeffs
 
-  double denom = pow(off-cut,5.0);
+  double denom = powint(off-cut,5);
   c0 = off*off2 * (off2 - 5.0*off*cut + 10.0*cut2) / denom;
   c1 = -30.0 * off2*cut2 / denom;
   c2 = 30.0 * (off2*cut+off*cut2) / denom;
@@ -2026,7 +2029,7 @@ void PairAmoeba::mix()
       } else if (epsilon_rule == HHG) {
         eij = 4.0 * (ei*ej) / ((sei+sej)*(sei+sej));
       } else if (epsilon_rule == W_H) {
-        eij = 2.0 * (sei*sej) * pow(ri*rj,3.0) / (pow(ri,6.0) + pow(rj,6.0));
+        eij = 2.0 * (sei*sej) * powint(ri*rj,3) / (powint(ri,6) + powint(rj,6));
       } else {
         eij = sei * sej;
       }