Unverified Commit 933cf92e authored by Steve Plimpton's avatar Steve Plimpton Committed by GitHub
Browse files

Merge pull request #923 from akohlmey/remove-register 

Remove deprecated 'register' keyword
parents 69903cb4 5daf1fe0

src/DIPOLE/pair_lj_long_dipole_long.cpp

+15 −15
Original line number Diff line number Diff line
@@ -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];

src/KSPACE/ewald_disp.cpp

+36 −36
Original line number Diff line number Diff line
@@ -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;

      }

    }

src/KSPACE/pair_buck_long_coul_long.cpp

+51 −51

File changed.

Preview size limit exceeded, changes collapsed.

src/KSPACE/pair_lj_long_coul_long.cpp

+50 −50

File changed.

Preview size limit exceeded, changes collapsed.

src/KSPACE/pair_lj_long_tip4p_long.cpp

+36 −36
Original line number Diff line number Diff line
@@ -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]);

CRA Git | Maintained and supported by SUSTech CRA and CCSE