src/contrib/gen_table.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.3.99_development_20071104
  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-2006, 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  * Groningen Machine for Chemical Simulation
  34  */
  35 #include <math.h>
  36 #include <string.h>
  37 #include <stdio.h>
  38
  39 #include "copyrite.h"
  40 #include "typedefs.h"
  41 #include "macros.h"
  42 #include "vec.h"
  43 #include "gromacs/commandline/pargs.h"
  44 #include "coulomb.h"
  45
  46 enum { mGuillot2001a, mAB1, mLjc, mMaaren, mGuillot_Maple, mHard_Wall, mGG, mGG_qd_q, mGG_qd_qd, mGG_q_q, mNR };
  47
  48 static double erf2(double x)
  49 {
  50   return -(2*x*M_2_SQRTPI)*exp(-x*x);
  51 }
  52
  53 static double erf1(double x)
  54 {
  55   return M_2_SQRTPI*exp(-x*x);
  56 }
  57
  58 static void do_hard(FILE *fp,int pts_nm,double efac,double delta)
  59 {
  60   int    i,k,imax;
  61   double x,vr,vr2,vc,vc2;
  62
  63   if (delta < 0)
  64     gmx_fatal(FARGS,"Delta should be >= 0 rather than %f\n",delta);
  65
  66   imax     = 3.0*pts_nm;
  67   for(i=0; (i<=imax); i++) {
  68     x   =  i*(1.0/pts_nm);
  69
  70     if (x < delta) {
  71       /* Avoid very high numbers */
  72       vc = vc2 = 1/delta;
  73     }
  74     else {
  75       vc  = 1/(x);
  76       vc2 = 2/pow(x,3);
  77     }
  78     vr  = erfc(efac*(x-delta))/2;
  79     vr2 = (1-erf2(efac*(x-delta)))/2;
  80     fprintf(fp,"%12.5e  %12.5e  %12.5e  %12.5e  %12.5e  %12.5e  %12.5e\n",
  81             x,vr,vr2,0.0,0.0,vc,vc2);
  82   }
  83
  84 }
  85
  86 static void do_AB1(FILE *fp,int eel,int pts_nm,int ndisp,int nrep)
  87 {
  88   int    i,k,imax;
  89   double myfac[3] = { 1, -1, 1 };
  90   double myexp[3] = { 1, 6, 0 };
  91   double x,v,v2;
  92
  93   myexp[1] = ndisp;
  94   myexp[2] = nrep;
  95   imax     = 3.0*pts_nm;
  96   for(i=0; (i<=imax); i++) {
  97     x   =  i*(1.0/pts_nm);
  98
  99     fprintf(fp,"%12.5e",x);
 100
 101     for(k=0; (k<3); k++) {
 102       if (x < 0.04) {
 103         /* Avoid very high numbers */
 104         v = v2 = 0;
 105       }
 106       else {
 107         v  =  myfac[k]*pow(x,-myexp[k]);
 108         v2 = (myexp[k]+1)*(myexp[k])*v/(x*x);
 109       }
 110       fprintf(fp,"  %12.5e  %12.5e",v,v2);
 111     }
 112     fprintf(fp,"\n");
 113   }
 114 }
 115
 116 static void lo_do_ljc(double r,
 117                       double *vc,double *fc,
 118                       double *vd,double *fd,
 119                       double *vr,double *fr)
 120 {
 121   double r2,r_6,r_12;
 122
 123   r2    = r*r;
 124   r_6   = 1.0/(r2*r2*r2);
 125   r_12  = r_6*r_6;
 126
 127   *vc   = 1.0/r;            /*  f(x)     Coulomb    */
 128   *fc   = 1.0/(r2);         /* -f'(x)               */
 129
 130   *vd   = -r_6;             /*  g(c)     Dispersion */
 131   *fd   =  6.0*(*vd)/r;     /* -g'(x)               */
 132
 133   *vr   = r_12;             /*  h(x)     Repulsion  */
 134   *fr   = 12.0*(*vr)/r;     /* -h'(x)               */
 135 }
 136
 137 /* use with coulombtype = user */
 138 static void lo_do_ljc_pme(double r,
 139                           double rcoulomb, double ewald_rtol,
 140                           double *vc,double *fc,
 141                           double *vd,double *fd,
 142                           double *vr,double *fr)
 143 {
 144   double r2,r_6,r_12;
 145   double ewc;
 146
 147   ewc  = calc_ewaldcoeff(rcoulomb,ewald_rtol);
 148
 149   r2   = r*r;
 150   r_6  = 1.0/(r2*r2*r2);
 151   r_12 = r_6*r_6;
 152
 153   *vc   = erfc(ewc*r)/r;
 154   /* *vc2  = 2*erfc(ewc*r)/(r*r2)+2*exp(-(ewc*ewc*r2))*ewc*M_2_SQRTPI/r2+
 155      2*ewc*ewc*ewc*exp(-(ewc*ewc*r2))*M_2_SQRTPI;*/
 156   *fc  = ewc*exp(-ewc*ewc*r2)*M_2_SQRTPI/r + erfc(ewc*r)/r2;
 157
 158   *vd  = -r_6;
 159   *fd  = -6.0*(*vd)/r;
 160
 161   *vr  = r_12;
 162   *fr  = 12.0*(*vr)/r;
 163 }
 164
 165 static void lo_do_guillot(double r,double xi, double xir,
 166                           double *vc,double *fc,
 167                           double *vd,double *fd,
 168                           double *vr,double *fr)
 169 {
 170   double qO     = -0.888;
 171   double qOd    =  0.226;
 172   double f0     = qOd/qO;
 173   double sqpi   = sqrt(M_PI);
 174   double rxi1,rxi2,z;
 175   double r2,r_6;
 176
 177   r2   = r*r;
 178   r_6  = 1.0/(r2*r2*r2);
 179
 180   rxi1    = r/(2*xi);
 181   rxi2    = r/(sqrt(2)*xi);
 182   *vc   = (1 + f0*f0*erf(r/(2*xi)) + 2*f0*erf(r/(sqrt(2)*xi)) )/r;
 183
 184   *fc   =  f0*f0*erf(r/(2*xi)) + 2*f0*erf(r/(sqrt(2)*xi));
 185     ;
 186  /* MuPad: Uc := erf(r/(2*xi))/r +
 187
 188      Mathematica:
 189      r1 := r/(2*xi);
 190      r2 := r/(Sqrt[2] * xi);
 191      Uc[r_] := (1 + f0 * f0 * Erf[r/(2*xi)] + 2 * f0 * Erf[r/(Sqrt[2]*xi)]) / r;
 192      -D[Uc[r],r]
 193      CForm=
 194      -(((2*f0*Sqrt(2/Pi))/(Power(E,Power(r,2)/(2.*Power(xi,2)))*xi) +
 195      Power(f0,2)/(Power(E,Power(r,2)/(4.*Power(xi,2)))*Sqrt(Pi)*xi))/r) +
 196      (1 + Power(f0,2)*Erf(r/(2.*xi)) + 2*f0*Erf(r/(Sqrt(2)*xi)))/Power(r,2)
 197
 198
 199 Uc1[r_] := 1/r;
 200 -D[Uc1[r],r]
 201           -2
 202 Out[20]= r
 203
 204 Uc2[r_] := f0^2 * Erf[r1] / r;
 205 -D[Uc2[r],r]
 206
 207
 208 Uc3[r_] := 2 * f0 * Erf[r2]/ r;
 209 -D[Uc3[r],r]
 210
 211 Uc[r_] := Erf[r/(2*xi)] / r
 212
 213 D[Uc[r],r]
 214
 215
 216 D[Erf[r],r]
 217
 218 */
 219     *vc   = (1 + sqr(f0)*erf(rxi1) + 2*f0*erf(rxi2))/r;
 220     *fc   =
 221       (1/r
 222         + (- f0 * (2 * sqrt(2) + exp(r2/4*xi*xi)*f0)/(exp(r2/(2*xi*xi))*sqrt(M_PI)*xi) + f0*f0*erf(r/(2*xi)) + 2 *f0 * erf(r/(sqrt(2 * xi)))  )/r2)
 223       ;
 224
 225
 226   /*  *vc2  = ((2/sqr(r))*(*vc -
 227                        sqr(f0)*erf1(r1)/(2*xi) -
 228                        4*f0*erf1(r2)/sqrt(2)*xi) +
 229                        (1/r)*(sqr(f0/(2.0*xi))*erf2(r1) + (2*f0/sqr(xi)))*erf2(r2)); */
 230
 231   *vd  = -r_6;
 232   *fd  = -6.0*(*vd)/r;
 233
 234   z     = r/(2.0*xir);
 235   *vr   = erfc(z)/z;
 236   *fr   = 0.0;
 237   //  *vr2  = (sqpi*(*vr)/(2.0*z*z)+(1.0/(z*z)+1)*exp(-z*z))/(sqpi*sqr(xir));
 238 }
 239
 240 void lo_do_guillot_maple(double r,double xi,double xir,
 241                          double *vc,double *vc2,
 242                          double *vd,double *vd2,
 243                          double *vr,double *vr2)
 244 {
 245   double qO     = -0.888;
 246   double qOd    = 0.226;
 247   double f0     = qOd/qO;
 248   double sqpi   = sqrt(M_PI);
 249
 250   *vc = pow(-f0/(1.0+f0)+1.0,2.0)/r+pow(-f0/(1.0+f0)+1.0,2.0)*f0*f0*erf(r/xi/2.0)/r+2.0*pow(-f0/(1.0+f0)+1.0,2.0)*f0*erf(r*sqrt(2.0)/xi/2.0)/r;
 251   *vc2 = 2.0*pow(-f0/(1.0+f0)+1.0,2.0)/(r*r*r)-pow(-f0/(1.0+f0)+1.0,2.0)*f0*f0/sqrt(M_PI)/(xi*xi*xi)*exp(-r*r/(xi*xi)/4.0)/2.0-2.0*pow(-f0/(1.0+f0)+1.0,2.0)*f0*f0/sqrt(M_PI)*exp(-r*r/(xi*xi)/4.0)/xi/(r*r)+2.0*pow(-f0/(1.0+f0)+1.0,2.0)*f0*f0*erf(r/xi/2.0)/(r*r*r)-2.0*pow(-f0/(1.0+f0)+1.0,2.0)*f0/sqrt(M_PI)/(xi*xi*xi)*exp(-r*r/(xi*xi)/2.0)*sqrt(2.0)-4.0*pow(-f0/(1.0+f0)+1.0,2.0)*f0/sqrt(M_PI)*exp(-r*r/(xi*xi)/2.0)*sqrt(2.0)/xi/(r*r)+4.0*pow(-f0/(1.0+f0)+1.0,2.0)*f0*erf(r*sqrt(2.0)/xi/2.0)/(r*r*r);
 252
 253   *vd   = -1.0/(r*r*r*r*r*r);
 254   *vd2  = -42.0/(r*r*r*r*r*r*r*r);
 255   *vr   = 2.0*erfc(r/xir/2.0)/r*xir;
 256   *vr2  = 1.0/sqrt(M_PI)/(xir*xir)*exp(-r*r/(xir*xir)/4.0)+4.0/sqrt(M_PI)*exp(-r*r/(xir*xir)/4.0)/(r*r)+4.0*erfc(r/xir/2.0)/(r*r*r)*xir;
 257 }
 258
 259 static void lo_do_GG(double r,double xi,double xir,
 260                      double *vc,double *fc,
 261                      double *vd,double *fd,
 262                      double *vr,double *fr)
 263 {
 264   double qO     = -0.888;
 265   double qOd    =  0.226;
 266   double f0     = qOd/qO;
 267   double sqpi   = sqrt(M_PI);
 268   double r2,xi2;
 269
 270   r2 = r*r;
 271   xi2 = xi*xi;
 272
 273   *vc = 1.0/r + f0*f0*erf(r/(2*xi))/r + 2*f0*erf(r/(sqrt(2)*xi))/r;
 274
 275   // -D[1/r,r] -D[f0*f0*Erf[r/(2*xi)]/r,r] -D[2*f0*Erf[r/(Sqrt[2]*xi)]/r,r]
 276   *fc  = (
 277     1.0/r2 +
 278     f0*f0*(-exp(-r2/(4*xi2))/(sqrt(M_PI) * r * xi) + erf(r/(2*xi))/r2) +
 279     2*f0*(-sqrt(2.0/M_PI)*exp(-r2/(2*xi2))/ (r*xi) + erf(r/(sqrt(2)*xi))/r2)
 280     );
 281
 282   // -D[1/r^6,r]
 283   *vd  = -1.0/(r*r*r*r*r*r);
 284   *fd  = 6.0*(*vd)/r;
 285
 286   //  -D[2*xir*Erfc[r/(2*xir)]/r,r]
 287   *vr  = 2.*xir*erfc(r/(2.*xir))/r;
 288   *fr  = -(-2.*exp(-r2/(4*xir*xir)) / (sqrt(M_PI)*r)  - 2*xir*erfc(r/(2*xir))/r2  );
 289
 290 }
 291
 292 /* Guillot2001 diffuse charge - diffuse charge interaction
 293    Mathematica
 294
 295 In[19]:= Uc[r_] := Erf[r/(2*xi)]/r
 296
 297 In[20]:= -D[Uc[r],r]
 298
 299                                              r
 300                                         Erf[----]
 301                        1                    2 xi
 302 Out[20]= -(-------------------------) + ---------
 303              2      2                       2
 304             r /(4 xi )                     r
 305            E           Sqrt[Pi] r xi
 306 */
 307 void lo_do_GG_qd_qd(double r,double xi,double xir,
 308                     double *vc,double *fc,
 309                     double *vd,double *fd,
 310                     double *vr,double *fr)
 311 {
 312   double sqpi   = sqrt(M_PI);
 313
 314   *vc = erf(r/(2*xi))/r;
 315     //erf((r*(1.0/2.0))/xi)/r;
 316   *fc = -(1.0/(exp(r*r/(4*xi*xi))*sqpi*r*xi)) + (erf(r/2*xi)/(r*r));
 317
 318     //2.0*pow(r, -3.0)*erf((r*(1.0/2.0))/xi) - (1.0/2.0)*pow(M_PI, -1.0/2.0)*pow(xi, -3.0)*exp((-1.0/4.0)*(r*r)*pow(xi, -2.0)) - (2.0*pow(M_PI, -1.0/2.0)*pow(r, -2.0)*exp((-1.0/4.0)*(r*r)*pow(xi, -2.0)))/xi ;
 319   *vd  = 0.0;
 320   *fd  = 0.0;
 321   *vr  = 0.0;
 322   *fr  = 0.0;
 323 }
 324
 325 /* Guillot2001 charge - diffuse charge interaction eqn 4 & 5
 326    Mathematica
 327 In[17]:= Uc[r_] := Erf[r/(Sqrt[2]*xi)]/r
 328
 329 In[18]:= -D[Uc[r],r]
 330
 331                     2                  r
 332                Sqrt[--]        Erf[----------]
 333                     Pi             Sqrt[2] xi
 334 Out[18]= -(----------------) + ---------------
 335              2      2                 2
 336             r /(2 xi )               r
 337            E           r xi
 338 */
 339 void lo_do_GG_q_qd(double r,double xi,double xir,
 340                    double *vc,double *fc,
 341                    double *vd,double *fd,
 342                    double *vr,double *fr)
 343 {
 344   double sqpi   = sqrt(M_PI);
 345
 346   *vc = erf(r/(sqrt(2)*xi)) / r;
 347     //erf(((1.0/2.0)*pow(2.0, 1.0/2.0)*r)/xi)/r ;
 348   *fc = -(sqrt(2/M_PI)/(exp(r*r/(2*xi*xi))*r*xi)) + (erf(r/(sqrt(2)*xi))/(r*r));
 349     //2.0*pow(r, -3.0)*erf(((1.0/2.0)*pow(2.0, 1.0/2.0)*r)/xi) - pow(2.0, 1.0/2.0)*pow(M_PI, -1.0/2.0)*pow(xi, -3.0)*exp((-1.0/2.0)*(r*r)*pow(xi, -2.0)) - (2.0*pow(2.0, 1.0/2.0)*pow(M_PI, -1.0/2.0)*pow(r, -2.0)*exp((-1.0/2.0)*(r*r)*pow(xi, -2.0)))/xi ;
 350
 351   *vd  = 0.0;
 352   *fd  = 0.0;
 353   *vr  = 0.0;
 354   *fr  = 0.0;
 355 }
 356
 357 /* Guillot2001 charge - charge interaction (normal coulomb), repulsion and dispersion
 358    Mathematica
 359
 360 In[6]:= Uc[r_] := 1.0/r
 361
 362 In[7]:= -D[Uc[r],r]
 363
 364         1.
 365 Out[7]= --
 366          2
 367         r
 368
 369 In[8]:= Ud[r_] := -1.0/r^6
 370
 371 In[9]:= -D[Ud[r],r]
 372
 373         -6.
 374 Out[9]= ---
 375          7
 376         r
 377
 378 In[13]:= Ur[r_] := (2*xir)*Erfc[r/(2*xir)]/r
 379
 380 In[14]:= -D[Ur[r],r]
 381                                                 r
 382                                    2 xir Erfc[-----]
 383                     2                         2 xir
 384 Out[16]= ----------------------- + -----------------
 385            2       2                       2
 386           r /(4 xir )                     r
 387          E            Sqrt[Pi] r
 388
 389
 390 */
 391 void lo_do_GG_q_q(double r,double xi,double xir,
 392                   double *vc,double *fc,
 393                   double *vd,double *fd,
 394                   double *vr,double *fr)
 395 {
 396   double sqpi   = sqrt(M_PI);
 397
 398   *vc  = 1.0/r;
 399   *fc  = 1.0/(r*r);
 400
 401   *vd  = -1.0/(r*r*r*r*r*r);
 402   *fd  = -6.0/(r*r*r*r*r*r*r);
 403
 404   *vr  = (2.0*xir*erfc(r/(2.0*xir)))/r;
 405   *fr  = 2.0/(exp((r*r)/(4*xir*xir)) * sqpi *r) + (2*xir*erfc((r*xir)/2.0))/(r*r);
 406     //4.0*pow(M_PI, -1.0/2.0)*pow(r, -2.0)*exp((-1.0/4.0)*(r*r)*pow(xir, -2.0)) + pow(M_PI, -1.0/2.0)*pow(xir, -2.0)*exp((-1.0/4.0)*(r*r)*pow(xir, -2.0)) + 4.0*pow(r, -3.0)*xir*erfc((r*(1.0/2.0))/xir);
 407 }
 408
 409 static void do_guillot(FILE *fp,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 410 {
 411   int    i,i0,imax;
 412   double r,vc,fc,vd,fd,vr,fr;
 413
 414   imax = 3*pts_nm;
 415   for(i=0; (i<=imax); i++) {
 416     r     = i*(1.0/pts_nm);
 417     /* Avoid very large numbers */
 418     if (r < 0.04) {
 419       vc = fc = vd = fd = vr = fr = 0;
 420     }
 421     else
 422       if (eel == eelPME) {
 423         fprintf(fp, "Not implemented\n");
 424       } else if (eel == eelCUT) {
 425         lo_do_guillot(r,xi,xir,&vc,&fc,&vd,&fd,&vr,&fr);
 426       }
 427     fprintf(fp,"%15.10e   %15.10e %15.10e   %15.10e %15.10e   %15.10e %15.10e\n",
 428             r,vc,fc,vd,fd,vr,fr);
 429
 430   }
 431 }
 432
 433 /* TODO:
 434    PvM: Everything is hardcoded, we should fix that. How?
 435 */
 436 static void do_guillot2001a(const char *file,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 437 {
 438   FILE *fp=NULL;
 439   static char buf[256];
 440   static char *atype[]   = { "HW", "OW", "HWd", "OWd", NULL };
 441   int    i,j,k,i0,imax,atypemax=4;
 442   double r,vc,fc,vd,fd,vr,fr;
 443
 444   /* For Guillot2001a we have four types: HW, OW, HWd and OWd. */
 445
 446   for (j=0;(j<atypemax);j++) {           /* loops over types */
 447     for (k=0; (k<=j); k++) {
 448       sprintf(buf,"table_%s_%s.xvg",atype[k],atype[j]);
 449
 450       printf("%d %d %s\n", j, k, buf);
 451       /* Guillot2001a eqn 2, 6 and 7 */
 452       if (((strcmp(atype[j],"HW") == 0) && (strcmp(atype[k],"HW") == 0)) ||
 453           ((strcmp(atype[j],"OW") == 0) && (strcmp(atype[k],"HW") == 0)) ||
 454           ((strcmp(atype[j],"OW") == 0) && (strcmp(atype[k],"OW") == 0))) {
 455
 456         fp = gmx_ffopen(buf,"w");
 457
 458         imax = 3*pts_nm;
 459         for(i=0; (i<=imax); i++) {
 460           r     = i*(1.0/pts_nm);
 461           /* Avoid very large numbers */
 462           if (r < 0.04) {
 463             vc = fc = vd = fd = vr = fr = 0;
 464           }
 465           else
 466             if (eel == eelPME || eel == eelRF) {
 467               fprintf(stderr, "Not implemented\n");
 468               exit(1);
 469             } else if (eel == eelCUT) {
 470               lo_do_GG_q_q(r,xi,xir,&vc,&fc,&vd,&fd,&vr,&fr);
 471             }
 472           fprintf(fp,"%15.10e   %15.10e %15.10e   %15.10e %15.10e   %15.10e %15.10e\n",
 473                   r,vc,fc,vd,fd,vr,fr);
 474
 475         }
 476         gmx_ffclose(fp);
 477
 478         /* Guillot eqn 4 and 5 */
 479       } else if (((strcmp(atype[j],"HWd") == 0) && (strcmp(atype[k],"HW") == 0)) ||
 480                  ((strcmp(atype[j],"HWd") == 0) && (strcmp(atype[k],"OW") == 0)) ||
 481                  ((strcmp(atype[j],"OWd") == 0) && (strcmp(atype[k],"HW") == 0)) ||
 482                  ((strcmp(atype[j],"OWd") == 0) && (strcmp(atype[k],"OW") == 0))) {
 483
 484         fp = gmx_ffopen(buf,"w");
 485
 486         imax = 3*pts_nm;
 487         for(i=0; (i<=imax); i++) {
 488           r     = i*(1.0/pts_nm);
 489           /* Avoid very large numbers */
 490           if (r < 0.04) {
 491             vc = fc = vd = fd = vr = fr = 0;
 492           }
 493           else
 494             if (eel == eelPME || eel == eelRF) {
 495               fprintf(stderr, "Not implemented\n");
 496               exit(1);
 497             } else if (eel == eelCUT) {
 498               lo_do_GG_q_qd(r,xi,xir,&vc,&fc,&vd,&fd,&vr,&fr);
 499             }
 500           fprintf(fp,"%15.10e   %15.10e %15.10e   %15.10e %15.10e   %15.10e %15.10e\n",
 501                   r,vc,fc,vd,fd,vr,fr);
 502
 503         }
 504         gmx_ffclose(fp);
 505
 506         /* Guillot2001a eqn 3 */
 507       } else if (((strcmp(atype[j],"HWd") == 0) && (strcmp(atype[k],"HWd") == 0)) ||
 508                  ((strcmp(atype[j],"OWd") == 0) && (strcmp(atype[k],"HWd") == 0)) ||
 509                  ((strcmp(atype[j],"OWd") == 0) && (strcmp(atype[k],"OWd") == 0))) {
 510
 511         fp = gmx_ffopen(buf,"w");
 512
 513         imax = 3*pts_nm;
 514         for(i=0; (i<=imax); i++) {
 515           r     = i*(1.0/pts_nm);
 516           /* Avoid very large numbers */
 517           if (r < 0.04) {
 518             vc = fc = vd = fd = vr = fr = 0;
 519           }
 520           else
 521             if (eel == eelPME || eel == eelRF) {
 522               fprintf(stderr, "Not implemented\n");
 523               exit(1);
 524             } else if (eel == eelCUT) {
 525               lo_do_GG_qd_qd(r,xi,xir,&vc,&fc,&vd,&fd,&vr,&fr);
 526             }
 527           fprintf(fp,"%15.10e   %15.10e %15.10e   %15.10e %15.10e   %15.10e %15.10e\n",
 528                   r,vc,fc,vd,fd,vr,fr);
 529
 530         }
 531         gmx_ffclose(fp);
 532
 533       } else
 534         gmx_fatal(FARGS,"Invalid atom type: %s %s", atype[j], atype[k]);
 535
 536
 537     }
 538   }
 539 }
 540
 541 static void do_ljc(FILE *fp,int eel,int pts_nm,real rc,real rtol)
 542 {
 543   int    i,i0,imax;
 544   double r,vc,fc,vd,fd,vr,fr;
 545
 546   imax = 3*pts_nm;
 547   for(i=0; (i<=imax); i++) {
 548     r     = i*(1.0/pts_nm);
 549     /* Avoid very large numbers */
 550     if (r < 0.04) {
 551       vc = fc = vd = fd = vr = fr = 0;
 552     } else {
 553       if (eel == eelPME) {
 554         lo_do_ljc_pme(r,rc,rtol,&vc,&fc,&vd,&fd,&vr,&fr);
 555       } else if (eel == eelCUT) {
 556         lo_do_ljc(r,&vc,&fc,&vd,&fd,&vr,&fr);
 557       }
 558     }
 559     fprintf(fp,"%15.10e   %15.10e %15.10e   %15.10e %15.10e   %15.10e %15.10e\n",
 560             r,vc,fc,vd,fd,vr,fr);
 561   }
 562 }
 563
 564 static void do_guillot_maple(FILE *fp,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 565 {
 566   int    i,i0,imax;
 567   /*  double xi     = 0.15;*/
 568   double r,vc,vc2,vd,vd2,vr,vr2;
 569
 570   imax = 3*pts_nm;
 571   for(i=0; (i<=imax); i++) {
 572     r     = i*(1.0/pts_nm);
 573     /* Avoid very large numbers */
 574     if (r < 0.04) {
 575       vc = vc2 = vd = vd2 = vr = vr2 = 0;
 576     }
 577     else
 578       if (eel == eelPME) {
 579         fprintf(fp, "Not implemented\n");
 580       } else if (eel == eelCUT) {
 581         lo_do_guillot_maple(r,xi,xir,&vc,&vc2,&vd,&vd2,&vr,&vr2);
 582       }
 583     fprintf(fp,"%12.5e  %12.5e  %12.5e   %12.5e  %12.5e  %12.5e  %12.5e\n",
 584             r,vc,vc2,vd,vd2,vr,vr2);
 585   }
 586 }
 587
 588 static void do_GG(FILE *fp,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 589 {
 590   int    i,i0,imax;
 591   double r,vc,vc2,vd,vd2,vr,vr2;
 592
 593   imax = 3*pts_nm;
 594   for(i=0; (i<=imax); i++) {
 595     r     = i*(1.0/pts_nm);
 596     /* Avoid very large numbers */
 597     if (r < 0.04) {
 598       vc = vc2 = vd = vd2 = vr = vr2 = 0;
 599     }
 600     else
 601       if (eel == eelPME) {
 602         fprintf(fp, "Not implemented\n");
 603       } else if (eel == eelCUT) {
 604         lo_do_GG(r,xi,xir,&vc,&vc2,&vd,&vd2,&vr,&vr2);
 605       }
 606     fprintf(fp,"%15.10e   %15.10e %15.10e   %15.10e %15.10e   %15.10e %15.10e\n",
 607             r,vc,vc2,vd,vd2,vr,vr2);
 608   }
 609 }
 610
 611 static void do_GG_q_q(FILE *fp,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 612 {
 613   int    i,i0,imax;
 614   double r,vc,vc2,vd,vd2,vr,vr2;
 615
 616   imax = 3*pts_nm;
 617   for(i=0; (i<=imax); i++) {
 618     r     = i*(1.0/pts_nm);
 619     /* Avoid very large numbers */
 620     if (r < 0.04) {
 621       vc = vc2 = vd = vd2 = vr = vr2 = 0;
 622     }
 623     else
 624       if (eel == eelPME) {
 625         fprintf(fp, "Not implemented\n");
 626       } else if (eel == eelCUT) {
 627         lo_do_GG_q_q(r,xi,xir,&vc,&vc2,&vd,&vd2,&vr,&vr2);
 628       }
 629     fprintf(fp,"%12.5e  %12.5e  %12.5e   %12.5e  %12.5e  %12.5e  %12.5e\n",
 630             r,vc,vc2,vd,vd2,vr,vr2);
 631   }
 632 }
 633
 634 static void do_GG_q_qd(FILE *fp,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 635 {
 636   int    i,i0,imax;
 637   /*  double xi     = 0.15;*/
 638   double r,vc,vc2,vd,vd2,vr,vr2;
 639
 640   imax = 3*pts_nm;
 641   for(i=0; (i<=imax); i++) {
 642     r     = i*(1.0/pts_nm);
 643     /* Avoid very large numbers */
 644     if (r < 0.04) {
 645       vc = vc2 = vd = vd2 = vr = vr2 = 0;
 646     }
 647     else
 648       if (eel == eelPME) {
 649         fprintf(fp, "Not implemented\n");
 650       } else if (eel == eelCUT) {
 651         lo_do_GG_q_qd(r,xi,xir,&vc,&vc2,&vd,&vd2,&vr,&vr2);
 652       }
 653     fprintf(fp,"%12.5e  %12.5e  %12.5e   %12.5e  %12.5e  %12.5e  %12.5e\n",
 654             r,vc,vc2,vd,vd2,vr,vr2);
 655   }
 656 }
 657
 658 static void do_GG_qd_qd(FILE *fp,int eel,int pts_nm,double rc,double rtol,double xi,double xir)
 659 {
 660   int    i,i0,imax;
 661   /*  double xi     = 0.15;*/
 662   double r,vc,vc2,vd,vd2,vr,vr2;
 663
 664   imax = 3*pts_nm;
 665   for(i=0; (i<=imax); i++) {
 666     r     = i*(1.0/pts_nm);
 667     /* Avoid very large numbers */
 668     if (r < 0.04) {
 669       vc = vc2 = vd = vd2 = vr = vr2 = 0;
 670     }
 671     else
 672       if (eel == eelPME) {
 673         fprintf(fp, "Not implemented\n");
 674       } else if (eel == eelCUT) {
 675         lo_do_GG_qd_qd(r,xi,xir,&vc,&vc2,&vd,&vd2,&vr,&vr2);
 676       }
 677     fprintf(fp,"%12.5e  %12.5e  %12.5e   %12.5e  %12.5e  %12.5e  %12.5e\n",
 678             r,vc,vc2,vd,vd2,vr,vr2);
 679   }
 680 }
 681
 682 static void do_maaren(FILE *fp,int eel,int pts_nm,int npow)
 683 {
 684   int    i,i0,imax;
 685   double xi     = 0.05;
 686   double xir     = 0.0615;
 687   double r,vc,vc2,vd,vd2,vr,vr2;
 688
 689   imax = 3*pts_nm;
 690   for(i=0; (i<=imax); i++) {
 691     r     = i*(1.0/pts_nm);
 692     /* Avoid very large numbers */
 693     if (r < 0.04) {
 694       vc = vc2 = vd = vd2 = vr = vr2 = 0;
 695     }
 696     else {
 697       lo_do_guillot_maple(r,xi,xir,&vc,&vc2,&vd,&vd2,&vr,&vr2);
 698       vr  =  pow(r,-1.0*npow);
 699       vr2 = (npow+1.0)*(npow)*vr/sqr(r);
 700     }
 701     fprintf(fp,"%12.5e  %12.5e  %12.5e   %12.5e  %12.5e  %12.5e  %12.5e\n",
 702             r,vc,vc2,vd,vd2,vr,vr2);
 703   }
 704 }
 705
 706 int main(int argc,char *argv[])
 707 {
 708   static char *desc[] = {
 709     "[TT]gen_table[tt] generates tables for [TT]mdrun[tt] for use with the USER defined",
 710     "potentials. Note that the format has been update for higher",
 711     "accuracy in the forces starting with version 4.0. Using older",
 712     "tables with 4.0 will silently crash your simulations, as will",
 713     "using new tables with an older GROMACS version. This is because in the",
 714     "old version the second derevative of the potential was specified",
 715     "whereas in the new version the first derivative of the potential",
 716     "is used instead.[PAR]"
 717   };
 718   static char *opt[]     = { NULL, "cut", "rf", "pme", NULL };
 719   /*  static char *model[]   = { NULL, "guillot", "AB1", "ljc", "maaren", "guillot_maple", "hard_wall", "gg_q_q", "gg_qd_q", "gg_qd_qd", NULL }; */
 720   static char *model[]   = { NULL, "ljc", "gg", "guillot2001a",
 721                              NULL };
 722   static real delta=0,efac=500,rc=0.9,rtol=1e-05,xi=0.15,xir=0.0615;
 723   static real w1=20,w2=20;
 724   static int  nrow1=1,nrow2=1;
 725   static int  nrep       = 12;
 726   static int  ndisp      = 6;
 727   static int  pts_nm     = 500;
 728   t_pargs pa[] = {
 729     { "-el",     FALSE, etENUM, {opt},
 730       "Electrostatics type: cut, rf or pme" },
 731     { "-rc",     FALSE, etREAL, {&rc},
 732       "Cut-off required for rf or pme" },
 733     { "-rtol",   FALSE, etREAL, {&rtol},
 734       "Ewald tolerance required for pme" },
 735     { "-xi",   FALSE, etREAL, {&xi},
 736       "Width of the Gaussian diffuse charge of the G&G model" },
 737     { "-xir",   FALSE, etREAL, {&xir},
 738       "Width of erfc(z)/z repulsion of the G&G model (z=0.5 rOO/xir)" },
 739     { "-m",      FALSE, etENUM, {model},
 740       "Model for the tables" },
 741     { "-resol",  FALSE, etINT,  {&pts_nm},
 742       "Resolution of the table (points per nm)" },
 743     { "-delta",  FALSE, etREAL, {&delta},
 744       "Displacement in the Coulomb functions (nm), used as 1/(r+delta). Only for hard wall potential." },
 745     { "-efac",   FALSE, etREAL, {&efac},
 746       "Number indicating the steepness of the hardwall potential." },
 747     { "-nrep",   FALSE, etINT,  {&nrep},
 748       "Power for the repulsion potential (with model AB1 or maaren)" },
 749     { "-ndisp",   FALSE, etINT,  {&ndisp},
 750       "Power for the dispersion potential (with model AB1 or maaren)" }
 751   };
 752 #define NPA asize(pa)
 753   t_filenm fnm[] = {
 754     { efXVG, "-o", "table", ffWRITE }
 755   };
 756 #define NFILE asize(fnm)
 757   FILE *fp;
 758   char *fn;
 759   int  eel=0,m=0;
 760
 761   CopyRight(stderr,argv[0]);
 762   parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
 763                     NFILE,fnm,NPA,pa,asize(desc),desc,0,NULL);
 764
 765   if (strcmp(opt[0],"cut") == 0)
 766     eel = eelCUT;
 767   else if (strcmp(opt[0],"rf") == 0)
 768     eel = eelRF;
 769   else if (strcmp(opt[0],"pme") == 0)
 770     eel = eelPME;
 771   else
 772     gmx_fatal(FARGS,"Invalid argument %s for option -e",opt[0]);
 773   if (strcmp(model[0],"maaren") == 0)
 774     m = mMaaren;
 775   else if (strcmp(model[0],"AB1") == 0)
 776     m = mAB1;
 777   else if (strcmp(model[0],"ljc") == 0)
 778     m = mLjc;
 779   else if (strcmp(model[0],"guillot2001a") == 0)
 780     m = mGuillot2001a;
 781   else if (strcmp(model[0],"guillot_maple") == 0)
 782     m = mGuillot_Maple;
 783   else if (strcmp(model[0],"hard_wall") == 0)
 784     m = mHard_Wall;
 785   else if (strcmp(model[0],"gg") == 0)
 786     m = mGG;
 787   else if (strcmp(model[0],"gg_qd_q") == 0)
 788     m = mGG_qd_q;
 789   else if (strcmp(model[0],"gg_qd_qd") == 0)
 790     m = mGG_qd_qd;
 791   else if (strcmp(model[0],"gg_q_q") == 0)
 792     m = mGG_q_q;
 793   else
 794     gmx_fatal(FARGS,"Invalid argument %s for option -m",opt[0]);
 795
 796   fn = opt2fn("-o",NFILE,fnm);
 797   if ((m != mGuillot2001a))
 798     fp = gmx_ffopen(fn,"w");
 799   switch (m) {
 800   case mGuillot2001a:
 801     do_guillot2001a(fn,eel,pts_nm,rc,rtol,xi,xir);
 802     break;
 803   case mGuillot_Maple:
 804     fprintf(fp, "#\n# Table Guillot_Maple: rc=%g, rtol=%g, xi=%g, xir=%g\n#\n",rc,rtol,xi,xir);
 805     do_guillot_maple(fp,eel,pts_nm,rc,rtol,xi,xir);
 806     break;
 807   case mGG_q_q:
 808     fprintf(fp, "#\n# Table GG_q_q: rc=%g, rtol=%g, xi=%g, xir=%g\n#\n",rc,rtol,xi,xir);
 809     do_GG_q_q(fp,eel,pts_nm,rc,rtol,xi,xir);
 810     break;
 811   case mGG:
 812     fprintf(fp, "#\n# Table GG: rc=%g, rtol=%g, xi=%g, xir=%g\n#\n",rc,rtol,xi,xir);
 813     do_GG(fp,eel,pts_nm,rc,rtol,xi,xir);
 814     break;
 815   case mGG_qd_q:
 816     fprintf(stdout, "case mGG_qd_q");
 817     fprintf(fp, "#\n# Table GG_qd_q: rc=%g, rtol=%g, xi=%g, xir=%g\n#\n",rc,rtol,xi,xir);
 818     do_GG_q_qd(fp,eel,pts_nm,rc,rtol,xi,xir);
 819     break;
 820   case mGG_qd_qd:
 821     fprintf(stdout, "case mGG_qd_qd");
 822     fprintf(fp, "#\n# Table GG_qd_qd: rc=%g, rtol=%g, xi=%g, xir=%g\n#\n",rc,rtol,xi,xir);
 823     do_GG_qd_qd(fp,eel,pts_nm,rc,rtol,xi,xir);
 824     break;
 825   case mMaaren:
 826     do_maaren(fp,eel,pts_nm,nrep);
 827     break;
 828   case mAB1:
 829     fprintf(fp, "#\n# Table AB1: ndisp=%d nrep=%d\n#\n",ndisp,nrep);
 830     do_AB1(fp,eel,pts_nm,ndisp,nrep);
 831     break;
 832   case mLjc:
 833     fprintf(fp, "#\n# Table LJC(12-6-1): rc=%g, rtol=%g\n#\n",rc,rtol);
 834     do_ljc(fp,eel,pts_nm,rc,rtol);
 835     break;
 836   case mHard_Wall:
 837     do_hard(fp,pts_nm,efac,delta);
 838     break;
 839   default:
 840     gmx_fatal(FARGS,"Model %s not supported yet",model[0]);
 841   }
 842   if ((m != mGuillot2001a))
 843     gmx_ffclose(fp);
 844
 845   gmx_thanx(stdout);
 846
 847   return 0;
 848 }