Merge pull request #2015 from jvita/fix-spline-meam-binning

Fix spline meam binning

src/USER-MISC/pair_meam_spline.cpp

@@ -717,6 +717,7 @@ void PairMEAMSpline::SplineFunction::prepareSpline(Error* error)
     Y2[i] /= h*6.0;
 #endif
   }
+  inv_h = (1/h);
   xmax_shifted = xmax - xmin;
 }
 
@@ -732,6 +733,7 @@ void PairMEAMSpline::SplineFunction::communicate(MPI_Comm& world, int me)
   MPI_Bcast(&isGridSpline, 1, MPI_INT, 0, world);
   MPI_Bcast(&h, 1, MPI_DOUBLE, 0, world);
   MPI_Bcast(&hsq, 1, MPI_DOUBLE, 0, world);
+  MPI_Bcast(&inv_h, 1, MPI_DOUBLE, 0, world);
   if(me != 0) {
     X = new double[N];
     Xs = new double[N];

src/USER-MISC/pair_meam_spline.h

@@ -142,7 +142,7 @@ protected:
           ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi])*(h*h)/6.0;
 #else
         // For a spline with regular grid, we directly calculate the interval X is in.
-        int klo = (int)(x / h);
+        int klo = (int)(x*inv_h);
         int khi = klo + 1;
         double a = Xs[khi] - x;
         double b = h - a;
@@ -186,7 +186,7 @@ protected:
            (b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
 #else
         // For a spline with regular grid, we directly calculate the interval X is in.
-        int klo = (int)(x / h);
+        int klo = (int)(x*inv_h);
         int khi = klo + 1;
         double a = Xs[khi] - x;
         double b = h - a;
@@ -224,6 +224,7 @@ protected:
     int isGridSpline;// Indicates that all spline knots are on a regular grid.
     double h;        // The distance between knots if this is a grid spline with equidistant knots.
     double hsq;      // The squared distance between knots if this is a grid spline with equidistant knots.
+    double inv_h;    // (1/h), used to avoid numerical errors in binnning for grid spline with equidistant knots.
     double xmax_shifted; // The end of the spline interval after it has been shifted to begin at X=0.
   };
 

src/USER-MISC/pair_meam_sw_spline.cpp

@@ -674,6 +674,7 @@ void PairMEAMSWSpline::SplineFunction::prepareSpline(Error* error)
                 Y2[i] /= h*6.0;
 #endif
         }
+        inv_h = 1/h;
         xmax_shifted = xmax - xmin;
 }
 
@@ -689,6 +690,7 @@ void PairMEAMSWSpline::SplineFunction::communicate(MPI_Comm& world, int me)
         MPI_Bcast(&isGridSpline, 1, MPI_INT, 0, world);
         MPI_Bcast(&h, 1, MPI_DOUBLE, 0, world);
         MPI_Bcast(&hsq, 1, MPI_DOUBLE, 0, world);
+        MPI_Bcast(&inv_h, 1, MPI_DOUBLE, 0, world);
         if(me != 0) {
                 X = new double[N];
                 Xs = new double[N];

src/USER-MISC/pair_meam_sw_spline.h

@@ -130,7 +130,7 @@ protected:
 #else
                                 // For a spline with grid points, we can directly calculate the interval X is in.
                                 //
-                                int klo = (int)(x / h);
+                                int klo = (int)(x*inv_h);
                                 if ( klo > N - 2 ) klo = N - 2;
                                 int khi = klo + 1;
                                 double a = Xs[khi] - x;
@@ -170,7 +170,7 @@ protected:
                                 return a * Y[klo] + b * Y[khi] + ((a*a*a - a) * Y2[klo] + (b*b*b - b) * Y2[khi]) * (h*h) / 6.0;
 #else
                                 // For a spline with grid points, we can directly calculate the interval X is in.
-                                int klo = (int)(x / h);
+                                int klo = (int)(x*inv_h);
                                 if ( klo > N - 2 ) klo = N - 2;
                                 int khi = klo + 1;
                                 double a = Xs[khi] - x;
@@ -207,6 +207,7 @@ protected:
                 int isGridSpline;                // Indicates that all spline knots are on a regular grid.
                 double h;                                // The distance between knots if this is a grid spline with equidistant knots.
                 double hsq;                                // The squared distance between knots if this is a grid spline with equidistant knots.
+                double inv_h;    // (1/h), used to avoid numerical errors in binnning for grid spline with equidistant knots.
                 double xmax_shifted;        // The end of the spline interval after it has been shifted to begin at X=0.
         };