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,2014, 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 #include "gmxpre.h"
  39
  40 #include "config.h"
  41
  42 #include "gromacs/legacyheaders/force.h"
  43 #include "gromacs/legacyheaders/macros.h"
  44
  45 #include "gromacs/math/vectypes.h"
  46 #include "gromacs/pbcutil/ishift.h"
  47 #include "gromacs/pbcutil/mshift.h"
  48 #include "gromacs/pbcutil/pbc.h"
  49
  50 #define XXXX    0
  51 #define XXYY    1
  52 #define XXZZ    2
  53 #define YYXX    3
  54 #define YYYY    4
  55 #define YYZZ    5
  56 #define ZZXX    6
  57 #define ZZYY    7
  58 #define ZZZZ    8
  59
  60 static void upd_vir(rvec vir, real dvx, real dvy, real dvz)
  61 {
  62     vir[XX] -= 0.5*dvx;
  63     vir[YY] -= 0.5*dvy;
  64     vir[ZZ] -= 0.5*dvz;
  65 }
  66
  67 void calc_vir(int nxf, rvec x[], rvec f[], tensor vir,
  68               gmx_bool bScrewPBC, matrix box)
  69 {
  70     int      i, isx;
  71     double   dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
  72
  73     for (i = 0; (i < nxf); i++)
  74     {
  75         dvxx += x[i][XX]*f[i][XX];
  76         dvxy += x[i][XX]*f[i][YY];
  77         dvxz += x[i][XX]*f[i][ZZ];
  78
  79         dvyx += x[i][YY]*f[i][XX];
  80         dvyy += x[i][YY]*f[i][YY];
  81         dvyz += x[i][YY]*f[i][ZZ];
  82
  83         dvzx += x[i][ZZ]*f[i][XX];
  84         dvzy += x[i][ZZ]*f[i][YY];
  85         dvzz += x[i][ZZ]*f[i][ZZ];
  86
  87         if (bScrewPBC)
  88         {
  89             isx = IS2X(i);
  90             /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
  91             if (isx == 1 || isx == -1)
  92             {
  93                 dvyy += box[YY][YY]*f[i][YY];
  94                 dvyz += box[YY][YY]*f[i][ZZ];
  95
  96                 dvzy += box[ZZ][ZZ]*f[i][YY];
  97                 dvzz += box[ZZ][ZZ]*f[i][ZZ];
  98             }
  99         }
 100     }
 101
 102     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 103     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 104     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 105 }
 106
 107
 108 static void lo_fcv(int i0, int i1,
 109                    real x[], real f[], tensor vir,
 110                    int is[], real box[], gmx_bool bTriclinic)
 111 {
 112     int      i, i3, tx, ty, tz;
 113     real     xx, yy, zz;
 114     real     dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
 115
 116     if (bTriclinic)
 117     {
 118         for (i = i0; (i < i1); i++)
 119         {
 120             i3 = DIM*i;
 121             tx = is[i3+XX];
 122             ty = is[i3+YY];
 123             tz = is[i3+ZZ];
 124
 125             xx    = x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
 126             dvxx += xx*f[i3+XX];
 127             dvxy += xx*f[i3+YY];
 128             dvxz += xx*f[i3+ZZ];
 129
 130             yy    = x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
 131             dvyx += yy*f[i3+XX];
 132             dvyy += yy*f[i3+YY];
 133             dvyz += yy*f[i3+ZZ];
 134
 135             zz    = x[i3+ZZ]-tz*box[ZZZZ];
 136             dvzx += zz*f[i3+XX];
 137             dvzy += zz*f[i3+YY];
 138             dvzz += zz*f[i3+ZZ];
 139         }
 140     }
 141     else
 142     {
 143         for (i = i0; (i < i1); i++)
 144         {
 145             i3 = DIM*i;
 146             tx = is[i3+XX];
 147             ty = is[i3+YY];
 148             tz = is[i3+ZZ];
 149
 150             xx    = x[i3+XX]-tx*box[XXXX];
 151             dvxx += xx*f[i3+XX];
 152             dvxy += xx*f[i3+YY];
 153             dvxz += xx*f[i3+ZZ];
 154
 155             yy    = x[i3+YY]-ty*box[YYYY];
 156             dvyx += yy*f[i3+XX];
 157             dvyy += yy*f[i3+YY];
 158             dvyz += yy*f[i3+ZZ];
 159
 160             zz    = x[i3+ZZ]-tz*box[ZZZZ];
 161             dvzx += zz*f[i3+XX];
 162             dvzy += zz*f[i3+YY];
 163             dvzz += zz*f[i3+ZZ];
 164         }
 165     }
 166
 167     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 168     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 169     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 170 }
 171
 172 static void lo_fcv2(int i0, int i1,
 173                     rvec x[], rvec f[], tensor vir,
 174                     ivec is[], matrix box, gmx_bool bTriclinic)
 175 {
 176     int      i, gg, tx, ty, tz;
 177     real     xx, yy, zz;
 178     real     dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
 179
 180     if (bTriclinic)
 181     {
 182         for (i = i0, gg = 0; (i < i1); i++, gg++)
 183         {
 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]-ty*box[YY][XX]-tz*box[ZZ][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]-tz*box[ZZ][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     else
 205     {
 206         for (i = i0, gg = 0; (i < i1); i++, gg++)
 207         {
 208             tx = is[gg][XX];
 209             ty = is[gg][YY];
 210             tz = is[gg][ZZ];
 211
 212             xx    = x[i][XX]-tx*box[XX][XX];
 213             dvxx += xx*f[i][XX];
 214             dvxy += xx*f[i][YY];
 215             dvxz += xx*f[i][ZZ];
 216
 217             yy    = x[i][YY]-ty*box[YY][YY];
 218             dvyx += yy*f[i][XX];
 219             dvyy += yy*f[i][YY];
 220             dvyz += yy*f[i][ZZ];
 221
 222             zz    = x[i][ZZ]-tz*box[ZZ][ZZ];
 223             dvzx += zz*f[i][XX];
 224             dvzy += zz*f[i][YY];
 225             dvzz += zz*f[i][ZZ];
 226         }
 227     }
 228
 229     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 230     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 231     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 232 }
 233
 234 void f_calc_vir(int i0, int i1, rvec x[], rvec f[], tensor vir,
 235                 t_graph *g, matrix box)
 236 {
 237     int start, end;
 238
 239     if (g && (g->nnodes > 0))
 240     {
 241         /* Calculate virial for bonded forces only when they belong to
 242          * this node.
 243          */
 244         start = max(i0, g->at_start);
 245         end   = min(i1, g->at_end);
 246 #ifdef SAFE
 247         lo_fcv2(start, end, x, f, vir, g->ishift, box, TRICLINIC(box));
 248 #else
 249         lo_fcv(start, end, x[0], f[0], vir, g->ishift[0], box[0], TRICLINIC(box));
 250 #endif
 251
 252         /* If not all atoms are bonded, calculate their virial contribution
 253          * anyway, without shifting back their coordinates.
 254          * Note the nifty pointer arithmetic...
 255          */
 256         if (start > i0)
 257         {
 258             calc_vir(start-i0, x + i0, f + i0, vir, FALSE, box);
 259         }
 260         if (end < i1)
 261         {
 262             calc_vir(i1-end, x + end, f + end, vir, FALSE, box);
 263         }
 264     }
 265     else
 266     {
 267         calc_vir(i1-i0, x + i0, f + i0, vir, FALSE, box);
 268     }
 269 }