src/gromacs/pbcutil/com.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020, 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 /*!\file
  36  * \internal
  37  * \brief
  38  * Implements helper methods to place particle COM in boxes.
  39  *
  40  * \author Paul Bauer <paul.bauer.q@gmail.com>
  41  * \ingroup module_pbcutil
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "com.h"
  46
  47 #include <algorithm>
  48 #include <vector>
  49
  50 #include "gromacs/pbcutil/pbc.h"
  51 #include "gromacs/topology/mtop_util.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/range.h"
  54
  55 namespace gmx
  56 {
  57
  58 namespace
  59 {
  60
  61 /*! \brief
  62  * Calculates shift to place COM into a box.
  63  *
  64  * \param[in] position     Unshifted COM position.
  65  * \param[in] box          The current box to place the COM in.
  66  * \param[in] pbcType      What kind of box is being used.
  67  * \param[in] unitCellType Type of unitcell being used.
  68  * \param[in] centerType   How things should be centered.
  69  * \returns The shift needed to place the COM into the box.
  70  */
  71 RVec evaluateShiftToBox(const RVec&          position,
  72                         const matrix         box,
  73                         const PbcType&       pbcType,
  74                         const UnitCellType&  unitCellType,
  75                         const CenteringType& centerType)
  76 {
  77     RVec newCOMPosition(position);
  78     auto comArrayRef = arrayRefFromArray(&newCOMPosition, 1);
  79     switch (unitCellType)
  80     {
  81         case (UnitCellType::Rectangular): put_atoms_in_box(pbcType, box, comArrayRef); break;
  82         case (UnitCellType::Triclinic):
  83             put_atoms_in_triclinic_unitcell(static_cast<int>(centerType), box, comArrayRef);
  84             break;
  85         case (UnitCellType::Compact):
  86             put_atoms_in_compact_unitcell(pbcType, static_cast<int>(centerType), box, comArrayRef);
  87             break;
  88         default: GMX_RELEASE_ASSERT(false, "Unhandled type of unit cell");
  89     }
  90     RVec shift(0, 0, 0);
  91     rvec_sub(newCOMPosition, position, shift);
  92     return shift;
  93 }
  94
  95 /*! \brief
  96  * Calculates the COM for each collection of atoms.
  97  *
  98  * \param[in] x       View on coordinates of the molecule.
  99  * \param[in] moltype Which molecule type to calculate for.
 100  * \param[in] atomOffset If needed, point from where to count the first atom to process.
 101  * \returns The center of mass for the molecule.
 102  */
 103 RVec calculateCOM(ArrayRef<const RVec> x, const gmx_moltype_t& moltype, const int atomOffset = 0)
 104 {
 105     RVec   com(0, 0, 0);
 106     double totalMass   = 0;
 107     int    currentAtom = atomOffset;
 108     for (const auto& coord : x)
 109     {
 110         real mass = moltype.atoms.atom[currentAtom].m;
 111         for (int d = 0; d < DIM; d++)
 112         {
 113             com[d] += mass * coord[d];
 114         }
 115         totalMass += mass;
 116         currentAtom++;
 117     }
 118     svmul(1.0 / totalMass, com, com);
 119     return com;
 120 }
 121
 122 } // namespace
 123
 124 void shiftAtoms(const RVec& shift, ArrayRef<RVec> x)
 125 {
 126     std::transform(std::begin(x), std::end(x), std::begin(x), [shift](RVec x) { return x + shift; });
 127 }
 128
 129 void placeCoordinatesWithCOMInBox(const PbcType&      pbcType,
 130                                   const UnitCellType  unitCellType,
 131                                   const CenteringType centerType,
 132                                   const matrix        box,
 133                                   ArrayRef<RVec>      x,
 134                                   const gmx_mtop_t&   mtop,
 135                                   const COMShiftType  comShiftType)
 136 {
 137     GMX_RELEASE_ASSERT(comShiftType != COMShiftType::Count, "Using COUNT of enumeration");
 138     // loop over all molecule blocks, then over all molecules in these blocks
 139     int molb = 0;
 140     for (const auto& molblock : mtop.molblock)
 141     {
 142         const MoleculeBlockIndices& ind              = mtop.moleculeBlockIndices[molb];
 143         const gmx_moltype_t&        moltype          = mtop.moltype[molblock.type];
 144         const int                   atomsPerMolecule = ind.numAtomsPerMolecule;
 145         int                         atomStart        = ind.globalAtomStart;
 146         for (int mol = 0; mol < molblock.nmol; mol++)
 147         {
 148             std::vector<Range<int>> atomRanges;
 149             if (comShiftType == COMShiftType::Molecule)
 150             {
 151                 atomRanges.emplace_back(gmx::Range<int>(0, atomsPerMolecule));
 152             }
 153             else if (comShiftType == COMShiftType::Residue)
 154             {
 155                 atomRanges = atomRangeOfEachResidue(moltype);
 156             }
 157             for (const auto& atomRange : atomRanges)
 158             {
 159                 const int      firstAtomGlobalPosition = atomStart + *atomRange.begin();
 160                 ArrayRef<RVec> coords = x.subArray(firstAtomGlobalPosition, atomRange.size());
 161                 const RVec     com    = calculateCOM(coords, moltype, *atomRange.begin());
 162                 const RVec shift = evaluateShiftToBox(com, box, pbcType, unitCellType, centerType);
 163                 shiftAtoms(shift, coords);
 164             }
 165             atomStart += atomsPerMolecule;
 166         }
 167         molb++;
 168     }
 169 }
 170
 171 } // namespace gmx