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