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