src/mdlib/ewald.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <stdio.h>
  43 #include <math.h>
  44 #include "typedefs.h"
  45 #include "vec.h"
  46 #include "gmxcomplex.h"
  47 #include "smalloc.h"
  48 #include "futil.h"
  49 #include "gmx_fatal.h"
  50 #include "physics.h"
  51 #include "coulomb.h"
  52
  53 #define TOL 2e-5
  54
  55 struct ewald_tab
  56 {
  57     int nx,ny,nz,kmax;
  58     cvec **eir;
  59     t_complex *tab_xy, *tab_qxyz;
  60 };
  61
  62
  63
  64 /* TODO: fix thread-safety */
  65
  66 /* the other routines are in complex.h */
  67 static t_complex conjmul(t_complex a,t_complex b)
  68 {
  69   t_complex c;
  70
  71   c.re = a.re*b.re + a.im*b.im;
  72   c.im = a.im*b.re - a.re*b.im;
  73
  74   return c;
  75 }
  76
  77
  78
  79
  80 static void tabulate_eir(int natom,rvec x[],int kmax,cvec **eir,rvec lll)
  81 {
  82   int  i,j,m;
  83
  84   if (kmax < 1) {
  85       printf("Go away! kmax = %d\n",kmax);
  86       exit(1);
  87   }
  88
  89   for(i=0; (i<natom); i++) {
  90     for(m=0; (m<3); m++) {
  91       eir[0][i][m].re = 1;
  92       eir[0][i][m].im = 0;
  93     }
  94
  95     for(m=0; (m<3); m++) {
  96       eir[1][i][m].re = cos(x[i][m]*lll[m]);
  97       eir[1][i][m].im = sin(x[i][m]*lll[m]);
  98     }
  99     for(j=2; (j<kmax); j++)
 100       for(m=0; (m<3); m++)
 101         eir[j][i][m] = cmul(eir[j-1][i][m],eir[1][i][m]);
 102   }
 103 }
 104
 105 void init_ewald_tab(ewald_tab_t *et, const t_commrec *cr, const t_inputrec *ir,
 106                     FILE *fp)
 107 {
 108     int n;
 109
 110     snew(*et, 1);
 111     if (fp)
 112         fprintf(fp,"Will do ordinary reciprocal space Ewald sum.\n");
 113
 114     (*et)->nx = ir->nkx+1;
 115     (*et)->ny = ir->nky+1;
 116     (*et)->nz = ir->nkz+1;
 117     (*et)->kmax = max((*et)->nx,max((*et)->ny,(*et)->nz));
 118     (*et)->eir = NULL;
 119     (*et)->tab_xy = NULL;
 120     (*et)->tab_qxyz = NULL;
 121 }
 122
 123
 124
 125 real do_ewald(FILE *log,       gmx_bool bVerbose,
 126               t_inputrec *ir,
 127               rvec x[],        rvec f[],
 128               real chargeA[],  real chargeB[],
 129               rvec box,
 130               t_commrec *cr,   int natoms,
 131               matrix lrvir,    real ewaldcoeff,
 132               real lambda,     real *dvdlambda,
 133               ewald_tab_t et)
 134 {
 135   real factor=-1.0/(4*ewaldcoeff*ewaldcoeff);
 136   real scaleRecip =4.0*M_PI/(box[XX]*box[YY]*box[ZZ])*ONE_4PI_EPS0/ir->epsilon_r; /* 1/(Vol*e0) */
 137   real *charge,energy_AB[2],energy;
 138   rvec lll;
 139   int  lowiy,lowiz,ix,iy,iz,n,q;
 140   real tmp,cs,ss,ak,akv,mx,my,mz,m2,scale;
 141   gmx_bool bFreeEnergy;
 142
 143     if (cr != NULL)
 144     {
 145         if (PAR(cr))
 146         {
 147             gmx_fatal(FARGS,"No parallel Ewald. Use PME instead.\n");
 148         }
 149     }
 150
 151
 152   if (!et->eir) /* allocate if we need to */
 153   {
 154       snew(et->eir,et->kmax);
 155       for(n=0;n<et->kmax;n++)
 156           snew(et->eir[n],natoms);
 157       snew(et->tab_xy,natoms);
 158       snew(et->tab_qxyz,natoms);
 159   }
 160
 161   bFreeEnergy = (ir->efep != efepNO);
 162
 163   clear_mat(lrvir);
 164
 165   calc_lll(box,lll);
 166   /* make tables for the structure factor parts */
 167   tabulate_eir(natoms,x,et->kmax,et->eir,lll);
 168
 169   for(q=0; q<(bFreeEnergy ? 2 : 1); q++) {
 170     if (!bFreeEnergy) {
 171       charge = chargeA;
 172       scale = 1.0;
 173     } else if (q==0) {
 174       charge = chargeA;
 175       scale = 1.0 - lambda;
 176     } else {
 177       charge = chargeB;
 178       scale = lambda;
 179     }
 180     lowiy=0;
 181     lowiz=1;
 182     energy_AB[q]=0;
 183     for(ix=0;ix<et->nx;ix++) {
 184       mx=ix*lll[XX];
 185       for(iy=lowiy;iy<et->ny;iy++) {
 186         my=iy*lll[YY];
 187         if(iy>=0)
 188           for(n=0;n<natoms;n++)
 189             et->tab_xy[n]=cmul(et->eir[ix][n][XX],et->eir[iy][n][YY]);
 190         else
 191           for(n=0;n<natoms;n++)
 192             et->tab_xy[n]=conjmul(et->eir[ix][n][XX],et->eir[-iy][n][YY]);
 193         for(iz=lowiz;iz<et->nz;iz++) {
 194           mz=iz*lll[ZZ];
 195           m2=mx*mx+my*my+mz*mz;
 196           ak=exp(m2*factor)/m2;
 197           akv=2.0*ak*(1.0/m2-factor);
 198           if(iz>=0)
 199             for(n=0;n<natoms;n++)
 200               et->tab_qxyz[n]=rcmul(charge[n],cmul(et->tab_xy[n],
 201                                                    et->eir[iz][n][ZZ]));
 202           else
 203             for(n=0;n<natoms;n++)
 204               et->tab_qxyz[n]=rcmul(charge[n],conjmul(et->tab_xy[n],
 205                                                       et->eir[-iz][n][ZZ]));
 206
 207           cs=ss=0;
 208           for(n=0;n<natoms;n++) {
 209             cs+=et->tab_qxyz[n].re;
 210             ss+=et->tab_qxyz[n].im;
 211           }
 212           energy_AB[q]+=ak*(cs*cs+ss*ss);
 213           tmp=scale*akv*(cs*cs+ss*ss);
 214           lrvir[XX][XX]-=tmp*mx*mx;
 215           lrvir[XX][YY]-=tmp*mx*my;
 216           lrvir[XX][ZZ]-=tmp*mx*mz;
 217           lrvir[YY][YY]-=tmp*my*my;
 218           lrvir[YY][ZZ]-=tmp*my*mz;
 219           lrvir[ZZ][ZZ]-=tmp*mz*mz;
 220           for(n=0;n<natoms;n++) {
 221             /*tmp=scale*ak*(cs*tab_qxyz[n].im-ss*tab_qxyz[n].re);*/
 222             tmp=scale*ak*(cs*et->tab_qxyz[n].im-ss*et->tab_qxyz[n].re);
 223             f[n][XX]+=tmp*mx*2*scaleRecip;
 224             f[n][YY]+=tmp*my*2*scaleRecip;
 225             f[n][ZZ]+=tmp*mz*2*scaleRecip;
 226 #if 0
 227             f[n][XX]+=tmp*mx;
 228             f[n][YY]+=tmp*my;
 229             f[n][ZZ]+=tmp*mz;
 230 #endif
 231           }
 232           lowiz=1-et->nz;
 233         }
 234         lowiy=1-et->ny;
 235       }
 236     }
 237   }
 238
 239   if (!bFreeEnergy) {
 240     energy = energy_AB[0];
 241   } else {
 242     energy = (1.0 - lambda)*energy_AB[0] + lambda*energy_AB[1];
 243     *dvdlambda += scaleRecip*(energy_AB[1] - energy_AB[0]);
 244   }
 245
 246   lrvir[XX][XX]=-0.5*scaleRecip*(lrvir[XX][XX]+energy);
 247   lrvir[XX][YY]=-0.5*scaleRecip*(lrvir[XX][YY]);
 248   lrvir[XX][ZZ]=-0.5*scaleRecip*(lrvir[XX][ZZ]);
 249   lrvir[YY][YY]=-0.5*scaleRecip*(lrvir[YY][YY]+energy);
 250   lrvir[YY][ZZ]=-0.5*scaleRecip*(lrvir[YY][ZZ]);
 251   lrvir[ZZ][ZZ]=-0.5*scaleRecip*(lrvir[ZZ][ZZ]+energy);
 252
 253   lrvir[YY][XX]=lrvir[XX][YY];
 254   lrvir[ZZ][XX]=lrvir[XX][ZZ];
 255   lrvir[ZZ][YY]=lrvir[YY][ZZ];
 256
 257   energy*=scaleRecip;
 258
 259   return energy;
 260 }