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