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