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