src/gromacs/domdec/domdec_box.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2014,2015,2017,2018, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  *
  38  * \brief This file defines functions used by the domdec module
  39  * for (bounding) box and pbc information generation.
  40  *
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_domdec
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "gromacs/domdec/domdec.h"
  48 #include "gromacs/domdec/domdec_network.h"
  49 #include "gromacs/domdec/domdec_struct.h"
  50 #include "gromacs/gmxlib/network.h"
  51 #include "gromacs/math/functions.h"
  52 #include "gromacs/math/vec.h"
  53 #include "gromacs/mdlib/nsgrid.h"
  54 #include "gromacs/mdtypes/commrec.h"
  55 #include "gromacs/mdtypes/inputrec.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/utility/fatalerror.h"
  58
  59 /*! \brief Calculates the average and standard deviation in 3D of n charge groups */
  60 static void calc_cgcm_av_stddev(const t_block *cgs, int n, const rvec *x,
  61                                 rvec av, rvec stddev,
  62                                 t_commrec *cr_sum)
  63 {
  64     int   *cgindex;
  65     dvec   s1, s2;
  66     double buf[7];
  67     int    cg, d, k0, k1, k, nrcg;
  68     real   inv_ncg;
  69     rvec   cg_cm;
  70
  71     clear_dvec(s1);
  72     clear_dvec(s2);
  73
  74     cgindex = cgs->index;
  75     for (cg = 0; cg < n; cg++)
  76     {
  77         k0      = cgindex[cg];
  78         k1      = cgindex[cg+1];
  79         nrcg    = k1 - k0;
  80         if (nrcg == 1)
  81         {
  82             copy_rvec(x[k0], cg_cm);
  83         }
  84         else
  85         {
  86             inv_ncg = 1.0/nrcg;
  87
  88             clear_rvec(cg_cm);
  89             for (k = k0; (k < k1); k++)
  90             {
  91                 rvec_inc(cg_cm, x[k]);
  92             }
  93             for (d = 0; (d < DIM); d++)
  94             {
  95                 cg_cm[d] *= inv_ncg;
  96             }
  97         }
  98         for (d = 0; d < DIM; d++)
  99         {
 100             s1[d] += cg_cm[d];
 101             s2[d] += cg_cm[d]*cg_cm[d];
 102         }
 103     }
 104
 105     if (cr_sum != nullptr)
 106     {
 107         for (d = 0; d < DIM; d++)
 108         {
 109             buf[d]     = s1[d];
 110             buf[DIM+d] = s2[d];
 111         }
 112         buf[6] = n;
 113         gmx_sumd(7, buf, cr_sum);
 114         for (d = 0; d < DIM; d++)
 115         {
 116             s1[d] = buf[d];
 117             s2[d] = buf[DIM+d];
 118         }
 119         n = (int)(buf[6] + 0.5);
 120     }
 121
 122     dsvmul(1.0/n, s1, s1);
 123     dsvmul(1.0/n, s2, s2);
 124
 125     for (d = 0; d < DIM; d++)
 126     {
 127         av[d]     = s1[d];
 128         stddev[d] = std::sqrt(s2[d] - s1[d]*s1[d]);
 129     }
 130 }
 131
 132 /*! \brief Determines if dimensions require triclinic treatment and stores this info in ddbox */
 133 static void set_tric_dir(const ivec *dd_nc, gmx_ddbox_t *ddbox, const matrix box)
 134 {
 135     int   npbcdim, d, i, j;
 136     rvec *v, *normal;
 137     real  dep, inv_skew_fac2;
 138
 139     npbcdim = ddbox->npbcdim;
 140     normal  = ddbox->normal;
 141     for (d = 0; d < DIM; d++)
 142     {
 143         ddbox->tric_dir[d] = 0;
 144         for (j = d+1; j < npbcdim; j++)
 145         {
 146             if (box[j][d] != 0)
 147             {
 148                 ddbox->tric_dir[d] = 1;
 149                 if (dd_nc != nullptr && (*dd_nc)[j] > 1 && (*dd_nc)[d] == 1)
 150                 {
 151                     gmx_fatal(FARGS, "Domain decomposition has not been implemented for box vectors that have non-zero components in directions that do not use domain decomposition: ncells = %d %d %d, box vector[%d] = %f %f %f",
 152                               (*dd_nc)[XX], (*dd_nc)[YY], (*dd_nc)[ZZ],
 153                               j+1, box[j][XX], box[j][YY], box[j][ZZ]);
 154                 }
 155             }
 156         }
 157
 158         /* Construct vectors v for dimension d that are perpendicular
 159          * to the triclinic plane of dimension d. Each vector v[i] has
 160          * v[i][i]=1 and v[i][d]!=0 for triclinic dimensions, while the third
 161          * component is zero. These are used for computing the distance
 162          * to a triclinic plane given the distance along dimension d.
 163          * Set the trilinic skewing factor that translates
 164          * the thickness of a slab perpendicular to this dimension
 165          * into the real thickness of the slab.
 166          */
 167         if (ddbox->tric_dir[d])
 168         {
 169             inv_skew_fac2 = 1;
 170             v             = ddbox->v[d];
 171             if (d == XX || d == YY)
 172             {
 173                 /* Normalize such that the "diagonal" is 1 */
 174                 svmul(1/box[d+1][d+1], box[d+1], v[d+1]);
 175                 for (i = 0; i < d; i++)
 176                 {
 177                     v[d+1][i] = 0;
 178                 }
 179                 inv_skew_fac2 += gmx::square(v[d+1][d]);
 180                 if (d == XX)
 181                 {
 182                     /* Normalize such that the "diagonal" is 1 */
 183                     svmul(1/box[d+2][d+2], box[d+2], v[d+2]);
 184                     /* Set v[d+2][d+1] to zero by shifting along v[d+1] */
 185                     dep = v[d+2][d+1]/v[d+1][d+1];
 186                     for (i = 0; i < DIM; i++)
 187                     {
 188                         v[d+2][i] -= dep*v[d+1][i];
 189                     }
 190                     inv_skew_fac2 += gmx::square(v[d+2][d]);
 191
 192                     cprod(v[d+1], v[d+2], normal[d]);
 193                 }
 194                 else
 195                 {
 196                     /* cross product with (1,0,0) */
 197                     normal[d][XX] =  0;
 198                     normal[d][YY] =  v[d+1][ZZ];
 199                     normal[d][ZZ] = -v[d+1][YY];
 200                 }
 201                 if (debug)
 202                 {
 203                     fprintf(debug, "box[%d]  %.3f %.3f %.3f\n",
 204                             d, box[d][XX], box[d][YY], box[d][ZZ]);
 205                     for (i = d+1; i < DIM; i++)
 206                     {
 207                         fprintf(debug, "  v[%d]  %.3f %.3f %.3f\n",
 208                                 i, v[i][XX], v[i][YY], v[i][ZZ]);
 209                     }
 210                 }
 211             }
 212             ddbox->skew_fac[d] = 1.0/std::sqrt(inv_skew_fac2);
 213             /* Set the normal vector length to skew_fac */
 214             dep = ddbox->skew_fac[d]/norm(normal[d]);
 215             svmul(dep, normal[d], normal[d]);
 216
 217             if (debug)
 218             {
 219                 fprintf(debug, "skew_fac[%d] = %f\n", d, ddbox->skew_fac[d]);
 220                 fprintf(debug, "normal[%d]  %.3f %.3f %.3f\n",
 221                         d, normal[d][XX], normal[d][YY], normal[d][ZZ]);
 222             }
 223         }
 224         else
 225         {
 226             ddbox->skew_fac[d] = 1;
 227
 228             for (i = 0; i < DIM; i++)
 229             {
 230                 clear_rvec(ddbox->v[d][i]);
 231                 ddbox->v[d][i][i] = 1;
 232             }
 233             clear_rvec(normal[d]);
 234             normal[d][d] = 1;
 235         }
 236     }
 237 }
 238
 239 /*! \brief This function calculates bounding box and pbc info and populates ddbox */
 240 static void low_set_ddbox(const t_inputrec *ir, const ivec *dd_nc, const matrix box,
 241                           gmx_bool bCalcUnboundedSize, int ncg, const t_block *cgs, const rvec *x,
 242                           t_commrec *cr_sum,
 243                           gmx_ddbox_t *ddbox)
 244 {
 245     rvec av, stddev;
 246     real b0, b1;
 247     int  d;
 248
 249     ddbox->npbcdim     = ePBC2npbcdim(ir->ePBC);
 250     ddbox->nboundeddim = inputrec2nboundeddim(ir);
 251
 252     for (d = 0; d < ddbox->nboundeddim; d++)
 253     {
 254         ddbox->box0[d]     = 0;
 255         ddbox->box_size[d] = box[d][d];
 256     }
 257
 258     if (ddbox->nboundeddim < DIM && bCalcUnboundedSize)
 259     {
 260         calc_cgcm_av_stddev(cgs, ncg, x, av, stddev, cr_sum);
 261
 262         /* GRID_STDDEV_FAC * stddev
 263          * gives a uniform load for a rectangular block of cg's.
 264          * For a sphere it is not a bad approximation for 4x1x1 up to 4x2x2.
 265          */
 266         for (d = ddbox->nboundeddim; d < DIM; d++)
 267         {
 268             b0 = av[d] - GRID_STDDEV_FAC*stddev[d];
 269             b1 = av[d] + GRID_STDDEV_FAC*stddev[d];
 270             if (debug)
 271             {
 272                 fprintf(debug, "Setting global DD grid boundaries to %f - %f\n",
 273                         b0, b1);
 274             }
 275             ddbox->box0[d]     = b0;
 276             ddbox->box_size[d] = b1 - b0;
 277         }
 278     }
 279
 280     set_tric_dir(dd_nc, ddbox, box);
 281 }
 282
 283 void set_ddbox(gmx_domdec_t *dd, gmx_bool bMasterState, t_commrec *cr_sum,
 284                const t_inputrec *ir, const matrix box,
 285                gmx_bool bCalcUnboundedSize, const t_block *cgs, const rvec *x,
 286                gmx_ddbox_t *ddbox)
 287 {
 288     if (!bMasterState || DDMASTER(dd))
 289     {
 290         low_set_ddbox(ir, &dd->nc, box, bCalcUnboundedSize,
 291                       bMasterState ? cgs->nr : dd->ncg_home, cgs, x,
 292                       bMasterState ? nullptr : cr_sum,
 293                       ddbox);
 294     }
 295
 296     if (bMasterState)
 297     {
 298         dd_bcast(dd, sizeof(gmx_ddbox_t), ddbox);
 299     }
 300 }
 301
 302 void set_ddbox_cr(t_commrec *cr, const ivec *dd_nc,
 303                   const t_inputrec *ir, const matrix box,
 304                   const t_block *cgs, const rvec *x,
 305                   gmx_ddbox_t *ddbox)
 306 {
 307     if (MASTER(cr))
 308     {
 309         low_set_ddbox(ir, dd_nc, box, TRUE, cgs->nr, cgs, x, nullptr, ddbox);
 310     }
 311
 312     gmx_bcast(sizeof(gmx_ddbox_t), ddbox, cr);
 313 }