src/gromacs/topology/mtop_lookup.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2016,2017,2018,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 /*! \libinternal \file
  36  *
  37  * \brief This file contains inline functions to look up atom information
  38  * using the global atom index.
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \inlibraryapi
  42  * \ingroup module_mtop
  43  */
  44
  45 #ifndef GMX_TOPOLOGY_MTOP_LOOKUP_H
  46 #define GMX_TOPOLOGY_MTOP_LOOKUP_H
  47
  48 #include "gromacs/topology/topology.h"
  49 #include "gromacs/utility/basedefinitions.h"
  50 #include "gromacs/utility/gmxassert.h"
  51
  52 struct t_atom;
  53
  54 // TODO All of the functions taking a const gmx_mtop * are deprecated
  55 // and should be replaced by versions taking const gmx_mtop & when
  56 // their callers are refactored similarly.
  57
  58 /*! \brief Look up the molecule block and other indices of a global atom index
  59  *
  60  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
  61  * The input value of moleculeBlock should be in range. Use 0 as starting value.
  62  * For subsequent calls to this function, e.g. in a loop, pass in the previously
  63  * returned value for best performance. Atoms in a group tend to be in the same
  64  * molecule(block), so this minimizes the search time.
  65  *
  66  * \param[in]     mtop                 The molecule topology
  67  * \param[in]     globalAtomIndex      The global atom index to look up
  68  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
  69  * \param[out]    moleculeIndex        The index of the molecule in the block, can be NULL
  70  * \param[out]    atomIndexInMolecule  The atom index in the molecule, can be NULL
  71  */
  72 static inline void
  73 mtopGetMolblockIndex(const gmx_mtop_t *mtop,
  74                      int               globalAtomIndex,
  75                      int              *moleculeBlock,
  76                      int              *moleculeIndex,
  77                      int              *atomIndexInMolecule)
  78 {
  79     GMX_ASSERT(globalAtomIndex >= 0 && globalAtomIndex < mtop->natoms, "The atom index to look up should be within range");
  80     GMX_ASSERT(moleculeBlock != nullptr, "molBlock can not be NULL");
  81     GMX_ASSERT(*moleculeBlock >= 0 && *moleculeBlock < gmx::ssize(mtop->moleculeBlockIndices), "The starting molecule block index for the search should be within range and moleculeBlockIndices should not be empty");
  82
  83     /* Search the molecule block index using bisection */
  84     int molBlock0 = -1;
  85     int molBlock1 = mtop->molblock.size();
  86
  87     int globalAtomStart;
  88     while (TRUE)
  89     {
  90         globalAtomStart = mtop->moleculeBlockIndices[*moleculeBlock].globalAtomStart;
  91         if (globalAtomIndex < globalAtomStart)
  92         {
  93             molBlock1 = *moleculeBlock;
  94         }
  95         else if (globalAtomIndex >= mtop->moleculeBlockIndices[*moleculeBlock].globalAtomEnd)
  96         {
  97             molBlock0 = *moleculeBlock;
  98         }
  99         else
 100         {
 101             break;
 102         }
 103         *moleculeBlock = ((molBlock0 + molBlock1 + 1) >> 1);
 104     }
 105
 106     int molIndex = (globalAtomIndex - globalAtomStart) / mtop->moleculeBlockIndices[*moleculeBlock].numAtomsPerMolecule;
 107     if (moleculeIndex != nullptr)
 108     {
 109         *moleculeIndex = molIndex;
 110     }
 111     if (atomIndexInMolecule != nullptr)
 112     {
 113         *atomIndexInMolecule = globalAtomIndex - globalAtomStart - molIndex*mtop->moleculeBlockIndices[*moleculeBlock].numAtomsPerMolecule;
 114     }
 115 }
 116
 117 /*! \brief Returns the global molecule index of a global atom index
 118  *
 119  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 120  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 121  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 122  * returned value for best performance. Atoms in a group tend to be in the same
 123  * molecule(block), so this minimizes the search time.
 124  *
 125  * \param[in]     mtop                 The molecule topology
 126  * \param[in]     globalAtomIndex      The global atom index to look up
 127  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 128  */
 129 static inline int
 130 mtopGetMoleculeIndex(const gmx_mtop_t *mtop,
 131                      int               globalAtomIndex,
 132                      int              *moleculeBlock)
 133 {
 134     int localMoleculeIndex;
 135     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, &localMoleculeIndex, nullptr);
 136
 137     return mtop->moleculeBlockIndices[*moleculeBlock].moleculeIndexStart + localMoleculeIndex;
 138 }
 139
 140 /*! \brief Returns the atom data for an atom based on global atom index
 141  *
 142  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 143  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 144  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 145  * returned value for best performance. Atoms in a group tend to be in the same
 146  * molecule(block), so this minimizes the search time.
 147  *
 148  * \param[in]     mtop                 The molecule topology
 149  * \param[in]     globalAtomIndex      The global atom index to look up
 150  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 151  */
 152 static inline const t_atom &
 153 mtopGetAtomParameters(const gmx_mtop_t  *mtop,
 154                       int                globalAtomIndex,
 155                       int               *moleculeBlock)
 156 {
 157     int atomIndexInMolecule;
 158     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock,
 159                          nullptr, &atomIndexInMolecule);
 160     const gmx_moltype_t &moltype = mtop->moltype[mtop->molblock[*moleculeBlock].type];
 161     return moltype.atoms.atom[atomIndexInMolecule];
 162 }
 163
 164 /*! \brief Returns the mass of an atom based on global atom index
 165  *
 166  * Returns that A-state mass of the atom with global index \p globalAtomIndex.
 167  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 168  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 169  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 170  * returned value for best performance. Atoms in a group tend to be in the same
 171  * molecule(block), so this minimizes the search time.
 172  *
 173  * \param[in]     mtop                 The molecule topology
 174  * \param[in]     globalAtomIndex      The global atom index to look up
 175  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 176  */
 177 static inline real
 178 mtopGetAtomMass(const gmx_mtop_t  *mtop,
 179                 int                globalAtomIndex,
 180                 int               *moleculeBlock)
 181 {
 182     const t_atom &atom = mtopGetAtomParameters(mtop, globalAtomIndex, moleculeBlock);
 183     return atom.m;
 184 }
 185
 186 /*! \brief Look up the atom and residue name and residue number and index of a global atom index
 187  *
 188  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 189  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 190  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 191  * returned value for best performance. Atoms in a group tend to be in the same
 192  * molecule(block), so this minimizes the search time.
 193  * Note that this function does a (somewhat expensive) lookup. If you want
 194  * to look up data sequentially for all atoms in a molecule or the system,
 195  * use one of the mtop loop functionalities.
 196  *
 197  * \param[in]     mtop                The molecule topology
 198  * \param[in]     globalAtomIndex     The global atom index to look up
 199  * \param[in,out] moleculeBlock       The molecule block index in \p mtop
 200  * \param[out]    atomName            The atom name, input can be NULL
 201  * \param[out]    residueNumber       The residue number, input can be NULL
 202  * \param[out]    residueName         The residue name, input can be NULL
 203  * \param[out]    globalResidueIndex  The gobal residue index, input can be NULL
 204  */
 205 static inline void
 206 mtopGetAtomAndResidueName(const gmx_mtop_t  *mtop,
 207                           int                globalAtomIndex,
 208                           int               *moleculeBlock,
 209                           const char       **atomName,
 210                           int               *residueNumber,
 211                           const char       **residueName,
 212                           int               *globalResidueIndex)
 213 {
 214     int moleculeIndex;
 215     int atomIndexInMolecule;
 216     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock,
 217                          &moleculeIndex, &atomIndexInMolecule);
 218
 219     const gmx_molblock_t       &molb    = mtop->molblock[*moleculeBlock];
 220     const t_atoms              &atoms   = mtop->moltype[molb.type].atoms;
 221     const MoleculeBlockIndices &indices = mtop->moleculeBlockIndices[*moleculeBlock];
 222     if (atomName != nullptr)
 223     {
 224         *atomName = *(atoms.atomname[atomIndexInMolecule]);
 225     }
 226     if (residueNumber != nullptr)
 227     {
 228         if (atoms.nres > mtop->maxres_renum)
 229         {
 230             *residueNumber = atoms.resinfo[atoms.atom[atomIndexInMolecule].resind].nr;
 231         }
 232         else
 233         {
 234             /* Single residue molecule, keep counting */
 235             *residueNumber = indices.residueNumberStart + moleculeIndex*atoms.nres + atoms.atom[atomIndexInMolecule].resind;
 236         }
 237     }
 238     if (residueName != nullptr)
 239     {
 240         *residueName = *(atoms.resinfo[atoms.atom[atomIndexInMolecule].resind].name);
 241     }
 242     if (globalResidueIndex != nullptr)
 243     {
 244         *globalResidueIndex = indices.globalResidueStart + moleculeIndex*atoms.nres + atoms.atom[atomIndexInMolecule].resind;
 245     }
 246 }
 247
 248 //! \copydoc mtopGetAtomAndResidueName()
 249 static inline void
 250 mtopGetAtomAndResidueName(const gmx_mtop_t  &mtop,
 251                           int                globalAtomIndex,
 252                           int               *moleculeBlock,
 253                           const char       **atomName,
 254                           int               *residueNumber,
 255                           const char       **residueName,
 256                           int               *globalResidueIndex)
 257 {
 258     mtopGetAtomAndResidueName(&mtop, globalAtomIndex, moleculeBlock,
 259                               atomName, residueNumber, residueName, globalResidueIndex);
 260 }
 261
 262 /*! \brief Returns residue information for an atom based on global atom index
 263  *
 264  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 265  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 266  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 267  * returned value for best performance. Atoms in a group tend to be in the same
 268  * molecule(block), so this minimizes the search time.
 269  *
 270  * \param[in]     mtop                 The molecule topology
 271  * \param[in]     globalAtomIndex      The global atom index to look up
 272  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 273  */
 274 static inline const t_resinfo &
 275 mtopGetResidueInfo(const gmx_mtop_t  *mtop,
 276                    int                globalAtomIndex,
 277                    int               *moleculeBlock)
 278 {
 279     int atomIndexInMolecule;
 280     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock,
 281                          nullptr, &atomIndexInMolecule);
 282     const gmx_moltype_t &moltype = mtop->moltype[mtop->molblock[*moleculeBlock].type];
 283     const int            resind  = moltype.atoms.atom[atomIndexInMolecule].resind;
 284     return moltype.atoms.resinfo[resind];
 285 }
 286
 287 /*! \brief Returns PDB information for an atom based on global atom index
 288  *
 289  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 290  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 291  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 292  * returned value for best performance. Atoms in a group tend to be in the same
 293  * molecule(block), so this minimizes the search time.
 294  *
 295  * \param[in]     mtop                 The molecule topology
 296  * \param[in]     globalAtomIndex      The global atom index to look up
 297  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 298  */
 299 static inline const t_pdbinfo &
 300 mtopGetAtomPdbInfo(const gmx_mtop_t  *mtop,
 301                    int                globalAtomIndex,
 302                    int               *moleculeBlock)
 303 {
 304     int atomIndexInMolecule;
 305     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock,
 306                          nullptr, &atomIndexInMolecule);
 307     const gmx_moltype_t &moltype = mtop->moltype[mtop->molblock[*moleculeBlock].type];
 308     GMX_ASSERT(moltype.atoms.havePdbInfo, "PDB information not present when requested");
 309     return moltype.atoms.pdbinfo[atomIndexInMolecule];
 310 }
 311
 312 #endif