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, 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.
42 #include "gromacs/math/utilities.h"
45 #include "gromacs/utility/smalloc.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/futil.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/math/vec.h"
51 #include "gromacs/math/units.h"
53 #include "gromacs/fileio/gmxfio.h"
57 /* All the possible (implemented) table functions */
83 /** Evaluates to true if the table type contains user data. */
84 #define ETAB_USER(e) ((e) == etabUSER || \
85 (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
92 /* This structure holds name and a flag that tells whether
93 this is a Coulomb type funtion */
94 static const t_tab_props tprops[etabNR] = {
97 { "LJ6Shift", FALSE },
98 { "LJ12Shift", FALSE },
104 { "Ewald-Switch", TRUE },
105 { "Ewald-User", TRUE },
106 { "Ewald-User-Switch", TRUE },
107 { "LJ6Ewald", FALSE },
108 { "LJ6Switch", FALSE },
109 { "LJ12Switch", FALSE },
110 { "COULSwitch", TRUE },
111 { "LJ6-Encad shift", FALSE },
112 { "LJ12-Encad shift", FALSE },
113 { "COUL-Encad shift", TRUE },
118 /* Index in the table that says which function to use */
120 etiCOUL, etiLJ6, etiLJ12, etiNR
129 #define pow2(x) ((x)*(x))
130 #define pow3(x) ((x)*(x)*(x))
131 #define pow4(x) ((x)*(x)*(x)*(x))
132 #define pow5(x) ((x)*(x)*(x)*(x)*(x))
134 double v_q_ewald_lr(double beta, double r)
138 return beta*2/sqrt(M_PI);
142 return gmx_erfd(beta*r)/r;
146 double v_lj_ewald_lr(double beta, double r)
148 double br, br2, br4, r6, factor;
151 return pow(beta, 6)/6;
159 factor = (1.0 - exp(-br2)*(1 + br2 + 0.5*br4))/r6;
164 void table_spline3_fill_ewald_lr(real *table_f,
170 real_space_grid_contribution_computer v_lr)
175 gmx_bool bOutOfRange;
176 double v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
179 /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument
180 * depending on wether we should create electrostatic or Lennard-Jones Ewald tables.
185 gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
188 /* We need some margin to be able to divide table values by r
189 * in the kernel and also to do the integration arithmetics
190 * without going out of range. Furthemore, we divide by dx below.
192 tab_max = GMX_REAL_MAX*0.0001;
194 /* This function produces a table with:
195 * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
196 * maximum force error: V'''/(6*4)*dx^2
197 * The rms force error is the max error times 1/sqrt(5)=0.45.
204 for (i = ntab-1; i >= 0; i--)
208 v_r0 = (*v_lr)(beta, x_r0);
219 /* Linear continuation for the last point in range */
220 vi = v_inrange - dc*(i - i_inrange)*dx;
233 /* Get the potential at table point i-1 */
234 v_r1 = (*v_lr)(beta, (i-1)*dx);
236 if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
243 /* Calculate the average second derivative times dx over interval i-1 to i.
244 * Using the function values at the end points and in the middle.
246 a2dx = (v_r0 + v_r1 - 2*(*v_lr)(beta, x_r0-0.5*dx))/(0.25*dx);
247 /* Set the derivative of the spline to match the difference in potential
248 * over the interval plus the average effect of the quadratic term.
249 * This is the essential step for minimizing the error in the force.
251 dc = (v_r0 - v_r1)/dx + 0.5*a2dx;
256 /* Fill the table with the force, minus the derivative of the spline */
261 /* tab[i] will contain the average of the splines over the two intervals */
262 table_f[i] += -0.5*dc;
267 /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
268 * matching the potential at the two end points
269 * and the derivative dc at the end point xr.
273 a2dx = (a1*dx + v_r1 - a0)*2/dx;
275 /* Set dc to the derivative at the next point */
278 if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
288 table_f[(i-1)] = -0.5*dc;
290 /* Currently the last value only contains half the force: double it */
293 if (table_v != NULL && table_fdv0 != NULL)
295 /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
296 * init_ewald_f_table().
298 for (i = 0; i < ntab-1; i++)
300 table_fdv0[4*i] = table_f[i];
301 table_fdv0[4*i+1] = table_f[i+1]-table_f[i];
302 table_fdv0[4*i+2] = table_v[i];
303 table_fdv0[4*i+3] = 0.0;
305 table_fdv0[4*(ntab-1)] = table_f[(ntab-1)];
306 table_fdv0[4*(ntab-1)+1] = -table_f[(ntab-1)];
307 table_fdv0[4*(ntab-1)+2] = table_v[(ntab-1)];
308 table_fdv0[4*(ntab-1)+3] = 0.0;
312 /* The scale (1/spacing) for third order spline interpolation
313 * of the Ewald mesh contribution which needs to be subtracted
314 * from the non-bonded interactions.
316 real ewald_spline3_table_scale(real ewaldcoeff, real rc)
318 double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
322 /* Force tolerance: single precision accuracy */
323 ftol = GMX_FLOAT_EPS;
324 sc_f = sqrt(erf_x_d3/(6*4*ftol*ewaldcoeff))*ewaldcoeff;
326 /* Energy tolerance: 10x more accurate than the cut-off jump */
327 etol = 0.1*gmx_erfc(ewaldcoeff*rc);
328 etol = max(etol, GMX_REAL_EPS);
329 sc_e = pow(erf_x_d3/(6*12*sqrt(3)*etol), 1.0/3.0)*ewaldcoeff;
331 return max(sc_f, sc_e);
334 /* Calculate the potential and force for an r value
335 * in exactly the same way it is done in the inner loop.
336 * VFtab is a pointer to the table data, offset is
337 * the point where we should begin and stride is
338 * 4 if we have a buckingham table, 3 otherwise.
339 * If you want to evaluate table no N, set offset to 4*N.
341 * We use normal precision here, since that is what we
342 * will use in the inner loops.
344 static void evaluate_table(real VFtab[], int offset, int stride,
345 real tabscale, real r, real *y, real *yp)
349 real Y, F, Geps, Heps2, Fp;
358 Geps = eps*VFtab[n+2];
359 Heps2 = eps2*VFtab[n+3];
362 *yp = (Fp+Geps+2.0*Heps2)*tabscale;
365 static void copy2table(int n, int offset, int stride,
366 double x[], double Vtab[], double Ftab[], real scalefactor,
369 /* Use double prec. for the intermediary variables
370 * and temporary x/vtab/vtab2 data to avoid unnecessary
377 for (i = 0; (i < n); i++)
383 G = 3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h;
384 H = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] + Ftab[i])*h;
388 /* Fill the last entry with a linear potential,
389 * this is mainly for rounding issues with angle and dihedral potentials.
395 nn0 = offset + i*stride;
396 dest[nn0] = scalefactor*Vtab[i];
397 dest[nn0+1] = scalefactor*F;
398 dest[nn0+2] = scalefactor*G;
399 dest[nn0+3] = scalefactor*H;
403 static void init_table(int n, int nx0,
404 double tabscale, t_tabledata *td, gmx_bool bAlloc)
410 td->tabscale = tabscale;
417 for (i = 0; (i < td->nx); i++)
419 td->x[i] = i/tabscale;
423 static void spline_forces(int nx, double h, double v[], gmx_bool bS3, gmx_bool bE3,
427 double v3, b_s, b_e, b;
430 /* Formulas can be found in:
431 * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
434 if (nx < 4 && (bS3 || bE3))
436 gmx_fatal(FARGS, "Can not generate splines with third derivative boundary conditions with less than 4 (%d) points", nx);
439 /* To make life easy we initially set the spacing to 1
440 * and correct for this at the end.
445 /* Fit V''' at the start */
446 v3 = v[3] - 3*v[2] + 3*v[1] - v[0];
449 fprintf(debug, "The left third derivative is %g\n", v3/(h*h*h));
451 b_s = 2*(v[1] - v[0]) + v3/6;
456 /* Fit V'' at the start */
459 v2 = -v[3] + 4*v[2] - 5*v[1] + 2*v[0];
460 /* v2 = v[2] - 2*v[1] + v[0]; */
463 fprintf(debug, "The left second derivative is %g\n", v2/(h*h));
465 b_s = 3*(v[1] - v[0]) - v2/2;
471 b_s = 3*(v[2] - v[0]) + f[0]*h;
476 /* Fit V''' at the end */
477 v3 = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4];
480 fprintf(debug, "The right third derivative is %g\n", v3/(h*h*h));
482 b_e = 2*(v[nx-1] - v[nx-2]) + v3/6;
487 /* V'=0 at the end */
488 b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h;
493 beta = (bS3 ? 1 : 4);
495 /* For V'' fitting */
496 /* beta = (bS3 ? 2 : 4); */
499 for (i = start+1; i < end; i++)
503 b = 3*(v[i+1] - v[i-1]);
504 f[i] = (b - f[i-1])/beta;
506 gamma[end-1] = 1/beta;
507 beta = (bE3 ? 1 : 4) - gamma[end-1];
508 f[end-1] = (b_e - f[end-2])/beta;
510 for (i = end-2; i >= start; i--)
512 f[i] -= gamma[i+1]*f[i+1];
516 /* Correct for the minus sign and the spacing */
517 for (i = start; i < end; i++)
523 static void set_forces(FILE *fp, int angle,
524 int nx, double h, double v[], double f[],
532 "Force generation for dihedral tables is not (yet) implemented");
536 while (v[start] == 0)
542 while (v[end-1] == 0)
557 fprintf(fp, "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
558 table+1, start*h, end == nx ? "V'''" : "V'=0", (end-1)*h);
560 spline_forces(end-start, h, v+start, TRUE, end == nx, f+start);
563 static void read_tables(FILE *fp, const char *fn,
564 int ntab, int angle, t_tabledata td[])
568 double **yy = NULL, start, end, dx0, dx1, ssd, vm, vp, f, numf;
569 int k, i, nx, nx0 = 0, ny, nny, ns;
570 gmx_bool bAllZero, bZeroV, bZeroF;
574 libfn = gmxlibfn(fn);
575 nx = read_xvg(libfn, &yy, &ny);
578 gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d",
586 "The first distance in file %s is %f nm instead of %f nm",
587 libfn, yy[0][0], 0.0);
601 if (yy[0][0] != start || yy[0][nx-1] != end)
603 gmx_fatal(FARGS, "The angles in file %s should go from %f to %f instead of %f to %f\n",
604 libfn, start, end, yy[0][0], yy[0][nx-1]);
608 tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
612 fprintf(fp, "Read user tables from %s with %d data points.\n", libfn, nx);
615 fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
620 for (k = 0; k < ntab; k++)
624 for (i = 0; (i < nx); i++)
628 dx0 = yy[0][i-1] - yy[0][i-2];
629 dx1 = yy[0][i] - yy[0][i-1];
630 /* Check for 1% deviation in spacing */
631 if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1)))
633 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]);
636 if (yy[1+k*2][i] != 0)
644 if (yy[1+k*2][i] > 0.01*GMX_REAL_MAX ||
645 yy[1+k*2][i] < -0.01*GMX_REAL_MAX)
647 gmx_fatal(FARGS, "Out of range potential value %g in file '%s'",
651 if (yy[1+k*2+1][i] != 0)
659 if (yy[1+k*2+1][i] > 0.01*GMX_REAL_MAX ||
660 yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX)
662 gmx_fatal(FARGS, "Out of range force value %g in file '%s'",
668 if (!bZeroV && bZeroF)
670 set_forces(fp, angle, nx, 1/tabscale, yy[1+k*2], yy[1+k*2+1], k);
674 /* Check if the second column is close to minus the numerical
675 * derivative of the first column.
679 for (i = 1; (i < nx-1); i++)
684 if (vm != 0 && vp != 0 && f != 0)
686 /* Take the centered difference */
687 numf = -(vp - vm)*0.5*tabscale;
688 ssd += fabs(2*(f - numf)/(f + numf));
695 sprintf(buf, "For the %d non-zero entries for table %d in %s the forces deviate on average %d%% from minus the numerical derivative of the potential\n", ns, k, libfn, (int)(100*ssd+0.5));
698 fprintf(debug, "%s", buf);
704 fprintf(fp, "\nWARNING: %s\n", buf);
706 fprintf(stderr, "\nWARNING: %s\n", buf);
713 fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn);
716 for (k = 0; (k < ntab); k++)
718 init_table(nx, nx0, tabscale, &(td[k]), TRUE);
719 for (i = 0; (i < nx); i++)
721 td[k].x[i] = yy[0][i];
722 td[k].v[i] = yy[2*k+1][i];
723 td[k].f[i] = yy[2*k+2][i];
726 for (i = 0; (i < ny); i++)
734 static void done_tabledata(t_tabledata *td)
748 static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
751 /* Fill the table according to the formulas in the manual.
752 * In principle, we only need the potential and the second
753 * derivative, but then we would have to do lots of calculations
754 * in the inner loop. By precalculating some terms (see manual)
755 * we get better eventual performance, despite a larger table.
757 * Since some of these higher-order terms are very small,
758 * we always use double precision to calculate them here, in order
759 * to avoid unnecessary loss of precision.
766 double r1, rc, r12, r13;
767 double r, r2, r6, rc2, rc6, rc12;
768 double expr, Vtab, Ftab;
769 /* Parameters for David's function */
770 double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
771 /* Parameters for the switching function */
772 double ksw, swi, swi1;
773 /* Temporary parameters */
774 gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
775 double ewc = fr->ewaldcoeff_q;
776 double ewclj = fr->ewaldcoeff_lj;
781 bPotentialSwitch = FALSE;
782 bForceSwitch = FALSE;
783 bPotentialShift = FALSE;
787 bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
788 (tp == etabCOULSwitch) ||
789 (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
790 (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
791 (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
792 bForceSwitch = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
794 (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodFORCESWITCH)) ||
795 (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodFORCESWITCH)));
796 bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
797 (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
802 if (tprops[tp].bCoulomb)
804 r1 = fr->rcoulomb_switch;
809 r1 = fr->rvdw_switch;
812 if (bPotentialSwitch)
814 ksw = 1.0/(pow5(rc-r1));
826 else if (tp == etabLJ6Shift)
835 A = p * ((p+1)*r1-(p+4)*rc)/(pow(rc, p+2)*pow2(rc-r1));
836 B = -p * ((p+1)*r1-(p+3)*rc)/(pow(rc, p+2)*pow3(rc-r1));
837 C = 1.0/pow(rc, p)-A/3.0*pow3(rc-r1)-B/4.0*pow4(rc-r1);
838 if (tp == etabLJ6Shift)
849 fprintf(debug, "Setting up tables\n"); fflush(debug);
853 fp = xvgropen("switch.xvg", "switch", "r", "s");
859 rc6 = 1.0/(rc2*rc2*rc2);
860 if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
866 rc12 = pow(rc, -reppow);
876 Vcut = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + pow4(ewclj)*rc2*rc2/2);
886 case etabEwaldSwitch:
887 Vtab = gmx_erfc(ewc*rc)/rc;
890 /* Only calculate minus the reciprocal space contribution */
891 Vtab = -gmx_erf(ewc*rc)/rc;
895 /* No need for preventing the usage of modifiers with RF */
902 gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
903 tprops[tp].name, __FILE__, __LINE__);
907 for (i = td->nx0; (i < td->nx); i++)
912 if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
918 r12 = pow(r, -reppow);
922 if (bPotentialSwitch)
924 /* swi is function, swi1 1st derivative and swi2 2nd derivative */
925 /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
926 * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
928 * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
942 swi = 1 - 10*pow3(r-r1)*ksw*pow2(rc-r1)
943 + 15*pow4(r-r1)*ksw*(rc-r1) - 6*pow5(r-r1)*ksw;
944 swi1 = -30*pow2(r-r1)*ksw*pow2(rc-r1)
945 + 60*pow3(r-r1)*ksw*(rc-r1) - 30*pow4(r-r1)*ksw;
948 else /* not really needed, but avoids compiler warnings... */
954 fprintf(fp, "%10g %10g %10g %10g\n", r, swi, swi1, swi2);
980 Ftab = reppow*Vtab/r;
988 Ftab = reppow*Vtab/r;
994 Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
995 Ftab = -(6.0*r6/r-6.0*rc6/rc);
1006 Vtab = -(r6-6.0*(rc-r)*rc6/rc-rc6);
1007 Ftab = -(6.0*r6/r-6.0*rc6/rc);
1019 case etabCOULSwitch:
1028 case etabEwaldSwitch:
1029 Vtab = gmx_erfc(ewc*r)/r;
1030 Ftab = gmx_erfc(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
1033 case etabEwaldUserSwitch:
1034 /* Only calculate the negative of the reciprocal space contribution */
1035 Vtab = -gmx_erf(ewc*r)/r;
1036 Ftab = -gmx_erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
1039 Vtab = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + pow4(ewclj)*r2*r2/2);
1040 Ftab = 6.0*Vtab/r - r6*exp(-ewclj*ewclj*r2)*pow5(ewclj)*ewclj*r2*r2*r;
1044 Vtab = 1.0/r + fr->k_rf*r2 - fr->c_rf;
1045 Ftab = 1.0/r2 - 2*fr->k_rf*r;
1046 if (tp == etabRF_ZERO && r >= rc)
1060 Vtab = 1.0/r-(rc-r)/(rc*rc)-1.0/rc;
1061 Ftab = 1.0/r2-1.0/(rc*rc);
1070 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
1071 tp, __FILE__, __LINE__);
1075 /* Normal coulomb with cut-off correction for potential */
1079 /* If in Shifting range add something to it */
1082 r12 = (r-r1)*(r-r1);
1084 Vtab += -A_3*r13 - B_4*r12*r12;
1085 Ftab += A*r12 + B*r13;
1090 /* Make sure interactions are zero outside cutoff with modifiers */
1095 if (bPotentialShift)
1103 /* Make sure interactions are zero outside cutoff with modifiers */
1115 if (bPotentialSwitch)
1119 /* Make sure interactions are zero outside cutoff with modifiers */
1125 Ftab = Ftab*swi - Vtab*swi1;
1129 /* Convert to single precision when we store to mem */
1134 /* Continue the table linearly from nx0 to 0.
1135 * These values are only required for energy minimization with overlap or TPI.
1137 for (i = td->nx0-1; i >= 0; i--)
1139 td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
1140 td->f[i] = td->f[i+1];
1148 static void set_table_type(int tabsel[], const t_forcerec *fr, gmx_bool b14only)
1150 int eltype, vdwtype;
1152 /* Set the different table indices.
1159 switch (fr->eeltype)
1166 case eelPMEUSERSWITCH:
1175 eltype = fr->eeltype;
1181 tabsel[etiCOUL] = etabCOUL;
1184 tabsel[etiCOUL] = etabShift;
1187 if (fr->rcoulomb > fr->rcoulomb_switch)
1189 tabsel[etiCOUL] = etabShift;
1193 tabsel[etiCOUL] = etabCOUL;
1199 tabsel[etiCOUL] = etabEwald;
1202 tabsel[etiCOUL] = etabEwaldSwitch;
1205 tabsel[etiCOUL] = etabEwaldUser;
1207 case eelPMEUSERSWITCH:
1208 tabsel[etiCOUL] = etabEwaldUserSwitch;
1213 tabsel[etiCOUL] = etabRF;
1216 tabsel[etiCOUL] = etabRF_ZERO;
1219 tabsel[etiCOUL] = etabCOULSwitch;
1222 tabsel[etiCOUL] = etabUSER;
1225 tabsel[etiCOUL] = etabCOULEncad;
1228 gmx_fatal(FARGS, "Invalid eeltype %d", eltype);
1231 /* Van der Waals time */
1232 if (fr->bBHAM && !b14only)
1234 tabsel[etiLJ6] = etabLJ6;
1235 tabsel[etiLJ12] = etabEXPMIN;
1239 if (b14only && fr->vdwtype != evdwUSER)
1245 vdwtype = fr->vdwtype;
1251 tabsel[etiLJ6] = etabLJ6Switch;
1252 tabsel[etiLJ12] = etabLJ12Switch;
1255 tabsel[etiLJ6] = etabLJ6Shift;
1256 tabsel[etiLJ12] = etabLJ12Shift;
1259 tabsel[etiLJ6] = etabUSER;
1260 tabsel[etiLJ12] = etabUSER;
1263 tabsel[etiLJ6] = etabLJ6;
1264 tabsel[etiLJ12] = etabLJ12;
1266 case evdwENCADSHIFT:
1267 tabsel[etiLJ6] = etabLJ6Encad;
1268 tabsel[etiLJ12] = etabLJ12Encad;
1271 tabsel[etiLJ6] = etabLJ6Ewald;
1272 tabsel[etiLJ12] = etabLJ12;
1275 gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype,
1276 __FILE__, __LINE__);
1279 if (!b14only && fr->vdw_modifier != eintmodNONE)
1281 if (fr->vdw_modifier != eintmodPOTSHIFT &&
1282 fr->vdwtype != evdwCUT)
1284 gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
1287 /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
1288 * to the original interaction forms when we fill the table, so we only check cutoffs here.
1290 if (fr->vdwtype == evdwCUT)
1292 switch (fr->vdw_modifier)
1295 case eintmodPOTSHIFT:
1296 case eintmodEXACTCUTOFF:
1297 /* No modification */
1299 case eintmodPOTSWITCH:
1300 tabsel[etiLJ6] = etabLJ6Switch;
1301 tabsel[etiLJ12] = etabLJ12Switch;
1303 case eintmodFORCESWITCH:
1304 tabsel[etiLJ6] = etabLJ6Shift;
1305 tabsel[etiLJ12] = etabLJ12Shift;
1308 gmx_incons("Unsupported vdw_modifier");
1315 t_forcetable make_tables(FILE *out, const output_env_t oenv,
1316 const t_forcerec *fr,
1317 gmx_bool bVerbose, const char *fn,
1318 real rtab, int flags)
1320 const char *fns[3] = { "ctab.xvg", "dtab.xvg", "rtab.xvg" };
1321 const char *fns14[3] = { "ctab14.xvg", "dtab14.xvg", "rtab14.xvg" };
1324 gmx_bool b14only, bReadTab, bGenTab;
1326 int i, j, k, nx, nx0, tabsel[etiNR];
1331 b14only = (flags & GMX_MAKETABLES_14ONLY);
1333 if (flags & GMX_MAKETABLES_FORCEUSER)
1335 tabsel[etiCOUL] = etabUSER;
1336 tabsel[etiLJ6] = etabUSER;
1337 tabsel[etiLJ12] = etabUSER;
1341 set_table_type(tabsel, fr, b14only);
1347 table.scale_exp = 0;
1351 table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
1352 table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1353 table.formatsize = 4;
1354 table.ninteractions = 3;
1355 table.stride = table.formatsize*table.ninteractions;
1357 /* Check whether we have to read or generate */
1360 for (i = 0; (i < etiNR); i++)
1362 if (ETAB_USER(tabsel[i]))
1366 if (tabsel[i] != etabUSER)
1373 read_tables(out, fn, etiNR, 0, td);
1374 if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1376 rtab = td[0].x[td[0].nx-1];
1382 if (td[0].x[td[0].nx-1] < rtab)
1384 gmx_fatal(FARGS, "Tables in file %s not long enough for cut-off:\n"
1385 "\tshould be at least %f nm\n", fn, rtab);
1387 nx = table.n = (int)(rtab*td[0].tabscale + 0.5);
1389 table.scale = td[0].tabscale;
1397 table.scale = 2000.0;
1399 table.scale = 500.0;
1401 nx = table.n = rtab*table.scale;
1406 if (fr->bham_b_max != 0)
1408 table.scale_exp = table.scale/fr->bham_b_max;
1412 table.scale_exp = table.scale;
1416 /* Each table type (e.g. coul,lj6,lj12) requires four
1417 * numbers per nx+1 data points. For performance reasons we want
1418 * the table data to be aligned to 16-byte.
1420 snew_aligned(table.data, 12*(nx+1)*sizeof(real), 32);
1422 for (k = 0; (k < etiNR); k++)
1424 if (tabsel[k] != etabUSER)
1427 (tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
1428 &(td[k]), !bReadTab);
1429 fill_table(&(td[k]), tabsel[k], fr, b14only);
1432 fprintf(out, "%s table with %d data points for %s%s.\n"
1433 "Tabscale = %g points/nm\n",
1434 ETAB_USER(tabsel[k]) ? "Modified" : "Generated",
1435 td[k].nx, b14only ? "1-4 " : "", tprops[tabsel[k]].name,
1440 /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
1441 * by including the derivative constants (6.0 or 12.0) in the parameters, since
1442 * we no longer calculate force in most steps. This means the c6/c12 parameters
1443 * have been scaled up, so we need to scale down the table interactions too.
1444 * It comes here since we need to scale user tables too.
1448 scalefactor = 1.0/6.0;
1450 else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1452 scalefactor = 1.0/12.0;
1459 copy2table(table.n, k*4, 12, td[k].x, td[k].v, td[k].f, scalefactor, table.data);
1461 if (bDebugMode() && bVerbose)
1465 fp = xvgropen(fns14[k], fns14[k], "r", "V", oenv);
1469 fp = xvgropen(fns[k], fns[k], "r", "V", oenv);
1471 /* plot the output 5 times denser than the table data */
1472 for (i = 5*((nx0+1)/2); i < 5*table.n; i++)
1474 x0 = i*table.r/(5*(table.n-1));
1475 evaluate_table(table.data, 4*k, 12, table.scale, x0, &y0, &yp);
1476 fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp);
1480 done_tabledata(&(td[k]));
1487 t_forcetable make_gb_table(const output_env_t oenv,
1488 const t_forcerec *fr)
1490 const char *fns[3] = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" };
1491 const char *fns14[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" };
1494 gmx_bool bReadTab, bGenTab;
1496 int i, j, k, nx, nx0, tabsel[etiNR];
1497 double r, r2, Vtab, Ftab, expterm;
1501 double abs_error_r, abs_error_r2;
1502 double rel_error_r, rel_error_r2;
1503 double rel_error_r_old = 0, rel_error_r2_old = 0;
1504 double x0_r_error, x0_r2_error;
1507 /* Only set a Coulomb table for GB */
1514 /* Set the table dimensions for GB, not really necessary to
1515 * use etiNR (since we only have one table, but ...)
1518 table.interaction = GMX_TABLE_INTERACTION_ELEC;
1519 table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1520 table.r = fr->gbtabr;
1521 table.scale = fr->gbtabscale;
1522 table.scale_exp = 0;
1523 table.n = table.scale*table.r;
1524 table.formatsize = 4;
1525 table.ninteractions = 1;
1526 table.stride = table.formatsize*table.ninteractions;
1528 nx = table.scale*table.r;
1530 /* Check whether we have to read or generate
1531 * We will always generate a table, so remove the read code
1532 * (Compare with original make_table function
1537 /* Each table type (e.g. coul,lj6,lj12) requires four
1538 * numbers per datapoint. For performance reasons we want
1539 * the table data to be aligned to 16-byte. This is accomplished
1540 * by allocating 16 bytes extra to a temporary pointer, and then
1541 * calculating an aligned pointer. This new pointer must not be
1542 * used in a free() call, but thankfully we're sloppy enough not
1546 snew_aligned(table.data, 4*nx, 32);
1548 init_table(nx, nx0, table.scale, &(td[0]), !bReadTab);
1550 /* Local implementation so we don't have to use the etabGB
1551 * enum above, which will cause problems later when
1552 * making the other tables (right now even though we are using
1553 * GB, the normal Coulomb tables will be created, but this
1554 * will cause a problem since fr->eeltype==etabGB which will not
1555 * be defined in fill_table and set_table_type
1558 for (i = nx0; i < nx; i++)
1562 expterm = exp(-0.25*r2);
1564 Vtab = 1/sqrt(r2+expterm);
1565 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1567 /* Convert to single precision when we store to mem */
1573 copy2table(table.n, 0, 4, td[0].x, td[0].v, td[0].f, 1.0, table.data);
1577 fp = xvgropen(fns[0], fns[0], "r", "V", oenv);
1578 /* plot the output 5 times denser than the table data */
1579 /* for(i=5*nx0;i<5*table.n;i++) */
1580 for (i = nx0; i < table.n; i++)
1582 /* x0=i*table.r/(5*table.n); */
1583 x0 = i*table.r/table.n;
1584 evaluate_table(table.data, 0, 4, table.scale, x0, &y0, &yp);
1585 fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp);
1592 for(i=100*nx0;i<99.81*table.n;i++)
1594 r = i*table.r/(100*table.n);
1596 expterm = exp(-0.25*r2);
1598 Vtab = 1/sqrt(r2+expterm);
1599 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1602 evaluate_table(table.data,0,4,table.scale,r,&y0,&yp);
1603 printf("gb: i=%d, x0=%g, y0=%15.15f, Vtab=%15.15f, yp=%15.15f, Ftab=%15.15f\n",i,r, y0, Vtab, yp, Ftab);
1605 abs_error_r=fabs(y0-Vtab);
1606 abs_error_r2=fabs(yp-(-1)*Ftab);
1608 rel_error_r=abs_error_r/y0;
1609 rel_error_r2=fabs(abs_error_r2/yp);
1612 if(rel_error_r>rel_error_r_old)
1614 rel_error_r_old=rel_error_r;
1618 if(rel_error_r2>rel_error_r2_old)
1620 rel_error_r2_old=rel_error_r2;
1625 printf("gb: MAX REL ERROR IN R=%15.15f, MAX REL ERROR IN R2=%15.15f\n",rel_error_r_old, rel_error_r2_old);
1626 printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error);
1629 done_tabledata(&(td[0]));
1637 t_forcetable make_atf_table(FILE *out, const output_env_t oenv,
1638 const t_forcerec *fr,
1642 const char *fns[3] = { "tf_tab.xvg", "atfdtab.xvg", "atfrtab.xvg" };
1645 real x0, y0, yp, rtab;
1647 real rx, ry, rz, box_r;
1652 /* Set the table dimensions for ATF, not really necessary to
1653 * use etiNR (since we only have one table, but ...)
1657 if (fr->adress_type == eAdressSphere)
1659 /* take half box diagonal direction as tab range */
1660 rx = 0.5*box[0][0]+0.5*box[1][0]+0.5*box[2][0];
1661 ry = 0.5*box[0][1]+0.5*box[1][1]+0.5*box[2][1];
1662 rz = 0.5*box[0][2]+0.5*box[1][2]+0.5*box[2][2];
1663 box_r = sqrt(rx*rx+ry*ry+rz*rz);
1668 /* xsplit: take half box x direction as tab range */
1669 box_r = box[0][0]/2;
1674 table.scale_exp = 0;
1678 read_tables(out, fn, 1, 0, td);
1679 rtab = td[0].x[td[0].nx-1];
1681 if (fr->adress_type == eAdressXSplit && (rtab < box[0][0]/2))
1683 gmx_fatal(FARGS, "AdResS full box therm force table in file %s extends to %f:\n"
1684 "\tshould extend to at least half the length of the box in x-direction"
1685 "%f\n", fn, rtab, box[0][0]/2);
1689 gmx_fatal(FARGS, "AdResS full box therm force table in file %s extends to %f:\n"
1690 "\tshould extend to at least for spherical adress"
1691 "%f (=distance from center to furthermost point in box \n", fn, rtab, box_r);
1697 table.scale = td[0].tabscale;
1700 /* Each table type (e.g. coul,lj6,lj12) requires four
1701 * numbers per datapoint. For performance reasons we want
1702 * the table data to be aligned to 16-byte. This is accomplished
1703 * by allocating 16 bytes extra to a temporary pointer, and then
1704 * calculating an aligned pointer. This new pointer must not be
1705 * used in a free() call, but thankfully we're sloppy enough not
1709 snew_aligned(table.data, 4*nx, 32);
1711 copy2table(table.n, 0, 4, td[0].x, td[0].v, td[0].f, 1.0, table.data);
1715 fp = xvgropen(fns[0], fns[0], "r", "V", oenv);
1716 /* plot the output 5 times denser than the table data */
1717 /* for(i=5*nx0;i<5*table.n;i++) */
1719 for (i = 5*((nx0+1)/2); i < 5*table.n; i++)
1721 /* x0=i*table.r/(5*table.n); */
1722 x0 = i*table.r/(5*(table.n-1));
1723 evaluate_table(table.data, 0, 4, table.scale, x0, &y0, &yp);
1724 fprintf(fp, "%15.10e %15.10e %15.10e\n", x0, y0, yp);
1730 done_tabledata(&(td[0]));
1733 table.interaction = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
1734 table.format = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1735 table.formatsize = 4;
1736 table.ninteractions = 3;
1737 table.stride = table.formatsize*table.ninteractions;
1743 bondedtable_t make_bonded_table(FILE *fplog, char *fn, int angle)
1758 read_tables(fplog, fn, 1, angle, &td);
1761 /* Convert the table from degrees to radians */
1762 for (i = 0; i < td.nx; i++)
1767 td.tabscale *= RAD2DEG;
1770 tab.scale = td.tabscale;
1771 snew(tab.data, tab.n*4);
1772 copy2table(tab.n, 0, 4, td.x, td.v, td.f, 1.0, tab.data);
1773 done_tabledata(&td);