src/mdlib/calcvir.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,2013, 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 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include "sysstuff.h"
  44 #include "force.h"
  45 #include "vec.h"
  46 #include "mshift.h"
  47 #include "macros.h"
  48
  49 static void upd_vir(rvec vir,real dvx,real dvy,real dvz)
  50 {
  51   vir[XX]-=0.5*dvx;
  52   vir[YY]-=0.5*dvy;
  53   vir[ZZ]-=0.5*dvz;
  54 }
  55
  56 void calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
  57               gmx_bool bScrewPBC,matrix box)
  58 {
  59   int      i,isx;
  60   double   dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
  61
  62   for(i=0; (i<nxf); i++) {
  63     dvxx+=x[i][XX]*f[i][XX];
  64     dvxy+=x[i][XX]*f[i][YY];
  65     dvxz+=x[i][XX]*f[i][ZZ];
  66
  67     dvyx+=x[i][YY]*f[i][XX];
  68     dvyy+=x[i][YY]*f[i][YY];
  69     dvyz+=x[i][YY]*f[i][ZZ];
  70
  71     dvzx+=x[i][ZZ]*f[i][XX];
  72     dvzy+=x[i][ZZ]*f[i][YY];
  73     dvzz+=x[i][ZZ]*f[i][ZZ];
  74
  75     if (bScrewPBC) {
  76       isx = IS2X(i);
  77       /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
  78       if (isx == 1 || isx == -1) {
  79         dvyy += box[YY][YY]*f[i][YY];
  80         dvyz += box[YY][YY]*f[i][ZZ];
  81
  82         dvzy += box[ZZ][ZZ]*f[i][YY];
  83         dvzz += box[ZZ][ZZ]*f[i][ZZ];
  84       }
  85     }
  86   }
  87
  88   upd_vir(vir[XX],dvxx,dvxy,dvxz);
  89   upd_vir(vir[YY],dvyx,dvyy,dvyz);
  90   upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
  91 }
  92
  93
  94 static void lo_fcv(int i0,int i1,
  95                    real x[],real f[],tensor vir,
  96                    int is[],real box[], gmx_bool bTriclinic)
  97 {
  98   int      i,i3,tx,ty,tz;
  99   real     xx,yy,zz;
 100   real     dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
 101
 102   if(bTriclinic) {
 103       for(i=i0; (i<i1); i++) {
 104           i3=DIM*i;
 105           tx=is[i3+XX];
 106           ty=is[i3+YY];
 107           tz=is[i3+ZZ];
 108
 109           xx=x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
 110           dvxx+=xx*f[i3+XX];
 111           dvxy+=xx*f[i3+YY];
 112           dvxz+=xx*f[i3+ZZ];
 113
 114           yy=x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
 115           dvyx+=yy*f[i3+XX];
 116           dvyy+=yy*f[i3+YY];
 117           dvyz+=yy*f[i3+ZZ];
 118
 119           zz=x[i3+ZZ]-tz*box[ZZZZ];
 120           dvzx+=zz*f[i3+XX];
 121           dvzy+=zz*f[i3+YY];
 122           dvzz+=zz*f[i3+ZZ];
 123       }
 124   } else {
 125       for(i=i0; (i<i1); i++) {
 126           i3=DIM*i;
 127           tx=is[i3+XX];
 128           ty=is[i3+YY];
 129           tz=is[i3+ZZ];
 130
 131           xx=x[i3+XX]-tx*box[XXXX];
 132           dvxx+=xx*f[i3+XX];
 133           dvxy+=xx*f[i3+YY];
 134           dvxz+=xx*f[i3+ZZ];
 135
 136           yy=x[i3+YY]-ty*box[YYYY];
 137           dvyx+=yy*f[i3+XX];
 138           dvyy+=yy*f[i3+YY];
 139           dvyz+=yy*f[i3+ZZ];
 140
 141           zz=x[i3+ZZ]-tz*box[ZZZZ];
 142           dvzx+=zz*f[i3+XX];
 143           dvzy+=zz*f[i3+YY];
 144           dvzz+=zz*f[i3+ZZ];
 145       }
 146   }
 147
 148   upd_vir(vir[XX],dvxx,dvxy,dvxz);
 149   upd_vir(vir[YY],dvyx,dvyy,dvyz);
 150   upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
 151 }
 152
 153 static void lo_fcv2(int i0,int i1,
 154                     rvec x[],rvec f[],tensor vir,
 155                     ivec is[],matrix box, gmx_bool bTriclinic)
 156 {
 157   int      i,gg,tx,ty,tz;
 158   real     xx,yy,zz;
 159   real     dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
 160
 161   if(bTriclinic) {
 162       for(i=i0,gg=0; (i<i1); i++,gg++) {
 163           tx=is[gg][XX];
 164           ty=is[gg][YY];
 165           tz=is[gg][ZZ];
 166
 167           xx=x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
 168           dvxx+=xx*f[i][XX];
 169           dvxy+=xx*f[i][YY];
 170           dvxz+=xx*f[i][ZZ];
 171
 172           yy=x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
 173           dvyx+=yy*f[i][XX];
 174           dvyy+=yy*f[i][YY];
 175           dvyz+=yy*f[i][ZZ];
 176
 177           zz=x[i][ZZ]-tz*box[ZZ][ZZ];
 178           dvzx+=zz*f[i][XX];
 179           dvzy+=zz*f[i][YY];
 180           dvzz+=zz*f[i][ZZ];
 181       }
 182   } else {
 183       for(i=i0,gg=0; (i<i1); i++,gg++) {
 184           tx=is[gg][XX];
 185           ty=is[gg][YY];
 186           tz=is[gg][ZZ];
 187
 188           xx=x[i][XX]-tx*box[XX][XX];
 189           dvxx+=xx*f[i][XX];
 190           dvxy+=xx*f[i][YY];
 191           dvxz+=xx*f[i][ZZ];
 192
 193           yy=x[i][YY]-ty*box[YY][YY];
 194           dvyx+=yy*f[i][XX];
 195           dvyy+=yy*f[i][YY];
 196           dvyz+=yy*f[i][ZZ];
 197
 198           zz=x[i][ZZ]-tz*box[ZZ][ZZ];
 199           dvzx+=zz*f[i][XX];
 200           dvzy+=zz*f[i][YY];
 201           dvzz+=zz*f[i][ZZ];
 202       }
 203   }
 204
 205   upd_vir(vir[XX],dvxx,dvxy,dvxz);
 206   upd_vir(vir[YY],dvyx,dvyy,dvyz);
 207   upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
 208 }
 209
 210 void f_calc_vir(FILE *log,int i0,int i1,rvec x[],rvec f[],tensor vir,
 211                 t_graph *g,matrix box)
 212 {
 213   int start,end;
 214
 215   if (g && (g->nnodes > 0)) {
 216     /* Calculate virial for bonded forces only when they belong to
 217      * this node.
 218      */
 219     start = max(i0,g->at_start);
 220     end   = min(i1,g->at_end);
 221 #ifdef SAFE
 222     lo_fcv2(start,end,x,f,vir,g->ishift,box,TRICLINIC(box));
 223 #else
 224     lo_fcv(start,end,x[0],f[0],vir,g->ishift[0],box[0],TRICLINIC(box));
 225 #endif
 226
 227     /* If not all atoms are bonded, calculate their virial contribution
 228      * anyway, without shifting back their coordinates.
 229      * Note the nifty pointer arithmetic...
 230      */
 231     if (start > i0)
 232       calc_vir(log,start-i0,x + i0,f + i0,vir,FALSE,box);
 233     if (end < i1)
 234       calc_vir(log,i1-end,x + end,f + end,vir,FALSE,box);
 235   }
 236   else
 237     calc_vir(log,i1-i0,x + i0,f + i0,vir,FALSE,box);
 238 }