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