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,2021, 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/pbc.h"
  52 #include "gromacs/utility/gmxassert.h"
  53
  54 static void upd_vir(rvec vir, real dvx, real dvy, real dvz)
  55 {
  56     vir[XX] -= 0.5 * dvx;
  57     vir[YY] -= 0.5 * dvy;
  58     vir[ZZ] -= 0.5 * dvz;
  59 }
  60
  61 static void calc_x_times_f(int nxf, const rvec x[], const rvec f[], gmx_bool bScrewPBC, const matrix box, matrix x_times_f)
  62 {
  63     clear_mat(x_times_f);
  64
  65     for (int i = 0; i < nxf; i++)
  66     {
  67         for (int d = 0; d < DIM; d++)
  68         {
  69             for (int n = 0; n < DIM; n++)
  70             {
  71                 x_times_f[d][n] += x[i][d] * f[i][n];
  72             }
  73         }
  74
  75         if (bScrewPBC)
  76         {
  77             int isx = gmx::shiftIndexToXDim(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                 for (int d = 0; d < DIM; d++)
  82                 {
  83                     for (int n = 0; n < DIM; n++)
  84                     {
  85                         x_times_f[d][n] += box[d][d] * f[i][n];
  86                     }
  87                 }
  88             }
  89         }
  90     }
  91 }
  92
  93 void calc_vir(int nxf, const rvec x[], const rvec f[], tensor vir, bool bScrewPBC, const matrix box)
  94 {
  95     matrix x_times_f;
  96
  97     int nthreads = gmx_omp_nthreads_get_simple_rvec_task(ModuleMultiThread::Default, nxf * 9);
  98
  99     GMX_ASSERT(nthreads >= 1, "Avoids uninitialized x_times_f (warning)");
 100
 101     if (nthreads == 1)
 102     {
 103         calc_x_times_f(nxf, x, f, bScrewPBC, box, x_times_f);
 104     }
 105     else
 106     {
 107         /* Use a buffer on the stack for storing thread-local results.
 108          * We use 2 extra elements (=18 reals) per thread to separate thread
 109          * local data by at least a cache line. Element 0 is not used.
 110          */
 111         matrix xf_buf[GMX_OPENMP_MAX_THREADS * 3];
 112
 113 #pragma omp parallel for num_threads(nthreads) schedule(static)
 114         for (int thread = 0; thread < nthreads; thread++)
 115         {
 116             int start = (nxf * thread) / nthreads;
 117             int end   = std::min(nxf * (thread + 1) / nthreads, nxf);
 118
 119             calc_x_times_f(end - start,
 120                            x + start,
 121                            f + start,
 122                            bScrewPBC,
 123                            box,
 124                            thread == 0 ? x_times_f : xf_buf[thread * 3]);
 125         }
 126
 127         for (int thread = 1; thread < nthreads; thread++)
 128         {
 129             m_add(x_times_f, xf_buf[thread * 3], x_times_f);
 130         }
 131     }
 132
 133     for (int d = 0; d < DIM; d++)
 134     {
 135         upd_vir(vir[d], x_times_f[d][XX], x_times_f[d][YY], x_times_f[d][ZZ]);
 136     }
 137 }
 138
 139 void f_calc_vir(int i0, int i1, const rvec x[], const rvec f[], tensor vir, const matrix box)
 140 {
 141     calc_vir(i1 - i0, x + i0, f + i0, vir, FALSE, box);
 142 }