src/mdlib/tables.c

   1 /*
   2  * $Id$
   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 #include <math.h>
  37 #include "typedefs.h"
  38 #include "names.h"
  39 #include "smalloc.h"
  40 #include "fatal.h"
  41 #include "futil.h"
  42 #include "xvgr.h"
  43 #include "vec.h"
  44 #include "main.h"
  45 #include "network.h"
  46
  47 /* All the possible (implemented) table functions */
  48 enum {
  49   etabLJ6,   etabLJ12, etabLJ6Shift, etabLJ12Shift, etabShift,
  50   etabRF,    etabCOUL, etabEwald, etabLJ6Switch, etabLJ12Switch,etabCOULSwitch,
  51   etabEXPMIN,etabUSER, etabNR
  52 };
  53 static const char *tabnm[etabNR] = {
  54   "LJ6",   "LJ12", "LJ6Shift", "LJ12Shift", "Shift",
  55   "RF",    "COUL", "Ewald", "LJ6Switch", "LJ12Switch","COULSwitch",
  56   "EXPMIN","USER"
  57 };
  58
  59 /* This flag tells whether this is a Coulomb type funtion */
  60 bool bCoulomb[etabNR] = { FALSE, FALSE, FALSE, FALSE, TRUE,
  61                           TRUE,  TRUE, TRUE, FALSE, FALSE, TRUE,
  62                           FALSE, FALSE };
  63
  64 /* Index in the table that says which function to use */
  65 enum { etiCOUL, etiLJ6, etiLJ12, etiNR };
  66
  67
  68 typedef struct {
  69   int  nx,nx0;
  70   double tabscale;
  71   double *x,*v,*v2;
  72 } t_tabledata;
  73
  74 #define pow2(x) ((x)*(x))
  75 #define pow3(x) ((x)*(x)*(x))
  76 #define pow4(x) ((x)*(x)*(x)*(x))
  77 #define pow5(x) ((x)*(x)*(x)*(x)*(x))
  78
  79 /* Calculate the potential and force for an r value
  80  * in exactly the same way it is done in the inner loop.
  81  * VFtab is a pointer to the table data, offset is
  82  * the point where we should begin and stride is
  83  * 4 if we have a buckingham table, 3 otherwise.
  84  * If you want to evaluate table no N, set offset to 4*N.
  85  *
  86  * We use normal precision here, since that is what we
  87  * will use in the inner loops.
  88  */
  89 static void evaluate_table(real VFtab[], int offset, int stride,
  90                            real tabscale, real r, real *y, real *yp)
  91 {
  92   int n;
  93   real rt,eps,eps2;
  94   real Y,F,Geps,Heps2,Fp;
  95
  96   rt       =  r*tabscale;
  97   n        =  (int)rt;
  98   eps      =  rt - n;
  99   eps2     =  eps*eps;
 100   n        =  offset+stride*n;
 101   Y        =  VFtab[n];
 102   F        =  VFtab[n+1];
 103   Geps     =  eps*VFtab[n+2];
 104   Heps2    =  eps2*VFtab[n+3];
 105   Fp       =  F+Geps+Heps2;
 106   *y       =  Y+eps*Fp;
 107   *yp      =  (Fp+Geps+2.0*Heps2)*tabscale;
 108 }
 109
 110
 111 static void splint(real xa[],real ya[],real y2a[],
 112                    int n,real x,real *y,real *yp)
 113 {
 114   int  klo,khi,k;
 115   real h,b,a,eps;
 116   real F,G,H;
 117
 118   klo=1;
 119   khi=n;
 120
 121   while ((khi-klo) > 1) {
 122     k=(khi+klo) >> 1;
 123     if (xa[k] > x)
 124       khi=k;
 125     else
 126       klo=k;
 127   }
 128   h = xa[khi]-xa[klo];
 129   if (h == 0.0)
 130     fatal_error(0,"Bad XA input to function splint");
 131   a   = (xa[khi]-x)/h;
 132   b   = (x-xa[klo])/h;
 133   *y  = a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
 134   *yp = (ya[khi]-ya[klo])/h+((3*a*a-1)*y2a[klo]-(3*b*b-1)*y2a[khi])*h/6.0;
 135
 136   eps = b;
 137   F   = (ya[khi]-ya[klo]-(h*h/6.0)*(2*y2a[klo]+y2a[khi]));
 138   G   = (h*h/2.0)*y2a[klo];
 139   H   = (h*h/6.0)*(y2a[khi]-y2a[klo]);
 140   *y  = ya[klo] + eps*F + eps*eps*G + eps*eps*eps*H;
 141   *yp = (F + 2*eps*G + 3*eps*eps*H)/h;
 142 }
 143
 144
 145 static void copy2table(int n,int offset,int stride,
 146                        double x[],double Vtab[],double Vtab2[],
 147                        real dest[],real r_zeros)
 148 {
 149 /* Use double prec. for the intermediary variables
 150  * and temporary x/vtab/vtab2 data to avoid unnecessary
 151  * loss of precision.
 152  */
 153   int  i,nn0;
 154   double F,G,H,h;
 155
 156   for(i=1; (i<n-1); i++) {
 157     h   = x[i+1]-x[i];
 158     F   = (Vtab[i+1]-Vtab[i]-(h*h/6.0)*(2*Vtab2[i]+Vtab2[i+1]));
 159     G   = (h*h/2.0)*Vtab2[i];
 160     H   = (h*h/6.0)*(Vtab2[i+1]-Vtab2[i]);
 161     nn0 = offset+i*stride;
 162     dest[nn0]   = Vtab[i];
 163     dest[nn0+1] = F;
 164     dest[nn0+2] = G;
 165     dest[nn0+3] = H;
 166   }
 167   if (r_zeros > 0.0) {
 168     for(i=1; (i<n-1); i++) {
 169       if (0.5*(x[i]+x[i+1]) >= r_zeros) {
 170         nn0 = offset+i*stride;
 171         dest[nn0]   = 0.0;
 172         dest[nn0+1] = 0.0;
 173         dest[nn0+2] = 0.0;
 174         dest[nn0+3] = 0.0;
 175       }
 176     }
 177   }
 178 }
 179
 180 static void init_table(FILE *fp,int n,int nx0,int tabsel,
 181                        real tabscale,t_tabledata *td,bool bAlloc)
 182 {
 183   int i;
 184
 185   td->nx  = n;
 186   td->nx0 = nx0;
 187   td->tabscale = tabscale;
 188   if (bAlloc) {
 189     snew(td->x,td->nx);
 190     snew(td->v,td->nx);
 191     snew(td->v2,td->nx);
 192   }
 193   for(i=td->nx0; (i<td->nx); i++)
 194     td->x[i] = i/tabscale;
 195 }
 196
 197 static void read_tables(FILE *fp,char *fn,t_tabledata td[])
 198 {
 199   char *libfn;
 200   real **yy=NULL;
 201   int  k,i,nx,nx0,ny,nny;
 202   bool bCont;
 203   real tabscale;
 204
 205   nny = 2*etiNR+1;
 206   libfn = low_libfn(fn,TRUE);
 207   nx  = read_xvg(libfn,&yy,&ny);
 208   if (ny != nny)
 209     fatal_error(0,"Trying to read file %s, but nr columns = %d, should be %d",
 210                 libfn,ny,nny);
 211   bCont = TRUE;
 212   for(nx0=0; bCont && (nx0 < nx); nx0++)
 213     for(k=1; (k<ny); k++)
 214       if (yy[k][nx0] != 0)
 215         bCont = FALSE;
 216   if (nx0 == nx)
 217     fatal_error(0,"All elements in table %s are zero!\n",libfn);
 218
 219   tabscale = (nx-1)/(yy[0][nx-1] - yy[0][0]);
 220   for(k=0; (k<etiNR); k++) {
 221     init_table(fp,nx,nx0,etabUSER,tabscale,&(td[k]),TRUE);
 222     for(i=0; (i<nx); i++) {
 223       td[k].x[i]  = yy[0][i];
 224       td[k].v[i]  = yy[2*k+1][i];
 225       td[k].v2[i] = yy[2*k+2][i];
 226     }
 227   }
 228   for(i=0; (i<ny); i++)
 229     sfree(yy[i]);
 230   sfree(yy);
 231
 232   if (fp)
 233     fprintf(fp,"Read user tables from %s with %d data points.\n"
 234             "Tabscale = %g points/nm\n",libfn,nx,tabscale);
 235 }
 236
 237 static void done_tabledata(t_tabledata *td)
 238 {
 239   int i;
 240
 241   if (!td)
 242     return;
 243
 244   sfree(td->x);
 245   sfree(td->v);
 246   sfree(td->v2);
 247 }
 248
 249 static void fill_table(t_tabledata *td,int tp,t_forcerec *fr)
 250 {
 251   /* Fill the table according to the formulas in the manual.
 252    * In principle, we only need the potential and the second
 253    * derivative, but then we would have to do lots of calculations
 254    * in the inner loop. By precalculating some terms (see manual)
 255    * we get better eventual performance, despite a larger table.
 256    *
 257    * Since some of these higher-order terms are very small,
 258    * we always use double precision to calculate them here, in order
 259    * to avoid unnecessary loss of precision.
 260    */
 261 #ifdef DEBUG_SWITCH
 262   FILE *fp;
 263 #endif
 264   int  i,p;
 265   double r1,rc,r12,r13;
 266   double r,r2,r6;
 267   double expr,Vtab,Ftab,Vtab2,Ftab2;
 268   /* Parameters for David's function */
 269   double A=0,B=0,C=0,A_3=0,B_4=0;
 270   /* Parameters for the switching function */
 271   double ksw,swi,swi1,swi2;
 272   /* Temporary parameters */
 273   bool bSwitch,bShift;
 274   double VtabT;
 275   double VtabT1;
 276   double VtabT2;
 277   double ewc=fr->ewaldcoeff;
 278   double isp= 0.564189583547756;
 279
 280   bSwitch = ((tp == etabLJ6Switch)    || (tp == etabLJ12Switch)    ||
 281              (tp == etabCOULSwitch));
 282   bShift  = ((tp == etabLJ6Shift) || (tp == etabLJ12Shift) ||
 283              (tp == etabShift));
 284
 285   if (bCoulomb[tp]) {
 286     r1 = fr->rcoulomb_switch;
 287     rc = fr->rcoulomb;
 288   }
 289   else {
 290     r1 = fr->rvdw_switch;
 291     rc = fr->rvdw;
 292   }
 293   if (bSwitch)
 294     ksw  = 1.0/(pow5(rc-r1));
 295   else
 296     ksw  = 0.0;
 297   if (bShift) {
 298     if (tp == etabShift)
 299       p=1;
 300     else if (tp == etabLJ6Shift)
 301       p=6;
 302     else
 303       p=12;
 304
 305     A = p * ((p+1)*r1-(p+4)*rc)/(pow(rc,p+2)*pow2(rc-r1));
 306     B = -p * ((p+1)*r1-(p+3)*rc)/(pow(rc,p+2)*pow3(rc-r1));
 307     C = 1.0/pow(rc,p)-A/3.0*pow3(rc-r1)-B/4.0*pow4(rc-r1);
 308     if (tp == etabLJ6Shift) {
 309       A=-A;
 310       B=-B;
 311       C=-C;
 312     }
 313     A_3=A/3.0;
 314     B_4=B/4.0;
 315   }
 316   if (debug) { fprintf(debug,"Further\n"); fflush(debug); }
 317
 318 #ifdef DEBUG_SWITCH
 319   fp=xvgropen("switch.xvg","switch","r","s");
 320 #endif
 321   for(i=td->nx0; (i<td->nx); i++) {
 322     r     = td->x[i];
 323     r2    = r*r;
 324     r6    = 1.0/(r2*r2*r2);
 325     r12   = r6*r6;
 326     Vtab  = 0.0;
 327     Ftab  = 0.0;
 328     Vtab2 = 0.0;
 329     Ftab2 = 0.0;
 330     if (bSwitch) {
 331       /* swi is function, swi1 1st derivative and swi2 2nd derivative */
 332       /* The switch function is 1 for r<r1, 0 for r>rc, and smooth for
 333        * r1<=r<=rc. The 1st and 2nd derivatives are both zero at
 334        * r1 and rc.
 335        * ksw is just the constant 1/(rc-r1)^5, to save some calculations...
 336        */
 337       if(r<=r1) {
 338         swi = 1.0;
 339         swi1 = swi2 = 0.0;
 340       } else if (r>=rc) {
 341         swi = swi1 = swi2 = 0.0;
 342       } else {
 343         swi      = 1 - 10*pow3(r-r1)*ksw*pow2(rc-r1)
 344           + 15*pow4(r-r1)*ksw*(rc-r1) - 6*pow5(r-r1)*ksw;
 345         swi1     = -30*pow2(r-r1)*ksw*pow2(rc-r1)
 346           + 60*pow3(r-r1)*ksw*(rc-r1) - 30*pow4(r-r1)*ksw;
 347         swi2     =  -60*(r-r1)*ksw*pow2(rc-r1)
 348           + 180*pow2(r-r1)*ksw*(rc-r1) - 120*pow3(r-r1)*ksw;
 349       }
 350     }
 351     else { /* not really needed, but avoids compiler warnings... */
 352       swi = 1.0;
 353       swi1 = swi2 = 0.0;
 354     }
 355 #ifdef DEBUG_SWITCH
 356     fprintf(fp,"%10g  %10g  %10g  %10g\n",r,swi,swi1,swi2);
 357 #endif
 358     switch (tp) {
 359     case etabLJ6:
 360       /* Dispersion */
 361       Vtab  = -r6;
 362       Ftab  = 6.0*Vtab/r;
 363       Vtab2 = 7.0*Ftab/r;
 364       Ftab2 = 8.0*Vtab2/r;
 365       break;
 366     case etabLJ6Switch:
 367     case etabLJ6Shift:
 368       /* Dispersion */
 369       if (r < rc) {
 370         Vtab  = -r6;
 371         Ftab  = 6.0*Vtab/r;
 372         Vtab2 = 7.0*Ftab/r;
 373         Ftab2 = 8.0*Vtab2/r;
 374       }
 375       break;
 376     case etabLJ12:
 377       /* Repulsion */
 378       Vtab  = r12;
 379       Ftab  = 12.0*Vtab/r;
 380       Vtab2 = 13.0*Ftab/r;
 381       Ftab2 = 14.0*Vtab2/r;
 382       break;
 383     case etabLJ12Switch:
 384     case etabLJ12Shift:
 385       /* Repulsion */
 386       if (r < rc) {
 387         Vtab  = r12;
 388         Ftab  = 12.0*Vtab/r;
 389         Vtab2 = 13.0*Ftab/r;
 390         Ftab2 = 14.0*Vtab2/r;
 391       }
 392       break;
 393     case etabCOUL:
 394       Vtab  = 1.0/r;
 395       Ftab  = 1.0/r2;
 396       Vtab2 = 2.0/(r*r2);
 397       Ftab2 = 6.0/(r2*r2);
 398       break;
 399     case etabCOULSwitch:
 400     case etabShift:
 401       if (r < rc) {
 402         Vtab  = 1.0/r;
 403         Ftab  = 1.0/r2;
 404         Vtab2 = 2.0/(r*r2);
 405         Ftab2 = 6.0/(r2*r2);
 406       }
 407       break;
 408     case etabEwald:
 409       Vtab  = erfc(ewc*r)/r;
 410       Ftab  = erfc(ewc*r)/r2+2*exp(-(ewc*ewc*r2))*ewc*isp/r;
 411       Vtab2 = 2*erfc(ewc*r)/(r*r2)+4*exp(-(ewc*ewc*r2))*ewc*isp/r2+
 412           4*ewc*ewc*ewc*exp(-(ewc*ewc*r2))*isp;
 413       Ftab2 = 6*erfc(ewc*r)/(r2*r2)+
 414           12*exp(-(ewc*ewc*r2))*ewc*isp/(r*r2)+
 415           8*ewc*ewc*ewc*exp(-(ewc*ewc*r2))*isp/r+
 416           8*ewc*ewc*ewc*ewc*ewc*r*exp(-(ewc*ewc*r2))*isp;
 417       break;
 418     case etabRF:
 419       Vtab  = 1.0/r      +   fr->k_rf*r2 - fr->c_rf;
 420       Ftab  = 1.0/r2     - 2*fr->k_rf*r;
 421       Vtab2 = 2.0/(r*r2) + 2*fr->k_rf;
 422       Ftab2 = 6.0/(r2*r2);
 423       break;
 424     case etabEXPMIN:
 425       expr  = exp(-r);
 426       Vtab  = expr;
 427       Ftab  = expr;
 428       Vtab2 = expr;
 429       Ftab2 = expr;
 430       break;
 431     default:
 432       fatal_error(0,"Table type %d not implemented yet. (%s,%d)",
 433                   tp,__FILE__,__LINE__);
 434     }
 435     if (bShift) {
 436       /* Normal coulomb with cut-off correction for potential */
 437       if (r < rc) {
 438         Vtab -= C;
 439         /* If in Shifting range add something to it */
 440         if (r > r1) {
 441           r12 = (r-r1)*(r-r1);
 442           r13 = (r-r1)*r12;
 443           Vtab  += - A_3*r13 - B_4*r12*r12;
 444           Ftab  +=   A*r12 + B*r13;
 445           Vtab2 += - 2.0*A*(r-r1) - 3.0*B*r12;
 446           Ftab2 +=   2.0*A + 6.0*B*(r-r1);
 447         }
 448       }
 449     }
 450
 451     if ((r > r1) && bSwitch) {
 452       VtabT     = Vtab;
 453       VtabT1    = -Ftab;
 454       VtabT2    = Vtab2;
 455       Vtab   = VtabT*swi;
 456       Vtab2  = VtabT2*swi + VtabT1*swi1 + VtabT1*swi1 + VtabT*swi2;
 457     }
 458
 459     /* Convert to single precision when we store to mem */
 460     td->v[i]  = Vtab;
 461     td->v2[i] = Vtab2;
 462   }
 463
 464 #ifdef DEBUG_SWITCH
 465   fclose(fp);
 466 #endif
 467 }
 468
 469 static void set_table_type(int tabsel[],t_forcerec *fr)
 470 {
 471   /* Set the different table indices.
 472    * Coulomb first.
 473    */
 474
 475   switch (fr->eeltype) {
 476   case eelCUT:
 477     tabsel[etiCOUL] = etabCOUL;
 478     break;
 479   case eelPPPM:
 480   case eelPOISSON:
 481     if ((fr->rcoulomb > fr->rcoulomb_switch) && fr->bcoultab)
 482         tabsel[etiCOUL] = etabShift;
 483     else
 484         tabsel[etiCOUL] = etabCOUL;  /* 1-4 */
 485     break;
 486   case eelSHIFT:
 487     if (fr->rcoulomb > fr->rcoulomb_switch)
 488       tabsel[etiCOUL] = etabShift;
 489     else
 490       tabsel[etiCOUL] = etabCOUL;
 491     break;
 492   case eelEWALD:
 493   case eelPME:
 494       if(fr->bcoultab)
 495           tabsel[etiCOUL] = etabEwald;
 496       else
 497           tabsel[etiCOUL] = etabCOUL; /* 1-4 */
 498       break;
 499   case eelRF:
 500   case eelGRF:
 501     tabsel[etiCOUL] = etabRF;
 502     break;
 503   case eelSWITCH:
 504     tabsel[etiCOUL] = etabCOULSwitch;
 505     break;
 506   case eelUSER:
 507     tabsel[etiCOUL] = etabUSER;
 508     break;
 509   default:
 510     fatal_error(0,"Invalid eeltype %d in %s line %d",fr->eeltype,
 511                 __FILE__,__LINE__);
 512   }
 513
 514   /* Van der Waals time */
 515   if (fr->bBHAM) {
 516     tabsel[etiLJ6]  = etabLJ6;
 517     tabsel[etiLJ12] = etabEXPMIN;
 518   }
 519   else {
 520     switch (fr->vdwtype) {
 521     case evdwSWITCH:
 522       tabsel[etiLJ6]  = etabLJ6Switch;
 523       tabsel[etiLJ12] = etabLJ12Switch;
 524       break;
 525     case evdwSHIFT:
 526       tabsel[etiLJ6]  = etabLJ6Shift;
 527       tabsel[etiLJ12] = etabLJ12Shift;
 528       break;
 529     case evdwUSER:
 530       tabsel[etiLJ6]  = etabUSER;
 531       tabsel[etiLJ12] = etabUSER;
 532       break;
 533     case evdwCUT:
 534       tabsel[etiLJ6]  = etabLJ6;
 535       tabsel[etiLJ12] = etabLJ12;
 536       break;
 537     default:
 538       fatal_error(0,"Invalid vdwtype %d in %s line %d",fr->vdwtype,
 539                   __FILE__,__LINE__);
 540     }
 541   }
 542 }
 543
 544 void make_tables(FILE *out,t_forcerec *fr,bool bVerbose,char *fn)
 545 {
 546   const char *fns[3] = { "ctab.xvg", "dtab.xvg", "rtab.xvg" };
 547   FILE        *fp;
 548   t_tabledata *td;
 549   bool        bReadTab,bGenTab;
 550   real        x0,y0,yp,rtab;
 551   int         i,j,k,nx,nx0,tabsel[etiNR];
 552
 553   set_table_type(tabsel,fr);
 554   snew(td,etiNR);
 555   fr->tabscale = 0;
 556   rtab         = fr->rtab;
 557   nx0          = 10;
 558   nx           = 0;
 559
 560   /* Check whether we have to read or generate */
 561   bReadTab = FALSE;
 562   bGenTab  = FALSE;
 563   for(i=0; (i<etiNR); i++) {
 564     if (tabsel[i] == etabUSER)
 565       bReadTab = TRUE;
 566     else
 567       bGenTab  = TRUE;
 568   }
 569   if (bReadTab) {
 570     read_tables(out,fn,td);
 571     fr->tabscale = td[0].tabscale;
 572     fr->rtab     = td[0].x[td[0].nx-1];
 573     nx0          = td[0].nx0;
 574     nx           = fr->ntab = fr->rtab*fr->tabscale;
 575     if (fr->rtab < rtab)
 576       fatal_error(0,"Tables in file %s not long enough for cut-off:\n"
 577                   "\tshould be at least %f nm\n",fn,rtab);
 578   }
 579   if (bGenTab) {
 580     if (!bReadTab) {
 581 #ifdef DOUBLE
 582       fr->tabscale = 2000.0;
 583 #else
 584       fr->tabscale = 500.0;
 585 #endif
 586       nx = fr->ntab = fr->rtab*fr->tabscale;
 587     }
 588   }
 589   /* Each table type (e.g. coul,lj6,lj12) requires four
 590    * numbers per datapoint. The exponential Buckingham table
 591    * requires an extra one for accuracy. Since we still have
 592    * the coulomb and dispersion, there will be four table
 593    * types with four elements each when we use Buckingham.
 594    */
 595
 596   snew(fr->coulvdwtab,( fr->bBHAM ? 16 :12 )*(nx+1)+1);
 597   for(k=0; (k<etiNR); k++) {
 598     if (tabsel[k] != etabUSER) {
 599       init_table(out,nx,nx0,tabsel[k],
 600                  (tabsel[k] == etabEXPMIN) ? fr->tabscale_exp : fr->tabscale,
 601                  &(td[k]),!bReadTab);
 602       fill_table(&(td[k]),tabsel[k],fr);
 603       if (out)
 604         fprintf(out,"Generated table with %d data points for %s.\n"
 605                 "Tabscale = %g points/nm\n",
 606                 td[k].nx,tabnm[tabsel[k]],td[k].tabscale);
 607
 608     }
 609     copy2table(td[k].nx,k*4,12,td[k].x,td[k].v,td[k].v2,fr->coulvdwtab,-1);
 610
 611     if (bDebugMode() && bVerbose) {
 612       fp=xvgropen(fns[k],fns[k],"r","V");
 613       /* plot the output 5 times denser than the table data */
 614       for(i=5*nx0;i<5*fr->ntab;i++) {
 615         x0=i*fr->rtab/(5*fr->ntab);
 616         evaluate_table(fr->coulvdwtab,4*k,12,fr->tabscale,x0,&y0,&yp);
 617         fprintf(fp,"%15.10e  %15.10e  %15.10e\n",x0,y0,yp);
 618       }
 619       ffclose(fp);
 620     }
 621     done_tabledata(&(td[k]));
 622   }
 623   sfree(td);
 624 }
 625