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