diff --git a/src/AMOEBA/amoeba_convolution.cpp b/src/AMOEBA/amoeba_convolution.cpp index 4759a047df..bb99f48422 100644 --- a/src/AMOEBA/amoeba_convolution.cpp +++ b/src/AMOEBA/amoeba_convolution.cpp @@ -72,9 +72,7 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair, flag3d = 1; if (which == POLAR_GRIDC || which == INDUCE_GRIDC) flag3d = 0; - // NOTE: worry about overflow - - nfft_global = nx * ny * nz; + nfft_global = (bingint) nx * ny * nz; // global indices of grid range from 0 to N-1 // nlo_in,nhi_in = lower/upper limits of the 3d sub-brick of @@ -111,7 +109,7 @@ AmoebaConvolution::AmoebaConvolution(LAMMPS *lmp, Pair *pair, // dist[3] = particle position bound = subbox + skin/2.0 // convert to triclinic if necessary // nlo_out,nhi_out = nlo,nhi + stencil size for particle mapping - // NOTE: this needs to be computed same as IGRID in amoeba + // NOTE: this needs to be computed same as IGRID in PairAmoeba double *prd,*boxlo,*sublo,*subhi; int triclinic = domain->triclinic; diff --git a/src/AMOEBA/amoeba_convolution.h b/src/AMOEBA/amoeba_convolution.h index d15fe8d198..270a501a71 100644 --- a/src/AMOEBA/amoeba_convolution.h +++ b/src/AMOEBA/amoeba_convolution.h @@ -33,11 +33,11 @@ class AmoebaConvolution : protected Pointers { public: int nx,ny,nz; int order; - int nfft_global; // nx * ny * nz int nfft_owned; // owned grid points in FFT decomp int nxlo_in,nxhi_in,nylo_in,nyhi_in,nzlo_in,nzhi_in; int nxlo_out,nxhi_out,nylo_out,nyhi_out,nzlo_out,nzhi_out; int nxlo_fft,nxhi_fft,nylo_fft,nyhi_fft,nzlo_fft,nzhi_fft; + bigint nfft_global; // nx * ny * nz double *grid_brick_start; // lower left corner of (c)grid_brick data AmoebaConvolution(class LAMMPS *, class Pair *, diff --git a/src/AMOEBA/amoeba_dispersion.cpp b/src/AMOEBA/amoeba_dispersion.cpp index 83405cfab9..905a67c92e 100644 --- a/src/AMOEBA/amoeba_dispersion.cpp +++ b/src/AMOEBA/amoeba_dispersion.cpp @@ -245,7 +245,7 @@ void PairAmoeba::dispersion_kspace() int it1,it2,it3; int j0,jgrd0; int k0,kgrd0; - double e,fi,denom; + double e,fi,denom,scale; double r1,r2,r3; double h1,h2,h3; double term,vterm; @@ -317,9 +317,6 @@ void PairAmoeba::dispersion_kspace() fac3 = -2.0*aewald*MY_PI*MY_PI; denom0 = (6.0*volbox)/pow(MY_PI,1.5); - //qgrid[0][0][0][0] = 0.0; // NOTE: why is this needed? - //qgrid[0][0][0][1] = 0.0; - n = 0; for (k = nzlo; k <= nzhi; k++) { for (j = nylo; j <= nyhi; j++) { @@ -355,9 +352,9 @@ void PairAmoeba::dispersion_kspace() virdisp[5] -= h2*h3*vterm; } } else term1 = 0.0; - // NOTE: pre-calc this division only once - gridfft[n] *= -(term1/denom); - gridfft[n+1] *= -(term1/denom); + scale = -term1 / denom; + gridfft[n] *= scale; + gridfft[n+1] *= scale; n += 2; } } diff --git a/src/AMOEBA/amoeba_file.cpp b/src/AMOEBA/amoeba_file.cpp index 34fe2f87a9..431edfd557 100644 --- a/src/AMOEBA/amoeba_file.cpp +++ b/src/AMOEBA/amoeba_file.cpp @@ -107,8 +107,6 @@ void PairAmoeba::read_prmfile(char *filename) if (n < 0) break; MPI_Bcast(line,n+1,MPI_CHAR,0,world); - //printf("Section: %s\n",line); - if (strstr(line,"Force Field") == line) section = FFIELD; else if (strstr(line,"Literature") == line) section = LITERATURE; else if (strstr(line,"Atom Type") == line) section = ATOMTYPE; diff --git a/src/AMOEBA/amoeba_induce.cpp b/src/AMOEBA/amoeba_induce.cpp index a4b7b2e1c6..1d911177ee 100644 --- a/src/AMOEBA/amoeba_induce.cpp +++ b/src/AMOEBA/amoeba_induce.cpp @@ -102,8 +102,8 @@ void PairAmoeba::induce() } // get induced dipoles via the OPT extrapolation method - // NOTE: any way to rewrite these loops to avoid allocating - // uopt,uoptp with a optorder+1 dimension, just optorder ?? + // NOTE: could rewrite these loops to avoid allocating + // uopt,uoptp with a optorder+1 dimension, just optorder // since no need to store optorder+1 values after these loops if (poltyp == OPT) { @@ -332,8 +332,6 @@ void PairAmoeba::induce() } } - // NOTE: comp of b,bp and allreduce only needed if pcgprec ? - reduce[0] = b; reduce[1] = bp; MPI_Allreduce(reduce,allreduce,4,MPI_DOUBLE,MPI_SUM,world); @@ -458,7 +456,6 @@ void PairAmoeba::ulspred() } // derive normal equations corresponding to least squares fit - // NOTE: check all N vs N-1 indices in code from here down } else if (polpred == LSQR) { double ***udalt = fixudalt->tstore; @@ -1201,8 +1198,6 @@ void PairAmoeba::udirect2b(double **field, double **fieldp) numneigh = list->numneigh; firstneigh = list->firstneigh; - // NOTE: does this have a problem if aewald is tiny ?? - aesq2 = 2.0 * aewald * aewald; aesq2n = 0.0; if (aewald > 0.0) aesq2n = 1.0 / (MY_PIS*aewald); diff --git a/src/AMOEBA/amoeba_utils.cpp b/src/AMOEBA/amoeba_utils.cpp index 059587a093..7bb95a26a2 100644 --- a/src/AMOEBA/amoeba_utils.cpp +++ b/src/AMOEBA/amoeba_utils.cpp @@ -34,7 +34,6 @@ enum{NOFRAME,ZONLY,ZTHENX,BISECTOR,ZBISECT,THREEFOLD}; any of the values can be 0 if not used yaxis can later be negative due to chkpole() also sets polaxe and pole[13] multipole for each owned atom - NOTE: may not always do this identically to Tinker b/c atom order matters ------------------------------------------------------------------------- */ void PairAmoeba::kmpole() @@ -71,6 +70,8 @@ void PairAmoeba::kmpole() flag = 0; // create a sorted version of bond/angle neighs from special[][] + // NOTE: this is try and do it identically to Tinker + // b/c I think in Tinker, atom order matters as to which case is seen fist for (j = 0; j < nspecial[i][0]; j++) bondneigh[j] = special[i][j]; diff --git a/src/AMOEBA/angle_amoeba.cpp b/src/AMOEBA/angle_amoeba.cpp index e6f866f63b..c81e22c8f3 100644 --- a/src/AMOEBA/angle_amoeba.cpp +++ b/src/AMOEBA/angle_amoeba.cpp @@ -315,7 +315,7 @@ void AngleAmoeba::tinker_anglep(int i1, int i2, int i3, int type, int eflag) rap2 = xap*xap + yap*yap + zap*zap; rcp2 = xcp*xcp + ycp*ycp + zcp*zcp; - // NOTE: can these be 0.0 ? what to do? + // Tinker just skips the computation in either is zero if (rap2 == 0.0 || rcp2 == 0.0) return; @@ -825,7 +825,10 @@ void AngleAmoeba::write_data(FILE *fp) fprintf(fp,"%d %g %g\n",i,ub_k[i],ub_r0[i]); } -/* ---------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + only computes tinker_angle() and tinker_bondangle() + does not compute tinker_anglep() and tinker_urey_bradley() +---------------------------------------------------------------------- */ double AngleAmoeba::single(int type, int i1, int i2, int i3) { @@ -866,7 +869,5 @@ double AngleAmoeba::single(int type, int i1, int i2, int i3) double dr2 = r2 - ba_r2[type]; energy += ba_k1[type]*dr1*dtheta + ba_k2[type]*dr2*dtheta; - // NOTE: add UB term - return energy; }