src/gromacs/mdlib/tables.c

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