address some NOTE comments

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;

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 *,

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

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;

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

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

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