diff --git a/src/MANYBODY/pair_comb3.cpp b/src/MANYBODY/pair_comb3.cpp
index 30a0e36c2d..7ba6605478 100644
--- a/src/MANYBODY/pair_comb3.cpp
+++ b/src/MANYBODY/pair_comb3.cpp
@@ -27,6 +27,7 @@
 #include "group.h"
 #include "math_const.h"
 #include "math_extra.h"
+#include "math_special.h"
 #include "memory.h"
 #include "my_page.h"
 #include "neigh_list.h"
@@ -42,6 +43,7 @@
 using namespace LAMMPS_NS;
 using namespace MathConst;
 using namespace MathExtra;
+using namespace MathSpecial;
 
 #define DELTA 4
 #define PGDELTA 1
@@ -666,7 +668,7 @@ void PairComb3::setup_params()
     params[m].Qo = (params[m].QU+params[m].QL)/2.0; // (A22)
     params[m].dQ = (params[m].QU-params[m].QL)/2.0; // (A21)
     params[m].aB = 1.0 /
-      (1.0-pow(fabs(params[m].Qo/params[m].dQ),10)); // (A20)
+      (1.0-powint(fabs(params[m].Qo/params[m].dQ),10)); // (A20)
     params[m].bB = pow(fabs(params[m].aB),0.1)/params[m].dQ; // (A19)
     params[m].nD = log(params[m].DU/(params[m].DU-params[m].DL))/
                     log(params[m].QU/(params[m].QU-params[m].QL));
@@ -1432,7 +1434,7 @@ void PairComb3::repulsive(Param *parami, Param *paramj, double rsq,
 
   vrcs = 1.0; fvrcs = 0.0;
   if (romi != 0.0 && r < addr) {
-    vrcs += romi * pow((1.0-r/addr),2.0);
+    vrcs += romi * square(1.0-r/addr);
     fvrcs = romi * 2.0 * (r/addr-1.0)/addr;
     fforce = fforce*vrcs - caj * tmp_fc * vrcs * fvrcs;
   }
@@ -1455,7 +1457,7 @@ double PairComb3::zeta(Param *parami, Param *paramj, double rsqij,
   costheta = dot3(delrij,delrik) / (rij*rik);
 
   rlm3 = parami->beta;
-  arg = pow(rlm3*(rij-rik),int(parami->powermint));
+  arg = powint(rlm3*(rij-rik),int(parami->powermint));
   if (arg > 69.0776) ex_delr = 1.e30;
   else if (arg < -69.0776) ex_delr = 0.0;
   else ex_delr = exp(arg);
@@ -1776,8 +1778,8 @@ double PairComb3::self(Param *param, double qi)
 
  self_tmp = qi*(s1+qi*(s2+qi*(s3+qi*s4)));
 
- if (qi < qmin) self_tmp += cmin * pow((qi-qmin),4);
- if (qi > qmax) self_tmp += cmax * pow((qi-qmax),4);
+ if (qi < qmin) self_tmp += cmin * powint((qi-qmin),4);
+ if (qi > qmax) self_tmp += cmax * powint((qi-qmax),4);
 
  return self_tmp;
 }
@@ -1802,8 +1804,8 @@ void PairComb3::comb_fa(double r, Param *parami, Param *paramj, double iq,
   Di = Dj = Bsi = 0.0;
   Di = parami->DU + pow(fabs(parami->bD*(parami->QU-qi)),parami->nD);
   Dj = paramj->DU + pow(fabs(paramj->bD*(paramj->QU-qj)),paramj->nD);
-  YYBn = (parami->aB-fabs(pow(parami->bB*(qi-parami->Qo),10)));
-  YYBj = (paramj->aB-fabs(pow(paramj->bB*(qj-paramj->Qo),10)));
+  YYBn = (parami->aB-fabs(powint(parami->bB*(qi-parami->Qo),10)));
+  YYBj = (paramj->aB-fabs(powint(paramj->bB*(qj-paramj->Qo),10)));
 
   if (YYBn*YYBj > 0.0) {
     AlfDiAlfDj = exp(0.5*(parami->alfi*Di+paramj->alfi*Dj));
@@ -2124,7 +2126,7 @@ void PairComb3::comb_zetaterm_d(double prefac_ij1, double prefac_ij2,
   fc_k = comb_fc(rik,paramk);
   dfc_k = comb_fc_d(rik,paramk);
   rlm3 = parami->beta;
-  tmp = pow(rlm3*(rij-rik),mint);
+  tmp = powint(rlm3*(rij-rik),mint);
 
   if (tmp > 69.0776) ex_delr = 1.e30;
   else if (tmp < -69.0776) ex_delr = 0.0;
@@ -2155,7 +2157,7 @@ void PairComb3::comb_zetaterm_d(double prefac_ij1, double prefac_ij2,
     com3k = 0.0;
   }
 
-  ex_delr_d = mint*pow(rlm3,mint)*pow((rij-rik),(mint-1))*ex_delr; // com3
+  ex_delr_d = mint*powint(rlm3,mint)*powint((rij-rik),(mint-1))*ex_delr; // com3
   scale3(-dfc_k*gijk*ex_delr,rik_hat,dri);          // com1
   scaleadd3(fc_k*gijk_d*ex_delr,dcosdri,dri,dri);   // com2
   scaleadd3(fc_k*gijk*ex_delr_d,rik_hat,dri,dri);   // com3 cont'd
@@ -2690,10 +2692,10 @@ void PairComb3::field(Param *parami, Param *paramj, double rsq, double iq,
   pcmi1 = parami->pcmn1;
   pcmi2 = parami->pcmn2;
 
-  rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
-  rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
-  rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
-  rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));
+  rf3i = r3/(square(r3)+cube(pcmi1));
+  rcf3i = rc3/(square(rc3)+cube(pcmi1));
+  rf5i = r5/(square(r5)+powint(pcmi2,5));
+  rcf5i = rc5/(square(rc5)+powint(pcmi2,5));
 
   drf3i = 3/r*rf3i-6*rsq*rf3i*rf3i;
   drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
@@ -2742,7 +2744,7 @@ void PairComb3::rad_calc(double r, Param *parami, Param *paramj,
   vrad = pradx = prady = pradz = 0.0;
   xrad = -comb_fc(r,parami)*parami->pcross + xcn;
   yrad = -comb_fc(r,paramj)*paramj->pcross + ycn;
-  zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
+  zcon = 1.0 + square(kconjug) + square(lconjug);
 
   if (xrad < 0.0) xrad = 0.0;
   if (yrad < 0.0) yrad = 0.0;
@@ -2893,14 +2895,14 @@ double PairComb3::bbtor1(int torindx, Param *paramk, Param *paraml,
     tork[2] = delrk[0]*delrj[1] - delrk[1]*delrj[0];
     torl[2] = delrj[0]*delrl[1] - delrj[1]*delrl[0];
     TT2 = dot3(tork,torl);
-    rmut = pow((TT2/TT1),2);
+    rmut = square(TT2/TT1);
     if (torindx>=1) {
       btt = 1.0 - rmut;
       return btt * fc1k * fc1l;
     }
     else {
       btt=paramk->ptork1-TT2/TT1;
-      btt=paramk->ptork2*pow(btt,2);
+      btt=paramk->ptork2*square(btt);
       return btt * fc1k * fc1l;
     }
 
@@ -2930,7 +2932,7 @@ void PairComb3::tor_calc(double r, Param *parami, Param *paramj,
   } else {
     xtor = -comb_fc(r, parami) * parami->pcross + xcn;
     ytor = -comb_fc(r, paramj) * paramj->pcross + ycn;
-    zcon = 1.0 + pow(kconjug,2) + pow(lconjug,2);
+    zcon = 1.0 + square(kconjug) + square(lconjug);
     if (xtor < 0.0) xtor = 0.0;
     if (ytor < 0.0) ytor = 0.0;
     if (zcon < 1.0) zcon = 1.0;
@@ -2981,9 +2983,9 @@ void PairComb3::tor_int(int torindx,double xtor, double ytor, double zcon, int l
           * pow(ytor,iin3[j][1]) * pow(zcon,iin3[j][2]);
     vtor += x;
 
-  if (xtor > 1.0e-8 ) dtorx += x*iin3[j][0]/xtor;
-  if (ytor > 1.0e-8 ) dtory += x*iin3[j][1]/ytor;
-  if (zcon > 1.0e-8 ) dtorz += x*iin3[j][2]/zcon;
+  if (xtor > 1.0e-8) dtorx += x*iin3[j][0]/xtor;
+  if (ytor > 1.0e-8) dtory += x*iin3[j][1]/ytor;
+  if (zcon > 1.0e-8) dtorz += x*iin3[j][2]/zcon;
   }
 }
 
@@ -3022,10 +3024,10 @@ void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
     fcp1k = comb_fc_d(rik,paramk);
     fc1l  = comb_fc(rjl,paraml);
     fcp1l = comb_fc_d(rjl,paraml);
-    srmul2 = pow(srmul,2);
+    srmul2 = square(srmul);
 
     TT1 = rik*rjl*rij*rij*srmu*srmul;
-    dt1dik = -rmu/pow(srmu,2);
+    dt1dik = -rmu/square(srmu);
     dt1djl = -rmul/srmul2;
     tork[0] = delrk[1]*delrj[2] - delrk[2]*delrj[1];
     torl[0] = delrj[1]*delrl[2] - delrj[2]*delrl[1];
@@ -3051,12 +3053,12 @@ void PairComb3::tor_force(int torindx, Param *paramk, Param *paraml,
     rmut = TT2/TT1;
 
     if (torindx>=1) {
-        btt = 1.0 - pow(rmut,2);
+        btt = 1.0 - square(rmut);
         AA = -2.0 * ptorr * rmut * fc1k * fc1l / TT1;
     }
     else {
         btt=paramk->ptork1-rmut;
-        btt=paramk->ptork2*pow(btt,2);
+        btt=paramk->ptork2*square(btt);
         AA = -2.0 * ptorr * paramk->ptork2 *
           (paramk->ptork1-rmut) * fc1k * fc1l /TT1;
    }
@@ -3300,8 +3302,8 @@ double PairComb3::qfo_self(Param *param, double qi)
   cmin = cmax = 100.0;
   self_d = s1+qi*(2.0*s2+qi*(3.0*s3+qi*4.0*s4));
 
-  if (qi < qmin) self_d += 4.0 * cmin * pow((qi-qmin),3);
-  if (qi > qmax) self_d += 4.0 * cmax * pow((qi-qmax),3);
+  if (qi < qmin) self_d += 4.0 * cmin * cube((qi-qmin));
+  if (qi > qmax) self_d += 4.0 * cmax * cube((qi-qmax));
 
   return self_d;
 }
@@ -3392,10 +3394,10 @@ void PairComb3::qfo_field(Param *parami, Param *paramj, double rsq,
   pcmi1 = parami->pcmn1;
   pcmi2 = parami->pcmn2;
 
-  rf3i = r3/(pow(r3,2)+pow(pcmi1,3));
-  rcf3i = rc3/(pow(rc3,2)+pow(pcmi1,3));
-  rf5i = r5/(pow(r5,2)+pow(pcmi2,5));
-  rcf5i = rc5/(pow(rc5,2)+pow(pcmi2,5));
+  rf3i = r3/(square(r3)+cube(pcmi1));
+  rcf3i = rc3/(square(rc3)+cube(pcmi1));
+  rf5i = r5/(square(r5)+powint(pcmi2,5));
+  rcf5i = rc5/(square(rc5)+powint(pcmi2,5));
 
   drcf3i = 3/rc*rcf3i-6*rc2*rcf3i*rcf3i;
   drcf5i = 5/rc*rcf5i-10*rc4*rcf5i*rcf5i;
@@ -3477,7 +3479,7 @@ void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
     Di = parami->DU + pow(QUchi,parami->nD);                            // YYDin
     dDi = -parami->nD * parami->bD * pow(QUchi,(parami->nD-1.0));       // YYDiqp
     Bsi = parami->aB - pow(QOchi,10);                                   // YYBsin
-    dBsi = -parami->bB * 10.0 * pow(QOchi,9.0);                         // YYBsiqp
+    dBsi = -parami->bB * 10.0 * powint(QOchi,9);                        // YYBsiqp
   }
 
   if (jq < paramj->QL-0.2) {
@@ -3492,7 +3494,7 @@ void PairComb3::qfo_short(Param *parami, Param *paramj, double rsq,
     Dj = paramj->DU + pow(QUchj,paramj->nD);                            // YYDij
     dDj = -paramj->nD * paramj->bD * pow(QUchj,(paramj->nD-1.0));       // YYDiqpj
     Bsj = paramj->aB - pow(QOchj,10);                                   // YYBsij
-    dBsj = -paramj->bB * 10.0 * pow(QOchj,9.0);                         // YYBsiqpj
+    dBsj = -paramj->bB * 10.0 * powint(QOchj,9);                        // YYBsiqpj
   }
 
   LamDiLamDj = exp(0.5*(parami->lami*Di+paramj->lami*Dj)-rlm1*r);