2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements simple keyword selection methods.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_selection
44 #include "gromacs/legacyheaders/macros.h"
46 #include "gromacs/topology/topology.h"
47 #include "gromacs/selection/position.h"
48 #include "gromacs/selection/selmethod.h"
49 #include "gromacs/utility/exceptions.h"
51 /** Evaluates the \p all selection keyword. */
53 evaluate_all(t_topology *top, t_trxframe *fr, t_pbc *pbc,
54 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
55 /** Evaluates the \p none selection keyword. */
57 evaluate_none(t_topology *top, t_trxframe *fr, t_pbc *pbc,
58 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
59 /** Evaluates the \p atomnr selection keyword. */
61 evaluate_atomnr(t_topology *top, t_trxframe *fr, t_pbc *pbc,
62 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
63 /** Evaluates the \p resnr selection keyword. */
65 evaluate_resnr(t_topology *top, t_trxframe *fr, t_pbc *pbc,
66 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
67 /** Evaluates the \p resindex selection keyword. */
69 evaluate_resindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
70 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
72 * Checks whether molecule information is present in the topology.
74 * \param[in] top Topology structure.
75 * \param npar Not used.
76 * \param param Not used.
77 * \param data Not used.
78 * \returns 0 if molecule info is present in the topology, -1 otherwise.
80 * If molecule information is not found, also prints an error message.
83 check_molecules(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
84 /** Evaluates the \p molindex selection keyword. */
86 evaluate_molindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
87 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
88 /** Evaluates the \p atomname selection keyword. */
90 evaluate_atomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
91 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
92 /** Evaluates the \p pdbatomname selection keyword. */
94 evaluate_pdbatomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
95 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
97 * Checks whether atom types are present in the topology.
99 * \param[in] top Topology structure.
100 * \param npar Not used.
101 * \param param Not used.
102 * \param data Not used.
103 * \returns 0 if atom types are present in the topology, -1 otherwise.
105 * If the atom types are not found, also prints an error message.
108 check_atomtype(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
109 /** Evaluates the \p atomtype selection keyword. */
111 evaluate_atomtype(t_topology *top, t_trxframe *fr, t_pbc *pbc,
112 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
113 /** Evaluates the \p insertcode selection keyword. */
115 evaluate_insertcode(t_topology *top, t_trxframe *fr, t_pbc *pbc,
116 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
117 /** Evaluates the \p chain selection keyword. */
119 evaluate_chain(t_topology *top, t_trxframe *fr, t_pbc *pbc,
120 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
121 /** Evaluates the \p mass selection keyword. */
123 evaluate_mass(t_topology *top, t_trxframe *fr, t_pbc *pbc,
124 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
125 /** Evaluates the \p charge selection keyword. */
127 evaluate_charge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
128 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
130 * Checks whether PDB info is present in the topology.
132 * \param[in] top Topology structure.
133 * \param npar Not used.
134 * \param param Not used.
135 * \param data Not used.
136 * \returns 0 if PDB info is present in the topology, -1 otherwise.
138 * If PDB info is not found, also prints an error message.
141 check_pdbinfo(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
142 /** Evaluates the \p altloc selection keyword. */
144 evaluate_altloc(t_topology *top, t_trxframe *fr, t_pbc *pbc,
145 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
146 /** Evaluates the \p occupancy selection keyword. */
148 evaluate_occupancy(t_topology *top, t_trxframe *fr, t_pbc *pbc,
149 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
150 /** Evaluates the \p betafactor selection keyword. */
152 evaluate_betafactor(t_topology *top, t_trxframe *fr, t_pbc *pbc,
153 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
154 /** Evaluates the \p resname selection keyword. */
156 evaluate_resname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
157 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
159 /** Evaluates the \p x selection keyword. */
161 evaluate_x(t_topology *top, t_trxframe *fr, t_pbc *pbc,
162 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
163 /** Evaluates the \p y selection keyword. */
165 evaluate_y(t_topology *top, t_trxframe *fr, t_pbc *pbc,
166 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
167 /** Evaluates the \p z selection keyword. */
169 evaluate_z(t_topology *top, t_trxframe *fr, t_pbc *pbc,
170 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
172 /** Help text for atom name selection keywords. */
173 static const char *help_atomname[] = {
174 "ATOM NAME SELECTION KEYWORDS[PAR]",
176 "[TT]name[tt] [TT]pdbname[tt] [TT]atomname[tt] [TT]pdbatomname[tt][PAR]",
178 "These keywords select atoms by name. [TT]name[tt] selects atoms using",
179 "the Gromacs atom naming convention.",
180 "For input formats other than PDB, the atom names are matched exactly",
181 "as they appear in the input file. For PDB files, 4 character atom names",
182 "that start with a digit are matched after moving the digit to the end",
183 "(e.g., to match 3HG2 from a PDB file, use [TT]name HG23[tt]).",
184 "[TT]pdbname[tt] can only be used with a PDB input file, and selects",
185 "atoms based on the exact name given in the input file, without the",
186 "transformation described above.[PAR]",
188 "[TT]atomname[tt] and [TT]pdbatomname[tt] are synonyms for the above two",
192 /** Selection method data for \p all selection keyword. */
193 gmx_ana_selmethod_t sm_all = {
194 "all", GROUP_VALUE, 0,
206 /** Selection method data for \p none selection keyword. */
207 gmx_ana_selmethod_t sm_none = {
208 "none", GROUP_VALUE, 0,
220 /** Selection method data for \p atomnr selection keyword. */
221 gmx_ana_selmethod_t sm_atomnr = {
222 "atomnr", INT_VALUE, 0,
234 /** Selection method data for \p resnr selection keyword. */
235 gmx_ana_selmethod_t sm_resnr = {
236 "resnr", INT_VALUE, SMETH_REQTOP,
248 /** Selection method data for \p resindex selection keyword. */
249 gmx_ana_selmethod_t sm_resindex = {
250 "resindex", INT_VALUE, SMETH_REQTOP,
262 /** Selection method data for \p molindex selection keyword. */
263 gmx_ana_selmethod_t sm_molindex = {
264 "molindex", INT_VALUE, SMETH_REQTOP,
276 /** Selection method data for \p atomname selection keyword. */
277 gmx_ana_selmethod_t sm_atomname = {
278 "atomname", STR_VALUE, SMETH_REQTOP,
288 {NULL, asize(help_atomname), help_atomname}
291 /** Selection method data for \p pdbatomname selection keyword. */
292 gmx_ana_selmethod_t sm_pdbatomname = {
293 "pdbatomname", STR_VALUE, SMETH_REQTOP,
301 &evaluate_pdbatomname,
303 {NULL, asize(help_atomname), help_atomname}
306 /** Selection method data for \p atomtype selection keyword. */
307 gmx_ana_selmethod_t sm_atomtype = {
308 "atomtype", STR_VALUE, SMETH_REQTOP,
320 /** Selection method data for \p resname selection keyword. */
321 gmx_ana_selmethod_t sm_resname = {
322 "resname", STR_VALUE, SMETH_REQTOP,
334 /** Selection method data for \p chain selection keyword. */
335 gmx_ana_selmethod_t sm_insertcode = {
336 "insertcode", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
344 &evaluate_insertcode,
348 /** Selection method data for \p chain selection keyword. */
349 gmx_ana_selmethod_t sm_chain = {
350 "chain", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
362 /** Selection method data for \p mass selection keyword. */
363 gmx_ana_selmethod_t sm_mass = {
364 "mass", REAL_VALUE, SMETH_REQTOP,
376 /** Selection method data for \p charge selection keyword. */
377 gmx_ana_selmethod_t sm_charge = {
378 "charge", REAL_VALUE, SMETH_REQTOP,
390 /** Selection method data for \p chain selection keyword. */
391 gmx_ana_selmethod_t sm_altloc = {
392 "altloc", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
404 /** Selection method data for \p occupancy selection keyword. */
405 gmx_ana_selmethod_t sm_occupancy = {
406 "occupancy", REAL_VALUE, SMETH_REQTOP,
418 /** Selection method data for \p betafactor selection keyword. */
419 gmx_ana_selmethod_t sm_betafactor = {
420 "betafactor", REAL_VALUE, SMETH_REQTOP,
428 &evaluate_betafactor,
432 /** Selection method data for \p x selection keyword. */
433 gmx_ana_selmethod_t sm_x = {
434 "x", REAL_VALUE, SMETH_DYNAMIC,
446 /** Selection method data for \p y selection keyword. */
447 gmx_ana_selmethod_t sm_y = {
448 "y", REAL_VALUE, SMETH_DYNAMIC,
460 /** Selection method data for \p z selection keyword. */
461 gmx_ana_selmethod_t sm_z = {
462 "z", REAL_VALUE, SMETH_DYNAMIC,
475 * See sel_updatefunc() for description of the parameters.
476 * \p data is not used.
478 * Copies \p g to \p out->u.g.
481 evaluate_all(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
482 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
484 gmx_ana_index_copy(out->u.g, g, false);
488 * See sel_updatefunc() for description of the parameters.
489 * \p data is not used.
491 * Returns an empty \p out->u.g.
494 evaluate_none(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
495 gmx_ana_index_t * /* g */, gmx_ana_selvalue_t *out, void * /* data */)
501 * See sel_updatefunc() for description of the parameters.
502 * \p data is not used.
504 * Returns the indices for each atom in \p out->u.i.
507 evaluate_atomnr(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
508 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
513 for (i = 0; i < g->isize; ++i)
515 out->u.i[i] = g->index[i] + 1;
520 * See sel_updatefunc() for description of the parameters.
521 * \p data is not used.
523 * Returns the residue numbers for each atom in \p out->u.i.
526 evaluate_resnr(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
527 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
533 for (i = 0; i < g->isize; ++i)
535 resind = top->atoms.atom[g->index[i]].resind;
536 out->u.i[i] = top->atoms.resinfo[resind].nr;
541 * See sel_updatefunc() for description of the parameters.
542 * \p data is not used.
544 * Returns the residue indices for each atom in \p out->u.i.
547 evaluate_resindex(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
548 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
553 for (i = 0; i < g->isize; ++i)
555 out->u.i[i] = top->atoms.atom[g->index[i]].resind + 1;
560 check_molecules(t_topology *top, int /* npar */, gmx_ana_selparam_t * /* param */, void * /* data */)
564 bOk = (top != NULL && top->mols.nr > 0);
567 GMX_THROW(gmx::InconsistentInputError("Molecule information not available in topology"));
572 * See sel_updatefunc() for description of the parameters.
573 * \p data is not used.
575 * Returns the molecule indices for each atom in \p out->u.i.
578 evaluate_molindex(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
579 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
584 for (i = j = 0; i < g->isize; ++i)
586 while (top->mols.index[j + 1] <= g->index[i])
595 * See sel_updatefunc() for description of the parameters.
596 * \p data is not used.
598 * Returns the atom name for each atom in \p out->u.s.
601 evaluate_atomname(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
602 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
607 for (i = 0; i < g->isize; ++i)
609 out->u.s[i] = *top->atoms.atomname[g->index[i]];
614 * See sel_updatefunc() for description of the parameters.
615 * \p data is not used.
617 * Returns the PDB atom name for each atom in \p out->u.s.
620 evaluate_pdbatomname(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
621 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
626 for (i = 0; i < g->isize; ++i)
628 char *s = top->atoms.pdbinfo[g->index[i]].atomnm;
629 while (std::isspace(*s))
638 check_atomtype(t_topology *top, int /* npar */, gmx_ana_selparam_t * /* param */, void * /* data */)
642 bOk = (top != NULL && top->atoms.atomtype != NULL);
645 GMX_THROW(gmx::InconsistentInputError("Atom types not available in topology"));
650 * See sel_updatefunc() for description of the parameters.
651 * \p data is not used.
653 * Returns the atom type for each atom in \p out->u.s.
654 * Segfaults if atom types are not found in the topology.
657 evaluate_atomtype(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
658 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
663 for (i = 0; i < g->isize; ++i)
665 out->u.s[i] = *top->atoms.atomtype[g->index[i]];
670 * See sel_updatefunc() for description of the parameters.
671 * \p data is not used.
673 * Returns the residue name for each atom in \p out->u.s.
676 evaluate_resname(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
677 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
683 for (i = 0; i < g->isize; ++i)
685 resind = top->atoms.atom[g->index[i]].resind;
686 out->u.s[i] = *top->atoms.resinfo[resind].name;
691 * See sel_updatefunc() for description of the parameters.
692 * \p data is not used.
694 * Returns the insertion code for each atom in \p out->u.s.
697 evaluate_insertcode(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
698 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
704 for (i = 0; i < g->isize; ++i)
706 resind = top->atoms.atom[g->index[i]].resind;
707 out->u.s[i][0] = top->atoms.resinfo[resind].ic;
712 * See sel_updatefunc() for description of the parameters.
713 * \p data is not used.
715 * Returns the chain for each atom in \p out->u.s.
718 evaluate_chain(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
719 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
725 for (i = 0; i < g->isize; ++i)
727 resind = top->atoms.atom[g->index[i]].resind;
728 out->u.s[i][0] = top->atoms.resinfo[resind].chainid;
733 * See sel_updatefunc() for description of the parameters.
734 * \p data is not used.
736 * Returns the mass for each atom in \p out->u.r.
739 evaluate_mass(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
740 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
745 for (i = 0; i < g->isize; ++i)
747 out->u.r[i] = top->atoms.atom[g->index[i]].m;
752 * See sel_updatefunc() for description of the parameters.
753 * \p data is not used.
755 * Returns the charge for each atom in \p out->u.r.
758 evaluate_charge(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
759 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
764 for (i = 0; i < g->isize; ++i)
766 out->u.r[i] = top->atoms.atom[g->index[i]].q;
771 check_pdbinfo(t_topology *top, int /* npar */, gmx_ana_selparam_t * /* param */, void * /* data */)
775 bOk = (top != NULL && top->atoms.pdbinfo != NULL);
778 GMX_THROW(gmx::InconsistentInputError("PDB info not available in topology"));
783 * See sel_updatefunc() for description of the parameters.
784 * \p data is not used.
786 * Returns the alternate location identifier for each atom in \p out->u.s.
789 evaluate_altloc(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
790 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
795 for (i = 0; i < g->isize; ++i)
797 out->u.s[i][0] = top->atoms.pdbinfo[g->index[i]].altloc;
802 * See sel_updatefunc() for description of the parameters.
803 * \p data is not used.
805 * Returns the occupancy numbers for each atom in \p out->u.r.
806 * Segfaults if PDB info is not found in the topology.
809 evaluate_occupancy(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
810 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
815 for (i = 0; i < g->isize; ++i)
817 out->u.r[i] = top->atoms.pdbinfo[g->index[i]].occup;
822 * See sel_updatefunc() for description of the parameters.
823 * \p data is not used.
825 * Returns the B-factors for each atom in \p out->u.r.
826 * Segfaults if PDB info is not found in the topology.
829 evaluate_betafactor(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
830 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
835 for (i = 0; i < g->isize; ++i)
837 out->u.r[i] = top->atoms.pdbinfo[g->index[i]].bfac;
842 * Internal utility function for position keyword evaluation.
844 * \param[out] out Output array.
845 * \param[in] pos Position data to use instead of atomic coordinates
847 * \param[in] d Coordinate index to evaluate (\p XX, \p YY or \p ZZ).
849 * This function is used internally by evaluate_x(), evaluate_y() and
850 * evaluate_z() to do the actual evaluation.
853 evaluate_coord(real out[], gmx_ana_pos_t *pos, int d)
855 for (int b = 0; b < pos->count(); ++b)
857 const real v = pos->x[b][d];
858 for (int i = pos->m.mapb.index[b]; i < pos->m.mapb.index[b+1]; ++i)
863 // TODO: Make this more efficient by directly extracting the coordinates
864 // from the frame coordinates for atomic positions instead of going through
865 // a position calculation.
869 * See sel_updatefunc_pos() for description of the parameters.
870 * \p data is not used.
872 * Returns the \p x coordinate for each atom in \p out->u.r.
875 evaluate_x(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
876 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
878 out->nr = pos->m.mapb.nra;
879 evaluate_coord(out->u.r, pos, XX);
883 * See sel_updatefunc() for description of the parameters.
884 * \p data is not used.
886 * Returns the \p y coordinate for each atom in \p out->u.r.
889 evaluate_y(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
890 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
892 out->nr = pos->m.mapb.nra;
893 evaluate_coord(out->u.r, pos, YY);
897 * See sel_updatefunc() for description of the parameters.
898 * \p data is not used.
900 * Returns the \p z coordinate for each atom in \p out->u.r.
903 evaluate_z(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
904 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
906 out->nr = pos->m.mapb.nra;
907 evaluate_coord(out->u.r, pos, ZZ);