src/gromacs/tables/forcetable.cpp

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