src/gromacs/domdec/nsgrid.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,2017,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 "nsgrid.h"
  42
  43 #include <cmath>
  44 #include <cstdio>
  45
  46 #include "gromacs/domdec/dlb.h"
  47 #include "gromacs/domdec/domdec.h"
  48 #include "gromacs/domdec/domdec_struct.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/utility/fatalerror.h"
  51
  52 /*! \brief The extent of the neighborsearch grid is a bit larger than sqrt(3)
  53  * to account for less dense regions at the edges of the system.
  54  */
  55 constexpr real c_stdDevFactor = 2.0;
  56
  57 /***********************************
  58  *         Grid Routines
  59  ***********************************/
  60
  61 static void calc_x_av_stddev(int n, rvec* x, rvec av, rvec stddev)
  62 {
  63     dvec s1, s2;
  64     int  i, d;
  65
  66     clear_dvec(s1);
  67     clear_dvec(s2);
  68
  69     for (i = 0; i < n; i++)
  70     {
  71         for (d = 0; d < DIM; d++)
  72         {
  73             s1[d] += x[i][d];
  74             s2[d] += x[i][d] * x[i][d];
  75         }
  76     }
  77
  78     dsvmul(1.0 / n, s1, s1);
  79     dsvmul(1.0 / n, s2, s2);
  80
  81     for (d = 0; d < DIM; d++)
  82     {
  83         av[d]     = s1[d];
  84         stddev[d] = std::sqrt(s2[d] - s1[d] * s1[d]);
  85     }
  86 }
  87
  88 static void get_nsgrid_boundaries_vac(real av, real stddev, real* bound0, real* bound1, real* bdens0, real* bdens1)
  89 {
  90     /* Set the grid to 2 times the standard deviation of
  91      * the charge group centers in both directions.
  92      * For a uniform bounded distribution the width is sqrt(3)*stddev,
  93      * so all charge groups fall within the width.
  94      * For a sphere stddev is r/sqrt(5): 99.2% falls within the width.
  95      * For a Gaussian distribution 98% fall within the width.
  96      */
  97     *bound0 = av - c_stdDevFactor * stddev;
  98     *bound1 = av + c_stdDevFactor * stddev;
  99
 100     *bdens0 = av - c_gridStdDevFactor * stddev;
 101     *bdens1 = av + c_gridStdDevFactor * stddev;
 102 }
 103
 104 static void dd_box_bounds_to_ns_bounds(real box0, real box_size, real* gr0, real* gr1)
 105 {
 106     real av, stddev;
 107
 108     /* Redetermine av and stddev from the DD box boundaries */
 109     av     = box0 + 0.5 * box_size;
 110     stddev = 0.5 * box_size / c_gridStdDevFactor;
 111
 112     *gr0 = av - c_stdDevFactor * stddev;
 113     *gr1 = av + c_stdDevFactor * stddev;
 114 }
 115
 116 void get_nsgrid_boundaries(int           nboundeddim,
 117                            matrix        box,
 118                            gmx_domdec_t* dd,
 119                            gmx_ddbox_t*  ddbox,
 120                            gmx::RVec*    gr0,
 121                            gmx::RVec*    gr1,
 122                            int           ncg,
 123                            rvec*         cgcm,
 124                            rvec          grid_x0,
 125                            rvec          grid_x1)
 126 {
 127     rvec av, stddev;
 128     real bdens0, bdens1;
 129     int  d;
 130
 131     if (nboundeddim < DIM)
 132     {
 133         calc_x_av_stddev(ncg, cgcm, av, stddev);
 134     }
 135
 136     for (d = 0; d < DIM; d++)
 137     {
 138         if (d < nboundeddim)
 139         {
 140             grid_x0[d] = (gr0 != nullptr ? (*gr0)[d] : 0);
 141             grid_x1[d] = (gr1 != nullptr ? (*gr1)[d] : box[d][d]);
 142         }
 143         else
 144         {
 145             if (ddbox == nullptr)
 146             {
 147                 get_nsgrid_boundaries_vac(av[d], stddev[d], &grid_x0[d], &grid_x1[d], &bdens0, &bdens1);
 148             }
 149             else
 150             {
 151                 /* Temporary fix which uses global ddbox boundaries
 152                  * for unbounded dimensions.
 153                  * Should be replaced by local boundaries, which makes
 154                  * the ns grid smaller and does not require global comm.
 155                  */
 156                 dd_box_bounds_to_ns_bounds(ddbox->box0[d], ddbox->box_size[d], &grid_x0[d], &grid_x1[d]);
 157                 bdens0 = grid_x0[d];
 158                 bdens1 = grid_x1[d];
 159             }
 160             /* Check for a DD cell not at a lower edge */
 161             if (dd != nullptr && gr0 != nullptr && dd->ci[d] > 0)
 162             {
 163                 grid_x0[d] = (*gr0)[d];
 164                 bdens0     = (*gr0)[d];
 165             }
 166             /* Check for a DD cell not at a higher edge */
 167             if (dd != nullptr && gr1 != nullptr && dd->ci[d] < dd->numCells[d] - 1)
 168             {
 169                 grid_x1[d] = (*gr1)[d];
 170                 bdens1     = (*gr1)[d];
 171             }
 172         }
 173
 174         if (debug)
 175         {
 176             fprintf(debug, "Set grid boundaries dim %d: %f %f\n", d, grid_x0[d], grid_x1[d]);
 177         }
 178     }
 179 }