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),
 436     nx0(nx0),
 437     tabscale(tabscale)
 438 {
 439     if (bAlloc)
 440     {
 441         x.resize(nx);
 442         v.resize(nx);
 443         f.resize(nx);
 444     }
 445 }
 446
 447 static void spline_forces(int nx, double h, const double v[], gmx_bool bS3, gmx_bool bE3, double f[])
 448 {
 449     int    start, end, i;
 450     double v3, b_s, b_e, b;
 451     double beta;
 452
 453     /* Formulas can be found in:
 454      * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
 455      */
 456
 457     if (nx < 4 && (bS3 || bE3))
 458     {
 459         gmx_fatal(FARGS,
 460                   "Can not generate splines with third derivative boundary conditions with less "
 461                   "than 4 (%d) points",
 462                   nx);
 463     }
 464
 465     /* To make life easy we initially set the spacing to 1
 466      * and correct for this at the end.
 467      */
 468     if (bS3)
 469     {
 470         /* Fit V''' at the start */
 471         v3 = v[3] - 3 * v[2] + 3 * v[1] - v[0];
 472         if (debug)
 473         {
 474             fprintf(debug, "The left third derivative is %g\n", v3 / (h * h * h));
 475         }
 476         b_s   = 2 * (v[1] - v[0]) + v3 / 6;
 477         start = 0;
 478
 479         if (FALSE)
 480         {
 481             /* Fit V'' at the start */
 482             real v2;
 483
 484             v2 = -v[3] + 4 * v[2] - 5 * v[1] + 2 * v[0];
 485             /* v2  = v[2] - 2*v[1] + v[0]; */
 486             if (debug)
 487             {
 488                 fprintf(debug, "The left second derivative is %g\n", v2 / (h * h));
 489             }
 490             b_s   = 3 * (v[1] - v[0]) - v2 / 2;
 491             start = 0;
 492         }
 493     }
 494     else
 495     {
 496         b_s   = 3 * (v[2] - v[0]) + f[0] * h;
 497         start = 1;
 498     }
 499     if (bE3)
 500     {
 501         /* Fit V''' at the end */
 502         v3 = v[nx - 1] - 3 * v[nx - 2] + 3 * v[nx - 3] - v[nx - 4];
 503         if (debug)
 504         {
 505             fprintf(debug, "The right third derivative is %g\n", v3 / (h * h * h));
 506         }
 507         b_e = 2 * (v[nx - 1] - v[nx - 2]) + v3 / 6;
 508         end = nx;
 509     }
 510     else
 511     {
 512         /* V'=0 at the end */
 513         b_e = 3 * (v[nx - 1] - v[nx - 3]) + f[nx - 1] * h;
 514         end = nx - 1;
 515     }
 516
 517     std::vector<double> gamma(nx);
 518     beta = (bS3 ? 1 : 4);
 519
 520     /* For V'' fitting */
 521     /* beta = (bS3 ? 2 : 4); */
 522
 523     f[start] = b_s / beta;
 524     for (i = start + 1; i < end; i++)
 525     {
 526         gamma[i] = 1 / beta;
 527         beta     = 4 - gamma[i];
 528         b        = 3 * (v[i + 1] - v[i - 1]);
 529         f[i]     = (b - f[i - 1]) / beta;
 530     }
 531     gamma[end - 1] = 1 / beta;
 532     beta           = (bE3 ? 1 : 4) - gamma[end - 1];
 533     f[end - 1]     = (b_e - f[end - 2]) / beta;
 534
 535     for (i = end - 2; i >= start; i--)
 536     {
 537         f[i] -= gamma[i + 1] * f[i + 1];
 538     }
 539
 540     /* Correct for the minus sign and the spacing */
 541     for (i = start; i < end; i++)
 542     {
 543         f[i] = -f[i] / h;
 544     }
 545 }
 546
 547 static void set_forces(FILE* fp, int angle, int nx, double h, double v[], double f[], int table)
 548 {
 549     int start, end;
 550
 551     if (angle == 2)
 552     {
 553         gmx_fatal(FARGS, "Force generation for dihedral tables is not (yet) implemented");
 554     }
 555
 556     start = 0;
 557     while (v[start] == 0)
 558     {
 559         start++;
 560     }
 561
 562     end = nx;
 563     while (v[end - 1] == 0)
 564     {
 565         end--;
 566     }
 567     if (end > nx - 2)
 568     {
 569         end = nx;
 570     }
 571     else
 572     {
 573         end++;
 574     }
 575
 576     if (fp)
 577     {
 578         fprintf(fp,
 579                 "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
 580                 table + 1,
 581                 start * h,
 582                 end == nx ? "V'''" : "V'=0",
 583                 (end - 1) * h);
 584     }
 585     spline_forces(end - start, h, v + start, TRUE, end == nx, f + start);
 586 }
 587
 588 static std::vector<t_tabledata> read_tables(FILE* fp, const char* filename, int ntab, int angle)
 589 {
 590     char     buf[STRLEN];
 591     double   start, end, dx0, dx1, ssd, vm, vp, f, numf;
 592     int      k, i, nx0 = 0, nny, ns;
 593     gmx_bool bAllZero, bZeroV, bZeroF;
 594     double   tabscale;
 595
 596     nny               = 2 * ntab + 1;
 597     std::string libfn = gmx::findLibraryFile(filename);
 598     gmx::MultiDimArray<std::vector<double>, gmx::dynamicExtents2D> xvgData    = readXvgData(libfn);
 599     int                                                            numColumns = xvgData.extent(0);
 600     if (numColumns != nny)
 601     {
 602         gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d", libfn.c_str(), numColumns, nny);
 603     }
 604     int numRows = xvgData.extent(1);
 605
 606     const auto& yy = xvgData.asView();
 607     if (angle == 0)
 608     {
 609         if (yy[0][0] != 0.0)
 610         {
 611             gmx_fatal(FARGS,
 612                       "The first distance in file %s is %f nm instead of %f nm",
 613                       libfn.c_str(),
 614                       yy[0][0],
 615                       0.0);
 616         }
 617     }
 618     else
 619     {
 620         if (angle == 1)
 621         {
 622             start = 0.0;
 623         }
 624         else
 625         {
 626             start = -180.0;
 627         }
 628         end = 180.0;
 629         if (yy[0][0] != start || yy[0][numRows - 1] != end)
 630         {
 631             gmx_fatal(FARGS,
 632                       "The angles in file %s should go from %f to %f instead of %f to %f\n",
 633                       libfn.c_str(),
 634                       start,
 635                       end,
 636                       yy[0][0],
 637                       yy[0][numRows - 1]);
 638         }
 639     }
 640
 641     tabscale = (numRows - 1) / (yy[0][numRows - 1] - yy[0][0]);
 642
 643     if (fp)
 644     {
 645         fprintf(fp, "Read user tables from %s with %d data points.\n", libfn.c_str(), numRows);
 646         if (angle == 0)
 647         {
 648             fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
 649         }
 650     }
 651
 652     bAllZero = TRUE;
 653     for (k = 0; k < ntab; k++)
 654     {
 655         bZeroV = TRUE;
 656         bZeroF = TRUE;
 657         for (i = 0; (i < numRows); i++)
 658         {
 659             if (i >= 2)
 660             {
 661                 dx0 = yy[0][i - 1] - yy[0][i - 2];
 662                 dx1 = yy[0][i] - yy[0][i - 1];
 663                 /* Check for 1% deviation in spacing */
 664                 if (fabs(dx1 - dx0) >= 0.005 * (fabs(dx0) + fabs(dx1)))
 665                 {
 666                     gmx_fatal(FARGS,
 667                               "In table file '%s' the x values are not equally spaced: %f %f %f",
 668                               filename,
 669                               yy[0][i - 2],
 670                               yy[0][i - 1],
 671                               yy[0][i]);
 672                 }
 673             }
 674             if (yy[1 + k * 2][i] != 0)
 675             {
 676                 bZeroV = FALSE;
 677                 if (bAllZero)
 678                 {
 679                     bAllZero = FALSE;
 680                     nx0      = i;
 681                 }
 682                 if (yy[1 + k * 2][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2][i] < -0.01 * GMX_REAL_MAX)
 683                 {
 684                     gmx_fatal(FARGS, "Out of range potential value %g in file '%s'", yy[1 + k * 2][i], filename);
 685                 }
 686             }
 687             if (yy[1 + k * 2 + 1][i] != 0)
 688             {
 689                 bZeroF = FALSE;
 690                 if (bAllZero)
 691                 {
 692                     bAllZero = FALSE;
 693                     nx0      = i;
 694                 }
 695                 if (yy[1 + k * 2 + 1][i] > 0.01 * GMX_REAL_MAX || yy[1 + k * 2 + 1][i] < -0.01 * GMX_REAL_MAX)
 696                 {
 697                     gmx_fatal(FARGS, "Out of range force value %g in file '%s'", yy[1 + k * 2 + 1][i], filename);
 698                 }
 699             }
 700         }
 701
 702         if (!bZeroV && bZeroF)
 703         {
 704             set_forces(fp, angle, numRows, 1 / tabscale, yy[1 + k * 2].data(), yy[1 + k * 2 + 1].data(), k);
 705         }
 706         else
 707         {
 708             /* Check if the second column is close to minus the numerical
 709              * derivative of the first column.
 710              */
 711             ssd = 0;
 712             ns  = 0;
 713             for (i = 1; (i < numRows - 1); i++)
 714             {
 715                 vm = yy[1 + 2 * k][i - 1];
 716                 vp = yy[1 + 2 * k][i + 1];
 717                 f  = yy[1 + 2 * k + 1][i];
 718                 if (vm != 0 && vp != 0 && f != 0)
 719                 {
 720                     /* Take the centered difference */
 721                     numf = -(vp - vm) * 0.5 * tabscale;
 722                     if (f + numf != 0)
 723                     {
 724                         ssd += fabs(2 * (f - numf) / (f + numf));
 725                     }
 726                     ns++;
 727                 }
 728             }
 729             if (ns > 0)
 730             {
 731                 ssd /= ns;
 732                 sprintf(buf,
 733                         "For the %d non-zero entries for table %d in %s the forces deviate on "
 734                         "average %" PRId64
 735                         "%% from minus the numerical derivative of the potential\n",
 736                         ns,
 737                         k,
 738                         libfn.c_str(),
 739                         gmx::roundToInt64(100 * ssd));
 740                 if (debug)
 741                 {
 742                     fprintf(debug, "%s", buf);
 743                 }
 744                 if (ssd > 0.2)
 745                 {
 746                     if (fp)
 747                     {
 748                         fprintf(fp, "\nWARNING: %s\n", buf);
 749                     }
 750                     fprintf(stderr, "\nWARNING: %s\n", buf);
 751                 }
 752             }
 753         }
 754     }
 755     if (bAllZero && fp)
 756     {
 757         fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn.c_str());
 758     }
 759
 760     std::vector<t_tabledata> td;
 761     for (k = 0; (k < ntab); k++)
 762     {
 763         td.emplace_back(t_tabledata(numRows, nx0, tabscale, true));
 764         for (i = 0; (i < numRows); i++)
 765         {
 766             td[k].x[i] = yy[0][i];
 767             td[k].v[i] = yy[2 * k + 1][i];
 768             td[k].f[i] = yy[2 * k + 2][i];
 769         }
 770     }
 771     return td;
 772 }
 773
 774 static void fill_table(t_tabledata* td, int tp, const interaction_const_t* ic, gmx_bool b14only)
 775 {
 776     /* Fill the table according to the formulas in the manual.
 777      * In principle, we only need the potential and the second
 778      * derivative, but then we would have to do lots of calculations
 779      * in the inner loop. By precalculating some terms (see manual)
 780      * we get better eventual performance, despite a larger table.
 781      *
 782      * Since some of these higher-order terms are very small,
 783      * we always use double precision to calculate them here, in order
 784      * to avoid unnecessary loss of precision.
 785      */
 786     int    i;
 787     double reppow, p;
 788     double r1, rc, r12, r13;
 789     double r, r2, r6, rc2;
 790     double expr, Vtab, Ftab;
 791     /* Parameters for David's function */
 792     double A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
 793     /* Parameters for the switching function */
 794     double ksw, swi, swi1;
 795     /* Temporary parameters */
 796     gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
 797     double   ewc   = ic->ewaldcoeff_q;
 798     double   ewclj = ic->ewaldcoeff_lj;
 799     double   Vcut  = 0;
 800
 801     if (b14only)
 802     {
 803         bPotentialSwitch = FALSE;
 804         bForceSwitch     = FALSE;
 805         bPotentialShift  = FALSE;
 806     }
 807     else
 808     {
 809         bPotentialSwitch =
 810                 ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) || (tp == etabCOULSwitch)
 811                  || (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch)
 812                  || (tprops[tp].bCoulomb && (ic->coulomb_modifier == InteractionModifiers::PotSwitch))
 813                  || (!tprops[tp].bCoulomb && (ic->vdw_modifier == InteractionModifiers::PotSwitch)));
 814         bForceSwitch =
 815                 ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) || (tp == etabShift)
 816                  || (tprops[tp].bCoulomb && (ic->coulomb_modifier == InteractionModifiers::ForceSwitch))
 817                  || (!tprops[tp].bCoulomb && (ic->vdw_modifier == InteractionModifiers::ForceSwitch)));
 818         bPotentialShift =
 819                 ((tprops[tp].bCoulomb && (ic->coulomb_modifier == InteractionModifiers::PotShift))
 820                  || (!tprops[tp].bCoulomb && (ic->vdw_modifier == InteractionModifiers::PotShift)));
 821     }
 822
 823     reppow = ic->reppow;
 824
 825     if (tprops[tp].bCoulomb)
 826     {
 827         r1 = ic->rcoulomb_switch;
 828         rc = ic->rcoulomb;
 829     }
 830     else
 831     {
 832         r1 = ic->rvdw_switch;
 833         rc = ic->rvdw;
 834     }
 835     if (bPotentialSwitch)
 836     {
 837         ksw = 1.0 / (gmx::power5(rc - r1));
 838     }
 839     else
 840     {
 841         ksw = 0.0;
 842     }
 843     if (bForceSwitch)
 844     {
 845         if (tp == etabShift)
 846         {
 847             p = 1;
 848         }
 849         else if (tp == etabLJ6Shift)
 850         {
 851             p = 6;
 852         }
 853         else
 854         {
 855             p = reppow;
 856         }
 857
 858         A = p * ((p + 1) * r1 - (p + 4) * rc) / (std::pow(rc, p + 2) * gmx::square(rc - r1));
 859         B = -p * ((p + 1) * r1 - (p + 3) * rc) / (std::pow(rc, p + 2) * gmx::power3(rc - r1));
 860         C = 1.0 / std::pow(rc, p) - A / 3.0 * gmx::power3(rc - r1) - B / 4.0 * gmx::power4(rc - r1);
 861         if (tp == etabLJ6Shift)
 862         {
 863             A = -A;
 864             B = -B;
 865             C = -C;
 866         }
 867         A_3 = A / 3.0;
 868         B_4 = B / 4.0;
 869     }
 870     if (debug)
 871     {
 872         fprintf(debug, "Setting up tables\n");
 873         fflush(debug);
 874     }
 875
 876     if (bPotentialShift)
 877     {
 878         rc2        = rc * rc;
 879         double rc6 = 1.0 / (rc2 * rc2 * rc2);
 880         double rc12;
 881         if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 882         {
 883             rc12 = rc6 * rc6;
 884         }
 885         else
 886         {
 887             rc12 = std::pow(rc, -reppow);
 888         }
 889
 890         switch (tp)
 891         {
 892             case etabLJ6:
 893                 /* Dispersion */
 894                 Vcut = -rc6;
 895                 break;
 896             case etabLJ6Ewald:
 897                 Vcut = -rc6 * exp(-ewclj * ewclj * rc2)
 898                        * (1 + ewclj * ewclj * rc2 + gmx::power4(ewclj) * rc2 * rc2 / 2);
 899                 break;
 900             case etabLJ12:
 901                 /* Repulsion */
 902                 Vcut = rc12;
 903                 break;
 904             case etabCOUL: Vcut = 1.0 / rc; break;
 905             case etabEwald:
 906             case etabEwaldSwitch: Vcut = std::erfc(ewc * rc) / rc; break;
 907             case etabEwaldUser:
 908                 /* Only calculate minus the reciprocal space contribution */
 909                 Vcut = -std::erf(ewc * rc) / rc;
 910                 break;
 911             case etabRF:
 912             case etabRF_ZERO:
 913                 /* No need for preventing the usage of modifiers with RF */
 914                 Vcut = 0.0;
 915                 break;
 916             case etabEXPMIN: Vcut = exp(-rc); break;
 917             default:
 918                 gmx_fatal(FARGS,
 919                           "Cannot apply new potential-shift modifier to interaction type '%s' yet. "
 920                           "(%s,%d)",
 921                           tprops[tp].name,
 922                           __FILE__,
 923                           __LINE__);
 924         }
 925     }
 926
 927     for (i = 0; (i < td->nx); i++)
 928     {
 929         td->x[i] = i / td->tabscale;
 930     }
 931     for (i = td->nx0; (i < td->nx); i++)
 932     {
 933         r  = td->x[i];
 934         r2 = r * r;
 935         r6 = 1.0 / (r2 * r2 * r2);
 936         if (gmx_within_tol(reppow, 12.0, 10 * GMX_DOUBLE_EPS))
 937         {
 938             r12 = r6 * r6;
 939         }
 940         else
 941         {
 942             r12 = std::pow(r, -reppow);
 943         }
 944         Vtab = 0.0;
 945         Ftab = 0.0;
 946         if (bPotentialSwitch)
 947         {
 948             /* swi is function, swi1 1st derivative and swi2 2nd derivative */
 949             /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
 950              * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
 951              * r1 and rc.
 952              * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
 953              */
 954             if (r <= r1)
 955             {
 956                 swi  = 1.0;
 957                 swi1 = 0.0;
 958             }
 959             else if (r >= rc)
 960             {
 961                 swi  = 0.0;
 962                 swi1 = 0.0;
 963             }
 964             else
 965             {
 966                 swi = 1 - 10 * gmx::power3(r - r1) * ksw * gmx::square(rc - r1)
 967                       + 15 * gmx::power4(r - r1) * ksw * (rc - r1) - 6 * gmx::power5(r - r1) * ksw;
 968                 swi1 = -30 * gmx::square(r - r1) * ksw * gmx::square(rc - r1)
 969                        + 60 * gmx::power3(r - r1) * ksw * (rc - r1) - 30 * gmx::power4(r - r1) * ksw;
 970             }
 971         }
 972         else /* not really needed, but avoids compiler warnings... */
 973         {
 974             swi  = 1.0;
 975             swi1 = 0.0;
 976         }
 977
 978         switch (tp)
 979         {
 980             case etabLJ6:
 981                 /* Dispersion */
 982                 Vtab = -r6;
 983                 Ftab = 6.0 * Vtab / r;
 984                 break;
 985             case etabLJ6Switch:
 986             case etabLJ6Shift:
 987                 /* Dispersion */
 988                 if (r < rc)
 989                 {
 990                     Vtab = -r6;
 991                     Ftab = 6.0 * Vtab / r;
 992                     break;
 993                 }
 994                 break;
 995             case etabLJ12:
 996                 /* Repulsion */
 997                 Vtab = r12;
 998                 Ftab = reppow * Vtab / r;
 999                 break;
1000             case etabLJ12Switch:
1001             case etabLJ12Shift:
1002                 /* Repulsion */
1003                 if (r < rc)
1004                 {
1005                     Vtab = r12;
1006                     Ftab = reppow * Vtab / r;
1007                 }
1008                 break;
1009             case etabCOUL:
1010                 Vtab = 1.0 / r;
1011                 Ftab = 1.0 / r2;
1012                 break;
1013             case etabCOULSwitch:
1014             case etabShift:
1015                 if (r < rc)
1016                 {
1017                     Vtab = 1.0 / r;
1018                     Ftab = 1.0 / r2;
1019                 }
1020                 break;
1021             case etabEwald:
1022             case etabEwaldSwitch:
1023                 Vtab = std::erfc(ewc * r) / r;
1024                 Ftab = std::erfc(ewc * r) / r2 + exp(-(ewc * ewc * r2)) * ewc * M_2_SQRTPI / r;
1025                 break;
1026             case etabEwaldUser:
1027             case etabEwaldUserSwitch:
1028                 /* Only calculate the negative of the reciprocal space contribution */
1029                 Vtab = -std::erf(ewc * r) / r;
1030                 Ftab = -std::erf(ewc * r) / r2 + exp(-(ewc * ewc * r2)) * ewc * M_2_SQRTPI / r;
1031                 break;
1032             case etabLJ6Ewald:
1033                 Vtab = -r6 * exp(-ewclj * ewclj * r2)
1034                        * (1 + ewclj * ewclj * r2 + gmx::power4(ewclj) * r2 * r2 / 2);
1035                 Ftab = 6.0 * Vtab / r
1036                        - r6 * exp(-ewclj * ewclj * r2) * gmx::power5(ewclj) * ewclj * r2 * r2 * r;
1037                 break;
1038             case etabRF:
1039             case etabRF_ZERO:
1040                 Vtab = 1.0 / r + ic->reactionFieldCoefficient * r2 - ic->reactionFieldShift;
1041                 Ftab = 1.0 / r2 - 2 * ic->reactionFieldCoefficient * r;
1042                 if (tp == etabRF_ZERO && r >= rc)
1043                 {
1044                     Vtab = 0;
1045                     Ftab = 0;
1046                 }
1047                 break;
1048             case etabEXPMIN:
1049                 expr = exp(-r);
1050                 Vtab = expr;
1051                 Ftab = expr;
1052                 break;
1053             default:
1054                 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)", tp, __FILE__, __LINE__);
1055         }
1056         if (bForceSwitch)
1057         {
1058             /* Normal coulomb with cut-off correction for potential */
1059             if (r < rc)
1060             {
1061                 Vtab -= C;
1062                 /* If in Shifting range add something to it */
1063                 if (r > r1)
1064                 {
1065                     r12 = (r - r1) * (r - r1);
1066                     r13 = (r - r1) * r12;
1067                     Vtab += -A_3 * r13 - B_4 * r12 * r12;
1068                     Ftab += A * r12 + B * r13;
1069                 }
1070             }
1071             else
1072             {
1073                 /* Make sure interactions are zero outside cutoff with modifiers */
1074                 Vtab = 0;
1075                 Ftab = 0;
1076             }
1077         }
1078         if (bPotentialShift)
1079         {
1080             if (r < rc)
1081             {
1082                 Vtab -= Vcut;
1083             }
1084             else
1085             {
1086                 /* Make sure interactions are zero outside cutoff with modifiers */
1087                 Vtab = 0;
1088                 Ftab = 0;
1089             }
1090         }
1091
1092         if (ETAB_USER(tp))
1093         {
1094             Vtab += td->v[i];
1095             Ftab += td->f[i];
1096         }
1097
1098         if (bPotentialSwitch)
1099         {
1100             if (r >= rc)
1101             {
1102                 /* Make sure interactions are zero outside cutoff with modifiers */
1103                 Vtab = 0;
1104                 Ftab = 0;
1105             }
1106             else if (r > r1)
1107             {
1108                 Ftab = Ftab * swi - Vtab * swi1;
1109                 Vtab = Vtab * swi;
1110             }
1111         }
1112         /* Convert to single precision when we store to mem */
1113         td->v[i] = Vtab;
1114         td->f[i] = Ftab;
1115     }
1116
1117     /* Continue the table linearly from nx0 to 0.
1118      * These values are only required for energy minimization with overlap or TPI.
1119      */
1120     for (i = td->nx0 - 1; i >= 0; i--)
1121     {
1122         td->v[i] = td->v[i + 1] + td->f[i + 1] * (td->x[i + 1] - td->x[i]);
1123         td->f[i] = td->f[i + 1];
1124     }
1125 }
1126
1127 static void set_table_type(int tabsel[], const interaction_const_t* ic, gmx_bool b14only)
1128 {
1129     /* Set the different table indices.
1130      * Coulomb first.
1131      */
1132
1133     CoulombInteractionType eltype;
1134     VanDerWaalsType        vdwtype;
1135
1136     if (b14only)
1137     {
1138         switch (ic->eeltype)
1139         {
1140             case CoulombInteractionType::User:
1141             case CoulombInteractionType::PmeUser:
1142             case CoulombInteractionType::PmeUserSwitch:
1143                 eltype = CoulombInteractionType::User;
1144                 break;
1145             default: eltype = CoulombInteractionType::Cut;
1146         }
1147     }
1148     else
1149     {
1150         eltype = ic->eeltype;
1151     }
1152
1153     switch (eltype)
1154     {
1155         case CoulombInteractionType::Cut: tabsel[etiCOUL] = etabCOUL; break;
1156         case CoulombInteractionType::Poisson: tabsel[etiCOUL] = etabShift; break;
1157         case CoulombInteractionType::Shift:
1158             if (ic->rcoulomb > ic->rcoulomb_switch)
1159             {
1160                 tabsel[etiCOUL] = etabShift;
1161             }
1162             else
1163             {
1164                 tabsel[etiCOUL] = etabCOUL;
1165             }
1166             break;
1167         case CoulombInteractionType::Ewald:
1168         case CoulombInteractionType::Pme:
1169         case CoulombInteractionType::P3mAD: tabsel[etiCOUL] = etabEwald; break;
1170         case CoulombInteractionType::PmeSwitch: tabsel[etiCOUL] = etabEwaldSwitch; break;
1171         case CoulombInteractionType::PmeUser: tabsel[etiCOUL] = etabEwaldUser; break;
1172         case CoulombInteractionType::PmeUserSwitch: tabsel[etiCOUL] = etabEwaldUserSwitch; break;
1173         case CoulombInteractionType::RF:
1174         case CoulombInteractionType::RFZero: tabsel[etiCOUL] = etabRF_ZERO; break;
1175         case CoulombInteractionType::Switch: tabsel[etiCOUL] = etabCOULSwitch; break;
1176         case CoulombInteractionType::User: tabsel[etiCOUL] = etabUSER; break;
1177         default: gmx_fatal(FARGS, "Invalid eeltype %s", enumValueToString(eltype));
1178     }
1179
1180     /* Van der Waals time */
1181     if (ic->useBuckingham && !b14only)
1182     {
1183         tabsel[etiLJ6]  = etabLJ6;
1184         tabsel[etiLJ12] = etabEXPMIN;
1185     }
1186     else
1187     {
1188         if (b14only && ic->vdwtype != VanDerWaalsType::User)
1189         {
1190             vdwtype = VanDerWaalsType::Cut;
1191         }
1192         else
1193         {
1194             vdwtype = ic->vdwtype;
1195         }
1196
1197         switch (vdwtype)
1198         {
1199             case VanDerWaalsType::Switch:
1200                 tabsel[etiLJ6]  = etabLJ6Switch;
1201                 tabsel[etiLJ12] = etabLJ12Switch;
1202                 break;
1203             case VanDerWaalsType::Shift:
1204                 tabsel[etiLJ6]  = etabLJ6Shift;
1205                 tabsel[etiLJ12] = etabLJ12Shift;
1206                 break;
1207             case VanDerWaalsType::User:
1208                 tabsel[etiLJ6]  = etabUSER;
1209                 tabsel[etiLJ12] = etabUSER;
1210                 break;
1211             case VanDerWaalsType::Cut:
1212                 tabsel[etiLJ6]  = etabLJ6;
1213                 tabsel[etiLJ12] = etabLJ12;
1214                 break;
1215             case VanDerWaalsType::Pme:
1216                 tabsel[etiLJ6]  = etabLJ6Ewald;
1217                 tabsel[etiLJ12] = etabLJ12;
1218                 break;
1219             default:
1220                 gmx_fatal(FARGS, "Invalid vdwtype %s in %s line %d", enumValueToString(vdwtype), __FILE__, __LINE__);
1221         }
1222
1223         if (!b14only && ic->vdw_modifier != InteractionModifiers::None)
1224         {
1225             if (ic->vdw_modifier != InteractionModifiers::PotShift && ic->vdwtype != VanDerWaalsType::Cut)
1226             {
1227                 gmx_incons(
1228                         "Potential modifiers other than potential-shift are only implemented for "
1229                         "LJ cut-off");
1230             }
1231
1232             /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
1233              * to the original interaction forms when we fill the table, so we only check cutoffs here.
1234              */
1235             if (ic->vdwtype == VanDerWaalsType::Cut)
1236             {
1237                 switch (ic->vdw_modifier)
1238                 {
1239                     case InteractionModifiers::None:
1240                     case InteractionModifiers::PotShift:
1241                     case InteractionModifiers::ExactCutoff:
1242                         /* No modification */
1243                         break;
1244                     case InteractionModifiers::PotSwitch:
1245                         tabsel[etiLJ6]  = etabLJ6Switch;
1246                         tabsel[etiLJ12] = etabLJ12Switch;
1247                         break;
1248                     case InteractionModifiers::ForceSwitch:
1249                         tabsel[etiLJ6]  = etabLJ6Shift;
1250                         tabsel[etiLJ12] = etabLJ12Shift;
1251                         break;
1252                     default: gmx_incons("Unsupported vdw_modifier");
1253                 }
1254             }
1255         }
1256     }
1257 }
1258
1259 std::unique_ptr<t_forcetable>
1260 make_tables(FILE* fp, const interaction_const_t* ic, const char* fn, real rtab, int flags)
1261 {
1262     gmx_bool b14only, useUserTable;
1263     int      nx0, tabsel[etiNR];
1264     real     scalefactor;
1265
1266     auto table = std::make_unique<t_forcetable>(GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP,
1267                                                 GMX_TABLE_FORMAT_CUBICSPLINE_YFGH);
1268
1269     b14only = ((flags & GMX_MAKETABLES_14ONLY) != 0);
1270
1271     if (flags & GMX_MAKETABLES_FORCEUSER)
1272     {
1273         tabsel[etiCOUL] = etabUSER;
1274         tabsel[etiLJ6]  = etabUSER;
1275         tabsel[etiLJ12] = etabUSER;
1276     }
1277     else
1278     {
1279         set_table_type(tabsel, ic, b14only);
1280     }
1281     std::vector<t_tabledata> td;
1282     table->r     = rtab;
1283     table->scale = 0;
1284     table->n     = 0;
1285
1286     table->formatsize    = 4;
1287     table->ninteractions = etiNR;
1288     table->stride        = table->formatsize * table->ninteractions;
1289
1290     /* Check whether we have to read or generate */
1291     useUserTable = FALSE;
1292     for (unsigned int i = 0; (i < etiNR); i++)
1293     {
1294         if (ETAB_USER(tabsel[i]))
1295         {
1296             useUserTable = TRUE;
1297         }
1298     }
1299     if (useUserTable)
1300     {
1301         td = read_tables(fp, fn, etiNR, 0);
1302         if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1303         {
1304             table->n = td[0].nx;
1305         }
1306         else
1307         {
1308             if (td[0].x[td[0].nx - 1] < rtab)
1309             {
1310                 gmx_fatal(FARGS,
1311                           "Tables in file %s not long enough for cut-off:\n"
1312                           "\tshould be at least %f nm\n",
1313                           fn,
1314                           rtab);
1315             }
1316             table->n = gmx::roundToInt(rtab * td[0].tabscale);
1317         }
1318         table->scale = td[0].tabscale;
1319         nx0          = td[0].nx0;
1320     }
1321     else
1322     {
1323         td.resize(etiNR);
1324         // No tables are read
1325 #if GMX_DOUBLE
1326         table->scale = 2000.0;
1327 #else
1328         table->scale = 500.0;
1329 #endif
1330         table->n = static_cast<int>(rtab * table->scale);
1331         nx0      = 10;
1332     }
1333
1334     /* Each table type (e.g. coul,lj6,lj12) requires four
1335      * numbers per table->n+1 data points. For performance reasons we want
1336      * the table data to be aligned to (at least) a 32-byte boundary.
1337      */
1338     table->data.resize(table->stride * (table->n + 1) * sizeof(real));
1339
1340     for (int k = 0; (k < etiNR); k++)
1341     {
1342         /* Now fill data for tables that have not been read
1343          * or add the Ewald long-range correction for Ewald user tables.
1344          */
1345         if (tabsel[k] != etabUSER)
1346         {
1347             real scale = table->scale;
1348             if (ic->useBuckingham && (ic->buckinghamBMax != 0) && tabsel[k] == etabEXPMIN)
1349             {
1350                 scale /= ic->buckinghamBMax;
1351             }
1352             td[k] = t_tabledata(table->n, nx0, scale, !useUserTable);
1353
1354             fill_table(&(td[k]), tabsel[k], ic, b14only);
1355             if (fp)
1356             {
1357                 fprintf(fp,
1358                         "Generated table with %d data points for %s%s.\n"
1359                         "Tabscale = %g points/nm\n",
1360                         td[k].nx,
1361                         b14only ? "1-4 " : "",
1362                         tprops[tabsel[k]].name,
1363                         td[k].tabscale);
1364             }
1365         }
1366
1367         /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
1368          * by including the derivative constants (6.0 or 12.0) in the parameters, since
1369          * we no longer calculate force in most steps. This means the c6/c12 parameters
1370          * have been scaled up, so we need to scale down the table interactions too.
1371          * It comes here since we need to scale user tables too.
1372          */
1373         if (k == etiLJ6)
1374         {
1375             scalefactor = 1.0 / 6.0;
1376         }
1377         else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1378         {
1379             scalefactor = 1.0 / 12.0;
1380         }
1381         else
1382         {
1383             scalefactor = 1.0;
1384         }
1385
1386         copy2table(table->n, k * table->formatsize, table->stride, td[k].x, td[k].v, td[k].f, scalefactor, table->data);
1387     }
1388
1389     return table;
1390 }
1391
1392 bondedtable_t make_bonded_table(FILE* fplog, const char* fn, int angle)
1393 {
1394     int           i;
1395     bondedtable_t tab;
1396     int           stride = 4;
1397
1398     t_tabledata td = read_tables(fplog, fn, 1, angle)[0];
1399     if (angle > 0)
1400     {
1401         /* Convert the table from degrees to radians */
1402         for (i = 0; i < td.nx; i++)
1403         {
1404             td.x[i] *= gmx::c_deg2Rad;
1405             td.f[i] *= gmx::c_rad2Deg;
1406         }
1407         td.tabscale *= gmx::c_rad2Deg;
1408     }
1409     tab.n     = td.nx;
1410     tab.scale = td.tabscale;
1411     tab.data.resize(tab.n * stride);
1412     copy2table(tab.n, 0, stride, td.x, td.v, td.f, 1.0, tab.data);
1413
1414     return tab;
1415 }
1416
1417 std::unique_ptr<t_forcetable>
1418 makeDispersionCorrectionTable(FILE* fp, const interaction_const_t* ic, real rtab, const char* tabfn)
1419 {
1420     GMX_RELEASE_ASSERT(ic->vdwtype != VanDerWaalsType::User || tabfn,
1421                        "With VdW user tables we need a table file name");
1422
1423     std::unique_ptr<t_forcetable> fullTable = make_tables(fp, ic, tabfn, rtab, 0);
1424     /* Copy the contents of the table to one that has just dispersion
1425      * and repulsion, to improve cache performance. We want the table
1426      * data to be aligned to 32-byte boundaries.
1427      */
1428     std::unique_ptr<t_forcetable> dispersionCorrectionTable =
1429             std::make_unique<t_forcetable>(GMX_TABLE_INTERACTION_VDWREP_VDWDISP, fullTable->format);
1430     dispersionCorrectionTable->r             = fullTable->r;
1431     dispersionCorrectionTable->n             = fullTable->n;
1432     dispersionCorrectionTable->scale         = fullTable->scale;
1433     dispersionCorrectionTable->formatsize    = fullTable->formatsize;
1434     dispersionCorrectionTable->ninteractions = 2;
1435     dispersionCorrectionTable->stride =
1436             dispersionCorrectionTable->formatsize * dispersionCorrectionTable->ninteractions;
1437     dispersionCorrectionTable->data.resize(dispersionCorrectionTable->stride
1438                                            * (dispersionCorrectionTable->n + 1));
1439
1440     for (int i = 0; i <= fullTable->n; i++)
1441     {
1442         for (int j = 0; j < 8; j++)
1443         {
1444             dispersionCorrectionTable->data[8 * i + j] = fullTable->data[12 * i + 4 + j];
1445         }
1446     }
1447
1448     return dispersionCorrectionTable;
1449 }
1450
1451 t_forcetable::t_forcetable(enum gmx_table_interaction interaction, enum gmx_table_format format) :
1452     interaction(interaction),
1453     format(format),
1454     r(0),
1455     n(0),
1456     scale(0),
1457     formatsize(0),
1458     ninteractions(0),
1459     stride(0)
1460 {
1461 }
1462
1463 t_forcetable::~t_forcetable() = default;