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