diff --git a/src/USER-MISC/pair_meam_spline.cpp b/src/USER-MISC/pair_meam_spline.cpp index 27deff9a6e..2b47b3dea6 100644 --- a/src/USER-MISC/pair_meam_spline.cpp +++ b/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]; diff --git a/src/USER-MISC/pair_meam_spline.h b/src/USER-MISC/pair_meam_spline.h index 9ef5b57381..2e6e7e00f3 100644 --- a/src/USER-MISC/pair_meam_spline.h +++ b/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*(1/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*(1/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. }; diff --git a/src/USER-MISC/pair_meam_sw_spline.cpp b/src/USER-MISC/pair_meam_sw_spline.cpp index 73d6c81004..8d33a6f04f 100644 --- a/src/USER-MISC/pair_meam_sw_spline.cpp +++ b/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]; diff --git a/src/USER-MISC/pair_meam_sw_spline.h b/src/USER-MISC/pair_meam_sw_spline.h index 9398619b63..b4480050db 100644 --- a/src/USER-MISC/pair_meam_sw_spline.h +++ b/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. };