src/gromacs/topology/atomsbuilder.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2016,2017,2019, 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 /*! \internal \file
  36  * \brief
  37  * Implements classes from atomsbuilder.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  */
  41 #include "gmxpre.h"
  42
  43 #include "atomsbuilder.h"
  44
  45 #include <algorithm>
  46
  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"
  52
  53 namespace gmx
  54 {
  55
  56 /********************************************************************
  57  * AtomsBuilder
  58  */
  59
  60 AtomsBuilder::AtomsBuilder(t_atoms* atoms, t_symtab* symtab) :
  61     atoms_(atoms),
  62     symtab_(symtab),
  63     nrAlloc_(atoms->nr),
  64     nresAlloc_(atoms->nres),
  65     currentResidueIndex_(atoms->nres),
  66     nextResidueNumber_(-1)
  67 {
  68     if (atoms->nres > 0)
  69     {
  70         nextResidueNumber_ = atoms->resinfo[atoms->nres - 1].nr + 1;
  71     }
  72 }
  73
  74 AtomsBuilder::~AtomsBuilder() {}
  75
  76 char** AtomsBuilder::symtabString(char** source)
  77 {
  78     if (symtab_ != nullptr)
  79     {
  80         return put_symtab(symtab_, *source);
  81     }
  82     return source;
  83 }
  84
  85 void AtomsBuilder::reserve(int atomCount, int residueCount)
  86 {
  87     srenew(atoms_->atom, atomCount);
  88     srenew(atoms_->atomname, atomCount);
  89     srenew(atoms_->resinfo, residueCount);
  90     if (atoms_->pdbinfo != nullptr)
  91     {
  92         srenew(atoms_->pdbinfo, atomCount);
  93     }
  94     nrAlloc_   = atomCount;
  95     nresAlloc_ = residueCount;
  96 }
  97
  98 void AtomsBuilder::clearAtoms()
  99 {
 100     atoms_->nr           = 0;
 101     atoms_->nres         = 0;
 102     currentResidueIndex_ = 0;
 103     nextResidueNumber_   = -1;
 104 }
 105
 106 int AtomsBuilder::currentAtomCount() const
 107 {
 108     return atoms_->nr;
 109 }
 110
 111 void AtomsBuilder::setNextResidueNumber(int number)
 112 {
 113     nextResidueNumber_ = number;
 114 }
 115
 116 void AtomsBuilder::addAtom(const t_atoms& atoms, int i)
 117 {
 118     const int index            = atoms_->nr;
 119     atoms_->atom[index]        = atoms.atom[i];
 120     atoms_->atomname[index]    = symtabString(atoms.atomname[i]);
 121     atoms_->atom[index].resind = currentResidueIndex_;
 122     if (atoms_->pdbinfo != nullptr)
 123     {
 124         if (atoms.pdbinfo != nullptr)
 125         {
 126             atoms_->pdbinfo[index] = atoms.pdbinfo[i];
 127         }
 128         else
 129         {
 130             gmx_pdbinfo_init_default(&atoms_->pdbinfo[index]);
 131         }
 132     }
 133     ++atoms_->nr;
 134 }
 135
 136 void AtomsBuilder::startResidue(const t_resinfo& resinfo)
 137 {
 138     if (nextResidueNumber_ == -1)
 139     {
 140         nextResidueNumber_ = resinfo.nr;
 141     }
 142     const int index             = atoms_->nres;
 143     atoms_->resinfo[index]      = resinfo;
 144     atoms_->resinfo[index].nr   = nextResidueNumber_;
 145     atoms_->resinfo[index].name = symtabString(resinfo.name);
 146     ++nextResidueNumber_;
 147     currentResidueIndex_ = index;
 148     ++atoms_->nres;
 149 }
 150
 151 void AtomsBuilder::finishResidue(const t_resinfo& resinfo)
 152 {
 153     if (nextResidueNumber_ == -1)
 154     {
 155         nextResidueNumber_ = resinfo.nr;
 156     }
 157     const int index             = currentResidueIndex_;
 158     atoms_->resinfo[index]      = resinfo;
 159     atoms_->resinfo[index].nr   = nextResidueNumber_;
 160     atoms_->resinfo[index].name = symtabString(resinfo.name);
 161     ++nextResidueNumber_;
 162     currentResidueIndex_ = index + 1;
 163     if (index >= atoms_->nres)
 164     {
 165         ++atoms_->nres;
 166     }
 167 }
 168
 169 void AtomsBuilder::discardCurrentResidue()
 170 {
 171     int index = atoms_->nr - 1;
 172     while (index > 0 && atoms_->atom[index - 1].resind == currentResidueIndex_)
 173     {
 174         --index;
 175     }
 176     atoms_->nr   = index;
 177     atoms_->nres = currentResidueIndex_;
 178 }
 179
 180 void AtomsBuilder::mergeAtoms(const t_atoms& atoms)
 181 {
 182     if (atoms_->nr + atoms.nr > nrAlloc_ || atoms_->nres + atoms.nres > nresAlloc_)
 183     {
 184         reserve(atoms_->nr + atoms.nr, atoms_->nres + atoms.nres);
 185     }
 186     int prevResInd = -1;
 187     for (int i = 0; i < atoms.nr; ++i)
 188     {
 189         const int resind = atoms.atom[i].resind;
 190         if (resind != prevResInd)
 191         {
 192             startResidue(atoms.resinfo[resind]);
 193             prevResInd = resind;
 194         }
 195         addAtom(atoms, i);
 196     }
 197 }
 198
 199 /********************************************************************
 200  * AtomsRemover
 201  */
 202
 203 AtomsRemover::AtomsRemover(const t_atoms& atoms) : removed_(atoms.nr, 0) {}
 204
 205 AtomsRemover::~AtomsRemover() {}
 206
 207 void AtomsRemover::refreshAtomCount(const t_atoms& atoms)
 208 {
 209     removed_.resize(atoms.nr, 0);
 210 }
 211
 212 void AtomsRemover::markAll()
 213 {
 214     std::fill(removed_.begin(), removed_.end(), 1);
 215 }
 216
 217 void AtomsRemover::markResidue(const t_atoms& atoms, int atomIndex, bool bStatus)
 218 {
 219     const int resind = atoms.atom[atomIndex].resind;
 220     while (atomIndex > 0 && resind == atoms.atom[atomIndex - 1].resind)
 221     {
 222         --atomIndex;
 223     }
 224     while (atomIndex < atoms.nr && resind == atoms.atom[atomIndex].resind)
 225     {
 226         removed_[atomIndex] = (bStatus ? 1 : 0);
 227         ++atomIndex;
 228     }
 229 }
 230
 231 void AtomsRemover::removeMarkedElements(std::vector<RVec>* container) const
 232 {
 233     GMX_RELEASE_ASSERT(container->size() == removed_.size(),
 234                        "Mismatching contained passed for removing values");
 235     int j = 0;
 236     for (size_t i = 0; i < removed_.size(); ++i)
 237     {
 238         if (!removed_[i])
 239         {
 240             (*container)[j] = (*container)[i];
 241             ++j;
 242         }
 243     }
 244     container->resize(j);
 245 }
 246
 247 void AtomsRemover::removeMarkedElements(std::vector<real>* container) const
 248 {
 249     GMX_RELEASE_ASSERT(container->size() == removed_.size(),
 250                        "Mismatching contained passed for removing values");
 251     int j = 0;
 252     for (size_t i = 0; i < removed_.size(); ++i)
 253     {
 254         if (!removed_[i])
 255         {
 256             (*container)[j] = (*container)[i];
 257             ++j;
 258         }
 259     }
 260     container->resize(j);
 261 }
 262
 263 void AtomsRemover::removeMarkedAtoms(t_atoms* atoms) const
 264 {
 265     const int    originalAtomCount = atoms->nr;
 266     AtomsBuilder builder(atoms, nullptr);
 267     if (atoms->nr > 0)
 268     {
 269         builder.setNextResidueNumber(atoms->resinfo[0].nr);
 270     }
 271     builder.clearAtoms();
 272     int prevResInd = -1;
 273     for (int i = 0; i < originalAtomCount; ++i)
 274     {
 275         if (!removed_[i])
 276         {
 277             const int resind = atoms->atom[i].resind;
 278             if (resind != prevResInd)
 279             {
 280                 builder.startResidue(atoms->resinfo[resind]);
 281                 prevResInd = resind;
 282             }
 283             builder.addAtom(*atoms, i);
 284         }
 285     }
 286 }
 287
 288 } // namespace gmx