diff --git a/src/KSPACE/pair_born_coul_long.cpp b/src/KSPACE/pair_born_coul_long.cpp index 431d7bcf1b..5260e5f7e0 100644 --- a/src/KSPACE/pair_born_coul_long.cpp +++ b/src/KSPACE/pair_born_coul_long.cpp @@ -46,6 +46,7 @@ using namespace MathConst; PairBornCoulLong::PairBornCoulLong(LAMMPS *lmp) : Pair(lmp) { ewaldflag = pppmflag = 1; + ftable = NULL; } /* ---------------------------------------------------------------------- */ @@ -69,14 +70,16 @@ PairBornCoulLong::~PairBornCoulLong() memory->destroy(born3); memory->destroy(offset); } + if (ftable) free_tables(); } /* ---------------------------------------------------------------------- */ void PairBornCoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itable,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; double rsq,r2inv,r6inv,forcecoul,forceborn,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double r,rexp; @@ -127,16 +130,31 @@ void PairBornCoulLong::compute(int eflag, int vflag) if (rsq < cutsq[itype][jtype]) { r2inv = 1.0/rsq; - r = sqrt(rsq); if (rsq < cut_coulsq) { - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } } else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -159,7 +177,12 @@ void PairBornCoulLong::compute(int eflag, int vflag) if (eflag) { if (rsq < cut_coulsq) { - ecoul = prefactor*erfc; + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -356,6 +379,9 @@ void PairBornCoulLong::init_style() g_ewald = force->kspace->g_ewald; neighbor->request(this); + + // setup force tables + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -427,6 +453,8 @@ void PairBornCoulLong::write_restart_settings(FILE *fp) fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&ncoultablebits,sizeof(int),1,fp); + fwrite(&tabinner,sizeof(double),1,fp); } /* ---------------------------------------------------------------------- @@ -441,12 +469,16 @@ void PairBornCoulLong::read_restart_settings(FILE *fp) fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); fread(&tail_flag,sizeof(int),1,fp); + fread(&ncoultablebits,sizeof(int),1,fp); + fread(&tabinner,sizeof(double),1,fp); } MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); + MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ @@ -457,18 +489,34 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype, double &fforce) { double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor; - double forcecoul,forceborn,phicoul,phiborn; + double fraction,table,forcecoul,forceborn,phicoul,phiborn; + int itable; r2inv = 1.0/rsq; if (rsq < cut_coulsq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + 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) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } } else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { r6inv = r2inv*r2inv*r2inv; @@ -481,7 +529,12 @@ double PairBornCoulLong::single(int i, int j, int itype, int jtype, double eng = 0.0; if (rsq < cut_coulsq) { - phicoul = prefactor*erfc; + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; eng += phicoul; } diff --git a/src/KSPACE/pair_born_coul_long.h b/src/KSPACE/pair_born_coul_long.h index dc60fb13f6..0cea32d0b9 100644 --- a/src/KSPACE/pair_born_coul_long.h +++ b/src/KSPACE/pair_born_coul_long.h @@ -46,6 +46,7 @@ class PairBornCoulLong : public Pair { double cut_coul,cut_coulsq; double **a,**rho,**sigma,**c,**d; double **rhoinv,**born1,**born2,**born3,**offset; + double *cut_respa; double g_ewald; void allocate(); diff --git a/src/KSPACE/pair_born_coul_msm.cpp b/src/KSPACE/pair_born_coul_msm.cpp index 3590265d8b..a5c977fe43 100644 --- a/src/KSPACE/pair_born_coul_msm.cpp +++ b/src/KSPACE/pair_born_coul_msm.cpp @@ -99,7 +99,7 @@ void PairBornCoulMSM::compute(int eflag, int vflag) r2inv = 1.0/rsq; if (rsq < cut_coulsq) { - r = sqrt(rsq); + r = sqrt(rsq); prefactor = qqrd2e * qtmp*q[j]/r; egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); fgamma = 1.0 + (rsq/cut_coulsq)*force->kspace->dgamma(r/cut_coul); diff --git a/src/KSPACE/pair_buck_coul_long.cpp b/src/KSPACE/pair_buck_coul_long.cpp index 9e5f8660ff..6f0f223701 100644 --- a/src/KSPACE/pair_buck_coul_long.cpp +++ b/src/KSPACE/pair_buck_coul_long.cpp @@ -42,6 +42,7 @@ using namespace MathConst; PairBuckCoulLong::PairBuckCoulLong(LAMMPS *lmp) : Pair(lmp) { ewaldflag = pppmflag = 1; + ftable = NULL; } /* ---------------------------------------------------------------------- */ @@ -62,14 +63,16 @@ PairBuckCoulLong::~PairBuckCoulLong() memory->destroy(buck2); memory->destroy(offset); } + if (ftable) free_tables(); } /* ---------------------------------------------------------------------- */ void PairBuckCoulLong::compute(int eflag, int vflag) { - int i,j,ii,jj,inum,jnum,itype,jtype; + int i,j,ii,jj,inum,jnum,itable,itype,jtype; double qtmp,xtmp,ytmp,ztmp,delx,dely,delz,evdwl,ecoul,fpair; + double fraction,table; double rsq,r2inv,r6inv,forcecoul,forcebuck,factor_coul,factor_lj; double grij,expm2,prefactor,t,erfc; double r,rexp; @@ -119,17 +122,31 @@ void PairBuckCoulLong::compute(int eflag, int vflag) jtype = type[j]; if (rsq < cutsq[itype][jtype]) { - r2inv = 1.0/rsq; - r = sqrt(rsq); - + r2inv = 1.0/rsq; if (rsq < cut_coulsq) { - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = qqrd2e * qtmp*q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = qqrd2e * qtmp*q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + fraction = (rsq_lookup.f - rtable[itable]) * drtable[itable]; + table = ftable[itable] + fraction*dftable[itable]; + forcecoul = qtmp*q[j] * table; + if (factor_coul < 1.0) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = qtmp*q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } } else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -151,7 +168,12 @@ void PairBuckCoulLong::compute(int eflag, int vflag) if (eflag) { if (rsq < cut_coulsq) { - ecoul = prefactor*erfc; + if (!ncoultablebits || rsq <= tabinnersq) + ecoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + ecoul = qtmp*q[j] * table; + } if (factor_coul < 1.0) ecoul -= (1.0-factor_coul)*prefactor; } else ecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { @@ -334,6 +356,9 @@ void PairBuckCoulLong::init_style() g_ewald = force->kspace->g_ewald; neighbor->request(this); + + // setup force tables + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -399,6 +424,8 @@ void PairBuckCoulLong::write_restart_settings(FILE *fp) fwrite(&offset_flag,sizeof(int),1,fp); fwrite(&mix_flag,sizeof(int),1,fp); fwrite(&tail_flag,sizeof(int),1,fp); + fwrite(&ncoultablebits,sizeof(int),1,fp); + fwrite(&tabinner,sizeof(double),1,fp); } /* ---------------------------------------------------------------------- @@ -413,12 +440,16 @@ void PairBuckCoulLong::read_restart_settings(FILE *fp) fread(&offset_flag,sizeof(int),1,fp); fread(&mix_flag,sizeof(int),1,fp); fread(&tail_flag,sizeof(int),1,fp); + fread(&ncoultablebits,sizeof(int),1,fp); + fread(&tabinner,sizeof(double),1,fp); } MPI_Bcast(&cut_lj_global,1,MPI_DOUBLE,0,world); MPI_Bcast(&cut_coul,1,MPI_DOUBLE,0,world); MPI_Bcast(&offset_flag,1,MPI_INT,0,world); MPI_Bcast(&mix_flag,1,MPI_INT,0,world); MPI_Bcast(&tail_flag,1,MPI_INT,0,world); + MPI_Bcast(&ncoultablebits,1,MPI_INT,0,world); + MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); } /* ---------------------------------------------------------------------- */ @@ -429,18 +460,34 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype, double &fforce) { double r2inv,r6inv,r,rexp,grij,expm2,t,erfc,prefactor; - double forcecoul,forcebuck,phicoul,phibuck; + double fraction,table,forcecoul,forcebuck,phicoul,phibuck; + int itable; r2inv = 1.0/rsq; if (rsq < cut_coulsq) { - r = sqrt(rsq); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - t = 1.0 / (1.0 + EWALD_P*grij); - erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; - prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; - forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); - if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + if (!ncoultablebits || rsq <= tabinnersq) { + r = sqrt(rsq); + grij = g_ewald * r; + expm2 = exp(-grij*grij); + t = 1.0 / (1.0 + EWALD_P*grij); + erfc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * expm2; + prefactor = force->qqrd2e * atom->q[i]*atom->q[j]/r; + forcecoul = prefactor * (erfc + EWALD_F*grij*expm2); + if (factor_coul < 1.0) forcecoul -= (1.0-factor_coul)*prefactor; + } else { + union_int_float_t rsq_lookup; + rsq_lookup.f = rsq; + itable = rsq_lookup.i & ncoulmask; + itable >>= ncoulshiftbits; + 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) { + table = ctable[itable] + fraction*dctable[itable]; + prefactor = atom->q[i]*atom->q[j] * table; + forcecoul -= (1.0-factor_coul)*prefactor; + } + } } else forcecoul = 0.0; if (rsq < cut_ljsq[itype][jtype]) { r6inv = r2inv*r2inv*r2inv; @@ -452,7 +499,12 @@ double PairBuckCoulLong::single(int i, int j, int itype, int jtype, double eng = 0.0; if (rsq < cut_coulsq) { - phicoul = prefactor*erfc; + if (!ncoultablebits || rsq <= tabinnersq) + phicoul = prefactor*erfc; + else { + table = etable[itable] + fraction*detable[itable]; + phicoul = atom->q[i]*atom->q[j] * table; + } if (factor_coul < 1.0) phicoul -= (1.0-factor_coul)*prefactor; eng += phicoul; } diff --git a/src/KSPACE/pair_buck_coul_long.h b/src/KSPACE/pair_buck_coul_long.h index 05903b6870..3dbcf224f5 100644 --- a/src/KSPACE/pair_buck_coul_long.h +++ b/src/KSPACE/pair_buck_coul_long.h @@ -46,6 +46,8 @@ class PairBuckCoulLong : public Pair { double cut_coul,cut_coulsq; double **a,**rho,**c; double **rhoinv,**buck1,**buck2,**offset; + + double *cut_respa; double g_ewald; void allocate(); diff --git a/src/KSPACE/pair_buck_long_coul_long.cpp b/src/KSPACE/pair_buck_long_coul_long.cpp index c354b3d28f..4506c989cc 100644 --- a/src/KSPACE/pair_buck_long_coul_long.cpp +++ b/src/KSPACE/pair_buck_long_coul_long.cpp @@ -290,7 +290,7 @@ void PairBuckLongCoulLong::init_style() // setup force tables - if (ncoultablebits) init_tables(); + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -929,195 +929,6 @@ void PairBuckLongCoulLong::compute_outer(int eflag, int vflag) } } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairBuckLongCoulLong::init_tables() -{ - int masklo,maskhi; - double r,grij,expm2,derfc,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - 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_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_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); - } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = 0.0,v_tmp = 0.0; - itablemin = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - 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); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * derfc; - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * derfc; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - 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); - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairBuckLongCoulLong::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - /* ---------------------------------------------------------------------- */ double PairBuckLongCoulLong::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_buck_long_coul_long.h b/src/KSPACE/pair_buck_long_coul_long.h index 7c36af1761..62858b487f 100644 --- a/src/KSPACE/pair_buck_long_coul_long.h +++ b/src/KSPACE/pair_buck_long_coul_long.h @@ -60,15 +60,8 @@ class PairBuckLongCoulLong : public Pair { double g_ewald_6; int ewald_order, ewald_off; - double tabinnersq; - double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable; - double *etable, *detable, *ptable, *dptable, *vtable, *dvtable; - int ncoulshiftbits, ncoulmask; - void options(char **arg, int order); void allocate(); - void init_tables(); - void free_tables(); }; } diff --git a/src/KSPACE/pair_coul_long.cpp b/src/KSPACE/pair_coul_long.cpp index 1ec2544703..7d2e6d8a95 100644 --- a/src/KSPACE/pair_coul_long.cpp +++ b/src/KSPACE/pair_coul_long.cpp @@ -257,7 +257,7 @@ void PairCoulLong::init_style() // setup force tables - if (ncoultablebits) init_tables(); + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -271,175 +271,6 @@ double PairCoulLong::init_one(int i, int j) return cut_coul; } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairCoulLong::init_tables() -{ - int masklo,maskhi; - double r,grij,expm2,derfc,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - 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_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_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); - } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - 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); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * derfc; - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * derfc; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - 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); - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -493,26 +324,6 @@ void PairCoulLong::read_restart_settings(FILE *fp) MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); } -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairCoulLong::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - /* ---------------------------------------------------------------------- */ double PairCoulLong::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_coul_long.h b/src/KSPACE/pair_coul_long.h index 70d97c1897..f899706a37 100644 --- a/src/KSPACE/pair_coul_long.h +++ b/src/KSPACE/pair_coul_long.h @@ -46,14 +46,7 @@ class PairCoulLong : public Pair { double g_ewald; double **scale; - double tabinnersq; - double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable; - double *etable,*detable,*ptable,*dptable,*vtable,*dvtable; - int ncoulshiftbits,ncoulmask; - void allocate(); - virtual void init_tables(); - void free_tables(); }; } diff --git a/src/KSPACE/pair_coul_msm.cpp b/src/KSPACE/pair_coul_msm.cpp index 8fc9ef1701..7e7ac3ecc6 100644 --- a/src/KSPACE/pair_coul_msm.cpp +++ b/src/KSPACE/pair_coul_msm.cpp @@ -94,7 +94,7 @@ void PairCoulMSM::compute(int eflag, int vflag) delz = ztmp - x[j][2]; rsq = delx*delx + dely*dely + delz*delz; jtype = type[j]; - + if (rsq < cut_coulsq) { r2inv = 1.0/rsq; if (!ncoultablebits || rsq <= tabinnersq) { @@ -149,173 +149,6 @@ void PairCoulMSM::compute(int eflag, int vflag) if (vflag_fdotr) virial_fdotr_compute(); } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairCoulMSM::init_tables() -{ - int masklo,maskhi; - double r,egamma,fgamma,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); - fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); - if (cut_respa == NULL) { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * fgamma; - ctable[i] = qqrd2e/r; - etable[i] = qqrd2e/r * egamma; - } else { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (fgamma - 1.0); - ctable[i] = 0.0; - etable[i] = qqrd2e/r * egamma; - ptable[i] = qqrd2e/r; - vtable[i] = qqrd2e/r * fgamma; - 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); - } else { - ftable[i] = qqrd2e/r * fgamma; - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - if (rsq_lookup.f < cut_coulsq) { - rsq_lookup.f = cut_coulsq; - r = sqrtf(rsq_lookup.f); - egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); - fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * fgamma; - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * egamma; - } else { - f_tmp = qqrd2e/r * (fgamma - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * egamma; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * fgamma; - 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); - } else { - f_tmp = qqrd2e/r * fgamma; - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - /* ---------------------------------------------------------------------- */ double PairCoulMSM::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_coul_msm.h b/src/KSPACE/pair_coul_msm.h index 2a9bc340fd..1dad2b4852 100644 --- a/src/KSPACE/pair_coul_msm.h +++ b/src/KSPACE/pair_coul_msm.h @@ -31,9 +31,6 @@ class PairCoulMSM : public PairCoulLong { virtual void compute(int, int); virtual double single(int, int, int, int, double, double, double, double &); virtual void *extract(const char *, int &); - - protected: - virtual void init_tables(); }; } diff --git a/src/KSPACE/pair_lj_charmm_coul_long.cpp b/src/KSPACE/pair_lj_charmm_coul_long.cpp index c11ab3bd21..8170483ecf 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_long.cpp @@ -775,7 +775,7 @@ void PairLJCharmmCoulLong::init_style() // setup force tables - if (ncoultablebits) init_tables(); + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -829,175 +829,6 @@ double PairLJCharmmCoulLong::init_one(int i, int j) return cut; } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJCharmmCoulLong::init_tables() -{ - int masklo,maskhi; - double r,grij,expm2,derfc,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - 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_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_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); - } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - 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); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * derfc; - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * derfc; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - 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); - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -1089,26 +920,6 @@ void PairLJCharmmCoulLong::read_restart_settings(FILE *fp) MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); } -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairLJCharmmCoulLong::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - /* ---------------------------------------------------------------------- */ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_lj_charmm_coul_long.h b/src/KSPACE/pair_lj_charmm_coul_long.h index 4c28126bf2..889f632833 100644 --- a/src/KSPACE/pair_lj_charmm_coul_long.h +++ b/src/KSPACE/pair_lj_charmm_coul_long.h @@ -59,14 +59,7 @@ class PairLJCharmmCoulLong : public Pair { double *cut_respa; double g_ewald; - double tabinnersq; - double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable; - double *etable,*detable,*ptable,*dptable,*vtable,*dvtable; - int ncoulshiftbits,ncoulmask; - void allocate(); - virtual void init_tables(); - void free_tables(); }; } diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.cpp b/src/KSPACE/pair_lj_charmm_coul_msm.cpp index 7e777605eb..fb81550a0b 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.cpp +++ b/src/KSPACE/pair_lj_charmm_coul_msm.cpp @@ -390,174 +390,6 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag) } } - -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJCharmmCoulMSM::init_tables() -{ - int masklo,maskhi; - double r,egamma,fgamma,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); - fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); - if (cut_respa == NULL) { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * fgamma; - ctable[i] = qqrd2e/r; - etable[i] = qqrd2e/r * egamma; - } else { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (fgamma - 1.0); - ctable[i] = 0.0; - etable[i] = qqrd2e/r * egamma; - ptable[i] = qqrd2e/r; - vtable[i] = qqrd2e/r * fgamma; - 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); - } else { - ftable[i] = qqrd2e/r * fgamma; - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - if (rsq_lookup.f < cut_coulsq) { - rsq_lookup.f = cut_coulsq; - r = sqrtf(rsq_lookup.f); - egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); - fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * fgamma; - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * egamma; - } else { - f_tmp = qqrd2e/r * (fgamma - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * egamma; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * fgamma; - 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); - } else { - f_tmp = qqrd2e/r * fgamma; - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - /* ---------------------------------------------------------------------- */ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_lj_charmm_coul_msm.h b/src/KSPACE/pair_lj_charmm_coul_msm.h index aa437ef51e..2efe34777b 100644 --- a/src/KSPACE/pair_lj_charmm_coul_msm.h +++ b/src/KSPACE/pair_lj_charmm_coul_msm.h @@ -32,9 +32,6 @@ class PairLJCharmmCoulMSM : public PairLJCharmmCoulLong { virtual double single(int, int, int, int, double, double, double, double &); virtual void compute_outer(int, int); virtual void *extract(const char *, int &); - - protected: - virtual void init_tables(); }; } diff --git a/src/KSPACE/pair_lj_cut_coul_long.cpp b/src/KSPACE/pair_lj_cut_coul_long.cpp index ae4c7c069b..97772aadd4 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.cpp +++ b/src/KSPACE/pair_lj_cut_coul_long.cpp @@ -707,7 +707,7 @@ void PairLJCutCoulLong::init_style() // setup force tables - if (ncoultablebits) init_tables(); + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -792,175 +792,6 @@ double PairLJCutCoulLong::init_one(int i, int j) return cut; } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::init_tables() -{ - int masklo,maskhi; - double r,grij,expm2,derfc,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - 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_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_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); - } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - 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); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * derfc; - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * derfc; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - 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); - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -1049,26 +880,6 @@ void PairLJCutCoulLong::read_restart_settings(FILE *fp) MPI_Bcast(&tabinner,1,MPI_DOUBLE,0,world); } -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairLJCutCoulLong::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - /* ---------------------------------------------------------------------- */ double PairLJCutCoulLong::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_lj_cut_coul_long.h b/src/KSPACE/pair_lj_cut_coul_long.h index 6f65573a50..8204823bb1 100644 --- a/src/KSPACE/pair_lj_cut_coul_long.h +++ b/src/KSPACE/pair_lj_cut_coul_long.h @@ -56,14 +56,7 @@ class PairLJCutCoulLong : public Pair { double qdist; // TIP4P distance from O site to negative charge double g_ewald; - double tabinnersq; - double *rtable,*drtable,*ftable,*dftable,*ctable,*dctable; - double *etable,*detable,*ptable,*dptable,*vtable,*dvtable; - int ncoulshiftbits,ncoulmask; - void allocate(); - virtual void init_tables(); - void free_tables(); }; } diff --git a/src/KSPACE/pair_lj_cut_coul_msm.cpp b/src/KSPACE/pair_lj_cut_coul_msm.cpp index b39eabf4c6..262184f9f3 100644 --- a/src/KSPACE/pair_lj_cut_coul_msm.cpp +++ b/src/KSPACE/pair_lj_cut_coul_msm.cpp @@ -347,173 +347,6 @@ void PairLJCutCoulMSM::compute_outer(int eflag, int vflag) } } } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJCutCoulMSM::init_tables() -{ - int masklo,maskhi; - double r,egamma,fgamma,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); - fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); - - if (cut_respa == NULL) { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * fgamma; - ctable[i] = qqrd2e/r; - etable[i] = qqrd2e/r * egamma; - } else { - rtable[i] = rsq_lookup.f; - ftable[i] = qqrd2e/r * (fgamma - 1.0); - ctable[i] = 0.0; - etable[i] = qqrd2e/r * egamma; - ptable[i] = qqrd2e/r; - vtable[i] = qqrd2e/r * fgamma; - 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); - } else { - ftable[i] = qqrd2e/r * fgamma; - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - if (rsq_lookup.f < cut_coulsq) { - rsq_lookup.f = cut_coulsq; - r = sqrtf(rsq_lookup.f); - egamma = 1.0 - (r/cut_coul)*force->kspace->gamma(r/cut_coul); - fgamma = 1.0 + (rsq_lookup.f/cut_coulsq)*force->kspace->dgamma(r/cut_coul); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * fgamma; - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * egamma; - } else { - f_tmp = qqrd2e/r * (fgamma - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * egamma; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * fgamma; - 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); - } else { - f_tmp = qqrd2e/r * fgamma; - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} /* ---------------------------------------------------------------------- */ diff --git a/src/KSPACE/pair_lj_cut_coul_msm.h b/src/KSPACE/pair_lj_cut_coul_msm.h index 6e69c3f727..3c88e22896 100644 --- a/src/KSPACE/pair_lj_cut_coul_msm.h +++ b/src/KSPACE/pair_lj_cut_coul_msm.h @@ -32,9 +32,6 @@ class PairLJCutCoulMSM : public PairLJCutCoulLong { virtual double single(int, int, int, int, double, double, double, double &); virtual void compute_outer(int, int); virtual void *extract(const char *, int &); - - protected: - virtual void init_tables(); }; } diff --git a/src/KSPACE/pair_lj_long_coul_long.cpp b/src/KSPACE/pair_lj_long_coul_long.cpp index 61c068f7f0..8f7eb79c8d 100644 --- a/src/KSPACE/pair_lj_long_coul_long.cpp +++ b/src/KSPACE/pair_lj_long_coul_long.cpp @@ -286,7 +286,7 @@ void PairLJLongCoulLong::init_style() // setup force tables - if (ncoultablebits) init_tables(); + if (ncoultablebits) init_tables(cut_coul,cut_respa); } /* ---------------------------------------------------------------------- @@ -919,195 +919,6 @@ void PairLJLongCoulLong::compute_outer(int eflag, int vflag) } } -/* ---------------------------------------------------------------------- - setup force tables used in compute routines -------------------------------------------------------------------------- */ - -void PairLJLongCoulLong::init_tables() -{ - int masklo,maskhi; - double r,grij,expm2,derfc,rsw; - double qqrd2e = force->qqrd2e; - - tabinnersq = tabinner*tabinner; - init_bitmap(tabinner,cut_coul,ncoultablebits, - masklo,maskhi,ncoulmask,ncoulshiftbits); - - int ntable = 1; - for (int i = 0; i < ncoultablebits; i++) ntable *= 2; - - // linear lookup tables of length N = 2^ncoultablebits - // stored value = value at lower edge of bin - // d values = delta from lower edge to upper edge of bin - - if (ftable) free_tables(); - - memory->create(rtable,ntable,"pair:rtable"); - memory->create(ftable,ntable,"pair:ftable"); - memory->create(ctable,ntable,"pair:ctable"); - memory->create(etable,ntable,"pair:etable"); - memory->create(drtable,ntable,"pair:drtable"); - memory->create(dftable,ntable,"pair:dftable"); - memory->create(dctable,ntable,"pair:dctable"); - memory->create(detable,ntable,"pair:detable"); - - if (cut_respa == NULL) { - vtable = ptable = dvtable = dptable = NULL; - } else { - memory->create(vtable,ntable,"pair:vtable"); - memory->create(ptable,ntable,"pair:ptable"); - memory->create(dvtable,ntable,"pair:dvtable"); - memory->create(dptable,ntable,"pair:dptable"); - } - - union_int_float_t rsq_lookup; - union_int_float_t minrsq_lookup; - int itablemin; - minrsq_lookup.i = 0 << ncoulshiftbits; - minrsq_lookup.i |= maskhi; - - for (int i = 0; i < ntable; i++) { - 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_lookup.f); - grij = g_ewald * r; - expm2 = exp(-grij*grij); - derfc = erfc(grij); - if (cut_respa == NULL) { - 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_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_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); - } else { - ftable[i] = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - ctable[i] = qqrd2e/r; - } - } - } - minrsq_lookup.f = MIN(minrsq_lookup.f,rsq_lookup.f); - } - - tabinnersq = minrsq_lookup.f; - - int ntablem1 = ntable - 1; - - for (int i = 0; i < ntablem1; i++) { - drtable[i] = 1.0/(rtable[i+1] - rtable[i]); - dftable[i] = ftable[i+1] - ftable[i]; - dctable[i] = ctable[i+1] - ctable[i]; - detable[i] = etable[i+1] - etable[i]; - } - - if (cut_respa) { - for (int i = 0; i < ntablem1; i++) { - dvtable[i] = vtable[i+1] - vtable[i]; - dptable[i] = ptable[i+1] - ptable[i]; - } - } - - // get the delta values for the last table entries - // tables are connected periodically between 0 and ntablem1 - - drtable[ntablem1] = 1.0/(rtable[0] - rtable[ntablem1]); - dftable[ntablem1] = ftable[0] - ftable[ntablem1]; - dctable[ntablem1] = ctable[0] - ctable[ntablem1]; - detable[ntablem1] = etable[0] - etable[ntablem1]; - if (cut_respa) { - dvtable[ntablem1] = vtable[0] - vtable[ntablem1]; - dptable[ntablem1] = ptable[0] - ptable[ntablem1]; - } - - // get the correct delta values at itablemax - // smallest r is in bin itablemin - // largest r is in bin itablemax, which is itablemin-1, - // 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 = 0.0,v_tmp = 0.0; - itablemin = minrsq_lookup.i & ncoulmask; - itablemin >>= ncoulshiftbits; - int itablemax = itablemin - 1; - if (itablemin == 0) itablemax = ntablem1; - rsq_lookup.i = itablemax << ncoulshiftbits; - rsq_lookup.i |= maskhi; - - 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); - - if (cut_respa == NULL) { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - e_tmp = qqrd2e/r * derfc; - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2 - 1.0); - c_tmp = 0.0; - e_tmp = qqrd2e/r * derfc; - p_tmp = qqrd2e/r; - v_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - 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); - } else { - f_tmp = qqrd2e/r * (derfc + EWALD_F*grij*expm2); - c_tmp = qqrd2e/r; - } - } - } - - 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]; - if (cut_respa) { - dvtable[itablemax] = v_tmp - vtable[itablemax]; - dptable[itablemax] = p_tmp - ptable[itablemax]; - } - } -} - -/* ---------------------------------------------------------------------- - free memory for tables used in pair computations -------------------------------------------------------------------------- */ - -void PairLJLongCoulLong::free_tables() -{ - memory->destroy(rtable); - memory->destroy(drtable); - memory->destroy(ftable); - memory->destroy(dftable); - memory->destroy(ctable); - memory->destroy(dctable); - memory->destroy(etable); - memory->destroy(detable); - memory->destroy(vtable); - memory->destroy(dvtable); - memory->destroy(ptable); - memory->destroy(dptable); -} - /* ---------------------------------------------------------------------- */ double PairLJLongCoulLong::single(int i, int j, int itype, int jtype, diff --git a/src/KSPACE/pair_lj_long_coul_long.h b/src/KSPACE/pair_lj_long_coul_long.h index 4d35f078b7..a8304ca21c 100644 --- a/src/KSPACE/pair_lj_long_coul_long.h +++ b/src/KSPACE/pair_lj_long_coul_long.h @@ -60,15 +60,8 @@ class PairLJLongCoulLong : public Pair { double g_ewald_6; int ewald_order, ewald_off; - double tabinnersq; - double *rtable, *drtable, *ftable, *dftable, *ctable, *dctable; - double *etable, *detable, *ptable, *dptable, *vtable, *dvtable; - int ncoulshiftbits, ncoulmask; - void options(char **arg, int order); void allocate(); - void init_tables(); - void free_tables(); }; }