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,2020,2021, 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 mtopGetMolblockIndex(const gmx_mtop_t* mtop,
  73                                         int               globalAtomIndex,
  74                                         int*              moleculeBlock,
  75                                         int*              moleculeIndex,
  76                                         int*              atomIndexInMolecule)
  77 {
  78     GMX_ASSERT(globalAtomIndex >= 0, "The atom index to look up should not be negative");
  79     GMX_ASSERT(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(!mtop->moleculeBlockIndices.empty(), "The moleculeBlockIndices should not be empty");
  82     GMX_ASSERT(*moleculeBlock >= 0,
  83                "The starting molecule block index for the search should not be negative");
  84     GMX_ASSERT(*moleculeBlock < gmx::ssize(mtop->moleculeBlockIndices),
  85                "The starting molecule block index for the search should be within range");
  86
  87     /* Search the molecule block index using bisection */
  88     int molBlock0 = -1;
  89     int molBlock1 = mtop->molblock.size();
  90
  91     int globalAtomStart = 0;
  92     while (TRUE)
  93     {
  94         globalAtomStart = mtop->moleculeBlockIndices[*moleculeBlock].globalAtomStart;
  95         if (globalAtomIndex < globalAtomStart)
  96         {
  97             molBlock1 = *moleculeBlock;
  98         }
  99         else if (globalAtomIndex >= mtop->moleculeBlockIndices[*moleculeBlock].globalAtomEnd)
 100         {
 101             molBlock0 = *moleculeBlock;
 102         }
 103         else
 104         {
 105             break;
 106         }
 107         *moleculeBlock = ((molBlock0 + molBlock1 + 1) >> 1);
 108     }
 109
 110     int molIndex = (globalAtomIndex - globalAtomStart)
 111                    / mtop->moleculeBlockIndices[*moleculeBlock].numAtomsPerMolecule;
 112     if (moleculeIndex != nullptr)
 113     {
 114         *moleculeIndex = molIndex;
 115     }
 116     if (atomIndexInMolecule != nullptr)
 117     {
 118         *atomIndexInMolecule = globalAtomIndex - globalAtomStart
 119                                - molIndex * mtop->moleculeBlockIndices[*moleculeBlock].numAtomsPerMolecule;
 120     }
 121 }
 122
 123 /*! \brief Returns the global molecule index of a global atom index
 124  *
 125  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 126  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 127  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 128  * returned value for best performance. Atoms in a group tend to be in the same
 129  * molecule(block), so this minimizes the search time.
 130  *
 131  * \param[in]     mtop                 The molecule topology
 132  * \param[in]     globalAtomIndex      The global atom index to look up
 133  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 134  */
 135 static inline int mtopGetMoleculeIndex(const gmx_mtop_t* mtop, int globalAtomIndex, int* moleculeBlock)
 136 {
 137     int localMoleculeIndex = 0;
 138     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, &localMoleculeIndex, nullptr);
 139
 140     return mtop->moleculeBlockIndices[*moleculeBlock].moleculeIndexStart + localMoleculeIndex;
 141 }
 142
 143 /*! \brief Returns the atom data for an atom based on global atom index
 144  *
 145  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 146  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 147  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 148  * returned value for best performance. Atoms in a group tend to be in the same
 149  * molecule(block), so this minimizes the search time.
 150  *
 151  * \param[in]     mtop                 The molecule topology
 152  * \param[in]     globalAtomIndex      The global atom index to look up
 153  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 154  */
 155 static inline const t_atom& mtopGetAtomParameters(const gmx_mtop_t* mtop, int globalAtomIndex, int* moleculeBlock)
 156 {
 157     int atomIndexInMolecule = 0;
 158     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, nullptr, &atomIndexInMolecule);
 159     const gmx_moltype_t& moltype = mtop->moltype[mtop->molblock[*moleculeBlock].type];
 160     return moltype.atoms.atom[atomIndexInMolecule];
 161 }
 162
 163 /*! \brief Returns the mass of an atom based on global atom index
 164  *
 165  * Returns that A-state mass of the atom with global index \p globalAtomIndex.
 166  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 167  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 168  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 169  * returned value for best performance. Atoms in a group tend to be in the same
 170  * molecule(block), so this minimizes the search time.
 171  *
 172  * \param[in]     mtop                 The molecule topology
 173  * \param[in]     globalAtomIndex      The global atom index to look up
 174  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 175  */
 176 static inline real mtopGetAtomMass(const gmx_mtop_t* mtop, int globalAtomIndex, int* moleculeBlock)
 177 {
 178     const t_atom& atom = mtopGetAtomParameters(mtop, globalAtomIndex, moleculeBlock);
 179     return atom.m;
 180 }
 181
 182 /*! \brief Look up the atom and residue name and residue number and index of a global atom index
 183  *
 184  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 185  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 186  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 187  * returned value for best performance. Atoms in a group tend to be in the same
 188  * molecule(block), so this minimizes the search time.
 189  * Note that this function does a (somewhat expensive) lookup. If you want
 190  * to look up data sequentially for all atoms in a molecule or the system,
 191  * use one of the mtop loop functionalities.
 192  *
 193  * \param[in]     mtop                The molecule topology
 194  * \param[in]     globalAtomIndex     The global atom index to look up
 195  * \param[in,out] moleculeBlock       The molecule block index in \p mtop
 196  * \param[out]    atomName            The atom name, input can be NULL
 197  * \param[out]    residueNumber       The residue number, input can be NULL
 198  * \param[out]    residueName         The residue name, input can be NULL
 199  * \param[out]    globalResidueIndex  The gobal residue index, input can be NULL
 200  */
 201 static inline void mtopGetAtomAndResidueName(const gmx_mtop_t* mtop,
 202                                              int               globalAtomIndex,
 203                                              int*              moleculeBlock,
 204                                              const char**      atomName,
 205                                              int*              residueNumber,
 206                                              const char**      residueName,
 207                                              int*              globalResidueIndex)
 208 {
 209     int moleculeIndex       = 0;
 210     int atomIndexInMolecule = 0;
 211     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, &moleculeIndex, &atomIndexInMolecule);
 212
 213     const gmx_molblock_t&       molb    = mtop->molblock[*moleculeBlock];
 214     const t_atoms&              atoms   = mtop->moltype[molb.type].atoms;
 215     const MoleculeBlockIndices& indices = mtop->moleculeBlockIndices[*moleculeBlock];
 216     if (atomName != nullptr)
 217     {
 218         *atomName = *(atoms.atomname[atomIndexInMolecule]);
 219     }
 220     if (residueNumber != nullptr)
 221     {
 222         if (atoms.nres > mtop->maxResiduesPerMoleculeToTriggerRenumber())
 223         {
 224             *residueNumber = atoms.resinfo[atoms.atom[atomIndexInMolecule].resind].nr;
 225         }
 226         else
 227         {
 228             /* Single residue molecule, keep counting */
 229             *residueNumber = indices.residueNumberStart + moleculeIndex * atoms.nres
 230                              + atoms.atom[atomIndexInMolecule].resind;
 231         }
 232     }
 233     if (residueName != nullptr)
 234     {
 235         *residueName = *(atoms.resinfo[atoms.atom[atomIndexInMolecule].resind].name);
 236     }
 237     if (globalResidueIndex != nullptr)
 238     {
 239         *globalResidueIndex = indices.globalResidueStart + moleculeIndex * atoms.nres
 240                               + atoms.atom[atomIndexInMolecule].resind;
 241     }
 242 }
 243
 244 //! \copydoc mtopGetAtomAndResidueName()
 245 static inline void mtopGetAtomAndResidueName(const gmx_mtop_t& mtop,
 246                                              int               globalAtomIndex,
 247                                              int*              moleculeBlock,
 248                                              const char**      atomName,
 249                                              int*              residueNumber,
 250                                              const char**      residueName,
 251                                              int*              globalResidueIndex)
 252 {
 253     mtopGetAtomAndResidueName(
 254             &mtop, globalAtomIndex, moleculeBlock, atomName, residueNumber, residueName, globalResidueIndex);
 255 }
 256
 257 /*! \brief Returns residue information for an atom based on global atom index
 258  *
 259  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 260  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 261  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 262  * returned value for best performance. Atoms in a group tend to be in the same
 263  * molecule(block), so this minimizes the search time.
 264  *
 265  * \param[in]     mtop                 The molecule topology
 266  * \param[in]     globalAtomIndex      The global atom index to look up
 267  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 268  */
 269 static inline const t_resinfo& mtopGetResidueInfo(const gmx_mtop_t* mtop, int globalAtomIndex, int* moleculeBlock)
 270 {
 271     int atomIndexInMolecule = 0;
 272     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, nullptr, &atomIndexInMolecule);
 273     const gmx_moltype_t& moltype = mtop->moltype[mtop->molblock[*moleculeBlock].type];
 274     const int            resind  = moltype.atoms.atom[atomIndexInMolecule].resind;
 275     return moltype.atoms.resinfo[resind];
 276 }
 277
 278 /*! \brief Returns PDB information for an atom based on global atom index
 279  *
 280  * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms.
 281  * The input value of moleculeBlock should be in range. Use 0 as starting value.
 282  * For subsequent calls to this function, e.g. in a loop, pass in the previously
 283  * returned value for best performance. Atoms in a group tend to be in the same
 284  * molecule(block), so this minimizes the search time.
 285  *
 286  * \param[in]     mtop                 The molecule topology
 287  * \param[in]     globalAtomIndex      The global atom index to look up
 288  * \param[in,out] moleculeBlock        The molecule block index in \p mtop
 289  */
 290 static inline const t_pdbinfo& mtopGetAtomPdbInfo(const gmx_mtop_t* mtop, int globalAtomIndex, int* moleculeBlock)
 291 {
 292     int atomIndexInMolecule = 0;
 293     mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, nullptr, &atomIndexInMolecule);
 294     const gmx_moltype_t& moltype = mtop->moltype[mtop->molblock[*moleculeBlock].type];
 295     GMX_ASSERT(moltype.atoms.havePdbInfo, "PDB information not present when requested");
 296     return moltype.atoms.pdbinfo[atomIndexInMolecule];
 297 }
 298
 299 #endif