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