src/gromacs/mdlib/calcvir.cpp

   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,2015,2016, 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" /* for GMX_MAX_OPENMP_THREADS */
  41
  42 #include <algorithm>
  43
  44 #include "gromacs/math/vec.h"
  45 #include "gromacs/math/vectypes.h"
  46 #include "gromacs/mdlib/force.h"
  47 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  48 #include "gromacs/pbcutil/ishift.h"
  49 #include "gromacs/pbcutil/mshift.h"
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/utility/gmxassert.h"
  52
  53 #define XXXX    0
  54 #define XXYY    1
  55 #define XXZZ    2
  56 #define YYXX    3
  57 #define YYYY    4
  58 #define YYZZ    5
  59 #define ZZXX    6
  60 #define ZZYY    7
  61 #define ZZZZ    8
  62
  63 static void upd_vir(rvec vir, real dvx, real dvy, real dvz)
  64 {
  65     vir[XX] -= 0.5*dvx;
  66     vir[YY] -= 0.5*dvy;
  67     vir[ZZ] -= 0.5*dvz;
  68 }
  69
  70 static void calc_x_times_f(int nxf, const rvec x[], const rvec f[],
  71                            gmx_bool bScrewPBC, const matrix box,
  72                            matrix x_times_f)
  73 {
  74     clear_mat(x_times_f);
  75
  76     for (int i = 0; i < nxf; i++)
  77     {
  78         for (int d = 0; d < DIM; d++)
  79         {
  80             for (int n = 0; n < DIM; n++)
  81             {
  82                 x_times_f[d][n] += x[i][d]*f[i][n];
  83             }
  84         }
  85
  86         if (bScrewPBC)
  87         {
  88             int isx = IS2X(i);
  89             /* We should correct all odd x-shifts, but the range of isx is -2 to 2 */
  90             if (isx == 1 || isx == -1)
  91             {
  92                 for (int d = 0; d < DIM; d++)
  93                 {
  94                     for (int n = 0; n < DIM; n++)
  95                     {
  96                         x_times_f[d][n] += box[d][d]*f[i][n];
  97                     }
  98                 }
  99             }
 100         }
 101     }
 102 }
 103
 104 void calc_vir(int nxf, rvec x[], rvec f[], tensor vir,
 105               gmx_bool bScrewPBC, matrix box)
 106 {
 107     matrix x_times_f;
 108
 109     int    nthreads = gmx_omp_nthreads_get_simple_rvec_task(emntDefault, nxf*9);
 110
 111     GMX_ASSERT(nthreads >= 1, "Avoids uninitialized x_times_f (warning)");
 112
 113     if (nthreads == 1)
 114     {
 115         calc_x_times_f(nxf, x, f, bScrewPBC, box, x_times_f);
 116     }
 117     else
 118     {
 119         /* Use a buffer on the stack for storing thread-local results.
 120          * We use 2 extra elements (=18 reals) per thread to separate thread
 121          * local data by at least a cache line. Element 0 is not used.
 122          */
 123         matrix xf_buf[GMX_OPENMP_MAX_THREADS*3];
 124
 125 #pragma omp parallel for num_threads(nthreads) schedule(static)
 126         for (int thread = 0; thread < nthreads; thread++)
 127         {
 128             int start = (nxf*thread)/nthreads;
 129             int end   = std::min(nxf*(thread + 1)/nthreads, nxf);
 130
 131             calc_x_times_f(end - start, x + start, f + start, bScrewPBC, box,
 132                            // cppcheck-suppress uninitvar
 133                            thread == 0 ? x_times_f : xf_buf[thread*3]);
 134         }
 135
 136         for (int thread = 1; thread < nthreads; thread++)
 137         {
 138             m_add(x_times_f, xf_buf[thread*3], x_times_f);
 139         }
 140     }
 141
 142     for (int d = 0; d < DIM; d++)
 143     {
 144         upd_vir(vir[d], x_times_f[d][XX], x_times_f[d][YY], x_times_f[d][ZZ]);
 145     }
 146 }
 147
 148
 149 static void lo_fcv(int i0, int i1,
 150                    real x[], real f[], tensor vir,
 151                    int is[], real box[], gmx_bool bTriclinic)
 152 {
 153     int      i, i3, tx, ty, tz;
 154     real     xx, yy, zz;
 155     real     dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
 156
 157     if (bTriclinic)
 158     {
 159         for (i = i0; (i < i1); i++)
 160         {
 161             i3 = DIM*i;
 162             tx = is[i3+XX];
 163             ty = is[i3+YY];
 164             tz = is[i3+ZZ];
 165
 166             xx    = x[i3+XX]-tx*box[XXXX]-ty*box[YYXX]-tz*box[ZZXX];
 167             dvxx += xx*f[i3+XX];
 168             dvxy += xx*f[i3+YY];
 169             dvxz += xx*f[i3+ZZ];
 170
 171             yy    = x[i3+YY]-ty*box[YYYY]-tz*box[ZZYY];
 172             dvyx += yy*f[i3+XX];
 173             dvyy += yy*f[i3+YY];
 174             dvyz += yy*f[i3+ZZ];
 175
 176             zz    = x[i3+ZZ]-tz*box[ZZZZ];
 177             dvzx += zz*f[i3+XX];
 178             dvzy += zz*f[i3+YY];
 179             dvzz += zz*f[i3+ZZ];
 180         }
 181     }
 182     else
 183     {
 184         for (i = i0; (i < i1); i++)
 185         {
 186             i3 = DIM*i;
 187             tx = is[i3+XX];
 188             ty = is[i3+YY];
 189             tz = is[i3+ZZ];
 190
 191             xx    = x[i3+XX]-tx*box[XXXX];
 192             dvxx += xx*f[i3+XX];
 193             dvxy += xx*f[i3+YY];
 194             dvxz += xx*f[i3+ZZ];
 195
 196             yy    = x[i3+YY]-ty*box[YYYY];
 197             dvyx += yy*f[i3+XX];
 198             dvyy += yy*f[i3+YY];
 199             dvyz += yy*f[i3+ZZ];
 200
 201             zz    = x[i3+ZZ]-tz*box[ZZZZ];
 202             dvzx += zz*f[i3+XX];
 203             dvzy += zz*f[i3+YY];
 204             dvzz += zz*f[i3+ZZ];
 205         }
 206     }
 207
 208     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 209     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 210     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 211 }
 212
 213 void f_calc_vir(int i0, int i1, rvec x[], rvec f[], tensor vir,
 214                 t_graph *g, matrix box)
 215 {
 216     int start, end;
 217
 218     if (g && (g->nnodes > 0))
 219     {
 220         /* Calculate virial for bonded forces only when they belong to
 221          * this node.
 222          */
 223         start = std::max(i0, g->at_start);
 224         end   = std::min(i1, g->at_end);
 225         lo_fcv(start, end, x[0], f[0], vir, g->ishift[0], box[0], TRICLINIC(box));
 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         {
 233             calc_vir(start-i0, x + i0, f + i0, vir, FALSE, box);
 234         }
 235         if (end < i1)
 236         {
 237             calc_vir(i1-end, x + end, f + end, vir, FALSE, box);
 238         }
 239     }
 240     else
 241     {
 242         calc_vir(i1-i0, x + i0, f + i0, vir, FALSE, box);
 243     }
 244 }