2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "forcetable.h"
45 #include "gromacs/fileio/xvgr.h"
46 #include "gromacs/math/functions.h"
47 #include "gromacs/math/units.h"
48 #include "gromacs/math/utilities.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/mdtypes/fcdata.h"
51 #include "gromacs/mdtypes/md_enums.h"
52 #include "gromacs/mdtypes/nblist.h"
53 #include "gromacs/utility/arrayref.h"
54 #include "gromacs/utility/cstringutil.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/futil.h"
57 #include "gromacs/utility/smalloc.h"
59 /* All the possible (implemented) table functions */
86 /** Evaluates to true if the table type contains user data. */
87 #define ETAB_USER(e) ((e) == etabUSER || (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
95 /* This structure holds name and a flag that tells whether
96 this is a Coulomb type funtion */
97 static const t_tab_props tprops[etabNR] = {
100 { "LJ6Shift", FALSE },
101 { "LJ12Shift", FALSE },
107 { "Ewald-Switch", TRUE },
108 { "Ewald-User", TRUE },
109 { "Ewald-User-Switch", TRUE },
110 { "LJ6Ewald", FALSE },
111 { "LJ6Switch", FALSE },
112 { "LJ12Switch", FALSE },
113 { "COULSwitch", TRUE },
114 { "LJ6-Encad shift", FALSE },
115 { "LJ12-Encad shift", FALSE },
116 { "COUL-Encad shift", TRUE },
128 double v_q_ewald_lr(double beta, double r)
132 return beta * 2 / sqrt(M_PI);
136 return std::erf(beta * r) / r;
140 double v_lj_ewald_lr(double beta, double r)
142 double br, br2, br4, r6, factor;
145 return gmx::power6(beta) / 6;
153 factor = (1.0 - std::exp(-br2) * (1 + br2 + 0.5 * br4)) / r6;
158 EwaldCorrectionTables generateEwaldCorrectionTables(const int numPoints,
159 const double tableScaling,
161 real_space_grid_contribution_computer v_lr)
166 gmx_bool bOutOfRange;
167 double v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
170 /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument
171 * depending on wether we should create electrostatic or Lennard-Jones Ewald tables.
176 gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
179 const double dx = 1 / tableScaling;
181 EwaldCorrectionTables tables;
182 tables.scale = tableScaling;
183 tables.tableF.resize(numPoints);
184 tables.tableV.resize(numPoints);
185 tables.tableFDV0.resize(numPoints * 4);
186 gmx::ArrayRef<real> table_f = tables.tableF;
187 gmx::ArrayRef<real> table_v = tables.tableV;
188 gmx::ArrayRef<real> table_fdv0 = tables.tableFDV0;
190 /* We need some margin to be able to divide table values by r
191 * in the kernel and also to do the integration arithmetics
192 * without going out of range. Furthemore, we divide by dx below.
194 tab_max = GMX_REAL_MAX * 0.0001;
196 /* This function produces a table with:
197 * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
198 * maximum force error: V'''/(6*4)*dx^2
199 * The rms force error is the max error times 1/sqrt(5)=0.45.
203 i_inrange = numPoints;
206 for (i = numPoints - 1; i >= 0; i--)
210 v_r0 = (*v_lr)(beta, x_r0);
221 /* Linear continuation for the last point in range */
222 vi = v_inrange - dc * (i - i_inrange) * dx;
232 /* Get the potential at table point i-1 */
233 v_r1 = (*v_lr)(beta, (i - 1) * dx);
235 if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
242 /* Calculate the average second derivative times dx over interval i-1 to i.
243 * Using the function values at the end points and in the middle.
245 a2dx = (v_r0 + v_r1 - 2 * (*v_lr)(beta, x_r0 - 0.5 * dx)) / (0.25 * dx);
246 /* Set the derivative of the spline to match the difference in potential
247 * over the interval plus the average effect of the quadratic term.
248 * This is the essential step for minimizing the error in the force.
250 dc = (v_r0 - v_r1) / dx + 0.5 * a2dx;
253 if (i == numPoints - 1)
255 /* Fill the table with the force, minus the derivative of the spline */
260 /* tab[i] will contain the average of the splines over the two intervals */
261 table_f[i] += -0.5 * dc;
266 /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
267 * matching the potential at the two end points
268 * and the derivative dc at the end point xr.
272 a2dx = (a1 * dx + v_r1 - a0) * 2 / dx;
274 /* Set dc to the derivative at the next point */
277 if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
287 table_f[(i - 1)] = -0.5 * dc;
289 /* Currently the last value only contains half the force: double it */
292 if (!table_fdv0.empty())
294 /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
295 * init_ewald_f_table().
297 for (i = 0; i < numPoints - 1; i++)
299 table_fdv0[4 * i] = table_f[i];
300 table_fdv0[4 * i + 1] = table_f[i + 1] - table_f[i];
301 table_fdv0[4 * i + 2] = table_v[i];
302 table_fdv0[4 * i + 3] = 0.0;
304 const int lastPoint = numPoints - 1;
305 table_fdv0[4 * lastPoint] = table_f[lastPoint];
306 table_fdv0[4 * lastPoint + 1] = -table_f[lastPoint];
307 table_fdv0[4 * lastPoint + 2] = table_v[lastPoint];
308 table_fdv0[4 * lastPoint + 3] = 0.0;
314 /* Returns the spacing for a function using the maximum of
315 * the third derivative, x_scale (unit 1/length)
316 * and function tolerance.
318 static double spline3_table_scale(double third_deriv_max, double x_scale, double func_tol)
321 double sc_deriv, sc_func;
323 /* Force tolerance: single precision accuracy */
324 deriv_tol = GMX_FLOAT_EPS;
325 sc_deriv = sqrt(third_deriv_max / (6 * 4 * deriv_tol * x_scale)) * x_scale;
327 /* Don't try to be more accurate on energy than the precision */
328 func_tol = std::max(func_tol, static_cast<double>(GMX_REAL_EPS));
329 sc_func = std::cbrt(third_deriv_max / (6 * 12 * std::sqrt(3.0) * func_tol)) * x_scale;
331 return std::max(sc_deriv, sc_func);
334 /* The scale (1/spacing) for third order spline interpolation
335 * of the Ewald mesh contribution which needs to be subtracted
336 * from the non-bonded interactions.
337 * Since there is currently only one spacing for Coulomb and LJ,
338 * the finest spacing is used if both Ewald types are present.
340 * Note that we could also implement completely separate tables
341 * for Coulomb and LJ Ewald, each with their own spacing.
342 * The current setup with the same spacing can provide slightly
343 * faster kernels with both Coulomb and LJ Ewald, especially
344 * when interleaving both tables (currently not implemented).
346 real ewald_spline3_table_scale(const interaction_const_t& ic,
347 const bool generateCoulombTables,
348 const bool generateVdwTables)
350 GMX_RELEASE_ASSERT(!generateCoulombTables || EEL_PME_EWALD(ic.eeltype),
351 "Can only use tables with Ewald");
352 GMX_RELEASE_ASSERT(!generateVdwTables || EVDW_PME(ic.vdwtype),
353 "Can only use tables with Ewald");
357 if (generateCoulombTables)
359 GMX_RELEASE_ASSERT(ic.ewaldcoeff_q > 0, "The Ewald coefficient shoule be positive");
361 double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
365 /* Energy tolerance: 0.1 times the cut-off jump */
366 etol = 0.1 * std::erfc(ic.ewaldcoeff_q * ic.rcoulomb);
368 sc_q = spline3_table_scale(erf_x_d3, ic.ewaldcoeff_q, etol);
372 fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f nm\n", 1 / sc_q);
375 sc = std::max(sc, sc_q);
378 if (generateVdwTables)
380 GMX_RELEASE_ASSERT(ic.ewaldcoeff_lj > 0, "The Ewald coefficient shoule be positive");
382 double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
386 /* Energy tolerance: 0.1 times the cut-off jump */
387 xrc2 = gmx::square(ic.ewaldcoeff_lj * ic.rvdw);
388 etol = 0.1 * std::exp(-xrc2) * (1 + xrc2 + xrc2 * xrc2 / 2.0);
390 sc_lj = spline3_table_scale(func_d3, ic.ewaldcoeff_lj, etol);
394 fprintf(debug, "Ewald LJ quadratic spline table spacing: %f nm\n", 1 / sc_lj);
397 sc = std::max(sc, sc_lj);
403 static void copy2table(int n,
412 /* Use double prec. for the intermediary variables
413 * and temporary x/vtab/vtab2 data to avoid unnecessary
420 for (i = 0; (i < n); i++)
426 G = 3 * (Vtab[i + 1] - Vtab[i]) + (Ftab[i + 1] + 2 * Ftab[i]) * h;
427 H = -2 * (Vtab[i + 1] - Vtab[i]) - (Ftab[i + 1] + Ftab[i]) * h;
431 /* Fill the last entry with a linear potential,
432 * this is mainly for rounding issues with angle and dihedral potentials.
438 nn0 = offset + i * stride;
439 dest[nn0] = scalefactor * Vtab[i];
440 dest[nn0 + 1] = scalefactor * F;
441 dest[nn0 + 2] = scalefactor * G;
442 dest[nn0 + 3] = scalefactor * H;
446 static void init_table(int n, int nx0, double tabscale, t_tabledata* td, gmx_bool bAlloc)
450 td->tabscale = tabscale;
459 static void spline_forces(int nx, double h, const double v[], gmx_bool bS3, gmx_bool bE3, double f[])
462 double v3, b_s, b_e, b;
465 /* Formulas can be found in:
466 * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
469 if (nx < 4 && (bS3 || bE3))
472 "Can not generate splines with third derivative boundary conditions with less "
473 "than 4 (%d) points",
477 /* To make life easy we initially set the spacing to 1
478 * and correct for this at the end.
482 /* Fit V''' at the start */
483 v3 = v[3] - 3 * v[2] + 3 * v[1] - v[0];
486 fprintf(debug, "The left third derivative is %g\n", v3 / (h * h * h));
488 b_s = 2 * (v[1] - v[0]) + v3 / 6;
493 /* Fit V'' at the start */
496 v2 = -v[3] + 4 * v[2] - 5 * v[1] + 2 * v[0];
497 /* v2 = v[2] - 2*v[1] + v[0]; */
500 fprintf(debug, "The left second derivative is %g\n", v2 / (h * h));
502 b_s = 3 * (v[1] - v[0]) - v2 / 2;
508 b_s = 3 * (v[2] - v[0]) + f[0] * h;
513 /* Fit V''' at the end */
514 v3 = v[nx - 1] - 3 * v[nx - 2] + 3 * v[nx - 3] - v[nx - 4];
517 fprintf(debug, "The right third derivative is %g\n", v3 / (h * h * h));
519 b_e = 2 * (v[nx - 1] - v[nx - 2]) + v3 / 6;
524 /* V'=0 at the end */
525 b_e = 3 * (v[nx - 1] - v[nx - 3]) + f[nx - 1] * h;
530 beta = (bS3 ? 1 : 4);
532 /* For V'' fitting */
533 /* beta = (bS3 ? 2 : 4); */
535 f[start] = b_s / beta;
536 for (i = start + 1; i < end; i++)
540 b = 3 * (v[i + 1] - v[i - 1]);
541 f[i] = (b - f[i - 1]) / beta;
543 gamma[end - 1] = 1 / beta;
544 beta = (bE3 ? 1 : 4) - gamma[end - 1];
545 f[end - 1] = (b_e - f[end - 2]) / beta;
547 for (i = end - 2; i >= start; i--)
549 f[i] -= gamma[i + 1] * f[i + 1];
553 /* Correct for the minus sign and the spacing */
554 for (i = start; i < end; i++)
560 static void set_forces(FILE* fp, int angle, int nx, double h, double v[], double f[], int table)
566 gmx_fatal(FARGS, "Force generation for dihedral tables is not (yet) implemented");
570 while (v[start] == 0)
576 while (v[end - 1] == 0)
591 fprintf(fp, "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
592 table + 1, start * h, end == nx ? "V'''" : "V'=0", (end - 1) * h);
594 spline_forces(end - start, h, v + start, TRUE, end == nx, f + start);
597 static void read_tables(FILE* fp, const char* filename, int ntab, int angle, t_tabledata td[])
600 double **yy = nullptr, start, end, dx0, dx1, ssd, vm, vp, f, numf;
601 int k, i, nx, nx0 = 0, ny, nny, ns;
602 gmx_bool bAllZero, bZeroV, bZeroF;
606 std::string libfn = gmx::findLibraryFile(filename);
607 nx = read_xvg(libfn.c_str(), &yy, &ny);
610 gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d", libfn.c_str(),
617 gmx_fatal(FARGS, "The first distance in file %s is %f nm instead of %f nm",
618 libfn.c_str(), yy[0][0], 0.0);
632 if (yy[0][0] != start || yy[0][nx - 1] != end)
634 gmx_fatal(FARGS, "The angles in file %s should go from %f to %f instead of %f to %f\n",
635 libfn.c_str(), start, end, yy[0][0], yy[0][nx - 1]);
639 tabscale = (nx - 1) / (yy[0][nx - 1] - yy[0][0]);
643 fprintf(fp, "Read user tables from %s with %d data points.\n", libfn.c_str(), nx);
646 fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
651 for (k = 0; k < ntab; k++)
655 for (i = 0; (i < nx); i++)
659 dx0 = yy[0][i - 1] - yy[0][i - 2];
660 dx1 = yy[0][i] - yy[0][i - 1];
661 /* Check for 1% deviation in spacing */
662 if (fabs(dx1 - dx0) >= 0.005 * (fabs(dx0) + fabs(dx1)))
665 "In table file '%s' the x values are not equally spaced: %f %f %f",
666 filename, yy[0][i - 2], yy[0][i - 1], yy[0][i]);
669 if (yy[1 + k * 2][i] != 0)
677 if (yy[1 + k * 2][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2][i] < -0.01 * GMX_REAL_MAX)
679 gmx_fatal(FARGS, "Out of range potential value %g in file '%s'",
680 yy[1 + k * 2][i], filename);
683 if (yy[1 + k * 2 + 1][i] != 0)
691 if (yy[1 + k * 2 + 1][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2 + 1][i] < -0.01 * GMX_REAL_MAX)
693 gmx_fatal(FARGS, "Out of range force value %g in file '%s'",
694 yy[1 + k * 2 + 1][i], filename);
699 if (!bZeroV && bZeroF)
701 set_forces(fp, angle, nx, 1 / tabscale, yy[1 + k * 2], yy[1 + k * 2 + 1], k);
705 /* Check if the second column is close to minus the numerical
706 * derivative of the first column.
710 for (i = 1; (i < nx - 1); i++)
712 vm = yy[1 + 2 * k][i - 1];
713 vp = yy[1 + 2 * k][i + 1];
714 f = yy[1 + 2 * k + 1][i];
715 if (vm != 0 && vp != 0 && f != 0)
717 /* Take the centered difference */
718 numf = -(vp - vm) * 0.5 * tabscale;
721 ssd += fabs(2 * (f - numf) / (f + numf));
730 "For the %d non-zero entries for table %d in %s the forces deviate on "
732 "%% from minus the numerical derivative of the potential\n",
733 ns, k, libfn.c_str(), gmx::roundToInt64(100 * ssd));
736 fprintf(debug, "%s", buf);
742 fprintf(fp, "\nWARNING: %s\n", buf);
744 fprintf(stderr, "\nWARNING: %s\n", buf);
751 fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn.c_str());
754 for (k = 0; (k < ntab); k++)
756 init_table(nx, nx0, tabscale, &(td[k]), TRUE);
757 for (i = 0; (i < nx); i++)
759 td[k].x[i] = yy[0][i];
760 td[k].v[i] = yy[2 * k + 1][i];
761 td[k].f[i] = yy[2 * k + 2][i];
764 for (i = 0; (i < ny); i++)
771 static void done_tabledata(t_tabledata* td)
783 static void fill_table(t_tabledata* td, int tp, const interaction_const_t* ic, gmx_bool b14only)
785 /* Fill the table according to the formulas in the manual.
786 * In principle, we only need the potential and the second
787 * derivative, but then we would have to do lots of calculations
788 * in the inner loop. By precalculating some terms (see manual)
789 * we get better eventual performance, despite a larger table.
791 * Since some of these higher-order terms are very small,
792 * we always use double precision to calculate them here, in order
793 * to avoid unnecessary loss of precision.
797 double r1, rc, r12, r13;
798 double r, r2, r6, rc2, rc6, rc12;
799 double expr, Vtab, Ftab;
800 /* Parameters for David's function */
801 double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
802 /* Parameters for the switching function */
803 double ksw, swi, swi1;
804 /* Temporary parameters */
805 gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
806 double ewc = ic->ewaldcoeff_q;
807 double ewclj = ic->ewaldcoeff_lj;
812 bPotentialSwitch = FALSE;
813 bForceSwitch = FALSE;
814 bPotentialShift = FALSE;
818 bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) || (tp == etabCOULSwitch)
819 || (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch)
820 || (tprops[tp].bCoulomb && (ic->coulomb_modifier == eintmodPOTSWITCH))
821 || (!tprops[tp].bCoulomb && (ic->vdw_modifier == eintmodPOTSWITCH)));
822 bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) || (tp == etabShift)
823 || (tprops[tp].bCoulomb && (ic->coulomb_modifier == eintmodFORCESWITCH))
824 || (!tprops[tp].bCoulomb && (ic->vdw_modifier == eintmodFORCESWITCH)));
825 bPotentialShift = ((tprops[tp].bCoulomb && (ic->coulomb_modifier == eintmodPOTSHIFT))
826 || (!tprops[tp].bCoulomb && (ic->vdw_modifier == eintmodPOTSHIFT)));
831 if (tprops[tp].bCoulomb)
833 r1 = ic->rcoulomb_switch;
838 r1 = ic->rvdw_switch;
841 if (bPotentialSwitch)
843 ksw = 1.0 / (gmx::power5(rc - r1));
855 else if (tp == etabLJ6Shift)
864 A = p * ((p + 1) * r1 - (p + 4) * rc) / (std::pow(rc, p + 2) * gmx::square(rc - r1));
865 B = -p * ((p + 1) * r1 - (p + 3) * rc) / (std::pow(rc, p + 2) * gmx::power3(rc - r1));
866 C = 1.0 / std::pow(rc, p) - A / 3.0 * gmx::power3(rc - r1) - B / 4.0 * gmx::power4(rc - r1);
867 if (tp == etabLJ6Shift)
878 fprintf(debug, "Setting up tables\n");
885 rc6 = 1.0 / (rc2 * rc2 * rc2);
886 if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
892 rc12 = std::pow(rc, -reppow);
902 Vcut = -rc6 * exp(-ewclj * ewclj * rc2)
903 * (1 + ewclj * ewclj * rc2 + gmx::power4(ewclj) * rc2 * rc2 / 2);
909 case etabCOUL: Vcut = 1.0 / rc; break;
911 case etabEwaldSwitch: Vcut = std::erfc(ewc * rc) / rc; break;
913 /* Only calculate minus the reciprocal space contribution */
914 Vcut = -std::erf(ewc * rc) / rc;
918 /* No need for preventing the usage of modifiers with RF */
921 case etabEXPMIN: Vcut = exp(-rc); break;
924 "Cannot apply new potential-shift modifier to interaction type '%s' yet. "
926 tprops[tp].name, __FILE__, __LINE__);
930 for (i = 0; (i < td->nx); i++)
932 td->x[i] = i / td->tabscale;
934 for (i = td->nx0; (i < td->nx); i++)
938 r6 = 1.0 / (r2 * r2 * r2);
939 if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
945 r12 = std::pow(r, -reppow);
949 if (bPotentialSwitch)
951 /* swi is function, swi1 1st derivative and swi2 2nd derivative */
952 /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
953 * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
955 * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
969 swi = 1 - 10 * gmx::power3(r - r1) * ksw * gmx::square(rc - r1)
970 + 15 * gmx::power4(r - r1) * ksw * (rc - r1) - 6 * gmx::power5(r - r1) * ksw;
971 swi1 = -30 * gmx::square(r - r1) * ksw * gmx::square(rc - r1)
972 + 60 * gmx::power3(r - r1) * ksw * (rc - r1) - 30 * gmx::power4(r - r1) * ksw;
975 else /* not really needed, but avoids compiler warnings... */
982 rc6 = 1.0 / (rc6 * rc6);
989 Ftab = 6.0 * Vtab / r;
997 Ftab = 6.0 * Vtab / r;
1004 Ftab = reppow * Vtab / r;
1006 case etabLJ12Switch:
1012 Ftab = reppow * Vtab / r;
1018 Vtab = -(r6 - 6.0 * (rc - r) * rc6 / rc - rc6);
1019 Ftab = -(6.0 * r6 / r - 6.0 * rc6 / rc);
1030 Vtab = -(r6 - 6.0 * (rc - r) * rc6 / rc - rc6);
1031 Ftab = -(6.0 * r6 / r - 6.0 * rc6 / rc);
1043 case etabCOULSwitch:
1052 case etabEwaldSwitch:
1053 Vtab = std::erfc(ewc * r) / r;
1054 Ftab = std::erfc(ewc * r) / r2 + exp(-(ewc * ewc * r2)) * ewc * M_2_SQRTPI / r;
1057 case etabEwaldUserSwitch:
1058 /* Only calculate the negative of the reciprocal space contribution */
1059 Vtab = -std::erf(ewc * r) / r;
1060 Ftab = -std::erf(ewc * r) / r2 + exp(-(ewc * ewc * r2)) * ewc * M_2_SQRTPI / r;
1063 Vtab = -r6 * exp(-ewclj * ewclj * r2)
1064 * (1 + ewclj * ewclj * r2 + gmx::power4(ewclj) * r2 * r2 / 2);
1065 Ftab = 6.0 * Vtab / r
1066 - r6 * exp(-ewclj * ewclj * r2) * gmx::power5(ewclj) * ewclj * r2 * r2 * r;
1070 Vtab = 1.0 / r + ic->k_rf * r2 - ic->c_rf;
1071 Ftab = 1.0 / r2 - 2 * ic->k_rf * r;
1072 if (tp == etabRF_ZERO && r >= rc)
1086 Vtab = 1.0 / r - (rc - r) / (rc * rc) - 1.0 / rc;
1087 Ftab = 1.0 / r2 - 1.0 / (rc * rc);
1096 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)", tp, __FILE__, __LINE__);
1100 /* Normal coulomb with cut-off correction for potential */
1104 /* If in Shifting range add something to it */
1107 r12 = (r - r1) * (r - r1);
1108 r13 = (r - r1) * r12;
1109 Vtab += -A_3 * r13 - B_4 * r12 * r12;
1110 Ftab += A * r12 + B * r13;
1115 /* Make sure interactions are zero outside cutoff with modifiers */
1120 if (bPotentialShift)
1128 /* Make sure interactions are zero outside cutoff with modifiers */
1140 if (bPotentialSwitch)
1144 /* Make sure interactions are zero outside cutoff with modifiers */
1150 Ftab = Ftab * swi - Vtab * swi1;
1154 /* Convert to single precision when we store to mem */
1159 /* Continue the table linearly from nx0 to 0.
1160 * These values are only required for energy minimization with overlap or TPI.
1162 for (i = td->nx0 - 1; i >= 0; i--)
1164 td->v[i] = td->v[i + 1] + td->f[i + 1] * (td->x[i + 1] - td->x[i]);
1165 td->f[i] = td->f[i + 1];
1169 static void set_table_type(int tabsel[], const interaction_const_t* ic, gmx_bool b14only)
1171 int eltype, vdwtype;
1173 /* Set the different table indices.
1180 switch (ic->eeltype)
1184 case eelPMEUSERSWITCH: eltype = eelUSER; break;
1185 default: eltype = eelCUT;
1190 eltype = ic->eeltype;
1195 case eelCUT: tabsel[etiCOUL] = etabCOUL; break;
1196 case eelPOISSON: tabsel[etiCOUL] = etabShift; break;
1198 if (ic->rcoulomb > ic->rcoulomb_switch)
1200 tabsel[etiCOUL] = etabShift;
1204 tabsel[etiCOUL] = etabCOUL;
1209 case eelP3M_AD: tabsel[etiCOUL] = etabEwald; break;
1210 case eelPMESWITCH: tabsel[etiCOUL] = etabEwaldSwitch; break;
1211 case eelPMEUSER: tabsel[etiCOUL] = etabEwaldUser; break;
1212 case eelPMEUSERSWITCH: tabsel[etiCOUL] = etabEwaldUserSwitch; break;
1214 case eelRF_ZERO: tabsel[etiCOUL] = etabRF_ZERO; break;
1215 case eelSWITCH: tabsel[etiCOUL] = etabCOULSwitch; break;
1216 case eelUSER: tabsel[etiCOUL] = etabUSER; break;
1217 case eelENCADSHIFT: tabsel[etiCOUL] = etabCOULEncad; break;
1218 default: gmx_fatal(FARGS, "Invalid eeltype %d", eltype);
1221 /* Van der Waals time */
1222 if (ic->useBuckingham && !b14only)
1224 tabsel[etiLJ6] = etabLJ6;
1225 tabsel[etiLJ12] = etabEXPMIN;
1229 if (b14only && ic->vdwtype != evdwUSER)
1235 vdwtype = ic->vdwtype;
1241 tabsel[etiLJ6] = etabLJ6Switch;
1242 tabsel[etiLJ12] = etabLJ12Switch;
1245 tabsel[etiLJ6] = etabLJ6Shift;
1246 tabsel[etiLJ12] = etabLJ12Shift;
1249 tabsel[etiLJ6] = etabUSER;
1250 tabsel[etiLJ12] = etabUSER;
1253 tabsel[etiLJ6] = etabLJ6;
1254 tabsel[etiLJ12] = etabLJ12;
1256 case evdwENCADSHIFT:
1257 tabsel[etiLJ6] = etabLJ6Encad;
1258 tabsel[etiLJ12] = etabLJ12Encad;
1261 tabsel[etiLJ6] = etabLJ6Ewald;
1262 tabsel[etiLJ12] = etabLJ12;
1265 gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype, __FILE__, __LINE__);
1268 if (!b14only && ic->vdw_modifier != eintmodNONE)
1270 if (ic->vdw_modifier != eintmodPOTSHIFT && ic->vdwtype != evdwCUT)
1273 "Potential modifiers other than potential-shift are only implemented for "
1277 /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
1278 * to the original interaction forms when we fill the table, so we only check cutoffs here.
1280 if (ic->vdwtype == evdwCUT)
1282 switch (ic->vdw_modifier)
1285 case eintmodPOTSHIFT:
1286 case eintmodEXACTCUTOFF:
1287 /* No modification */
1289 case eintmodPOTSWITCH:
1290 tabsel[etiLJ6] = etabLJ6Switch;
1291 tabsel[etiLJ12] = etabLJ12Switch;
1293 case eintmodFORCESWITCH:
1294 tabsel[etiLJ6] = etabLJ6Shift;
1295 tabsel[etiLJ12] = etabLJ12Shift;
1297 default: gmx_incons("Unsupported vdw_modifier");
1304 t_forcetable* make_tables(FILE* out, const interaction_const_t* ic, const char* fn, real rtab, int flags)
1307 gmx_bool b14only, useUserTable;
1308 int nx0, tabsel[etiNR];
1311 t_forcetable* table = new t_forcetable(GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
1312 GMX_TABLE_FORMAT_CUBICSPLINE_YFGH);
1314 b14only = ((flags & GMX_MAKETABLES_14ONLY) != 0);
1316 if (flags & GMX_MAKETABLES_FORCEUSER)
1318 tabsel[etiCOUL] = etabUSER;
1319 tabsel[etiLJ6] = etabUSER;
1320 tabsel[etiLJ12] = etabUSER;
1324 set_table_type(tabsel, ic, b14only);
1331 table->formatsize = 4;
1332 table->ninteractions = etiNR;
1333 table->stride = table->formatsize * table->ninteractions;
1335 /* Check whether we have to read or generate */
1336 useUserTable = FALSE;
1337 for (unsigned int i = 0; (i < etiNR); i++)
1339 if (ETAB_USER(tabsel[i]))
1341 useUserTable = TRUE;
1346 read_tables(out, fn, etiNR, 0, td);
1347 if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1349 table->n = td[0].nx;
1353 if (td[0].x[td[0].nx - 1] < rtab)
1356 "Tables in file %s not long enough for cut-off:\n"
1357 "\tshould be at least %f nm\n",
1360 table->n = gmx::roundToInt(rtab * td[0].tabscale);
1362 table->scale = td[0].tabscale;
1367 // No tables are read
1369 table->scale = 2000.0;
1371 table->scale = 500.0;
1373 table->n = static_cast<int>(rtab * table->scale);
1377 /* Each table type (e.g. coul,lj6,lj12) requires four
1378 * numbers per table->n+1 data points. For performance reasons we want
1379 * the table data to be aligned to a 32-byte boundary.
1381 snew_aligned(table->data, table->stride * (table->n + 1) * sizeof(real), 32);
1383 for (int k = 0; (k < etiNR); k++)
1385 /* Now fill data for tables that have not been read
1386 * or add the Ewald long-range correction for Ewald user tables.
1388 if (tabsel[k] != etabUSER)
1390 real scale = table->scale;
1391 if (ic->useBuckingham && (ic->buckinghamBMax != 0) && tabsel[k] == etabEXPMIN)
1393 scale /= ic->buckinghamBMax;
1395 init_table(table->n, nx0, scale, &(td[k]), !useUserTable);
1397 fill_table(&(td[k]), tabsel[k], ic, b14only);
1401 "Generated table with %d data points for %s%s.\n"
1402 "Tabscale = %g points/nm\n",
1403 td[k].nx, b14only ? "1-4 " : "", tprops[tabsel[k]].name, td[k].tabscale);
1407 /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
1408 * by including the derivative constants (6.0 or 12.0) in the parameters, since
1409 * we no longer calculate force in most steps. This means the c6/c12 parameters
1410 * have been scaled up, so we need to scale down the table interactions too.
1411 * It comes here since we need to scale user tables too.
1415 scalefactor = 1.0 / 6.0;
1417 else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1419 scalefactor = 1.0 / 12.0;
1426 copy2table(table->n, k * table->formatsize, table->stride, td[k].x, td[k].v, td[k].f,
1427 scalefactor, table->data);
1429 done_tabledata(&(td[k]));
1436 bondedtable_t make_bonded_table(FILE* fplog, const char* fn, int angle)
1443 read_tables(fplog, fn, 1, angle, &td);
1446 /* Convert the table from degrees to radians */
1447 for (i = 0; i < td.nx; i++)
1452 td.tabscale *= RAD2DEG;
1455 tab.scale = td.tabscale;
1456 snew(tab.data, tab.n * stride);
1457 copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data);
1458 done_tabledata(&td);
1463 std::unique_ptr<t_forcetable>
1464 makeDispersionCorrectionTable(FILE* fp, const interaction_const_t* ic, real rtab, const char* tabfn)
1466 GMX_RELEASE_ASSERT(ic->vdwtype != evdwUSER || tabfn,
1467 "With VdW user tables we need a table file name");
1469 if (tabfn == nullptr)
1471 return std::unique_ptr<t_forcetable>(nullptr);
1474 t_forcetable* fullTable = make_tables(fp, ic, tabfn, rtab, 0);
1475 /* Copy the contents of the table to one that has just dispersion
1476 * and repulsion, to improve cache performance. We want the table
1477 * data to be aligned to 32-byte boundaries.
1479 std::unique_ptr<t_forcetable> dispersionCorrectionTable =
1480 std::make_unique<t_forcetable>(GMX_TABLE_INTERACTION_VDWREP_VDWDISP, fullTable->format);
1481 dispersionCorrectionTable->r = fullTable->r;
1482 dispersionCorrectionTable->n = fullTable->n;
1483 dispersionCorrectionTable->scale = fullTable->scale;
1484 dispersionCorrectionTable->formatsize = fullTable->formatsize;
1485 dispersionCorrectionTable->ninteractions = 2;
1486 dispersionCorrectionTable->stride =
1487 dispersionCorrectionTable->formatsize * dispersionCorrectionTable->ninteractions;
1488 snew_aligned(dispersionCorrectionTable->data,
1489 dispersionCorrectionTable->stride * (dispersionCorrectionTable->n + 1), 32);
1491 for (int i = 0; i <= fullTable->n; i++)
1493 for (int j = 0; j < 8; j++)
1495 dispersionCorrectionTable->data[8 * i + j] = fullTable->data[12 * i + 4 + j];
1500 return dispersionCorrectionTable;
1503 t_forcetable::t_forcetable(enum gmx_table_interaction interaction, enum gmx_table_format format) :
1504 interaction(interaction),
1516 t_forcetable::~t_forcetable()
1518 sfree_aligned(data);