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