2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016, 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.
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.
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.
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.
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.
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.
37 * Implements classes from atomsbuilder.h.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
43 #include "atomsbuilder.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/topology/atoms.h"
49 #include "gromacs/topology/symtab.h"
50 #include "gromacs/utility/gmxassert.h"
51 #include "gromacs/utility/smalloc.h"
56 /********************************************************************
60 AtomsBuilder::AtomsBuilder(t_atoms *atoms, t_symtab *symtab)
61 : atoms_(atoms), symtab_(symtab),
62 nrAlloc_(atoms->nr), nresAlloc_(atoms->nres),
63 currentResidueIndex_(atoms->nres), nextResidueNumber_(-1)
67 nextResidueNumber_ = atoms->resinfo[atoms->nres - 1].nr + 1;
71 AtomsBuilder::~AtomsBuilder()
75 char **AtomsBuilder::symtabString(char **source)
79 return put_symtab(symtab_, *source);
84 void AtomsBuilder::reserve(int atomCount, int residueCount)
86 srenew(atoms_->atom, atomCount);
87 srenew(atoms_->atomname, atomCount);
88 srenew(atoms_->resinfo, residueCount);
89 if (atoms_->pdbinfo != nullptr)
91 srenew(atoms_->pdbinfo, atomCount);
94 nresAlloc_ = residueCount;
97 void AtomsBuilder::clearAtoms()
101 currentResidueIndex_ = 0;
102 nextResidueNumber_ = -1;
105 int AtomsBuilder::currentAtomCount() const
110 void AtomsBuilder::setNextResidueNumber(int number)
112 nextResidueNumber_ = number;
115 void AtomsBuilder::addAtom(const t_atoms &atoms, int i)
117 const int index = atoms_->nr;
118 atoms_->atom[index] = atoms.atom[i];
119 atoms_->atomname[index] = symtabString(atoms.atomname[i]);
120 atoms_->atom[index].resind = currentResidueIndex_;
121 if (atoms_->pdbinfo != nullptr)
123 if (atoms.pdbinfo != nullptr)
125 atoms_->pdbinfo[index] = atoms.pdbinfo[i];
129 gmx_pdbinfo_init_default(&atoms_->pdbinfo[index]);
135 void AtomsBuilder::startResidue(const t_resinfo &resinfo)
137 if (nextResidueNumber_ == -1)
139 nextResidueNumber_ = resinfo.nr;
141 const int index = atoms_->nres;
142 atoms_->resinfo[index] = resinfo;
143 atoms_->resinfo[index].nr = nextResidueNumber_;
144 atoms_->resinfo[index].name = symtabString(resinfo.name);
145 ++nextResidueNumber_;
146 currentResidueIndex_ = index;
150 void AtomsBuilder::finishResidue(const t_resinfo &resinfo)
152 if (nextResidueNumber_ == -1)
154 nextResidueNumber_ = resinfo.nr;
156 const int index = currentResidueIndex_;
157 atoms_->resinfo[index] = resinfo;
158 atoms_->resinfo[index].nr = nextResidueNumber_;
159 atoms_->resinfo[index].name = symtabString(resinfo.name);
160 ++nextResidueNumber_;
161 currentResidueIndex_ = index + 1;
162 if (index >= atoms_->nres)
168 void AtomsBuilder::discardCurrentResidue()
170 int index = atoms_->nr - 1;
171 while (index > 0 && atoms_->atom[index - 1].resind == currentResidueIndex_)
176 atoms_->nres = currentResidueIndex_;
179 void AtomsBuilder::mergeAtoms(const t_atoms &atoms)
181 if (atoms_->nr + atoms.nr > nrAlloc_
182 || atoms_->nres + atoms.nres > nresAlloc_)
184 reserve(atoms_->nr + atoms.nr, atoms_->nres + atoms.nres);
187 for (int i = 0; i < atoms.nr; ++i)
189 const int resind = atoms.atom[i].resind;
190 if (resind != prevResInd)
192 startResidue(atoms.resinfo[resind]);
199 /********************************************************************
203 AtomsRemover::AtomsRemover(const t_atoms &atoms)
204 : removed_(atoms.nr, 0)
208 AtomsRemover::~AtomsRemover()
212 void AtomsRemover::refreshAtomCount(const t_atoms &atoms)
214 removed_.resize(atoms.nr, 0);
217 void AtomsRemover::markAll()
219 std::fill(removed_.begin(), removed_.end(), 1);
222 void AtomsRemover::markResidue(const t_atoms &atoms, int atomIndex, bool bStatus)
224 const int resind = atoms.atom[atomIndex].resind;
225 while (atomIndex > 0 && resind == atoms.atom[atomIndex - 1].resind)
229 while (atomIndex < atoms.nr && resind == atoms.atom[atomIndex].resind)
231 removed_[atomIndex] = (bStatus ? 1 : 0);
236 void AtomsRemover::removeMarkedElements(std::vector<RVec> *container) const
238 GMX_RELEASE_ASSERT(container->size() == removed_.size(),
239 "Mismatching contained passed for removing values");
241 for (size_t i = 0; i < removed_.size(); ++i)
245 (*container)[j] = (*container)[i];
249 container->resize(j);
252 void AtomsRemover::removeMarkedElements(std::vector<real> *container) const
254 GMX_RELEASE_ASSERT(container->size() == removed_.size(),
255 "Mismatching contained passed for removing values");
257 for (size_t i = 0; i < removed_.size(); ++i)
261 (*container)[j] = (*container)[i];
265 container->resize(j);
268 void AtomsRemover::removeMarkedAtoms(t_atoms *atoms) const
270 const int originalAtomCount = atoms->nr;
271 AtomsBuilder builder(atoms, NULL);
274 builder.setNextResidueNumber(atoms->resinfo[0].nr);
276 builder.clearAtoms();
278 for (int i = 0; i < originalAtomCount; ++i)
282 const int resind = atoms->atom[i].resind;
283 if (resind != prevResInd)
285 builder.startResidue(atoms->resinfo[resind]);
288 builder.addAtom(*atoms, i);