src/gromacs/mdlib/tables.c

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