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.2.0
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * GROwing Monsters And Cloning Shrimps
  35  */
  36 #include <stdio.h>
  37 #include <math.h>
  38 #include "typedefs.h"
  39 #include "vec.h"
  40 #include "gmxcomplex.h"
  41 #include "smalloc.h"
  42 #include "futil.h"
  43 #include "ewald_util.h"
  44 #include "shift_util.h"
  45 #include "fftgrid.h"
  46 #include "fatal.h"
  47 #include "physics.h"
  48
  49 #define TOL 2e-5
  50
  51
  52
  53
  54
  55 /* the other routines are in complex.h */
  56 static t_complex conjmul(t_complex a,t_complex b)
  57 {
  58   t_complex c;
  59
  60   c.re = a.re*b.re + a.im*b.im;
  61   c.im = a.im*b.re - a.re*b.im;
  62
  63   return c;
  64 }
  65
  66
  67
  68
  69 static void tabulate_eir(int natom,rvec x[],int kmax,cvec **eir,rvec lll)
  70 {
  71   int  i,j,m;
  72
  73   if (kmax < 1) {
  74       printf("Go away! kmax = %d\n",kmax);
  75       exit(1);
  76   }
  77
  78   for(i=0; (i<natom); i++) {
  79     for(m=0; (m<3); m++) {
  80       eir[0][i][m].re = 1;
  81       eir[0][i][m].im = 0;
  82     }
  83
  84     for(m=0; (m<3); m++) {
  85       eir[1][i][m].re = cos(x[i][m]*lll[m]);
  86       eir[1][i][m].im = sin(x[i][m]*lll[m]);
  87     }
  88     for(j=2; (j<kmax); j++)
  89       for(m=0; (m<3); m++)
  90         eir[j][i][m] = cmul(eir[j-1][i][m],eir[1][i][m]);
  91   }
  92 }
  93
  94
  95
  96
  97 real do_ewald(FILE *log,       bool bVerbose,
  98               t_inputrec *ir,
  99               rvec x[],        rvec f[],
 100                   real charge[],   rvec box,
 101               t_commrec *cr,          t_nsborder *nsb,
 102               matrix lrvir, real ewaldcoeff)
 103 {
 104   static    bool bFirst = TRUE;
 105   static    int       nx,ny,nz,kmax;
 106   static    cvec      **eir;
 107   static    t_complex  *tab_xy,*tab_qxyz;
 108   real factor=-1.0/(4*ewaldcoeff*ewaldcoeff);
 109   real energy;
 110   rvec lll;
 111   int  lowiy,lowiz,ix,iy,iz,n;
 112   real tmp,cs,ss,ak,akv,mx,my,mz,m2;
 113
 114   if (bFirst) {
 115       if (bVerbose)
 116     fprintf(log,"Will do ordinary reciprocal space Ewald sum.\n");
 117
 118     if (cr != NULL) {
 119       if (cr->nnodes > 1 || cr->nthreads>1)
 120         fatal_error(0,"No parallel Ewald. Use PME instead.\n");
 121     }
 122
 123     nx = ir->nkx+1;
 124     ny = ir->nky+1;
 125     nz = ir->nkz+1;
 126     kmax = max(nx,max(ny,nz));
 127     snew(eir,kmax);
 128     for(n=0;n<kmax;n++)
 129       snew(eir[n],HOMENR(nsb));
 130     snew(tab_xy,HOMENR(nsb));
 131     snew(tab_qxyz,HOMENR(nsb));
 132     bFirst = FALSE;
 133   }
 134   clear_mat(lrvir);
 135
 136   calc_lll(box,lll);
 137   /* make tables for the structure factor parts */
 138   tabulate_eir(HOMENR(nsb),x,kmax,eir,lll);
 139
 140   lowiy=0;
 141   lowiz=1;
 142   energy=0;
 143   for(ix=0;ix<nx;ix++) {
 144     mx=ix*lll[XX];
 145     for(iy=lowiy;iy<ny;iy++) {
 146       my=iy*lll[YY];
 147       if(iy>=0)
 148         for(n=0;n<HOMENR(nsb);n++)
 149           tab_xy[n]=cmul(eir[ix][n][XX],eir[iy][n][YY]);
 150       else
 151         for(n=0;n<HOMENR(nsb);n++)
 152           tab_xy[n]=conjmul(eir[ix][n][XX],eir[-iy][n][YY]);
 153       for(iz=lowiz;iz<nz;iz++) {
 154         mz=iz*lll[ZZ];
 155         m2=mx*mx+my*my+mz*mz;
 156         ak=exp(m2*factor)/m2;
 157         akv=2.0*ak*(1.0/m2-factor);
 158         if(iz>=0)
 159           for(n=0;n<HOMENR(nsb);n++)
 160             tab_qxyz[n]=rcmul(charge[n],cmul(tab_xy[n],eir[iz][n][ZZ]));
 161         else
 162           for(n=0;n<HOMENR(nsb);n++)
 163             tab_qxyz[n]=rcmul(charge[n],conjmul(tab_xy[n],eir[-iz][n][ZZ]));
 164
 165         cs=ss=0;
 166         for(n=0;n<HOMENR(nsb);n++) {
 167           cs+=tab_qxyz[n].re;
 168           ss+=tab_qxyz[n].im;
 169         }
 170         energy+=ak*(cs*cs+ss*ss);
 171         tmp=akv*(cs*cs+ss*ss);
 172         lrvir[XX][XX]-=tmp*mx*mx;
 173         lrvir[XX][YY]-=tmp*mx*my;
 174         lrvir[XX][ZZ]-=tmp*mx*mz;
 175         lrvir[YY][YY]-=tmp*my*my;
 176         lrvir[YY][ZZ]-=tmp*my*mz;
 177         lrvir[ZZ][ZZ]-=tmp*mz*mz;
 178         for(n=0;n<HOMENR(nsb);n++) {
 179           tmp=ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);
 180           f[n][XX]+=tmp*mx;
 181           f[n][YY]+=tmp*my;
 182           f[n][ZZ]+=tmp*mz;
 183         }
 184         lowiz=1-nz;
 185       }
 186       lowiy=1-ny;
 187     }
 188   }
 189   tmp=4.0*M_PI/(box[XX]*box[YY]*box[ZZ])*ONE_4PI_EPS0;
 190   for(n=0;n<HOMENR(nsb);n++) {
 191     f[n][XX]*=2*tmp;
 192     f[n][YY]*=2*tmp;
 193     f[n][ZZ]*=2*tmp;
 194   }
 195   lrvir[XX][XX]=-0.5*tmp*(lrvir[XX][XX]+energy);
 196   lrvir[XX][YY]=-0.5*tmp*(lrvir[XX][YY]);
 197   lrvir[XX][ZZ]=-0.5*tmp*(lrvir[XX][ZZ]);
 198   lrvir[YY][YY]=-0.5*tmp*(lrvir[YY][YY]+energy);
 199   lrvir[YY][ZZ]=-0.5*tmp*(lrvir[YY][ZZ]);
 200   lrvir[ZZ][ZZ]=-0.5*tmp*(lrvir[ZZ][ZZ]+energy);
 201
 202   lrvir[YY][XX]=lrvir[XX][YY];
 203   lrvir[ZZ][XX]=lrvir[XX][ZZ];
 204   lrvir[ZZ][YY]=lrvir[YY][ZZ];
 205
 206   energy*=tmp;
 207
 208   return energy;
 209 }