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