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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "forcetable.h"
  41
  42 #include <cmath>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/fileio/xvgr.h"
  47 #include "gromacs/math/functions.h"
  48 #include "gromacs/math/multidimarray.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/math/utilities.h"
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdspan/extensions.h"
  53 #include "gromacs/mdtypes/fcdata.h"
  54 #include "gromacs/mdtypes/interaction_const.h"
  55 #include "gromacs/mdtypes/md_enums.h"
  56 #include "gromacs/mdtypes/nblist.h"
  57 #include "gromacs/utility/arrayref.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/futil.h"
  61 #include "gromacs/utility/smalloc.h"
  62
  63 /* All the possible (implemented) table functions */
  64 enum
  65 {
  66     etabLJ6,
  67     etabLJ12,
  68     etabLJ6Shift,
  69     etabLJ12Shift,
  70     etabShift,
  71     etabRF,
  72     etabRF_ZERO,
  73     etabCOUL,
  74     etabEwald,
  75     etabEwaldSwitch,
  76     etabEwaldUser,
  77     etabEwaldUserSwitch,
  78     etabLJ6Ewald,
  79     etabLJ6Switch,
  80     etabLJ12Switch,
  81     etabCOULSwitch,
  82     etabEXPMIN,
  83     etabUSER,
  84     etabNR
  85 };
  86
  87 /** Evaluates to true if the table type contains user data. */
  88 #define ETAB_USER(e) ((e) == etabUSER || (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
  89
  90 typedef struct
  91 {
  92     const char* name;
  93     gmx_bool    bCoulomb;
  94 } t_tab_props;
  95
  96 /* This structure holds name and a flag that tells whether
  97    this is a Coulomb type funtion */
  98 static const t_tab_props tprops[etabNR] = {
  99     { "LJ6", FALSE },         { "LJ12", FALSE },      { "LJ6Shift", FALSE },
 100     { "LJ12Shift", FALSE },   { "Shift", TRUE },      { "RF", TRUE },
 101     { "RF-zero", TRUE },      { "COUL", TRUE },       { "Ewald", TRUE },
 102     { "Ewald-Switch", TRUE }, { "Ewald-User", TRUE }, { "Ewald-User-Switch", TRUE },
 103     { "LJ6Ewald", FALSE },    { "LJ6Switch", FALSE }, { "LJ12Switch", FALSE },
 104     { "COULSwitch", TRUE },   { "EXPMIN", FALSE },    { "USER", FALSE },
 105 };
 106
 107 typedef struct
 108 {
 109     int     nx, nx0;
 110     double  tabscale;
 111     double *x, *v, *f;
 112 } t_tabledata;
 113
 114 double v_q_ewald_lr(double beta, double r)
 115 {
 116     if (r == 0)
 117     {
 118         return beta * 2 / sqrt(M_PI);
 119     }
 120     else
 121     {
 122         return std::erf(beta * r) / r;
 123     }
 124 }
 125
 126 double v_lj_ewald_lr(double beta, double r)
 127 {
 128     double br, br2, br4, r6, factor;
 129     if (r == 0)
 130     {
 131         return gmx::power6(beta) / 6;
 132     }
 133     else
 134     {
 135         br     = beta * r;
 136         br2    = br * br;
 137         br4    = br2 * br2;
 138         r6     = gmx::power6(r);
 139         factor = (1.0 - std::exp(-br2) * (1 + br2 + 0.5 * br4)) / r6;
 140         return factor;
 141     }
 142 }
 143
 144 EwaldCorrectionTables generateEwaldCorrectionTables(const int    numPoints,
 145                                                     const double tableScaling,
 146                                                     const real   beta,
 147                                                     real_space_grid_contribution_computer v_lr)
 148 {
 149     real     tab_max;
 150     int      i, i_inrange;
 151     double   dc, dc_new;
 152     gmx_bool bOutOfRange;
 153     double   v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
 154     double   x_r0;
 155
 156     /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument
 157      * depending on wether we should create electrostatic or Lennard-Jones Ewald tables.
 158      */
 159
 160     if (numPoints < 2)
 161     {
 162         gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
 163     }
 164
 165     const double dx = 1 / tableScaling;
 166
 167     EwaldCorrectionTables tables;
 168     tables.scale = tableScaling;
 169     tables.tableF.resize(numPoints);
 170     tables.tableV.resize(numPoints);
 171     tables.tableFDV0.resize(numPoints * 4);
 172     gmx::ArrayRef<real> table_f    = tables.tableF;
 173     gmx::ArrayRef<real> table_v    = tables.tableV;
 174     gmx::ArrayRef<real> table_fdv0 = tables.tableFDV0;
 175
 176     /* We need some margin to be able to divide table values by r
 177      * in the kernel and also to do the integration arithmetics
 178      * without going out of range. Furthemore, we divide by dx below.
 179      */
 180     tab_max = GMX_REAL_MAX * 0.0001;
 181
 182     /* This function produces a table with:
 183      * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
 184      * maximum force error:  V'''/(6*4)*dx^2
 185      * The rms force error is the max error times 1/sqrt(5)=0.45.
 186      */
 187
 188     bOutOfRange = FALSE;
 189     i_inrange   = numPoints;
 190     v_inrange   = 0;
 191     dc          = 0;
 192     for (i = numPoints - 1; i >= 0; i--)
 193     {
 194         x_r0 = i * dx;
 195
 196         v_r0 = (*v_lr)(beta, x_r0);
 197
 198         if (!bOutOfRange)
 199         {
 200             i_inrange = i;
 201             v_inrange = v_r0;
 202
 203             vi = v_r0;
 204         }
 205         else
 206         {
 207             /* Linear continuation for the last point in range */
 208             vi = v_inrange - dc * (i - i_inrange) * dx;
 209         }
 210
 211         table_v[i] = vi;
 212
 213         if (i == 0)
 214         {
 215             continue;
 216         }
 217
 218         /* Get the potential at table point i-1 */
 219         v_r1 = (*v_lr)(beta, (i - 1) * dx);
 220
 221         if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
 222         {
 223             bOutOfRange = TRUE;
 224         }
 225
 226         if (!bOutOfRange)
 227         {
 228             /* Calculate the average second derivative times dx over interval i-1 to i.
 229              * Using the function values at the end points and in the middle.
 230              */
 231             a2dx = (v_r0 + v_r1 - 2 * (*v_lr)(beta, x_r0 - 0.5 * dx)) / (0.25 * dx);
 232             /* Set the derivative of the spline to match the difference in potential
 233              * over the interval plus the average effect of the quadratic term.
 234              * This is the essential step for minimizing the error in the force.
 235              */
 236             dc = (v_r0 - v_r1) / dx + 0.5 * a2dx;
 237         }
 238
 239         if (i == numPoints - 1)
 240         {
 241             /* Fill the table with the force, minus the derivative of the spline */
 242             table_f[i] = -dc;
 243         }
 244         else
 245         {
 246             /* tab[i] will contain the average of the splines over the two intervals */
 247             table_f[i] += -0.5 * dc;
 248         }
 249
 250         if (!bOutOfRange)
 251         {
 252             /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
 253              * matching the potential at the two end points
 254              * and the derivative dc at the end point xr.
 255              */
 256             a0   = v_r0;
 257             a1   = dc;
 258             a2dx = (a1 * dx + v_r1 - a0) * 2 / dx;
 259
 260             /* Set dc to the derivative at the next point */
 261             dc_new = a1 - a2dx;
 262
 263             if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
 264             {
 265                 bOutOfRange = TRUE;
 266             }
 267             else
 268             {
 269                 dc = dc_new;
 270             }
 271         }
 272
 273         table_f[(i - 1)] = -0.5 * dc;
 274     }
 275     /* Currently the last value only contains half the force: double it */
 276     table_f[0] *= 2;
 277
 278     if (!table_fdv0.empty())
 279     {
 280         /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
 281          * init_ewald_f_table().
 282          */
 283         for (i = 0; i < numPoints - 1; i++)
 284         {
 285             table_fdv0[4 * i]     = table_f[i];
 286             table_fdv0[4 * i + 1] = table_f[i + 1] - table_f[i];
 287             table_fdv0[4 * i + 2] = table_v[i];
 288             table_fdv0[4 * i + 3] = 0.0;
 289         }
 290         const int lastPoint           = numPoints - 1;
 291         table_fdv0[4 * lastPoint]     = table_f[lastPoint];
 292         table_fdv0[4 * lastPoint + 1] = -table_f[lastPoint];
 293         table_fdv0[4 * lastPoint + 2] = table_v[lastPoint];
 294         table_fdv0[4 * lastPoint + 3] = 0.0;
 295     }
 296
 297     return tables;
 298 }
 299
 300 /* Returns the spacing for a function using the maximum of
 301  * the third derivative, x_scale (unit 1/length)
 302  * and function tolerance.
 303  */
 304 static double spline3_table_scale(double third_deriv_max, double x_scale, double func_tol)
 305 {
 306     double deriv_tol;
 307     double sc_deriv, sc_func;
 308
 309     /* Force tolerance: single precision accuracy */
 310     deriv_tol = GMX_FLOAT_EPS;
 311     sc_deriv  = sqrt(third_deriv_max / (6 * 4 * deriv_tol * x_scale)) * x_scale;
 312
 313     /* Don't try to be more accurate on energy than the precision */
 314     func_tol = std::max(func_tol, static_cast<double>(GMX_REAL_EPS));
 315     sc_func  = std::cbrt(third_deriv_max / (6 * 12 * std::sqrt(3.0) * func_tol)) * x_scale;
 316
 317     return std::max(sc_deriv, sc_func);
 318 }
 319
 320 /* The scale (1/spacing) for third order spline interpolation
 321  * of the Ewald mesh contribution which needs to be subtracted
 322  * from the non-bonded interactions.
 323  * Since there is currently only one spacing for Coulomb and LJ,
 324  * the finest spacing is used if both Ewald types are present.
 325  *
 326  * Note that we could also implement completely separate tables
 327  * for Coulomb and LJ Ewald, each with their own spacing.
 328  * The current setup with the same spacing can provide slightly
 329  * faster kernels with both Coulomb and LJ Ewald, especially
 330  * when interleaving both tables (currently not implemented).
 331  */
 332 real ewald_spline3_table_scale(const interaction_const_t& ic,
 333                                const bool                 generateCoulombTables,
 334                                const bool                 generateVdwTables)
 335 {
 336     GMX_RELEASE_ASSERT(!generateCoulombTables || EEL_PME_EWALD(ic.eeltype),
 337                        "Can only use tables with Ewald");
 338     GMX_RELEASE_ASSERT(!generateVdwTables || EVDW_PME(ic.vdwtype), "Can only use tables with Ewald");
 339
 340     real sc = 0;
 341
 342     if (generateCoulombTables)
 343     {
 344         GMX_RELEASE_ASSERT(ic.ewaldcoeff_q > 0, "The Ewald coefficient shoule be positive");
 345
 346         double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
 347         double etol;
 348         real   sc_q;
 349
 350         /* Energy tolerance: 0.1 times the cut-off jump */
 351         etol = 0.1 * std::erfc(ic.ewaldcoeff_q * ic.rcoulomb);
 352
 353         sc_q = spline3_table_scale(erf_x_d3, ic.ewaldcoeff_q, etol);
 354
 355         if (debug)
 356         {
 357             fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f nm\n", 1 / sc_q);
 358         }
 359
 360         sc = std::max(sc, sc_q);
 361     }
 362
 363     if (generateVdwTables)
 364     {
 365         GMX_RELEASE_ASSERT(ic.ewaldcoeff_lj > 0, "The Ewald coefficient shoule be positive");
 366
 367         double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
 368         double xrc2, etol;
 369         real   sc_lj;
 370
 371         /* Energy tolerance: 0.1 times the cut-off jump */
 372         xrc2 = gmx::square(ic.ewaldcoeff_lj * ic.rvdw);
 373         etol = 0.1 * std::exp(-xrc2) * (1 + xrc2 + xrc2 * xrc2 / 2.0);
 374
 375         sc_lj = spline3_table_scale(func_d3, ic.ewaldcoeff_lj, etol);
 376
 377         if (debug)
 378         {
 379             fprintf(debug, "Ewald LJ quadratic spline table spacing: %f nm\n", 1 / sc_lj);
 380         }
 381
 382         sc = std::max(sc, sc_lj);
 383     }
 384
 385     return sc;
 386 }
 387
 388 static void copy2table(int          n,
 389                        int          offset,
 390                        int          stride,
 391                        const double x[],
 392                        const double Vtab[],
 393                        const double Ftab[],
 394                        real         scalefactor,
 395                        real         dest[])
 396 {
 397     /* Use double prec. for the intermediary variables
 398      * and temporary x/vtab/vtab2 data to avoid unnecessary
 399      * loss of precision.
 400      */
 401     int    i, nn0;
 402     double F, G, H, h;
 403
 404     h = 0;
 405     for (i = 0; (i < n); i++)
 406     {
 407         if (i < n - 1)
 408         {
 409             h = x[i + 1] - x[i];
 410             F = -Ftab[i] * h;
 411             G = 3 * (Vtab[i + 1] - Vtab[i]) + (Ftab[i + 1] + 2 * Ftab[i]) * h;
 412             H = -2 * (Vtab[i + 1] - Vtab[i]) - (Ftab[i + 1] + Ftab[i]) * h;
 413         }
 414         else
 415         {
 416             /* Fill the last entry with a linear potential,
 417              * this is mainly for rounding issues with angle and dihedral potentials.
 418              */
 419             F = -Ftab[i] * h;
 420             G = 0;
 421             H = 0;
 422         }
 423         nn0           = offset + i * stride;
 424         dest[nn0]     = scalefactor * Vtab[i];
 425         dest[nn0 + 1] = scalefactor * F;
 426         dest[nn0 + 2] = scalefactor * G;
 427         dest[nn0 + 3] = scalefactor * H;
 428     }
 429 }
 430
 431 static void init_table(int n, int nx0, double tabscale, t_tabledata* td, gmx_bool bAlloc)
 432 {
 433     td->nx       = n;
 434     td->nx0      = nx0;
 435     td->tabscale = tabscale;
 436     if (bAlloc)
 437     {
 438         snew(td->x, td->nx);
 439         snew(td->v, td->nx);
 440         snew(td->f, td->nx);
 441     }
 442 }
 443
 444 static void spline_forces(int nx, double h, const double v[], gmx_bool bS3, gmx_bool bE3, double f[])
 445 {
 446     int    start, end, i;
 447     double v3, b_s, b_e, b;
 448     double beta, *gamma;
 449
 450     /* Formulas can be found in:
 451      * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
 452      */
 453
 454     if (nx < 4 && (bS3 || bE3))
 455     {
 456         gmx_fatal(FARGS,
 457                   "Can not generate splines with third derivative boundary conditions with less "
 458                   "than 4 (%d) points",
 459                   nx);
 460     }
 461
 462     /* To make life easy we initially set the spacing to 1
 463      * and correct for this at the end.
 464      */
 465     if (bS3)
 466     {
 467         /* Fit V''' at the start */
 468         v3 = v[3] - 3 * v[2] + 3 * v[1] - v[0];
 469         if (debug)
 470         {
 471             fprintf(debug, "The left third derivative is %g\n", v3 / (h * h * h));
 472         }
 473         b_s   = 2 * (v[1] - v[0]) + v3 / 6;
 474         start = 0;
 475
 476         if (FALSE)
 477         {
 478             /* Fit V'' at the start */
 479             real v2;
 480
 481             v2 = -v[3] + 4 * v[2] - 5 * v[1] + 2 * v[0];
 482             /* v2  = v[2] - 2*v[1] + v[0]; */
 483             if (debug)
 484             {
 485                 fprintf(debug, "The left second derivative is %g\n", v2 / (h * h));
 486             }
 487             b_s   = 3 * (v[1] - v[0]) - v2 / 2;
 488             start = 0;
 489         }
 490     }
 491     else
 492     {
 493         b_s   = 3 * (v[2] - v[0]) + f[0] * h;
 494         start = 1;
 495     }
 496     if (bE3)
 497     {
 498         /* Fit V''' at the end */
 499         v3 = v[nx - 1] - 3 * v[nx - 2] + 3 * v[nx - 3] - v[nx - 4];
 500         if (debug)
 501         {
 502             fprintf(debug, "The right third derivative is %g\n", v3 / (h * h * h));
 503         }
 504         b_e = 2 * (v[nx - 1] - v[nx - 2]) + v3 / 6;
 505         end = nx;
 506     }
 507     else
 508     {
 509         /* V'=0 at the end */
 510         b_e = 3 * (v[nx - 1] - v[nx - 3]) + f[nx - 1] * h;
 511         end = nx - 1;
 512     }
 513
 514     snew(gamma, nx);
 515     beta = (bS3 ? 1 : 4);
 516
 517     /* For V'' fitting */
 518     /* beta = (bS3 ? 2 : 4); */
 519
 520     f[start] = b_s / beta;
 521     for (i = start + 1; i < end; i++)
 522     {
 523         gamma[i] = 1 / beta;
 524         beta     = 4 - gamma[i];
 525         b        = 3 * (v[i + 1] - v[i - 1]);
 526         f[i]     = (b - f[i - 1]) / beta;
 527     }
 528     gamma[end - 1] = 1 / beta;
 529     beta           = (bE3 ? 1 : 4) - gamma[end - 1];
 530     f[end - 1]     = (b_e - f[end - 2]) / beta;
 531
 532     for (i = end - 2; i >= start; i--)
 533     {
 534         f[i] -= gamma[i + 1] * f[i + 1];
 535     }
 536     sfree(gamma);
 537
 538     /* Correct for the minus sign and the spacing */
 539     for (i = start; i < end; i++)
 540     {
 541         f[i] = -f[i] / h;
 542     }
 543 }
 544
 545 static void set_forces(FILE* fp, int angle, int nx, double h, double v[], double f[], int table)
 546 {
 547     int start, end;
 548
 549     if (angle == 2)
 550     {
 551         gmx_fatal(FARGS, "Force generation for dihedral tables is not (yet) implemented");
 552     }
 553
 554     start = 0;
 555     while (v[start] == 0)
 556     {
 557         start++;
 558     }
 559
 560     end = nx;
 561     while (v[end - 1] == 0)
 562     {
 563         end--;
 564     }
 565     if (end > nx - 2)
 566     {
 567         end = nx;
 568     }
 569     else
 570     {
 571         end++;
 572     }
 573
 574     if (fp)
 575     {
 576         fprintf(fp,
 577                 "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
 578                 table + 1,
 579                 start * h,
 580                 end == nx ? "V'''" : "V'=0",
 581                 (end - 1) * h);
 582     }
 583     spline_forces(end - start, h, v + start, TRUE, end == nx, f + start);
 584 }
 585
 586 static void read_tables(FILE* fp, const char* filename, int ntab, int angle, t_tabledata td[])
 587 {
 588     char     buf[STRLEN];
 589     double   start, end, dx0, dx1, ssd, vm, vp, f, numf;
 590     int      k, i, nx0 = 0, nny, ns;
 591     gmx_bool bAllZero, bZeroV, bZeroF;
 592     double   tabscale;
 593
 594     nny               = 2 * ntab + 1;
 595     std::string libfn = gmx::findLibraryFile(filename);
 596     gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgData    = readXvgData(libfn);
 597     int                                                            numColumns = xvgData.extent(0);
 598     if (numColumns != nny)
 599     {
 600         gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d", libfn.c_str(), numColumns, nny);
 601     }
 602     int numRows = xvgData.extent(1);
 603
 604     const auto& yy = xvgData.asView();
 605     if (angle == 0)
 606     {
 607         if (yy[0][0] != 0.0)
 608         {
 609             gmx_fatal(FARGS,
 610                       "The first distance in file %s is %f nm instead of %f nm",
 611                       libfn.c_str(),
 612                       yy[0][0],
 613                       0.0);
 614         }
 615     }
 616     else
 617     {
 618         if (angle == 1)
 619         {
 620             start = 0.0;
 621         }
 622         else
 623         {
 624             start = -180.0;
 625         }
 626         end = 180.0;
 627         if (yy[0][0] != start || yy[0][numRows - 1] != end)
 628         {
 629             gmx_fatal(FARGS,
 630                       "The angles in file %s should go from %f to %f instead of %f to %f\n",
 631                       libfn.c_str(),
 632                       start,
 633                       end,
 634                       yy[0][0],
 635                       yy[0][numRows - 1]);
 636         }
 637     }
 638
 639     tabscale = (numRows - 1) / (yy[0][numRows - 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(), numRows);
 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 < numRows); 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,
 667                               yy[0][i - 2],
 668                               yy[0][i - 1],
 669                               yy[0][i]);
 670                 }
 671             }
 672             if (yy[1 + k * 2][i] != 0)
 673             {
 674                 bZeroV = FALSE;
 675                 if (bAllZero)
 676                 {
 677                     bAllZero = FALSE;
 678                     nx0      = i;
 679                 }
 680                 if (yy[1 + k * 2][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2][i] < -0.01 * GMX_REAL_MAX)
 681                 {
 682                     gmx_fatal(FARGS, "Out of range potential value %g in file '%s'", yy[1 + k * 2][i], filename);
 683                 }
 684             }
 685             if (yy[1 + k * 2 + 1][i] != 0)
 686             {
 687                 bZeroF = FALSE;
 688                 if (bAllZero)
 689                 {
 690                     bAllZero = FALSE;
 691                     nx0      = i;
 692                 }
 693                 if (yy[1 + k * 2 + 1][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2 + 1][i] < -0.01 * GMX_REAL_MAX)
 694                 {
 695                     gmx_fatal(FARGS, "Out of range force value %g in file '%s'", yy[1 + k * 2 + 1][i], filename);
 696                 }
 697             }
 698         }
 699
 700         if (!bZeroV && bZeroF)
 701         {
 702             set_forces(fp, angle, numRows, 1 / tabscale, yy[1 + k * 2].data(), yy[1 + k * 2 + 1].data(), k);
 703         }
 704         else
 705         {
 706             /* Check if the second column is close to minus the numerical
 707              * derivative of the first column.
 708              */
 709             ssd = 0;
 710             ns  = 0;
 711             for (i = 1; (i < numRows - 1); i++)
 712             {
 713                 vm = yy[1 + 2 * k][i - 1];
 714                 vp = yy[1 + 2 * k][i + 1];
 715                 f  = yy[1 + 2 * k + 1][i];
 716                 if (vm != 0 && vp != 0 && f != 0)
 717                 {
 718                     /* Take the centered difference */
 719                     numf = -(vp - vm) * 0.5 * tabscale;
 720                     if (f + numf != 0)
 721                     {
 722                         ssd += fabs(2 * (f - numf) / (f + numf));
 723                     }
 724                     ns++;
 725                 }
 726             }
 727             if (ns > 0)
 728             {
 729                 ssd /= ns;
 730                 sprintf(buf,
 731                         "For the %d non-zero entries for table %d in %s the forces deviate on "
 732                         "average %" PRId64
 733                         "%% from minus the numerical derivative of the potential\n",
 734                         ns,
 735                         k,
 736                         libfn.c_str(),
 737                         gmx::roundToInt64(100 * ssd));
 738                 if (debug)
 739                 {
 740                     fprintf(debug, "%s", buf);
 741                 }
 742                 if (ssd > 0.2)
 743                 {
 744                     if (fp)
 745                     {
 746                         fprintf(fp, "\nWARNING: %s\n", buf);
 747                     }
 748                     fprintf(stderr, "\nWARNING: %s\n", buf);
 749                 }
 750             }
 751         }
 752     }
 753     if (bAllZero && fp)
 754     {
 755         fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn.c_str());
 756     }
 757
 758     for (k = 0; (k < ntab); k++)
 759     {
 760         init_table(numRows, nx0, tabscale, &(td[k]), TRUE);
 761         for (i = 0; (i < numRows); i++)
 762         {
 763             td[k].x[i] = yy[0][i];
 764             td[k].v[i] = yy[2 * k + 1][i];
 765             td[k].f[i] = yy[2 * k + 2][i];
 766         }
 767     }
 768 }
 769
 770 static void done_tabledata(t_tabledata* td)
 771 {
 772     if (!td)
 773     {
 774         return;
 775     }
 776
 777     sfree(td->x);
 778     sfree(td->v);
 779     sfree(td->f);
 780 }
 781
 782 static void fill_table(t_tabledata* td, int tp, const interaction_const_t* ic, gmx_bool b14only)
 783 {
 784     /* Fill the table according to the formulas in the manual.
 785      * In principle, we only need the potential and the second
 786      * derivative, but then we would have to do lots of calculations
 787      * in the inner loop. By precalculating some terms (see manual)
 788      * we get better eventual performance, despite a larger table.
 789      *
 790      * Since some of these higher-order terms are very small,
 791      * we always use double precision to calculate them here, in order
 792      * to avoid unnecessary loss of precision.
 793      */
 794     int    i;
 795     double reppow, p;
 796     double r1, rc, r12, r13;
 797     double r, r2, r6, rc2;
 798     double expr, Vtab, Ftab;
 799     /* Parameters for David's function */
 800     double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
 801     /* Parameters for the switching function */
 802     double ksw, swi, swi1;
 803     /* Temporary parameters */
 804     gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
 805     double   ewc   = ic->ewaldcoeff_q;
 806     double   ewclj = ic->ewaldcoeff_lj;
 807     double   Vcut  = 0;
 808
 809     if (b14only)
 810     {
 811         bPotentialSwitch = FALSE;
 812         bForceSwitch     = FALSE;
 813         bPotentialShift  = FALSE;
 814     }
 815     else
 816     {
 817         bPotentialSwitch =
 818                 ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) || (tp == etabCOULSwitch)
 819                  || (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch)
 820                  || (tprops[tp].bCoulomb && (ic->coulomb_modifier == InteractionModifiers::PotSwitch))
 821                  || (!tprops[tp].bCoulomb && (ic->vdw_modifier == InteractionModifiers::PotSwitch)));
 822         bForceSwitch =
 823                 ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) || (tp == etabShift)
 824                  || (tprops[tp].bCoulomb && (ic->coulomb_modifier == InteractionModifiers::ForceSwitch))
 825                  || (!tprops[tp].bCoulomb && (ic->vdw_modifier == InteractionModifiers::ForceSwitch)));
 826         bPotentialShift =
 827                 ((tprops[tp].bCoulomb && (ic->coulomb_modifier == InteractionModifiers::PotShift))
 828                  || (!tprops[tp].bCoulomb && (ic->vdw_modifier == InteractionModifiers::PotShift)));
 829     }
 830
 831     reppow = ic->reppow;
 832
 833     if (tprops[tp].bCoulomb)
 834     {
 835         r1 = ic->rcoulomb_switch;
 836         rc = ic->rcoulomb;
 837     }
 838     else
 839     {
 840         r1 = ic->rvdw_switch;
 841         rc = ic->rvdw;
 842     }
 843     if (bPotentialSwitch)
 844     {
 845         ksw = 1.0 / (gmx::power5(rc - r1));
 846     }
 847     else
 848     {
 849         ksw = 0.0;
 850     }
 851     if (bForceSwitch)
 852     {
 853         if (tp == etabShift)
 854         {
 855             p = 1;
 856         }
 857         else if (tp == etabLJ6Shift)
 858         {
 859             p = 6;
 860         }
 861         else
 862         {
 863             p = reppow;
 864         }
 865
 866         A = p * ((p + 1) * r1 - (p + 4) * rc) / (std::pow(rc, p + 2) * gmx::square(rc - r1));
 867         B = -p * ((p + 1) * r1 - (p + 3) * rc) / (std::pow(rc, p + 2) * gmx::power3(rc - r1));
 868         C = 1.0 / std::pow(rc, p) - A / 3.0 * gmx::power3(rc - r1) - B / 4.0 * gmx::power4(rc - r1);
 869         if (tp == etabLJ6Shift)
 870         {
 871             A = -A;
 872             B = -B;
 873             C = -C;
 874         }
 875         A_3 = A / 3.0;
 876         B_4 = B / 4.0;
 877     }
 878     if (debug)
 879     {
 880         fprintf(debug, "Setting up tables\n");
 881         fflush(debug);
 882     }
 883
 884     if (bPotentialShift)
 885     {
 886         rc2        = rc * rc;
 887         double rc6 = 1.0 / (rc2 * rc2 * rc2);
 888         double rc12;
 889         if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 890         {
 891             rc12 = rc6 * rc6;
 892         }
 893         else
 894         {
 895             rc12 = std::pow(rc, -reppow);
 896         }
 897
 898         switch (tp)
 899         {
 900             case etabLJ6:
 901                 /* Dispersion */
 902                 Vcut = -rc6;
 903                 break;
 904             case etabLJ6Ewald:
 905                 Vcut = -rc6 * exp(-ewclj * ewclj * rc2)
 906                        * (1 + ewclj * ewclj * rc2 + gmx::power4(ewclj) * rc2 * rc2 / 2);
 907                 break;
 908             case etabLJ12:
 909                 /* Repulsion */
 910                 Vcut = rc12;
 911                 break;
 912             case etabCOUL: Vcut = 1.0 / rc; break;
 913             case etabEwald:
 914             case etabEwaldSwitch: Vcut = std::erfc(ewc * rc) / rc; break;
 915             case etabEwaldUser:
 916                 /* Only calculate minus the reciprocal space contribution */
 917                 Vcut = -std::erf(ewc * rc) / rc;
 918                 break;
 919             case etabRF:
 920             case etabRF_ZERO:
 921                 /* No need for preventing the usage of modifiers with RF */
 922                 Vcut = 0.0;
 923                 break;
 924             case etabEXPMIN: Vcut = exp(-rc); break;
 925             default:
 926                 gmx_fatal(FARGS,
 927                           "Cannot apply new potential-shift modifier to interaction type '%s' yet. "
 928                           "(%s,%d)",
 929                           tprops[tp].name,
 930                           __FILE__,
 931                           __LINE__);
 932         }
 933     }
 934
 935     for (i = 0; (i < td->nx); i++)
 936     {
 937         td->x[i] = i / td->tabscale;
 938     }
 939     for (i = td->nx0; (i < td->nx); i++)
 940     {
 941         r  = td->x[i];
 942         r2 = r * r;
 943         r6 = 1.0 / (r2 * r2 * r2);
 944         if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 945         {
 946             r12 = r6 * r6;
 947         }
 948         else
 949         {
 950             r12 = std::pow(r, -reppow);
 951         }
 952         Vtab = 0.0;
 953         Ftab = 0.0;
 954         if (bPotentialSwitch)
 955         {
 956             /* swi is function, swi1 1st derivative and swi2 2nd derivative */
 957             /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
 958              * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
 959              * r1 and rc.
 960              * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
 961              */
 962             if (r <= r1)
 963             {
 964                 swi  = 1.0;
 965                 swi1 = 0.0;
 966             }
 967             else if (r >= rc)
 968             {
 969                 swi  = 0.0;
 970                 swi1 = 0.0;
 971             }
 972             else
 973             {
 974                 swi = 1 - 10 * gmx::power3(r - r1) * ksw * gmx::square(rc - r1)
 975                       + 15 * gmx::power4(r - r1) * ksw * (rc - r1) - 6 * gmx::power5(r - r1) * ksw;
 976                 swi1 = -30 * gmx::square(r - r1) * ksw * gmx::square(rc - r1)
 977                        + 60 * gmx::power3(r - r1) * ksw * (rc - r1) - 30 * gmx::power4(r - r1) * ksw;
 978             }
 979         }
 980         else /* not really needed, but avoids compiler warnings... */
 981         {
 982             swi  = 1.0;
 983             swi1 = 0.0;
 984         }
 985
 986         switch (tp)
 987         {
 988             case etabLJ6:
 989                 /* Dispersion */
 990                 Vtab = -r6;
 991                 Ftab = 6.0 * Vtab / r;
 992                 break;
 993             case etabLJ6Switch:
 994             case etabLJ6Shift:
 995                 /* Dispersion */
 996                 if (r < rc)
 997                 {
 998                     Vtab = -r6;
 999                     Ftab = 6.0 * Vtab / r;
1000                     break;
1001                 }
1002                 break;
1003             case etabLJ12:
1004                 /* Repulsion */
1005                 Vtab = r12;
1006                 Ftab = reppow * Vtab / r;
1007                 break;
1008             case etabLJ12Switch:
1009             case etabLJ12Shift:
1010                 /* Repulsion */
1011                 if (r < rc)
1012                 {
1013                     Vtab = r12;
1014                     Ftab = reppow * Vtab / r;
1015                 }
1016                 break;
1017             case etabCOUL:
1018                 Vtab = 1.0 / r;
1019                 Ftab = 1.0 / r2;
1020                 break;
1021             case etabCOULSwitch:
1022             case etabShift:
1023                 if (r < rc)
1024                 {
1025                     Vtab = 1.0 / r;
1026                     Ftab = 1.0 / r2;
1027                 }
1028                 break;
1029             case etabEwald:
1030             case etabEwaldSwitch:
1031                 Vtab = std::erfc(ewc * r) / r;
1032                 Ftab = std::erfc(ewc * r) / r2 + exp(-(ewc * ewc * r2)) * ewc * M_2_SQRTPI / r;
1033                 break;
1034             case etabEwaldUser:
1035             case etabEwaldUserSwitch:
1036                 /* Only calculate the negative of the reciprocal space contribution */
1037                 Vtab = -std::erf(ewc * r) / r;
1038                 Ftab = -std::erf(ewc * r) / r2 + exp(-(ewc * ewc * r2)) * ewc * M_2_SQRTPI / r;
1039                 break;
1040             case etabLJ6Ewald:
1041                 Vtab = -r6 * exp(-ewclj * ewclj * r2)
1042                        * (1 + ewclj * ewclj * r2 + gmx::power4(ewclj) * r2 * r2 / 2);
1043                 Ftab = 6.0 * Vtab / r
1044                        - r6 * exp(-ewclj * ewclj * r2) * gmx::power5(ewclj) * ewclj * r2 * r2 * r;
1045                 break;
1046             case etabRF:
1047             case etabRF_ZERO:
1048                 Vtab = 1.0 / r + ic->reactionFieldCoefficient * r2 - ic->reactionFieldShift;
1049                 Ftab = 1.0 / r2 - 2 * ic->reactionFieldCoefficient * r;
1050                 if (tp == etabRF_ZERO && r >= rc)
1051                 {
1052                     Vtab = 0;
1053                     Ftab = 0;
1054                 }
1055                 break;
1056             case etabEXPMIN:
1057                 expr = exp(-r);
1058                 Vtab = expr;
1059                 Ftab = expr;
1060                 break;
1061             default:
1062                 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)", tp, __FILE__, __LINE__);
1063         }
1064         if (bForceSwitch)
1065         {
1066             /* Normal coulomb with cut-off correction for potential */
1067             if (r < rc)
1068             {
1069                 Vtab -= C;
1070                 /* If in Shifting range add something to it */
1071                 if (r > r1)
1072                 {
1073                     r12 = (r - r1) * (r - r1);
1074                     r13 = (r - r1) * r12;
1075                     Vtab += -A_3 * r13 - B_4 * r12 * r12;
1076                     Ftab += A * r12 + B * r13;
1077                 }
1078             }
1079             else
1080             {
1081                 /* Make sure interactions are zero outside cutoff with modifiers */
1082                 Vtab = 0;
1083                 Ftab = 0;
1084             }
1085         }
1086         if (bPotentialShift)
1087         {
1088             if (r < rc)
1089             {
1090                 Vtab -= Vcut;
1091             }
1092             else
1093             {
1094                 /* Make sure interactions are zero outside cutoff with modifiers */
1095                 Vtab = 0;
1096                 Ftab = 0;
1097             }
1098         }
1099
1100         if (ETAB_USER(tp))
1101         {
1102             Vtab += td->v[i];
1103             Ftab += td->f[i];
1104         }
1105
1106         if (bPotentialSwitch)
1107         {
1108             if (r >= rc)
1109             {
1110                 /* Make sure interactions are zero outside cutoff with modifiers */
1111                 Vtab = 0;
1112                 Ftab = 0;
1113             }
1114             else if (r > r1)
1115             {
1116                 Ftab = Ftab * swi - Vtab * swi1;
1117                 Vtab = Vtab * swi;
1118             }
1119         }
1120         /* Convert to single precision when we store to mem */
1121         td->v[i] = Vtab;
1122         td->f[i] = Ftab;
1123     }
1124
1125     /* Continue the table linearly from nx0 to 0.
1126      * These values are only required for energy minimization with overlap or TPI.
1127      */
1128     for (i = td->nx0 - 1; i >= 0; i--)
1129     {
1130         td->v[i] = td->v[i + 1] + td->f[i + 1] * (td->x[i + 1] - td->x[i]);
1131         td->f[i] = td->f[i + 1];
1132     }
1133 }
1134
1135 static void set_table_type(int tabsel[], const interaction_const_t* ic, gmx_bool b14only)
1136 {
1137     /* Set the different table indices.
1138      * Coulomb first.
1139      */
1140
1141     CoulombInteractionType eltype;
1142     VanDerWaalsType        vdwtype;
1143
1144     if (b14only)
1145     {
1146         switch (ic->eeltype)
1147         {
1148             case CoulombInteractionType::User:
1149             case CoulombInteractionType::PmeUser:
1150             case CoulombInteractionType::PmeUserSwitch:
1151                 eltype = CoulombInteractionType::User;
1152                 break;
1153             default: eltype = CoulombInteractionType::Cut;
1154         }
1155     }
1156     else
1157     {
1158         eltype = ic->eeltype;
1159     }
1160
1161     switch (eltype)
1162     {
1163         case CoulombInteractionType::Cut: tabsel[etiCOUL] = etabCOUL; break;
1164         case CoulombInteractionType::Poisson: tabsel[etiCOUL] = etabShift; break;
1165         case CoulombInteractionType::Shift:
1166             if (ic->rcoulomb > ic->rcoulomb_switch)
1167             {
1168                 tabsel[etiCOUL] = etabShift;
1169             }
1170             else
1171             {
1172                 tabsel[etiCOUL] = etabCOUL;
1173             }
1174             break;
1175         case CoulombInteractionType::Ewald:
1176         case CoulombInteractionType::Pme:
1177         case CoulombInteractionType::P3mAD: tabsel[etiCOUL] = etabEwald; break;
1178         case CoulombInteractionType::PmeSwitch: tabsel[etiCOUL] = etabEwaldSwitch; break;
1179         case CoulombInteractionType::PmeUser: tabsel[etiCOUL] = etabEwaldUser; break;
1180         case CoulombInteractionType::PmeUserSwitch: tabsel[etiCOUL] = etabEwaldUserSwitch; break;
1181         case CoulombInteractionType::RF:
1182         case CoulombInteractionType::RFZero: tabsel[etiCOUL] = etabRF_ZERO; break;
1183         case CoulombInteractionType::Switch: tabsel[etiCOUL] = etabCOULSwitch; break;
1184         case CoulombInteractionType::User: tabsel[etiCOUL] = etabUSER; break;
1185         default: gmx_fatal(FARGS, "Invalid eeltype %s", enumValueToString(eltype));
1186     }
1187
1188     /* Van der Waals time */
1189     if (ic->useBuckingham && !b14only)
1190     {
1191         tabsel[etiLJ6]  = etabLJ6;
1192         tabsel[etiLJ12] = etabEXPMIN;
1193     }
1194     else
1195     {
1196         if (b14only && ic->vdwtype != VanDerWaalsType::User)
1197         {
1198             vdwtype = VanDerWaalsType::Cut;
1199         }
1200         else
1201         {
1202             vdwtype = ic->vdwtype;
1203         }
1204
1205         switch (vdwtype)
1206         {
1207             case VanDerWaalsType::Switch:
1208                 tabsel[etiLJ6]  = etabLJ6Switch;
1209                 tabsel[etiLJ12] = etabLJ12Switch;
1210                 break;
1211             case VanDerWaalsType::Shift:
1212                 tabsel[etiLJ6]  = etabLJ6Shift;
1213                 tabsel[etiLJ12] = etabLJ12Shift;
1214                 break;
1215             case VanDerWaalsType::User:
1216                 tabsel[etiLJ6]  = etabUSER;
1217                 tabsel[etiLJ12] = etabUSER;
1218                 break;
1219             case VanDerWaalsType::Cut:
1220                 tabsel[etiLJ6]  = etabLJ6;
1221                 tabsel[etiLJ12] = etabLJ12;
1222                 break;
1223             case VanDerWaalsType::Pme:
1224                 tabsel[etiLJ6]  = etabLJ6Ewald;
1225                 tabsel[etiLJ12] = etabLJ12;
1226                 break;
1227             default:
1228                 gmx_fatal(FARGS, "Invalid vdwtype %s in %s line %d", enumValueToString(vdwtype), __FILE__, __LINE__);
1229         }
1230
1231         if (!b14only && ic->vdw_modifier != InteractionModifiers::None)
1232         {
1233             if (ic->vdw_modifier != InteractionModifiers::PotShift && ic->vdwtype != VanDerWaalsType::Cut)
1234             {
1235                 gmx_incons(
1236                         "Potential modifiers other than potential-shift are only implemented for "
1237                         "LJ cut-off");
1238             }
1239
1240             /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
1241              * to the original interaction forms when we fill the table, so we only check cutoffs here.
1242              */
1243             if (ic->vdwtype == VanDerWaalsType::Cut)
1244             {
1245                 switch (ic->vdw_modifier)
1246                 {
1247                     case InteractionModifiers::None:
1248                     case InteractionModifiers::PotShift:
1249                     case InteractionModifiers::ExactCutoff:
1250                         /* No modification */
1251                         break;
1252                     case InteractionModifiers::PotSwitch:
1253                         tabsel[etiLJ6]  = etabLJ6Switch;
1254                         tabsel[etiLJ12] = etabLJ12Switch;
1255                         break;
1256                     case InteractionModifiers::ForceSwitch:
1257                         tabsel[etiLJ6]  = etabLJ6Shift;
1258                         tabsel[etiLJ12] = etabLJ12Shift;
1259                         break;
1260                     default: gmx_incons("Unsupported vdw_modifier");
1261                 }
1262             }
1263         }
1264     }
1265 }
1266
1267 std::unique_ptr<t_forcetable>
1268 make_tables(FILE* fp, const interaction_const_t* ic, const char* fn, real rtab, int flags)
1269 {
1270     t_tabledata* td;
1271     gmx_bool     b14only, useUserTable;
1272     int          nx0, tabsel[etiNR];
1273     real         scalefactor;
1274
1275     auto table = std::make_unique<t_forcetable>(GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
1276                                                 GMX_TABLE_FORMAT_CUBICSPLINE_YFGH);
1277
1278     b14only = ((flags & GMX_MAKETABLES_14ONLY) != 0);
1279
1280     if (flags & GMX_MAKETABLES_FORCEUSER)
1281     {
1282         tabsel[etiCOUL] = etabUSER;
1283         tabsel[etiLJ6]  = etabUSER;
1284         tabsel[etiLJ12] = etabUSER;
1285     }
1286     else
1287     {
1288         set_table_type(tabsel, ic, b14only);
1289     }
1290     snew(td, etiNR);
1291     table->r     = rtab;
1292     table->scale = 0;
1293     table->n     = 0;
1294
1295     table->formatsize    = 4;
1296     table->ninteractions = etiNR;
1297     table->stride        = table->formatsize * table->ninteractions;
1298
1299     /* Check whether we have to read or generate */
1300     useUserTable = FALSE;
1301     for (unsigned int i = 0; (i < etiNR); i++)
1302     {
1303         if (ETAB_USER(tabsel[i]))
1304         {
1305             useUserTable = TRUE;
1306         }
1307     }
1308     if (useUserTable)
1309     {
1310         read_tables(fp, fn, etiNR, 0, td);
1311         if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1312         {
1313             table->n = td[0].nx;
1314         }
1315         else
1316         {
1317             if (td[0].x[td[0].nx - 1] < rtab)
1318             {
1319                 gmx_fatal(FARGS,
1320                           "Tables in file %s not long enough for cut-off:\n"
1321                           "\tshould be at least %f nm\n",
1322                           fn,
1323                           rtab);
1324             }
1325             table->n = gmx::roundToInt(rtab * td[0].tabscale);
1326         }
1327         table->scale = td[0].tabscale;
1328         nx0          = td[0].nx0;
1329     }
1330     else
1331     {
1332         // No tables are read
1333 #if GMX_DOUBLE
1334         table->scale = 2000.0;
1335 #else
1336         table->scale = 500.0;
1337 #endif
1338         table->n = static_cast<int>(rtab * table->scale);
1339         nx0      = 10;
1340     }
1341
1342     /* Each table type (e.g. coul,lj6,lj12) requires four
1343      * numbers per table->n+1 data points. For performance reasons we want
1344      * the table data to be aligned to (at least) a 32-byte boundary.
1345      */
1346     table->data.resize(table->stride * (table->n + 1) * sizeof(real));
1347
1348     for (int k = 0; (k < etiNR); k++)
1349     {
1350         /* Now fill data for tables that have not been read
1351          * or add the Ewald long-range correction for Ewald user tables.
1352          */
1353         if (tabsel[k] != etabUSER)
1354         {
1355             real scale = table->scale;
1356             if (ic->useBuckingham && (ic->buckinghamBMax != 0) && tabsel[k] == etabEXPMIN)
1357             {
1358                 scale /= ic->buckinghamBMax;
1359             }
1360             init_table(table->n, nx0, scale, &(td[k]), !useUserTable);
1361
1362             fill_table(&(td[k]), tabsel[k], ic, b14only);
1363             if (fp)
1364             {
1365                 fprintf(fp,
1366                         "Generated table with %d data points for %s%s.\n"
1367                         "Tabscale = %g points/nm\n",
1368                         td[k].nx,
1369                         b14only ? "1-4 " : "",
1370                         tprops[tabsel[k]].name,
1371                         td[k].tabscale);
1372             }
1373         }
1374
1375         /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
1376          * by including the derivative constants (6.0 or 12.0) in the parameters, since
1377          * we no longer calculate force in most steps. This means the c6/c12 parameters
1378          * have been scaled up, so we need to scale down the table interactions too.
1379          * It comes here since we need to scale user tables too.
1380          */
1381         if (k == etiLJ6)
1382         {
1383             scalefactor = 1.0 / 6.0;
1384         }
1385         else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1386         {
1387             scalefactor = 1.0 / 12.0;
1388         }
1389         else
1390         {
1391             scalefactor = 1.0;
1392         }
1393
1394         copy2table(table->n,
1395                    k * table->formatsize,
1396                    table->stride,
1397                    td[k].x,
1398                    td[k].v,
1399                    td[k].f,
1400                    scalefactor,
1401                    table->data.data());
1402
1403         done_tabledata(&(td[k]));
1404     }
1405     sfree(td);
1406
1407     return table;
1408 }
1409
1410 bondedtable_t make_bonded_table(FILE* fplog, const char* fn, int angle)
1411 {
1412     t_tabledata   td;
1413     int           i;
1414     bondedtable_t tab;
1415     int           stride = 4;
1416
1417     read_tables(fplog, fn, 1, angle, &td);
1418     if (angle > 0)
1419     {
1420         /* Convert the table from degrees to radians */
1421         for (i = 0; i < td.nx; i++)
1422         {
1423             td.x[i] *= gmx::c_deg2Rad;
1424             td.f[i] *= gmx::c_rad2Deg;
1425         }
1426         td.tabscale *= gmx::c_rad2Deg;
1427     }
1428     tab.n     = td.nx;
1429     tab.scale = td.tabscale;
1430     tab.data.resize(tab.n * stride);
1431     copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data.data());
1432     done_tabledata(&td);
1433
1434     return tab;
1435 }
1436
1437 std::unique_ptr<t_forcetable>
1438 makeDispersionCorrectionTable(FILE* fp, const interaction_const_t* ic, real rtab, const char* tabfn)
1439 {
1440     GMX_RELEASE_ASSERT(ic->vdwtype != VanDerWaalsType::User || tabfn,
1441                        "With VdW user tables we need a table file name");
1442
1443     std::unique_ptr<t_forcetable> fullTable = make_tables(fp, ic, tabfn, rtab, 0);
1444     /* Copy the contents of the table to one that has just dispersion
1445      * and repulsion, to improve cache performance. We want the table
1446      * data to be aligned to 32-byte boundaries.
1447      */
1448     std::unique_ptr<t_forcetable> dispersionCorrectionTable =
1449             std::make_unique<t_forcetable>(GMX_TABLE_INTERACTION_VDWREP_VDWDISP, fullTable->format);
1450     dispersionCorrectionTable->r             = fullTable->r;
1451     dispersionCorrectionTable->n             = fullTable->n;
1452     dispersionCorrectionTable->scale         = fullTable->scale;
1453     dispersionCorrectionTable->formatsize    = fullTable->formatsize;
1454     dispersionCorrectionTable->ninteractions = 2;
1455     dispersionCorrectionTable->stride =
1456             dispersionCorrectionTable->formatsize * dispersionCorrectionTable->ninteractions;
1457     dispersionCorrectionTable->data.resize(dispersionCorrectionTable->stride
1458                                            * (dispersionCorrectionTable->n + 1));
1459
1460     for (int i = 0; i <= fullTable->n; i++)
1461     {
1462         for (int j = 0; j < 8; j++)
1463         {
1464             dispersionCorrectionTable->data[8 * i + j] = fullTable->data[12 * i + 4 + j];
1465         }
1466     }
1467
1468     return dispersionCorrectionTable;
1469 }
1470
1471 t_forcetable::t_forcetable(enum gmx_table_interaction interaction, enum gmx_table_format format) :
1472     interaction(interaction),
1473     format(format),
1474     r(0),
1475     n(0),
1476     scale(0),
1477     formatsize(0),
1478     ninteractions(0),
1479     stride(0)
1480 {
1481 }
1482
1483 t_forcetable::~t_forcetable() = default;