src/gromacs/mdlib/tables.c

   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, 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 <math.h>
  40 #include "gromacs/math/utilities.h"
  41 #include "gromacs/legacyheaders/typedefs.h"
  42 #include "gromacs/legacyheaders/names.h"
  43 #include "gromacs/utility/smalloc.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "gromacs/utility/futil.h"
  46 #include "gromacs/fileio/xvgr.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/legacyheaders/network.h"
  49 #include "gromacs/math/units.h"
  50 #include "gromacs/legacyheaders/force.h"
  51 #include "gromacs/fileio/gmxfio.h"
  52 #include "gromacs/legacyheaders/macros.h"
  53 #include "gromacs/legacyheaders/tables.h"
  54
  55 /* All the possible (implemented) table functions */
  56 enum {
  57     etabLJ6,
  58     etabLJ12,
  59     etabLJ6Shift,
  60     etabLJ12Shift,
  61     etabShift,
  62     etabRF,
  63     etabRF_ZERO,
  64     etabCOUL,
  65     etabEwald,
  66     etabEwaldSwitch,
  67     etabEwaldUser,
  68     etabEwaldUserSwitch,
  69     etabLJ6Ewald,
  70     etabLJ6Switch,
  71     etabLJ12Switch,
  72     etabCOULSwitch,
  73     etabLJ6Encad,
  74     etabLJ12Encad,
  75     etabCOULEncad,
  76     etabEXPMIN,
  77     etabUSER,
  78     etabNR
  79 };
  80
  81 /** Evaluates to true if the table type contains user data. */
  82 #define ETAB_USER(e)  ((e) == etabUSER || \
  83                        (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
  84
  85 typedef struct {
  86     const char *name;
  87     gmx_bool    bCoulomb;
  88 } t_tab_props;
  89
  90 /* This structure holds name and a flag that tells whether
  91    this is a Coulomb type funtion */
  92 static const t_tab_props tprops[etabNR] = {
  93     { "LJ6",  FALSE },
  94     { "LJ12", FALSE },
  95     { "LJ6Shift", FALSE },
  96     { "LJ12Shift", FALSE },
  97     { "Shift", TRUE },
  98     { "RF", TRUE },
  99     { "RF-zero", TRUE },
 100     { "COUL", TRUE },
 101     { "Ewald", TRUE },
 102     { "Ewald-Switch", TRUE },
 103     { "Ewald-User", TRUE },
 104     { "Ewald-User-Switch", TRUE },
 105     { "LJ6Ewald", FALSE },
 106     { "LJ6Switch", FALSE },
 107     { "LJ12Switch", FALSE },
 108     { "COULSwitch", TRUE },
 109     { "LJ6-Encad shift", FALSE },
 110     { "LJ12-Encad shift", FALSE },
 111     { "COUL-Encad shift",  TRUE },
 112     { "EXPMIN", FALSE },
 113     { "USER", FALSE },
 114 };
 115
 116 /* Index in the table that says which function to use */
 117 enum {
 118     etiCOUL, etiLJ6, etiLJ12, etiNR
 119 };
 120
 121 typedef struct {
 122     int     nx, nx0;
 123     double  tabscale;
 124     double *x, *v, *f;
 125 } t_tabledata;
 126
 127 #define pow2(x) ((x)*(x))
 128 #define pow3(x) ((x)*(x)*(x))
 129 #define pow4(x) ((x)*(x)*(x)*(x))
 130 #define pow5(x) ((x)*(x)*(x)*(x)*(x))
 131
 132 double v_q_ewald_lr(double beta, double r)
 133 {
 134     if (r == 0)
 135     {
 136         return beta*2/sqrt(M_PI);
 137     }
 138     else
 139     {
 140         return gmx_erfd(beta*r)/r;
 141     }
 142 }
 143
 144 double v_lj_ewald_lr(double beta, double r)
 145 {
 146     double br, br2, br4, r6, factor;
 147     if (r == 0)
 148     {
 149         return pow(beta, 6)/6;
 150     }
 151     else
 152     {
 153         br     = beta*r;
 154         br2    = br*br;
 155         br4    = br2*br2;
 156         r6     = pow(r, 6.0);
 157         factor = (1.0 - exp(-br2)*(1 + br2 + 0.5*br4))/r6;
 158         return factor;
 159     }
 160 }
 161
 162 void table_spline3_fill_ewald_lr(real                                 *table_f,
 163                                  real                                 *table_v,
 164                                  real                                 *table_fdv0,
 165                                  int                                   ntab,
 166                                  double                                dx,
 167                                  real                                  beta,
 168                                  real_space_grid_contribution_computer v_lr)
 169 {
 170     real     tab_max;
 171     int      i, i_inrange;
 172     double   dc, dc_new;
 173     gmx_bool bOutOfRange;
 174     double   v_r0, v_r1, v_inrange, vi, a0, a1, a2dx;
 175     double   x_r0;
 176
 177     /* This function is called using either v_ewald_lr or v_lj_ewald_lr as a function argument
 178      * depending on wether we should create electrostatic or Lennard-Jones Ewald tables.
 179      */
 180
 181     if (ntab < 2)
 182     {
 183         gmx_fatal(FARGS, "Can not make a spline table with less than 2 points");
 184     }
 185
 186     /* We need some margin to be able to divide table values by r
 187      * in the kernel and also to do the integration arithmetics
 188      * without going out of range. Furthemore, we divide by dx below.
 189      */
 190     tab_max = GMX_REAL_MAX*0.0001;
 191
 192     /* This function produces a table with:
 193      * maximum energy error: V'''/(6*12*sqrt(3))*dx^3
 194      * maximum force error:  V'''/(6*4)*dx^2
 195      * The rms force error is the max error times 1/sqrt(5)=0.45.
 196      */
 197
 198     bOutOfRange = FALSE;
 199     i_inrange   = ntab;
 200     v_inrange   = 0;
 201     dc          = 0;
 202     for (i = ntab-1; i >= 0; i--)
 203     {
 204         x_r0 = i*dx;
 205
 206         v_r0 = (*v_lr)(beta, x_r0);
 207
 208         if (!bOutOfRange)
 209         {
 210             i_inrange = i;
 211             v_inrange = v_r0;
 212
 213             vi = v_r0;
 214         }
 215         else
 216         {
 217             /* Linear continuation for the last point in range */
 218             vi = v_inrange - dc*(i - i_inrange)*dx;
 219         }
 220
 221         if (table_v != NULL)
 222         {
 223             table_v[i] = vi;
 224         }
 225
 226         if (i == 0)
 227         {
 228             continue;
 229         }
 230
 231         /* Get the potential at table point i-1 */
 232         v_r1 = (*v_lr)(beta, (i-1)*dx);
 233
 234         if (v_r1 != v_r1 || v_r1 < -tab_max || v_r1 > tab_max)
 235         {
 236             bOutOfRange = TRUE;
 237         }
 238
 239         if (!bOutOfRange)
 240         {
 241             /* Calculate the average second derivative times dx over interval i-1 to i.
 242              * Using the function values at the end points and in the middle.
 243              */
 244             a2dx = (v_r0 + v_r1 - 2*(*v_lr)(beta, x_r0-0.5*dx))/(0.25*dx);
 245             /* Set the derivative of the spline to match the difference in potential
 246              * over the interval plus the average effect of the quadratic term.
 247              * This is the essential step for minimizing the error in the force.
 248              */
 249             dc = (v_r0 - v_r1)/dx + 0.5*a2dx;
 250         }
 251
 252         if (i == ntab - 1)
 253         {
 254             /* Fill the table with the force, minus the derivative of the spline */
 255             table_f[i] = -dc;
 256         }
 257         else
 258         {
 259             /* tab[i] will contain the average of the splines over the two intervals */
 260             table_f[i] += -0.5*dc;
 261         }
 262
 263         if (!bOutOfRange)
 264         {
 265             /* Make spline s(x) = a0 + a1*(x - xr) + 0.5*a2*(x - xr)^2
 266              * matching the potential at the two end points
 267              * and the derivative dc at the end point xr.
 268              */
 269             a0   = v_r0;
 270             a1   = dc;
 271             a2dx = (a1*dx + v_r1 - a0)*2/dx;
 272
 273             /* Set dc to the derivative at the next point */
 274             dc_new = a1 - a2dx;
 275
 276             if (dc_new != dc_new || dc_new < -tab_max || dc_new > tab_max)
 277             {
 278                 bOutOfRange = TRUE;
 279             }
 280             else
 281             {
 282                 dc = dc_new;
 283             }
 284         }
 285
 286         table_f[(i-1)] = -0.5*dc;
 287     }
 288     /* Currently the last value only contains half the force: double it */
 289     table_f[0] *= 2;
 290
 291     if (table_v != NULL && table_fdv0 != NULL)
 292     {
 293         /* Copy to FDV0 table too. Allocation occurs in forcerec.c,
 294          * init_ewald_f_table().
 295          */
 296         for (i = 0; i < ntab-1; i++)
 297         {
 298             table_fdv0[4*i]     = table_f[i];
 299             table_fdv0[4*i+1]   = table_f[i+1]-table_f[i];
 300             table_fdv0[4*i+2]   = table_v[i];
 301             table_fdv0[4*i+3]   = 0.0;
 302         }
 303         table_fdv0[4*(ntab-1)]    = table_f[(ntab-1)];
 304         table_fdv0[4*(ntab-1)+1]  = -table_f[(ntab-1)];
 305         table_fdv0[4*(ntab-1)+2]  = table_v[(ntab-1)];
 306         table_fdv0[4*(ntab-1)+3]  = 0.0;
 307     }
 308 }
 309
 310 /* Returns the spacing for a function using the maximum of
 311  * the third derivative, x_scale (unit 1/length)
 312  * and function tolerance.
 313  */
 314 static double spline3_table_scale(double third_deriv_max,
 315                                   double x_scale,
 316                                   double func_tol)
 317 {
 318     double deriv_tol;
 319     double sc_deriv, sc_func;
 320
 321     /* Force tolerance: single precision accuracy */
 322     deriv_tol = GMX_FLOAT_EPS;
 323     sc_deriv  = sqrt(third_deriv_max/(6*4*deriv_tol*x_scale))*x_scale;
 324
 325     /* Don't try to be more accurate on energy than the precision */
 326     func_tol  = max(func_tol, GMX_REAL_EPS);
 327     sc_func   = pow(third_deriv_max/(6*12*sqrt(3)*func_tol), 1.0/3.0)*x_scale;
 328
 329     return max(sc_deriv, sc_func);
 330 }
 331
 332 /* The scale (1/spacing) for third order spline interpolation
 333  * of the Ewald mesh contribution which needs to be subtracted
 334  * from the non-bonded interactions.
 335  * Since there is currently only one spacing for Coulomb and LJ,
 336  * the finest spacing is used if both Ewald types are present.
 337  *
 338  * Note that we could also implement completely separate tables
 339  * for Coulomb and LJ Ewald, each with their own spacing.
 340  * The current setup with the same spacing can provide slightly
 341  * faster kernels with both Coulomb and LJ Ewald, especially
 342  * when interleaving both tables (currently not implemented).
 343  */
 344 real ewald_spline3_table_scale(const interaction_const_t *ic)
 345 {
 346     real sc;
 347
 348     sc = 0;
 349
 350     if (ic->eeltype == eelEWALD || EEL_PME(ic->eeltype))
 351     {
 352         double erf_x_d3 = 1.0522; /* max of (erf(x)/x)''' */
 353         double etol;
 354         real   sc_q;
 355
 356         /* Energy tolerance: 0.1 times the cut-off jump */
 357         etol  = 0.1*gmx_erfc(ic->ewaldcoeff_q*ic->rcoulomb);
 358
 359         sc_q  = spline3_table_scale(erf_x_d3, ic->ewaldcoeff_q, etol);
 360
 361         if (debug)
 362         {
 363             fprintf(debug, "Ewald Coulomb quadratic spline table spacing: %f 1/nm\n", 1/sc_q);
 364         }
 365
 366         sc    = max(sc, sc_q);
 367     }
 368
 369     if (EVDW_PME(ic->vdwtype))
 370     {
 371         double func_d3 = 0.42888; /* max of (x^-6 (1 - exp(-x^2)(1+x^2+x^4/2)))''' */
 372         double xrc2, etol;
 373         real   sc_lj;
 374
 375         /* Energy tolerance: 0.1 times the cut-off jump */
 376         xrc2  = sqr(ic->ewaldcoeff_lj*ic->rvdw);
 377         etol  = 0.1*exp(-xrc2)*(1 + xrc2 + xrc2*xrc2/2.0);
 378
 379         sc_lj = spline3_table_scale(func_d3, ic->ewaldcoeff_lj, etol);
 380
 381         if (debug)
 382         {
 383             fprintf(debug, "Ewald LJ quadratic spline table spacing: %f 1/nm\n", 1/sc_lj);
 384         }
 385
 386         sc = max(sc, sc_lj);
 387     }
 388
 389     return sc;
 390 }
 391
 392 /* Calculate the potential and force for an r value
 393  * in exactly the same way it is done in the inner loop.
 394  * VFtab is a pointer to the table data, offset is
 395  * the point where we should begin and stride is
 396  * 4 if we have a buckingham table, 3 otherwise.
 397  * If you want to evaluate table no N, set offset to 4*N.
 398  *
 399  * We use normal precision here, since that is what we
 400  * will use in the inner loops.
 401  */
 402 static void evaluate_table(real VFtab[], int offset, int stride,
 403                            real tabscale, real r, real *y, real *yp)
 404 {
 405     int  n;
 406     real rt, eps, eps2;
 407     real Y, F, Geps, Heps2, Fp;
 408
 409     rt       =  r*tabscale;
 410     n        =  (int)rt;
 411     eps      =  rt - n;
 412     eps2     =  eps*eps;
 413     n        =  offset+stride*n;
 414     Y        =  VFtab[n];
 415     F        =  VFtab[n+1];
 416     Geps     =  eps*VFtab[n+2];
 417     Heps2    =  eps2*VFtab[n+3];
 418     Fp       =  F+Geps+Heps2;
 419     *y       =  Y+eps*Fp;
 420     *yp      =  (Fp+Geps+2.0*Heps2)*tabscale;
 421 }
 422
 423 static void copy2table(int n, int offset, int stride,
 424                        double x[], double Vtab[], double Ftab[], real scalefactor,
 425                        real dest[])
 426 {
 427 /* Use double prec. for the intermediary variables
 428  * and temporary x/vtab/vtab2 data to avoid unnecessary
 429  * loss of precision.
 430  */
 431     int    i, nn0;
 432     double F, G, H, h;
 433
 434     h = 0;
 435     for (i = 0; (i < n); i++)
 436     {
 437         if (i < n-1)
 438         {
 439             h   = x[i+1] - x[i];
 440             F   = -Ftab[i]*h;
 441             G   =  3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h;
 442             H   = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] +   Ftab[i])*h;
 443         }
 444         else
 445         {
 446             /* Fill the last entry with a linear potential,
 447              * this is mainly for rounding issues with angle and dihedral potentials.
 448              */
 449             F   = -Ftab[i]*h;
 450             G   = 0;
 451             H   = 0;
 452         }
 453         nn0         = offset + i*stride;
 454         dest[nn0]   = scalefactor*Vtab[i];
 455         dest[nn0+1] = scalefactor*F;
 456         dest[nn0+2] = scalefactor*G;
 457         dest[nn0+3] = scalefactor*H;
 458     }
 459 }
 460
 461 static void init_table(int n, int nx0,
 462                        double tabscale, t_tabledata *td, gmx_bool bAlloc)
 463 {
 464     int i;
 465
 466     td->nx       = n;
 467     td->nx0      = nx0;
 468     td->tabscale = tabscale;
 469     if (bAlloc)
 470     {
 471         snew(td->x, td->nx);
 472         snew(td->v, td->nx);
 473         snew(td->f, td->nx);
 474     }
 475     for (i = 0; (i < td->nx); i++)
 476     {
 477         td->x[i] = i/tabscale;
 478     }
 479 }
 480
 481 static void spline_forces(int nx, double h, double v[], gmx_bool bS3, gmx_bool bE3,
 482                           double f[])
 483 {
 484     int    start, end, i;
 485     double v3, b_s, b_e, b;
 486     double beta, *gamma;
 487
 488     /* Formulas can be found in:
 489      * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
 490      */
 491
 492     if (nx < 4 && (bS3 || bE3))
 493     {
 494         gmx_fatal(FARGS, "Can not generate splines with third derivative boundary conditions with less than 4 (%d) points", nx);
 495     }
 496
 497     /* To make life easy we initially set the spacing to 1
 498      * and correct for this at the end.
 499      */
 500     beta = 2;
 501     if (bS3)
 502     {
 503         /* Fit V''' at the start */
 504         v3  = v[3] - 3*v[2] + 3*v[1] - v[0];
 505         if (debug)
 506         {
 507             fprintf(debug, "The left third derivative is %g\n", v3/(h*h*h));
 508         }
 509         b_s   = 2*(v[1] - v[0]) + v3/6;
 510         start = 0;
 511
 512         if (FALSE)
 513         {
 514             /* Fit V'' at the start */
 515             real v2;
 516
 517             v2  = -v[3] + 4*v[2] - 5*v[1] + 2*v[0];
 518             /* v2  = v[2] - 2*v[1] + v[0]; */
 519             if (debug)
 520             {
 521                 fprintf(debug, "The left second derivative is %g\n", v2/(h*h));
 522             }
 523             b_s   = 3*(v[1] - v[0]) - v2/2;
 524             start = 0;
 525         }
 526     }
 527     else
 528     {
 529         b_s   = 3*(v[2] - v[0]) + f[0]*h;
 530         start = 1;
 531     }
 532     if (bE3)
 533     {
 534         /* Fit V''' at the end */
 535         v3  = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4];
 536         if (debug)
 537         {
 538             fprintf(debug, "The right third derivative is %g\n", v3/(h*h*h));
 539         }
 540         b_e = 2*(v[nx-1] - v[nx-2]) + v3/6;
 541         end = nx;
 542     }
 543     else
 544     {
 545         /* V'=0 at the end */
 546         b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h;
 547         end = nx - 1;
 548     }
 549
 550     snew(gamma, nx);
 551     beta = (bS3 ? 1 : 4);
 552
 553     /* For V'' fitting */
 554     /* beta = (bS3 ? 2 : 4); */
 555
 556     f[start] = b_s/beta;
 557     for (i = start+1; i < end; i++)
 558     {
 559         gamma[i] = 1/beta;
 560         beta     = 4 - gamma[i];
 561         b        =  3*(v[i+1] - v[i-1]);
 562         f[i]     = (b - f[i-1])/beta;
 563     }
 564     gamma[end-1] = 1/beta;
 565     beta         = (bE3 ? 1 : 4) - gamma[end-1];
 566     f[end-1]     = (b_e - f[end-2])/beta;
 567
 568     for (i = end-2; i >= start; i--)
 569     {
 570         f[i] -= gamma[i+1]*f[i+1];
 571     }
 572     sfree(gamma);
 573
 574     /* Correct for the minus sign and the spacing */
 575     for (i = start; i < end; i++)
 576     {
 577         f[i] = -f[i]/h;
 578     }
 579 }
 580
 581 static void set_forces(FILE *fp, int angle,
 582                        int nx, double h, double v[], double f[],
 583                        int table)
 584 {
 585     int start, end;
 586
 587     if (angle == 2)
 588     {
 589         gmx_fatal(FARGS,
 590                   "Force generation for dihedral tables is not (yet) implemented");
 591     }
 592
 593     start = 0;
 594     while (v[start] == 0)
 595     {
 596         start++;
 597     }
 598
 599     end = nx;
 600     while (v[end-1] == 0)
 601     {
 602         end--;
 603     }
 604     if (end > nx - 2)
 605     {
 606         end = nx;
 607     }
 608     else
 609     {
 610         end++;
 611     }
 612
 613     if (fp)
 614     {
 615         fprintf(fp, "Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
 616                 table+1, start*h, end == nx ? "V'''" : "V'=0", (end-1)*h);
 617     }
 618     spline_forces(end-start, h, v+start, TRUE, end == nx, f+start);
 619 }
 620
 621 static void read_tables(FILE *fp, const char *fn,
 622                         int ntab, int angle, t_tabledata td[])
 623 {
 624     char    *libfn;
 625     char     buf[STRLEN];
 626     double **yy = NULL, start, end, dx0, dx1, ssd, vm, vp, f, numf;
 627     int      k, i, nx, nx0 = 0, ny, nny, ns;
 628     gmx_bool bAllZero, bZeroV, bZeroF;
 629     double   tabscale;
 630
 631     nny   = 2*ntab+1;
 632     libfn = gmxlibfn(fn);
 633     nx    = read_xvg(libfn, &yy, &ny);
 634     if (ny != nny)
 635     {
 636         gmx_fatal(FARGS, "Trying to read file %s, but nr columns = %d, should be %d",
 637                   libfn, ny, nny);
 638     }
 639     if (angle == 0)
 640     {
 641         if (yy[0][0] != 0.0)
 642         {
 643             gmx_fatal(FARGS,
 644                       "The first distance in file %s is %f nm instead of %f nm",
 645                       libfn, yy[0][0], 0.0);
 646         }
 647     }
 648     else
 649     {
 650         if (angle == 1)
 651         {
 652             start = 0.0;
 653         }
 654         else
 655         {
 656             start = -180.0;
 657         }
 658         end = 180.0;
 659         if (yy[0][0] != start || yy[0][nx-1] != end)
 660         {
 661             gmx_fatal(FARGS, "The angles in file %s should go from %f to %f instead of %f to %f\n",
 662                       libfn, start, end, yy[0][0], yy[0][nx-1]);
 663         }
 664     }
 665
 666     tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
 667
 668     if (fp)
 669     {
 670         fprintf(fp, "Read user tables from %s with %d data points.\n", libfn, nx);
 671         if (angle == 0)
 672         {
 673             fprintf(fp, "Tabscale = %g points/nm\n", tabscale);
 674         }
 675     }
 676
 677     bAllZero = TRUE;
 678     for (k = 0; k < ntab; k++)
 679     {
 680         bZeroV = TRUE;
 681         bZeroF = TRUE;
 682         for (i = 0; (i < nx); i++)
 683         {
 684             if (i >= 2)
 685             {
 686                 dx0 = yy[0][i-1] - yy[0][i-2];
 687                 dx1 = yy[0][i]   - yy[0][i-1];
 688                 /* Check for 1% deviation in spacing */
 689                 if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1)))
 690                 {
 691                     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]);
 692                 }
 693             }
 694             if (yy[1+k*2][i] != 0)
 695             {
 696                 bZeroV = FALSE;
 697                 if (bAllZero)
 698                 {
 699                     bAllZero = FALSE;
 700                     nx0      = i;
 701                 }
 702                 if (yy[1+k*2][i] >  0.01*GMX_REAL_MAX ||
 703                     yy[1+k*2][i] < -0.01*GMX_REAL_MAX)
 704                 {
 705                     gmx_fatal(FARGS, "Out of range potential value %g in file '%s'",
 706                               yy[1+k*2][i], fn);
 707                 }
 708             }
 709             if (yy[1+k*2+1][i] != 0)
 710             {
 711                 bZeroF = FALSE;
 712                 if (bAllZero)
 713                 {
 714                     bAllZero = FALSE;
 715                     nx0      = i;
 716                 }
 717                 if (yy[1+k*2+1][i] >  0.01*GMX_REAL_MAX ||
 718                     yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX)
 719                 {
 720                     gmx_fatal(FARGS, "Out of range force value %g in file '%s'",
 721                               yy[1+k*2+1][i], fn);
 722                 }
 723             }
 724         }
 725
 726         if (!bZeroV && bZeroF)
 727         {
 728             set_forces(fp, angle, nx, 1/tabscale, yy[1+k*2], yy[1+k*2+1], k);
 729         }
 730         else
 731         {
 732             /* Check if the second column is close to minus the numerical
 733              * derivative of the first column.
 734              */
 735             ssd = 0;
 736             ns  = 0;
 737             for (i = 1; (i < nx-1); i++)
 738             {
 739                 vm = yy[1+2*k][i-1];
 740                 vp = yy[1+2*k][i+1];
 741                 f  = yy[1+2*k+1][i];
 742                 if (vm != 0 && vp != 0 && f != 0)
 743                 {
 744                     /* Take the centered difference */
 745                     numf = -(vp - vm)*0.5*tabscale;
 746                     ssd += fabs(2*(f - numf)/(f + numf));
 747                     ns++;
 748                 }
 749             }
 750             if (ns > 0)
 751             {
 752                 ssd /= ns;
 753                 sprintf(buf, "For the %d non-zero entries for table %d in %s the forces deviate on average %d%% from minus the numerical derivative of the potential\n", ns, k, libfn, (int)(100*ssd+0.5));
 754                 if (debug)
 755                 {
 756                     fprintf(debug, "%s", buf);
 757                 }
 758                 if (ssd > 0.2)
 759                 {
 760                     if (fp)
 761                     {
 762                         fprintf(fp, "\nWARNING: %s\n", buf);
 763                     }
 764                     fprintf(stderr, "\nWARNING: %s\n", buf);
 765                 }
 766             }
 767         }
 768     }
 769     if (bAllZero && fp)
 770     {
 771         fprintf(fp, "\nNOTE: All elements in table %s are zero\n\n", libfn);
 772     }
 773
 774     for (k = 0; (k < ntab); k++)
 775     {
 776         init_table(nx, nx0, tabscale, &(td[k]), TRUE);
 777         for (i = 0; (i < nx); i++)
 778         {
 779             td[k].x[i] = yy[0][i];
 780             td[k].v[i] = yy[2*k+1][i];
 781             td[k].f[i] = yy[2*k+2][i];
 782         }
 783     }
 784     for (i = 0; (i < ny); i++)
 785     {
 786         sfree(yy[i]);
 787     }
 788     sfree(yy);
 789     sfree(libfn);
 790 }
 791
 792 static void done_tabledata(t_tabledata *td)
 793 {
 794     int i;
 795
 796     if (!td)
 797     {
 798         return;
 799     }
 800
 801     sfree(td->x);
 802     sfree(td->v);
 803     sfree(td->f);
 804 }
 805
 806 static void fill_table(t_tabledata *td, int tp, const t_forcerec *fr,
 807                        gmx_bool b14only)
 808 {
 809     /* Fill the table according to the formulas in the manual.
 810      * In principle, we only need the potential and the second
 811      * derivative, but then we would have to do lots of calculations
 812      * in the inner loop. By precalculating some terms (see manual)
 813      * we get better eventual performance, despite a larger table.
 814      *
 815      * Since some of these higher-order terms are very small,
 816      * we always use double precision to calculate them here, in order
 817      * to avoid unnecessary loss of precision.
 818      */
 819 #ifdef DEBUG_SWITCH
 820     FILE    *fp;
 821 #endif
 822     int      i;
 823     double   reppow, p;
 824     double   r1, rc, r12, r13;
 825     double   r, r2, r6, rc2, rc6, rc12;
 826     double   expr, Vtab, Ftab;
 827     /* Parameters for David's function */
 828     double   A = 0, B = 0, C = 0, A_3 = 0, B_4 = 0;
 829     /* Parameters for the switching function */
 830     double   ksw, swi, swi1;
 831     /* Temporary parameters */
 832     gmx_bool bPotentialSwitch, bForceSwitch, bPotentialShift;
 833     double   ewc   = fr->ewaldcoeff_q;
 834     double   ewclj = fr->ewaldcoeff_lj;
 835     double   Vcut  = 0;
 836
 837     if (b14only)
 838     {
 839         bPotentialSwitch = FALSE;
 840         bForceSwitch     = FALSE;
 841         bPotentialShift  = FALSE;
 842     }
 843     else
 844     {
 845         bPotentialSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
 846                             (tp == etabCOULSwitch) ||
 847                             (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch) ||
 848                             (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSWITCH)) ||
 849                             (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSWITCH)));
 850         bForceSwitch  = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
 851                          (tp == etabShift) ||
 852                          (tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodFORCESWITCH)) ||
 853                          (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodFORCESWITCH)));
 854         bPotentialShift = ((tprops[tp].bCoulomb && (fr->coulomb_modifier == eintmodPOTSHIFT)) ||
 855                            (!tprops[tp].bCoulomb && (fr->vdw_modifier == eintmodPOTSHIFT)));
 856     }
 857
 858     reppow = fr->reppow;
 859
 860     if (tprops[tp].bCoulomb)
 861     {
 862         r1 = fr->rcoulomb_switch;
 863         rc = fr->rcoulomb;
 864     }
 865     else
 866     {
 867         r1 = fr->rvdw_switch;
 868         rc = fr->rvdw;
 869     }
 870     if (bPotentialSwitch)
 871     {
 872         ksw  = 1.0/(pow5(rc-r1));
 873     }
 874     else
 875     {
 876         ksw  = 0.0;
 877     }
 878     if (bForceSwitch)
 879     {
 880         if (tp == etabShift)
 881         {
 882             p = 1;
 883         }
 884         else if (tp == etabLJ6Shift)
 885         {
 886             p = 6;
 887         }
 888         else
 889         {
 890             p = reppow;
 891         }
 892
 893         A = p * ((p+1)*r1-(p+4)*rc)/(pow(rc, p+2)*pow2(rc-r1));
 894         B = -p * ((p+1)*r1-(p+3)*rc)/(pow(rc, p+2)*pow3(rc-r1));
 895         C = 1.0/pow(rc, p)-A/3.0*pow3(rc-r1)-B/4.0*pow4(rc-r1);
 896         if (tp == etabLJ6Shift)
 897         {
 898             A = -A;
 899             B = -B;
 900             C = -C;
 901         }
 902         A_3 = A/3.0;
 903         B_4 = B/4.0;
 904     }
 905     if (debug)
 906     {
 907         fprintf(debug, "Setting up tables\n"); fflush(debug);
 908     }
 909
 910 #ifdef DEBUG_SWITCH
 911     fp = xvgropen("switch.xvg", "switch", "r", "s");
 912 #endif
 913
 914     if (bPotentialShift)
 915     {
 916         rc2   = rc*rc;
 917         rc6   = 1.0/(rc2*rc2*rc2);
 918         if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
 919         {
 920             rc12 = rc6*rc6;
 921         }
 922         else
 923         {
 924             rc12 = pow(rc, -reppow);
 925         }
 926
 927         switch (tp)
 928         {
 929             case etabLJ6:
 930                 /* Dispersion */
 931                 Vcut = -rc6;
 932                 break;
 933             case etabLJ6Ewald:
 934                 Vcut  = -rc6*exp(-ewclj*ewclj*rc2)*(1 + ewclj*ewclj*rc2 + pow4(ewclj)*rc2*rc2/2);
 935                 break;
 936             case etabLJ12:
 937                 /* Repulsion */
 938                 Vcut  = rc12;
 939                 break;
 940             case etabCOUL:
 941                 Vcut  = 1.0/rc;
 942                 break;
 943             case etabEwald:
 944             case etabEwaldSwitch:
 945                 Vtab  = gmx_erfc(ewc*rc)/rc;
 946                 break;
 947             case etabEwaldUser:
 948                 /* Only calculate minus the reciprocal space contribution */
 949                 Vtab  = -gmx_erf(ewc*rc)/rc;
 950                 break;
 951             case etabRF:
 952             case etabRF_ZERO:
 953                 /* No need for preventing the usage of modifiers with RF */
 954                 Vcut  = 0.0;
 955                 break;
 956             case etabEXPMIN:
 957                 Vcut  = exp(-rc);
 958                 break;
 959             default:
 960                 gmx_fatal(FARGS, "Cannot apply new potential-shift modifier to interaction type '%s' yet. (%s,%d)",
 961                           tprops[tp].name, __FILE__, __LINE__);
 962         }
 963     }
 964
 965     for (i = td->nx0; (i < td->nx); i++)
 966     {
 967         r     = td->x[i];
 968         r2    = r*r;
 969         r6    = 1.0/(r2*r2*r2);
 970         if (gmx_within_tol(reppow, 12.0, 10*GMX_DOUBLE_EPS))
 971         {
 972             r12 = r6*r6;
 973         }
 974         else
 975         {
 976             r12 = pow(r, -reppow);
 977         }
 978         Vtab  = 0.0;
 979         Ftab  = 0.0;
 980         if (bPotentialSwitch)
 981         {
 982             /* swi is function, swi1 1st derivative and swi2 2nd derivative */
 983             /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
 984              * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
 985              * r1 and rc.
 986              * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
 987              */
 988             if (r <= r1)
 989             {
 990                 swi  = 1.0;
 991                 swi1 = 0.0;
 992             }
 993             else if (r >= rc)
 994             {
 995                 swi  = 0.0;
 996                 swi1 = 0.0;
 997             }
 998             else
 999             {
1000                 swi      = 1 - 10*pow3(r-r1)*ksw*pow2(rc-r1)
1001                     + 15*pow4(r-r1)*ksw*(rc-r1) - 6*pow5(r-r1)*ksw;
1002                 swi1     = -30*pow2(r-r1)*ksw*pow2(rc-r1)
1003                     + 60*pow3(r-r1)*ksw*(rc-r1) - 30*pow4(r-r1)*ksw;
1004             }
1005         }
1006         else /* not really needed, but avoids compiler warnings... */
1007         {
1008             swi  = 1.0;
1009             swi1 = 0.0;
1010         }
1011 #ifdef DEBUG_SWITCH
1012         fprintf(fp, "%10g  %10g  %10g  %10g\n", r, swi, swi1, swi2);
1013 #endif
1014
1015         rc6 = rc*rc*rc;
1016         rc6 = 1.0/(rc6*rc6);
1017
1018         switch (tp)
1019         {
1020             case etabLJ6:
1021                 /* Dispersion */
1022                 Vtab = -r6;
1023                 Ftab = 6.0*Vtab/r;
1024                 break;
1025             case etabLJ6Switch:
1026             case etabLJ6Shift:
1027                 /* Dispersion */
1028                 if (r < rc)
1029                 {
1030                     Vtab = -r6;
1031                     Ftab = 6.0*Vtab/r;
1032                     break;
1033                 }
1034                 break;
1035             case etabLJ12:
1036                 /* Repulsion */
1037                 Vtab  = r12;
1038                 Ftab  = reppow*Vtab/r;
1039                 break;
1040             case etabLJ12Switch:
1041             case etabLJ12Shift:
1042                 /* Repulsion */
1043                 if (r < rc)
1044                 {
1045                     Vtab  = r12;
1046                     Ftab  = reppow*Vtab/r;
1047                 }
1048                 break;
1049             case etabLJ6Encad:
1050                 if (r < rc)
1051                 {
1052                     Vtab  = -(r6-6.0*(rc-r)*rc6/rc-rc6);
1053                     Ftab  = -(6.0*r6/r-6.0*rc6/rc);
1054                 }
1055                 else /* r>rc */
1056                 {
1057                     Vtab  = 0;
1058                     Ftab  = 0;
1059                 }
1060                 break;
1061             case etabLJ12Encad:
1062                 if (r < rc)
1063                 {
1064                     Vtab  = -(r6-6.0*(rc-r)*rc6/rc-rc6);
1065                     Ftab  = -(6.0*r6/r-6.0*rc6/rc);
1066                 }
1067                 else /* r>rc */
1068                 {
1069                     Vtab  = 0;
1070                     Ftab  = 0;
1071                 }
1072                 break;
1073             case etabCOUL:
1074                 Vtab  = 1.0/r;
1075                 Ftab  = 1.0/r2;
1076                 break;
1077             case etabCOULSwitch:
1078             case etabShift:
1079                 if (r < rc)
1080                 {
1081                     Vtab  = 1.0/r;
1082                     Ftab  = 1.0/r2;
1083                 }
1084                 break;
1085             case etabEwald:
1086             case etabEwaldSwitch:
1087                 Vtab  = gmx_erfc(ewc*r)/r;
1088                 Ftab  = gmx_erfc(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
1089                 break;
1090             case etabEwaldUser:
1091             case etabEwaldUserSwitch:
1092                 /* Only calculate the negative of the reciprocal space contribution */
1093                 Vtab  = -gmx_erf(ewc*r)/r;
1094                 Ftab  = -gmx_erf(ewc*r)/r2+exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r;
1095                 break;
1096             case etabLJ6Ewald:
1097                 Vtab  = -r6*exp(-ewclj*ewclj*r2)*(1 + ewclj*ewclj*r2 + pow4(ewclj)*r2*r2/2);
1098                 Ftab  = 6.0*Vtab/r - r6*exp(-ewclj*ewclj*r2)*pow5(ewclj)*ewclj*r2*r2*r;
1099                 break;
1100             case etabRF:
1101             case etabRF_ZERO:
1102                 Vtab  = 1.0/r      +   fr->k_rf*r2 - fr->c_rf;
1103                 Ftab  = 1.0/r2     - 2*fr->k_rf*r;
1104                 if (tp == etabRF_ZERO && r >= rc)
1105                 {
1106                     Vtab = 0;
1107                     Ftab = 0;
1108                 }
1109                 break;
1110             case etabEXPMIN:
1111                 expr  = exp(-r);
1112                 Vtab  = expr;
1113                 Ftab  = expr;
1114                 break;
1115             case etabCOULEncad:
1116                 if (r < rc)
1117                 {
1118                     Vtab  = 1.0/r-(rc-r)/(rc*rc)-1.0/rc;
1119                     Ftab  = 1.0/r2-1.0/(rc*rc);
1120                 }
1121                 else /* r>rc */
1122                 {
1123                     Vtab  = 0;
1124                     Ftab  = 0;
1125                 }
1126                 break;
1127             default:
1128                 gmx_fatal(FARGS, "Table type %d not implemented yet. (%s,%d)",
1129                           tp, __FILE__, __LINE__);
1130         }
1131         if (bForceSwitch)
1132         {
1133             /* Normal coulomb with cut-off correction for potential */
1134             if (r < rc)
1135             {
1136                 Vtab -= C;
1137                 /* If in Shifting range add something to it */
1138                 if (r > r1)
1139                 {
1140                     r12    = (r-r1)*(r-r1);
1141                     r13    = (r-r1)*r12;
1142                     Vtab  += -A_3*r13 - B_4*r12*r12;
1143                     Ftab  +=   A*r12 + B*r13;
1144                 }
1145             }
1146             else
1147             {
1148                 /* Make sure interactions are zero outside cutoff with modifiers */
1149                 Vtab = 0;
1150                 Ftab = 0;
1151             }
1152         }
1153         if (bPotentialShift)
1154         {
1155             if (r < rc)
1156             {
1157                 Vtab -= Vcut;
1158             }
1159             else
1160             {
1161                 /* Make sure interactions are zero outside cutoff with modifiers */
1162                 Vtab = 0;
1163                 Ftab = 0;
1164             }
1165         }
1166
1167         if (ETAB_USER(tp))
1168         {
1169             Vtab += td->v[i];
1170             Ftab += td->f[i];
1171         }
1172
1173         if (bPotentialSwitch)
1174         {
1175             if (r >= rc)
1176             {
1177                 /* Make sure interactions are zero outside cutoff with modifiers */
1178                 Vtab = 0;
1179                 Ftab = 0;
1180             }
1181             else if (r > r1)
1182             {
1183                 Ftab = Ftab*swi - Vtab*swi1;
1184                 Vtab = Vtab*swi;
1185             }
1186         }
1187         /* Convert to single precision when we store to mem */
1188         td->v[i]  = Vtab;
1189         td->f[i]  = Ftab;
1190     }
1191
1192     /* Continue the table linearly from nx0 to 0.
1193      * These values are only required for energy minimization with overlap or TPI.
1194      */
1195     for (i = td->nx0-1; i >= 0; i--)
1196     {
1197         td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
1198         td->f[i] = td->f[i+1];
1199     }
1200
1201 #ifdef DEBUG_SWITCH
1202     gmx_fio_fclose(fp);
1203 #endif
1204 }
1205
1206 static void set_table_type(int tabsel[], const t_forcerec *fr, gmx_bool b14only)
1207 {
1208     int eltype, vdwtype;
1209
1210     /* Set the different table indices.
1211      * Coulomb first.
1212      */
1213
1214
1215     if (b14only)
1216     {
1217         switch (fr->eeltype)
1218         {
1219             case eelRF_NEC:
1220                 eltype = eelRF;
1221                 break;
1222             case eelUSER:
1223             case eelPMEUSER:
1224             case eelPMEUSERSWITCH:
1225                 eltype = eelUSER;
1226                 break;
1227             default:
1228                 eltype = eelCUT;
1229         }
1230     }
1231     else
1232     {
1233         eltype = fr->eeltype;
1234     }
1235
1236     switch (eltype)
1237     {
1238         case eelCUT:
1239             tabsel[etiCOUL] = etabCOUL;
1240             break;
1241         case eelPOISSON:
1242             tabsel[etiCOUL] = etabShift;
1243             break;
1244         case eelSHIFT:
1245             if (fr->rcoulomb > fr->rcoulomb_switch)
1246             {
1247                 tabsel[etiCOUL] = etabShift;
1248             }
1249             else
1250             {
1251                 tabsel[etiCOUL] = etabCOUL;
1252             }
1253             break;
1254         case eelEWALD:
1255         case eelPME:
1256         case eelP3M_AD:
1257             tabsel[etiCOUL] = etabEwald;
1258             break;
1259         case eelPMESWITCH:
1260             tabsel[etiCOUL] = etabEwaldSwitch;
1261             break;
1262         case eelPMEUSER:
1263             tabsel[etiCOUL] = etabEwaldUser;
1264             break;
1265         case eelPMEUSERSWITCH:
1266             tabsel[etiCOUL] = etabEwaldUserSwitch;
1267             break;
1268         case eelRF:
1269         case eelGRF:
1270         case eelRF_NEC:
1271             tabsel[etiCOUL] = etabRF;
1272             break;
1273         case eelRF_ZERO:
1274             tabsel[etiCOUL] = etabRF_ZERO;
1275             break;
1276         case eelSWITCH:
1277             tabsel[etiCOUL] = etabCOULSwitch;
1278             break;
1279         case eelUSER:
1280             tabsel[etiCOUL] = etabUSER;
1281             break;
1282         case eelENCADSHIFT:
1283             tabsel[etiCOUL] = etabCOULEncad;
1284             break;
1285         default:
1286             gmx_fatal(FARGS, "Invalid eeltype %d", eltype);
1287     }
1288
1289     /* Van der Waals time */
1290     if (fr->bBHAM && !b14only)
1291     {
1292         tabsel[etiLJ6]  = etabLJ6;
1293         tabsel[etiLJ12] = etabEXPMIN;
1294     }
1295     else
1296     {
1297         if (b14only && fr->vdwtype != evdwUSER)
1298         {
1299             vdwtype = evdwCUT;
1300         }
1301         else
1302         {
1303             vdwtype = fr->vdwtype;
1304         }
1305
1306         switch (vdwtype)
1307         {
1308             case evdwSWITCH:
1309                 tabsel[etiLJ6]  = etabLJ6Switch;
1310                 tabsel[etiLJ12] = etabLJ12Switch;
1311                 break;
1312             case evdwSHIFT:
1313                 tabsel[etiLJ6]  = etabLJ6Shift;
1314                 tabsel[etiLJ12] = etabLJ12Shift;
1315                 break;
1316             case evdwUSER:
1317                 tabsel[etiLJ6]  = etabUSER;
1318                 tabsel[etiLJ12] = etabUSER;
1319                 break;
1320             case evdwCUT:
1321                 tabsel[etiLJ6]  = etabLJ6;
1322                 tabsel[etiLJ12] = etabLJ12;
1323                 break;
1324             case evdwENCADSHIFT:
1325                 tabsel[etiLJ6]  = etabLJ6Encad;
1326                 tabsel[etiLJ12] = etabLJ12Encad;
1327                 break;
1328             case evdwPME:
1329                 tabsel[etiLJ6]  = etabLJ6Ewald;
1330                 tabsel[etiLJ12] = etabLJ12;
1331                 break;
1332             default:
1333                 gmx_fatal(FARGS, "Invalid vdwtype %d in %s line %d", vdwtype,
1334                           __FILE__, __LINE__);
1335         }
1336
1337         if (!b14only && fr->vdw_modifier != eintmodNONE)
1338         {
1339             if (fr->vdw_modifier != eintmodPOTSHIFT &&
1340                 fr->vdwtype != evdwCUT)
1341             {
1342                 gmx_incons("Potential modifiers other than potential-shift are only implemented for LJ cut-off");
1343             }
1344
1345             /* LJ-PME and other (shift-only) modifiers are handled by applying the modifiers
1346              * to the original interaction forms when we fill the table, so we only check cutoffs here.
1347              */
1348             if (fr->vdwtype == evdwCUT)
1349             {
1350                 switch (fr->vdw_modifier)
1351                 {
1352                     case eintmodNONE:
1353                     case eintmodPOTSHIFT:
1354                     case eintmodEXACTCUTOFF:
1355                         /* No modification */
1356                         break;
1357                     case eintmodPOTSWITCH:
1358                         tabsel[etiLJ6]  = etabLJ6Switch;
1359                         tabsel[etiLJ12] = etabLJ12Switch;
1360                         break;
1361                     case eintmodFORCESWITCH:
1362                         tabsel[etiLJ6]  = etabLJ6Shift;
1363                         tabsel[etiLJ12] = etabLJ12Shift;
1364                         break;
1365                     default:
1366                         gmx_incons("Unsupported vdw_modifier");
1367                 }
1368             }
1369         }
1370     }
1371 }
1372
1373 t_forcetable make_tables(FILE *out, const output_env_t oenv,
1374                          const t_forcerec *fr,
1375                          gmx_bool bVerbose, const char *fn,
1376                          real rtab, int flags)
1377 {
1378     const char     *fns[3]   = { "ctab.xvg", "dtab.xvg", "rtab.xvg" };
1379     const char     *fns14[3] = { "ctab14.xvg", "dtab14.xvg", "rtab14.xvg" };
1380     FILE           *fp;
1381     t_tabledata    *td;
1382     gmx_bool        b14only, bReadTab, bGenTab;
1383     real            x0, y0, yp;
1384     int             i, j, k, nx, nx0, tabsel[etiNR];
1385     real            scalefactor;
1386
1387     t_forcetable    table;
1388
1389     b14only = (flags & GMX_MAKETABLES_14ONLY);
1390
1391     if (flags & GMX_MAKETABLES_FORCEUSER)
1392     {
1393         tabsel[etiCOUL] = etabUSER;
1394         tabsel[etiLJ6]  = etabUSER;
1395         tabsel[etiLJ12] = etabUSER;
1396     }
1397     else
1398     {
1399         set_table_type(tabsel, fr, b14only);
1400     }
1401     snew(td, etiNR);
1402     table.r         = rtab;
1403     table.scale     = 0;
1404     table.n         = 0;
1405     table.scale_exp = 0;
1406     nx0             = 10;
1407     nx              = 0;
1408
1409     table.interaction   = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
1410     table.format        = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1411     table.formatsize    = 4;
1412     table.ninteractions = 3;
1413     table.stride        = table.formatsize*table.ninteractions;
1414
1415     /* Check whether we have to read or generate */
1416     bReadTab = FALSE;
1417     bGenTab  = FALSE;
1418     for (i = 0; (i < etiNR); i++)
1419     {
1420         if (ETAB_USER(tabsel[i]))
1421         {
1422             bReadTab = TRUE;
1423         }
1424         if (tabsel[i] != etabUSER)
1425         {
1426             bGenTab  = TRUE;
1427         }
1428     }
1429     if (bReadTab)
1430     {
1431         read_tables(out, fn, etiNR, 0, td);
1432         if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY))
1433         {
1434             rtab      = td[0].x[td[0].nx-1];
1435             table.n   = td[0].nx;
1436             nx        = table.n;
1437         }
1438         else
1439         {
1440             if (td[0].x[td[0].nx-1] < rtab)
1441             {
1442                 gmx_fatal(FARGS, "Tables in file %s not long enough for cut-off:\n"
1443                           "\tshould be at least %f nm\n", fn, rtab);
1444             }
1445             nx        = table.n = (int)(rtab*td[0].tabscale + 0.5);
1446         }
1447         table.scale = td[0].tabscale;
1448         nx0         = td[0].nx0;
1449     }
1450     if (bGenTab)
1451     {
1452         if (!bReadTab)
1453         {
1454 #ifdef GMX_DOUBLE
1455             table.scale = 2000.0;
1456 #else
1457             table.scale = 500.0;
1458 #endif
1459             nx = table.n = rtab*table.scale;
1460         }
1461     }
1462     if (fr->bBHAM)
1463     {
1464         if (fr->bham_b_max != 0)
1465         {
1466             table.scale_exp = table.scale/fr->bham_b_max;
1467         }
1468         else
1469         {
1470             table.scale_exp = table.scale;
1471         }
1472     }
1473
1474     /* Each table type (e.g. coul,lj6,lj12) requires four
1475      * numbers per nx+1 data points. For performance reasons we want
1476      * the table data to be aligned to 16-byte.
1477      */
1478     snew_aligned(table.data, 12*(nx+1)*sizeof(real), 32);
1479
1480     for (k = 0; (k < etiNR); k++)
1481     {
1482         if (tabsel[k] != etabUSER)
1483         {
1484             init_table(nx, nx0,
1485                        (tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
1486                        &(td[k]), !bReadTab);
1487             fill_table(&(td[k]), tabsel[k], fr, b14only);
1488             if (out)
1489             {
1490                 fprintf(out, "%s table with %d data points for %s%s.\n"
1491                         "Tabscale = %g points/nm\n",
1492                         ETAB_USER(tabsel[k]) ? "Modified" : "Generated",
1493                         td[k].nx, b14only ? "1-4 " : "", tprops[tabsel[k]].name,
1494                         td[k].tabscale);
1495             }
1496         }
1497
1498         /* Set scalefactor for c6/c12 tables. This is because we save flops in the non-table kernels
1499          * by including the derivative constants (6.0 or 12.0) in the parameters, since
1500          * we no longer calculate force in most steps. This means the c6/c12 parameters
1501          * have been scaled up, so we need to scale down the table interactions too.
1502          * It comes here since we need to scale user tables too.
1503          */
1504         if (k == etiLJ6)
1505         {
1506             scalefactor = 1.0/6.0;
1507         }
1508         else if (k == etiLJ12 && tabsel[k] != etabEXPMIN)
1509         {
1510             scalefactor = 1.0/12.0;
1511         }
1512         else
1513         {
1514             scalefactor = 1.0;
1515         }
1516
1517         copy2table(table.n, k*4, 12, td[k].x, td[k].v, td[k].f, scalefactor, table.data);
1518
1519         if (bDebugMode() && bVerbose)
1520         {
1521             if (b14only)
1522             {
1523                 fp = xvgropen(fns14[k], fns14[k], "r", "V", oenv);
1524             }
1525             else
1526             {
1527                 fp = xvgropen(fns[k], fns[k], "r", "V", oenv);
1528             }
1529             /* plot the output 5 times denser than the table data */
1530             for (i = 5*((nx0+1)/2); i < 5*table.n; i++)
1531             {
1532                 x0 = i*table.r/(5*(table.n-1));
1533                 evaluate_table(table.data, 4*k, 12, table.scale, x0, &y0, &yp);
1534                 fprintf(fp, "%15.10e  %15.10e  %15.10e\n", x0, y0, yp);
1535             }
1536             gmx_fio_fclose(fp);
1537         }
1538         done_tabledata(&(td[k]));
1539     }
1540     sfree(td);
1541
1542     return table;
1543 }
1544
1545 t_forcetable make_gb_table(const output_env_t oenv,
1546                            const t_forcerec  *fr)
1547 {
1548     const char     *fns[3]   = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" };
1549     const char     *fns14[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" };
1550     FILE           *fp;
1551     t_tabledata    *td;
1552     gmx_bool        bReadTab, bGenTab;
1553     real            x0, y0, yp;
1554     int             i, j, k, nx, nx0, tabsel[etiNR];
1555     double          r, r2, Vtab, Ftab, expterm;
1556
1557     t_forcetable    table;
1558
1559     double          abs_error_r, abs_error_r2;
1560     double          rel_error_r, rel_error_r2;
1561     double          rel_error_r_old = 0, rel_error_r2_old = 0;
1562     double          x0_r_error, x0_r2_error;
1563
1564
1565     /* Only set a Coulomb table for GB */
1566     /*
1567        tabsel[0]=etabGB;
1568        tabsel[1]=-1;
1569        tabsel[2]=-1;
1570      */
1571
1572     /* Set the table dimensions for GB, not really necessary to
1573      * use etiNR (since we only have one table, but ...)
1574      */
1575     snew(td, 1);
1576     table.interaction   = GMX_TABLE_INTERACTION_ELEC;
1577     table.format        = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1578     table.r             = fr->gbtabr;
1579     table.scale         = fr->gbtabscale;
1580     table.scale_exp     = 0;
1581     table.n             = table.scale*table.r;
1582     table.formatsize    = 4;
1583     table.ninteractions = 1;
1584     table.stride        = table.formatsize*table.ninteractions;
1585     nx0                 = 0;
1586     nx                  = table.scale*table.r;
1587
1588     /* Check whether we have to read or generate
1589      * We will always generate a table, so remove the read code
1590      * (Compare with original make_table function
1591      */
1592     bReadTab = FALSE;
1593     bGenTab  = TRUE;
1594
1595     /* Each table type (e.g. coul,lj6,lj12) requires four
1596      * numbers per datapoint. For performance reasons we want
1597      * the table data to be aligned to 16-byte. This is accomplished
1598      * by allocating 16 bytes extra to a temporary pointer, and then
1599      * calculating an aligned pointer. This new pointer must not be
1600      * used in a free() call, but thankfully we're sloppy enough not
1601      * to do this :-)
1602      */
1603
1604     snew_aligned(table.data, 4*nx, 32);
1605
1606     init_table(nx, nx0, table.scale, &(td[0]), !bReadTab);
1607
1608     /* Local implementation so we don't have to use the etabGB
1609      * enum above, which will cause problems later when
1610      * making the other tables (right now even though we are using
1611      * GB, the normal Coulomb tables will be created, but this
1612      * will cause a problem since fr->eeltype==etabGB which will not
1613      * be defined in fill_table and set_table_type
1614      */
1615
1616     for (i = nx0; i < nx; i++)
1617     {
1618         r       = td->x[i];
1619         r2      = r*r;
1620         expterm = exp(-0.25*r2);
1621
1622         Vtab = 1/sqrt(r2+expterm);
1623         Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1624
1625         /* Convert to single precision when we store to mem */
1626         td->v[i]  = Vtab;
1627         td->f[i]  = Ftab;
1628
1629     }
1630
1631     copy2table(table.n, 0, 4, td[0].x, td[0].v, td[0].f, 1.0, table.data);
1632
1633     if (bDebugMode())
1634     {
1635         fp = xvgropen(fns[0], fns[0], "r", "V", oenv);
1636         /* plot the output 5 times denser than the table data */
1637         /* for(i=5*nx0;i<5*table.n;i++) */
1638         for (i = nx0; i < table.n; i++)
1639         {
1640             /* x0=i*table.r/(5*table.n); */
1641             x0 = i*table.r/table.n;
1642             evaluate_table(table.data, 0, 4, table.scale, x0, &y0, &yp);
1643             fprintf(fp, "%15.10e  %15.10e  %15.10e\n", x0, y0, yp);
1644
1645         }
1646         gmx_fio_fclose(fp);
1647     }
1648
1649     /*
1650        for(i=100*nx0;i<99.81*table.n;i++)
1651        {
1652        r = i*table.r/(100*table.n);
1653        r2      = r*r;
1654        expterm = exp(-0.25*r2);
1655
1656        Vtab = 1/sqrt(r2+expterm);
1657        Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1658
1659
1660        evaluate_table(table.data,0,4,table.scale,r,&y0,&yp);
1661        printf("gb: i=%d, x0=%g, y0=%15.15f, Vtab=%15.15f, yp=%15.15f, Ftab=%15.15f\n",i,r, y0, Vtab, yp, Ftab);
1662
1663        abs_error_r=fabs(y0-Vtab);
1664        abs_error_r2=fabs(yp-(-1)*Ftab);
1665
1666        rel_error_r=abs_error_r/y0;
1667        rel_error_r2=fabs(abs_error_r2/yp);
1668
1669
1670        if(rel_error_r>rel_error_r_old)
1671        {
1672        rel_error_r_old=rel_error_r;
1673        x0_r_error=x0;
1674        }
1675
1676        if(rel_error_r2>rel_error_r2_old)
1677        {
1678        rel_error_r2_old=rel_error_r2;
1679        x0_r2_error=x0;
1680        }
1681        }
1682
1683        printf("gb: MAX REL ERROR IN R=%15.15f, MAX REL ERROR IN R2=%15.15f\n",rel_error_r_old, rel_error_r2_old);
1684        printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error);
1685
1686        exit(1); */
1687     done_tabledata(&(td[0]));
1688     sfree(td);
1689
1690     return table;
1691
1692
1693 }
1694
1695 t_forcetable make_atf_table(FILE *out, const output_env_t oenv,
1696                             const t_forcerec *fr,
1697                             const char *fn,
1698                             matrix box)
1699 {
1700     const char  *fns[3] = { "tf_tab.xvg", "atfdtab.xvg", "atfrtab.xvg" };
1701     FILE        *fp;
1702     t_tabledata *td;
1703     real         x0, y0, yp, rtab;
1704     int          i, nx, nx0;
1705     real         rx, ry, rz, box_r;
1706
1707     t_forcetable table;
1708
1709
1710     /* Set the table dimensions for ATF, not really necessary to
1711      * use etiNR (since we only have one table, but ...)
1712      */
1713     snew(td, 1);
1714
1715     if (fr->adress_type == eAdressSphere)
1716     {
1717         /* take half box diagonal direction as tab range */
1718         rx    = 0.5*box[0][0]+0.5*box[1][0]+0.5*box[2][0];
1719         ry    = 0.5*box[0][1]+0.5*box[1][1]+0.5*box[2][1];
1720         rz    = 0.5*box[0][2]+0.5*box[1][2]+0.5*box[2][2];
1721         box_r = sqrt(rx*rx+ry*ry+rz*rz);
1722
1723     }
1724     else
1725     {
1726         /* xsplit: take half box x direction as tab range */
1727         box_r        = box[0][0]/2;
1728     }
1729     table.r         = box_r;
1730     table.scale     = 0;
1731     table.n         = 0;
1732     table.scale_exp = 0;
1733     nx0             = 10;
1734     nx              = 0;
1735
1736     read_tables(out, fn, 1, 0, td);
1737     rtab      = td[0].x[td[0].nx-1];
1738
1739     if (fr->adress_type == eAdressXSplit && (rtab < box[0][0]/2))
1740     {
1741         gmx_fatal(FARGS, "AdResS full box therm force table in file %s extends to %f:\n"
1742                   "\tshould extend to at least half the length of the box in x-direction"
1743                   "%f\n", fn, rtab, box[0][0]/2);
1744     }
1745     if (rtab < box_r)
1746     {
1747         gmx_fatal(FARGS, "AdResS full box therm force table in file %s extends to %f:\n"
1748                   "\tshould extend to at least for spherical adress"
1749                   "%f (=distance from center to furthermost point in box \n", fn, rtab, box_r);
1750     }
1751
1752
1753     table.n     = td[0].nx;
1754     nx          = table.n;
1755     table.scale = td[0].tabscale;
1756     nx0         = td[0].nx0;
1757
1758     /* Each table type (e.g. coul,lj6,lj12) requires four
1759      * numbers per datapoint. For performance reasons we want
1760      * the table data to be aligned to 16-byte. This is accomplished
1761      * by allocating 16 bytes extra to a temporary pointer, and then
1762      * calculating an aligned pointer. This new pointer must not be
1763      * used in a free() call, but thankfully we're sloppy enough not
1764      * to do this :-)
1765      */
1766
1767     snew_aligned(table.data, 4*nx, 32);
1768
1769     copy2table(table.n, 0, 4, td[0].x, td[0].v, td[0].f, 1.0, table.data);
1770
1771     if (bDebugMode())
1772     {
1773         fp = xvgropen(fns[0], fns[0], "r", "V", oenv);
1774         /* plot the output 5 times denser than the table data */
1775         /* for(i=5*nx0;i<5*table.n;i++) */
1776
1777         for (i = 5*((nx0+1)/2); i < 5*table.n; i++)
1778         {
1779             /* x0=i*table.r/(5*table.n); */
1780             x0 = i*table.r/(5*(table.n-1));
1781             evaluate_table(table.data, 0, 4, table.scale, x0, &y0, &yp);
1782             fprintf(fp, "%15.10e  %15.10e  %15.10e\n", x0, y0, yp);
1783
1784         }
1785         gmx_ffclose(fp);
1786     }
1787
1788     done_tabledata(&(td[0]));
1789     sfree(td);
1790
1791     table.interaction   = GMX_TABLE_INTERACTION_ELEC_VDWREP_VDWDISP;
1792     table.format        = GMX_TABLE_FORMAT_CUBICSPLINE_YFGH;
1793     table.formatsize    = 4;
1794     table.ninteractions = 3;
1795     table.stride        = table.formatsize*table.ninteractions;
1796
1797
1798     return table;
1799 }
1800
1801 bondedtable_t make_bonded_table(FILE *fplog, char *fn, int angle)
1802 {
1803     t_tabledata   td;
1804     double        start;
1805     int           i;
1806     bondedtable_t tab;
1807
1808     if (angle < 2)
1809     {
1810         start = 0;
1811     }
1812     else
1813     {
1814         start = -180.0;
1815     }
1816     read_tables(fplog, fn, 1, angle, &td);
1817     if (angle > 0)
1818     {
1819         /* Convert the table from degrees to radians */
1820         for (i = 0; i < td.nx; i++)
1821         {
1822             td.x[i] *= DEG2RAD;
1823             td.f[i] *= RAD2DEG;
1824         }
1825         td.tabscale *= RAD2DEG;
1826     }
1827     tab.n     = td.nx;
1828     tab.scale = td.tabscale;
1829     snew(tab.data, tab.n*4);
1830     copy2table(tab.n, 0, 4, td.x, td.v, td.f, 1.0, tab.data);
1831     done_tabledata(&td);
1832
1833     return tab;
1834 }