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