diff --git a/src/DIPOLE/pair_lj_long_dipole_long.cpp b/src/DIPOLE/pair_lj_long_dipole_long.cpp
index 1635bd4f636f33d3b4f678aef45a4d2b4727fe2f..21fc035854d7a2fabeba3b6fdf54eb24055c4aab 100644
--- a/src/DIPOLE/pair_lj_long_dipole_long.cpp
+++ b/src/DIPOLE/pair_lj_long_dipole_long.cpp
@@ -452,7 +452,7 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
       ni = sbmask(j);                                   // special index
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                           // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -463,9 +463,9 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
       if (order3 && (rsq < cut_coulsq)) {               // dipole
         memcpy(muj, jmu = mu0+(j<<2), sizeof(vector));
         {                                               // series real space
-          register double r = sqrt(rsq);
-          register double x = g_ewald*r;
-          register double f = exp(-x*x)*qqrd2e;
+          double r = sqrt(rsq);
+          double x = g_ewald*r;
+          double f = exp(-x*x)*qqrd2e;
 
           B0 = 1.0/(1.0+EWALD_P*x);                     // eqn 2.8
           B0 *= ((((A5*B0+A4)*B0+A3)*B0+A2)*B0+A1)*f/r;
@@ -524,8 +524,8 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
 
       if (rsq < cut_ljsqi[typej]) {                     // lj
         if (order6) {                                   // long-range lj
-          register double rn = r2inv*r2inv*r2inv;
-          register double x2 = g2*rsq, a2 = 1.0/x2;
+          double rn = r2inv*r2inv*r2inv;
+          double x2 = g2*rsq, a2 = 1.0/x2;
           x2 = a2*exp(-x2)*lj4i[typej];
           if (ni < 0) {
             force_lj =
@@ -533,7 +533,7 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
             if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
           }
           else {                                        // special case
-            register double f = special_lj[ni], t = rn*(1.0-f);
+            double f = special_lj[ni], t = rn*(1.0-f);
             force_lj = f*(rn *= rn)*lj1i[typej]-
               g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
             if (eflag) evdwl =
@@ -541,13 +541,13 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
           }
         }
         else {                                          // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni < 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag) evdwl = f*(
                 rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -559,14 +559,14 @@ void PairLJLongDipoleLong::compute(int eflag, int vflag)
 
       fpair = force_coul+force_lj;                      // force
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1));
+        double *fj = f0+(j+(j<<1));
         fi[0] += fx = d[0]*fpair+force_d[0]; fj[0] -= fx;
         fi[1] += fy = d[1]*fpair+force_d[1]; fj[1] -= fy;
         fi[2] += fz = d[2]*fpair+force_d[2]; fj[2] -= fz;
         tqi[0] += mui[1]*ti[2]-mui[2]*ti[1];            // torque
         tqi[1] += mui[2]*ti[0]-mui[0]*ti[2];
         tqi[2] += mui[0]*ti[1]-mui[1]*ti[0];
-        register double *tqj = tq0+(j+(j<<1));
+        double *tqj = tq0+(j+(j<<1));
         tqj[0] += muj[1]*tj[2]-muj[2]*tj[1];
         tqj[1] += muj[2]*tj[0]-muj[0]*tj[2];
         tqj[2] += muj[0]*tj[1]-muj[1]*tj[0];
@@ -608,9 +608,9 @@ double PairLJLongDipoleLong::single(int i, int j, int itype, int jtype,
     double G0, G1, G2, B0, B1, B2, B3, mudi, mudj, muij;
     vector d = {xi[0]-xj[0], xi[1]-xj[1], xi[2]-xj[2]};
     {                                                   // series real space
-      register double r = sqrt(rsq);
-      register double x = g_ewald*r;
-      register double f = exp(-x*x)*qqrd2e;
+      double r = sqrt(rsq);
+      double x = g_ewald*r;
+      double f = exp(-x*x)*qqrd2e;
 
       B0 = 1.0/(1.0+EWALD_P*x);                 // eqn 2.8
       B0 *= ((((A5*B0+A4)*B0+A3)*B0+A2)*B0+A1)*f/r;
@@ -644,7 +644,7 @@ double PairLJLongDipoleLong::single(int i, int j, int itype, int jtype,
   if (rsq < cut_ljsq[itype][jtype]) {                   // lennard-jones
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&0x40) {                             // long-range
-      register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
+      double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
       x2 = a2*exp(-x2)*lj4[itype][jtype];
       force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
         g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp
index 8cdaa044aa737917ebfa9c04ff26391cbb2eb3eb..78d7b38e4d424a990b687dca472cb294b318b64a 100644
--- a/src/KSPACE/ewald_disp.cpp
+++ b/src/KSPACE/ewald_disp.cpp
@@ -778,7 +778,7 @@ void EwaldDisp::compute_ek()
         cek->re += zxyz.re*ci[i]; (cek++)->im += zxyz.im*ci[i];
       }
       if (func[3]) {
-        register double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h;
+        double muk = mui[0]*h->x+mui[1]*h->y+mui[2]*h->z; ++h;
         cek->re += zxyz.re*muk; (cek++)->im += zxyz.im*muk;
       }
     }
@@ -819,7 +819,7 @@ void EwaldDisp::compute_force()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -830,33 +830,33 @@ void EwaldDisp::compute_force()
       }
       C_CRMULT(zc, z[k->z].z, zxy);
       if (func[0]) {                                        // 1/r
-        register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re);
+        double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re);
         if (func[3]) cek_coul = cek;
         ++cek;
         sum[0][0] += h->x*im; sum[0][1] += h->y*im; sum[0][2] += h->z*im;
       }
       if (func[1]) {                                        // geometric 1/r^6
-        register double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek;
+        double im = *(ke++)*(zc.im*cek->re+cek->im*zc.re); ++cek;
         sum[1][0] += h->x*im; sum[1][1] += h->y*im; sum[1][2] += h->z*im;
       }
       if (func[2]) {                                        // arithmetic 1/r^6
-        register double im, c = *(ke++);
+        double im, c = *(ke++);
         for (i=2; i<9; ++i) {
           im = c*(zc.im*cek->re+cek->im*zc.re); ++cek;
           sum[i][0] += h->x*im; sum[i][1] += h->y*im; sum[i][2] += h->z*im;
         }
       }
       if (func[3]) {                                        // dipole
-        register double im = *(ke)*(zc.im*cek->re+
+        double im = *(ke)*(zc.im*cek->re+
             cek->im*zc.re)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
-        register double im2 = *(ke)*(zc.re*cek->re-
+        double im2 = *(ke)*(zc.re*cek->re-
             cek->im*zc.im);
         sum[9][0] += h->x*im; sum[9][1] += h->y*im; sum[9][2] += h->z*im;
         t[0] += -mui[1]*h->z*im2 + mui[2]*h->y*im2;        // torque
         t[1] += -mui[2]*h->x*im2 + mui[0]*h->z*im2;
         t[2] += -mui[0]*h->y*im2 + mui[1]*h->x*im2;
         if (func[0]) {                                      // charge-dipole
-          register double qi = *(q)*c[0];
+          double qi = *(q)*c[0];
           im = - *(ke)*(zc.re*cek_coul->re -
               cek_coul->im*zc.im)*(mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
           im += *(ke)*(zc.re*cek->re - cek->im*zc.im)*qi;
@@ -873,17 +873,17 @@ void EwaldDisp::compute_force()
       }
     }
     if (func[0]) {                                        // 1/r
-      register double qi = *(q++)*c[0];
+      double qi = *(q++)*c[0];
       f[0] -= sum[0][0]*qi; f[1] -= sum[0][1]*qi; f[2] -= sum[0][2]*qi;
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double bi = B[*type]*c[1];
+      double bi = B[*type]*c[1];
       f[0] -= sum[1][0]*bi; f[1] -= sum[1][1]*bi; f[2] -= sum[1][2]*bi;
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double *bi = B+7*type[0]+7;
+      double *bi = B+7*type[0]+7;
       for (i=2; i<9; ++i) {
-        register double c2 = (--bi)[0]*c[2];
+        double c2 = (--bi)[0]*c[2];
         f[0] -= sum[i][0]*c2; f[1] -= sum[i][1]*c2; f[2] -= sum[i][2]*c2;
       }
     }
@@ -962,7 +962,7 @@ void EwaldDisp::compute_energy()
     if (func[1]) {                                        // geometric 1/r^6
       sum[1] += *(ke++)*(cek->re*cek->re+cek->im*cek->im); ++cek; }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double r =
+      double r =
             (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+
             (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+
             (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+
@@ -1013,7 +1013,7 @@ void EwaldDisp::compute_energy_peratom()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(double));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -1031,7 +1031,7 @@ void EwaldDisp::compute_energy_peratom()
       if (func[1]) {                                        // geometric 1/r^6
         sum[1] += *(ke++)*(cek->re*zc.re - cek->im*zc.im); ++cek; }
       if (func[2]) {                                        // arithmetic 1/r^6
-        register double im, c = *(ke++);
+        double im, c = *(ke++);
         for (i=2; i<9; ++i) {
           im = c*(cek->re*zc.re - cek->im*zc.im); ++cek;
           sum[i] += im;
@@ -1041,7 +1041,7 @@ void EwaldDisp::compute_energy_peratom()
         double muk = (mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
         sum[9] += *(ke)*(cek->re*zc.re - cek->im*zc.im)*muk;
         if (func[0]) {                                      // charge-dipole
-          register double qj = *(q)*c[0];
+          double qj = *(q)*c[0];
           sum[9] += *(ke)*(cek_coul->im*zc.re + cek_coul->re*zc.im)*muk;
           sum[9] -= *(ke)*(cek->re*zc.im + cek->im*zc.re)*qj;
         }
@@ -1051,17 +1051,17 @@ void EwaldDisp::compute_energy_peratom()
     }
 
     if (func[0]) {                                        // 1/r
-      register double qj = *(q++)*c[0];
+      double qj = *(q++)*c[0];
       *eatomj += sum[0]*qj - energy_self_peratom[j][0];
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double bj = B[*type]*c[1];
+      double bj = B[*type]*c[1];
       *eatomj += sum[1]*bj - energy_self_peratom[j][1];
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double *bj = B+7*type[0]+7;
+      double *bj = B+7*type[0]+7;
       for (i=2; i<9; ++i) {
-        register double c2 = (--bj)[0]*c[2];
+        double c2 = (--bj)[0]*c[2];
         *eatomj += 0.5*sum[i]*c2;
       }
       *eatomj -= energy_self_peratom[j][2];
@@ -1097,19 +1097,19 @@ void EwaldDisp::compute_virial()
   memset(sum, 0, EWALD_NFUNCS*sizeof(shape));
   for (int k=0; k<nkvec; ++k) {                      // sum over k vectors
     if (func[0]) {                                         // 1/r
-      register double r = cek->re*cek->re+cek->im*cek->im;
+      double r = cek->re*cek->re+cek->im*cek->im;
       if (func[3]) cek_coul = cek;
       ++cek;
       sum[0][0] += *(kv++)*r; sum[0][1] += *(kv++)*r; sum[0][2] += *(kv++)*r;
       sum[0][3] += *(kv++)*r; sum[0][4] += *(kv++)*r; sum[0][5] += *(kv++)*r;
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double r = cek->re*cek->re+cek->im*cek->im; ++cek;
+      double r = cek->re*cek->re+cek->im*cek->im; ++cek;
       sum[1][0] += *(kv++)*r; sum[1][1] += *(kv++)*r; sum[1][2] += *(kv++)*r;
       sum[1][3] += *(kv++)*r; sum[1][4] += *(kv++)*r; sum[1][5] += *(kv++)*r;
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double r =
+      double r =
             (cek[0].re*cek[6].re+cek[0].im*cek[6].im)+
             (cek[1].re*cek[5].re+cek[1].im*cek[5].im)+
             (cek[2].re*cek[4].re+cek[2].im*cek[4].im)+
@@ -1118,12 +1118,12 @@ void EwaldDisp::compute_virial()
       sum[2][3] += *(kv++)*r; sum[2][4] += *(kv++)*r; sum[2][5] += *(kv++)*r;
     }
     if (func[3]) {
-      register double r = cek->re*cek->re+cek->im*cek->im;
+      double r = cek->re*cek->re+cek->im*cek->im;
       sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r;
       sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r;
       if (func[0]) {                                      // charge-dipole
         kv -= 6;
-        register double r = 2.0*(cek->re*cek_coul->im - cek->im*cek_coul->re);
+        double r = 2.0*(cek->re*cek_coul->im - cek->im*cek_coul->re);
         sum[3][0] += *(kv++)*r; sum[3][1] += *(kv++)*r; sum[3][2] += *(kv++)*r;
         sum[3][3] += *(kv++)*r; sum[3][4] += *(kv++)*r; sum[3][5] += *(kv++)*r;
       }
@@ -1173,7 +1173,7 @@ void EwaldDisp::compute_virial_dipole()
     cek = cek_global;
     memset(&sum[0], 0, 6*sizeof(double));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -1267,7 +1267,7 @@ void EwaldDisp::compute_virial_peratom()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(shape));
     if (func[3]) {
-      register double di = c[3];
+      double di = c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[0]; mui[2] = di*(mu++)[0];
       mu++;
     }
@@ -1279,7 +1279,7 @@ void EwaldDisp::compute_virial_peratom()
       C_CRMULT(zc, z[k->z].z, zxy);
       if (func[0]) {                                        // 1/r
           if (func[3]) cek_coul = cek;
-          register double r = cek->re*zc.re - cek->im*zc.im; ++cek;
+          double r = cek->re*zc.re - cek->im*zc.im; ++cek;
           sum[0][0] += *(kv++)*r;
           sum[0][1] += *(kv++)*r;
           sum[0][2] += *(kv++)*r;
@@ -1288,7 +1288,7 @@ void EwaldDisp::compute_virial_peratom()
           sum[0][5] += *(kv++)*r;
       }
       if (func[1]) {                                        // geometric 1/r^6
-          register double r = cek->re*zc.re - cek->im*zc.im; ++cek;
+          double r = cek->re*zc.re - cek->im*zc.im; ++cek;
           sum[1][0] += *(kv++)*r;
           sum[1][1] += *(kv++)*r;
           sum[1][2] += *(kv++)*r;
@@ -1297,7 +1297,7 @@ void EwaldDisp::compute_virial_peratom()
           sum[1][5] += *(kv++)*r;
       }
       if (func[2]) {                                        // arithmetic 1/r^6
-        register double r;
+        double r;
         for (i=2; i<9; ++i) {
           r = cek->re*zc.re - cek->im*zc.im; ++cek;
           sum[i][0] += *(kv++)*r;
@@ -1312,7 +1312,7 @@ void EwaldDisp::compute_virial_peratom()
       }
       if (func[3]) {                                        // dipole
          double muk = (mui[0]*h->x+mui[1]*h->y+mui[2]*h->z);
-         register double
+         double
            r = (cek->re*zc.re - cek->im*zc.im)*muk;
          sum[9][0] += *(kv++)*r;
          sum[9][1] += *(kv++)*r;
@@ -1322,7 +1322,7 @@ void EwaldDisp::compute_virial_peratom()
          sum[9][5] += *(kv++)*r;
          if (func[0]) {                                      // charge-dipole
            kv -= 6;
-           register double qj = *(q)*c[0];
+           double qj = *(q)*c[0];
            r = (cek_coul->im*zc.re + cek_coul->re*zc.im)*muk;
            r += -(cek->re*zc.im + cek->im*zc.re)*qj;
            sum[9][0] += *(kv++)*r; sum[9][1] += *(kv++)*r; sum[9][2] += *(kv++)*r;
@@ -1333,17 +1333,17 @@ void EwaldDisp::compute_virial_peratom()
     }
 
     if (func[0]) {                                        // 1/r
-      register double qi = *(q++)*c[0];
+      double qi = *(q++)*c[0];
       for (int n = 0; n < 6; n++) vatomj[n] += sum[0][n]*qi;
     }
     if (func[1]) {                                        // geometric 1/r^6
-      register double bi = B[*type]*c[1];
+      double bi = B[*type]*c[1];
       for (int n = 0; n < 6; n++) vatomj[n] += sum[1][n]*bi;
     }
     if (func[2]) {                                        // arithmetic 1/r^6
-      register double *bj = B+7*type[0]+7;
+      double *bj = B+7*type[0]+7;
       for (i=2; i<9; ++i) {
-        register double c2 = (--bj)[0]*c[2];
+        double c2 = (--bj)[0]*c[2];
         for (int n = 0; n < 6; n++) vatomj[n] += 0.5*sum[i][n]*c2;
       }
     }
diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp
index 55825b56571c492eca8c3262a75916e29f82ae60..d5bb7b6c5b6106686e2bd1f641b2a17de7118c05 100644
--- a/src/KSPACE/pair_buck_long_coul_long.cpp
+++ b/src/KSPACE/pair_buck_long_coul_long.cpp
@@ -483,7 +483,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -494,25 +494,25 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
             if (eflag) ecoul = t;
           }
           else {                                        // special case
-            register double f = s*(1.0-special_coul[ni])/r;
+            double f = s*(1.0-special_coul[ni])/r;
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
             if (eflag) ecoul = t-f;
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -527,11 +527,11 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
       else force_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (order6) {                                        // long-range
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -539,7 +539,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
               if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
               if (eflag) evdwl = f*expr*buckai[typej] -
@@ -547,16 +547,16 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
             }
           }
           else {                                              //table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
               if (eflag) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                                             //speial case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
               if (eflag) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -569,7 +569,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
                          rn*buckci[typej]-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
             if (eflag)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -581,7 +581,7 @@ void PairBuckLongCoulLong::compute(int eflag, int vflag)
       fpair = (force_coul+force_buck)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -640,7 +640,7 @@ void PairBuckLongCoulLong::compute_inner()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -654,7 +654,7 @@ void PairBuckLongCoulLong::compute_inner()
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -665,12 +665,12 @@ void PairBuckLongCoulLong::compute_inner()
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -730,7 +730,7 @@ void PairBuckLongCoulLong::compute_middle()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -745,7 +745,7 @@ void PairBuckLongCoulLong::compute_middle()
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -756,16 +756,16 @@ void PairBuckLongCoulLong::compute_middle()
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -831,7 +831,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -845,36 +845,36 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       respa_buck = 0.0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (r-cut_in_off)/cut_in_diff;
+        double rsw = (r-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double s = qri*q[j];
+          double s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
             if (eflag) ecoul = t;
           }
           else {                                        // correct for special
-            register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
+            double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul;
             if (eflag) ecoul = t-ri;
           }
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double s = qri*q[j];
+            double s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -892,14 +892,14 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (respa_flag) respa_buck = ni == 0 ?                 // correct for respa
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) :
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
         if (order6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -907,7 +907,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
               if (eflag) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck;
               if (eflag) evdwl = f*expr*buckai[typej] -
@@ -915,17 +915,17 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
             }
           }
           else {          // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck;
               if (eflag) evdwl =  expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                             //special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck;
               if (eflag) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -938,7 +938,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
               evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck;
             if (eflag)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -950,7 +950,7 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag)
       fpair = (force_coul+force_buck)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -985,17 +985,17 @@ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype,
 
   if ((ewald_order&2) && (rsq < cut_coulsq)) {                // coulombic
     if (!ncoultablebits || rsq <= tabinnersq) {                // series real space
-      register double x = g_ewald*r;
-      register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
+      double x = g_ewald*r;
+      double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
       f = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
       force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
       eng += t-f;
     }
     else {                                                // table real space
-      register union_int_float_t t;
+      union_int_float_t t;
       t.f = rsq;
-      register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-      register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
+      const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+      double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
       t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
       force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
       eng += qiqj*(etable[k]+f*detable[k]-t.f);
@@ -1003,10 +1003,10 @@ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype,
   } else force_coul = 0.0;
 
   if (rsq < cut_bucksq[itype][jtype]) {                        // buckingham
-    register double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]);
+    double expr = factor_buck*exp(-sqrt(rsq)*rhoinv[itype][jtype]);
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&64) {                                // long-range
-      register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_buck);
+      double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_buck);
       x2 = a2*exp(-x2)*buck_c[itype][jtype];
       force_buck = buck1[itype][jtype]*r*expr-
                g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*buck2[itype][jtype];
diff --git a/src/KSPACE/pair_lj_long_coul_long.cpp b/src/KSPACE/pair_lj_long_coul_long.cpp
index 32d6d9e0a3e82e602d838beade4b656496e689b3..2de527ab396e8f0a0a27bcdda23f62707b563181 100644
--- a/src/KSPACE/pair_lj_long_coul_long.cpp
+++ b/src/KSPACE/pair_lj_long_coul_long.cpp
@@ -478,7 +478,7 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -488,8 +488,8 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double r = sqrt(rsq), x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
@@ -502,10 +502,10 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask)>>ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask)>>ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -522,8 +522,8 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
       if (rsq < cut_ljsqi[typej]) {                        // lj
         if (order6) {                                        // long-range lj
           if(!ndisptablebits || rsq <= tabinnerdispsq) {                // series real space
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -532,7 +532,7 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
                 evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
               g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
               if (eflag)
@@ -540,30 +540,30 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
               if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
               if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
           }
         }
         else {                                                // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (eflag)
               evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -576,7 +576,7 @@ void PairLJLongCoulLong::compute(int eflag, int vflag)
       fpair = (force_coul+force_lj)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -635,7 +635,7 @@ void PairLJLongCoulLong::compute_inner()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -650,7 +650,7 @@ void PairLJLongCoulLong::compute_inner()
       }
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -660,12 +660,12 @@ void PairLJLongCoulLong::compute_inner()
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -725,7 +725,7 @@ void PairLJLongCoulLong::compute_middle()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -739,7 +739,7 @@ void PairLJLongCoulLong::compute_middle()
           qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni];
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -749,16 +749,16 @@ void PairLJLongCoulLong::compute_middle()
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -823,7 +823,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -836,16 +836,16 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       respa_lj = 0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (order1 && (rsq < cut_coulsq)) {                // coulombic
         if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), s = qri*q[j];
+          double r = sqrt(rsq), s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -859,13 +859,13 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -884,13 +884,13 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_ljsqi[typej]) {                        // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
             frespa*rn*(rn*lj1i[typej]-lj2i[typej]) :
             frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
         if (order6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -898,7 +898,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
               if (eflag) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj;
               if (eflag)
@@ -906,17 +906,17 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj;
               if (eflag) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj;
               if (eflag) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
@@ -928,7 +928,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
             if (eflag) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj;
             if (eflag)
               evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -940,7 +940,7 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag)
       fpair = (force_coul+force_lj)*r2inv;
 
       if (newton_pair || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -974,17 +974,17 @@ double PairLJLongCoulLong::single(int i, int j, int itype, int jtype,
   r2inv = 1.0/rsq;
   if ((ewald_order&2) && (rsq < cut_coulsq)) {                // coulombic
     if (!ncoultablebits || rsq <= tabinnersq) {                // series real space
-      register double r = sqrt(rsq), x = g_ewald*r;
-      register double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
+      double r = sqrt(rsq), x = g_ewald*r;
+      double s = force->qqrd2e*q[i]*q[j], t = 1.0/(1.0+EWALD_P*x);
       r = s*(1.0-factor_coul)/r; s *= g_ewald*exp(-x*x);
       force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-r;
       eng += t-r;
     }
     else {                                                // table real space
-      register union_int_float_t t;
+      union_int_float_t t;
       t.f = rsq;
-      register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-      register double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
+      const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+      double f = (rsq-rtable[k])*drtable[k], qiqj = q[i]*q[j];
       t.f = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
       force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
       eng += qiqj*(etable[k]+f*detable[k]-t.f);
@@ -994,7 +994,7 @@ double PairLJLongCoulLong::single(int i, int j, int itype, int jtype,
   if (rsq < cut_ljsq[itype][jtype]) {                        // lennard-jones
     r6inv = r2inv*r2inv*r2inv;
     if (ewald_order&64) {                                // long-range
-      register double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
+      double x2 = g2*rsq, a2 = 1.0/x2, t = r6inv*(1.0-factor_lj);
       x2 = a2*exp(-x2)*lj4[itype][jtype];
       force_lj = factor_lj*(r6inv *= r6inv)*lj1[itype][jtype]-
                g8*(((6.0*a2+6.0)*a2+3.0)*a2+a2)*x2*rsq+t*lj2[itype][jtype];
diff --git a/src/KSPACE/pair_lj_long_tip4p_long.cpp b/src/KSPACE/pair_lj_long_tip4p_long.cpp
index 78316b9daad147c3fb7b236760c13a99a08c89dd..a571e30fc377b6bafb48126e149576882aa4d27d 100644
--- a/src/KSPACE/pair_lj_long_tip4p_long.cpp
+++ b/src/KSPACE/pair_lj_long_tip4p_long.cpp
@@ -190,8 +190,8 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
         r2inv = 1.0/rsq;
         if (order6) {                   // long-range lj
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -200,7 +200,7 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
                 evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype];
               if (eflag)
@@ -208,30 +208,30 @@ void PairLJLongTIP4PLong::compute(int eflag, int vflag)
             }
           }
           else {                                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype];
               if (eflag) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype];
               if (eflag) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
           }
         }
         else {                      // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (eflag) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                    // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (eflag)
               evdwl = f * (rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -573,15 +573,15 @@ void PairLJLongTIP4PLong::compute_inner()
 
       if (rsq < cut_ljsq[itype][jtype] && rsq < cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq > cut_out_on_sq) {                        // switching
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -642,7 +642,7 @@ void PairLJLongTIP4PLong::compute_inner()
           }
 
           if (rsq > cut_out_on_sq) {                        // switching
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -826,19 +826,19 @@ void PairLJLongTIP4PLong::compute_middle()
 
       if (rsq < cut_ljsq[itype][jtype] && rsq >= cut_in_off_sq && rsq <= cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq < cut_in_on_sq) {                                // switching
-          register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
           forcelj  *= rsw*rsw*(3.0 - 2.0*rsw);
         }
         if (rsq > cut_out_on_sq) {
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -899,11 +899,11 @@ void PairLJLongTIP4PLong::compute_middle()
           }
 
           if (rsq < cut_in_on_sq) {                                // switching
-            register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
             forcecoul  *= rsw*rsw*(3.0 - 2.0*rsw);
           }
           if (rsq > cut_out_on_sq) {
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -1112,18 +1112,18 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
         frespa = 1.0;                                       // check whether and how to compute respa corrections
         respa_flag = rsq < cut_in_on_sq ? 1 : 0;
         if (respa_flag && (rsq > cut_in_off_sq)) {
-          register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
           frespa = 1-rsw*rsw*(3.0-2.0*rsw);
         }
 
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype]) :
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype])*special_lj[ni];
         if (order6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -1131,7 +1131,7 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
               if (eflag) evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype]-respa_lj;
               if (eflag)
@@ -1139,16 +1139,16 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]-respa_lj;
               if (eflag) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype]-respa_lj;
               if (eflag) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
@@ -1160,7 +1160,7 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
             if (eflag) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype])-respa_lj;
             if (eflag)
               evdwl = f*(rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -1225,16 +1225,16 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
           frespa = 1.0;                                       // check whether and how to compute respa corrections
           respa_flag = rsq < cut_in_on_sq ? 1 : 0;
           if (respa_flag && (rsq > cut_in_off_sq)) {
-            register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
             frespa = 1-rsw*rsw*(3.0-2.0*rsw);
           }
 
           r2inv = 1.0 / rsq;
           if (!ncoultablebits || rsq <= tabinnersq) {        // series real space
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             if (respa_flag)                                // correct for respa
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-            register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+            double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
             if (ni == 0) {
               s *= g_ewald*exp(-x*x);
               forcecoul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -1248,13 +1248,13 @@ void PairLJLongTIP4PLong::compute_outer(int eflag, int vflag)
           }                                                // table real space
           else {
             if (respa_flag) {
-              register double r = sqrt(rsq), s = qri*q[j];
+              double r = sqrt(rsq), s = qri*q[j];
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
             }
-            register union_int_float_t t;
+            union_int_float_t t;
             t.f = rsq;
-            register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-            register double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
+            const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+            double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
             if (ni == 0) {
               forcecoul = qiqj*(ftable[k]+f*dftable[k]);
               if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
diff --git a/src/OPT/pair_lj_long_coul_long_opt.cpp b/src/OPT/pair_lj_long_coul_long_opt.cpp
index 91ffc2c01ed20e96f317cf82cd48d778cb1233a0..6bd1e8cb9682411d4b3e2b7034beace732b95c2b 100644
--- a/src/OPT/pair_lj_long_coul_long_opt.cpp
+++ b/src/OPT/pair_lj_long_coul_long_opt.cpp
@@ -574,7 +574,7 @@ void PairLJLongCoulLongOpt::eval()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -584,8 +584,8 @@ void PairLJLongCoulLongOpt::eval()
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double r = sqrt(rsq), x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
@@ -598,10 +598,10 @@ void PairLJLongCoulLongOpt::eval()
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask)>>ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask)>>ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -618,8 +618,8 @@ void PairLJLongCoulLongOpt::eval()
       if (rsq < cut_ljsqi[typej]) {                        // lj
         if (ORDER6) {                                        // long-range lj
           if(!LJTABLE || rsq <= tabinnerdispsq) {               // series real space
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -628,7 +628,7 @@ void PairLJLongCoulLongOpt::eval()
                 evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
               g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
               if (EFLAG)
@@ -636,30 +636,30 @@ void PairLJLongCoulLongOpt::eval()
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
           }
         }
         else {                                                // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG)
               evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -671,7 +671,7 @@ void PairLJLongCoulLongOpt::eval()
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -742,7 +742,7 @@ void PairLJLongCoulLongOpt::eval_outer()
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register double *xj = x0+(j+(j<<1));
+      { double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -755,16 +755,16 @@ void PairLJLongCoulLongOpt::eval_outer()
       respa_lj = 0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), s = qri*q[j];
+          double r = sqrt(rsq), s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -778,13 +778,13 @@ void PairLJLongCoulLongOpt::eval_outer()
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -803,13 +803,13 @@ void PairLJLongCoulLongOpt::eval_outer()
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_ljsqi[typej]) {                        // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
             frespa*rn*(rn*lj1i[typej]-lj2i[typej]) :
             frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!LJTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -817,7 +817,7 @@ void PairLJLongCoulLongOpt::eval_outer()
               if (EFLAG) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj;
               if (EFLAG)
@@ -825,17 +825,17 @@ void PairLJLongCoulLongOpt::eval_outer()
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj;
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj;
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
@@ -847,7 +847,7 @@ void PairLJLongCoulLongOpt::eval_outer()
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj;
             if (EFLAG)
               evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -859,7 +859,7 @@ void PairLJLongCoulLongOpt::eval_outer()
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
diff --git a/src/USER-MISC/pair_edip.cpp b/src/USER-MISC/pair_edip.cpp
index f556a099c97370ddd350196a613e31f524d9087f..2976d7ad73c994fcb4237e6d05edd642f0f2c749 100644
--- a/src/USER-MISC/pair_edip.cpp
+++ b/src/USER-MISC/pair_edip.cpp
@@ -94,7 +94,7 @@ void PairEDIP::compute(int eflag, int vflag)
   int itype,jtype,ktype,ijparam,ikparam;
   double xtmp,ytmp,ztmp,evdwl;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  register int preForceCoord_counter;
+  int preForceCoord_counter;
 
   double invR_ij;
   double invR_ik;
diff --git a/src/USER-MISC/pair_edip_multi.cpp b/src/USER-MISC/pair_edip_multi.cpp
index 7a4193f55cbebb25944b8266a71444d3b265af0c..6be57eee74c9c32ffb92f28f27391b225e10d6ed 100644
--- a/src/USER-MISC/pair_edip_multi.cpp
+++ b/src/USER-MISC/pair_edip_multi.cpp
@@ -106,7 +106,7 @@ void PairEDIPMulti::compute(int eflag, int vflag)
   int itype,jtype,ktype,ijparam,ikparam,ijkparam;
   double xtmp,ytmp,ztmp,evdwl;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  register int preForceCoord_counter;
+  int preForceCoord_counter;
 
   double zeta_i;
   double dzetair;
diff --git a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
index 4a67cf9544fda593f82862e976498f1537d3a3c5..acb95ab76cd8c49f63940edb5c2300e9ae910aa1 100644
--- a/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_buck_long_coul_long_omp.cpp
@@ -678,7 +678,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -689,25 +689,25 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
             if (EFLAG) ecoul = t;
           }
           else {                                        // special case
-            register double f = s*(1.0-special_coul[ni])/r;
+            double f = s*(1.0-special_coul[ni])/r;
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-f;
             if (EFLAG) ecoul = t-f;
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -722,11 +722,11 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       else force_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (ORDER6) {                                        // long-range
           if (!DISPTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -734,7 +734,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
               if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej];
               if (EFLAG) evdwl = f*expr*buckai[typej] -
@@ -742,16 +742,16 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
             }
           }
           else {                                              //table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej];
               if (EFLAG) evdwl = expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                                             //speial case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej] -(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej] +t*buck2i[typej];
               if (EFLAG) evdwl = f*expr*buckai[typej] -(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -764,7 +764,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
                          rn*buckci[typej]-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej]);
             if (EFLAG)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -776,7 +776,7 @@ void PairBuckLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       fpair = (force_coul+force_buck)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -842,7 +842,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -856,7 +856,7 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -867,12 +867,12 @@ void PairBuckLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const t
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -939,7 +939,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -954,7 +954,7 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
           qri*q[j]/r : qri*q[j]/r*special_coul[ni];
 
       if (rsq < cut_bucksqi[typej = type[j]]) {                // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         force_buck = ni == 0 ?
           (r*expr*buck1i[typej]-rn*buck2i[typej]) :
@@ -965,16 +965,16 @@ void PairBuckLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const
       fpair = (force_coul + force_buck) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -1042,7 +1042,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -1056,36 +1056,36 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       respa_buck = 0.0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (r-cut_in_off)/cut_in_diff;
+        double rsw = (r-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double s = qri*q[j];
+          double s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
             if (EFLAG) ecoul = t;
           }
           else {                                        // correct for special
-            register double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
+            double ri = s*(1.0-special_coul[ni])/r; s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-ri-respa_coul;
             if (EFLAG) ecoul = t-ri;
           }
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double s = qri*q[j];
+            double s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -1103,14 +1103,14 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_bucksqi[typej]) {                        // buckingham
-        register double rn = r2inv*r2inv*r2inv,
+        double rn = r2inv*r2inv*r2inv,
                         expr = exp(-r*rhoinvi[typej]);
         if (respa_flag) respa_buck = ni == 0 ?                 // correct for respa
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej]) :
             frespa*(r*expr*buck1i[typej]-rn*buck2i[typej])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!DISPTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*buckci[typej];
             if (ni == 0) {
               force_buck =
@@ -1118,7 +1118,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
               if (EFLAG) evdwl = expr*buckai[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*buck2i[typej]-respa_buck;
               if (EFLAG) evdwl = f*expr*buckai[typej] -
@@ -1126,17 +1126,17 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
             }
           }
           else {          // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_buck = r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]-respa_buck;
               if (EFLAG) evdwl =  expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej];
             }
             else {                             //special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_buck = f*r*expr*buck1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*buckci[typej]+t*buck2i[typej]-respa_buck;
               if (EFLAG) evdwl = f*expr*buckai[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*buckci[typej]+t*buckci[typej];
             }
@@ -1149,7 +1149,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
               evdwl = expr*buckai[typej]-rn*buckci[typej]-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_buck = f*(r*expr*buck1i[typej]-rn*buck2i[typej])-respa_buck;
             if (EFLAG)
               evdwl = f*(expr*buckai[typej]-rn*buckci[typej]-offseti[typej]);
@@ -1161,7 +1161,7 @@ void PairBuckLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const t
       fpair = (force_coul+force_buck)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
diff --git a/src/USER-OMP/pair_edip_omp.cpp b/src/USER-OMP/pair_edip_omp.cpp
index 5ad8c808ba5e39a0ce728f1375065f296b247048..3392160d925f0e504b8a3fb8450f0202f681dd37 100644
--- a/src/USER-OMP/pair_edip_omp.cpp
+++ b/src/USER-OMP/pair_edip_omp.cpp
@@ -84,7 +84,7 @@ void PairEDIPOMP::eval(int iifrom, int iito, ThrData * const thr)
   int itype,jtype,ktype,ijparam,ikparam;
   double xtmp,ytmp,ztmp,evdwl;
   int *ilist,*jlist,*numneigh,**firstneigh;
-  register int preForceCoord_counter;
+  int preForceCoord_counter;
 
   double invR_ij;
   double invR_ik;
diff --git a/src/USER-OMP/pair_lj_long_coul_long_omp.cpp b/src/USER-OMP/pair_lj_long_coul_long_omp.cpp
index 0afc230e92c54c40f1a7578353ee367a80ef82b2..391840285b94e1d54d1d7ab5c98cb16285302954 100644
--- a/src/USER-OMP/pair_lj_long_coul_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_long_coul_long_omp.cpp
@@ -673,7 +673,7 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -683,8 +683,8 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), x = g_ewald*r;
-          register double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
+          double r = sqrt(rsq), x = g_ewald*r;
+          double s = qri*q[j], t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s;
@@ -697,10 +697,10 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
           }
         }                                                // table real space
         else {
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask)>>ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask)>>ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -717,8 +717,8 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       if (rsq < cut_ljsqi[typej]) {                        // lj
         if (ORDER6) {                                        // long-range lj
           if(!LJTABLE || rsq <= tabinnerdispsq) {            //series real space
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -727,7 +727,7 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
                 evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej];
               if (EFLAG)
@@ -735,30 +735,30 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
             }
           }
           else {                                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej];
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                                      // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej];
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
           }
         }
         else {                                                // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             force_lj = rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej]);
             if (EFLAG)
               evdwl = f * (rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -770,7 +770,7 @@ void PairLJLongCoulLongOMP::eval(int iifrom, int iito, ThrData * const thr)
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -834,7 +834,7 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -849,7 +849,7 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr
       }
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -859,12 +859,12 @@ void PairLJLongCoulLongOMP::eval_inner(int iifrom, int iito, ThrData * const thr
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq > cut_out_on_sq) {                        // switching
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -932,7 +932,7 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -946,7 +946,7 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
           qri*q[j]*sqrt(r2inv) : qri*q[j]*sqrt(r2inv)*special_coul[ni];
 
       if (rsq < cut_ljsqi[typej = type[j]]) {                // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         force_lj = ni == 0 ?
           rn*(rn*lj1i[typej]-lj2i[typej]) :
           rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
@@ -956,16 +956,16 @@ void PairLJLongCoulLongOMP::eval_middle(int iifrom, int iito, ThrData * const th
       fpair = (force_coul + force_lj) * r2inv;
 
       if (rsq < cut_in_on_sq) {                                // switching
-        register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
         fpair  *= rsw*rsw*(3.0 - 2.0*rsw);
       }
       if (rsq > cut_out_on_sq) {
-        register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+        double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
         fpair  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
       }
 
       if (newton_pair || j < nlocal) {                        // force update
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
@@ -1034,7 +1034,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       ni = sbmask(j);
       j &= NEIGHMASK;
 
-      { register const double *xj = x0+(j+(j<<1));
+      { const double *xj = x0+(j+(j<<1));
         d[0] = xi[0] - xj[0];                                // pair vector
         d[1] = xi[1] - xj[1];
         d[2] = xi[2] - xj[2]; }
@@ -1047,16 +1047,16 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       respa_lj = 0;
       respa_flag = rsq < cut_in_on_sq ? 1 : 0;
       if (respa_flag && (rsq > cut_in_off_sq)) {
-        register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+        double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
         frespa = 1-rsw*rsw*(3.0-2.0*rsw);
       }
 
       if (ORDER1 && (rsq < cut_coulsq)) {                // coulombic
         if (!CTABLE || rsq <= tabinnersq) {        // series real space
-          register double r = sqrt(rsq), s = qri*q[j];
+          double r = sqrt(rsq), s = qri*q[j];
           if (respa_flag)                                // correct for respa
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-          register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+          double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
           if (ni == 0) {
             s *= g_ewald*exp(-x*x);
             force_coul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -1070,13 +1070,13 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
         }                                                // table real space
         else {
           if (respa_flag) {
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
           }
-          register union_int_float_t t;
+          union_int_float_t t;
           t.f = rsq;
-          register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-          register double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
+          const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+          double f = (rsq-rtable[k])*drtable[k], qiqj = qi*q[j];
           if (ni == 0) {
             force_coul = qiqj*(ftable[k]+f*dftable[k]);
             if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
@@ -1095,13 +1095,13 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       else force_coul = respa_coul = ecoul = 0.0;
 
       if (rsq < cut_ljsqi[typej]) {                        // lennard-jones
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
             frespa*rn*(rn*lj1i[typej]-lj2i[typej]) :
             frespa*rn*(rn*lj1i[typej]-lj2i[typej])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!LJTABLE || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[typej];
             if (ni == 0) {
               force_lj =
@@ -1109,7 +1109,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
               if (EFLAG) evdwl = rn*lj3i[typej]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[typej]-respa_lj;
               if (EFLAG)
@@ -1117,17 +1117,17 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
             }
           }
           else {                                                // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               force_lj = (rn*=rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]-respa_lj;
               if (EFLAG) evdwl = rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej];
             }
             else {                                      // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               force_lj = f*(rn *= rn)*lj1i[typej]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[typej]+t*lj2i[typej]-respa_lj;
               if (EFLAG) evdwl = f*rn*lj3i[typej]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[typej]+t*lj4i[typej];
             }
@@ -1139,7 +1139,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
             if (EFLAG) evdwl = rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             force_lj = f*rn*(rn*lj1i[typej]-lj2i[typej])-respa_lj;
             if (EFLAG)
               evdwl = f*(rn*(rn*lj3i[typej]-lj4i[typej])-offseti[typej]);
@@ -1151,7 +1151,7 @@ void PairLJLongCoulLongOMP::eval_outer(int iiform, int iito, ThrData * const thr
       fpair = (force_coul+force_lj)*r2inv;
 
       if (NEWTON_PAIR || j < nlocal) {
-        register double *fj = f0+(j+(j<<1)), f;
+        double *fj = f0+(j+(j<<1)), f;
         fi[0] += f = d[0]*fpair; fj[0] -= f;
         fi[1] += f = d[1]*fpair; fj[1] -= f;
         fi[2] += f = d[2]*fpair; fj[2] -= f;
diff --git a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
index 8c0d72c74180a5dc61a32f6da8576196c5ee825e..19278d19bb6fbbe7c6fd70da5d94f49d4bffac14 100644
--- a/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
+++ b/src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
@@ -812,8 +812,8 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
         r2inv = 1.0/rsq;
         if (ORDER6) {                   // long-range lj
           if (!LJTABLE || rsq <= tabinnerdispsq) {
-            register double rn = r2inv*r2inv*r2inv;
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double rn = r2inv*r2inv*r2inv;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -822,7 +822,7 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
                 evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype];
               if (EFLAG)
@@ -830,30 +830,30 @@ void PairLJLongTIP4PLongOMP::eval(int iifrom, int iito, ThrData * const thr)
             }
           }
           else {                                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
-            register double rn = r2inv*r2inv*r2inv;
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            double rn = r2inv*r2inv*r2inv;
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype];
               if (EFLAG) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype];
               if (EFLAG) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
           }
         }
         else {                      // cut lj
-          register double rn = r2inv*r2inv*r2inv;
+          double rn = r2inv*r2inv*r2inv;
           if (ni == 0) {
             forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (EFLAG) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                    // special case
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
             if (EFLAG)
               evdwl = f * (rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -1191,15 +1191,15 @@ void PairLJLongTIP4PLongOMP::eval_inner(int iifrom, int iito, ThrData * const th
 
       if (rsq < cut_ljsq[itype][jtype] && rsq < cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq > cut_out_on_sq) {                        // switching
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -1260,7 +1260,7 @@ void PairLJLongTIP4PLongOMP::eval_inner(int iifrom, int iito, ThrData * const th
           }
 
           if (rsq > cut_out_on_sq) {                        // switching
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -1445,19 +1445,19 @@ void PairLJLongTIP4PLongOMP::eval_middle(int iifrom, int iito, ThrData * const t
 
       if (rsq < cut_ljsq[itype][jtype] && rsq >= cut_in_off_sq && rsq <= cut_out_off_sq ) {  // lj
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (ni == 0) forcelj = rn*(rn*lj1i[jtype]-lj2i[jtype]);
         else {                  // special case
-          register double f = special_lj[ni];
+          double f = special_lj[ni];
           forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype]);
         }
 
         if (rsq < cut_in_on_sq) {                                // switching
-          register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
           forcelj  *= rsw*rsw*(3.0 - 2.0*rsw);
         }
         if (rsq > cut_out_on_sq) {
-          register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+          double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
           forcelj  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
         }
 
@@ -1518,11 +1518,11 @@ void PairLJLongTIP4PLongOMP::eval_middle(int iifrom, int iito, ThrData * const t
           }
 
           if (rsq < cut_in_on_sq) {                                // switching
-            register double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
             forcecoul  *= rsw*rsw*(3.0 - 2.0*rsw);
           }
           if (rsq > cut_out_on_sq) {
-            register double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
+            double rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
             forcecoul  *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
           }
 
@@ -1719,18 +1719,18 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
         frespa = 1.0;                                       // check whether and how to compute respa corrections
         respa_flag = rsq < cut_in_on_sq ? 1 : 0;
         if (respa_flag && (rsq > cut_in_off_sq)) {
-          register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+          double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
           frespa = 1-rsw*rsw*(3.0-2.0*rsw);
         }
 
         r2inv = 1.0/rsq;
-        register double rn = r2inv*r2inv*r2inv;
+        double rn = r2inv*r2inv*r2inv;
         if (respa_flag) respa_lj = ni == 0 ?                 // correct for respa
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype]) :
                           frespa*rn*(rn*lj1i[jtype]-lj2i[jtype])*special_lj[ni];
         if (ORDER6) {                                        // long-range form
           if (!ndisptablebits || rsq <= tabinnerdispsq) {
-            register double x2 = g2*rsq, a2 = 1.0/x2;
+            double x2 = g2*rsq, a2 = 1.0/x2;
             x2 = a2*exp(-x2)*lj4i[jtype];
             if (ni == 0) {
               forcelj =
@@ -1738,7 +1738,7 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
               if (EFLAG) evdwl = rn*lj3i[jtype]-g6*((a2+1.0)*a2+0.5)*x2;
             }
             else {                                        // correct for special
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-
                 g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq+t*lj2i[jtype]-respa_lj;
               if (EFLAG)
@@ -1746,16 +1746,16 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
             }
           }
           else {                        // table real space
-            register union_int_float_t disp_t;
+            union_int_float_t disp_t;
             disp_t.f = rsq;
-            register const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
-            register double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
+            const int disp_k = (disp_t.i & ndispmask)>>ndispshiftbits;
+            double f_disp = (rsq-rdisptable[disp_k])*drdisptable[disp_k];
             if (ni == 0) {
               forcelj = (rn*=rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]-respa_lj;
               if (EFLAG) evdwl = rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype];
             }
             else {                  // special case
-              register double f = special_lj[ni], t = rn*(1.0-f);
+              double f = special_lj[ni], t = rn*(1.0-f);
               forcelj = f*(rn *= rn)*lj1i[jtype]-(fdisptable[disp_k]+f_disp*dfdisptable[disp_k])*lj4i[jtype]+t*lj2i[jtype]-respa_lj;
               if (EFLAG) evdwl = f*rn*lj3i[jtype]-(edisptable[disp_k]+f_disp*dedisptable[disp_k])*lj4i[jtype]+t*lj4i[jtype];
             }
@@ -1767,7 +1767,7 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
             if (EFLAG) evdwl = rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype];
           }
           else {                                        // correct for special
-            register double f = special_lj[ni];
+            double f = special_lj[ni];
             forcelj = f*rn*(rn*lj1i[jtype]-lj2i[jtype])-respa_lj;
             if (EFLAG)
               evdwl = f*(rn*(rn*lj3i[jtype]-lj4i[jtype])-offseti[jtype]);
@@ -1832,16 +1832,16 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
           frespa = 1.0;                                       // check whether and how to compute respa corrections
           respa_flag = rsq < cut_in_on_sq ? 1 : 0;
           if (respa_flag && (rsq > cut_in_off_sq)) {
-            register double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
+            double rsw = (sqrt(rsq)-cut_in_off)/cut_in_diff;
             frespa = 1-rsw*rsw*(3.0-2.0*rsw);
           }
 
           r2inv = 1.0 / rsq;
           if (!CTABLE || rsq <= tabinnersq) {        // series real space
-            register double r = sqrt(rsq), s = qri*q[j];
+            double r = sqrt(rsq), s = qri*q[j];
             if (respa_flag)                                // correct for respa
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
-            register double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
+            double x = g_ewald*r, t = 1.0/(1.0+EWALD_P*x);
             if (ni == 0) {
               s *= g_ewald*exp(-x*x);
               forcecoul = (t *= ((((t*A5+A4)*t+A3)*t+A2)*t+A1)*s/x)+EWALD_F*s-respa_coul;
@@ -1855,13 +1855,13 @@ void PairLJLongTIP4PLongOMP::eval_outer(int iifrom, int iito, ThrData * const th
           }                                                // table real space
           else {
             if (respa_flag) {
-              register double r = sqrt(rsq), s = qri*q[j];
+              double r = sqrt(rsq), s = qri*q[j];
               respa_coul = ni == 0 ? frespa*s/r : frespa*s/r*special_coul[ni];
             }
-            register union_int_float_t t;
+            union_int_float_t t;
             t.f = rsq;
-            register const int k = (t.i & ncoulmask) >> ncoulshiftbits;
-            register double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
+            const int k = (t.i & ncoulmask) >> ncoulshiftbits;
+            double f = (t.f-rtable[k])*drtable[k], qiqj = qtmp*q[j];
             if (ni == 0) {
               forcecoul = qiqj*(ftable[k]+f*dftable[k]);
               if (EFLAG) ecoul = qiqj*(etable[k]+f*detable[k]);
diff --git a/src/math_complex.h b/src/math_complex.h
index e4d889c4e6fa438a5febe706a1041af9d85ec648..4cf25b6dacaa83cdccfd2c9c732607b1ca4ca84c 100644
--- a/src/math_complex.h
+++ b/src/math_complex.h
@@ -32,12 +32,12 @@ typedef struct complex {
   d.im = x.re*y.im+x.im*y.re; }
 
 #define C_RMULT(d, x, y) { \
-  register complex t = x; \
+  complex t = x; \
   d.re = t.re*y.re-t.im*y.im; \
   d.im = t.re*y.im+t.im*y.re; }
 
 #define C_CRMULT(d, x, y) { \
-  register complex t = x; \
+  complex t = x; \
   d.re = t.re*y.re-t.im*y.im; \
   d.im = -t.re*y.im-t.im*y.re; }
 
@@ -62,7 +62,7 @@ typedef struct complex {
   d.im = y; }
 
 #define C_ANGLE(d, angle) { \
-  register double a = angle; \
+  double a = angle; \
   d.re = cos(a); \
   d.im = sin(a); }
 
diff --git a/src/math_vector.h b/src/math_vector.h
index c2227572729a96367f3ca35d1d11fd12c1391592..ff260b794d19a7ee4c599420e4a88ca278465e84 100644
--- a/src/math_vector.h
+++ b/src/math_vector.h
@@ -62,7 +62,7 @@ inline void vec_neg(vector &dest) {                                // -a
   dest[2] = -dest[2]; }
 
 inline void vec_norm(vector &dest) {                                 // a/|a|
-  register double f = sqrt(vec_dot(dest, dest));
+  double f = sqrt(vec_dot(dest, dest));
   dest[0] /= f;
   dest[1] /= f;
   dest[2] /= f; }
@@ -190,7 +190,7 @@ inline void form_subtr(shape &dest, form &src) {                // m_a-m_b
   dest[3] -= src[3]; dest[4] -= src[4]; dest[5] -= src[5]; }
 
 inline int form_inv(form &m_inv, form &m) {                        // m^-1
-  register double det = form_det(m);
+  double det = form_det(m);
   if (fzero(det)) return 0;
   m_inv[0] = (m[1]*m[2]-m[3]*m[3])/det;
   m_inv[1] = (m[0]*m[2]-m[4]*m[4])/det;
@@ -345,7 +345,7 @@ inline void form4_unit(form4 &dest) {
   dest[0] = dest[1] = dest[2] = dest[3] = 1.0; }
 
 inline double form4_det(form4 &m) {
-  register double f = m[6]*m[7]-m[5]*m[8];
+  double f = m[6]*m[7]-m[5]*m[8];
   return m[0]*(
       m[1]*(m[2]*m[3]-m[4]*m[4])+
       m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])+f*f+
@@ -355,7 +355,7 @@ inline double form4_det(form4 &m) {
         m[9]*(m[4]*m[4]-m[2]*m[3])); }
 
 inline int form4_inv(form4 &m_inv, form4 &m) {
-  register double det = form4_det(m);
+  double det = form4_det(m);
   if (fzero(det)) return 0;
   m_inv[0] = (m[1]*(m[2]*m[3]-m[4]*m[4])+
       m[5]*(2.0*m[4]*m[7]-m[2]*m[5])-m[3]*m[7]*m[7])/det;
diff --git a/src/pack.h b/src/pack.h
index 5db35117c55162ed2045ca37f1b2850a1c8922e3..066535f5c9fc576094ef9a059e090a2d10393286 100644
--- a/src/pack.h
+++ b/src/pack.h
@@ -55,8 +55,8 @@ struct pack_plan_3d {
 
 static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -81,8 +81,8 @@ static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 
 static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -107,8 +107,8 @@ static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan
 
 static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -133,8 +133,8 @@ static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 
 static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -162,8 +162,8 @@ static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,iqty,instart,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
+  int in,out,iqty,instart,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -192,8 +192,8 @@ static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -218,8 +218,8 @@ static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  int in,out,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -246,8 +246,8 @@ static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register int in,out,iqty,instart,fast,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
+  int in,out,iqty,instart,fast,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -283,9 +283,9 @@ static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -312,9 +312,9 @@ static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -341,9 +341,9 @@ static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan
 static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -370,9 +370,9 @@ static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -401,9 +401,9 @@ static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -433,9 +433,9 @@ static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -461,9 +461,9 @@ static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -491,9 +491,9 @@ static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -532,9 +532,9 @@ static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out;
-  register int mid,slow,size;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
+  PACK_DATA *in,*out;
+  int mid,slow,size;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -561,9 +561,9 @@ static void pack_3d(PACK_DATA *data, PACK_DATA *buf, struct pack_plan_3d *plan)
 static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out;
-  register int mid,slow,size;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
+  PACK_DATA *in,*out;
+  int mid,slow,size;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,upto;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -590,9 +590,9 @@ static void unpack_3d(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan
 static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -619,9 +619,9 @@ static void unpack_3d_permute1_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -650,9 +650,9 @@ static void unpack_3d_permute1_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -682,9 +682,9 @@ static void unpack_3d_permute1_n(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -710,9 +710,9 @@ static void unpack_3d_permute2_1(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*begin,*end;
-  register int mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane;
+  PACK_DATA *in,*out,*begin,*end;
+  int mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
@@ -740,9 +740,9 @@ static void unpack_3d_permute2_2(PACK_DATA *buf, PACK_DATA *data, struct pack_pl
 static void unpack_3d_permute2_n(PACK_DATA *buf, PACK_DATA *data, struct pack_plan_3d *plan)
 
 {
-  register PACK_DATA *in,*out,*instart,*begin,*end;
-  register int iqty,mid,slow;
-  register int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
+  PACK_DATA *in,*out,*instart,*begin,*end;
+  int iqty,mid,slow;
+  int nfast,nmid,nslow,nstride_line,nstride_plane,nqty;
 
   nfast = plan->nfast;
   nmid = plan->nmid;
diff --git a/src/pair.cpp b/src/pair.cpp
index 3b839dff497808a033a97f9c551381b80658e267..73dbb9428e37939a7a407ad70d8b9ea516a8a133 100644
--- a/src/pair.cpp
+++ b/src/pair.cpp
@@ -563,7 +563,7 @@ void Pair::init_tables_disp(double cut_lj_global)
       rsq_lookup.i |= maskhi;
     }
     rsq = rsq_lookup.f;
-    register double x2 = g2*rsq, a2 = 1.0/x2;
+    double x2 = g2*rsq, a2 = 1.0/x2;
     x2 = a2*exp(-x2);
 
     rdisptable[i] = rsq_lookup.f;
@@ -609,7 +609,7 @@ void Pair::init_tables_disp(double cut_lj_global)
   if (rsq_lookup.f < (cut_lj_globalsq = cut_lj_global * cut_lj_global)) {
     rsq_lookup.f = cut_lj_globalsq;
 
-    register double x2 = g2*rsq, a2 = 1.0/x2;
+    double x2 = g2*rsq, a2 = 1.0/x2;
     x2 = a2*exp(-x2);
     f_tmp = g8*(((6.0*a2+6.0)*a2+3.0)*a2+1.0)*x2*rsq;
     e_tmp = g6*((a2+1.0)*a2+0.5)*x2;