src/mdlib/tables.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROwing Monsters And Cloning Shrimps
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "maths.h"
  41 #include "typedefs.h"
  42 #include "names.h"
  43 #include "smalloc.h"
  44 #include "gmx_fatal.h"
  45 #include "futil.h"
  46 #include "xvgr.h"
  47 #include "vec.h"
  48 #include "main.h"
  49 #include "network.h"
  50 #include "physics.h"
  51 #include "force.h"
  52 #include "gmxfio.h"
  53
  54 /* All the possible (implemented) table functions */
  55 enum {
  56   etabLJ6,
  57   etabLJ12,
  58   etabLJ6Shift,
  59   etabLJ12Shift,
  60   etabShift,
  61   etabRF,
  62   etabRF_ZERO,
  63   etabCOUL,
  64   etabEwald,
  65   etabEwaldSwitch,
  66   etabEwaldUser,
  67   etabEwaldUserSwitch,
  68   etabLJ6Switch,
  69   etabLJ12Switch,
  70   etabCOULSwitch,
  71   etabLJ6Encad,
  72   etabLJ12Encad,
  73   etabCOULEncad,
  74   etabEXPMIN,
  75   etabUSER,
  76   etabNR
  77 };
  78
  79 /** Evaluates to true if the table type contains user data. */
  80 #define ETAB_USER(e)  ((e) == etabUSER || \
  81                        (e) == etabEwaldUser || (e) == etabEwaldUserSwitch)
  82
  83 typedef struct {
  84   const char *name;
  85   gmx_bool bCoulomb;
  86 } t_tab_props;
  87
  88 /* This structure holds name and a flag that tells whether
  89    this is a Coulomb type funtion */
  90 static const t_tab_props tprops[etabNR] = {
  91   { "LJ6",  FALSE },
  92   { "LJ12", FALSE },
  93   { "LJ6Shift", FALSE },
  94   { "LJ12Shift", FALSE },
  95   { "Shift", TRUE },
  96   { "RF", TRUE },
  97   { "RF-zero", TRUE },
  98   { "COUL", TRUE },
  99   { "Ewald", TRUE },
 100   { "Ewald-Switch", TRUE },
 101   { "Ewald-User", TRUE },
 102   { "Ewald-User-Switch", TRUE },
 103   { "LJ6Switch", FALSE },
 104   { "LJ12Switch", FALSE },
 105   { "COULSwitch", TRUE },
 106   { "LJ6-Encad shift", FALSE },
 107   { "LJ12-Encad shift", FALSE },
 108   { "COUL-Encad shift",  TRUE },
 109   { "EXPMIN", FALSE },
 110   { "USER", FALSE }
 111 };
 112
 113 /* Index in the table that says which function to use */
 114 enum { etiCOUL, etiLJ6, etiLJ12, etiNR };
 115
 116 typedef struct {
 117   int  nx,nx0;
 118   double tabscale;
 119   double *x,*v,*f;
 120 } t_tabledata;
 121
 122 #define pow2(x) ((x)*(x))
 123 #define pow3(x) ((x)*(x)*(x))
 124 #define pow4(x) ((x)*(x)*(x)*(x))
 125 #define pow5(x) ((x)*(x)*(x)*(x)*(x))
 126
 127 /* Calculate the potential and force for an r value
 128  * in exactly the same way it is done in the inner loop.
 129  * VFtab is a pointer to the table data, offset is
 130  * the point where we should begin and stride is
 131  * 4 if we have a buckingham table, 3 otherwise.
 132  * If you want to evaluate table no N, set offset to 4*N.
 133  *
 134  * We use normal precision here, since that is what we
 135  * will use in the inner loops.
 136  */
 137 static void evaluate_table(real VFtab[], int offset, int stride,
 138                            real tabscale, real r, real *y, real *yp)
 139 {
 140   int n;
 141   real rt,eps,eps2;
 142   real Y,F,Geps,Heps2,Fp;
 143
 144   rt       =  r*tabscale;
 145   n        =  (int)rt;
 146   eps      =  rt - n;
 147   eps2     =  eps*eps;
 148   n        =  offset+stride*n;
 149   Y        =  VFtab[n];
 150   F        =  VFtab[n+1];
 151   Geps     =  eps*VFtab[n+2];
 152   Heps2    =  eps2*VFtab[n+3];
 153   Fp       =  F+Geps+Heps2;
 154   *y       =  Y+eps*Fp;
 155   *yp      =  (Fp+Geps+2.0*Heps2)*tabscale;
 156 }
 157
 158
 159 static void splint(real xa[],real ya[],real y2a[],
 160                    int n,real x,real *y,real *yp)
 161 {
 162   int  klo,khi,k;
 163   real h,b,a,eps;
 164   real F,G,H;
 165
 166   klo=1;
 167   khi=n;
 168
 169   while ((khi-klo) > 1) {
 170     k=(khi+klo) >> 1;
 171     if (xa[k] > x)
 172       khi=k;
 173     else
 174       klo=k;
 175   }
 176   h = xa[khi]-xa[klo];
 177   if (h == 0.0)
 178     gmx_fatal(FARGS,"Bad XA input to function splint");
 179   a   = (xa[khi]-x)/h;
 180   b   = (x-xa[klo])/h;
 181   *y  = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
 182   *yp = (ya[khi]-ya[klo])/h+((3*a*a-1)*y2a[klo]-(3*b*b-1)*y2a[khi])*h/6.0;
 183
 184   eps = b;
 185   F   = (ya[khi]-ya[klo]-(h*h/6.0)*(2*y2a[klo]+y2a[khi]));
 186   G   = (h*h/2.0)*y2a[klo];
 187   H   = (h*h/6.0)*(y2a[khi]-y2a[klo]);
 188   *y  = ya[klo] + eps*F + eps*eps*G + eps*eps*eps*H;
 189   *yp = (F + 2*eps*G + 3*eps*eps*H)/h;
 190 }
 191
 192
 193 static void copy2table(int n,int offset,int stride,
 194                        double x[],double Vtab[],double Ftab[],
 195                        real dest[])
 196 {
 197 /* Use double prec. for the intermediary variables
 198  * and temporary x/vtab/vtab2 data to avoid unnecessary
 199  * loss of precision.
 200  */
 201   int  i,nn0;
 202   double F,G,H,h;
 203
 204   h = 0;
 205   for(i=0; (i<n); i++) {
 206     if (i < n-1) {
 207       h   = x[i+1] - x[i];
 208       F   = -Ftab[i]*h;
 209       G   =  3*(Vtab[i+1] - Vtab[i]) + (Ftab[i+1] + 2*Ftab[i])*h;
 210       H   = -2*(Vtab[i+1] - Vtab[i]) - (Ftab[i+1] +   Ftab[i])*h;
 211     } else {
 212       /* Fill the last entry with a linear potential,
 213        * this is mainly for rounding issues with angle and dihedral potentials.
 214        */
 215       F   = -Ftab[i]*h;
 216       G   = 0;
 217       H   = 0;
 218     }
 219     nn0 = offset + i*stride;
 220     dest[nn0]   = Vtab[i];
 221     dest[nn0+1] = F;
 222     dest[nn0+2] = G;
 223     dest[nn0+3] = H;
 224   }
 225 }
 226
 227 static void init_table(FILE *fp,int n,int nx0,
 228                        double tabscale,t_tabledata *td,gmx_bool bAlloc)
 229 {
 230   int i;
 231
 232   td->nx  = n;
 233   td->nx0 = nx0;
 234   td->tabscale = tabscale;
 235   if (bAlloc) {
 236     snew(td->x,td->nx);
 237     snew(td->v,td->nx);
 238     snew(td->f,td->nx);
 239   }
 240   for(i=0; (i<td->nx); i++)
 241     td->x[i] = i/tabscale;
 242 }
 243
 244 static void spline_forces(int nx,double h,double v[],gmx_bool bS3,gmx_bool bE3,
 245                           double f[])
 246 {
 247   int    start,end,i;
 248   double v3,b_s,b_e,b;
 249   double beta,*gamma;
 250
 251   /* Formulas can be found in:
 252    * H.J.C. Berendsen, Simulating the Physical World, Cambridge 2007
 253    */
 254
 255   if (nx < 4 && (bS3 || bE3))
 256     gmx_fatal(FARGS,"Can not generate splines with third derivative boundary conditions with less than 4 (%d) points",nx);
 257
 258   /* To make life easy we initially set the spacing to 1
 259    * and correct for this at the end.
 260    */
 261   beta = 2;
 262   if (bS3) {
 263     /* Fit V''' at the start */
 264     v3  = v[3] - 3*v[2] + 3*v[1] - v[0];
 265     if (debug)
 266       fprintf(debug,"The left third derivative is %g\n",v3/(h*h*h));
 267     b_s = 2*(v[1] - v[0]) + v3/6;
 268     start = 0;
 269
 270     if (FALSE) {
 271       /* Fit V'' at the start */
 272       real v2;
 273
 274       v2  = -v[3] + 4*v[2] - 5*v[1] + 2*v[0];
 275       /* v2  = v[2] - 2*v[1] + v[0]; */
 276       if (debug)
 277         fprintf(debug,"The left second derivative is %g\n",v2/(h*h));
 278       b_s = 3*(v[1] - v[0]) - v2/2;
 279       start = 0;
 280     }
 281   } else {
 282     b_s = 3*(v[2] - v[0]) + f[0]*h;
 283     start = 1;
 284   }
 285   if (bE3) {
 286     /* Fit V''' at the end */
 287     v3  = v[nx-1] - 3*v[nx-2] + 3*v[nx-3] - v[nx-4];
 288     if (debug)
 289       fprintf(debug,"The right third derivative is %g\n",v3/(h*h*h));
 290     b_e = 2*(v[nx-1] - v[nx-2]) + v3/6;
 291     end = nx;
 292   } else {
 293     /* V'=0 at the end */
 294     b_e = 3*(v[nx-1] - v[nx-3]) + f[nx-1]*h;
 295     end = nx - 1;
 296   }
 297
 298   snew(gamma,nx);
 299   beta = (bS3 ? 1 : 4);
 300
 301   /* For V'' fitting */
 302   /* beta = (bS3 ? 2 : 4); */
 303
 304   f[start] = b_s/beta;
 305   for(i=start+1; i<end; i++) {
 306     gamma[i] = 1/beta;
 307     beta = 4 - gamma[i];
 308     b    =  3*(v[i+1] - v[i-1]);
 309     f[i] = (b - f[i-1])/beta;
 310   }
 311   gamma[end-1] = 1/beta;
 312   beta = (bE3 ? 1 : 4) - gamma[end-1];
 313   f[end-1] = (b_e - f[end-2])/beta;
 314
 315   for(i=end-2; i>=start; i--)
 316     f[i] -= gamma[i+1]*f[i+1];
 317   sfree(gamma);
 318
 319   /* Correct for the minus sign and the spacing */
 320   for(i=start; i<end; i++)
 321     f[i] = -f[i]/h;
 322 }
 323
 324 static void set_forces(FILE *fp,int angle,
 325                        int nx,double h,double v[],double f[],
 326                        int table)
 327 {
 328   int start,end;
 329
 330   if (angle == 2)
 331     gmx_fatal(FARGS,
 332               "Force generation for dihedral tables is not (yet) implemented");
 333
 334   start = 0;
 335   while (v[start] == 0)
 336     start++;
 337
 338   end = nx;
 339   while(v[end-1] == 0)
 340     end--;
 341   if (end > nx - 2)
 342     end = nx;
 343   else
 344     end++;
 345
 346   if (fp)
 347     fprintf(fp,"Generating forces for table %d, boundary conditions: V''' at %g, %s at %g\n",
 348             table+1,start*h,end==nx ? "V'''" : "V'=0",(end-1)*h);
 349   spline_forces(end-start,h,v+start,TRUE,end==nx,f+start);
 350 }
 351
 352 static void read_tables(FILE *fp,const char *fn,
 353                         int ntab,int angle,t_tabledata td[])
 354 {
 355   char *libfn;
 356   char buf[STRLEN];
 357   double **yy=NULL,start,end,dx0,dx1,ssd,vm,vp,f,numf;
 358   int  k,i,nx,nx0=0,ny,nny,ns;
 359   gmx_bool bAllZero,bZeroV,bZeroF;
 360   double tabscale;
 361
 362   nny = 2*ntab+1;
 363   libfn = gmxlibfn(fn);
 364   nx  = read_xvg(libfn,&yy,&ny);
 365   if (ny != nny)
 366     gmx_fatal(FARGS,"Trying to read file %s, but nr columns = %d, should be %d",
 367                 libfn,ny,nny);
 368   if (angle == 0) {
 369     if (yy[0][0] != 0.0)
 370       gmx_fatal(FARGS,
 371                 "The first distance in file %s is %f nm instead of %f nm",
 372                 libfn,yy[0][0],0.0);
 373   } else {
 374     if (angle == 1)
 375       start = 0.0;
 376     else
 377       start = -180.0;
 378     end = 180.0;
 379     if (yy[0][0] != start || yy[0][nx-1] != end)
 380       gmx_fatal(FARGS,"The angles in file %s should go from %f to %f instead of %f to %f\n",
 381                 libfn,start,end,yy[0][0],yy[0][nx-1]);
 382   }
 383
 384   tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
 385
 386   if (fp) {
 387     fprintf(fp,"Read user tables from %s with %d data points.\n",libfn,nx);
 388     if (angle == 0)
 389       fprintf(fp,"Tabscale = %g points/nm\n",tabscale);
 390   }
 391
 392   bAllZero = TRUE;
 393   for(k=0; k<ntab; k++) {
 394     bZeroV = TRUE;
 395     bZeroF = TRUE;
 396     for(i=0; (i < nx); i++) {
 397       if (i >= 2) {
 398         dx0 = yy[0][i-1] - yy[0][i-2];
 399         dx1 = yy[0][i]   - yy[0][i-1];
 400         /* Check for 1% deviation in spacing */
 401         if (fabs(dx1 - dx0) >= 0.005*(fabs(dx0) + fabs(dx1))) {
 402           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]);
 403         }
 404       }
 405       if (yy[1+k*2][i] != 0) {
 406         bZeroV = FALSE;
 407         if (bAllZero) {
 408           bAllZero = FALSE;
 409           nx0 = i;
 410         }
 411         if (yy[1+k*2][i] >  0.01*GMX_REAL_MAX ||
 412             yy[1+k*2][i] < -0.01*GMX_REAL_MAX) {
 413           gmx_fatal(FARGS,"Out of range potential value %g in file '%s'",
 414                     yy[1+k*2][i],fn);
 415         }
 416       }
 417       if (yy[1+k*2+1][i] != 0) {
 418         bZeroF = FALSE;
 419         if (bAllZero) {
 420           bAllZero = FALSE;
 421           nx0 = i;
 422         }
 423         if (yy[1+k*2+1][i] >  0.01*GMX_REAL_MAX ||
 424             yy[1+k*2+1][i] < -0.01*GMX_REAL_MAX) {
 425           gmx_fatal(FARGS,"Out of range force value %g in file '%s'",
 426                     yy[1+k*2+1][i],fn);
 427         }
 428       }
 429     }
 430
 431     if (!bZeroV && bZeroF) {
 432       set_forces(fp,angle,nx,1/tabscale,yy[1+k*2],yy[1+k*2+1],k);
 433     } else {
 434       /* Check if the second column is close to minus the numerical
 435        * derivative of the first column.
 436        */
 437       ssd = 0;
 438       ns = 0;
 439       for(i=1; (i < nx-1); i++) {
 440         vm = yy[1+2*k][i-1];
 441         vp = yy[1+2*k][i+1];
 442         f  = yy[1+2*k+1][i];
 443         if (vm != 0 && vp != 0 && f != 0) {
 444           /* Take the centered difference */
 445           numf = -(vp - vm)*0.5*tabscale;
 446           ssd += fabs(2*(f - numf)/(f + numf));
 447           ns++;
 448         }
 449       }
 450       if (ns > 0) {
 451         ssd /= ns;
 452         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));
 453         if (debug)
 454           fprintf(debug,"%s",buf);
 455         if (ssd > 0.2) {
 456           if (fp)
 457             fprintf(fp,"\nWARNING: %s\n",buf);
 458           fprintf(stderr,"\nWARNING: %s\n",buf);
 459         }
 460       }
 461     }
 462   }
 463   if (bAllZero && fp) {
 464     fprintf(fp,"\nNOTE: All elements in table %s are zero\n\n",libfn);
 465   }
 466
 467   for(k=0; (k<ntab); k++) {
 468     init_table(fp,nx,nx0,tabscale,&(td[k]),TRUE);
 469     for(i=0; (i<nx); i++) {
 470       td[k].x[i] = yy[0][i];
 471       td[k].v[i] = yy[2*k+1][i];
 472       td[k].f[i] = yy[2*k+2][i];
 473     }
 474   }
 475   for(i=0; (i<ny); i++)
 476     sfree(yy[i]);
 477   sfree(yy);
 478   sfree(libfn);
 479 }
 480
 481 static void done_tabledata(t_tabledata *td)
 482 {
 483   int i;
 484
 485   if (!td)
 486     return;
 487
 488   sfree(td->x);
 489   sfree(td->v);
 490   sfree(td->f);
 491 }
 492
 493 static void fill_table(t_tabledata *td,int tp,const t_forcerec *fr)
 494 {
 495   /* Fill the table according to the formulas in the manual.
 496    * In principle, we only need the potential and the second
 497    * derivative, but then we would have to do lots of calculations
 498    * in the inner loop. By precalculating some terms (see manual)
 499    * we get better eventual performance, despite a larger table.
 500    *
 501    * Since some of these higher-order terms are very small,
 502    * we always use double precision to calculate them here, in order
 503    * to avoid unnecessary loss of precision.
 504    */
 505 #ifdef DEBUG_SWITCH
 506   FILE *fp;
 507 #endif
 508   int  i;
 509   double reppow,p;
 510   double r1,rc,r12,r13;
 511   double r,r2,r6,rc6;
 512   double expr,Vtab,Ftab;
 513   /* Parameters for David's function */
 514   double A=0,B=0,C=0,A_3=0,B_4=0;
 515   /* Parameters for the switching function */
 516   double ksw,swi,swi1;
 517   /* Temporary parameters */
 518   gmx_bool bSwitch,bShift;
 519   double ewc=fr->ewaldcoeff;
 520   double isp= 0.564189583547756;
 521
 522   bSwitch = ((tp == etabLJ6Switch) || (tp == etabLJ12Switch) ||
 523              (tp == etabCOULSwitch) ||
 524              (tp == etabEwaldSwitch) || (tp == etabEwaldUserSwitch));
 525   bShift  = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
 526              (tp == etabShift));
 527
 528   reppow = fr->reppow;
 529
 530   if (tprops[tp].bCoulomb) {
 531     r1 = fr->rcoulomb_switch;
 532     rc = fr->rcoulomb;
 533   }
 534   else {
 535     r1 = fr->rvdw_switch;
 536     rc = fr->rvdw;
 537   }
 538   if (bSwitch)
 539     ksw  = 1.0/(pow5(rc-r1));
 540   else
 541     ksw  = 0.0;
 542   if (bShift) {
 543     if (tp == etabShift)
 544       p = 1;
 545     else if (tp == etabLJ6Shift)
 546       p = 6;
 547     else
 548       p = reppow;
 549
 550     A = p * ((p+1)*r1-(p+4)*rc)/(pow(rc,p+2)*pow2(rc-r1));
 551     B = -p * ((p+1)*r1-(p+3)*rc)/(pow(rc,p+2)*pow3(rc-r1));
 552     C = 1.0/pow(rc,p)-A/3.0*pow3(rc-r1)-B/4.0*pow4(rc-r1);
 553     if (tp == etabLJ6Shift) {
 554       A=-A;
 555       B=-B;
 556       C=-C;
 557     }
 558     A_3=A/3.0;
 559     B_4=B/4.0;
 560   }
 561   if (debug) { fprintf(debug,"Setting up tables\n"); fflush(debug); }
 562
 563 #ifdef DEBUG_SWITCH
 564   fp=xvgropen("switch.xvg","switch","r","s");
 565 #endif
 566
 567   for(i=td->nx0; (i<td->nx); i++) {
 568     r     = td->x[i];
 569     r2    = r*r;
 570     r6    = 1.0/(r2*r2*r2);
 571     if (gmx_within_tol(reppow,12.0,10*GMX_DOUBLE_EPS)) {
 572       r12 = r6*r6;
 573     } else {
 574       r12 = pow(r,-reppow);
 575     }
 576     Vtab  = 0.0;
 577     Ftab  = 0.0;
 578     if (bSwitch) {
 579       /* swi is function, swi1 1st derivative and swi2 2nd derivative */
 580       /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
 581        * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
 582        * r1 and rc.
 583        * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
 584        */
 585       if(r<=r1) {
 586         swi  = 1.0;
 587         swi1 = 0.0;
 588       } else if (r>=rc) {
 589         swi  = 0.0;
 590         swi1 = 0.0;
 591       } else {
 592         swi      = 1 - 10*pow3(r-r1)*ksw*pow2(rc-r1)
 593           + 15*pow4(r-r1)*ksw*(rc-r1) - 6*pow5(r-r1)*ksw;
 594         swi1     = -30*pow2(r-r1)*ksw*pow2(rc-r1)
 595           + 60*pow3(r-r1)*ksw*(rc-r1) - 30*pow4(r-r1)*ksw;
 596       }
 597     }
 598     else { /* not really needed, but avoids compiler warnings... */
 599       swi  = 1.0;
 600       swi1 = 0.0;
 601     }
 602 #ifdef DEBUG_SWITCH
 603     fprintf(fp,"%10g  %10g  %10g  %10g\n",r,swi,swi1,swi2);
 604 #endif
 605
 606     rc6 = rc*rc*rc;
 607     rc6 = 1.0/(rc6*rc6);
 608
 609     switch (tp) {
 610     case etabLJ6:
 611       /* Dispersion */
 612       Vtab  = -r6;
 613       Ftab  = 6.0*Vtab/r;
 614       break;
 615     case etabLJ6Switch:
 616     case etabLJ6Shift:
 617       /* Dispersion */
 618       if (r < rc) {
 619         Vtab  = -r6;
 620         Ftab  = 6.0*Vtab/r;
 621       }
 622       break;
 623     case etabLJ12:
 624       /* Repulsion */
 625       Vtab  = r12;
 626       Ftab  = reppow*Vtab/r;
 627       break;
 628     case etabLJ12Switch:
 629     case etabLJ12Shift:
 630       /* Repulsion */
 631       if (r < rc) {
 632         Vtab  = r12;
 633         Ftab  = reppow*Vtab/r;
 634       }
 635       break;
 636         case etabLJ6Encad:
 637         if(r < rc) {
 638             Vtab  = -(r6-6.0*(rc-r)*rc6/rc-rc6);
 639             Ftab  = -(6.0*r6/r-6.0*rc6/rc);
 640         } else { /* r>rc */
 641             Vtab  = 0;
 642             Ftab  = 0;
 643         }
 644         break;
 645     case etabLJ12Encad:
 646         if(r < rc) {
 647             Vtab  = r12-12.0*(rc-r)*rc6*rc6/rc-1.0*rc6*rc6;
 648             Ftab  = 12.0*r12/r-12.0*rc6*rc6/rc;
 649         } else { /* r>rc */
 650             Vtab  = 0;
 651             Ftab  = 0;
 652         }
 653         break;
 654     case etabCOUL:
 655       Vtab  = 1.0/r;
 656       Ftab  = 1.0/r2;
 657       break;
 658     case etabCOULSwitch:
 659     case etabShift:
 660       if (r < rc) {
 661         Vtab  = 1.0/r;
 662         Ftab  = 1.0/r2;
 663       }
 664       break;
 665     case etabEwald:
 666     case etabEwaldSwitch:
 667       Vtab  = gmx_erfc(ewc*r)/r;
 668       Ftab  = gmx_erfc(ewc*r)/r2+2*exp(-(ewc*ewc*r2))*ewc*isp/r;
 669       break;
 670     case etabEwaldUser:
 671     case etabEwaldUserSwitch:
 672       /* Only calculate minus the reciprocal space contribution */
 673       Vtab  = -gmx_erf(ewc*r)/r;
 674       Ftab  = -gmx_erf(ewc*r)/r2+2*exp(-(ewc*ewc*r2))*ewc*isp/r;
 675       break;
 676     case etabRF:
 677     case etabRF_ZERO:
 678       Vtab  = 1.0/r      +   fr->k_rf*r2 - fr->c_rf;
 679       Ftab  = 1.0/r2     - 2*fr->k_rf*r;
 680       if (tp == etabRF_ZERO && r >= rc) {
 681         Vtab = 0;
 682         Ftab = 0;
 683       }
 684       break;
 685     case etabEXPMIN:
 686       expr  = exp(-r);
 687       Vtab  = expr;
 688       Ftab  = expr;
 689       break;
 690     case etabCOULEncad:
 691         if(r < rc) {
 692             Vtab  = 1.0/r-(rc-r)/(rc*rc)-1.0/rc;
 693             Ftab  = 1.0/r2-1.0/(rc*rc);
 694         } else { /* r>rc */
 695             Vtab  = 0;
 696             Ftab  = 0;
 697         }
 698         break;
 699     default:
 700       gmx_fatal(FARGS,"Table type %d not implemented yet. (%s,%d)",
 701                   tp,__FILE__,__LINE__);
 702     }
 703     if (bShift) {
 704       /* Normal coulomb with cut-off correction for potential */
 705       if (r < rc) {
 706         Vtab -= C;
 707         /* If in Shifting range add something to it */
 708         if (r > r1) {
 709           r12 = (r-r1)*(r-r1);
 710           r13 = (r-r1)*r12;
 711           Vtab  += - A_3*r13 - B_4*r12*r12;
 712           Ftab  +=   A*r12 + B*r13;
 713         }
 714       }
 715     }
 716
 717     if (ETAB_USER(tp)) {
 718       Vtab += td->v[i];
 719       Ftab += td->f[i];
 720     }
 721
 722     if ((r > r1) && bSwitch) {
 723       Ftab = Ftab*swi - Vtab*swi1;
 724       Vtab = Vtab*swi;
 725     }
 726
 727     /* Convert to single precision when we store to mem */
 728     td->v[i]  = Vtab;
 729     td->f[i]  = Ftab;
 730   }
 731
 732   /* Continue the table linearly from nx0 to 0.
 733    * These values are only required for energy minimization with overlap or TPI.
 734    */
 735   for(i=td->nx0-1; i>=0; i--) {
 736     td->v[i] = td->v[i+1] + td->f[i+1]*(td->x[i+1] - td->x[i]);
 737     td->f[i] = td->f[i+1];
 738   }
 739
 740 #ifdef DEBUG_SWITCH
 741   gmx_fio_fclose(fp);
 742 #endif
 743 }
 744
 745 static void set_table_type(int tabsel[],const t_forcerec *fr,gmx_bool b14only)
 746 {
 747   int eltype,vdwtype;
 748
 749   /* Set the different table indices.
 750    * Coulomb first.
 751    */
 752
 753
 754   if (b14only) {
 755     switch (fr->eeltype) {
 756     case eelRF_NEC:
 757       eltype = eelRF;
 758       break;
 759     case eelUSER:
 760     case eelPMEUSER:
 761     case eelPMEUSERSWITCH:
 762       eltype = eelUSER;
 763       break;
 764     default:
 765       eltype = eelCUT;
 766     }
 767   } else {
 768     eltype = fr->eeltype;
 769   }
 770
 771   switch (eltype) {
 772   case eelCUT:
 773     tabsel[etiCOUL] = etabCOUL;
 774     break;
 775   case eelPPPM:
 776   case eelPOISSON:
 777     tabsel[etiCOUL] = etabShift;
 778     break;
 779   case eelSHIFT:
 780     if (fr->rcoulomb > fr->rcoulomb_switch)
 781       tabsel[etiCOUL] = etabShift;
 782     else
 783       tabsel[etiCOUL] = etabCOUL;
 784     break;
 785   case eelEWALD:
 786   case eelPME:
 787     tabsel[etiCOUL] = etabEwald;
 788     break;
 789   case eelPMESWITCH:
 790     tabsel[etiCOUL] = etabEwaldSwitch;
 791     break;
 792   case eelPMEUSER:
 793     tabsel[etiCOUL] = etabEwaldUser;
 794     break;
 795   case eelPMEUSERSWITCH:
 796     tabsel[etiCOUL] = etabEwaldUserSwitch;
 797     break;
 798   case eelRF:
 799   case eelGRF:
 800   case eelRF_NEC:
 801     tabsel[etiCOUL] = etabRF;
 802     break;
 803   case eelRF_ZERO:
 804     tabsel[etiCOUL] = etabRF_ZERO;
 805     break;
 806   case eelSWITCH:
 807     tabsel[etiCOUL] = etabCOULSwitch;
 808     break;
 809   case eelUSER:
 810     tabsel[etiCOUL] = etabUSER;
 811     break;
 812   case eelENCADSHIFT:
 813     tabsel[etiCOUL] = etabCOULEncad;
 814     break;
 815   default:
 816     gmx_fatal(FARGS,"Invalid eeltype %d",eltype);
 817   }
 818
 819   /* Van der Waals time */
 820   if (fr->bBHAM && !b14only) {
 821     tabsel[etiLJ6]  = etabLJ6;
 822     tabsel[etiLJ12] = etabEXPMIN;
 823   } else {
 824     if (b14only && fr->vdwtype != evdwUSER)
 825       vdwtype = evdwCUT;
 826     else
 827       vdwtype = fr->vdwtype;
 828
 829     switch (vdwtype) {
 830     case evdwSWITCH:
 831       tabsel[etiLJ6]  = etabLJ6Switch;
 832       tabsel[etiLJ12] = etabLJ12Switch;
 833       break;
 834     case evdwSHIFT:
 835       tabsel[etiLJ6]  = etabLJ6Shift;
 836       tabsel[etiLJ12] = etabLJ12Shift;
 837       break;
 838     case evdwUSER:
 839       tabsel[etiLJ6]  = etabUSER;
 840       tabsel[etiLJ12] = etabUSER;
 841       break;
 842     case evdwCUT:
 843       tabsel[etiLJ6]  = etabLJ6;
 844       tabsel[etiLJ12] = etabLJ12;
 845       break;
 846     case evdwENCADSHIFT:
 847       tabsel[etiLJ6]  = etabLJ6Encad;
 848       tabsel[etiLJ12] = etabLJ12Encad;
 849       break;
 850     default:
 851       gmx_fatal(FARGS,"Invalid vdwtype %d in %s line %d",vdwtype,
 852                   __FILE__,__LINE__);
 853     }
 854   }
 855 }
 856
 857 t_forcetable make_tables(FILE *out,const output_env_t oenv,
 858                          const t_forcerec *fr,
 859                          gmx_bool bVerbose,const char *fn,
 860                          real rtab,int flags)
 861 {
 862   const char *fns[3] = { "ctab.xvg", "dtab.xvg", "rtab.xvg" };
 863   const char *fns14[3] = { "ctab14.xvg", "dtab14.xvg", "rtab14.xvg" };
 864   FILE        *fp;
 865   t_tabledata *td;
 866   gmx_bool        b14only,bReadTab,bGenTab;
 867   real        x0,y0,yp;
 868   int         i,j,k,nx,nx0,tabsel[etiNR];
 869
 870   t_forcetable table;
 871
 872   b14only = (flags & GMX_MAKETABLES_14ONLY);
 873
 874   if (flags & GMX_MAKETABLES_FORCEUSER) {
 875     tabsel[etiCOUL] = etabUSER;
 876     tabsel[etiLJ6]  = etabUSER;
 877     tabsel[etiLJ12] = etabUSER;
 878   } else {
 879     set_table_type(tabsel,fr,b14only);
 880   }
 881   snew(td,etiNR);
 882   table.r         = rtab;
 883   table.scale     = 0;
 884   table.n         = 0;
 885   table.scale_exp = 0;
 886   nx0             = 10;
 887   nx              = 0;
 888
 889   /* Check whether we have to read or generate */
 890   bReadTab = FALSE;
 891   bGenTab  = FALSE;
 892   for(i=0; (i<etiNR); i++) {
 893     if (ETAB_USER(tabsel[i]))
 894       bReadTab = TRUE;
 895     if (tabsel[i] != etabUSER)
 896       bGenTab  = TRUE;
 897   }
 898   if (bReadTab) {
 899     read_tables(out,fn,etiNR,0,td);
 900     if (rtab == 0 || (flags & GMX_MAKETABLES_14ONLY)) {
 901       rtab      = td[0].x[td[0].nx-1];
 902       table.n   = td[0].nx;
 903       nx        = table.n;
 904     } else {
 905       if (td[0].x[td[0].nx-1] < rtab)
 906         gmx_fatal(FARGS,"Tables in file %s not long enough for cut-off:\n"
 907                   "\tshould be at least %f nm\n",fn,rtab);
 908       nx        = table.n = (int)(rtab*td[0].tabscale + 0.5);
 909     }
 910     table.scale = td[0].tabscale;
 911     nx0         = td[0].nx0;
 912   }
 913   if (bGenTab) {
 914     if (!bReadTab) {
 915 #ifdef GMX_DOUBLE
 916       table.scale = 2000.0;
 917 #else
 918       table.scale = 500.0;
 919 #endif
 920       nx = table.n = rtab*table.scale;
 921     }
 922   }
 923   if (fr->bBHAM) {
 924     if(fr->bham_b_max!=0)
 925       table.scale_exp = table.scale/fr->bham_b_max;
 926     else
 927       table.scale_exp = table.scale;
 928   }
 929
 930   /* Each table type (e.g. coul,lj6,lj12) requires four
 931    * numbers per nx+1 data points. For performance reasons we want
 932    * the table data to be aligned to 16-byte.
 933    */
 934   snew_aligned(table.tab, 12*(nx+1)*sizeof(real),16);
 935
 936   for(k=0; (k<etiNR); k++) {
 937     if (tabsel[k] != etabUSER) {
 938       init_table(out,nx,nx0,
 939                  (tabsel[k] == etabEXPMIN) ? table.scale_exp : table.scale,
 940                  &(td[k]),!bReadTab);
 941       fill_table(&(td[k]),tabsel[k],fr);
 942       if (out)
 943         fprintf(out,"%s table with %d data points for %s%s.\n"
 944                 "Tabscale = %g points/nm\n",
 945                 ETAB_USER(tabsel[k]) ? "Modified" : "Generated",
 946                 td[k].nx,b14only?"1-4 ":"",tprops[tabsel[k]].name,
 947                 td[k].tabscale);
 948     }
 949     copy2table(table.n,k*4,12,td[k].x,td[k].v,td[k].f,table.tab);
 950
 951     if (bDebugMode() && bVerbose) {
 952       if (b14only)
 953         fp=xvgropen(fns14[k],fns14[k],"r","V",oenv);
 954       else
 955         fp=xvgropen(fns[k],fns[k],"r","V",oenv);
 956       /* plot the output 5 times denser than the table data */
 957       for(i=5*((nx0+1)/2); i<5*table.n; i++) {
 958         x0 = i*table.r/(5*(table.n-1));
 959         evaluate_table(table.tab,4*k,12,table.scale,x0,&y0,&yp);
 960         fprintf(fp,"%15.10e  %15.10e  %15.10e\n",x0,y0,yp);
 961       }
 962       gmx_fio_fclose(fp);
 963     }
 964     done_tabledata(&(td[k]));
 965   }
 966   sfree(td);
 967
 968   return table;
 969 }
 970
 971 t_forcetable make_gb_table(FILE *out,const output_env_t oenv,
 972                            const t_forcerec *fr,
 973                            const char *fn,
 974                            real rtab)
 975 {
 976         const char *fns[3] = { "gbctab.xvg", "gbdtab.xvg", "gbrtab.xvg" };
 977         const char *fns14[3] = { "gbctab14.xvg", "gbdtab14.xvg", "gbrtab14.xvg" };
 978         FILE        *fp;
 979         t_tabledata *td;
 980         gmx_bool        bReadTab,bGenTab;
 981         real        x0,y0,yp;
 982         int         i,j,k,nx,nx0,tabsel[etiNR];
 983         void *      p_tmp;
 984         double      r,r2,Vtab,Ftab,expterm;
 985
 986         t_forcetable table;
 987
 988         double abs_error_r, abs_error_r2;
 989         double rel_error_r, rel_error_r2;
 990         double rel_error_r_old=0, rel_error_r2_old=0;
 991         double x0_r_error, x0_r2_error;
 992
 993
 994         /* Only set a Coulomb table for GB */
 995         /*
 996          tabsel[0]=etabGB;
 997          tabsel[1]=-1;
 998          tabsel[2]=-1;
 999         */
1000
1001         /* Set the table dimensions for GB, not really necessary to
1002          * use etiNR (since we only have one table, but ...)
1003          */
1004         snew(td,1);
1005         table.r         = fr->gbtabr;
1006         table.scale     = fr->gbtabscale;
1007         table.scale_exp = 0;
1008         table.n         = table.scale*table.r;
1009         nx0             = 0;
1010         nx              = table.scale*table.r;
1011
1012         /* Check whether we have to read or generate
1013          * We will always generate a table, so remove the read code
1014          * (Compare with original make_table function
1015          */
1016         bReadTab = FALSE;
1017         bGenTab  = TRUE;
1018
1019         /* Each table type (e.g. coul,lj6,lj12) requires four
1020          * numbers per datapoint. For performance reasons we want
1021          * the table data to be aligned to 16-byte. This is accomplished
1022          * by allocating 16 bytes extra to a temporary pointer, and then
1023          * calculating an aligned pointer. This new pointer must not be
1024          * used in a free() call, but thankfully we're sloppy enough not
1025          * to do this :-)
1026          */
1027
1028         /* 4 fp entries per table point, nx+1 points, and 16 bytes extra
1029            to align it. */
1030         p_tmp = malloc(4*(nx+1)*sizeof(real)+16);
1031
1032         /* align it - size_t has the same same as a pointer */
1033         table.tab = (real *) (((size_t) p_tmp + 16) & (~((size_t) 15)));
1034
1035         init_table(out,nx,nx0,table.scale,&(td[0]),!bReadTab);
1036
1037         /* Local implementation so we don't have to use the etabGB
1038          * enum above, which will cause problems later when
1039          * making the other tables (right now even though we are using
1040          * GB, the normal Coulomb tables will be created, but this
1041          * will cause a problem since fr->eeltype==etabGB which will not
1042          * be defined in fill_table and set_table_type
1043          */
1044
1045         for(i=nx0;i<nx;i++)
1046     {
1047                 Vtab    = 0.0;
1048                 Ftab    = 0.0;
1049                 r       = td->x[i];
1050                 r2      = r*r;
1051                 expterm = exp(-0.25*r2);
1052
1053                 Vtab = 1/sqrt(r2+expterm);
1054                 Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1055
1056                 /* Convert to single precision when we store to mem */
1057                 td->v[i]  = Vtab;
1058                 td->f[i]  = Ftab;
1059
1060     }
1061
1062         copy2table(table.n,0,4,td[0].x,td[0].v,td[0].f,table.tab);
1063
1064         if(bDebugMode())
1065     {
1066                 fp=xvgropen(fns[0],fns[0],"r","V",oenv);
1067                 /* plot the output 5 times denser than the table data */
1068                 /* for(i=5*nx0;i<5*table.n;i++) */
1069                 for(i=nx0;i<table.n;i++)
1070                 {
1071                         /* x0=i*table.r/(5*table.n); */
1072                         x0=i*table.r/table.n;
1073                         evaluate_table(table.tab,0,4,table.scale,x0,&y0,&yp);
1074                         fprintf(fp,"%15.10e  %15.10e  %15.10e\n",x0,y0,yp);
1075
1076                 }
1077                 ffclose(fp);
1078     }
1079
1080         /*
1081          for(i=100*nx0;i<99.81*table.n;i++)
1082          {
1083          r = i*table.r/(100*table.n);
1084          r2      = r*r;
1085          expterm = exp(-0.25*r2);
1086
1087          Vtab = 1/sqrt(r2+expterm);
1088          Ftab = (r-0.25*r*expterm)/((r2+expterm)*sqrt(r2+expterm));
1089
1090
1091          evaluate_table(table.tab,0,4,table.scale,r,&y0,&yp);
1092          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);
1093
1094          abs_error_r=fabs(y0-Vtab);
1095          abs_error_r2=fabs(yp-(-1)*Ftab);
1096
1097          rel_error_r=abs_error_r/y0;
1098          rel_error_r2=fabs(abs_error_r2/yp);
1099
1100
1101          if(rel_error_r>rel_error_r_old)
1102          {
1103          rel_error_r_old=rel_error_r;
1104          x0_r_error=x0;
1105          }
1106
1107          if(rel_error_r2>rel_error_r2_old)
1108          {
1109          rel_error_r2_old=rel_error_r2;
1110          x0_r2_error=x0;
1111          }
1112          }
1113
1114          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);
1115          printf("gb: XO_R=%g, X0_R2=%g\n",x0_r_error, x0_r2_error);
1116
1117          exit(1); */
1118         done_tabledata(&(td[0]));
1119         sfree(td);
1120
1121         return table;
1122
1123
1124 }
1125
1126 bondedtable_t make_bonded_table(FILE *fplog,char *fn,int angle)
1127 {
1128   t_tabledata td;
1129   double start;
1130   int    i;
1131   bondedtable_t tab;
1132
1133   if (angle < 2)
1134     start = 0;
1135   else
1136     start = -180.0;
1137   read_tables(fplog,fn,1,angle,&td);
1138   if (angle > 0) {
1139     /* Convert the table from degrees to radians */
1140     for(i=0; i<td.nx; i++) {
1141       td.x[i] *= DEG2RAD;
1142       td.f[i] *= RAD2DEG;
1143     }
1144     td.tabscale *= RAD2DEG;
1145   }
1146   tab.n = td.nx;
1147   tab.scale = td.tabscale;
1148   snew(tab.tab,tab.n*4);
1149   copy2table(tab.n,0,4,td.x,td.v,td.f,tab.tab);
1150   done_tabledata(&td);
1151
1152   return tab;
1153 }
1154
1155