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,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #include "gmxpre.h"
  40
  41 #include "calcvir.h"
  42
  43 #include "config.h" /* for GMX_MAX_OPENMP_THREADS */
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/math/vectypes.h"
  49 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  50 #include "gromacs/pbcutil/ishift.h"
  51 #include "gromacs/pbcutil/mshift.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/utility/gmxassert.h"
  54
  55 #define XXXX 0
  56 #define XXYY 1
  57 #define XXZZ 2
  58 #define YYXX 3
  59 #define YYYY 4
  60 #define YYZZ 5
  61 #define ZZXX 6
  62 #define ZZYY 7
  63 #define ZZZZ 8
  64
  65 static void upd_vir(rvec vir, real dvx, real dvy, real dvz)
  66 {
  67     vir[XX] -= 0.5 * dvx;
  68     vir[YY] -= 0.5 * dvy;
  69     vir[ZZ] -= 0.5 * dvz;
  70 }
  71
  72 static void calc_x_times_f(int nxf, const rvec x[], const rvec f[], gmx_bool bScrewPBC, const matrix box, 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, const rvec x[], const rvec f[], tensor vir, bool bScrewPBC, const matrix box)
 105 {
 106     matrix x_times_f;
 107
 108     int nthreads = gmx_omp_nthreads_get_simple_rvec_task(emntDefault, nxf * 9);
 109
 110     GMX_ASSERT(nthreads >= 1, "Avoids uninitialized x_times_f (warning)");
 111
 112     if (nthreads == 1)
 113     {
 114         calc_x_times_f(nxf, x, f, bScrewPBC, box, x_times_f);
 115     }
 116     else
 117     {
 118         /* Use a buffer on the stack for storing thread-local results.
 119          * We use 2 extra elements (=18 reals) per thread to separate thread
 120          * local data by at least a cache line. Element 0 is not used.
 121          */
 122         matrix xf_buf[GMX_OPENMP_MAX_THREADS * 3];
 123
 124 #pragma omp parallel for num_threads(nthreads) schedule(static)
 125         for (int thread = 0; thread < nthreads; thread++)
 126         {
 127             int start = (nxf * thread) / nthreads;
 128             int end   = std::min(nxf * (thread + 1) / nthreads, nxf);
 129
 130             calc_x_times_f(end - start, x + start, f + start, bScrewPBC, box,
 131                            thread == 0 ? x_times_f : xf_buf[thread * 3]);
 132         }
 133
 134         for (int thread = 1; thread < nthreads; thread++)
 135         {
 136             m_add(x_times_f, xf_buf[thread * 3], x_times_f);
 137         }
 138     }
 139
 140     for (int d = 0; d < DIM; d++)
 141     {
 142         upd_vir(vir[d], x_times_f[d][XX], x_times_f[d][YY], x_times_f[d][ZZ]);
 143     }
 144 }
 145
 146
 147 static void
 148 lo_fcv(int i0, int i1, const real x[], const real f[], tensor vir, const int is[], const real box[], gmx_bool bTriclinic)
 149 {
 150     int  i, i3, tx, ty, tz;
 151     real xx, yy, zz;
 152     real dvxx = 0, dvxy = 0, dvxz = 0, dvyx = 0, dvyy = 0, dvyz = 0, dvzx = 0, dvzy = 0, dvzz = 0;
 153
 154     if (bTriclinic)
 155     {
 156         for (i = i0; (i < i1); i++)
 157         {
 158             i3 = DIM * i;
 159             tx = is[i3 + XX];
 160             ty = is[i3 + YY];
 161             tz = is[i3 + ZZ];
 162
 163             xx = x[i3 + XX] - tx * box[XXXX] - ty * box[YYXX] - tz * box[ZZXX];
 164             dvxx += xx * f[i3 + XX];
 165             dvxy += xx * f[i3 + YY];
 166             dvxz += xx * f[i3 + ZZ];
 167
 168             yy = x[i3 + YY] - ty * box[YYYY] - tz * box[ZZYY];
 169             dvyx += yy * f[i3 + XX];
 170             dvyy += yy * f[i3 + YY];
 171             dvyz += yy * f[i3 + ZZ];
 172
 173             zz = x[i3 + ZZ] - tz * box[ZZZZ];
 174             dvzx += zz * f[i3 + XX];
 175             dvzy += zz * f[i3 + YY];
 176             dvzz += zz * f[i3 + ZZ];
 177         }
 178     }
 179     else
 180     {
 181         for (i = i0; (i < i1); i++)
 182         {
 183             i3 = DIM * i;
 184             tx = is[i3 + XX];
 185             ty = is[i3 + YY];
 186             tz = is[i3 + ZZ];
 187
 188             xx = x[i3 + XX] - tx * box[XXXX];
 189             dvxx += xx * f[i3 + XX];
 190             dvxy += xx * f[i3 + YY];
 191             dvxz += xx * f[i3 + ZZ];
 192
 193             yy = x[i3 + YY] - ty * box[YYYY];
 194             dvyx += yy * f[i3 + XX];
 195             dvyy += yy * f[i3 + YY];
 196             dvyz += yy * f[i3 + ZZ];
 197
 198             zz = x[i3 + ZZ] - tz * box[ZZZZ];
 199             dvzx += zz * f[i3 + XX];
 200             dvzy += zz * f[i3 + YY];
 201             dvzz += zz * f[i3 + ZZ];
 202         }
 203     }
 204
 205     upd_vir(vir[XX], dvxx, dvxy, dvxz);
 206     upd_vir(vir[YY], dvyx, dvyy, dvyz);
 207     upd_vir(vir[ZZ], dvzx, dvzy, dvzz);
 208 }
 209
 210 void f_calc_vir(int i0, int i1, const rvec x[], const rvec f[], tensor vir, const t_graph* g, const matrix box)
 211 {
 212     int start, end;
 213
 214     if (g && (g->nnodes > 0))
 215     {
 216         /* Calculate virial for bonded forces only when they belong to
 217          * this node.
 218          */
 219         start = std::max(i0, g->at_start);
 220         end   = std::min(i1, g->at_end);
 221         lo_fcv(start, end, x[0], f[0], vir, g->ishift[0], box[0], TRICLINIC(box));
 222
 223         /* If not all atoms are bonded, calculate their virial contribution
 224          * anyway, without shifting back their coordinates.
 225          * Note the nifty pointer arithmetic...
 226          */
 227         if (start > i0)
 228         {
 229             calc_vir(start - i0, x + i0, f + i0, vir, FALSE, box);
 230         }
 231         if (end < i1)
 232         {
 233             calc_vir(i1 - end, x + end, f + end, vir, FALSE, box);
 234         }
 235     }
 236     else
 237     {
 238         calc_vir(i1 - i0, x + i0, f + i0, vir, FALSE, box);
 239     }
 240 }