jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

remove deprecated 'register' keyword. closes #922

This commit is contained in:
Axel Kohlmeyer 2018-05-18 17:43:18 -04:00
parent 94c6d2d546
commit 5daf1fe0d4
16 changed files with 423 additions and 423 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
src/DIPOLE/pair_lj_long_dipole_long.cpp
View File

@ -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];

72
src/KSPACE/ewald_disp.cpp
View File

@ -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;
}
}

102
src/KSPACE/pair_buck_long_coul_long.cpp
View File

@ -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];

100
src/KSPACE/pair_lj_long_coul_long.cpp
View File

@ -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];

72
src/KSPACE/pair_lj_long_tip4p_long.cpp
View File

@ -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]);

70
src/OPT/pair_lj_long_coul_long_opt.cpp
View File

@ -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;

2
src/USER-MISC/pair_edip.cpp
View File

@ -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;

2
src/USER-MISC/pair_edip_multi.cpp
View File

@ -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;

88
src/USER-OMP/pair_buck_long_coul_long_omp.cpp
View File

@ -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;

2
src/USER-OMP/pair_edip_omp.cpp
View File

@ -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;

88
src/USER-OMP/pair_lj_long_coul_long_omp.cpp
View File

@ -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;

72
src/USER-OMP/pair_lj_long_tip4p_long_omp.cpp
View File

@ -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]);

6
src/math_complex.h
View File

@ -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); }

8
src/math_vector.h
View File

@ -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;

128
src/pack.h
View File

@ -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;

4
src/pair.cpp
View File

@ -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;