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