src/mdlib/psgather.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 3.1
  11  * Copyright (c) 1991-2001, University of Groningen, The Netherlands
  12  * This program is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU General Public License
  14  * as published by the Free Software Foundation; either version 2
  15  * of the License, or (at your option) any later version.
  16  *
  17  * If you want to redistribute modifications, please consider that
  18  * scientific software is very special. Version control is crucial -
  19  * bugs must be traceable. We will be happy to consider code for
  20  * inclusion in the official distribution, but derived work must not
  21  * be called official GROMACS. Details are found in the README & COPYING
  22  * files - if they are missing, get the official version at www.gromacs.org.
  23  *
  24  * To help us fund GROMACS development, we humbly ask that you cite
  25  * the papers on the package - you can find them in the top README file.
  26  *
  27  * For more info, check our website at http://www.gromacs.org
  28  *
  29  * And Hey:
  30  * GROup of MAchos and Cynical Suckers
  31  */
  32 static char *SRCID_psgather_c = "$Id$";
  33 #include <math.h>
  34 #include "poisson.h"
  35 #include "nrnb.h"
  36 #include "shift_util.h"
  37
  38 real ps_gather_inner(int JCXYZ[],real WXYZ[],int ixw[],int iyw[],int izw[],
  39                      real c1x,real c1y,real c1z,real c2x,real c2y,real c2z,
  40                      real qi,rvec f,real ***ptr)
  41 {
  42   real pi,fX,fY,fZ,weight;
  43   int  jxyz,m,jcx,jcy,jcz;
  44   int  jcx0,jcy0,jcz0;
  45
  46   pi = 0.0;
  47   fX = 0.0;
  48   fY = 0.0;
  49   fZ = 0.0;
  50
  51   /* Now loop over 27 surrounding vectors */
  52   for(jxyz=m=0; (jxyz < 27); jxyz++,m+=3) {
  53     jcx    = JCXYZ[m];
  54     jcy    = JCXYZ[m+1];
  55     jcz    = JCXYZ[m+2];
  56     weight = WXYZ[jxyz];
  57
  58     jcx0   = ixw[jcx];
  59     jcy0   = iyw[jcy];
  60     jcz0   = izw[jcz];
  61
  62     /* Electrostatic Potential ! */
  63     pi += weight * ptr[jcx0][jcy0][jcz0];
  64
  65     /* Forces */
  66     fX += weight * ((c1x*(ptr[ixw[jcx-1]] [jcy0]       [jcz0] -
  67                           ptr[ixw[jcx+1]] [jcy0]       [jcz0] )) +
  68                     (c2x*(ptr[ixw[jcx-2]] [jcy0]       [jcz0] -
  69                           ptr[ixw[jcx+2]] [jcy0]       [jcz0] )));
  70     fY += weight * ((c1y*(ptr[jcx0]       [iyw[jcy-1]] [jcz0] -
  71                           ptr[jcx0]       [iyw[jcy+1]] [jcz0] ))  +
  72                     (c2y*(ptr[jcx0]       [iyw[jcy-2]] [jcz0] -
  73                           ptr[jcx0]       [iyw[jcy+2]] [jcz0] )));
  74     fZ += weight * ((c1z*(ptr[jcx0]       [jcy0]       [izw[jcz-1]] -
  75                           ptr[jcx0]       [jcy0]       [izw[jcz+1]] ))  +
  76                     (c2z*(ptr[jcx0]       [jcy0]       [izw[jcz-2]] -
  77                           ptr[jcx0]       [jcy0]       [izw[jcz+2]] )));
  78   }
  79   f[XX] += qi*fX;
  80   f[YY] += qi*fY;
  81   f[ZZ] += qi*fZ;
  82
  83   return pi;
  84 }
  85
  86 real ps_gather_f(FILE *log,bool bVerbose,
  87                  int natoms,rvec x[],rvec f[],real charge[],rvec box,
  88                  real pot[],t_PSgrid *grid,rvec beta,t_nrnb *nrnb)
  89 {
  90   static bool bFirst=TRUE;
  91   static int  *nnx,*nny,*nnz;
  92   static int  JCXYZ[81];
  93   int    i,m;
  94   real   energy;
  95   real   qi,pi;
  96   ivec   ixyz;
  97   rvec   invh,h,c1,c2;
  98   real   WXYZ[27];
  99   real   c1x,c1y,c1z,c2x,c2y,c2z;
 100   int    ixw[7],iyw[7],izw[7];
 101   int    ll;
 102   int    nx,ny,nz;
 103   real   ***ptr;
 104
 105   unpack_PSgrid(grid,&nx,&ny,&nz,&ptr);
 106
 107   calc_invh_h(box,nx,ny,nz,invh,h);
 108
 109   for(m=0; (m<DIM); m++) {
 110     c1[m] = (beta[m]/2.0)*invh[m];
 111     c2[m] = ((1.0-beta[m])/4.0)*invh[m];
 112   }
 113   c1x = c1[XX];
 114   c1y = c1[YY];
 115   c1z = c1[ZZ];
 116   c2x = c2[XX];
 117   c2y = c2[YY];
 118   c2z = c2[ZZ];
 119
 120   if (bFirst) {
 121     fprintf(log,"Gathering Forces using Triangle Shaped on %dx%dx%d grid\n",
 122             nx,ny,nz);
 123     fprintf(log,"beta = %10g,%10g,%10g\n",beta[XX],beta[YY],beta[ZZ]);
 124     fprintf(log,"c1   = %10g,%10g,%10g\n",c1[XX],c1[YY],c1[ZZ]);
 125     fprintf(log,"c2   = %10g,%10g,%10g\n",c2[XX],c2[YY],c2[ZZ]);
 126     fprintf(log,"invh = %10g,%10g,%10g\n",invh[XX],invh[YY],invh[ZZ]);
 127
 128     calc_nxyz(nx,ny,nz,&nnx,&nny,&nnz);
 129
 130     for(i=0; (i<27); i++) {
 131       JCXYZ[3*i]   = 2 + (i/9);
 132       JCXYZ[3*i+1] = 2 + (i/3) % 3;
 133       JCXYZ[3*i+2] = 2 + (i % 3);
 134     }
 135
 136     bFirst = FALSE;
 137   }
 138
 139   energy=0.0;
 140   for(i=0; (i<natoms); i++) {
 141     /* Each charge is spread over the nearest 27 grid cells,
 142      * thus we loop over -1..1 in X,Y and Z direction
 143      * We apply the TSC (triangle shaped charge)
 144      * see Luty et. al, JCP 103 (1995) 3014
 145      */
 146
 147     calc_weights(i,nx,ny,nz,x[i],box,invh,ixyz,WXYZ);
 148
 149     for(ll=llim2; (ll<=ulim2); ll++) {
 150       ixw[ll-llim2] = nnx[ixyz[XX]+ll+nx];
 151       iyw[ll-llim2] = nny[ixyz[YY]+ll+ny];
 152       izw[ll-llim2] = nnz[ixyz[ZZ]+ll+nz];
 153     }
 154
 155     qi      = charge[i];
 156     pi      = ps_gather_inner(JCXYZ,WXYZ,ixw,iyw,izw,
 157                               c1x,c1y,c1z,c2x,c2y,c2z,
 158                               qi,f[i],ptr);
 159
 160     energy += pi*qi;
 161     pot[i]  = pi;
 162   }
 163
 164   inc_nrnb(nrnb,eNR_GATHERF,27*natoms);
 165   inc_nrnb(nrnb,eNR_WEIGHTS,3*natoms);
 166
 167   return energy*0.5;
 168 }
 169