src/mdlib/calcvir.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROwing Monsters And Cloning Shrimps
  34  */
  35 /* This file is completely threadsafe - keep it that way! */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include "sysstuff.h"
  41 #include "force.h"
  42 #include "vec.h"
  43 #include "mshift.h"
  44 #include "macros.h"
  45
  46 static void upd_vir(rvec vir,real dvx,real dvy,real dvz)
  47 {
  48   vir[XX]-=0.5*dvx;
  49   vir[YY]-=0.5*dvy;
  50   vir[ZZ]-=0.5*dvz;
  51 }
  52
  53 void calc_vir(FILE *log,int nxf,rvec x[],rvec f[],tensor vir,
  54               gmx_bool bScrewPBC,matrix box)
  55 {
  56   int      i,isx;
  57   double   dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
  58
  59   for(i=0; (i<nxf); i++) {
  60     dvxx+=x[i][XX]*f[i][XX];
  61     dvxy+=x[i][XX]*f[i][YY];
  62     dvxz+=x[i][XX]*f[i][ZZ];
  63
  64     dvyx+=x[i][YY]*f[i][XX];
  65     dvyy+=x[i][YY]*f[i][YY];
  66     dvyz+=x[i][YY]*f[i][ZZ];
  67
  68     dvzx+=x[i][ZZ]*f[i][XX];
  69     dvzy+=x[i][ZZ]*f[i][YY];
  70     dvzz+=x[i][ZZ]*f[i][ZZ];
  71
  72     if (bScrewPBC) {
  73       isx = IS2X(i);
  74       /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
  75       if (isx == 1 || isx == -1) {
  76         dvyy += box[YY][YY]*f[i][YY];
  77         dvyz += box[YY][YY]*f[i][ZZ];
  78
  79         dvzy += box[ZZ][ZZ]*f[i][YY];
  80         dvzz += box[ZZ][ZZ]*f[i][ZZ];
  81       }
  82     }
  83   }
  84
  85   upd_vir(vir[XX],dvxx,dvxy,dvxz);
  86   upd_vir(vir[YY],dvyx,dvyy,dvyz);
  87   upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
  88 }
  89
  90
  91 static void lo_fcv(int i0,int i1,
  92                    real x[],real f[],tensor vir,
  93                    int is[],real box[], gmx_bool bTriclinic)
  94 {
  95   int      i,i3,tx,ty,tz;
  96   real     xx,yy,zz;
  97   real     dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
  98
  99   if(bTriclinic) {
 100       for(i=i0; (i<i1); i++) {
 101           i3=DIM*i;
 102           tx=is[i3+XX];
 103           ty=is[i3+YY];
 104           tz=is[i3+ZZ];
 105
 106           xx=x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
 107           dvxx+=xx*f[i3+XX];
 108           dvxy+=xx*f[i3+YY];
 109           dvxz+=xx*f[i3+ZZ];
 110
 111           yy=x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
 112           dvyx+=yy*f[i3+XX];
 113           dvyy+=yy*f[i3+YY];
 114           dvyz+=yy*f[i3+ZZ];
 115
 116           zz=x[i3+ZZ]-tz*box[ZZZZ];
 117           dvzx+=zz*f[i3+XX];
 118           dvzy+=zz*f[i3+YY];
 119           dvzz+=zz*f[i3+ZZ];
 120       }
 121   } else {
 122       for(i=i0; (i<i1); i++) {
 123           i3=DIM*i;
 124           tx=is[i3+XX];
 125           ty=is[i3+YY];
 126           tz=is[i3+ZZ];
 127
 128           xx=x[i3+XX]-tx*box[XXXX];
 129           dvxx+=xx*f[i3+XX];
 130           dvxy+=xx*f[i3+YY];
 131           dvxz+=xx*f[i3+ZZ];
 132
 133           yy=x[i3+YY]-ty*box[YYYY];
 134           dvyx+=yy*f[i3+XX];
 135           dvyy+=yy*f[i3+YY];
 136           dvyz+=yy*f[i3+ZZ];
 137
 138           zz=x[i3+ZZ]-tz*box[ZZZZ];
 139           dvzx+=zz*f[i3+XX];
 140           dvzy+=zz*f[i3+YY];
 141           dvzz+=zz*f[i3+ZZ];
 142       }
 143   }
 144
 145   upd_vir(vir[XX],dvxx,dvxy,dvxz);
 146   upd_vir(vir[YY],dvyx,dvyy,dvyz);
 147   upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
 148 }
 149
 150 static void lo_fcv2(int i0,int i1,
 151                     rvec x[],rvec f[],tensor vir,
 152                     ivec is[],matrix box, gmx_bool bTriclinic)
 153 {
 154   int      i,gg,tx,ty,tz;
 155   real     xx,yy,zz;
 156   real     dvxx=0,dvxy=0,dvxz=0,dvyx=0,dvyy=0,dvyz=0,dvzx=0,dvzy=0,dvzz=0;
 157
 158   if(bTriclinic) {
 159       for(i=i0,gg=0; (i<i1); i++,gg++) {
 160           tx=is[gg][XX];
 161           ty=is[gg][YY];
 162           tz=is[gg][ZZ];
 163
 164           xx=x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
 165           dvxx+=xx*f[i][XX];
 166           dvxy+=xx*f[i][YY];
 167           dvxz+=xx*f[i][ZZ];
 168
 169           yy=x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
 170           dvyx+=yy*f[i][XX];
 171           dvyy+=yy*f[i][YY];
 172           dvyz+=yy*f[i][ZZ];
 173
 174           zz=x[i][ZZ]-tz*box[ZZ][ZZ];
 175           dvzx+=zz*f[i][XX];
 176           dvzy+=zz*f[i][YY];
 177           dvzz+=zz*f[i][ZZ];
 178       }
 179   } else {
 180       for(i=i0,gg=0; (i<i1); i++,gg++) {
 181           tx=is[gg][XX];
 182           ty=is[gg][YY];
 183           tz=is[gg][ZZ];
 184
 185           xx=x[i][XX]-tx*box[XX][XX];
 186           dvxx+=xx*f[i][XX];
 187           dvxy+=xx*f[i][YY];
 188           dvxz+=xx*f[i][ZZ];
 189
 190           yy=x[i][YY]-ty*box[YY][YY];
 191           dvyx+=yy*f[i][XX];
 192           dvyy+=yy*f[i][YY];
 193           dvyz+=yy*f[i][ZZ];
 194
 195           zz=x[i][ZZ]-tz*box[ZZ][ZZ];
 196           dvzx+=zz*f[i][XX];
 197           dvzy+=zz*f[i][YY];
 198           dvzz+=zz*f[i][ZZ];
 199       }
 200   }
 201
 202   upd_vir(vir[XX],dvxx,dvxy,dvxz);
 203   upd_vir(vir[YY],dvyx,dvyy,dvyz);
 204   upd_vir(vir[ZZ],dvzx,dvzy,dvzz);
 205 }
 206
 207 void f_calc_vir(FILE *log,int i0,int i1,rvec x[],rvec f[],tensor vir,
 208                 t_graph *g,matrix box)
 209 {
 210   int start,end;
 211
 212   if (g && (g->nnodes > 0)) {
 213     /* Calculate virial for bonded forces only when they belong to
 214      * this node.
 215      */
 216     start = max(i0,g->at_start);
 217     end   = min(i1,g->at_end);
 218 #ifdef SAFE
 219     lo_fcv2(start,end,x,f,vir,g->ishift,box,TRICLINIC(box));
 220 #else
 221     lo_fcv(start,end,x[0],f[0],vir,g->ishift[0],box[0],TRICLINIC(box));
 222 #endif
 223
 224     /* If not all atoms are bonded, calculate their virial contribution
 225      * anyway, without shifting back their coordinates.
 226      * Note the nifty pointer arithmetic...
 227      */
 228     if (start > i0)
 229       calc_vir(log,start-i0,x + i0,f + i0,vir,FALSE,box);
 230     if (end < i1)
 231       calc_vir(log,i1-end,x + end,f + end,vir,FALSE,box);
 232   }
 233   else
 234     calc_vir(log,i1-i0,x + i0,f + i0,vir,FALSE,box);
 235 }