src/gromacs/selection/selelem.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014, 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 /*! \internal \file
  36  * \brief
  37  * Implements functions in selelem.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "selelem.h"
  45
  46 #include <cstring>
  47
  48 #include "gromacs/selection/indexutil.h"
  49 #include "gromacs/selection/position.h"
  50 #include "gromacs/utility/exceptions.h"
  51 #include "gromacs/utility/gmxassert.h"
  52 #include "gromacs/utility/smalloc.h"
  53 #include "gromacs/utility/stringutil.h"
  54
  55 #include "keywords.h"
  56 #include "mempool.h"
  57 #include "poscalc.h"
  58 #include "selmethod.h"
  59
  60 /*!
  61  * \param[in] sel Selection for which the string is requested
  62  * \returns   Pointer to a string that corresponds to \p sel->type.
  63  *
  64  * The return value points to a string constant and should not be \p free'd.
  65  *
  66  * The function returns NULL if \p sel->type is not one of the valid values.
  67  */
  68 const char *
  69 _gmx_selelem_type_str(const gmx::SelectionTreeElement &sel)
  70 {
  71     const char *p = NULL;
  72     switch (sel.type)
  73     {
  74         case SEL_CONST:      p = "CONST";    break;
  75         case SEL_EXPRESSION: p = "EXPR";     break;
  76         case SEL_BOOLEAN:    p = "BOOL";     break;
  77         case SEL_ARITHMETIC: p = "ARITH";    break;
  78         case SEL_ROOT:       p = "ROOT";     break;
  79         case SEL_SUBEXPR:    p = "SUBEXPR";  break;
  80         case SEL_SUBEXPRREF: p = "REF";      break;
  81         case SEL_GROUPREF:   p = "GROUPREF"; break;
  82         case SEL_MODIFIER:   p = "MODIFIER"; break;
  83             // No default clause so we intentionally get compiler errors
  84             // if new selection choices are added later.
  85     }
  86     return p;
  87 }
  88
  89 /*!
  90  * \param[in] val Value structore for which the string is requested.
  91  * \returns   Pointer to a string that corresponds to \p val->type,
  92  *   NULL if the type value is invalid.
  93  *
  94  * The return value points to a string constant and should not be \p free'd.
  95  */
  96 const char *
  97 _gmx_sel_value_type_str(const gmx_ana_selvalue_t *val)
  98 {
  99     const char *p = NULL;
 100     switch (val->type)
 101     {
 102         case NO_VALUE:       p = "NONE";  break;
 103         case INT_VALUE:      p = "INT";   break;
 104         case REAL_VALUE:     p = "REAL";  break;
 105         case STR_VALUE:      p = "STR";   break;
 106         case POS_VALUE:      p = "VEC";   break;
 107         case GROUP_VALUE:    p = "GROUP"; break;
 108             // No default clause so we intentionally get compiler errors
 109             // if new selection choices are added later.
 110     }
 111     return p;
 112 }
 113
 114 /*! \copydoc _gmx_selelem_type_str() */
 115 const char *
 116 _gmx_selelem_boolean_type_str(const gmx::SelectionTreeElement &sel)
 117 {
 118     const char *p = NULL;
 119     switch (sel.u.boolt)
 120     {
 121         case BOOL_NOT:  p = "NOT"; break;
 122         case BOOL_AND:  p = "AND"; break;
 123         case BOOL_OR:   p = "OR";  break;
 124         case BOOL_XOR:  p = "XOR"; break;
 125             // No default clause so we intentionally get compiler errors
 126             // if new selection choices are added later.
 127     }
 128     return p;
 129 }
 130
 131
 132 namespace gmx
 133 {
 134
 135 SelectionTreeElement::SelectionTreeElement(e_selelem_t type)
 136 {
 137     this->type       = type;
 138     this->flags      = (type != SEL_ROOT) ? SEL_ALLOCVAL : 0;
 139     if (type == SEL_BOOLEAN)
 140     {
 141         this->v.type = GROUP_VALUE;
 142         this->flags |= SEL_ALLOCDATA;
 143     }
 144     else
 145     {
 146         this->v.type = NO_VALUE;
 147     }
 148     _gmx_selvalue_clear(&this->v);
 149     std::memset(&this->u, 0, sizeof(this->u));
 150     this->evaluate   = NULL;
 151     this->mempool    = NULL;
 152     this->cdata      = NULL;
 153 }
 154
 155 SelectionTreeElement::~SelectionTreeElement()
 156 {
 157     /* Free the children.
 158      * Must be done before freeing other data, because the children may hold
 159      * references to data in this element. */
 160     child.reset();
 161
 162     freeValues();
 163     freeExpressionData();
 164     freeCompilerData();
 165 }
 166
 167 void SelectionTreeElement::freeValues()
 168 {
 169     mempoolRelease();
 170     if ((flags & SEL_ALLOCDATA) && v.u.ptr)
 171     {
 172         /* The number of position/group structures is constant, so the
 173          * backup of using sel->v.nr should work for them.
 174          * For strings, we report an error if we don't know the allocation
 175          * size here. */
 176         int n = (v.nalloc > 0) ? v.nalloc : v.nr;
 177         switch (v.type)
 178         {
 179             case STR_VALUE:
 180                 GMX_RELEASE_ASSERT(v.nalloc != 0,
 181                                    "SEL_ALLOCDATA should only be set for allocated "
 182                                    "STR_VALUE values");
 183                 for (int i = 0; i < n; ++i)
 184                 {
 185                     sfree(v.u.s[i]);
 186                 }
 187                 break;
 188             case GROUP_VALUE:
 189                 for (int i = 0; i < n; ++i)
 190                 {
 191                     gmx_ana_index_deinit(&v.u.g[i]);
 192                 }
 193                 break;
 194             default: /* No special handling for other types */
 195                 break;
 196         }
 197     }
 198     if (v.nalloc > 0)
 199     {
 200         if (v.type == POS_VALUE)
 201         {
 202             delete [] v.u.p;
 203         }
 204         else
 205         {
 206             sfree(v.u.ptr);
 207         }
 208     }
 209     _gmx_selvalue_setstore(&v, NULL);
 210     if (type == SEL_SUBEXPRREF && u.param)
 211     {
 212         u.param->val.u.ptr = NULL;
 213     }
 214 }
 215
 216 void
 217 SelectionTreeElement::freeExpressionData()
 218 {
 219     if (type == SEL_EXPRESSION || type == SEL_MODIFIER)
 220     {
 221         _gmx_selelem_free_method(u.expr.method, u.expr.mdata);
 222         u.expr.mdata  = NULL;
 223         u.expr.method = NULL;
 224         /* Free position data */
 225         delete u.expr.pos;
 226         u.expr.pos = NULL;
 227         /* Free position calculation data */
 228         if (u.expr.pc)
 229         {
 230             gmx_ana_poscalc_free(u.expr.pc);
 231             u.expr.pc = NULL;
 232         }
 233     }
 234     if (type == SEL_ARITHMETIC)
 235     {
 236         sfree(u.arith.opstr);
 237         u.arith.opstr = NULL;
 238     }
 239     if (type == SEL_SUBEXPR || type == SEL_ROOT
 240         || (type == SEL_CONST && v.type == GROUP_VALUE))
 241     {
 242         gmx_ana_index_deinit(&u.cgrp);
 243     }
 244     if (type == SEL_GROUPREF)
 245     {
 246         sfree(u.gref.name);
 247     }
 248 }
 249
 250 void SelectionTreeElement::mempoolReserve(int count)
 251 {
 252     if (!mempool)
 253     {
 254         return;
 255     }
 256     switch (v.type)
 257     {
 258         case INT_VALUE:
 259             v.u.i = static_cast<int *>(
 260                     _gmx_sel_mempool_alloc(mempool, sizeof(*v.u.i)*count));
 261             break;
 262
 263         case REAL_VALUE:
 264             v.u.r = static_cast<real *>(
 265                     _gmx_sel_mempool_alloc(mempool, sizeof(*v.u.r)*count));
 266             break;
 267
 268         case GROUP_VALUE:
 269             _gmx_sel_mempool_alloc_group(mempool, v.u.g, count);
 270             break;
 271
 272         default:
 273             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
 274     }
 275 }
 276
 277 void SelectionTreeElement::mempoolRelease()
 278 {
 279     if (!mempool)
 280     {
 281         return;
 282     }
 283     switch (v.type)
 284     {
 285         case INT_VALUE:
 286         case REAL_VALUE:
 287             _gmx_sel_mempool_free(mempool, v.u.ptr);
 288             _gmx_selvalue_setstore(&v, NULL);
 289             break;
 290
 291         case GROUP_VALUE:
 292             if (v.u.g)
 293             {
 294                 _gmx_sel_mempool_free_group(mempool, v.u.g);
 295             }
 296             break;
 297
 298         default:
 299             GMX_THROW(gmx::InternalError("Memory pooling not implemented for requested type"));
 300     }
 301 }
 302
 303 void SelectionTreeElement::fillNameIfMissing(const char *selectionText)
 304 {
 305     GMX_RELEASE_ASSERT(type == SEL_ROOT,
 306                        "Should not be called for non-root elements");
 307     if (name().empty())
 308     {
 309         // Check whether the actual selection given was from an external group,
 310         // and if so, use the name of the external group.
 311         SelectionTreeElementPointer child = this->child;
 312         while (child->type == SEL_MODIFIER)
 313         {
 314             if (!child->child || child->child->type != SEL_SUBEXPRREF
 315                 || !child->child->child)
 316             {
 317                 break;
 318             }
 319             child = child->child->child;
 320         }
 321         if (child->type == SEL_EXPRESSION
 322             && child->child && child->child->type == SEL_SUBEXPRREF
 323             && child->child->child)
 324         {
 325             if (child->child->child->type == SEL_CONST
 326                 && child->child->child->v.type == GROUP_VALUE)
 327             {
 328                 setName(child->child->child->name());
 329                 return;
 330             }
 331             // If the group reference is still unresolved, leave the name empty
 332             // and fill it later.
 333             if (child->child->child->type == SEL_GROUPREF)
 334             {
 335                 return;
 336             }
 337         }
 338         // If there still is no name, use the selection string.
 339         setName(selectionText);
 340     }
 341 }
 342
 343 void SelectionTreeElement::checkUnsortedAtoms(
 344         bool bUnsortedAllowed, ExceptionInitializer *errors) const
 345 {
 346     const bool bUnsortedSupported
 347         = (type == SEL_CONST && v.type == GROUP_VALUE)
 348             || type == SEL_ROOT || type == SEL_SUBEXPR || type == SEL_SUBEXPRREF
 349             // TODO: Consolidate.
 350             || type == SEL_MODIFIER
 351             || (type == SEL_EXPRESSION && (u.expr.method->flags & SMETH_ALLOW_UNSORTED));
 352
 353     // TODO: For some complicated selections, this may result in the same
 354     // index group reference being flagged as an error multiple times for the
 355     // same selection.
 356     SelectionTreeElementPointer child = this->child;
 357     while (child)
 358     {
 359         child->checkUnsortedAtoms(bUnsortedAllowed && bUnsortedSupported,
 360                                   errors);
 361         child = child->next;
 362     }
 363
 364     // The logic here is simplified by the fact that only constant groups can
 365     // currently be the root cause of SEL_UNSORTED being set, so only those
 366     // need to be considered in triggering the error.
 367     if (!bUnsortedAllowed && (flags & SEL_UNSORTED)
 368         && type == SEL_CONST && v.type == GROUP_VALUE)
 369     {
 370         std::string message = formatString(
 371                     "Group '%s' cannot be used in selections except "
 372                     "as a full value of the selection, "
 373                     "because atom indices in it are not sorted and/or "
 374                     "it contains duplicate atoms.",
 375                     name().c_str());
 376         errors->addNested(InconsistentInputError(message));
 377     }
 378 }
 379
 380 void SelectionTreeElement::resolveIndexGroupReference(
 381         gmx_ana_indexgrps_t *grps, int natoms)
 382 {
 383     GMX_RELEASE_ASSERT(type == SEL_GROUPREF,
 384                        "Should only be called for index group reference elements");
 385     if (grps == NULL)
 386     {
 387         std::string message = formatString(
 388                     "Cannot match '%s', because index groups are not available.",
 389                     name().c_str());
 390         GMX_THROW(InconsistentInputError(message));
 391     }
 392
 393     gmx_ana_index_t foundGroup;
 394     std::string     foundName;
 395     if (u.gref.name != NULL)
 396     {
 397         if (!gmx_ana_indexgrps_find(&foundGroup, &foundName, grps, u.gref.name))
 398         {
 399             std::string message = formatString(
 400                         "Cannot match '%s', because no such index group can be found.",
 401                         name().c_str());
 402             GMX_THROW(InconsistentInputError(message));
 403         }
 404     }
 405     else
 406     {
 407         if (!gmx_ana_indexgrps_extract(&foundGroup, &foundName, grps, u.gref.id))
 408         {
 409             std::string message = formatString(
 410                         "Cannot match '%s', because no such index group can be found.",
 411                         name().c_str());
 412             GMX_THROW(InconsistentInputError(message));
 413         }
 414     }
 415
 416     if (!gmx_ana_index_check_sorted(&foundGroup))
 417     {
 418         flags |= SEL_UNSORTED;
 419     }
 420
 421     sfree(u.gref.name);
 422     type = SEL_CONST;
 423     gmx_ana_index_set(&u.cgrp, foundGroup.isize, foundGroup.index,
 424                       foundGroup.nalloc_index);
 425     setName(foundName);
 426
 427     if (natoms > 0)
 428     {
 429         checkIndexGroup(natoms);
 430     }
 431 }
 432
 433 void SelectionTreeElement::checkIndexGroup(int natoms)
 434 {
 435     GMX_RELEASE_ASSERT(type == SEL_CONST && v.type == GROUP_VALUE,
 436                        "Should only be called for index group elements");
 437     if (!gmx_ana_index_check_range(&u.cgrp, natoms))
 438     {
 439         std::string message = formatString(
 440                     "Group '%s' cannot be used in selections, because it "
 441                     "contains negative atom indices and/or references atoms "
 442                     "not present (largest allowed atom index is %d).",
 443                     name().c_str(), natoms);
 444         GMX_THROW(InconsistentInputError(message));
 445     }
 446 }
 447
 448 } // namespace gmx
 449
 450 /*!
 451  * \param[in,out] sel   Selection element to set the type for.
 452  * \param[in]     vtype Value type for the selection element.
 453  *
 454  * If the new type is \ref GROUP_VALUE or \ref POS_VALUE, the
 455  * \ref SEL_ALLOCDATA flag is also set.
 456  *
 457  * This function should only be called at most once for each element,
 458  * preferably right after calling _gmx_selelem_create().
 459  */
 460 void
 461 _gmx_selelem_set_vtype(const gmx::SelectionTreeElementPointer &sel,
 462                        e_selvalue_t                            vtype)
 463 {
 464     GMX_RELEASE_ASSERT(sel->type != SEL_BOOLEAN || vtype == GROUP_VALUE,
 465                        "Boolean elements must have a group value");
 466     GMX_RELEASE_ASSERT(sel->v.type == NO_VALUE || vtype == sel->v.type,
 467                        "_gmx_selelem_set_vtype() called more than once");
 468     sel->v.type = vtype;
 469     if (vtype == GROUP_VALUE || vtype == POS_VALUE)
 470     {
 471         sel->flags |= SEL_ALLOCDATA;
 472     }
 473 }
 474
 475 /*!
 476  * \param[in] method Method to free.
 477  * \param[in] mdata  Method data to free.
 478  */
 479 void
 480 _gmx_selelem_free_method(gmx_ana_selmethod_t *method, void *mdata)
 481 {
 482     sel_freefunc free_func = NULL;
 483
 484     /* Save the pointer to the free function. */
 485     if (method && method->free)
 486     {
 487         free_func = method->free;
 488     }
 489
 490     /* Free the method itself.
 491      * Has to be done before freeing the method data, because parameter
 492      * values are typically stored in the method data, and here we may
 493      * access them. */
 494     if (method)
 495     {
 496         int  i, j;
 497
 498         /* Free the memory allocated for the parameters that are not managed
 499          * by the selection method itself. */
 500         for (i = 0; i < method->nparams; ++i)
 501         {
 502             gmx_ana_selparam_t *param = &method->param[i];
 503
 504             if (param->val.u.ptr)
 505             {
 506                 if (param->val.type == GROUP_VALUE)
 507                 {
 508                     for (j = 0; j < param->val.nr; ++j)
 509                     {
 510                         gmx_ana_index_deinit(&param->val.u.g[j]);
 511                     }
 512                 }
 513                 else if (param->val.type == POS_VALUE)
 514                 {
 515                     if (param->val.nalloc > 0)
 516                     {
 517                         delete[] param->val.u.p;
 518                         _gmx_selvalue_setstore(&param->val, NULL);
 519                     }
 520                 }
 521
 522                 if (param->val.nalloc > 0)
 523                 {
 524                     sfree(param->val.u.ptr);
 525                 }
 526             }
 527         }
 528         sfree(method->param);
 529         sfree(method);
 530     }
 531     /* Free method data. */
 532     if (mdata)
 533     {
 534         if (free_func)
 535         {
 536             free_func(mdata);
 537         }
 538         else
 539         {
 540             sfree(mdata);
 541         }
 542     }
 543 }
 544
 545 /*!
 546  * \param[in] fp      File handle to receive the output.
 547  * \param[in] sel     Root of the selection subtree to print.
 548  * \param[in] bValues If true, the evaluated values of selection elements
 549  *   are printed as well.
 550  * \param[in] level   Indentation level, starting from zero.
 551  */
 552 void
 553 _gmx_selelem_print_tree(FILE *fp, const gmx::SelectionTreeElement &sel,
 554                         bool bValues, int level)
 555 {
 556     int          i;
 557
 558     fprintf(fp, "%*c %s %s", level*2+1, '*',
 559             _gmx_selelem_type_str(sel), _gmx_sel_value_type_str(&sel.v));
 560     if (!sel.name().empty())
 561     {
 562         fprintf(fp, " \"%s\"", sel.name().c_str());
 563     }
 564     fprintf(fp, " flg=");
 565     if (sel.flags & SEL_FLAGSSET)
 566     {
 567         fprintf(fp, "s");
 568     }
 569     if (sel.flags & SEL_SINGLEVAL)
 570     {
 571         fprintf(fp, "S");
 572     }
 573     if (sel.flags & SEL_ATOMVAL)
 574     {
 575         fprintf(fp, "A");
 576     }
 577     if (sel.flags & SEL_VARNUMVAL)
 578     {
 579         fprintf(fp, "V");
 580     }
 581     if (sel.flags & SEL_DYNAMIC)
 582     {
 583         fprintf(fp, "D");
 584     }
 585     if (!(sel.flags & SEL_VALFLAGMASK))
 586     {
 587         fprintf(fp, "0");
 588     }
 589     if (sel.flags & SEL_ALLOCVAL)
 590     {
 591         fprintf(fp, "Av");
 592     }
 593     if (sel.flags & SEL_ALLOCDATA)
 594     {
 595         fprintf(fp, "Ad");
 596     }
 597     if (sel.mempool)
 598     {
 599         fprintf(fp, "P");
 600     }
 601     if (sel.type == SEL_CONST)
 602     {
 603         if (sel.v.type == INT_VALUE)
 604         {
 605             fprintf(fp, " %d", sel.v.u.i[0]);
 606         }
 607         else if (sel.v.type == REAL_VALUE)
 608         {
 609             fprintf(fp, " %f", sel.v.u.r[0]);
 610         }
 611         else if (sel.v.type == GROUP_VALUE)
 612         {
 613             const gmx_ana_index_t *g = sel.v.u.g;
 614             if (!g || g->isize == 0)
 615             {
 616                 g = &sel.u.cgrp;
 617             }
 618             fprintf(fp, " (%d atoms)", g->isize);
 619         }
 620     }
 621     else if (sel.type == SEL_BOOLEAN)
 622     {
 623         fprintf(fp, " %s", _gmx_selelem_boolean_type_str(sel));
 624     }
 625     else if (sel.type == SEL_EXPRESSION
 626              && sel.u.expr.method->name == sm_compare.name)
 627     {
 628         _gmx_selelem_print_compare_info(fp, sel.u.expr.mdata);
 629     }
 630     if (sel.evaluate)
 631     {
 632         fprintf(fp, " eval=");
 633         _gmx_sel_print_evalfunc_name(fp, sel.evaluate);
 634     }
 635     if (!(sel.flags & SEL_ALLOCVAL))
 636     {
 637         fprintf(fp, " (ext)");
 638     }
 639     fprintf(fp, "\n");
 640
 641     if ((sel.type == SEL_CONST && sel.v.type == GROUP_VALUE) || sel.type == SEL_ROOT)
 642     {
 643         const gmx_ana_index_t *g = sel.v.u.g;
 644         if (!g || g->isize == 0 || sel.evaluate != NULL)
 645         {
 646             g = &sel.u.cgrp;
 647         }
 648         if (g->isize < 0)
 649         {
 650             fprintf(fp, "%*c group: (null)\n", level*2+1, ' ');
 651         }
 652         else if (g->isize > 0)
 653         {
 654             fprintf(fp, "%*c group:", level*2+1, ' ');
 655             if (g->isize <= 20)
 656             {
 657                 for (i = 0; i < g->isize; ++i)
 658                 {
 659                     fprintf(fp, " %d", g->index[i] + 1);
 660                 }
 661             }
 662             else
 663             {
 664                 fprintf(fp, " %d atoms", g->isize);
 665             }
 666             fprintf(fp, "\n");
 667         }
 668     }
 669     else if (sel.type == SEL_EXPRESSION)
 670     {
 671         if (sel.u.expr.pc)
 672         {
 673             fprintf(fp, "%*c COM", level*2+3, '*');
 674             fprintf(fp, "\n");
 675         }
 676     }
 677
 678     if (sel.cdata)
 679     {
 680         _gmx_selelem_print_compiler_info(fp, sel, level);
 681     }
 682
 683     if (bValues && sel.type != SEL_CONST && sel.type != SEL_ROOT && sel.v.u.ptr)
 684     {
 685         fprintf(fp, "%*c value: ", level*2+1, ' ');
 686         switch (sel.v.type)
 687         {
 688             case POS_VALUE:
 689                 /* In normal use, the pointer should never be NULL, but it's
 690                  * useful to have the check for debugging to avoid accidental
 691                  * segfaults when printing the selection tree. */
 692                 if (sel.v.u.p->x)
 693                 {
 694                     fprintf(fp, "(%f, %f, %f)",
 695                             sel.v.u.p->x[0][XX], sel.v.u.p->x[0][YY],
 696                             sel.v.u.p->x[0][ZZ]);
 697                 }
 698                 else
 699                 {
 700                     fprintf(fp, "(null)");
 701                 }
 702                 break;
 703             case GROUP_VALUE:
 704                 fprintf(fp, "%d atoms", sel.v.u.g->isize);
 705                 if (sel.v.u.g->isize < 20)
 706                 {
 707                     if (sel.v.u.g->isize > 0)
 708                     {
 709                         fprintf(fp, ":");
 710                     }
 711                     for (i = 0; i < sel.v.u.g->isize; ++i)
 712                     {
 713                         fprintf(fp, " %d", sel.v.u.g->index[i] + 1);
 714                     }
 715                 }
 716                 break;
 717             default:
 718                 fprintf(fp, "???");
 719                 break;
 720         }
 721         fprintf(fp, "\n");
 722     }
 723
 724     /* Print the subexpressions with one more level of indentation */
 725     gmx::SelectionTreeElementPointer child = sel.child;
 726     while (child)
 727     {
 728         if (!(sel.type == SEL_SUBEXPRREF && child->type == SEL_SUBEXPR))
 729         {
 730             _gmx_selelem_print_tree(fp, *child, bValues, level+1);
 731         }
 732         child = child->next;
 733     }
 734 }
 735
 736 /*!
 737  * \param[in] root Root of the subtree to query.
 738  * \returns true if \p root or any any of its elements require topology
 739  *   information, false otherwise.
 740  */
 741 bool
 742 _gmx_selelem_requires_top(const gmx::SelectionTreeElement &root)
 743 {
 744     if (root.type == SEL_EXPRESSION || root.type == SEL_MODIFIER)
 745     {
 746         if (root.u.expr.method && (root.u.expr.method->flags & SMETH_REQTOP))
 747         {
 748             return true;
 749         }
 750         if (root.u.expr.pc && gmx_ana_poscalc_requires_top(root.u.expr.pc))
 751         {
 752             return true;
 753         }
 754     }
 755     gmx::SelectionTreeElementPointer child = root.child;
 756     while (child)
 757     {
 758         if (_gmx_selelem_requires_top(*child))
 759         {
 760             return true;
 761         }
 762         child = child->next;
 763     }
 764     return false;
 765 }