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