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