src/gromacs/selection/selection.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements classes in selection.h.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #include "selection.h"
  39
  40 #include "position.h"
  41 #include "selelem.h"
  42 #include "selvalue.h"
  43
  44 namespace gmx
  45 {
  46
  47 namespace internal
  48 {
  49
  50 SelectionData::SelectionData(t_selelem *elem, const char *selstr)
  51     : name_(elem->name), selectionText_(selstr),
  52       rootElement_(elem), coveredFractionType_(CFRAC_NONE),
  53       coveredFraction_(1.0), averageCoveredFraction_(1.0),
  54       bDynamic_(false), bDynamicCoveredFraction_(false)
  55 {
  56     gmx_ana_pos_clear(&rawPositions_);
  57
  58     if (elem->child->type == SEL_CONST)
  59     {
  60         // TODO: This is not exception-safe if any called function throws.
  61         gmx_ana_pos_copy(&rawPositions_, elem->child->v.u.p, true);
  62     }
  63     else
  64     {
  65         t_selelem *child;
  66
  67         child = elem->child;
  68         child->flags     &= ~SEL_ALLOCVAL;
  69         _gmx_selvalue_setstore(&child->v, &rawPositions_);
  70         /* We should also skip any modifiers to determine the dynamic
  71          * status. */
  72         while (child->type == SEL_MODIFIER)
  73         {
  74             child = child->child;
  75             if (child->type == SEL_SUBEXPRREF)
  76             {
  77                 child = child->child;
  78                 /* Because most subexpression elements are created
  79                  * during compilation, we need to check for them
  80                  * explicitly here.
  81                  */
  82                 if (child->type == SEL_SUBEXPR)
  83                 {
  84                     child = child->child;
  85                 }
  86             }
  87         }
  88         /* For variable references, we should skip the
  89          * SEL_SUBEXPRREF and SEL_SUBEXPR elements. */
  90         if (child->type == SEL_SUBEXPRREF)
  91         {
  92             child = child->child->child;
  93         }
  94         bDynamic_ = (child->child->flags & SEL_DYNAMIC);
  95     }
  96     initCoveredFraction(CFRAC_NONE);
  97 }
  98
  99
 100 SelectionData::~SelectionData()
 101 {
 102     gmx_ana_pos_deinit(&rawPositions_);
 103 }
 104
 105
 106 bool
 107 SelectionData::initCoveredFraction(e_coverfrac_t type)
 108 {
 109     coveredFractionType_ = type;
 110     if (type == CFRAC_NONE || rootElement_ == NULL)
 111     {
 112         bDynamicCoveredFraction_ = false;
 113     }
 114     else if (!_gmx_selelem_can_estimate_cover(rootElement_))
 115     {
 116         coveredFractionType_ = CFRAC_NONE;
 117         bDynamicCoveredFraction_ = false;
 118     }
 119     else
 120     {
 121         bDynamicCoveredFraction_ = true;
 122     }
 123     coveredFraction_ = bDynamicCoveredFraction_ ? 0.0 : 1.0;
 124     averageCoveredFraction_ = coveredFraction_;
 125     return type == CFRAC_NONE || coveredFractionType_ != CFRAC_NONE;
 126 }
 127
 128
 129 void
 130 SelectionData::initializeMassesAndCharges(const t_topology *top)
 131 {
 132     posInfo_.reserve(posCount());
 133     for (int b = 0; b < posCount(); ++b)
 134     {
 135         real mass   = 1.0;
 136         real charge = 0.0;
 137         if (top != NULL)
 138         {
 139             mass = 0.0;
 140             for (int i = rawPositions_.m.mapb.index[b];
 141                      i < rawPositions_.m.mapb.index[b+1];
 142                      ++i)
 143             {
 144                 int index = rawPositions_.g->index[i];
 145                 mass   += top->atoms.atom[index].m;
 146                 charge += top->atoms.atom[index].q;
 147             }
 148         }
 149         posInfo_.push_back(PositionInfo(mass, charge));
 150     }
 151     if (isDynamic() && !hasFlag(efSelection_DynamicMask))
 152     {
 153         originalPosInfo_ = posInfo_;
 154     }
 155 }
 156
 157
 158 void
 159 SelectionData::refreshMassesAndCharges()
 160 {
 161     if (!originalPosInfo_.empty())
 162     {
 163         posInfo_.clear();
 164         for (int i = 0; i < posCount(); ++i)
 165         {
 166             int refid  = rawPositions_.m.refid[i];
 167             posInfo_.push_back(originalPosInfo_[refid]);
 168         }
 169     }
 170 }
 171
 172
 173 void
 174 SelectionData::updateCoveredFractionForFrame()
 175 {
 176     if (isCoveredFractionDynamic())
 177     {
 178         real cfrac = _gmx_selelem_estimate_coverfrac(rootElement());
 179         coveredFraction_ = cfrac;
 180         averageCoveredFraction_ += cfrac;
 181     }
 182 }
 183
 184
 185 void
 186 SelectionData::computeAverageCoveredFraction(int nframes)
 187 {
 188     if (isCoveredFractionDynamic() && nframes > 0)
 189     {
 190         averageCoveredFraction_ /= nframes;
 191     }
 192 }
 193
 194
 195 void
 196 SelectionData::restoreOriginalPositions()
 197 {
 198     if (isDynamic())
 199     {
 200         gmx_ana_pos_t &p = rawPositions_;
 201         gmx_ana_index_copy(p.g, rootElement_->v.u.g, false);
 202         p.g->name = NULL;
 203         gmx_ana_indexmap_update(&p.m, p.g, hasFlag(gmx::efSelection_DynamicMask));
 204         p.nr = p.m.nr;
 205         refreshMassesAndCharges();
 206     }
 207 }
 208
 209 } // namespace internal
 210
 211
 212 void
 213 Selection::printInfo(FILE *fp) const
 214 {
 215     fprintf(fp, "\"%s\" (%d position%s, %d atom%s%s)", name(),
 216             posCount(),  posCount()  == 1 ? "" : "s",
 217             atomCount(), atomCount() == 1 ? "" : "s",
 218             isDynamic() ? ", dynamic" : "");
 219     fprintf(fp, "\n");
 220 }
 221
 222
 223 void
 224 Selection::printDebugInfo(FILE *fp, int nmaxind) const
 225 {
 226     const gmx_ana_pos_t &p = data().rawPositions_;
 227
 228     fprintf(fp, "  ");
 229     printInfo(fp);
 230     fprintf(fp, "    ");
 231     gmx_ana_index_dump(fp, p.g, -1, nmaxind);
 232
 233     fprintf(fp, "    Block (size=%d):", p.m.mapb.nr);
 234     if (!p.m.mapb.index)
 235     {
 236         fprintf(fp, " (null)");
 237     }
 238     else
 239     {
 240         int n = p.m.mapb.nr;
 241         if (nmaxind >= 0 && n > nmaxind)
 242             n = nmaxind;
 243         for (int i = 0; i <= n; ++i)
 244             fprintf(fp, " %d", p.m.mapb.index[i]);
 245         if (n < p.m.mapb.nr)
 246             fprintf(fp, " ...");
 247     }
 248     fprintf(fp, "\n");
 249
 250     int n = posCount();
 251     if (nmaxind >= 0 && n > nmaxind)
 252         n = nmaxind;
 253     fprintf(fp, "    RefId:");
 254     if (!p.m.refid)
 255     {
 256         fprintf(fp, " (null)");
 257     }
 258     else
 259     {
 260         for (int i = 0; i < n; ++i)
 261             fprintf(fp, " %d", p.m.refid[i]);
 262         if (n < posCount())
 263             fprintf(fp, " ...");
 264     }
 265     fprintf(fp, "\n");
 266
 267     fprintf(fp, "    MapId:");
 268     if (!p.m.mapid)
 269     {
 270         fprintf(fp, " (null)");
 271     }
 272     else
 273     {
 274         for (int i = 0; i < n; ++i)
 275             fprintf(fp, " %d", p.m.mapid[i]);
 276         if (n < posCount())
 277             fprintf(fp, " ...");
 278     }
 279     fprintf(fp, "\n");
 280 }
 281
 282 } // namespace gmx