src/gromacs/mdlib/tables.c

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