src/gromacs/tables/forcetable.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
  10  *
  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.
  15  *
  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.
  20  *
  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.
  25  *
  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.
  33  *
  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.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "forcetable.h"
  40
  41 #include <cmath>
  42
  43 #include <algorithm>
  44
  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"
  58
  59 /* All the possible (implemented) table functions */
  60 enum
  61 {
  62     etabLJ6,
  63     etabLJ12,
  64     etabLJ6Shift,
  65     etabLJ12Shift,
  66     etabShift,
  67     etabRF,
  68     etabRF_ZERO,
  69     etabCOUL,
  70     etabEwald,
  71     etabEwaldSwitch,
  72     etabEwaldUser,
  73     etabEwaldUserSwitch,
  74     etabLJ6Ewald,
  75     etabLJ6Switch,
  76     etabLJ12Switch,
  77     etabCOULSwitch,
  78     etabLJ6Encad,
  79     etabLJ12Encad,
  80     etabCOULEncad,
  81     etabEXPMIN,
  82     etabUSER,
  83     etabNR
  84 };
  85
  86 /** Evaluates to true if the table type contains user data. */
  87 #define ETAB_USER(e) ((e) == etabUSER || (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
  88
  89 typedef struct
  90 {
  91     const char* name;
  92     gmx_bool    bCoulomb;
  93 } t_tab_props;
  94
  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] = {
  98     { "LJ6", FALSE },
  99     { "LJ12", FALSE },
 100     { "LJ6Shift", FALSE },
 101     { "LJ12Shift", FALSE },
 102     { "Shift", TRUE },
 103     { "RF", TRUE },
 104     { "RF-zero", TRUE },
 105     { "COUL", TRUE },
 106     { "Ewald", TRUE },
 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 },
 117     { "EXPMIN", FALSE },
 118     { "USER", FALSE },
 119 };
 120
 121 typedef struct
 122 {
 123     int     nx, nx0;
 124     double  tabscale;
 125     double *x, *v, *f;
 126 } t_tabledata;
 127
 128 double v_q_ewald_lr(double beta, double r)
 129 {
 130     if (r == 0)
 131     {
 132         return beta * 2 / sqrt(M_PI);
 133     }
 134     else
 135     {
 136         return std::erf(beta * r) / r;
 137     }
 138 }
 139
 140 double v_lj_ewald_lr(double beta, double r)
 141 {
 142     double br, br2, br4, r6, factor;
 143     if (r == 0)
 144     {
 145         return gmx::power6(beta) / 6;
 146     }
 147     else
 148     {
 149         br     = beta * r;
 150         br2    = br * br;
 151         br4    = br2 * br2;
 152         r6     = gmx::power6(r);
 153         factor = (1.0 - std::exp(-br2) * (1 + br2 + 0.5 * br4)) / r6;
 154         return factor;
 155     }
 156 }
 157
 158 EwaldCorrectionTables generateEwaldCorrectionTables(const int    numPoints,
 159                                                     const double tableScaling,
 160                                                     const real   beta,
 161                                                     real_space_grid_contribution_computer v_lr)
 162 {
 163     real     tab_max;
 164     int      i, i_inrange;
 165     double   dc, dc_new;
 166     gmx_bool bOutOfRange;
 167     double   v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
 168     double   x_r0;
 169
 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.
 172      */
 173
 174     if (numPoints < 2)
 175     {
 176         gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
 177     }
 178
 179     const double dx = 1 / tableScaling;
 180
 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;
 189
 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.
 193      */
 194     tab_max = GMX_REAL_MAX * 0.0001;
 195
 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.
 200      */
 201
 202     bOutOfRange = FALSE;
 203     i_inrange   = numPoints;
 204     v_inrange   = 0;
 205     dc          = 0;
 206     for (i = numPoints - 1; i >= 0; i--)
 207     {
 208         x_r0 = i * dx;
 209
 210         v_r0 = (*v_lr)(beta, x_r0);
 211
 212         if (!bOutOfRange)
 213         {
 214             i_inrange = i;
 215             v_inrange = v_r0;
 216
 217             vi = v_r0;
 218         }
 219         else
 220         {
 221             /* Linear continuation for the last point in range */
 222             vi = v_inrange - dc * (i - i_inrange) * dx;
 223         }
 224
 225         table_v[i] = vi;
 226
 227         if (i == 0)
 228         {
 229             continue;
 230         }
 231
 232         /* Get the potential at table point i-1 */
 233         v_r1 = (*v_lr)(beta, (i - 1) * dx);
 234
 235         if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
 236         {
 237             bOutOfRange = TRUE;
 238         }
 239
 240         if (!bOutOfRange)
 241         {
 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.
 244              */
 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.
 249              */
 250             dc = (v_r0 - v_r1) / dx + 0.5 * a2dx;
 251         }
 252
 253         if (i == numPoints - 1)
 254         {
 255             /* Fill the table with the force, minus the derivative of the spline */
 256             table_f[i] = -dc;
 257         }
 258         else
 259         {
 260             /* tab[i] will contain the average of the splines over the two intervals */
 261             table_f[i] += -0.5 * dc;
 262         }
 263
 264         if (!bOutOfRange)
 265         {
 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.
 269              */
 270             a0   = v_r0;
 271             a1   = dc;
 272             a2dx = (a1 * dx + v_r1 - a0) * 2 / dx;
 273
 274             /* Set dc to the derivative at the next point */
 275             dc_new = a1 - a2dx;
 276
 277             if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
 278             {
 279                 bOutOfRange = TRUE;
 280             }
 281             else
 282             {
 283                 dc = dc_new;
 284             }
 285         }
 286
 287         table_f[(i - 1)] = -0.5 * dc;
 288     }
 289     /* Currently the last value only contains half the force: double it */
 290     table_f[0] *= 2;
 291
 292     if (!table_fdv0.empty())
 293     {
 294         /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
 295          * init_ewald_f_table().
 296          */
 297         for (i = 0; i < numPoints - 1; i++)
 298         {
 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;
 303         }
 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;
 309     }
 310
 311     return tables;
 312 }
 313
 314 /* Returns the spacing for a function using the maximum of
 315  * the third derivative, x_scale (unit 1/length)
 316  * and function tolerance.
 317  */
 318 static double spline3_table_scale(double third_deriv_max, double x_scale, double func_tol)
 319 {
 320     double deriv_tol;
 321     double sc_deriv, sc_func;
 322
 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;
 326
 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;
 330
 331     return std::max(sc_deriv, sc_func);
 332 }
 333
 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.
 339  *
 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).
 345  */
 346 real ewald_spline3_table_scale(const interaction_const_t& ic,
 347                                const bool                 generateCoulombTables,
 348                                const bool                 generateVdwTables)
 349 {
 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");
 354
 355     real sc = 0;
 356
 357     if (generateCoulombTables)
 358     {
 359         GMX_RELEASE_ASSERT(ic.ewaldcoeff_q > 0, "The Ewald coefficient shoule be positive");
 360
 361         double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
 362         double etol;
 363         real   sc_q;
 364
 365         /* Energy tolerance: 0.1 times the cut-off jump */
 366         etol = 0.1 * std::erfc(ic.ewaldcoeff_q * ic.rcoulomb);
 367
 368         sc_q = spline3_table_scale(erf_x_d3, ic.ewaldcoeff_q, etol);
 369
 370         if (debug)
 371         {
 372             fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f nm\n", 1 / sc_q);
 373         }
 374
 375         sc = std::max(sc, sc_q);
 376     }
 377
 378     if (generateVdwTables)
 379     {
 380         GMX_RELEASE_ASSERT(ic.ewaldcoeff_lj > 0, "The Ewald coefficient shoule be positive");
 381
 382         double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
 383         double xrc2, etol;
 384         real   sc_lj;
 385
 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);
 389
 390         sc_lj = spline3_table_scale(func_d3, ic.ewaldcoeff_lj, etol);
 391
 392         if (debug)
 393         {
 394             fprintf(debug, "Ewald LJ quadratic spline table spacing: %f nm\n", 1 / sc_lj);
 395         }
 396
 397         sc = std::max(sc, sc_lj);
 398     }
 399
 400     return sc;
 401 }
 402
 403 static void copy2table(int          n,
 404                        int          offset,
 405                        int          stride,
 406                        const double x[],
 407                        const double Vtab[],
 408                        const double Ftab[],
 409                        real         scalefactor,
 410                        real         dest[])
 411 {
 412     /* Use double prec. for the intermediary variables
 413      * and temporary x/vtab/vtab2 data to avoid unnecessary
 414      * loss of precision.
 415      */
 416     int    i, nn0;
 417     double F, G, H, h;
 418
 419     h = 0;
 420     for (i = 0; (i < n); i++)
 421     {
 422         if (i < n - 1)
 423         {
 424             h = x[i + 1] - x[i];
 425             F = -Ftab[i] * h;
 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;
 428         }
 429         else
 430         {
 431             /* Fill the last entry with a linear potential,
 432              * this is mainly for rounding issues with angle and dihedral potentials.
 433              */
 434             F = -Ftab[i] * h;
 435             G = 0;
 436             H = 0;
 437         }
 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;
 443     }
 444 }
 445
 446 static void init_table(int n, int nx0, double tabscale, t_tabledata* td, gmx_bool bAlloc)
 447 {
 448     td->nx       = n;
 449     td->nx0      = nx0;
 450     td->tabscale = tabscale;
 451     if (bAlloc)
 452     {
 453         snew(td->x, td->nx);
 454         snew(td->v, td->nx);
 455         snew(td->f, td->nx);
 456     }
 457 }
 458
 459 static void spline_forces(int nx, double h, const double v[], gmx_bool bS3, gmx_bool bE3, double f[])
 460 {
 461     int    start, end, i;
 462     double v3, b_s, b_e, b;
 463     double beta, *gamma;
 464
 465     /* Formulas can be found in:
 466      * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
 467      */
 468
 469     if (nx < 4 && (bS3 || bE3))
 470     {
 471         gmx_fatal(FARGS,
 472                   "Can not generate splines with third derivative boundary conditions with less "
 473                   "than 4 (%d) points",
 474                   nx);
 475     }
 476
 477     /* To make life easy we initially set the spacing to 1
 478      * and correct for this at the end.
 479      */
 480     if (bS3)
 481     {
 482         /* Fit V''' at the start */
 483         v3 = v[3] - 3 * v[2] + 3 * v[1] - v[0];
 484         if (debug)
 485         {
 486             fprintf(debug, "The left third derivative is %g\n", v3 / (h * h * h));
 487         }
 488         b_s   = 2 * (v[1] - v[0]) + v3 / 6;
 489         start = 0;
 490
 491         if (FALSE)
 492         {
 493             /* Fit V'' at the start */
 494             real v2;
 495
 496             v2 = -v[3] + 4 * v[2] - 5 * v[1] + 2 * v[0];
 497             /* v2  = v[2] - 2*v[1] + v[0]; */
 498             if (debug)
 499             {
 500                 fprintf(debug, "The left second derivative is %g\n", v2 / (h * h));
 501             }
 502             b_s   = 3 * (v[1] - v[0]) - v2 / 2;
 503             start = 0;
 504         }
 505     }
 506     else
 507     {
 508         b_s   = 3 * (v[2] - v[0]) + f[0] * h;
 509         start = 1;
 510     }
 511     if (bE3)
 512     {
 513         /* Fit V''' at the end */
 514         v3 = v[nx - 1] - 3 * v[nx - 2] + 3 * v[nx - 3] - v[nx - 4];
 515         if (debug)
 516         {
 517             fprintf(debug, "The right third derivative is %g\n", v3 / (h * h * h));
 518         }
 519         b_e = 2 * (v[nx - 1] - v[nx - 2]) + v3 / 6;
 520         end = nx;
 521     }
 522     else
 523     {
 524         /* V'=0 at the end */
 525         b_e = 3 * (v[nx - 1] - v[nx - 3]) + f[nx - 1] * h;
 526         end = nx - 1;
 527     }
 528
 529     snew(gamma, nx);
 530     beta = (bS3 ? 1 : 4);
 531
 532     /* For V'' fitting */
 533     /* beta = (bS3 ? 2 : 4); */
 534
 535     f[start] = b_s / beta;
 536     for (i = start + 1; i < end; i++)
 537     {
 538         gamma[i] = 1 / beta;
 539         beta     = 4 - gamma[i];
 540         b        = 3 * (v[i + 1] - v[i - 1]);
 541         f[i]     = (b - f[i - 1]) / beta;
 542     }
 543     gamma[end - 1] = 1 / beta;
 544     beta           = (bE3 ? 1 : 4) - gamma[end - 1];
 545     f[end - 1]     = (b_e - f[end - 2]) / beta;
 546
 547     for (i = end - 2; i >= start; i--)
 548     {
 549         f[i] -= gamma[i + 1] * f[i + 1];
 550     }
 551     sfree(gamma);
 552
 553     /* Correct for the minus sign and the spacing */
 554     for (i = start; i < end; i++)
 555     {
 556         f[i] = -f[i] / h;
 557     }
 558 }
 559
 560 static void set_forces(FILE* fp, int angle, int nx, double h, double v[], double f[], int table)
 561 {
 562     int start, end;
 563
 564     if (angle == 2)
 565     {
 566         gmx_fatal(FARGS, "Force generation for dihedral tables is not (yet) implemented");
 567     }
 568
 569     start = 0;
 570     while (v[start] == 0)
 571     {
 572         start++;
 573     }
 574
 575     end = nx;
 576     while (v[end - 1] == 0)
 577     {
 578         end--;
 579     }
 580     if (end > nx - 2)
 581     {
 582         end = nx;
 583     }
 584     else
 585     {
 586         end++;
 587     }
 588
 589     if (fp)
 590     {
 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);
 593     }
 594     spline_forces(end - start, h, v + start, TRUE, end == nx, f + start);
 595 }
 596
 597 static void read_tables(FILE* fp, const char* filename, int ntab, int angle, t_tabledata td[])
 598 {
 599     char     buf[STRLEN];
 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;
 603     double   tabscale;
 604
 605     nny               = 2 * ntab + 1;
 606     std::string libfn = gmx::findLibraryFile(filename);
 607     nx                = read_xvg(libfn.c_str(), &yy, &ny);
 608     if (ny != nny)
 609     {
 610         gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d", libfn.c_str(),
 611                   ny, nny);
 612     }
 613     if (angle == 0)
 614     {
 615         if (yy[0][0] != 0.0)
 616         {
 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);
 619         }
 620     }
 621     else
 622     {
 623         if (angle == 1)
 624         {
 625             start = 0.0;
 626         }
 627         else
 628         {
 629             start = -180.0;
 630         }
 631         end = 180.0;
 632         if (yy[0][0] != start || yy[0][nx - 1] != end)
 633         {
 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]);
 636         }
 637     }
 638
 639     tabscale = (nx - 1) / (yy[0][nx - 1] - yy[0][0]);
 640
 641     if (fp)
 642     {
 643         fprintf(fp, "Read user tables from %s with %d data points.\n", libfn.c_str(), nx);
 644         if (angle == 0)
 645         {
 646             fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
 647         }
 648     }
 649
 650     bAllZero = TRUE;
 651     for (k = 0; k < ntab; k++)
 652     {
 653         bZeroV = TRUE;
 654         bZeroF = TRUE;
 655         for (i = 0; (i < nx); i++)
 656         {
 657             if (i >= 2)
 658             {
 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)))
 663                 {
 664                     gmx_fatal(FARGS,
 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]);
 667                 }
 668             }
 669             if (yy[1 + k * 2][i] != 0)
 670             {
 671                 bZeroV = FALSE;
 672                 if (bAllZero)
 673                 {
 674                     bAllZero = FALSE;
 675                     nx0      = i;
 676                 }
 677                 if (yy[1 + k * 2][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2][i] < -0.01 * GMX_REAL_MAX)
 678                 {
 679                     gmx_fatal(FARGS, "Out of range potential value %g in file '%s'",
 680                               yy[1 + k * 2][i], filename);
 681                 }
 682             }
 683             if (yy[1 + k * 2 + 1][i] != 0)
 684             {
 685                 bZeroF = FALSE;
 686                 if (bAllZero)
 687                 {
 688                     bAllZero = FALSE;
 689                     nx0      = i;
 690                 }
 691                 if (yy[1 + k * 2 + 1][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2 + 1][i] < -0.01 * GMX_REAL_MAX)
 692                 {
 693                     gmx_fatal(FARGS, "Out of range force value %g in file '%s'",
 694                               yy[1 + k * 2 + 1][i], filename);
 695                 }
 696             }
 697         }
 698
 699         if (!bZeroV && bZeroF)
 700         {
 701             set_forces(fp, angle, nx, 1 / tabscale, yy[1 + k * 2], yy[1 + k * 2 + 1], k);
 702         }
 703         else
 704         {
 705             /* Check if the second column is close to minus the numerical
 706              * derivative of the first column.
 707              */
 708             ssd = 0;
 709             ns  = 0;
 710             for (i = 1; (i < nx - 1); i++)
 711             {
 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)
 716                 {
 717                     /* Take the centered difference */
 718                     numf = -(vp - vm) * 0.5 * tabscale;
 719                     if (f + numf != 0)
 720                     {
 721                         ssd += fabs(2 * (f - numf) / (f + numf));
 722                     }
 723                     ns++;
 724                 }
 725             }
 726             if (ns > 0)
 727             {
 728                 ssd /= ns;
 729                 sprintf(buf,
 730                         "For the %d non-zero entries for table %d in %s the forces deviate on "
 731                         "average %" PRId64
 732                         "%% from minus the numerical derivative of the potential\n",
 733                         ns, k, libfn.c_str(), gmx::roundToInt64(100 * ssd));
 734                 if (debug)
 735                 {
 736                     fprintf(debug, "%s", buf);
 737                 }
 738                 if (ssd > 0.2)
 739                 {
 740                     if (fp)
 741                     {
 742                         fprintf(fp, "\nWARNING: %s\n", buf);
 743                     }
 744                     fprintf(stderr, "\nWARNING: %s\n", buf);
 745                 }
 746             }
 747         }
 748     }
 749     if (bAllZero && fp)
 750     {
 751         fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn.c_str());
 752     }
 753
 754     for (k = 0; (k < ntab); k++)
 755     {
 756         init_table(nx, nx0, tabscale, &(td[k]), TRUE);
 757         for (i = 0; (i < nx); i++)
 758         {
 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];
 762         }
 763     }
 764     for (i = 0; (i < ny); i++)
 765     {
 766         sfree(yy[i]);
 767     }
 768     sfree(yy);
 769 }
 770
 771 static void done_tabledata(t_tabledata* td)
 772 {
 773     if (!td)
 774     {
 775         return;
 776     }
 777
 778     sfree(td->x);
 779     sfree(td->v);
 780     sfree(td->f);
 781 }
 782
 783 static void fill_table(t_tabledata* td, int tp, const interaction_const_t* ic, gmx_bool b14only)
 784 {
 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.
 790      *
 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.
 794      */
 795     int    i;
 796     double reppow, p;
 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;
 808     double   Vcut  = 0;
 809
 810     if (b14only)
 811     {
 812         bPotentialSwitch = FALSE;
 813         bForceSwitch     = FALSE;
 814         bPotentialShift  = FALSE;
 815     }
 816     else
 817     {
 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)));
 827     }
 828
 829     reppow = ic->reppow;
 830
 831     if (tprops[tp].bCoulomb)
 832     {
 833         r1 = ic->rcoulomb_switch;
 834         rc = ic->rcoulomb;
 835     }
 836     else
 837     {
 838         r1 = ic->rvdw_switch;
 839         rc = ic->rvdw;
 840     }
 841     if (bPotentialSwitch)
 842     {
 843         ksw = 1.0 / (gmx::power5(rc - r1));
 844     }
 845     else
 846     {
 847         ksw = 0.0;
 848     }
 849     if (bForceSwitch)
 850     {
 851         if (tp == etabShift)
 852         {
 853             p = 1;
 854         }
 855         else if (tp == etabLJ6Shift)
 856         {
 857             p = 6;
 858         }
 859         else
 860         {
 861             p = reppow;
 862         }
 863
 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)
 868         {
 869             A = -A;
 870             B = -B;
 871             C = -C;
 872         }
 873         A_3 = A / 3.0;
 874         B_4 = B / 4.0;
 875     }
 876     if (debug)
 877     {
 878         fprintf(debug, "Setting up tables\n");
 879         fflush(debug);
 880     }
 881
 882     if (bPotentialShift)
 883     {
 884         rc2 = rc * rc;
 885         rc6 = 1.0 / (rc2 * rc2 * rc2);
 886         if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 887         {
 888             rc12 = rc6 * rc6;
 889         }
 890         else
 891         {
 892             rc12 = std::pow(rc, -reppow);
 893         }
 894
 895         switch (tp)
 896         {
 897             case etabLJ6:
 898                 /* Dispersion */
 899                 Vcut = -rc6;
 900                 break;
 901             case etabLJ6Ewald:
 902                 Vcut = -rc6 * exp(-ewclj * ewclj * rc2)
 903                        * (1 + ewclj * ewclj * rc2 + gmx::power4(ewclj) * rc2 * rc2 / 2);
 904                 break;
 905             case etabLJ12:
 906                 /* Repulsion */
 907                 Vcut = rc12;
 908                 break;
 909             case etabCOUL: Vcut = 1.0 / rc; break;
 910             case etabEwald:
 911             case etabEwaldSwitch: Vcut = std::erfc(ewc * rc) / rc; break;
 912             case etabEwaldUser:
 913                 /* Only calculate minus the reciprocal space contribution */
 914                 Vcut = -std::erf(ewc * rc) / rc;
 915                 break;
 916             case etabRF:
 917             case etabRF_ZERO:
 918                 /* No need for preventing the usage of modifiers with RF */
 919                 Vcut = 0.0;
 920                 break;
 921             case etabEXPMIN: Vcut = exp(-rc); break;
 922             default:
 923                 gmx_fatal(FARGS,
 924                           "Cannot apply new potential-shift modifier to interaction type '%s' yet. "
 925                           "(%s,%d)",
 926                           tprops[tp].name, __FILE__, __LINE__);
 927         }
 928     }
 929
 930     for (i = 0; (i < td->nx); i++)
 931     {
 932         td->x[i] = i / td->tabscale;
 933     }
 934     for (i = td->nx0; (i < td->nx); i++)
 935     {
 936         r  = td->x[i];
 937         r2 = r * r;
 938         r6 = 1.0 / (r2 * r2 * r2);
 939         if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 940         {
 941             r12 = r6 * r6;
 942         }
 943         else
 944         {
 945             r12 = std::pow(r, -reppow);
 946         }
 947         Vtab = 0.0;
 948         Ftab = 0.0;
 949         if (bPotentialSwitch)
 950         {
 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
 954              * r1 and rc.
 955              * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
 956              */
 957             if (r <= r1)
 958             {
 959                 swi  = 1.0;
 960                 swi1 = 0.0;
 961             }
 962             else if (r >= rc)
 963             {
 964                 swi  = 0.0;
 965                 swi1 = 0.0;
 966             }
 967             else
 968             {
 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;
 973             }
 974         }
 975         else /* not really needed, but avoids compiler warnings... */
 976         {
 977             swi  = 1.0;
 978             swi1 = 0.0;
 979         }
 980
 981         rc6 = rc * rc * rc;
 982         rc6 = 1.0 / (rc6 * rc6);
 983
 984         switch (tp)
 985         {
 986             case etabLJ6:
 987                 /* Dispersion */
 988                 Vtab = -r6;
 989                 Ftab = 6.0 * Vtab / r;
 990                 break;
 991             case etabLJ6Switch:
 992             case etabLJ6Shift:
 993                 /* Dispersion */
 994                 if (r < rc)
 995                 {
 996                     Vtab = -r6;
 997                     Ftab = 6.0 * Vtab / r;
 998                     break;
 999                 }
1000                 break;
1001             case etabLJ12:
1002                 /* Repulsion */
1003                 Vtab = r12;
1004                 Ftab = reppow * Vtab / r;
1005                 break;
1006             case etabLJ12Switch:
1007             case etabLJ12Shift:
1008                 /* Repulsion */
1009                 if (r < rc)
1010                 {
1011                     Vtab = r12;
1012                     Ftab = reppow * Vtab / r;
1013                 }
1014                 break;
1015             case etabLJ6Encad:
1016                 if (r < rc)
1017                 {
1018                     Vtab = -(r6 - 6.0 * (rc - r) * rc6 / rc - rc6);
1019                     Ftab = -(6.0 * r6 / r - 6.0 * rc6 / rc);
1020                 }
1021                 else /* r>rc */
1022                 {
1023                     Vtab = 0;
1024                     Ftab = 0;
1025                 }
1026                 break;
1027             case etabLJ12Encad:
1028                 if (r < rc)
1029                 {
1030                     Vtab = -(r6 - 6.0 * (rc - r) * rc6 / rc - rc6);
1031                     Ftab = -(6.0 * r6 / r - 6.0 * rc6 / rc);
1032                 }
1033                 else /* r>rc */
1034                 {
1035                     Vtab = 0;
1036                     Ftab = 0;
1037                 }
1038                 break;
1039             case etabCOUL:
1040                 Vtab = 1.0 / r;
1041                 Ftab = 1.0 / r2;
1042                 break;
1043             case etabCOULSwitch:
1044             case etabShift:
1045                 if (r < rc)
1046                 {
1047                     Vtab = 1.0 / r;
1048                     Ftab = 1.0 / r2;
1049                 }
1050                 break;
1051             case etabEwald:
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;
1055                 break;
1056             case etabEwaldUser:
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;
1061                 break;
1062             case etabLJ6Ewald:
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;
1067                 break;
1068             case etabRF:
1069             case etabRF_ZERO:
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)
1073                 {
1074                     Vtab = 0;
1075                     Ftab = 0;
1076                 }
1077                 break;
1078             case etabEXPMIN:
1079                 expr = exp(-r);
1080                 Vtab = expr;
1081                 Ftab = expr;
1082                 break;
1083             case etabCOULEncad:
1084                 if (r < rc)
1085                 {
1086                     Vtab = 1.0 / r - (rc - r) / (rc * rc) - 1.0 / rc;
1087                     Ftab = 1.0 / r2 - 1.0 / (rc * rc);
1088                 }
1089                 else /* r>rc */
1090                 {
1091                     Vtab = 0;
1092                     Ftab = 0;
1093                 }
1094                 break;
1095             default:
1096                 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)", tp, __FILE__, __LINE__);
1097         }
1098         if (bForceSwitch)
1099         {
1100             /* Normal coulomb with cut-off correction for potential */
1101             if (r < rc)
1102             {
1103                 Vtab -= C;
1104                 /* If in Shifting range add something to it */
1105                 if (r > r1)
1106                 {
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;
1111                 }
1112             }
1113             else
1114             {
1115                 /* Make sure interactions are zero outside cutoff with modifiers */
1116                 Vtab = 0;
1117                 Ftab = 0;
1118             }
1119         }
1120         if (bPotentialShift)
1121         {
1122             if (r < rc)
1123             {
1124                 Vtab -= Vcut;
1125             }
1126             else
1127             {
1128                 /* Make sure interactions are zero outside cutoff with modifiers */
1129                 Vtab = 0;
1130                 Ftab = 0;
1131             }
1132         }
1133
1134         if (ETAB_USER(tp))
1135         {
1136             Vtab += td->v[i];
1137             Ftab += td->f[i];
1138         }
1139
1140         if (bPotentialSwitch)
1141         {
1142             if (r >= rc)
1143             {
1144                 /* Make sure interactions are zero outside cutoff with modifiers */
1145                 Vtab = 0;
1146                 Ftab = 0;
1147             }
1148             else if (r > r1)
1149             {
1150                 Ftab = Ftab * swi - Vtab * swi1;
1151                 Vtab = Vtab * swi;
1152             }
1153         }
1154         /* Convert to single precision when we store to mem */
1155         td->v[i] = Vtab;
1156         td->f[i] = Ftab;
1157     }
1158
1159     /* Continue the table linearly from nx0 to 0.
1160      * These values are only required for energy minimization with overlap or TPI.
1161      */
1162     for (i = td->nx0 - 1; i >= 0; i--)
1163     {
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];
1166     }
1167 }
1168
1169 static void set_table_type(int tabsel[], const interaction_const_t* ic, gmx_bool b14only)
1170 {
1171     int eltype, vdwtype;
1172
1173     /* Set the different table indices.
1174      * Coulomb first.
1175      */
1176
1177
1178     if (b14only)
1179     {
1180         switch (ic->eeltype)
1181         {
1182             case eelUSER:
1183             case eelPMEUSER:
1184             case eelPMEUSERSWITCH: eltype = eelUSER; break;
1185             default: eltype = eelCUT;
1186         }
1187     }
1188     else
1189     {
1190         eltype = ic->eeltype;
1191     }
1192
1193     switch (eltype)
1194     {
1195         case eelCUT: tabsel[etiCOUL] = etabCOUL; break;
1196         case eelPOISSON: tabsel[etiCOUL] = etabShift; break;
1197         case eelSHIFT:
1198             if (ic->rcoulomb > ic->rcoulomb_switch)
1199             {
1200                 tabsel[etiCOUL] = etabShift;
1201             }
1202             else
1203             {
1204                 tabsel[etiCOUL] = etabCOUL;
1205             }
1206             break;
1207         case eelEWALD:
1208         case eelPME:
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;
1213         case eelRF:
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);
1219     }
1220
1221     /* Van der Waals time */
1222     if (ic->useBuckingham && !b14only)
1223     {
1224         tabsel[etiLJ6]  = etabLJ6;
1225         tabsel[etiLJ12] = etabEXPMIN;
1226     }
1227     else
1228     {
1229         if (b14only && ic->vdwtype != evdwUSER)
1230         {
1231             vdwtype = evdwCUT;
1232         }
1233         else
1234         {
1235             vdwtype = ic->vdwtype;
1236         }
1237
1238         switch (vdwtype)
1239         {
1240             case evdwSWITCH:
1241                 tabsel[etiLJ6]  = etabLJ6Switch;
1242                 tabsel[etiLJ12] = etabLJ12Switch;
1243                 break;
1244             case evdwSHIFT:
1245                 tabsel[etiLJ6]  = etabLJ6Shift;
1246                 tabsel[etiLJ12] = etabLJ12Shift;
1247                 break;
1248             case evdwUSER:
1249                 tabsel[etiLJ6]  = etabUSER;
1250                 tabsel[etiLJ12] = etabUSER;
1251                 break;
1252             case evdwCUT:
1253                 tabsel[etiLJ6]  = etabLJ6;
1254                 tabsel[etiLJ12] = etabLJ12;
1255                 break;
1256             case evdwENCADSHIFT:
1257                 tabsel[etiLJ6]  = etabLJ6Encad;
1258                 tabsel[etiLJ12] = etabLJ12Encad;
1259                 break;
1260             case evdwPME:
1261                 tabsel[etiLJ6]  = etabLJ6Ewald;
1262                 tabsel[etiLJ12] = etabLJ12;
1263                 break;
1264             default:
1265                 gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype, __FILE__, __LINE__);
1266         }
1267
1268         if (!b14only && ic->vdw_modifier != eintmodNONE)
1269         {
1270             if (ic->vdw_modifier != eintmodPOTSHIFT && ic->vdwtype != evdwCUT)
1271             {
1272                 gmx_incons(
1273                         "Potential modifiers other than potential-shift are only implemented for "
1274                         "LJ cut-off");
1275             }
1276
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.
1279              */
1280             if (ic->vdwtype == evdwCUT)
1281             {
1282                 switch (ic->vdw_modifier)
1283                 {
1284                     case eintmodNONE:
1285                     case eintmodPOTSHIFT:
1286                     case eintmodEXACTCUTOFF:
1287                         /* No modification */
1288                         break;
1289                     case eintmodPOTSWITCH:
1290                         tabsel[etiLJ6]  = etabLJ6Switch;
1291                         tabsel[etiLJ12] = etabLJ12Switch;
1292                         break;
1293                     case eintmodFORCESWITCH:
1294                         tabsel[etiLJ6]  = etabLJ6Shift;
1295                         tabsel[etiLJ12] = etabLJ12Shift;
1296                         break;
1297                     default: gmx_incons("Unsupported vdw_modifier");
1298                 }
1299             }
1300         }
1301     }
1302 }
1303
1304 t_forcetable* make_tables(FILE* out, const interaction_const_t* ic, const char* fn, real rtab, int flags)
1305 {
1306     t_tabledata* td;
1307     gmx_bool     b14only, useUserTable;
1308     int          nx0, tabsel[etiNR];
1309     real         scalefactor;
1310
1311     t_forcetable* table = new t_forcetable(GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
1312                                            GMX_TABLE_FORMAT_CUBICSPLINE_YFGH);
1313
1314     b14only = ((flags & GMX_MAKETABLES_14ONLY) != 0);
1315
1316     if (flags & GMX_MAKETABLES_FORCEUSER)
1317     {
1318         tabsel[etiCOUL] = etabUSER;
1319         tabsel[etiLJ6]  = etabUSER;
1320         tabsel[etiLJ12] = etabUSER;
1321     }
1322     else
1323     {
1324         set_table_type(tabsel, ic, b14only);
1325     }
1326     snew(td, etiNR);
1327     table->r     = rtab;
1328     table->scale = 0;
1329     table->n     = 0;
1330
1331     table->formatsize    = 4;
1332     table->ninteractions = etiNR;
1333     table->stride        = table->formatsize * table->ninteractions;
1334
1335     /* Check whether we have to read or generate */
1336     useUserTable = FALSE;
1337     for (unsigned int i = 0; (i < etiNR); i++)
1338     {
1339         if (ETAB_USER(tabsel[i]))
1340         {
1341             useUserTable = TRUE;
1342         }
1343     }
1344     if (useUserTable)
1345     {
1346         read_tables(out, fn, etiNR, 0, td);
1347         if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1348         {
1349             table->n = td[0].nx;
1350         }
1351         else
1352         {
1353             if (td[0].x[td[0].nx - 1] < rtab)
1354             {
1355                 gmx_fatal(FARGS,
1356                           "Tables in file %s not long enough for cut-off:\n"
1357                           "\tshould be at least %f nm\n",
1358                           fn, rtab);
1359             }
1360             table->n = gmx::roundToInt(rtab * td[0].tabscale);
1361         }
1362         table->scale = td[0].tabscale;
1363         nx0          = td[0].nx0;
1364     }
1365     else
1366     {
1367         // No tables are read
1368 #if GMX_DOUBLE
1369         table->scale = 2000.0;
1370 #else
1371         table->scale = 500.0;
1372 #endif
1373         table->n = static_cast<int>(rtab * table->scale);
1374         nx0      = 10;
1375     }
1376
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.
1380      */
1381     snew_aligned(table->data, table->stride * (table->n + 1) * sizeof(real), 32);
1382
1383     for (int k = 0; (k < etiNR); k++)
1384     {
1385         /* Now fill data for tables that have not been read
1386          * or add the Ewald long-range correction for Ewald user tables.
1387          */
1388         if (tabsel[k] != etabUSER)
1389         {
1390             real scale = table->scale;
1391             if (ic->useBuckingham && (ic->buckinghamBMax != 0) && tabsel[k] == etabEXPMIN)
1392             {
1393                 scale /= ic->buckinghamBMax;
1394             }
1395             init_table(table->n, nx0, scale, &(td[k]), !useUserTable);
1396
1397             fill_table(&(td[k]), tabsel[k], ic, b14only);
1398             if (out)
1399             {
1400                 fprintf(out,
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);
1404             }
1405         }
1406
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.
1412          */
1413         if (k == etiLJ6)
1414         {
1415             scalefactor = 1.0 / 6.0;
1416         }
1417         else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1418         {
1419             scalefactor = 1.0 / 12.0;
1420         }
1421         else
1422         {
1423             scalefactor = 1.0;
1424         }
1425
1426         copy2table(table->n, k * table->formatsize, table->stride, td[k].x, td[k].v, td[k].f,
1427                    scalefactor, table->data);
1428
1429         done_tabledata(&(td[k]));
1430     }
1431     sfree(td);
1432
1433     return table;
1434 }
1435
1436 bondedtable_t make_bonded_table(FILE* fplog, const char* fn, int angle)
1437 {
1438     t_tabledata   td;
1439     int           i;
1440     bondedtable_t tab;
1441     int           stride = 4;
1442
1443     read_tables(fplog, fn, 1, angle, &td);
1444     if (angle > 0)
1445     {
1446         /* Convert the table from degrees to radians */
1447         for (i = 0; i < td.nx; i++)
1448         {
1449             td.x[i] *= DEG2RAD;
1450             td.f[i] *= RAD2DEG;
1451         }
1452         td.tabscale *= RAD2DEG;
1453     }
1454     tab.n     = td.nx;
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);
1459
1460     return tab;
1461 }
1462
1463 std::unique_ptr<t_forcetable>
1464 makeDispersionCorrectionTable(FILE* fp, const interaction_const_t* ic, real rtab, const char* tabfn)
1465 {
1466     GMX_RELEASE_ASSERT(ic->vdwtype != evdwUSER || tabfn,
1467                        "With VdW user tables we need a table file name");
1468
1469     if (tabfn == nullptr)
1470     {
1471         return std::unique_ptr<t_forcetable>(nullptr);
1472     }
1473
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.
1478      */
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);
1490
1491     for (int i = 0; i <= fullTable->n; i++)
1492     {
1493         for (int j = 0; j < 8; j++)
1494         {
1495             dispersionCorrectionTable->data[8 * i + j] = fullTable->data[12 * i + 4 + j];
1496         }
1497     }
1498     delete fullTable;
1499
1500     return dispersionCorrectionTable;
1501 }
1502
1503 t_forcetable::t_forcetable(enum gmx_table_interaction interaction, enum gmx_table_format format) :
1504     interaction(interaction),
1505     format(format),
1506     r(0),
1507     n(0),
1508     scale(0),
1509     data(nullptr),
1510     formatsize(0),
1511     ninteractions(0),
1512     stride(0)
1513 {
1514 }
1515
1516 t_forcetable::~t_forcetable()
1517 {
1518     sfree_aligned(data);
1519 }