src/gromacs/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  * Copyright (c) 2013, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /* This file is completely threadsafe - keep it that way! */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "sysstuff.h"
  43 #include "force.h"
  44 #include "vec.h"
  45 #include "mshift.h"
  46 #include "macros.h"
  47
  48 static void upd_vir(rvec vir, real dvx, real dvy, real dvz)
  49 {
  50     vir[XX] -= 0.5*dvx;
  51     vir[YY] -= 0.5*dvy;
  52     vir[ZZ] -= 0.5*dvz;
  53 }
  54
  55 void calc_vir(int nxf, rvec x[], rvec f[], tensor vir,
  56               gmx_bool bScrewPBC, matrix box)
  57 {
  58     int      i, isx;
  59     double   dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
  60
  61     for (i = 0; (i < nxf); i++)
  62     {
  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         {
  77             isx = IS2X(i);
  78             /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
  79             if (isx == 1 || isx == -1)
  80             {
  81                 dvyy += box[YY][YY]*f[i][YY];
  82                 dvyz += box[YY][YY]*f[i][ZZ];
  83
  84                 dvzy += box[ZZ][ZZ]*f[i][YY];
  85                 dvzz += box[ZZ][ZZ]*f[i][ZZ];
  86             }
  87         }
  88     }
  89
  90     upd_vir(vir[XX], dvxx, dvxy, dvxz);
  91     upd_vir(vir[YY], dvyx, dvyy, dvyz);
  92     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
  93 }
  94
  95
  96 static void lo_fcv(int i0, int i1,
  97                    real x[], real f[], tensor vir,
  98                    int is[], real box[], gmx_bool bTriclinic)
  99 {
 100     int      i, i3, tx, ty, tz;
 101     real     xx, yy, zz;
 102     real     dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
 103
 104     if (bTriclinic)
 105     {
 106         for (i = i0; (i < i1); i++)
 107         {
 108             i3 = DIM*i;
 109             tx = is[i3+XX];
 110             ty = is[i3+YY];
 111             tz = is[i3+ZZ];
 112
 113             xx    = x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
 114             dvxx += xx*f[i3+XX];
 115             dvxy += xx*f[i3+YY];
 116             dvxz += xx*f[i3+ZZ];
 117
 118             yy    = x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
 119             dvyx += yy*f[i3+XX];
 120             dvyy += yy*f[i3+YY];
 121             dvyz += yy*f[i3+ZZ];
 122
 123             zz    = x[i3+ZZ]-tz*box[ZZZZ];
 124             dvzx += zz*f[i3+XX];
 125             dvzy += zz*f[i3+YY];
 126             dvzz += zz*f[i3+ZZ];
 127         }
 128     }
 129     else
 130     {
 131         for (i = i0; (i < i1); i++)
 132         {
 133             i3 = DIM*i;
 134             tx = is[i3+XX];
 135             ty = is[i3+YY];
 136             tz = is[i3+ZZ];
 137
 138             xx    = x[i3+XX]-tx*box[XXXX];
 139             dvxx += xx*f[i3+XX];
 140             dvxy += xx*f[i3+YY];
 141             dvxz += xx*f[i3+ZZ];
 142
 143             yy    = x[i3+YY]-ty*box[YYYY];
 144             dvyx += yy*f[i3+XX];
 145             dvyy += yy*f[i3+YY];
 146             dvyz += yy*f[i3+ZZ];
 147
 148             zz    = x[i3+ZZ]-tz*box[ZZZZ];
 149             dvzx += zz*f[i3+XX];
 150             dvzy += zz*f[i3+YY];
 151             dvzz += zz*f[i3+ZZ];
 152         }
 153     }
 154
 155     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 156     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 157     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 158 }
 159
 160 static void lo_fcv2(int i0, int i1,
 161                     rvec x[], rvec f[], tensor vir,
 162                     ivec is[], matrix box, gmx_bool bTriclinic)
 163 {
 164     int      i, gg, tx, ty, tz;
 165     real     xx, yy, zz;
 166     real     dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
 167
 168     if (bTriclinic)
 169     {
 170         for (i = i0, gg = 0; (i < i1); i++, gg++)
 171         {
 172             tx = is[gg][XX];
 173             ty = is[gg][YY];
 174             tz = is[gg][ZZ];
 175
 176             xx    = x[i][XX]-tx*box[XX][XX]-ty*box[YY][XX]-tz*box[ZZ][XX];
 177             dvxx += xx*f[i][XX];
 178             dvxy += xx*f[i][YY];
 179             dvxz += xx*f[i][ZZ];
 180
 181             yy    = x[i][YY]-ty*box[YY][YY]-tz*box[ZZ][YY];
 182             dvyx += yy*f[i][XX];
 183             dvyy += yy*f[i][YY];
 184             dvyz += yy*f[i][ZZ];
 185
 186             zz    = x[i][ZZ]-tz*box[ZZ][ZZ];
 187             dvzx += zz*f[i][XX];
 188             dvzy += zz*f[i][YY];
 189             dvzz += zz*f[i][ZZ];
 190         }
 191     }
 192     else
 193     {
 194         for (i = i0, gg = 0; (i < i1); i++, gg++)
 195         {
 196             tx = is[gg][XX];
 197             ty = is[gg][YY];
 198             tz = is[gg][ZZ];
 199
 200             xx    = x[i][XX]-tx*box[XX][XX];
 201             dvxx += xx*f[i][XX];
 202             dvxy += xx*f[i][YY];
 203             dvxz += xx*f[i][ZZ];
 204
 205             yy    = x[i][YY]-ty*box[YY][YY];
 206             dvyx += yy*f[i][XX];
 207             dvyy += yy*f[i][YY];
 208             dvyz += yy*f[i][ZZ];
 209
 210             zz    = x[i][ZZ]-tz*box[ZZ][ZZ];
 211             dvzx += zz*f[i][XX];
 212             dvzy += zz*f[i][YY];
 213             dvzz += zz*f[i][ZZ];
 214         }
 215     }
 216
 217     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 218     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 219     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 220 }
 221
 222 void f_calc_vir(int i0, int i1, rvec x[], rvec f[], tensor vir,
 223                 t_graph *g, matrix box)
 224 {
 225     int start, end;
 226
 227     if (g && (g->nnodes > 0))
 228     {
 229         /* Calculate virial for bonded forces only when they belong to
 230          * this node.
 231          */
 232         start = max(i0, g->at_start);
 233         end   = min(i1, g->at_end);
 234 #ifdef SAFE
 235         lo_fcv2(start, end, x, f, vir, g->ishift, box, TRICLINIC(box));
 236 #else
 237         lo_fcv(start, end, x[0], f[0], vir, g->ishift[0], box[0], TRICLINIC(box));
 238 #endif
 239
 240         /* If not all atoms are bonded, calculate their virial contribution
 241          * anyway, without shifting back their coordinates.
 242          * Note the nifty pointer arithmetic...
 243          */
 244         if (start > i0)
 245         {
 246             calc_vir(start-i0, x + i0, f + i0, vir, FALSE, box);
 247         }
 248         if (end < i1)
 249         {
 250             calc_vir(i1-end, x + end, f + end, vir, FALSE, box);
 251         }
 252     }
 253     else
 254     {
 255         calc_vir(i1-i0, x + i0, f + i0, vir, FALSE, box);
 256     }
 257 }