2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
5 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * Implements test helper routines from toputils.h.
40 * \author Teemu Murtola <teemu.murtola@gmail.com>
41 * \ingroup module_selection
53 #include "gromacs/fileio/confio.h"
54 #include "gromacs/fileio/trxio.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/topology/atoms.h"
57 #include "gromacs/topology/mtop_util.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/trajectory/trajectoryframe.h"
60 #include "gromacs/utility/arrayref.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/smalloc.h"
65 #include "testutils/testfilemanager.h"
72 TopologyManager::TopologyManager() : frame_(nullptr) {}
74 TopologyManager::~TopologyManager()
76 if (frame_
!= nullptr)
85 for (char* atomtype
: atomtypes_
)
91 void TopologyManager::requestFrame()
93 GMX_RELEASE_ASSERT(mtop_
== nullptr, "Frame must be requested before initializing topology");
94 if (frame_
== nullptr)
100 void TopologyManager::requestVelocities()
102 GMX_RELEASE_ASSERT(frame_
!= nullptr, "Velocities requested before requesting a frame");
104 if (frame_
->natoms
> 0)
106 snew(frame_
->v
, frame_
->natoms
);
110 void TopologyManager::requestForces()
112 GMX_RELEASE_ASSERT(frame_
!= nullptr, "Forces requested before requesting a frame");
114 if (frame_
->natoms
> 0)
116 snew(frame_
->f
, frame_
->natoms
);
120 void TopologyManager::loadTopology(const char* filename
)
124 rvec
* xtop
= nullptr;
127 GMX_RELEASE_ASSERT(mtop_
== nullptr, "Topology initialized more than once");
128 mtop_
= std::make_unique
<gmx_mtop_t
>();
129 readConfAndTopology(gmx::test::TestFileManager::getInputFilePath(filename
).c_str(), &fullTopology
,
130 mtop_
.get(), &pbcType
, frame_
!= nullptr ? &xtop
: nullptr, nullptr, box
);
132 if (frame_
!= nullptr)
134 GMX_ASSERT(xtop
!= nullptr, "Keep the static analyzer happy");
135 frame_
->natoms
= mtop_
->natoms
;
137 snew(frame_
->x
, frame_
->natoms
);
138 std::memcpy(frame_
->x
, xtop
, sizeof(*frame_
->x
) * frame_
->natoms
);
140 copy_mat(box
, frame_
->box
);
146 void TopologyManager::initAtoms(int count
)
148 GMX_RELEASE_ASSERT(mtop_
== nullptr, "Topology initialized more than once");
149 mtop_
= std::make_unique
<gmx_mtop_t
>();
150 mtop_
->moltype
.resize(1);
151 init_t_atoms(&mtop_
->moltype
[0].atoms
, count
, FALSE
);
152 mtop_
->molblock
.resize(1);
153 mtop_
->molblock
[0].type
= 0;
154 mtop_
->molblock
[0].nmol
= 1;
155 mtop_
->natoms
= count
;
157 GMX_RELEASE_ASSERT(mtop_
->maxResiduesPerMoleculeToTriggerRenumber() == 0,
158 "maxres_renum in mtop can be modified by an env.var., that is not supported "
160 t_atoms
& atoms
= this->atoms();
161 for (int i
= 0; i
< count
; ++i
)
163 atoms
.atom
[i
].m
= (i
% 3 == 0 ? 2.0 : 1.0);
165 atoms
.haveMass
= TRUE
;
166 if (frame_
!= nullptr)
168 frame_
->natoms
= count
;
170 snew(frame_
->x
, count
);
173 snew(frame_
->v
, count
);
177 snew(frame_
->f
, count
);
182 void TopologyManager::initAtomTypes(const ArrayRef
<const char* const>& types
)
184 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
185 atomtypes_
.reserve(types
.size());
186 for (const char* type
: types
)
188 atomtypes_
.push_back(gmx_strdup(type
));
190 t_atoms
& atoms
= this->atoms();
191 snew(atoms
.atomtype
, atoms
.nr
);
193 for (int i
= 0; i
< atoms
.nr
; ++i
, ++j
)
195 if (j
== types
.ssize())
199 atoms
.atomtype
[i
] = &atomtypes_
[j
];
201 atoms
.haveType
= TRUE
;
204 void TopologyManager::initTopology(const int numMoleculeTypes
, const int numMoleculeBlocks
)
206 GMX_RELEASE_ASSERT(mtop_
== nullptr, "Topology initialized more than once");
207 mtop_
= std::make_unique
<gmx_mtop_t
>();
208 mtop_
->moltype
.resize(numMoleculeTypes
);
209 mtop_
->molblock
.resize(numMoleculeBlocks
);
212 void TopologyManager::setMoleculeType(const int moleculeTypeIndex
, const ArrayRef
<const int> residueSizes
)
214 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
216 // Make a molecule block that refers to a new molecule type
217 auto& newMoleculeType
= mtop_
->moltype
[moleculeTypeIndex
];
218 auto& atoms
= newMoleculeType
.atoms
;
220 auto numAtomsInMolecule
= std::accumulate(residueSizes
.begin(), residueSizes
.end(), 0);
221 init_t_atoms(&atoms
, numAtomsInMolecule
, FALSE
);
222 atoms
.nres
= residueSizes
.size();
224 // Fill the residues in the molecule type
225 int residueIndex
= 0;
226 auto residueSize
= std::begin(residueSizes
);
227 for (int i
= 0; i
< atoms
.nr
&& residueSize
!= std::end(residueSizes
); ++residueSize
, ++residueIndex
)
229 for (int j
= 0; i
< atoms
.nr
&& j
< *residueSize
; ++i
, ++j
)
231 atoms
.atom
[i
].resind
= residueIndex
;
232 atoms
.atom
[i
].m
= (i
% 3 == 0 ? 2.0 : 1.0);
235 atoms
.nres
= residueIndex
;
236 atoms
.haveMass
= true;
239 void TopologyManager::setMoleculeBlock(const int moleculeBlockIndex
,
240 const int moleculeTypeIndex
,
241 const int numMoleculesToAdd
)
243 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
245 auto& newMoleculeBlock
= mtop_
->molblock
[moleculeBlockIndex
];
246 newMoleculeBlock
.type
= moleculeTypeIndex
;
247 newMoleculeBlock
.nmol
= numMoleculesToAdd
;
249 mtop_
->natoms
+= numMoleculesToAdd
* mtop_
->moltype
[moleculeTypeIndex
].atoms
.nr
;
252 void TopologyManager::finalizeTopology()
254 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
256 mtop_
->haveMoleculeIndices
= true;
260 void TopologyManager::initUniformResidues(int residueSize
)
262 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
263 t_atoms
& atoms
= this->atoms();
264 int residueIndex
= -1;
265 for (int i
= 0; i
< atoms
.nr
; ++i
)
267 if (i
% residueSize
== 0)
271 atoms
.atom
[i
].resind
= residueIndex
;
273 atoms
.nres
= residueIndex
;
276 void TopologyManager::initUniformMolecules(int moleculeSize
)
278 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
279 GMX_RELEASE_ASSERT(mtop_
->molblock
.size() == 1,
280 "initUniformMolecules only implemented for a single molblock");
281 gmx_molblock_t
& molblock
= mtop_
->molblock
[0];
282 t_atoms
& atoms
= mtop_
->moltype
[molblock
.type
].atoms
;
283 GMX_RELEASE_ASSERT(atoms
.nr
% moleculeSize
== 0,
284 "The number of atoms should be a multiple of moleculeSize");
285 molblock
.nmol
= atoms
.nr
/ moleculeSize
;
286 atoms
.nr
= moleculeSize
;
287 const int nres
= atoms
.atom
[atoms
.nr
].resind
;
288 GMX_RELEASE_ASSERT(atoms
.atom
[atoms
.nr
- 1].resind
!= nres
,
289 "The residues should break at molecule boundaries");
291 mtop_
->haveMoleculeIndices
= true;
295 void TopologyManager::initFrameIndices(const ArrayRef
<const int>& index
)
297 GMX_RELEASE_ASSERT(frame_
!= nullptr, "Frame not initialized");
298 GMX_RELEASE_ASSERT(!frame_
->bIndex
, "Frame atom indices can only be set once");
300 frame_
->bIndex
= TRUE
;
301 snew(frame_
->index
, index
.size());
302 std::copy(index
.begin(), index
.end(), frame_
->index
);
304 frame_
->natoms
= index
.size();
307 t_atoms
& TopologyManager::atoms()
309 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
310 GMX_RELEASE_ASSERT(mtop_
->natoms
== mtop_
->moltype
[0].atoms
.nr
,
311 "Test setup assumes all atoms in a single molecule type");
312 return mtop_
->moltype
[0].atoms
;