src/kernel/mk_ghat.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include <stdio.h>
  41 #include "copyrite.h"
  42 #include "macros.h"
  43 #include "smalloc.h"
  44 #include "typedefs.h"
  45 #include "futil.h"
  46 #include "xvgr.h"
  47 #include "vec.h"
  48 #include "maths.h"
  49 #include "coulomb.h"
  50 #include "physics.h"
  51 #include "statutil.h"
  52 #include "tpxio.h"
  53 #include "copyrite.h"
  54
  55 const real tol = 1e-8;
  56
  57 real sinx_x(real x)
  58 {
  59   if (x == 0)
  60     return 1.0;
  61   else
  62     return sin(x)/x;
  63 }
  64
  65 real  uhat(int porder,real k1,real k2,real k3,real h1,real h2,real h3)
  66 {
  67   real fac1,fac2,fac3;
  68   real f123,fff;
  69   int  i;
  70
  71   fac1 = sinx_x(k1*h1*0.5);
  72   fac2 = sinx_x(k2*h2*0.5);
  73   fac3 = sinx_x(k3*h3*0.5);
  74
  75   fff  = 1;
  76   f123 = fac1*fac2*fac3;
  77   for(i=1; (i<=porder+1); i++)
  78     fff *= f123;
  79
  80   if (fabs(fff) < tol)
  81     return 0.0;
  82   else
  83     return fff;
  84 }
  85
  86 real shat(real acut,real kmag,real r1)
  87 {
  88   return gk(kmag,acut,r1);
  89 }
  90
  91 real dhat(real alpha,real k,real h)
  92 {
  93   real dh;
  94
  95   dh = alpha*(sin(k*h)/h) + (1.0-alpha)*(sin(2.0*k*h)/(2.0*h) );
  96
  97   if (fabs(dh) < tol)
  98     return 0.0;
  99   else
 100     return dh;
 101 }
 102
 103 real cufrth(int porder,real k1,real k2,real k3,
 104             real h1,real h2,real h3,int nalias)
 105 {
 106   real kn1,kn2,kn3,tmp;
 107   int  n1,n2,n3;
 108   real ufrth;
 109
 110   real twopi=2*M_PI;
 111
 112   if (nalias == 0) {
 113     tmp   = uhat(porder,k1,k2,k3,h1,h2,h3);
 114     ufrth = tmp*tmp*tmp*tmp;
 115   }
 116   else {
 117     ufrth = 0.0;
 118     for( n1 = -nalias; (n1<= nalias); n1++) {
 119       for( n2 = -nalias; (n2<= nalias); n2++) {
 120         for( n3 = -nalias; (n3<= nalias); n3++) {
 121           kn1 = k1 + n1*twopi/h1;
 122           kn2 = k2 + n2*twopi/h2;
 123           kn3 = k3 + n3*twopi/h3;
 124           tmp = uhat(porder,kn1,kn2,kn3,h1,h2,h3);
 125           ufrth += tmp*tmp*tmp*tmp;
 126         }
 127       }
 128     }
 129   }
 130   if (fabs(ufrth) < tol)
 131     return 0.0;
 132   else
 133     return ufrth;
 134 }
 135
 136 real crsqal(real acut,real r1,real k1,real k2,real k3,
 137             real h1,real h2,real h3,int nalias)
 138 {
 139   real kn1,kn2,kn3;
 140   real h_1,h_2,h_3;
 141   real ksq,kmag,tmp,Rsqal;
 142   int  n1,n2,n3;
 143
 144   real twopi=2*M_PI;
 145
 146   h_1=twopi/h1;
 147   h_2=twopi/h2;
 148   h_3=twopi/h3;
 149   if (nalias==0) {
 150     ksq   = k1*k1 + k2*k2 + k3*k3;
 151     kmag  = sqrt(ksq);
 152     tmp   = shat(acut,kmag,r1);
 153     Rsqal = tmp*tmp/(ksq);
 154   }
 155   else {
 156     Rsqal = 0.0;
 157     for(n1 = -nalias; (n1<= nalias); n1++) {
 158       for( n2 = -nalias; (n2<= nalias); n2++) {
 159         for( n3 = -nalias; (n3<= nalias); n3++) {
 160           kn1   = k1 + n1*h_1;
 161           kn2   = k2 + n2*h_2;
 162           kn3   = k3 + n3*h_3;
 163           ksq   = kn1*kn1 + kn2*kn2 + kn3*kn3;
 164           kmag  = sqrt(ksq);
 165           tmp   = shat(acut,kmag,r1);
 166           Rsqal = Rsqal + tmp*tmp/ksq;
 167         }
 168       }
 169     }
 170   }
 171   return Rsqal;
 172 }
 173
 174
 175 real usqsq(int porder,real k1,real k2,real k3,real h1,real h2,real h3)
 176 {
 177   const real tt=2.0/3.0;
 178   const real tx=2.0/15.0;
 179   real t1,t2,t3,t12,t22,t32,tmp;
 180
 181   t1 = sin(k1*h1*0.5);
 182   t2 = sin(k2*h2*0.5);
 183   t3 = sin(k3*h3*0.5);
 184   t12 = t1*t1;
 185   t22 = t2*t2;
 186   t32 = t3*t3;
 187
 188   if (porder == 1)
 189     tmp = (1.0-tt*t12)*(1-tt*t22)*(1-tt*t32);
 190   else if (porder == 2)
 191     tmp = ( (1.0 - t12 + tx*t12*t12)*
 192             (1.0 - t22 + tx*t22*t22)*
 193             (1.0 - t32 + tx*t32*t32) );
 194   else
 195     gmx_fatal(FARGS,"porder = %d in usqsq",porder);
 196
 197   return tmp*tmp;
 198 }
 199
 200 real  ursum(int term,int porder,real acut,real r1,
 201             real k1,real k2,real k3,real h1,real h2,real h3,int nalias)
 202 {
 203   real kt,ksq,kmag;
 204   /*   real kcutsq; */
 205   real kn1,kn2,kn3,urs,tmp;
 206   real h_1,h_2,h_3;
 207   int  n1,n2,n3;
 208
 209   real twopi=2*M_PI;
 210   h_1=twopi/h1;
 211   h_2=twopi/h2;
 212   h_3=twopi/h3;
 213   /*
 214     c
 215     c     for large enough values of k, the terms become negligable
 216     c     if shat(k) = exp(-k^2/4*acut) < eps
 217     c     kcutsq = 4*alpha* (-ln(eps))
 218     c     eps = 10^-6, -ln(eps) = 14
 219     c     eps = 10^-10, -ln(eps) = 23
 220     c     eps = 10^-20, -ln(eps) = 46
 221     c
 222     c     kcutsq = 4.0*acut*115;
 223   */
 224
 225   if (nalias==0) {
 226     if (term==XX) kt = k1;
 227     if (term==YY) kt = k2;
 228     if (term==ZZ) kt = k3;
 229     ksq = k1*k1 + k2*k2 + k3*k3;
 230     kmag = sqrt(ksq);
 231     tmp = uhat(porder,k1,k2,k3,h1,h2,h3);
 232     urs = tmp*tmp*kt*shat(acut,kmag,r1)/(ksq);
 233   }
 234   else {
 235     urs = 0.0;
 236     for(n1 = -nalias; (n1<= nalias); n1++) {
 237       for( n2 = -nalias; (n2<= nalias); n2++) {
 238         for( n3 = -nalias; (n3<= nalias); n3++) {
 239           kn1 = k1 + n1*h_1;
 240           kn2 = k2 + n2*h_2;
 241           kn3 = k3 + n3*h_3;
 242           ksq = kn1*kn1 + kn2*kn2 + kn3*kn3;
 243
 244           if (term==XX) kt = kn1;
 245           if (term==YY) kt = kn2;
 246           if (term==ZZ) kt = kn3;
 247           if (kt != 0.0) {
 248             kmag = sqrt(ksq);
 249             tmp = uhat(porder,kn1,kn2,kn3,h1,h2,h3);
 250             if (tmp != 0.0)
 251               urs = urs + tmp*tmp*kt*shat(acut,kmag,r1)/ksq;
 252           }
 253         }
 254       }
 255     }
 256   }
 257   return urs;
 258 }
 259
 260 real sym(int indx,int maxind)
 261 {
 262   if ( (indx == 0 ) || (indx == maxind/2) )
 263     return 1.0;
 264   else
 265     return 2.0;
 266 }
 267
 268 void calc(gmx_bool bSym,gmx_bool bVerbose,
 269           const int n1max,const int n2max,const int n3max,
 270           const real h1,const real h2,const real h3,
 271           int nalias,int porder,real acut,real r1,const real alpha,
 272           const gmx_bool bSearch,
 273           real ***ghat,real *ppval,real *zzval,real *eeref,real *qqopt)
 274 {
 275   real box1,box2,box3;
 276   real k1,k2,k3,ksq,kmag;
 277   real gnumer,dsq,gdenom,gsq;
 278   real ufrth,rsqal,rsq;
 279   real symfac;
 280   int  l1,l2,l3;
 281   real twopi=2*M_PI;
 282   real d1,d2,d3,u1,u2,u3,ss,gg;
 283   real pval,zval,eref,qopt;
 284   int  N1MAX,N2MAX,N3MAX;
 285
 286   if (bSym) {
 287     N1MAX = n1max/2+1;
 288     N2MAX = n2max/2+1;
 289     N3MAX = n3max/2+1;
 290   }
 291   else {
 292     N1MAX = n1max;
 293     N2MAX = n2max;
 294     N3MAX = n3max;
 295   }
 296
 297   box1 = n1max*h1;
 298   box2 = n2max*h2;
 299   box3 = n3max*h3;
 300
 301   pval = 0.0;
 302   zval = 0.0;
 303   eref = 0.0;
 304   qopt = 0.0;
 305
 306   for(l1=0; (l1<N1MAX); l1++) {
 307     if (bVerbose)
 308       fprintf(stderr,"\rl1=%5d  qopt=%12.6e",l1,qopt);
 309
 310     k1   = twopi*l1/box1;
 311     d1   = dhat(alpha,k1,h1);
 312
 313     for(l2=0; (l2<N2MAX); l2++) {
 314       k2   = twopi*l2/box2;
 315       d2   = dhat(alpha,k2,h2);
 316
 317       for(l3=0; (l3<N3MAX); l3++) {
 318         if (((l1+l2+l3) == 0) /*|| (l1 == n1max/2) || (l2 == n2max/2) ||
 319             (l3 == n3max/2)*/)
 320           ghat[0][0][0] = 0.0;
 321         else {
 322           k3   = twopi*l3/box3;
 323           ksq  = k1*k1 + k2*k2 + k3*k3;
 324           kmag = sqrt(ksq);
 325
 326           d3   = dhat(alpha,k3,h3);
 327           u1   = ursum(XX,porder,acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 328           u2   = ursum(YY,porder,acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 329           u3   = ursum(ZZ,porder,acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 330
 331           gnumer = d1*u1+d2*u2+d3*u3;
 332           dsq    = d1*d1+d2*d2+d3*d3;
 333           gdenom = dsq*usqsq(porder,k1,k2,k3,h1,h2,h3);
 334           if (bSym)
 335             symfac = sym(l1,n1max)*sym(l2,n2max)*sym(l3,n3max);
 336           else
 337             symfac = 1.0;
 338
 339           rsqal  = crsqal(acut,r1,k1,k2,k3,h1,h2,h3,nalias);
 340
 341           if (gdenom != 0)
 342             qopt  += (symfac*(rsqal - sqr(gnumer)/gdenom));
 343           if (debug)
 344             fprintf(debug,"rsqal: %10.3e, gnumer: %10.3e, gdenom: %10.3e, ratio: %10.3e\n",
 345                     rsqal,gnumer,gdenom,gnumer/gdenom);
 346
 347 #ifdef DEBUG
 348           if ((l1 == n1max/2) || (l2 == n2max/2) || (l3 == n3max/2))
 349             printf("L(%2d,%2d,%2d)  D(%10.3e,%10.3e,%10.3e) U(%10.3e,%10.3e,%10.3e) gnumer=%10.3em dsq=%10.3e, gdenom=%10.3e, ghat=%10.3e\n",
 350                    l1,l2,l3,d1,d2,d3,u1,u2,u3,gnumer,dsq,gdenom,
 351                    (gdenom == 0.0) ? 0 : gnumer/gdenom);
 352 #endif
 353           if (!bSearch) {
 354             if (gdenom != 0)
 355               gg = gnumer/gdenom;
 356             else
 357               gg = 0.0;
 358             ghat[l1][l2][l3] = gg/EPSILON0;
 359             gsq = gg*gg;
 360
 361             ufrth = cufrth(porder,k1,k2,k3,h1,h2,h3,nalias);
 362             pval  = pval + symfac*
 363               (dsq*gsq*(usqsq(porder,k1,k2,k3,h1,h2,h3)-ufrth));
 364             if (gdenom != 0)
 365               zval  = zval + symfac*
 366                 (dsq*gsq*ufrth - 2.0*gnumer*gnumer/gdenom + rsqal);
 367             ss    = shat(acut,kmag,r1);
 368             rsq   = ss*ss/ksq;
 369             eref  = eref + symfac* (rsqal - rsq);
 370           }
 371         }
 372       }
 373     }
 374   }
 375   if (bVerbose)
 376     fprintf(stderr,"\n");
 377   *ppval = pval/(box1*box2*box3);
 378   *zzval = zval/(box1*box2*box3);
 379   *eeref = eref/(box1*box2*box3);
 380   *qqopt = qopt/(EPSILON0*box1*box2*box3);
 381 }
 382
 383 void read_params(char *fn,t_inputrec *ir,rvec boxs)
 384 {
 385   real   t,lambda;
 386   int    step,natoms,m;
 387   matrix box;
 388
 389   /* Read topology and coordinates */
 390   read_tpx(fn,&step,&t,&lambda,ir,box,&natoms,NULL,NULL,NULL,NULL);
 391   for(m=0; (m<DIM); m++)
 392     boxs[m] = box[m][m];
 393 }
 394
 395 int main(int argc,char *argv[])
 396 {
 397   /* FILE   *fp; */
 398   const  real gold=0.38197;
 399   const  int  mxiter=12;
 400   int    n1max,n2max,n3max;
 401   real   h1,h2,h3;
 402   real   box[3];
 403   int    nalias,porder;
 404   real   acut,alpha,r1;
 405   rvec   beta;
 406   gmx_bool   bSearch,bConv;
 407   /* gmx_bool   bNL=FALSE; */
 408   real   ***ghat;
 409   real   pval,zval,eref,qopt,norm;
 410   real   alpha0,alpha1,alpha2,alpha3,alptol;
 411   real   q0,q1,q2,q3;
 412   int    niter;
 413   /* int    ii,jj,kk,nn; */
 414   t_inputrec ir;
 415   t_filenm fnm[] = {
 416     { efTPX, NULL, NULL,   ffREAD },
 417     { efHAT, "-o", "ghat", ffWRITE }
 418   };
 419 #define NFILE asize(fnm)
 420   static gmx_bool bVerbose=FALSE,bCubic=TRUE,bSym=TRUE;
 421   static t_pargs pa[] = {
 422     { "-v",     FALSE, etBOOL, &bVerbose, "Verbose on"},
 423     { "-cubic", FALSE, etBOOL, &bCubic,   "Force beta to be the same in all directions" },
 424     { "-sym",   FALSE, etBOOL, &bSym,     "HIDDENUse symmetry for the generation of ghat function (turn off for debugging only!)" }
 425   };
 426
 427   CopyRight(stdout,argv[0]);
 428   parse_common_args(&argc,argv,0,TRUE,NFILE,fnm,asize(pa),pa,0,NULL,0,NULL);
 429
 430   read_params(ftp2fn(efTPX,NFILE,fnm),&ir,box);
 431   acut   = ir.rcoulomb;
 432   r1     = ir.rcoulomb_switch;
 433   n1max  = ir.nkx;
 434   n2max  = ir.nky;
 435   n3max  = ir.nkz;
 436   h1     = box[XX]/n1max;
 437   h2     = box[YY]/n2max;
 438   h3     = box[ZZ]/n3max;
 439
 440   /* These are not parameters. They are fixed.
 441    * Luty et al. determined their optimal values to be 2.
 442    */
 443   nalias = 2;
 444   porder = 2;
 445   /*bSym   = TRUE;*/
 446
 447   set_LRconsts(stdout,r1,acut,box,NULL);
 448
 449   printf("Grid cell size is %8.4f x %8.4f x %8.4f nm\n",h1,h2,h3);
 450
 451   ghat=mk_rgrid(n1max,n2max,n3max);
 452
 453   bSearch = TRUE;
 454
 455   niter = 0;
 456   if (!bSearch) {
 457     /* c          alpha = 1.33*/
 458     alpha = 1.0;
 459     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 460          nalias,porder,acut,r1,alpha,bSearch,
 461          ghat,&pval,&zval,&eref,&qopt);
 462   }
 463   else /* Elsje */ {
 464     alpha0 = 1.0;
 465     alpha1 = 4.0/3.0;
 466     alpha3 = 2.0;
 467     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 468          nalias,porder,acut,r1,alpha0,bSearch,
 469          ghat,&pval,&zval,&eref,&q0);
 470     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 471          nalias,porder,acut,r1,alpha1,bSearch,
 472          ghat,&pval,&zval,&eref,&q1);
 473     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 474          nalias,porder,acut,r1,alpha3,bSearch,
 475          ghat,&pval,&zval,&eref,&q3);
 476
 477     /* if ( (q1 > q0) || (q1 > q3) )
 478       gmx_fatal(FARGS,"oops q1=%f,q0=%f,q3=%f",q1,q0,q3); */
 479     alpha2 = alpha1 + gold*(alpha3-alpha1);
 480     calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 481          nalias,porder,acut,r1,alpha2,bSearch,
 482          ghat,&pval,&zval,&eref,&q2);
 483     bConv  = FALSE;
 484     alptol = 0.05;
 485
 486     fprintf(stderr,"q0 = %10g, q1= %10g, q2 = %10g, q3 = %10g\n",
 487             q0,q1,q2,q3);
 488
 489     while ((!bConv) && (niter < mxiter)) {
 490       fprintf(stderr,"%2d  %10.4f  %10.4f  %10.4f  %10.4f\n",
 491               niter,alpha0,alpha1,alpha2,alpha3);
 492       niter = niter + 1;
 493       if (q2 < q1) {
 494         alpha0 = alpha1;
 495         q0     = q1;
 496         alpha1 = alpha2;
 497         q1     = q2;
 498         alpha2 = alpha1 + gold*(alpha3-alpha1);
 499         calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 500              nalias,porder,acut,r1,alpha2,bSearch,
 501              ghat,&pval,&zval,&eref,&q2);
 502       }
 503       else {
 504         alpha3 = alpha2;
 505         q3     = q2;
 506         alpha2 = alpha1;
 507         q2     = q1;
 508         alpha1 = alpha2 - gold*(alpha2-alpha0);
 509         calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 510              nalias,porder,acut,r1,alpha1,bSearch,
 511              ghat,&pval,&zval,&eref,&q1);
 512       }
 513       if ((alpha3-alpha0) < alptol)
 514         bConv = TRUE;
 515     }
 516     bSearch = FALSE;
 517     if (q1 < q2) {
 518       alpha = alpha1;
 519       calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 520            nalias,porder,acut,r1,alpha,bSearch,
 521            ghat,&pval,&zval,&eref,&qopt);
 522     }
 523     else {
 524       alpha = alpha2;
 525       calc(bSym,bVerbose,n1max,n2max,n3max,h1,h2,h3,
 526            nalias,porder,acut,r1,alpha,bSearch,
 527            ghat,&pval,&zval,&eref,&qopt);
 528     }
 529
 530     fprintf(stderr,"%10g  %10g  %10g  %10g  %10g  %10g\n",
 531             acut,r1,pval,zval,eref,qopt);
 532     norm = sqr(1.0/(4.0*M_PI*h1*h1))*(4.0/3.0*M_PI*h1*h1*h1);
 533     if (norm != 0) {
 534       pval = sqrt(fabs(pval)/norm)*100.0;
 535       zval = sqrt(fabs(zval)/norm)*100.0;
 536       eref = sqrt(fabs(eref)/norm)*100.0;
 537       qopt = sqrt(fabs(qopt)/norm)*100.0;
 538     }
 539     beta[XX] = beta[YY] = beta[ZZ] = alpha;
 540     wr_ghat(ftp2fn(efHAT,NFILE,fnm),n1max,n2max,n3max,h1,h2,h3,
 541             ghat,nalias,porder,niter,bSym,beta,
 542             r1,acut,pval,zval,eref,qopt);
 543
 544     /*    fp=ftp2FILE(efHAT,NFILE,fnm,"w");
 545           fprintf(fp,"%8d  %8d  %8d  %15.10e  %15.10e %15.10e\n",
 546           n1max,n2max,n3max,h1,h2,h3);
 547           fprintf(fp,"%8d  %8d  %8d  %8d  %15.10e  %15.10e  %15.10e\n",
 548           nalias,porder,niter,bSym,alpha,alpha,alpha);
 549           fprintf(fp,"%10g  %10g  %10g  %10g  %10g  %10g\n",
 550           acut,r1,pval,zval,eref,qopt);
 551           {
 552           int  N1MAX,N2MAX,N3MAX;
 553
 554           if (bSym) {
 555           N1MAX = n1max/2+1;
 556           N2MAX = n2max/2+1;
 557           N3MAX = n3max/2+1;
 558           }
 559           else {
 560           N1MAX = n1max;
 561           N2MAX = n2max;
 562           N3MAX = n3max;
 563           }
 564           for(ii=0; (ii<N1MAX); ii++) {
 565           for(jj=0; (jj<N2MAX); jj++) {
 566           for(kk=0,nn=1; (kk<N3MAX); kk++,nn++) {
 567           bNL=FALSE;
 568           fprintf(fp,"  %12.5e",ghat[ii][jj][kk]);
 569           if ((nn % 6) == 0) {
 570           fprintf(fp,"\n");
 571           bNL=TRUE;
 572           }
 573           }
 574           if (!bNL)
 575           fprintf(fp,"\n");
 576           }
 577           }
 578           ffclose(fp);
 579           }
 580     */
 581     pr_scalar_gk("ghat.xvg",n1max,n2max,n3max,box,ghat);
 582   }
 583   return 0;
 584 }
 585