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  */
evaluateShiftToBox(const RVec & position,const matrix box,const PbcType & pbcType,const UnitCellType & unitCellType,const CenteringType & centerType)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  */
calculateCOM(ArrayRef<const RVec> x,const gmx_moltype_t & moltype,const int atomOffset=0)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 
shiftAtoms(const RVec & shift,ArrayRef<RVec> x)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 
placeCoordinatesWithCOMInBox(const PbcType & pbcType,const UnitCellType unitCellType,const CenteringType centerType,const matrix box,ArrayRef<RVec> x,const gmx_mtop_t & mtop,const COMShiftType comShiftType)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
172