src/gromacs/selection/selection.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team.
   5  * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
   6  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief
  39  * Implements classes in selection.h.
  40  *
  41  * \author Teemu Murtola <teemu.murtola@gmail.com>
  42  * \ingroup module_selection
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "selection.h"
  47
  48 #include <string>
  49
  50 #include "gromacs/selection/nbsearch.h"
  51 #include "gromacs/topology/mtop_lookup.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/exceptions.h"
  54 #include "gromacs/utility/gmxassert.h"
  55 #include "gromacs/utility/stringutil.h"
  56 #include "gromacs/utility/textwriter.h"
  57
  58 #include "position.h"
  59 #include "selelem.h"
  60 #include "selvalue.h"
  61
  62 namespace gmx
  63 {
  64
  65 namespace internal
  66 {
  67
  68 /********************************************************************
  69  * SelectionData
  70  */
  71
  72 SelectionData::SelectionData(SelectionTreeElement* elem, const char* selstr) :
  73     name_(elem->name()),
  74     selectionText_(selstr),
  75     rootElement_(*elem),
  76     coveredFractionType_(CFRAC_NONE),
  77     coveredFraction_(1.0),
  78     averageCoveredFraction_(1.0),
  79     bDynamic_(false),
  80     bDynamicCoveredFraction_(false)
  81 {
  82     if (elem->child->type == SEL_CONST)
  83     {
  84         // TODO: This is not exception-safe if any called function throws.
  85         gmx_ana_pos_copy(&rawPositions_, elem->child->v.u.p, true);
  86     }
  87     else
  88     {
  89         SelectionTreeElementPointer child = elem->child;
  90         child->flags &= ~SEL_ALLOCVAL;
  91         _gmx_selvalue_setstore(&child->v, &rawPositions_);
  92         /* We should also skip any modifiers to determine the dynamic
  93          * status. */
  94         while (child->type == SEL_MODIFIER)
  95         {
  96             child = child->child;
  97             if (!child)
  98             {
  99                 break;
 100             }
 101             if (child->type == SEL_SUBEXPRREF)
 102             {
 103                 child = child->child;
 104                 /* Because most subexpression elements are created
 105                  * during compilation, we need to check for them
 106                  * explicitly here.
 107                  */
 108                 if (child->type == SEL_SUBEXPR)
 109                 {
 110                     child = child->child;
 111                 }
 112             }
 113         }
 114         if (child)
 115         {
 116             /* For variable references, we should skip the
 117              * SEL_SUBEXPRREF and SEL_SUBEXPR elements. */
 118             if (child->type == SEL_SUBEXPRREF)
 119             {
 120                 child = child->child->child;
 121             }
 122             bDynamic_ = ((child->child->flags & SEL_DYNAMIC) != 0);
 123         }
 124     }
 125     initCoveredFraction(CFRAC_NONE);
 126 }
 127
 128
 129 SelectionData::~SelectionData() {}
 130
 131
 132 bool SelectionData::initCoveredFraction(e_coverfrac_t type)
 133 {
 134     coveredFractionType_ = type;
 135     if (type == CFRAC_NONE)
 136     {
 137         bDynamicCoveredFraction_ = false;
 138     }
 139     else if (!_gmx_selelem_can_estimate_cover(rootElement()))
 140     {
 141         coveredFractionType_     = CFRAC_NONE;
 142         bDynamicCoveredFraction_ = false;
 143     }
 144     else
 145     {
 146         bDynamicCoveredFraction_ = true;
 147     }
 148     coveredFraction_        = bDynamicCoveredFraction_ ? 0.0 : 1.0;
 149     averageCoveredFraction_ = coveredFraction_;
 150     return type == CFRAC_NONE || coveredFractionType_ != CFRAC_NONE;
 151 }
 152
 153 namespace
 154 {
 155
 156 /*! \brief
 157  * Helper function to compute total masses and charges for positions.
 158  *
 159  * \param[in]  top     Topology to take atom masses from.
 160  * \param[in]  pos     Positions to compute masses and charges for.
 161  * \param[out] masses  Output masses.
 162  * \param[out] charges Output charges.
 163  *
 164  * Does not throw if enough space has been reserved for the output vectors.
 165  */
 166 void computeMassesAndCharges(const gmx_mtop_t*    top,
 167                              const gmx_ana_pos_t& pos,
 168                              std::vector<real>*   masses,
 169                              std::vector<real>*   charges)
 170 {
 171     GMX_ASSERT(top != nullptr, "Should not have been called with NULL topology");
 172     masses->clear();
 173     charges->clear();
 174     int molb = 0;
 175     for (int b = 0; b < pos.count(); ++b)
 176     {
 177         real mass   = 0.0;
 178         real charge = 0.0;
 179         for (int i = pos.m.mapb.index[b]; i < pos.m.mapb.index[b + 1]; ++i)
 180         {
 181             const int     index = pos.m.mapb.a[i];
 182             const t_atom& atom  = mtopGetAtomParameters(top, index, &molb);
 183             mass += atom.m;
 184             charge += atom.q;
 185         }
 186         masses->push_back(mass);
 187         charges->push_back(charge);
 188     }
 189 }
 190
 191 } // namespace
 192
 193 bool SelectionData::hasSortedAtomIndices() const
 194 {
 195     gmx_ana_index_t g;
 196     gmx_ana_index_set(&g, rawPositions_.m.mapb.nra, rawPositions_.m.mapb.a, -1);
 197     return gmx_ana_index_check_sorted(&g);
 198 }
 199
 200 void SelectionData::refreshName()
 201 {
 202     rootElement_.fillNameIfMissing(selectionText_.c_str());
 203     name_ = rootElement_.name();
 204 }
 205
 206 void SelectionData::initializeMassesAndCharges(const gmx_mtop_t* top)
 207 {
 208     GMX_ASSERT(posMass_.empty() && posCharge_.empty(), "Should not be called more than once");
 209     posMass_.reserve(posCount());
 210     posCharge_.reserve(posCount());
 211     if (top == nullptr)
 212     {
 213         posMass_.resize(posCount(), 1.0);
 214         posCharge_.resize(posCount(), 0.0);
 215     }
 216     else
 217     {
 218         computeMassesAndCharges(top, rawPositions_, &posMass_, &posCharge_);
 219     }
 220 }
 221
 222
 223 void SelectionData::refreshMassesAndCharges(const gmx_mtop_t* top)
 224 {
 225     if (top != nullptr && isDynamic() && !hasFlag(efSelection_DynamicMask))
 226     {
 227         computeMassesAndCharges(top, rawPositions_, &posMass_, &posCharge_);
 228     }
 229 }
 230
 231
 232 void SelectionData::updateCoveredFractionForFrame()
 233 {
 234     if (isCoveredFractionDynamic())
 235     {
 236         real cfrac       = _gmx_selelem_estimate_coverfrac(rootElement());
 237         coveredFraction_ = cfrac;
 238         averageCoveredFraction_ += cfrac;
 239     }
 240 }
 241
 242
 243 void SelectionData::computeAverageCoveredFraction(int nframes)
 244 {
 245     if (isCoveredFractionDynamic() && nframes > 0)
 246     {
 247         averageCoveredFraction_ /= nframes;
 248     }
 249 }
 250
 251
 252 void SelectionData::restoreOriginalPositions(const gmx_mtop_t* top)
 253 {
 254     if (isDynamic())
 255     {
 256         gmx_ana_pos_t& p = rawPositions_;
 257         gmx_ana_indexmap_update(&p.m, rootElement().v.u.g, hasFlag(gmx::efSelection_DynamicMask));
 258         refreshMassesAndCharges(top);
 259     }
 260 }
 261
 262 } // namespace internal
 263
 264 /********************************************************************
 265  * Selection
 266  */
 267
 268 Selection::operator AnalysisNeighborhoodPositions() const
 269 {
 270     AnalysisNeighborhoodPositions pos(data().rawPositions_.x, data().rawPositions_.count());
 271     if (hasOnlyAtoms())
 272     {
 273         pos.exclusionIds(atomIndices());
 274     }
 275     return pos;
 276 }
 277
 278
 279 void Selection::setOriginalId(int i, int id)
 280 {
 281     data().rawPositions_.m.mapid[i] = id;
 282     data().rawPositions_.m.orgid[i] = id;
 283 }
 284
 285
 286 int Selection::initOriginalIdsToGroup(const gmx_mtop_t* top, e_index_t type)
 287 {
 288     try
 289     {
 290         return gmx_ana_indexmap_init_orgid_group(&data().rawPositions_.m, top, type);
 291     }
 292     catch (const InconsistentInputError&)
 293     {
 294         GMX_ASSERT(type == INDEX_RES || type == INDEX_MOL,
 295                    "Expected that only grouping by residue/molecule would fail");
 296         std::string message = formatString(
 297                 "Cannot group selection '%s' into %s, because some "
 298                 "positions have atoms from more than one such group.",
 299                 name(), type == INDEX_MOL ? "molecules" : "residues");
 300         GMX_THROW(InconsistentInputError(message));
 301     }
 302 }
 303
 304
 305 void Selection::printInfo(FILE* fp) const
 306 {
 307     fprintf(fp, "\"%s\" (%d position%s, %d atom%s%s)", name(), posCount(), posCount() == 1 ? "" : "s",
 308             atomCount(), atomCount() == 1 ? "" : "s", isDynamic() ? ", dynamic" : "");
 309     fprintf(fp, "\n");
 310 }
 311
 312
 313 void Selection::printDebugInfo(FILE* fp, int nmaxind) const
 314 {
 315     const gmx_ana_pos_t& p = data().rawPositions_;
 316
 317     fprintf(fp, "  ");
 318     printInfo(fp);
 319     fprintf(fp, "    Group ");
 320     gmx_ana_index_t g;
 321     gmx_ana_index_set(&g, p.m.mapb.nra, p.m.mapb.a, 0);
 322     TextWriter writer(fp);
 323     gmx_ana_index_dump(&writer, &g, nmaxind);
 324
 325     fprintf(fp, "    Block (size=%d):", p.m.mapb.nr);
 326     if (!p.m.mapb.index)
 327     {
 328         fprintf(fp, " (null)");
 329     }
 330     else
 331     {
 332         int n = p.m.mapb.nr;
 333         if (nmaxind >= 0 && n > nmaxind)
 334         {
 335             n = nmaxind;
 336         }
 337         for (int i = 0; i <= n; ++i)
 338         {
 339             fprintf(fp, " %d", p.m.mapb.index[i]);
 340         }
 341         if (n < p.m.mapb.nr)
 342         {
 343             fprintf(fp, " ...");
 344         }
 345     }
 346     fprintf(fp, "\n");
 347
 348     int n = posCount();
 349     if (nmaxind >= 0 && n > nmaxind)
 350     {
 351         n = nmaxind;
 352     }
 353     fprintf(fp, "    RefId:");
 354     if (!p.m.refid)
 355     {
 356         fprintf(fp, " (null)");
 357     }
 358     else
 359     {
 360         for (int i = 0; i < n; ++i)
 361         {
 362             fprintf(fp, " %d", p.m.refid[i]);
 363         }
 364         if (n < posCount())
 365         {
 366             fprintf(fp, " ...");
 367         }
 368     }
 369     fprintf(fp, "\n");
 370
 371     fprintf(fp, "    MapId:");
 372     if (!p.m.mapid)
 373     {
 374         fprintf(fp, " (null)");
 375     }
 376     else
 377     {
 378         for (int i = 0; i < n; ++i)
 379         {
 380             fprintf(fp, " %d", p.m.mapid[i]);
 381         }
 382         if (n < posCount())
 383         {
 384             fprintf(fp, " ...");
 385         }
 386     }
 387     fprintf(fp, "\n");
 388 }
 389
 390
 391 /********************************************************************
 392  * SelectionPosition
 393  */
 394
 395 SelectionPosition::operator AnalysisNeighborhoodPositions() const
 396 {
 397     AnalysisNeighborhoodPositions pos(sel_->rawPositions_.x, sel_->rawPositions_.count());
 398     if (sel_->hasOnlyAtoms())
 399     {
 400         // TODO: Move atomIndices() such that it can be reused here as well.
 401         pos.exclusionIds(constArrayRefFromArray<int>(sel_->rawPositions_.m.mapb.a,
 402                                                      sel_->rawPositions_.m.mapb.nra));
 403     }
 404     return pos.selectSingleFromArray(i_);
 405 }
 406
 407 } // namespace gmx