src/gmxlib/ewald_util.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_ewald_util_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include "assert.h"
  34 #include "typedefs.h"
  35 #include "vec.h"
  36 #include "ewald_util.h"
  37 #include "smalloc.h"
  38 #include "physics.h"
  39 #include "txtdump.h"
  40 #include "futil.h"
  41 #include "names.h"
  42 #include "fftgrid.h"
  43 #include "writeps.h"
  44 #include "macros.h"
  45 #include "xvgr.h"
  46
  47 real calc_ewaldcoeff(real rc,real dtol)
  48 {
  49   real x=5,low,high;
  50   int n,i=0;
  51
  52
  53   do {
  54     i++;
  55     x*=2;
  56   } while((erfc(x*rc)/rc)>dtol);
  57
  58   n=i+60; /* search tolerance is 2^-60 */
  59   low=0;
  60   high=x;
  61   for(i=0;i<n;i++) {
  62     x=(low+high)/2;
  63     if((erfc(x*rc)/rc)>dtol)
  64       low=x;
  65     else
  66       high=x;
  67   }
  68   return x;
  69 }
  70
  71
  72
  73 real ewald_LRcorrection(FILE *fp,t_nsborder *nsb,t_commrec *cr,t_forcerec *fr,
  74                         real charge[],t_block *excl,rvec x[],
  75                         rvec box_size,matrix lr_vir)
  76 {
  77   static  bool bFirst=TRUE;
  78   static  real Vself;
  79
  80   int     i,i1,i2,j,k,m,iv,jv;
  81   atom_id *AA;
  82   double  qq; /* Necessary for precision */
  83   real    vc,qi,dr,ddd,dr2,rinv,fscal,Vexcl,rinv2,ewc=fr->ewaldcoeff;
  84   rvec    df,dx;
  85   rvec    *flr=fr->flr;
  86   /*#define TABLES*/
  87 #ifdef TABLES
  88   real    tabscale=fr->tabscale;
  89   real    eps,eps2,VV,FF,F,Y,Geps,Heps2,Fp,fijC,r1t;
  90   real    *VFtab=fr->VFtab;
  91   int     n0,n1,nnn;
  92 #else
  93   double  isp=0.564189583547756;
  94   real    dr_3;
  95 #endif
  96   int     start = START(nsb);
  97   int     end   = start+HOMENR(nsb);
  98
  99   if (bFirst) {
 100     qq =0;
 101     for(i=start; (i<end); i++)
 102       qq  += charge[i]*charge[i];
 103
 104     /* Charge and dipole correction remain to be implemented... */
 105
 106     Vself=ewc*ONE_4PI_EPS0*qq/sqrt(M_PI);
 107     if (debug) {
 108       fprintf(fp,"Long Range corrections for Ewald interactions:\n");
 109       fprintf(fp,"start=%d,natoms=%d\n",start,end-start);
 110       fprintf(fp,"qq = %g, Vself=%g\n",qq,Vself);
 111     }
 112   }
 113
 114   AA    = excl->a;
 115   Vexcl = 0;
 116
 117   for(i=start; (i<end); i++) {
 118       /* Initiate local variables (for this i-particle) to 0 */
 119       qi  = charge[i]*ONE_4PI_EPS0;
 120       i1  = excl->index[i];
 121       i2  = excl->index[i+1];
 122       /* Loop over excluded neighbours */
 123       for(j=i1; (j<i2); j++) {
 124         k = AA[j];
 125           /*
 126            * First we must test whether k <> i, and then, because the
 127            * exclusions are all listed twice i->k and k->i we must select
 128            * just one of the two.
 129        * As a minor optimization we only compute forces when the charges
 130        * are non-zero.
 131        */
 132       if (k > i) {
 133         qq = qi*charge[k];
 134         if (qq != 0.0) {
 135           /* Compute distance vector, no PBC check! */
 136           dr2 = 0;
 137           for(m=0; (m<DIM); m++) {
 138             ddd = x[i][m] - x[k][m];
 139               if(ddd>box_size[m]/2) {  /* ugly hack,   */
 140                 ddd-=box_size[m];      /* to fix pbc.. */
 141                 /* Can this be done better? */
 142               } else if (ddd<-box_size[m]/2)
 143                   ddd+=box_size[m];
 144
 145             dx[m] = ddd;
 146             dr2  += ddd*ddd;
 147           }
 148
 149           /* It might be possible to optimize this slightly
 150            * when using spc or similar water models:
 151            * Potential could be stored once, at the beginning,
 152            * and forces could use that bonds are constant
 153            */
 154           /* The Ewald interactions are tabulated. If you
 155            * remove the tabulation you must replace this
 156            * part too
 157            */
 158
 159           rinv              = invsqrt(dr2);
 160           rinv2             = rinv*rinv;
 161           dr                = 1.0/rinv;
 162 #ifdef TABLES
 163           r1t               = tabscale*dr;
 164           n0                = r1t;
 165           assert(n0 >= 3);
 166           n1                = 12*n0;
 167           eps               = r1t-n0;
 168           eps2              = eps*eps;
 169           nnn               = n1;
 170           Y                 = VFtab[nnn];
 171           F                 = VFtab[nnn+1];
 172           Geps              = eps*VFtab[nnn+2];
 173           Heps2             = eps2*VFtab[nnn+3];
 174           Fp                = F+Geps+Heps2;
 175           VV                = Y+eps*Fp;
 176           FF                = Fp+Geps+2.0*Heps2;
 177           vc                = qq*(rinv-VV);
 178           fijC              = qq*FF;
 179           Vexcl            += vc;
 180
 181           fscal             = vc*rinv2+fijC*tabscale*rinv;
 182           /* End of tabulated interaction part */
 183 #else
 184
 185           /* This is the code you would want instead if not using
 186            * tables:
 187            */
 188           dr_3    = rinv*rinv2;
 189           vc      = qq*erf(ewc*dr)*rinv;
 190           Vexcl  += vc;
 191           fscal   = rinv2*(vc-2.0*qq*ewc*isp*exp(-ewc*ewc*dr2));
 192 #endif
 193
 194           /* The force vector is obtained by multiplication with the
 195            * distance vector
 196            */
 197           svmul(fscal,dx,df);
 198           if (debug)
 199             fprintf(debug,"dr=%8.4f, fscal=%8.0f, df=%10.0f,%10.0f,%10.0f\n",
 200                     dr,fscal,df[XX],df[YY],df[ZZ]);
 201           rvec_inc(flr[k],df);
 202           rvec_dec(flr[i],df);
 203           for(iv=0; (iv<DIM); iv++)
 204             for(jv=0; (jv<DIM); jv++)
 205               lr_vir[iv][jv]+=0.5*dx[iv]*df[jv];
 206         }
 207       }
 208     }
 209   }
 210   if (bFirst && debug)
 211     fprintf(fp,"Long Range correction: Vexcl=%g\n",Vexcl);
 212
 213   bFirst = FALSE;
 214
 215   return (Vself+Vexcl);
 216 }
 217