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