src/gromacs/domdec/computemultibodycutoffs.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, 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  * while managing the construction, use and error checking for
  40  * topologies local to a DD rank.
  41  *
  42  * \author Berk Hess <hess@kth.se>
  43  * \ingroup module_domdec
  44  */
  45
  46 #include "gmxpre.h"
  47
  48 #include "gromacs/domdec/computemultibodycutoffs.h"
  49
  50 #include <vector>
  51
  52 #include "gromacs/domdec/options.h"
  53 #include "gromacs/domdec/reversetopology.h"
  54 #include "gromacs/math/vec.h"
  55 #include "gromacs/mdtypes/inputrec.h"
  56 #include "gromacs/pbcutil/mshift.h"
  57 #include "gromacs/pbcutil/pbc.h"
  58 #include "gromacs/utility/arrayref.h"
  59 #include "gromacs/utility/logger.h"
  60 #include "gromacs/topology/mtop_util.h"
  61
  62 using gmx::ArrayRef;
  63 using gmx::DDBondedChecking;
  64 using gmx::RVec;
  65
  66 typedef struct
  67 {
  68     real r2;
  69     int  ftype;
  70     int  a1;
  71     int  a2;
  72 } bonded_distance_t;
  73
  74 /*! \brief Compare distance^2 \p r2 against the distance in \p bd and if larger store it along with \p ftype and atom indices \p a1 and \p a2 */
  75 static void update_max_bonded_distance(real r2, int ftype, int a1, int a2, bonded_distance_t* bd)
  76 {
  77     if (r2 > bd->r2)
  78     {
  79         bd->r2    = r2;
  80         bd->ftype = ftype;
  81         bd->a1    = a1;
  82         bd->a2    = a2;
  83     }
  84 }
  85
  86 /*! \brief Set the distance, function type and atom indices for the longest distance between atoms of molecule type \p molt for two-body and multi-body bonded interactions */
  87 static void bonded_cg_distance_mol(const gmx_moltype_t*   molt,
  88                                    const DDBondedChecking ddBondedChecking,
  89                                    gmx_bool               bExcl,
  90                                    ArrayRef<const RVec>   x,
  91                                    bonded_distance_t*     bd_2b,
  92                                    bonded_distance_t*     bd_mb)
  93 {
  94     const ReverseTopOptions rtOptions(ddBondedChecking);
  95
  96     for (int ftype = 0; ftype < F_NRE; ftype++)
  97     {
  98         if (dd_check_ftype(ftype, rtOptions))
  99         {
 100             const auto& il   = molt->ilist[ftype];
 101             int         nral = NRAL(ftype);
 102             if (nral > 1)
 103             {
 104                 for (int i = 0; i < il.size(); i += 1 + nral)
 105                 {
 106                     for (int ai = 0; ai < nral; ai++)
 107                     {
 108                         int atomI = il.iatoms[i + 1 + ai];
 109                         for (int aj = ai + 1; aj < nral; aj++)
 110                         {
 111                             int atomJ = il.iatoms[i + 1 + aj];
 112                             if (atomI != atomJ)
 113                             {
 114                                 real rij2 = distance2(x[atomI], x[atomJ]);
 115
 116                                 update_max_bonded_distance(
 117                                         rij2, ftype, atomI, atomJ, (nral == 2) ? bd_2b : bd_mb);
 118                             }
 119                         }
 120                     }
 121                 }
 122             }
 123         }
 124     }
 125     if (bExcl)
 126     {
 127         const auto& excls = molt->excls;
 128         for (gmx::index ai = 0; ai < excls.ssize(); ai++)
 129         {
 130             for (const int aj : excls[ai])
 131             {
 132                 if (ai != aj)
 133                 {
 134                     real rij2 = distance2(x[ai], x[aj]);
 135
 136                     /* There is no function type for exclusions, use -1 */
 137                     update_max_bonded_distance(rij2, -1, ai, aj, bd_2b);
 138                 }
 139             }
 140         }
 141     }
 142 }
 143
 144 /*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */
 145 static void bonded_distance_intermol(const InteractionLists& ilists_intermol,
 146                                      const DDBondedChecking  ddBondedChecking,
 147                                      ArrayRef<const RVec>    x,
 148                                      PbcType                 pbcType,
 149                                      const matrix            box,
 150                                      bonded_distance_t*      bd_2b,
 151                                      bonded_distance_t*      bd_mb)
 152 {
 153     t_pbc pbc;
 154
 155     set_pbc(&pbc, pbcType, box);
 156
 157     const ReverseTopOptions rtOptions(ddBondedChecking);
 158
 159     for (int ftype = 0; ftype < F_NRE; ftype++)
 160     {
 161         if (dd_check_ftype(ftype, rtOptions))
 162         {
 163             const auto& il   = ilists_intermol[ftype];
 164             int         nral = NRAL(ftype);
 165
 166             /* No nral>1 check here, since intermol interactions always
 167              * have nral>=2 (and the code is also correct for nral=1).
 168              */
 169             for (int i = 0; i < il.size(); i += 1 + nral)
 170             {
 171                 for (int ai = 0; ai < nral; ai++)
 172                 {
 173                     int atom_i = il.iatoms[i + 1 + ai];
 174
 175                     for (int aj = ai + 1; aj < nral; aj++)
 176                     {
 177                         rvec dx;
 178
 179                         int atom_j = il.iatoms[i + 1 + aj];
 180
 181                         pbc_dx(&pbc, x[atom_i], x[atom_j], dx);
 182
 183                         const real rij2 = norm2(dx);
 184
 185                         update_max_bonded_distance(rij2, ftype, atom_i, atom_j, (nral == 2) ? bd_2b : bd_mb);
 186                     }
 187                 }
 188             }
 189         }
 190     }
 191 }
 192
 193 //! Returns whether \p molt has at least one virtual site
 194 static bool moltypeHasVsite(const gmx_moltype_t& molt)
 195 {
 196     bool hasVsite = false;
 197     for (int i = 0; i < F_NRE; i++)
 198     {
 199         if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty())
 200         {
 201             hasVsite = true;
 202         }
 203     }
 204
 205     return hasVsite;
 206 }
 207
 208 //! Returns coordinates not broken over PBC for a molecule
 209 static void getWholeMoleculeCoordinates(const gmx_moltype_t*  molt,
 210                                         const gmx_ffparams_t* ffparams,
 211                                         PbcType               pbcType,
 212                                         t_graph*              graph,
 213                                         const matrix          box,
 214                                         ArrayRef<const RVec>  x,
 215                                         ArrayRef<RVec>        xs)
 216 {
 217     if (pbcType != PbcType::No)
 218     {
 219         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(x.data()));
 220
 221         shift_x(graph, box, as_rvec_array(x.data()), as_rvec_array(xs.data()));
 222         /* By doing an extra mk_mshift the molecules that are broken
 223          * because they were e.g. imported from another software
 224          * will be made whole again. Such are the healing powers
 225          * of GROMACS.
 226          */
 227         mk_mshift(nullptr, graph, pbcType, box, as_rvec_array(xs.data()));
 228     }
 229     else
 230     {
 231         /* We copy the coordinates so the original coordinates remain
 232          * unchanged, just to be 100% sure that we do not affect
 233          * binary reproducibility of simulations.
 234          */
 235         for (int i = 0; i < molt->atoms.nr; i++)
 236         {
 237             copy_rvec(x[i], xs[i]);
 238         }
 239     }
 240
 241     if (moltypeHasVsite(*molt))
 242     {
 243         gmx::constructVirtualSites(xs, ffparams->iparams, molt->ilist);
 244     }
 245 }
 246
 247 void dd_bonded_cg_distance(const gmx::MDLogger&   mdlog,
 248                            const gmx_mtop_t&      mtop,
 249                            const t_inputrec&      inputrec,
 250                            ArrayRef<const RVec>   x,
 251                            const matrix           box,
 252                            const DDBondedChecking ddBondedChecking,
 253                            real*                  r_2b,
 254                            real*                  r_mb)
 255 {
 256     bonded_distance_t bd_2b = { 0, -1, -1, -1 };
 257     bonded_distance_t bd_mb = { 0, -1, -1, -1 };
 258
 259     bool bExclRequired = inputrecExclForces(&inputrec);
 260
 261     *r_2b         = 0;
 262     *r_mb         = 0;
 263     int at_offset = 0;
 264     for (const gmx_molblock_t& molb : mtop.molblock)
 265     {
 266         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 267         if (molt.atoms.nr == 1 || molb.nmol == 0)
 268         {
 269             at_offset += molb.nmol * molt.atoms.nr;
 270         }
 271         else
 272         {
 273             t_graph graph;
 274             if (inputrec.pbcType != PbcType::No)
 275             {
 276                 graph = mk_graph_moltype(molt);
 277             }
 278
 279             std::vector<RVec> xs(molt.atoms.nr);
 280             for (int mol = 0; mol < molb.nmol; mol++)
 281             {
 282                 getWholeMoleculeCoordinates(&molt,
 283                                             &mtop.ffparams,
 284                                             inputrec.pbcType,
 285                                             &graph,
 286                                             box,
 287                                             x.subArray(at_offset, molt.atoms.nr),
 288                                             xs);
 289
 290                 bonded_distance_t bd_mol_2b = { 0, -1, -1, -1 };
 291                 bonded_distance_t bd_mol_mb = { 0, -1, -1, -1 };
 292
 293                 bonded_cg_distance_mol(&molt, ddBondedChecking, bExclRequired, xs, &bd_mol_2b, &bd_mol_mb);
 294
 295                 /* Process the mol data adding the atom index offset */
 296                 update_max_bonded_distance(bd_mol_2b.r2,
 297                                            bd_mol_2b.ftype,
 298                                            at_offset + bd_mol_2b.a1,
 299                                            at_offset + bd_mol_2b.a2,
 300                                            &bd_2b);
 301                 update_max_bonded_distance(bd_mol_mb.r2,
 302                                            bd_mol_mb.ftype,
 303                                            at_offset + bd_mol_mb.a1,
 304                                            at_offset + bd_mol_mb.a2,
 305                                            &bd_mb);
 306
 307                 at_offset += molt.atoms.nr;
 308             }
 309         }
 310     }
 311
 312     if (mtop.bIntermolecularInteractions)
 313     {
 314         GMX_RELEASE_ASSERT(mtop.intermolecular_ilist,
 315                            "We should have an ilist when intermolecular interactions are on");
 316
 317         bonded_distance_intermol(
 318                 *mtop.intermolecular_ilist, ddBondedChecking, x, inputrec.pbcType, box, &bd_2b, &bd_mb);
 319     }
 320
 321     *r_2b = sqrt(bd_2b.r2);
 322     *r_mb = sqrt(bd_mb.r2);
 323
 324     if (*r_2b > 0 || *r_mb > 0)
 325     {
 326         GMX_LOG(mdlog.info).appendText("Initial maximum distances in bonded interactions:");
 327         if (*r_2b > 0)
 328         {
 329             GMX_LOG(mdlog.info)
 330                     .appendTextFormatted(
 331                             "    two-body bonded interactions: %5.3f nm, %s, atoms %d %d",
 332                             *r_2b,
 333                             (bd_2b.ftype >= 0) ? interaction_function[bd_2b.ftype].longname : "Exclusion",
 334                             bd_2b.a1 + 1,
 335                             bd_2b.a2 + 1);
 336         }
 337         if (*r_mb > 0)
 338         {
 339             GMX_LOG(mdlog.info)
 340                     .appendTextFormatted(
 341                             "  multi-body bonded interactions: %5.3f nm, %s, atoms %d %d",
 342                             *r_mb,
 343                             interaction_function[bd_mb.ftype].longname,
 344                             bd_mb.a1 + 1,
 345                             bd_mb.a2 + 1);
 346         }
 347     }
 348 }