src/gromacs/topology/mtop_lookup.h

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