src/mdlib/ewald.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  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_ewald_c = "$Id$";
  30
  31 #include <stdio.h>
  32 #include <math.h>
  33 #include "typedefs.h"
  34 #include "vec.h"
  35 #include "complex.h"
  36 #include "smalloc.h"
  37 #include "futil.h"
  38 #include "ewald_util.h"
  39 #include "shift_util.h"
  40 #include "fftgrid.h"
  41 #include "fatal.h"
  42 #include "physics.h"
  43
  44 #define TOL 2e-5
  45
  46
  47
  48
  49
  50 /* the other routines are in complex.h */
  51 static t_complex conjmul(t_complex a,t_complex b)
  52 {
  53   t_complex c;
  54
  55   c.re = a.re*b.re + a.im*b.im;
  56   c.im = a.im*b.re - a.re*b.im;
  57
  58   return c;
  59 }
  60
  61
  62
  63
  64 static void tabulate_eir(int natom,rvec x[],int kmax,cvec **eir,rvec lll)
  65 {
  66   int  i,j,m;
  67
  68   if (kmax < 1) {
  69       printf("Go away! kmax = %d\n",kmax);
  70       exit(1);
  71   }
  72
  73   for(i=0; (i<natom); i++) {
  74     for(m=0; (m<3); m++) {
  75       eir[0][i][m].re = 1;
  76       eir[0][i][m].im = 0;
  77     }
  78
  79     for(m=0; (m<3); m++) {
  80       eir[1][i][m].re = cos(x[i][m]*lll[m]);
  81       eir[1][i][m].im = sin(x[i][m]*lll[m]);
  82     }
  83     for(j=2; (j<kmax); j++)
  84       for(m=0; (m<3); m++)
  85         eir[j][i][m] = cmul(eir[j-1][i][m],eir[1][i][m]);
  86   }
  87 }
  88
  89
  90
  91
  92 real do_ewald(FILE *log,       bool bVerbose,
  93               t_inputrec *ir,
  94               rvec x[],        rvec f[],
  95                   real charge[],   rvec box,
  96               t_commrec *cr,          t_nsborder *nsb,
  97               matrix lrvir, real ewaldcoeff)
  98 {
  99   static    bool bFirst = TRUE;
 100   static    int       nx,ny,nz,kmax;
 101   static    cvec      **eir;
 102   static    t_complex  *tab_xy,*tab_qxyz;
 103   real factor=-1.0/(4*ewaldcoeff*ewaldcoeff);
 104   real energy;
 105   rvec lll;
 106   int  lowiy,lowiz,ix,iy,iz,n;
 107   real tmp,cs,ss,ak,akv,mx,my,mz,m2;
 108
 109   if (bFirst) {
 110       if (bVerbose)
 111     fprintf(log,"Will do ordinary reciprocal space Ewald sum.\n");
 112
 113     if (cr != NULL) {
 114       if (cr->nprocs > 1)
 115         fatal_error(0,"No parallel Ewald. Use PME instead.\n");
 116     }
 117
 118     nx = ir->nkx+1;
 119     ny = ir->nky+1;
 120     nz = ir->nkz+1;
 121     kmax = max(nx,max(ny,nz));
 122     snew(eir,kmax);
 123     for(n=0;n<kmax;n++)
 124       snew(eir[n],HOMENR(nsb));
 125     snew(tab_xy,HOMENR(nsb));
 126     snew(tab_qxyz,HOMENR(nsb));
 127     bFirst = FALSE;
 128   }
 129   clear_mat(lrvir);
 130
 131   calc_lll(box,lll);
 132   /* make tables for the structure factor parts */
 133   tabulate_eir(HOMENR(nsb),x,kmax,eir,lll);
 134
 135   lowiy=0;
 136   lowiz=1;
 137   energy=0;
 138   for(ix=0;ix<nx;ix++) {
 139     mx=ix*lll[XX];
 140     for(iy=lowiy;iy<ny;iy++) {
 141       my=iy*lll[YY];
 142       if(iy>=0)
 143         for(n=0;n<HOMENR(nsb);n++)
 144           tab_xy[n]=cmul(eir[ix][n][XX],eir[iy][n][YY]);
 145       else
 146         for(n=0;n<HOMENR(nsb);n++)
 147           tab_xy[n]=conjmul(eir[ix][n][XX],eir[-iy][n][YY]);
 148       for(iz=lowiz;iz<nz;iz++) {
 149         mz=iz*lll[ZZ];
 150         m2=mx*mx+my*my+mz*mz;
 151         ak=exp(m2*factor)/m2;
 152         akv=2.0*ak*(1.0/m2-factor);
 153         if(iz>=0)
 154           for(n=0;n<HOMENR(nsb);n++)
 155             tab_qxyz[n]=rcmul(charge[n],cmul(tab_xy[n],eir[iz][n][ZZ]));
 156         else
 157           for(n=0;n<HOMENR(nsb);n++)
 158             tab_qxyz[n]=rcmul(charge[n],conjmul(tab_xy[n],eir[-iz][n][ZZ]));
 159
 160         cs=ss=0;
 161         for(n=0;n<HOMENR(nsb);n++) {
 162           cs+=tab_qxyz[n].re;
 163           ss+=tab_qxyz[n].im;
 164         }
 165         energy+=ak*(cs*cs+ss*ss);
 166         tmp=akv*(cs*cs+ss*ss);
 167         lrvir[XX][XX]-=tmp*mx*mx;
 168         lrvir[XX][YY]-=tmp*mx*my;
 169         lrvir[XX][ZZ]-=tmp*mx*mz;
 170         lrvir[YY][YY]-=tmp*my*my;
 171         lrvir[YY][ZZ]-=tmp*my*mz;
 172         lrvir[ZZ][ZZ]-=tmp*mz*mz;
 173         for(n=0;n<HOMENR(nsb);n++) {
 174           tmp=ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);
 175           f[n][XX]+=tmp*mx;
 176           f[n][YY]+=tmp*my;
 177           f[n][ZZ]+=tmp*mz;
 178         }
 179         lowiz=1-nz;
 180       }
 181       lowiy=1-ny;
 182     }
 183   }
 184   tmp=4.0*M_PI/(box[XX]*box[YY]*box[ZZ])*ONE_4PI_EPS0;
 185   for(n=0;n<HOMENR(nsb);n++) {
 186     f[n][XX]*=2*tmp;
 187     f[n][YY]*=2*tmp;
 188     f[n][ZZ]*=2*tmp;
 189   }
 190   lrvir[XX][XX]=-0.5*tmp*(lrvir[XX][XX]+energy);
 191   lrvir[XX][YY]=-0.5*tmp*(lrvir[XX][YY]);
 192   lrvir[XX][ZZ]=-0.5*tmp*(lrvir[XX][ZZ]);
 193   lrvir[YY][YY]=-0.5*tmp*(lrvir[YY][YY]+energy);
 194   lrvir[YY][ZZ]=-0.5*tmp*(lrvir[YY][ZZ]);
 195   lrvir[ZZ][ZZ]=-0.5*tmp*(lrvir[ZZ][ZZ]+energy);
 196
 197   lrvir[YY][XX]=lrvir[XX][YY];
 198   lrvir[ZZ][XX]=lrvir[XX][ZZ];
 199   lrvir[ZZ][YY]=lrvir[YY][ZZ];
 200
 201   energy*=tmp;
 202
 203   return energy;
 204 }