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, 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/cstringutil.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/futil.h"
56 #include "gromacs/utility/smalloc.h"
58 /* All the possible (implemented) table functions */
84 /** Evaluates to true if the table type contains user data. */
85 #define ETAB_USER(e) ((e) == etabUSER || \
86 (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
93 /* This structure holds name and a flag that tells whether
94 this is a Coulomb type funtion */
95 static const t_tab_props tprops[etabNR] = {
98 { "LJ6Shift", FALSE },
99 { "LJ12Shift", FALSE },
105 { "Ewald-Switch", TRUE },
106 { "Ewald-User", TRUE },
107 { "Ewald-User-Switch", TRUE },
108 { "LJ6Ewald", FALSE },
109 { "LJ6Switch", FALSE },
110 { "LJ12Switch", FALSE },
111 { "COULSwitch", TRUE },
112 { "LJ6-Encad shift", FALSE },
113 { "LJ12-Encad shift", FALSE },
114 { "COUL-Encad shift", TRUE },
125 double v_q_ewald_lr(double beta, double r)
129 return beta*2/sqrt(M_PI);
133 return std::erf(beta*r)/r;
137 double v_lj_ewald_lr(double beta, double r)
139 double br, br2, br4, r6, factor;
142 return gmx::power6(beta)/6;
150 factor = (1.0 - std::exp(-br2)*(1 + br2 + 0.5*br4))/r6;
155 void table_spline3_fill_ewald_lr(real *table_f,
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 /* We need some margin to be able to divide table values by r
180 * in the kernel and also to do the integration arithmetics
181 * without going out of range. Furthemore, we divide by dx below.
183 tab_max = GMX_REAL_MAX*0.0001;
185 /* This function produces a table with:
186 * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
187 * maximum force error: V'''/(6*4)*dx^2
188 * The rms force error is the max error times 1/sqrt(5)=0.45.
195 for (i = ntab-1; i >= 0; i--)
199 v_r0 = (*v_lr)(beta, x_r0);
210 /* Linear continuation for the last point in range */
211 vi = v_inrange - dc*(i - i_inrange)*dx;
214 if (table_v != nullptr)
224 /* Get the potential at table point i-1 */
225 v_r1 = (*v_lr)(beta, (i-1)*dx);
227 if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
234 /* Calculate the average second derivative times dx over interval i-1 to i.
235 * Using the function values at the end points and in the middle.
237 a2dx = (v_r0 + v_r1 - 2*(*v_lr)(beta, x_r0-0.5*dx))/(0.25*dx);
238 /* Set the derivative of the spline to match the difference in potential
239 * over the interval plus the average effect of the quadratic term.
240 * This is the essential step for minimizing the error in the force.
242 dc = (v_r0 - v_r1)/dx + 0.5*a2dx;
247 /* Fill the table with the force, minus the derivative of the spline */
252 /* tab[i] will contain the average of the splines over the two intervals */
253 table_f[i] += -0.5*dc;
258 /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
259 * matching the potential at the two end points
260 * and the derivative dc at the end point xr.
264 a2dx = (a1*dx + v_r1 - a0)*2/dx;
266 /* Set dc to the derivative at the next point */
269 if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
279 table_f[(i-1)] = -0.5*dc;
281 /* Currently the last value only contains half the force: double it */
284 if (table_v != nullptr && table_fdv0 != nullptr)
286 /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
287 * init_ewald_f_table().
289 for (i = 0; i < ntab-1; i++)
291 table_fdv0[4*i] = table_f[i];
292 table_fdv0[4*i+1] = table_f[i+1]-table_f[i];
293 table_fdv0[4*i+2] = table_v[i];
294 table_fdv0[4*i+3] = 0.0;
296 table_fdv0[4*(ntab-1)] = table_f[(ntab-1)];
297 table_fdv0[4*(ntab-1)+1] = -table_f[(ntab-1)];
298 table_fdv0[4*(ntab-1)+2] = table_v[(ntab-1)];
299 table_fdv0[4*(ntab-1)+3] = 0.0;
303 /* Returns the spacing for a function using the maximum of
304 * the third derivative, x_scale (unit 1/length)
305 * and function tolerance.
307 static double spline3_table_scale(double third_deriv_max,
312 double sc_deriv, sc_func;
314 /* Force tolerance: single precision accuracy */
315 deriv_tol = GMX_FLOAT_EPS;
316 sc_deriv = sqrt(third_deriv_max/(6*4*deriv_tol*x_scale))*x_scale;
318 /* Don't try to be more accurate on energy than the precision */
319 func_tol = std::max(func_tol, static_cast<double>(GMX_REAL_EPS));
320 sc_func = std::cbrt(third_deriv_max/(6*12*std::sqrt(3.0)*func_tol))*x_scale;
322 return std::max(sc_deriv, sc_func);
325 /* The scale (1/spacing) for third order spline interpolation
326 * of the Ewald mesh contribution which needs to be subtracted
327 * from the non-bonded interactions.
328 * Since there is currently only one spacing for Coulomb and LJ,
329 * the finest spacing is used if both Ewald types are present.
331 * Note that we could also implement completely separate tables
332 * for Coulomb and LJ Ewald, each with their own spacing.
333 * The current setup with the same spacing can provide slightly
334 * faster kernels with both Coulomb and LJ Ewald, especially
335 * when interleaving both tables (currently not implemented).
337 real ewald_spline3_table_scale(const interaction_const_t *ic)
343 if (ic->eeltype == eelEWALD || EEL_PME(ic->eeltype))
345 double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
349 /* Energy tolerance: 0.1 times the cut-off jump */
350 etol = 0.1*std::erfc(ic->ewaldcoeff_q*ic->rcoulomb);
352 sc_q = spline3_table_scale(erf_x_d3, ic->ewaldcoeff_q, etol);
356 fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f 1/nm\n", 1/sc_q);
359 sc = std::max(sc, sc_q);
362 if (EVDW_PME(ic->vdwtype))
364 double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
368 /* Energy tolerance: 0.1 times the cut-off jump */
369 xrc2 = gmx::square(ic->ewaldcoeff_lj*ic->rvdw);
370 etol = 0.1*std::exp(-xrc2)*(1 + xrc2 + xrc2*xrc2/2.0);
372 sc_lj = spline3_table_scale(func_d3, ic->ewaldcoeff_lj, etol);
376 fprintf(debug, "Ewald LJ quadratic spline table spacing: %f 1/nm\n", 1/sc_lj);
379 sc = std::max(sc, sc_lj);
385 static void copy2table(int n, int offset, int stride,
386 const double x[], const double Vtab[], const double Ftab[], real scalefactor,
389 /* Use double prec. for the intermediary variables
390 * and temporary x/vtab/vtab2 data to avoid unnecessary
397 for (i = 0; (i < n); i++)
403 G = 3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h;
404 H = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] + Ftab[i])*h;
408 /* Fill the last entry with a linear potential,
409 * this is mainly for rounding issues with angle and dihedral potentials.
415 nn0 = offset + i*stride;
416 dest[nn0] = scalefactor*Vtab[i];
417 dest[nn0+1] = scalefactor*F;
418 dest[nn0+2] = scalefactor*G;
419 dest[nn0+3] = scalefactor*H;
423 static void init_table(int n, int nx0,
424 double tabscale, t_tabledata *td, gmx_bool bAlloc)
428 td->tabscale = tabscale;
437 static void spline_forces(int nx, double h, const double v[], gmx_bool bS3, gmx_bool bE3,
441 double v3, b_s, b_e, b;
444 /* Formulas can be found in:
445 * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
448 if (nx < 4 && (bS3 || bE3))
450 gmx_fatal(FARGS, "Can not generate splines with third derivative boundary conditions with less than 4 (%d) points", nx);
453 /* To make life easy we initially set the spacing to 1
454 * and correct for this at the end.
458 /* Fit V''' at the start */
459 v3 = v[3] - 3*v[2] + 3*v[1] - v[0];
462 fprintf(debug, "The left third derivative is %g\n", v3/(h*h*h));
464 b_s = 2*(v[1] - v[0]) + v3/6;
469 /* Fit V'' at the start */
472 v2 = -v[3] + 4*v[2] - 5*v[1] + 2*v[0];
473 /* v2 = v[2] - 2*v[1] + v[0]; */
476 fprintf(debug, "The left second derivative is %g\n", v2/(h*h));
478 b_s = 3*(v[1] - v[0]) - v2/2;
484 b_s = 3*(v[2] - v[0]) + f[0]*h;
489 /* Fit V''' at the end */
490 v3 = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4];
493 fprintf(debug, "The right third derivative is %g\n", v3/(h*h*h));
495 b_e = 2*(v[nx-1] - v[nx-2]) + v3/6;
500 /* V'=0 at the end */
501 b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h;
506 beta = (bS3 ? 1 : 4);
508 /* For V'' fitting */
509 /* beta = (bS3 ? 2 : 4); */
512 for (i = start+1; i < end; i++)
516 b = 3*(v[i+1] - v[i-1]);
517 f[i] = (b - f[i-1])/beta;
519 gamma[end-1] = 1/beta;
520 beta = (bE3 ? 1 : 4) - gamma[end-1];
521 f[end-1] = (b_e - f[end-2])/beta;
523 for (i = end-2; i >= start; i--)
525 f[i] -= gamma[i+1]*f[i+1];
529 /* Correct for the minus sign and the spacing */
530 for (i = start; i < end; i++)
536 static void set_forces(FILE *fp, int angle,
537 int nx, double h, double v[], double f[],
545 "Force generation for dihedral tables is not (yet) implemented");
549 while (v[start] == 0)
555 while (v[end-1] == 0)
570 fprintf(fp, "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
571 table+1, start*h, end == nx ? "V'''" : "V'=0", (end-1)*h);
573 spline_forces(end-start, h, v+start, TRUE, end == nx, f+start);
576 static void read_tables(FILE *fp, const char *fn,
577 int ntab, int angle, t_tabledata td[])
581 double **yy = nullptr, start, end, dx0, dx1, ssd, vm, vp, f, numf;
582 int k, i, nx, nx0 = 0, ny, nny, ns;
583 gmx_bool bAllZero, bZeroV, bZeroF;
587 libfn = gmxlibfn(fn);
588 nx = read_xvg(libfn, &yy, &ny);
591 gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d",
599 "The first distance in file %s is %f nm instead of %f nm",
600 libfn, yy[0][0], 0.0);
614 if (yy[0][0] != start || yy[0][nx-1] != end)
616 gmx_fatal(FARGS, "The angles in file %s should go from %f to %f instead of %f to %f\n",
617 libfn, start, end, yy[0][0], yy[0][nx-1]);
621 tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
625 fprintf(fp, "Read user tables from %s with %d data points.\n", libfn, nx);
628 fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
633 for (k = 0; k < ntab; k++)
637 for (i = 0; (i < nx); i++)
641 dx0 = yy[0][i-1] - yy[0][i-2];
642 dx1 = yy[0][i] - yy[0][i-1];
643 /* Check for 1% deviation in spacing */
644 if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1)))
646 gmx_fatal(FARGS, "In table file '%s' the x values are not equally spaced: %f %f %f", fn, yy[0][i-2], yy[0][i-1], yy[0][i]);
649 if (yy[1+k*2][i] != 0)
657 if (yy[1+k*2][i] > 0.01*GMX_REAL_MAX ||
658 yy[1+k*2][i] < -0.01*GMX_REAL_MAX)
660 gmx_fatal(FARGS, "Out of range potential value %g in file '%s'",
664 if (yy[1+k*2+1][i] != 0)
672 if (yy[1+k*2+1][i] > 0.01*GMX_REAL_MAX ||
673 yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX)
675 gmx_fatal(FARGS, "Out of range force value %g in file '%s'",
681 if (!bZeroV && bZeroF)
683 set_forces(fp, angle, nx, 1/tabscale, yy[1+k*2], yy[1+k*2+1], k);
687 /* Check if the second column is close to minus the numerical
688 * derivative of the first column.
692 for (i = 1; (i < nx-1); i++)
697 if (vm != 0 && vp != 0 && f != 0)
699 /* Take the centered difference */
700 numf = -(vp - vm)*0.5*tabscale;
703 ssd += fabs(2*(f - numf)/(f + numf));
711 sprintf(buf, "For the %d non-zero entries for table %d in %s the forces deviate on average %lld%% from minus the numerical derivative of the potential\n", ns, k, libfn, (long long int)(100*ssd+0.5));
714 fprintf(debug, "%s", buf);
720 fprintf(fp, "\nWARNING: %s\n", buf);
722 fprintf(stderr, "\nWARNING: %s\n", buf);
729 fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn);
732 for (k = 0; (k < ntab); k++)
734 init_table(nx, nx0, tabscale, &(td[k]), TRUE);
735 for (i = 0; (i < nx); i++)
737 td[k].x[i] = yy[0][i];
738 td[k].v[i] = yy[2*k+1][i];
739 td[k].f[i] = yy[2*k+2][i];
742 for (i = 0; (i < ny); i++)
750 static void done_tabledata(t_tabledata *td)
762 static void fill_table(t_tabledata *td, int tp, const interaction_const_t *ic,
765 /* Fill the table according to the formulas in the manual.
766 * In principle, we only need the potential and the second
767 * derivative, but then we would have to do lots of calculations
768 * in the inner loop. By precalculating some terms (see manual)
769 * we get better eventual performance, despite a larger table.
771 * Since some of these higher-order terms are very small,
772 * we always use double precision to calculate them here, in order
773 * to avoid unnecessary loss of precision.
777 double r1, rc, r12, r13;
778 double r, r2, r6, rc2, rc6, rc12;
779 double expr, Vtab, Ftab;
780 /* Parameters for David's function */
781 double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
782 /* Parameters for the switching function */
783 double ksw, swi, swi1;
784 /* Temporary parameters */
785 gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
786 double ewc = ic->ewaldcoeff_q;
787 double ewclj = ic->ewaldcoeff_lj;
792 bPotentialSwitch = FALSE;
793 bForceSwitch = FALSE;
794 bPotentialShift = FALSE;
798 bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
799 (tp == etabCOULSwitch) ||
800 (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
801 (tprops[tp].bCoulomb && (ic->coulomb_modifier == eintmodPOTSWITCH)) ||
802 (!tprops[tp].bCoulomb && (ic->vdw_modifier == eintmodPOTSWITCH)));
803 bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
805 (tprops[tp].bCoulomb && (ic->coulomb_modifier == eintmodFORCESWITCH)) ||
806 (!tprops[tp].bCoulomb && (ic->vdw_modifier == eintmodFORCESWITCH)));
807 bPotentialShift = ((tprops[tp].bCoulomb && (ic->coulomb_modifier == eintmodPOTSHIFT)) ||
808 (!tprops[tp].bCoulomb && (ic->vdw_modifier == eintmodPOTSHIFT)));
813 if (tprops[tp].bCoulomb)
815 r1 = ic->rcoulomb_switch;
820 r1 = ic->rvdw_switch;
823 if (bPotentialSwitch)
825 ksw = 1.0/(gmx::power5(rc-r1));
837 else if (tp == etabLJ6Shift)
846 A = p * ((p+1)*r1-(p+4)*rc)/(std::pow(rc, p+2)*gmx::square(rc-r1));
847 B = -p * ((p+1)*r1-(p+3)*rc)/(std::pow(rc, p+2)*gmx::power3(rc-r1));
848 C = 1.0/std::pow(rc, p)-A/3.0*gmx::power3(rc-r1)-B/4.0*gmx::power4(rc-r1);
849 if (tp == etabLJ6Shift)
860 fprintf(debug, "Setting up tables\n"); fflush(debug);
866 rc6 = 1.0/(rc2*rc2*rc2);
867 if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
873 rc12 = std::pow(rc, -reppow);
883 Vcut = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + gmx::power4(ewclj)*rc2*rc2/2);
893 case etabEwaldSwitch:
894 Vcut = std::erfc(ewc*rc)/rc;
897 /* Only calculate minus the reciprocal space contribution */
898 Vcut = -std::erf(ewc*rc)/rc;
902 /* No need for preventing the usage of modifiers with RF */
909 gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
910 tprops[tp].name, __FILE__, __LINE__);
914 for (i = 0; (i < td->nx); i++)
916 td->x[i] = i/td->tabscale;
918 for (i = td->nx0; (i < td->nx); i++)
923 if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
929 r12 = std::pow(r, -reppow);
933 if (bPotentialSwitch)
935 /* swi is function, swi1 1st derivative and swi2 2nd derivative */
936 /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
937 * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
939 * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
953 swi = 1 - 10*gmx::power3(r-r1)*ksw*gmx::square(rc-r1)
954 + 15*gmx::power4(r-r1)*ksw*(rc-r1) - 6*gmx::power5(r-r1)*ksw;
955 swi1 = -30*gmx::square(r-r1)*ksw*gmx::square(rc-r1)
956 + 60*gmx::power3(r-r1)*ksw*(rc-r1) - 30*gmx::power4(r-r1)*ksw;
959 else /* not really needed, but avoids compiler warnings... */
988 Ftab = reppow*Vtab/r;
996 Ftab = reppow*Vtab/r;
1002 Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
1003 Ftab = -(6.0*r6/r-6.0*rc6/rc);
1014 Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
1015 Ftab = -(6.0*r6/r-6.0*rc6/rc);
1027 case etabCOULSwitch:
1036 case etabEwaldSwitch:
1037 Vtab = std::erfc(ewc*r)/r;
1038 Ftab = std::erfc(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
1041 case etabEwaldUserSwitch:
1042 /* Only calculate the negative of the reciprocal space contribution */
1043 Vtab = -std::erf(ewc*r)/r;
1044 Ftab = -std::erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
1047 Vtab = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + gmx::power4(ewclj)*r2*r2/2);
1048 Ftab = 6.0*Vtab/r - r6*exp(-ewclj*ewclj*r2)*gmx::power5(ewclj)*ewclj*r2*r2*r;
1052 Vtab = 1.0/r + ic->k_rf*r2 - ic->c_rf;
1053 Ftab = 1.0/r2 - 2*ic->k_rf*r;
1054 if (tp == etabRF_ZERO && r >= rc)
1068 Vtab = 1.0/r-(rc-r)/(rc*rc)-1.0/rc;
1069 Ftab = 1.0/r2-1.0/(rc*rc);
1078 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
1079 tp, __FILE__, __LINE__);
1083 /* Normal coulomb with cut-off correction for potential */
1087 /* If in Shifting range add something to it */
1090 r12 = (r-r1)*(r-r1);
1092 Vtab += -A_3*r13 - B_4*r12*r12;
1093 Ftab += A*r12 + B*r13;
1098 /* Make sure interactions are zero outside cutoff with modifiers */
1103 if (bPotentialShift)
1111 /* Make sure interactions are zero outside cutoff with modifiers */
1123 if (bPotentialSwitch)
1127 /* Make sure interactions are zero outside cutoff with modifiers */
1133 Ftab = Ftab*swi - Vtab*swi1;
1137 /* Convert to single precision when we store to mem */
1142 /* Continue the table linearly from nx0 to 0.
1143 * These values are only required for energy minimization with overlap or TPI.
1145 for (i = td->nx0-1; i >= 0; i--)
1147 td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
1148 td->f[i] = td->f[i+1];
1152 static void set_table_type(int tabsel[], const interaction_const_t *ic, gmx_bool b14only)
1154 int eltype, vdwtype;
1156 /* Set the different table indices.
1163 switch (ic->eeltype)
1167 case eelPMEUSERSWITCH:
1176 eltype = ic->eeltype;
1182 tabsel[etiCOUL] = etabCOUL;
1185 tabsel[etiCOUL] = etabShift;
1188 if (ic->rcoulomb > ic->rcoulomb_switch)
1190 tabsel[etiCOUL] = etabShift;
1194 tabsel[etiCOUL] = etabCOUL;
1200 tabsel[etiCOUL] = etabEwald;
1203 tabsel[etiCOUL] = etabEwaldSwitch;
1206 tabsel[etiCOUL] = etabEwaldUser;
1208 case eelPMEUSERSWITCH:
1209 tabsel[etiCOUL] = etabEwaldUserSwitch;
1214 tabsel[etiCOUL] = etabRF_ZERO;
1217 tabsel[etiCOUL] = etabCOULSwitch;
1220 tabsel[etiCOUL] = etabUSER;
1223 tabsel[etiCOUL] = etabCOULEncad;
1226 gmx_fatal(FARGS, "Invalid eeltype %d", eltype);
1229 /* Van der Waals time */
1230 if (ic->useBuckingham && !b14only)
1232 tabsel[etiLJ6] = etabLJ6;
1233 tabsel[etiLJ12] = etabEXPMIN;
1237 if (b14only && ic->vdwtype != evdwUSER)
1243 vdwtype = ic->vdwtype;
1249 tabsel[etiLJ6] = etabLJ6Switch;
1250 tabsel[etiLJ12] = etabLJ12Switch;
1253 tabsel[etiLJ6] = etabLJ6Shift;
1254 tabsel[etiLJ12] = etabLJ12Shift;
1257 tabsel[etiLJ6] = etabUSER;
1258 tabsel[etiLJ12] = etabUSER;
1261 tabsel[etiLJ6] = etabLJ6;
1262 tabsel[etiLJ12] = etabLJ12;
1264 case evdwENCADSHIFT:
1265 tabsel[etiLJ6] = etabLJ6Encad;
1266 tabsel[etiLJ12] = etabLJ12Encad;
1269 tabsel[etiLJ6] = etabLJ6Ewald;
1270 tabsel[etiLJ12] = etabLJ12;
1273 gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype,
1274 __FILE__, __LINE__);
1277 if (!b14only && ic->vdw_modifier != eintmodNONE)
1279 if (ic->vdw_modifier != eintmodPOTSHIFT &&
1280 ic->vdwtype != evdwCUT)
1282 gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
1285 /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
1286 * to the original interaction forms when we fill the table, so we only check cutoffs here.
1288 if (ic->vdwtype == evdwCUT)
1290 switch (ic->vdw_modifier)
1293 case eintmodPOTSHIFT:
1294 case eintmodEXACTCUTOFF:
1295 /* No modification */
1297 case eintmodPOTSWITCH:
1298 tabsel[etiLJ6] = etabLJ6Switch;
1299 tabsel[etiLJ12] = etabLJ12Switch;
1301 case eintmodFORCESWITCH:
1302 tabsel[etiLJ6] = etabLJ6Shift;
1303 tabsel[etiLJ12] = etabLJ12Shift;
1306 gmx_incons("Unsupported vdw_modifier");
1313 t_forcetable *make_tables(FILE *out,
1314 const interaction_const_t *ic,
1316 real rtab, int flags)
1319 gmx_bool b14only, useUserTable;
1320 int nx0, tabsel[etiNR];
1323 t_forcetable *table;
1326 b14only = (flags & GMX_MAKETABLES_14ONLY);
1328 if (flags & GMX_MAKETABLES_FORCEUSER)
1330 tabsel[etiCOUL] = etabUSER;
1331 tabsel[etiLJ6] = etabUSER;
1332 tabsel[etiLJ12] = etabUSER;
1336 set_table_type(tabsel, ic, b14only);
1343 table->interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
1344 table->format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1345 table->formatsize = 4;
1346 table->ninteractions = etiNR;
1347 table->stride = table->formatsize*table->ninteractions;
1349 /* Check whether we have to read or generate */
1350 useUserTable = FALSE;
1351 for (unsigned int i = 0; (i < etiNR); i++)
1353 if (ETAB_USER(tabsel[i]))
1355 useUserTable = TRUE;
1360 read_tables(out, fn, etiNR, 0, td);
1361 if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1363 table->n = td[0].nx;
1367 if (td[0].x[td[0].nx-1] < rtab)
1369 gmx_fatal(FARGS, "Tables in file %s not long enough for cut-off:\n"
1370 "\tshould be at least %f nm\n", fn, rtab);
1372 table->n = (int)(rtab*td[0].tabscale + 0.5);
1374 table->scale = td[0].tabscale;
1379 // No tables are read
1381 table->scale = 2000.0;
1383 table->scale = 500.0;
1385 table->n = static_cast<int>(rtab*table->scale);
1389 /* Each table type (e.g. coul,lj6,lj12) requires four
1390 * numbers per table->n+1 data points. For performance reasons we want
1391 * the table data to be aligned to a 32-byte boundary.
1393 snew_aligned(table->data, table->stride*(table->n+1)*sizeof(real), 32);
1395 for (int k = 0; (k < etiNR); k++)
1397 /* Now fill data for tables that have not been read
1398 * or add the Ewald long-range correction for Ewald user tables.
1400 if (tabsel[k] != etabUSER)
1402 real scale = table->scale;
1403 if (ic->useBuckingham &&
1404 (ic->buckinghamBMax != 0) &&
1405 tabsel[k] == etabEXPMIN)
1407 scale /= ic->buckinghamBMax;
1409 init_table(table->n, nx0, scale, &(td[k]), !useUserTable);
1411 fill_table(&(td[k]), tabsel[k], ic, b14only);
1414 fprintf(out, "Generated table with %d data points for %s%s.\n"
1415 "Tabscale = %g points/nm\n",
1416 td[k].nx, b14only ? "1-4 " : "", tprops[tabsel[k]].name,
1421 /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
1422 * by including the derivative constants (6.0 or 12.0) in the parameters, since
1423 * we no longer calculate force in most steps. This means the c6/c12 parameters
1424 * have been scaled up, so we need to scale down the table interactions too.
1425 * It comes here since we need to scale user tables too.
1429 scalefactor = 1.0/6.0;
1431 else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1433 scalefactor = 1.0/12.0;
1440 copy2table(table->n, k*table->formatsize, table->stride, td[k].x, td[k].v, td[k].f, scalefactor, table->data);
1442 done_tabledata(&(td[k]));
1449 bondedtable_t make_bonded_table(FILE *fplog, const char *fn, int angle)
1456 read_tables(fplog, fn, 1, angle, &td);
1459 /* Convert the table from degrees to radians */
1460 for (i = 0; i < td.nx; i++)
1465 td.tabscale *= RAD2DEG;
1468 tab.scale = td.tabscale;
1469 snew(tab.data, tab.n*stride);
1470 copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data);
1471 done_tabledata(&td);
1476 t_forcetable *makeDispersionCorrectionTable(FILE *fp,
1477 const interaction_const_t *ic,
1481 t_forcetable *dispersionCorrectionTable = nullptr;
1483 if (tabfn == nullptr)
1487 fprintf(debug, "No table file name passed, can not read table, can not do non-bonded interactions\n");
1489 return dispersionCorrectionTable;
1492 t_forcetable *fullTable = make_tables(fp, ic, tabfn, rtab, 0);
1493 /* Copy the contents of the table to one that has just dispersion
1494 * and repulsion, to improve cache performance. We want the table
1495 * data to be aligned to 32-byte boundaries. The pointers could be
1496 * freed but currently aren't. */
1497 snew(dispersionCorrectionTable, 1);
1498 dispersionCorrectionTable->interaction = GMX_TABLE_INTERACTION_VDWREP_VDWDISP;
1499 dispersionCorrectionTable->format = fullTable->format;
1500 dispersionCorrectionTable->r = fullTable->r;
1501 dispersionCorrectionTable->n = fullTable->n;
1502 dispersionCorrectionTable->scale = fullTable->scale;
1503 dispersionCorrectionTable->formatsize = fullTable->formatsize;
1504 dispersionCorrectionTable->ninteractions = 2;
1505 dispersionCorrectionTable->stride = dispersionCorrectionTable->formatsize * dispersionCorrectionTable->ninteractions;
1506 snew_aligned(dispersionCorrectionTable->data, dispersionCorrectionTable->stride*(dispersionCorrectionTable->n+1), 32);
1508 for (int i = 0; i <= fullTable->n; i++)
1510 for (int j = 0; j < 8; j++)
1512 dispersionCorrectionTable->data[8*i+j] = fullTable->data[12*i+4+j];
1515 sfree_aligned(fullTable->data);
1518 return dispersionCorrectionTable;