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