src/contrib/gen_table.c

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