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(), &ePBC
, 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
;
156 mtop_
->maxres_renum
= 0;
157 gmx_mtop_finalize(mtop_
.get());
158 GMX_RELEASE_ASSERT(mtop_
->maxres_renum
== 0,
159 "maxres_renum in mtop can be modified by an env.var., that is not supported "
161 t_atoms
& atoms
= this->atoms();
162 for (int i
= 0; i
< count
; ++i
)
164 atoms
.atom
[i
].m
= (i
% 3 == 0 ? 2.0 : 1.0);
166 atoms
.haveMass
= TRUE
;
167 if (frame_
!= nullptr)
169 frame_
->natoms
= count
;
171 snew(frame_
->x
, count
);
174 snew(frame_
->v
, count
);
178 snew(frame_
->f
, count
);
183 void TopologyManager::initAtomTypes(const ArrayRef
<const char* const>& types
)
185 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
186 atomtypes_
.reserve(types
.size());
187 for (const char* type
: types
)
189 atomtypes_
.push_back(gmx_strdup(type
));
191 t_atoms
& atoms
= this->atoms();
192 snew(atoms
.atomtype
, atoms
.nr
);
194 for (int i
= 0; i
< atoms
.nr
; ++i
, ++j
)
196 if (j
== types
.ssize())
200 atoms
.atomtype
[i
] = &atomtypes_
[j
];
202 atoms
.haveType
= TRUE
;
205 void TopologyManager::initTopology(const int numMoleculeTypes
, const int numMoleculeBlocks
)
207 GMX_RELEASE_ASSERT(mtop_
== nullptr, "Topology initialized more than once");
208 mtop_
= std::make_unique
<gmx_mtop_t
>();
209 mtop_
->moltype
.resize(numMoleculeTypes
);
210 mtop_
->molblock
.resize(numMoleculeBlocks
);
213 void TopologyManager::setMoleculeType(const int moleculeTypeIndex
, const ArrayRef
<const int> residueSizes
)
215 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
217 // Make a molecule block that refers to a new molecule type
218 auto& newMoleculeType
= mtop_
->moltype
[moleculeTypeIndex
];
219 auto& atoms
= newMoleculeType
.atoms
;
221 auto numAtomsInMolecule
= std::accumulate(residueSizes
.begin(), residueSizes
.end(), 0);
222 init_t_atoms(&atoms
, numAtomsInMolecule
, FALSE
);
223 atoms
.nres
= residueSizes
.size();
225 // Fill the residues in the molecule type
226 int residueIndex
= 0;
227 auto residueSize
= std::begin(residueSizes
);
228 for (int i
= 0; i
< atoms
.nr
&& residueSize
!= std::end(residueSizes
); ++residueSize
, ++residueIndex
)
230 for (int j
= 0; i
< atoms
.nr
&& j
< *residueSize
; ++i
, ++j
)
232 atoms
.atom
[i
].resind
= residueIndex
;
233 atoms
.atom
[i
].m
= (i
% 3 == 0 ? 2.0 : 1.0);
236 atoms
.nres
= residueIndex
;
237 atoms
.haveMass
= true;
240 void TopologyManager::setMoleculeBlock(const int moleculeBlockIndex
,
241 const int moleculeTypeIndex
,
242 const int numMoleculesToAdd
)
244 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
246 auto& newMoleculeBlock
= mtop_
->molblock
[moleculeBlockIndex
];
247 newMoleculeBlock
.type
= moleculeTypeIndex
;
248 newMoleculeBlock
.nmol
= numMoleculesToAdd
;
250 mtop_
->natoms
+= numMoleculesToAdd
* mtop_
->moltype
[moleculeTypeIndex
].atoms
.nr
;
253 void TopologyManager::finalizeTopology()
255 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
257 mtop_
->maxres_renum
= 0;
258 mtop_
->haveMoleculeIndices
= true;
259 gmx_mtop_finalize(mtop_
.get());
262 void TopologyManager::initUniformResidues(int residueSize
)
264 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
265 t_atoms
& atoms
= this->atoms();
266 int residueIndex
= -1;
267 for (int i
= 0; i
< atoms
.nr
; ++i
)
269 if (i
% residueSize
== 0)
273 atoms
.atom
[i
].resind
= residueIndex
;
275 atoms
.nres
= residueIndex
;
278 void TopologyManager::initUniformMolecules(int moleculeSize
)
280 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
281 GMX_RELEASE_ASSERT(mtop_
->molblock
.size() == 1,
282 "initUniformMolecules only implemented for a single molblock");
283 gmx_molblock_t
& molblock
= mtop_
->molblock
[0];
284 t_atoms
& atoms
= mtop_
->moltype
[molblock
.type
].atoms
;
285 GMX_RELEASE_ASSERT(atoms
.nr
% moleculeSize
== 0,
286 "The number of atoms should be a multiple of moleculeSize");
287 molblock
.nmol
= atoms
.nr
/ moleculeSize
;
288 atoms
.nr
= moleculeSize
;
289 const int nres
= atoms
.atom
[atoms
.nr
].resind
;
290 GMX_RELEASE_ASSERT(atoms
.atom
[atoms
.nr
- 1].resind
!= nres
,
291 "The residues should break at molecule boundaries");
293 mtop_
->haveMoleculeIndices
= true;
294 gmx_mtop_finalize(mtop_
.get());
297 void TopologyManager::initFrameIndices(const ArrayRef
<const int>& index
)
299 GMX_RELEASE_ASSERT(frame_
!= nullptr, "Frame not initialized");
300 GMX_RELEASE_ASSERT(!frame_
->bIndex
, "Frame atom indices can only be set once");
302 frame_
->bIndex
= TRUE
;
303 snew(frame_
->index
, index
.size());
304 std::copy(index
.begin(), index
.end(), frame_
->index
);
306 frame_
->natoms
= index
.size();
309 t_atoms
& TopologyManager::atoms()
311 GMX_RELEASE_ASSERT(mtop_
!= nullptr, "Topology not initialized");
312 GMX_RELEASE_ASSERT(mtop_
->natoms
== mtop_
->moltype
[0].atoms
.nr
,
313 "Test setup assumes all atoms in a single molecule type");
314 return mtop_
->moltype
[0].atoms
;