lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3245 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2009-10-29 21:50:07 +00:00
parent ac33e812e8
commit 64c45b8d56
12 changed files with 317 additions and 320 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
src/KSPACE/pair_coul_long.cpp
View File

@ -73,8 +73,7 @@ void PairCoulLong::compute(int eflag, int vflag)
double r,r2inv,forcecoul,factor_coul;
double grij,expm2,prefactor,t,erfc;
int *ilist,*jlist,*numneigh,**firstneigh;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -118,11 +117,11 @@ void PairCoulLong::compute(int eflag, int vflag)
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < cut_coulsq) {
r2inv = 1.0/rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r2inv = 1.0/rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -131,9 +130,11 @@ void PairCoulLong::compute(int eflag, int vflag)
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
@ -310,39 +311,37 @@ void PairCoulLong::init_tables()
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
table_lookup_t rsq_lookup;
table_lookup_t minrsq_lookup;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq);
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -352,10 +351,10 @@ void PairCoulLong::init_tables()
}
}
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq;
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
@ -391,18 +390,18 @@ void PairCoulLong::init_tables()
// or ntablem1 if itablemin=0
// deltas at itablemax only needed if corresponding rsq < cut*cut
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq < cut_coulsq) {
rsq = cut_coulsq;
r = sqrtf(rsq);
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
@ -417,8 +416,8 @@ void PairCoulLong::init_tables()
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -429,7 +428,7 @@ void PairCoulLong::init_tables()
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
@ -529,11 +528,11 @@ double PairCoulLong::single(int i, int j, int itype, int jtype,
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
float rsq_single = rsq;
int *int_rsq = (int *) &rsq_single;
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_single - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = atom->q[i]*atom->q[j] * table;
if (factor_coul < 1.0) {

84
src/KSPACE/pair_lj_charmm_coul_long.cpp
View File

@ -90,8 +90,7 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
double grij,expm2,prefactor,t,erfc;
double philj,switch1,switch2;
int *ilist,*jlist,*numneigh,**firstneigh;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -145,7 +144,7 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -154,9 +153,11 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
@ -431,8 +432,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
double philj,switch1,switch2;
double rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -494,7 +494,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -515,9 +515,11 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
}
}
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
@ -540,7 +542,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = forcelj*switch1 + philj*switch2;
}
if (rsq < cut_in_on_sq) {
rsw = (sqrtf(rsq) - cut_in_off)/cut_in_diff;
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
}
} else forcelj = 0.0;
@ -888,39 +890,37 @@ void PairLJCharmmCoulLong::init_tables()
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
table_lookup_t rsq_lookup;
table_lookup_t minrsq_lookup;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq);
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -930,10 +930,10 @@ void PairLJCharmmCoulLong::init_tables()
}
}
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq;
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
@ -971,16 +971,16 @@ void PairLJCharmmCoulLong::init_tables()
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq < cut_coulsq) {
rsq = cut_coulsq;
r = sqrtf(rsq);
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
@ -995,8 +995,8 @@ void PairLJCharmmCoulLong::init_tables()
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -1007,7 +1007,7 @@ void PairLJCharmmCoulLong::init_tables()
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
@ -1146,11 +1146,11 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
float rsq_single = rsq;
int *int_rsq = (int *) &rsq_single;
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_single - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = atom->q[i]*atom->q[j] * table;
if (factor_coul < 1.0) {

92
src/KSPACE/pair_lj_cut_coul_long.cpp
View File

@ -85,8 +85,7 @@ void PairLJCutCoulLong::compute(int eflag, int vflag)
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
int *ilist,*jlist,*numneigh,**firstneigh;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -141,7 +140,7 @@ void PairLJCutCoulLong::compute(int eflag, int vflag)
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -150,9 +149,11 @@ void PairLJCutCoulLong::compute(int eflag, int vflag)
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
@ -406,8 +407,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
double grij,expm2,prefactor,t,erfc;
double rsw;
int *ilist,*jlist,*numneigh,**firstneigh;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -469,7 +469,7 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -479,9 +479,9 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
if (rsq > cut_in_off_sq) {
if (rsq < cut_in_on_sq) {
rsw = (r - cut_in_off)/cut_in_diff;
forcecoul += prefactor*rsw*rsw*(3 - 2*rsw);
forcecoul += prefactor*rsw*rsw*(3.0 - 2.0*rsw);
if (factor_coul < 1.0)
forcecoul -= (1.0-factor_coul)*prefactor*rsw*rsw*(3 - 2*rsw);
forcecoul -= (1.0-factor_coul)*prefactor*rsw*rsw*(3.0 - 2.0*rsw);
} else {
forcecoul += prefactor;
if (factor_coul < 1.0)
@ -489,9 +489,11 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
}
}
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {
@ -506,8 +508,8 @@ void PairLJCutCoulLong::compute_outer(int eflag, int vflag)
r6inv = r2inv*r2inv*r2inv;
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq < cut_in_on_sq) {
rsw = (sqrtf(rsq) - cut_in_off)/cut_in_diff;
forcelj *= rsw*rsw*(3 - 2*rsw);
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
forcelj *= rsw*rsw*(3.0 - 2.0*rsw);
}
} else forcelj = 0.0;
@ -847,39 +849,37 @@ void PairLJCutCoulLong::init_tables()
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
table_lookup_t rsq_lookup;
table_lookup_t minrsq_lookup;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq);
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -889,10 +889,10 @@ void PairLJCutCoulLong::init_tables()
}
}
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq;
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
@ -928,18 +928,18 @@ void PairLJCutCoulLong::init_tables()
// or ntablem1 if itablemin=0
// deltas at itablemax only needed if corresponding rsq < cut*cut
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq < cut_coulsq) {
rsq = cut_coulsq;
r = sqrtf(rsq);
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
@ -954,8 +954,8 @@ void PairLJCutCoulLong::init_tables()
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -966,7 +966,7 @@ void PairLJCutCoulLong::init_tables()
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
@ -1099,11 +1099,11 @@ double PairLJCutCoulLong::single(int i, int j, int itype, int jtype,
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor;
} else {
float rsq_single = rsq;
int *int_rsq = (int *) &rsq_single;
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup_single;
rsq_lookup_single.f = rsq;
itable = rsq_lookup_single.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq_single - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup_single.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = atom->q[i]*atom->q[j] * table;
if (factor_coul < 1.0) {

14
src/KSPACE/pair_lj_cut_coul_long_tip4p.cpp
View File

@ -77,8 +77,7 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
double xiM[3],xjM[3],fO[3],fH[3],v[6];
double *x1,*x2;
int *ilist,*jlist,*numneigh,**firstneigh;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
evdwl = ecoul = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -174,9 +173,9 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
// test current rsq against cutoff and compute Coulombic force
if (rsq < cut_coulsq) {
r2inv = 1 / rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r2inv = 1 / rsq;
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -187,10 +186,11 @@ void PairLJCutCoulLongTIP4P::compute(int eflag, int vflag)
forcecoul -= (1.0-factor_coul)*prefactor;
}
} else {
r2inv = 1 / rsq;
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (factor_coul < 1.0) {

21
src/OPT/pair_lj_charmm_coul_long_opt.h
View File

@ -65,9 +65,8 @@ void PairLJCharmmCoulLongOpt::eval()
double r,r2inv,r6inv,forcecoul,forcelj,factor_coul,factor_lj;
double grij,expm2,prefactor,t,erfc;
double philj,switch1,switch2;
float rsq;
int *int_rsq = (int *) &rsq;
double rsq;
double evdwl = 0.0;
double ecoul = 0.0;
@ -142,7 +141,7 @@ void PairLJCharmmCoulLongOpt::eval()
forcecoul = 0.0;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -152,9 +151,11 @@ void PairLJCharmmCoulLongOpt::eval()
prefactor = qqrd2e * tmp_coef3/r;
forcecoul = prefactor * (erfc + EWALD_F*grij*expm2);
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = tmp_coef3 * table;
}
@ -228,7 +229,7 @@ void PairLJCharmmCoulLongOpt::eval()
forcecoul = 0.0;
if (rsq < cut_coulsq) {
if (!ncoultablebits || rsq <= tabinnersq) {
r = sqrtf(rsq);
r = sqrt(rsq);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
t = 1.0 / (1.0 + EWALD_P*grij);
@ -241,9 +242,11 @@ void PairLJCharmmCoulLongOpt::eval()
forcecoul -= (1.0-factor_coul)*prefactor;
}
} else {
itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
fraction = (rsq - rtable[itable]) * drtable[itable];
fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
table = ftable[itable] + fraction*dftable[itable];
forcecoul = tmp_coef3 * table;
if (factor_coul < 1.0) {

52
src/USER-CG-CMM/pair_cg_cmm_coul_long.cpp
View File

@ -163,38 +163,36 @@ void PairCGCMMCoulLong::init_tables()
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
table_lookup_t rsq_lookup;
table_lookup_t minrsq_lookup;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq);
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -204,9 +202,9 @@ void PairCGCMMCoulLong::init_tables()
}
}
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq;
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
@ -244,15 +242,15 @@ void PairCGCMMCoulLong::init_tables()
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
if (rsq < cut_coulsq_global) {
rsq = cut_coulsq_global;
r = sqrtf(rsq);
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq_lookup.f < cut_coulsq_global) {
rsq_lookup.f = cut_coulsq_global;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
@ -267,8 +265,8 @@ void PairCGCMMCoulLong::init_tables()
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -279,7 +277,7 @@ void PairCGCMMCoulLong::init_tables()
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];

20
src/USER-CG-CMM/pair_cmm_common.h
View File

@ -257,9 +257,8 @@ namespace LAMMPS_NS {
if (COUL_TYPE == CG_COUL_LONG) {
if (rsq < cut_coulsq_global) {
const float rsqf = rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
const float r = sqrtf(rsq);
const double r = sqrt(rsq);
const double grij = g_ewald * r;
const double expm2 = exp(-grij*grij);
@ -273,10 +272,11 @@ namespace LAMMPS_NS {
if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
}
} else {
const int *int_rsq = (int *) &rsqf;
int itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
int itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
const double fraction = (rsq - rtable[itable]) * drtable[itable];
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
const double table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (EFLAG) {
@ -654,9 +654,8 @@ namespace LAMMPS_NS {
if (COUL_TYPE == CG_COUL_LONG) {
if (rsq < cut_coulsq_global) {
const float rsqf = rsq;
if (!ncoultablebits || rsq <= tabinnersq) {
const float r = sqrtf(rsq);
const double r = sqrt(rsq);
const double grij = g_ewald * r;
const double expm2 = exp(-grij*grij);
@ -670,10 +669,11 @@ namespace LAMMPS_NS {
if (EFLAG) ecoul -= (1.0-factor_coul)*prefactor;
}
} else {
const int *int_rsq = (int *) &rsqf;
int itable = *int_rsq & ncoulmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
int itable = rsq_lookup.i & ncoulmask;
itable >>= ncoulshiftbits;
const double fraction = (rsq - rtable[itable]) * drtable[itable];
const double fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable];
const double table = ftable[itable] + fraction*dftable[itable];
forcecoul = qtmp*q[j] * table;
if (EFLAG) {

86
src/USER-EWALDN/pair_buck_coul.cpp
View File

@ -224,7 +224,6 @@ void PairBuckCoul::coeff(int narg, char **arg)
void PairBuckCoul::init_style()
{
int i,j;
// require an atom style with charge defined
@ -508,17 +507,18 @@ void PairBuckCoul::compute(int eflag, int vflag)
}
} // table real space
else {
register float t = rsq;
register int k = (((int *) &t)[0] & ncoulmask)>>ncoulshiftbits;
register table_lookup_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];
if (ni < 0) {
force_coul = qiqj*(ftable[k]+f*dftable[k]);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
}
else { // special case
t = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t);
t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
}
}
}
@ -848,17 +848,18 @@ void PairBuckCoul::compute_outer(int eflag, int vflag)
if (respa_flag) respa_coul = ni<0 ? // correct for respa
frespa*qri*q[j]/r :
frespa*qri*q[j]/r*special_coul[ni];
register float t = rsq;
register int k = (((int *) &t)[0] & ncoulmask)>>ncoulshiftbits;
register table_lookup_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];
if (ni < 0) {
force_coul = qiqj*(ftable[k]+f*dftable[k]);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
}
else { // correct for special
t = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t);
t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
}
}
}
@ -964,39 +965,37 @@ void PairBuckCoul::init_tables()
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
table_lookup_t rsq_lookup;
table_lookup_t minrsq_lookup;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrt(rsq);
r = sqrt(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -1006,10 +1005,10 @@ void PairBuckCoul::init_tables()
}
}
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq;
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
@ -1047,16 +1046,16 @@ void PairBuckCoul::init_tables()
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq < cut_coulsq) {
rsq = cut_coulsq;
r = sqrt(rsq);
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
@ -1071,8 +1070,8 @@ void PairBuckCoul::init_tables()
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -1083,7 +1082,7 @@ void PairBuckCoul::init_tables()
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
@ -1136,12 +1135,13 @@ double PairBuckCoul::single(int i, int j, int itype, int jtype,
eng += t-f;
}
else { // table real space
register float t = rsq;
register int k = (((int *) &t)[0] & ncoulmask)>>ncoulshiftbits;
register table_lookup_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];
t = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
eng += qiqj*(etable[k]+f*detable[k]-t);
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);
}
} else force_coul = 0.0;

87
src/USER-EWALDN/pair_lj_coul.cpp
View File

@ -220,7 +220,7 @@ void PairLJCoul::init_style()
{
char *style1[] = {"ewald", "ewald/n", "pppm", NULL};
char *style6[] = {"ewald/n", NULL};
int i,j;
int i;
// require an atom style with charge defined
@ -508,17 +508,18 @@ void PairLJCoul::compute(int eflag, int vflag)
}
} // table real space
else {
register float t = rsq;
register int k = (((int *) &t)[0] & ncoulmask)>>ncoulshiftbits;
register table_lookup_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];
if (ni < 0) {
force_coul = qiqj*(ftable[k]+f*dftable[k]);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
}
else { // special case
t = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t);
t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
}
}
}
@ -841,17 +842,18 @@ void PairLJCoul::compute_outer(int eflag, int vflag)
if (respa_flag) respa_coul = ni<0 ? // correct for respa
frespa*qri*q[j]/sqrt(rsq) :
frespa*qri*q[j]/sqrt(rsq)*special_coul[ni];
register float t = rsq;
register int k = (((int *) &t)[0] & ncoulmask)>>ncoulshiftbits;
register table_lookup_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];
if (ni < 0) {
force_coul = qiqj*(ftable[k]+f*dftable[k]);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]);
}
else { // correct for special
t = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t);
t.f = (1.0-special_coul[ni])*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t.f);
if (eflag) ecoul = qiqj*(etable[k]+f*detable[k]-t.f);
}
}
}
@ -955,39 +957,37 @@ void PairLJCoul::init_tables()
dptable = (double *) memory->smalloc(ntable*sizeof(double),"pair:dptable");
}
float rsq;
int *int_rsq = (int *) &rsq;
float minrsq;
int *int_minrsq = (int *) &minrsq;
table_lookup_t rsq_lookup;
table_lookup_t minrsq_lookup;
int itablemin;
*int_minrsq = 0 << ncoulshiftbits;
*int_minrsq = *int_minrsq | maskhi;
minrsq_lookup.i = 0 << ncoulshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tabinnersq) {
*int_rsq = i << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tabinnersq) {
rsq_lookup.i = i << ncoulshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrt(rsq);
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
if (cut_respa == NULL) {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
ctable[i] = qqrd2e/r;
etable[i] = qqrd2e/r * derfc;
} else {
rtable[i] = rsq;
rtable[i] = rsq_lookup.f;
ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0);
ctable[i] = 0.0;
etable[i] = qqrd2e/r * derfc;
ptable[i] = qqrd2e/r;
vtable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
ftable[i] += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
ctable[i] = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -997,10 +997,10 @@ void PairLJCoul::init_tables()
}
}
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tabinnersq = minrsq;
tabinnersq = minrsq_lookup.f;
int ntablem1 = ntable - 1;
@ -1038,16 +1038,16 @@ void PairLJCoul::init_tables()
// if so, compute deltas between rsq and cut*cut
double f_tmp,c_tmp,e_tmp,p_tmp,v_tmp;
itablemin = *int_minrsq & ncoulmask;
itablemin = minrsq_lookup.i & ncoulmask;
itablemin >>= ncoulshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
*int_rsq = itablemax << ncoulshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = itablemax << ncoulshiftbits;
rsq_lookup.i |= maskhi;
if (rsq < cut_coulsq) {
rsq = cut_coulsq;
r = sqrt(rsq);
if (rsq_lookup.f < cut_coulsq) {
rsq_lookup.f = cut_coulsq;
r = sqrtf(rsq_lookup.f);
grij = g_ewald * r;
expm2 = exp(-grij*grij);
derfc = erfc(grij);
@ -1062,8 +1062,8 @@ void PairLJCoul::init_tables()
e_tmp = qqrd2e/r * derfc;
p_tmp = qqrd2e/r;
v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2);
if (rsq > cut_respa[2]*cut_respa[2]) {
if (rsq < cut_respa[3]*cut_respa[3]) {
if (rsq_lookup.f > cut_respa[2]*cut_respa[2]) {
if (rsq_lookup.f < cut_respa[3]*cut_respa[3]) {
rsw = (r - cut_respa[2])/(cut_respa[3] - cut_respa[2]);
f_tmp += qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
c_tmp = qqrd2e/r * rsw*rsw*(3.0 - 2.0*rsw);
@ -1074,7 +1074,7 @@ void PairLJCoul::init_tables()
}
}
drtable[itablemax] = 1.0/(rsq - rtable[itablemax]);
drtable[itablemax] = 1.0/(rsq_lookup.f - rtable[itablemax]);
dftable[itablemax] = f_tmp - ftable[itablemax];
dctable[itablemax] = c_tmp - ctable[itablemax];
detable[itablemax] = e_tmp - etable[itablemax];
@ -1126,12 +1126,13 @@ double PairLJCoul::single(int i, int j, int itype, int jtype,
eng += t-r;
}
else { // table real space
register float t = rsq;
register int k = (((int *) &t)[0] & ncoulmask)>>ncoulshiftbits;
register table_lookup_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];
t = (1.0-factor_coul)*(ctable[k]+f*dctable[k]);
force_coul = qiqj*(ftable[k]+f*dftable[k]-t);
eng += qiqj*(etable[k]+f*detable[k]-t);
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);
}
} else force_coul = 0.0;

26
src/pair.cpp
View File

@ -899,8 +899,7 @@ void Pair::write_file(int narg, char **arg)
}
double r,e,f,rsq;
float rsq_float;
int *int_rsq = (int *) &rsq_float;
table_lookup_t rsq_lookup;
for (int i = 0; i < n; i++) {
if (style == R) {
@ -910,13 +909,13 @@ void Pair::write_file(int narg, char **arg)
rsq = inner*inner + (outer*outer - inner*inner) * i/(n-1);
r = sqrt(rsq);
} else if (style == BMP) {
*int_rsq = i << nshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq_float < inner*inner) {
*int_rsq = i << nshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << nshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < inner*inner) {
rsq_lookup.i = i << nshiftbits;
rsq_lookup.i |= maskhi;
}
rsq = rsq_float;
rsq = rsq_lookup.f;
r = sqrt(rsq);
}
@ -981,12 +980,11 @@ void Pair::init_bitmap(double inner, double outer, int ntablebits,
for (int j = 0; j < ntablebits+nshiftbits; j++) nmask *= 2;
nmask -= 1;
float rsq;
int *int_rsq = (int *) &rsq;
rsq = outer*outer;
maskhi = *int_rsq & ~(nmask);
rsq = inner*inner;
masklo = *int_rsq & ~(nmask);
table_lookup_t rsq_lookup;
rsq_lookup.f = outer*outer;
maskhi = rsq_lookup.i & ~(nmask);
rsq_lookup.f = inner*inner;
masklo = rsq_lookup.i & ~(nmask);
}
/* ---------------------------------------------------------------------- */

1
src/pair.h
View File

@ -105,6 +105,7 @@ class Pair : protected Pointers {
int offset_flag,mix_flag; // flags for offset and mixing
int ncoultablebits; // size of Coulomb table
double tabinner; // inner cutoff for Coulomb table
typedef union { int i; float f;} table_lookup_t; // custom data type for accessing Coulomb tables.
double THIRD;

75
src/pair_table.cpp
View File

@ -74,8 +74,8 @@ void PairTable::compute(int eflag, int vflag)
double rsq,factor_lj,fraction,value,a,b;
int *ilist,*jlist,*numneigh,**firstneigh;
Table *tb;
float rsq_single;
int *int_rsq = (int *) &rsq_single;
table_lookup_t rsq_lookup;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
@ -148,10 +148,10 @@ void PairTable::compute(int eflag, int vflag)
tb->deltasq6;
fpair = factor_lj * value;
} else {
rsq_single = rsq;
itable = *int_rsq & tb->nmask;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & tb->nmask;
itable >>= tb->nshiftbits;
fraction = (rsq_single - tb->rsq[itable]) * tb->drsq[itable];
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
value = tb->f[itable] + fraction*tb->df[itable];
fpair = factor_lj * value;
}
@ -394,8 +394,7 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
int itmp;
double rtmp;
float rsq;
int *int_rsq = (int *) &rsq;
table_lookup_t rsq_lookup;
fgets(line,MAXLINE,fp);
for (int i = 0; i < tb->ninput; i++) {
@ -409,13 +408,13 @@ void PairTable::read_table(Table *tb, char *file, char *keyword)
(tb->rhi*tb->rhi - tb->rlo*tb->rlo)*i/(tb->ninput-1);
rtmp = sqrt(rtmp);
} else if (tb->rflag == BMP) {
*int_rsq = i << nshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tb->rlo*tb->rlo) {
*int_rsq = i << nshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << nshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tb->rlo*tb->rlo) {
rsq_lookup.i = i << nshiftbits;
rsq_lookup.i |= maskhi;
}
rtmp = sqrtf(rsq);
rtmp = sqrtf(rsq_lookup.f);
}
tb->rfile[i] = rtmp;
@ -678,8 +677,7 @@ void PairTable::compute_table(Table *tb)
if (tabstyle == BITMAP) {
double r;
float rsq;
int *int_rsq = (int *) &rsq;
table_lookup_t rsq_lookup;
int masklo,maskhi;
// linear lookup tables of length ntable = 2^n
@ -696,20 +694,19 @@ void PairTable::compute_table(Table *tb)
tb->df = (double *) memory->smalloc(ntable*sizeof(double),"pair:df");
tb->drsq = (double *) memory->smalloc(ntable*sizeof(double),"pair:drsq");
float minrsq;
int *int_minrsq = (int *) &minrsq;
*int_minrsq = 0 << tb->nshiftbits;
*int_minrsq = *int_minrsq | maskhi;
table_lookup_t minrsq_lookup;
minrsq_lookup.i = 0 << tb->nshiftbits;
minrsq_lookup.i |= maskhi;
for (int i = 0; i < ntable; i++) {
*int_rsq = i << tb->nshiftbits;
*int_rsq = *int_rsq | masklo;
if (rsq < tb->innersq) {
*int_rsq = i << tb->nshiftbits;
*int_rsq = *int_rsq | maskhi;
rsq_lookup.i = i << tb->nshiftbits;
rsq_lookup.i |= masklo;
if (rsq_lookup.f < tb->innersq) {
rsq_lookup.i = i << tb->nshiftbits;
rsq_lookup.i |= maskhi;
}
r = sqrtf(rsq);
tb->rsq[i] = rsq;
r = sqrtf(rsq_lookup.f);
tb->rsq[i] = rsq_lookup.f;
if (tb->match) {
tb->e[i] = tb->efile[i];
tb->f[i] = tb->ffile[i]/r;
@ -717,10 +714,10 @@ void PairTable::compute_table(Table *tb)
tb->e[i] = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r);
tb->f[i] = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,r)/r;
}
minrsq = MIN(minrsq,rsq);
minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f);
}
tb->innersq = minrsq;
tb->innersq = minrsq_lookup.f;
for (int i = 0; i < ntablem1; i++) {
tb->de[i] = tb->e[i+1] - tb->e[i];
@ -746,27 +743,27 @@ void PairTable::compute_table(Table *tb)
// deltas at itablemax-1 as a good approximation
double e_tmp,f_tmp;
int itablemin = *int_minrsq & tb->nmask;
int itablemin = minrsq_lookup.i & tb->nmask;
itablemin >>= tb->nshiftbits;
int itablemax = itablemin - 1;
if (itablemin == 0) itablemax = ntablem1;
int itablemaxm1 = itablemax - 1;
if (itablemax == 0) itablemaxm1 = ntablem1;
*int_rsq = itablemax << tb->nshiftbits;
*int_rsq = *int_rsq | maskhi;
if (rsq < tb->cut*tb->cut) {
rsq_lookup.i = itablemax << tb->nshiftbits;
rsq_lookup.i |= maskhi;
if (rsq_lookup.f < tb->cut*tb->cut) {
if (tb->match) {
tb->de[itablemax] = tb->de[itablemaxm1];
tb->df[itablemax] = tb->df[itablemaxm1];
tb->drsq[itablemax] = tb->drsq[itablemaxm1];
} else {
rsq = tb->cut*tb->cut;
r = sqrtf(rsq);
rsq_lookup.f = tb->cut*tb->cut;
r = sqrtf(rsq_lookup.f);
e_tmp = splint(tb->rfile,tb->efile,tb->e2file,tb->ninput,r);
f_tmp = splint(tb->rfile,tb->ffile,tb->f2file,tb->ninput,r)/r;
tb->de[itablemax] = e_tmp - tb->e[itablemax];
tb->df[itablemax] = f_tmp - tb->f[itablemax];
tb->drsq[itablemax] = 1.0/(rsq - tb->rsq[itablemax]);
tb->drsq[itablemax] = 1.0/(rsq_lookup.f - tb->rsq[itablemax]);
}
}
}
@ -938,11 +935,11 @@ double PairTable::single(int i, int j, int itype, int jtype, double rsq,
tb->deltasq6;
fforce = factor_lj * value;
} else {
float rsq_single = rsq;
int *int_rsq = (int *) &rsq_single;
itable = *int_rsq & tb->nmask;
table_lookup_t rsq_lookup;
rsq_lookup.f = rsq;
itable = rsq_lookup.i & tb->nmask;
itable >>= tb->nshiftbits;
fraction = (rsq_single - tb->rsq[itable]) * tb->drsq[itable];
fraction = (rsq_lookup.f - tb->rsq[itable]) * tb->drsq[itable];
value = tb->f[itable] + fraction*tb->df[itable];
fforce = factor_lj * value;
}