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