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/selection/position.h"
47 #include "gromacs/selection/selmethod.h"
48 #include "gromacs/utility/exceptions.h"
50 /** Evaluates the \p all selection keyword. */
52 evaluate_all(t_topology *top, t_trxframe *fr, t_pbc *pbc,
53 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
54 /** Evaluates the \p none selection keyword. */
56 evaluate_none(t_topology *top, t_trxframe *fr, t_pbc *pbc,
57 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
58 /** Evaluates the \p atomnr selection keyword. */
60 evaluate_atomnr(t_topology *top, t_trxframe *fr, t_pbc *pbc,
61 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
62 /** Evaluates the \p resnr selection keyword. */
64 evaluate_resnr(t_topology *top, t_trxframe *fr, t_pbc *pbc,
65 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
66 /** Evaluates the \p resindex selection keyword. */
68 evaluate_resindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
69 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
71 * Checks whether molecule information is present in the topology.
73 * \param[in] top Topology structure.
74 * \param npar Not used.
75 * \param param Not used.
76 * \param data Not used.
77 * \returns 0 if molecule info is present in the topology, -1 otherwise.
79 * If molecule information is not found, also prints an error message.
82 check_molecules(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
83 /** Evaluates the \p molindex selection keyword. */
85 evaluate_molindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
86 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
87 /** Evaluates the \p atomname selection keyword. */
89 evaluate_atomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
90 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
91 /** Evaluates the \p pdbatomname selection keyword. */
93 evaluate_pdbatomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
94 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
96 * Checks whether atom types are present in the topology.
98 * \param[in] top Topology structure.
99 * \param npar Not used.
100 * \param param Not used.
101 * \param data Not used.
102 * \returns 0 if atom types are present in the topology, -1 otherwise.
104 * If the atom types are not found, also prints an error message.
107 check_atomtype(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
108 /** Evaluates the \p atomtype selection keyword. */
110 evaluate_atomtype(t_topology *top, t_trxframe *fr, t_pbc *pbc,
111 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
112 /** Evaluates the \p insertcode selection keyword. */
114 evaluate_insertcode(t_topology *top, t_trxframe *fr, t_pbc *pbc,
115 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
116 /** Evaluates the \p chain selection keyword. */
118 evaluate_chain(t_topology *top, t_trxframe *fr, t_pbc *pbc,
119 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
120 /** Evaluates the \p mass selection keyword. */
122 evaluate_mass(t_topology *top, t_trxframe *fr, t_pbc *pbc,
123 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
124 /** Evaluates the \p charge selection keyword. */
126 evaluate_charge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
127 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
129 * Checks whether PDB info is present in the topology.
131 * \param[in] top Topology structure.
132 * \param npar Not used.
133 * \param param Not used.
134 * \param data Not used.
135 * \returns 0 if PDB info is present in the topology, -1 otherwise.
137 * If PDB info is not found, also prints an error message.
140 check_pdbinfo(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
141 /** Evaluates the \p altloc selection keyword. */
143 evaluate_altloc(t_topology *top, t_trxframe *fr, t_pbc *pbc,
144 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
145 /** Evaluates the \p occupancy selection keyword. */
147 evaluate_occupancy(t_topology *top, t_trxframe *fr, t_pbc *pbc,
148 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
149 /** Evaluates the \p betafactor selection keyword. */
151 evaluate_betafactor(t_topology *top, t_trxframe *fr, t_pbc *pbc,
152 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
153 /** Evaluates the \p resname selection keyword. */
155 evaluate_resname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
156 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
158 /** Evaluates the \p x selection keyword. */
160 evaluate_x(t_topology *top, t_trxframe *fr, t_pbc *pbc,
161 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
162 /** Evaluates the \p y selection keyword. */
164 evaluate_y(t_topology *top, t_trxframe *fr, t_pbc *pbc,
165 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
166 /** Evaluates the \p z selection keyword. */
168 evaluate_z(t_topology *top, t_trxframe *fr, t_pbc *pbc,
169 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
171 /** Help text for atom name selection keywords. */
172 static const char *help_atomname[] = {
173 "ATOM NAME SELECTION KEYWORDS[PAR]",
175 "[TT]name[tt] [TT]pdbname[tt] [TT]atomname[tt] [TT]pdbatomname[tt][PAR]",
177 "These keywords select atoms by name. [TT]name[tt] selects atoms using",
178 "the Gromacs atom naming convention.",
179 "For input formats other than PDB, the atom names are matched exactly",
180 "as they appear in the input file. For PDB files, 4 character atom names",
181 "that start with a digit are matched after moving the digit to the end",
182 "(e.g., to match 3HG2 from a PDB file, use [TT]name HG23[tt]).",
183 "[TT]pdbname[tt] can only be used with a PDB input file, and selects",
184 "atoms based on the exact name given in the input file, without the",
185 "transformation described above.[PAR]",
187 "[TT]atomname[tt] and [TT]pdbatomname[tt] are synonyms for the above two",
191 /** Selection method data for \p all selection keyword. */
192 gmx_ana_selmethod_t sm_all = {
193 "all", GROUP_VALUE, 0,
205 /** Selection method data for \p none selection keyword. */
206 gmx_ana_selmethod_t sm_none = {
207 "none", GROUP_VALUE, 0,
219 /** Selection method data for \p atomnr selection keyword. */
220 gmx_ana_selmethod_t sm_atomnr = {
221 "atomnr", INT_VALUE, 0,
233 /** Selection method data for \p resnr selection keyword. */
234 gmx_ana_selmethod_t sm_resnr = {
235 "resnr", INT_VALUE, SMETH_REQTOP,
247 /** Selection method data for \p resindex selection keyword. */
248 gmx_ana_selmethod_t sm_resindex = {
249 "resindex", INT_VALUE, SMETH_REQTOP,
261 /** Selection method data for \p molindex selection keyword. */
262 gmx_ana_selmethod_t sm_molindex = {
263 "molindex", INT_VALUE, SMETH_REQTOP,
275 /** Selection method data for \p atomname selection keyword. */
276 gmx_ana_selmethod_t sm_atomname = {
277 "atomname", STR_VALUE, SMETH_REQTOP,
287 {NULL, asize(help_atomname), help_atomname}
290 /** Selection method data for \p pdbatomname selection keyword. */
291 gmx_ana_selmethod_t sm_pdbatomname = {
292 "pdbatomname", STR_VALUE, SMETH_REQTOP,
300 &evaluate_pdbatomname,
302 {NULL, asize(help_atomname), help_atomname}
305 /** Selection method data for \p atomtype selection keyword. */
306 gmx_ana_selmethod_t sm_atomtype = {
307 "atomtype", STR_VALUE, SMETH_REQTOP,
319 /** Selection method data for \p resname selection keyword. */
320 gmx_ana_selmethod_t sm_resname = {
321 "resname", STR_VALUE, SMETH_REQTOP,
333 /** Selection method data for \p chain selection keyword. */
334 gmx_ana_selmethod_t sm_insertcode = {
335 "insertcode", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
343 &evaluate_insertcode,
347 /** Selection method data for \p chain selection keyword. */
348 gmx_ana_selmethod_t sm_chain = {
349 "chain", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
361 /** Selection method data for \p mass selection keyword. */
362 gmx_ana_selmethod_t sm_mass = {
363 "mass", REAL_VALUE, SMETH_REQTOP,
375 /** Selection method data for \p charge selection keyword. */
376 gmx_ana_selmethod_t sm_charge = {
377 "charge", REAL_VALUE, SMETH_REQTOP,
389 /** Selection method data for \p chain selection keyword. */
390 gmx_ana_selmethod_t sm_altloc = {
391 "altloc", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
403 /** Selection method data for \p occupancy selection keyword. */
404 gmx_ana_selmethod_t sm_occupancy = {
405 "occupancy", REAL_VALUE, SMETH_REQTOP,
417 /** Selection method data for \p betafactor selection keyword. */
418 gmx_ana_selmethod_t sm_betafactor = {
419 "betafactor", REAL_VALUE, SMETH_REQTOP,
427 &evaluate_betafactor,
431 /** Selection method data for \p x selection keyword. */
432 gmx_ana_selmethod_t sm_x = {
433 "x", REAL_VALUE, SMETH_DYNAMIC,
445 /** Selection method data for \p y selection keyword. */
446 gmx_ana_selmethod_t sm_y = {
447 "y", REAL_VALUE, SMETH_DYNAMIC,
459 /** Selection method data for \p z selection keyword. */
460 gmx_ana_selmethod_t sm_z = {
461 "z", REAL_VALUE, SMETH_DYNAMIC,
474 * See sel_updatefunc() for description of the parameters.
475 * \p data is not used.
477 * Copies \p g to \p out->u.g.
480 evaluate_all(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
481 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
483 gmx_ana_index_copy(out->u.g, g, false);
487 * See sel_updatefunc() for description of the parameters.
488 * \p data is not used.
490 * Returns an empty \p out->u.g.
493 evaluate_none(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
494 gmx_ana_index_t * /* g */, gmx_ana_selvalue_t *out, void * /* data */)
500 * See sel_updatefunc() for description of the parameters.
501 * \p data is not used.
503 * Returns the indices for each atom in \p out->u.i.
506 evaluate_atomnr(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
507 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
512 for (i = 0; i < g->isize; ++i)
514 out->u.i[i] = g->index[i] + 1;
519 * See sel_updatefunc() for description of the parameters.
520 * \p data is not used.
522 * Returns the residue numbers for each atom in \p out->u.i.
525 evaluate_resnr(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
526 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
532 for (i = 0; i < g->isize; ++i)
534 resind = top->atoms.atom[g->index[i]].resind;
535 out->u.i[i] = top->atoms.resinfo[resind].nr;
540 * See sel_updatefunc() for description of the parameters.
541 * \p data is not used.
543 * Returns the residue indices for each atom in \p out->u.i.
546 evaluate_resindex(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
547 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
552 for (i = 0; i < g->isize; ++i)
554 out->u.i[i] = top->atoms.atom[g->index[i]].resind + 1;
559 check_molecules(t_topology *top, int /* npar */, gmx_ana_selparam_t * /* param */, void * /* data */)
563 bOk = (top != NULL && top->mols.nr > 0);
566 GMX_THROW(gmx::InconsistentInputError("Molecule information not available in topology"));
571 * See sel_updatefunc() for description of the parameters.
572 * \p data is not used.
574 * Returns the molecule indices for each atom in \p out->u.i.
577 evaluate_molindex(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
578 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
583 for (i = j = 0; i < g->isize; ++i)
585 while (top->mols.index[j + 1] <= g->index[i])
594 * See sel_updatefunc() for description of the parameters.
595 * \p data is not used.
597 * Returns the atom name for each atom in \p out->u.s.
600 evaluate_atomname(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
601 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
606 for (i = 0; i < g->isize; ++i)
608 out->u.s[i] = *top->atoms.atomname[g->index[i]];
613 * See sel_updatefunc() for description of the parameters.
614 * \p data is not used.
616 * Returns the PDB atom name for each atom in \p out->u.s.
619 evaluate_pdbatomname(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
620 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
625 for (i = 0; i < g->isize; ++i)
627 char *s = top->atoms.pdbinfo[g->index[i]].atomnm;
628 while (std::isspace(*s))
637 check_atomtype(t_topology *top, int /* npar */, gmx_ana_selparam_t * /* param */, void * /* data */)
641 bOk = (top != NULL && top->atoms.atomtype != NULL);
644 GMX_THROW(gmx::InconsistentInputError("Atom types not available in topology"));
649 * See sel_updatefunc() for description of the parameters.
650 * \p data is not used.
652 * Returns the atom type for each atom in \p out->u.s.
653 * Segfaults if atom types are not found in the topology.
656 evaluate_atomtype(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
657 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
662 for (i = 0; i < g->isize; ++i)
664 out->u.s[i] = *top->atoms.atomtype[g->index[i]];
669 * See sel_updatefunc() for description of the parameters.
670 * \p data is not used.
672 * Returns the residue name for each atom in \p out->u.s.
675 evaluate_resname(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
676 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
682 for (i = 0; i < g->isize; ++i)
684 resind = top->atoms.atom[g->index[i]].resind;
685 out->u.s[i] = *top->atoms.resinfo[resind].name;
690 * See sel_updatefunc() for description of the parameters.
691 * \p data is not used.
693 * Returns the insertion code for each atom in \p out->u.s.
696 evaluate_insertcode(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
697 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
703 for (i = 0; i < g->isize; ++i)
705 resind = top->atoms.atom[g->index[i]].resind;
706 out->u.s[i][0] = top->atoms.resinfo[resind].ic;
711 * See sel_updatefunc() for description of the parameters.
712 * \p data is not used.
714 * Returns the chain for each atom in \p out->u.s.
717 evaluate_chain(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
718 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
724 for (i = 0; i < g->isize; ++i)
726 resind = top->atoms.atom[g->index[i]].resind;
727 out->u.s[i][0] = top->atoms.resinfo[resind].chainid;
732 * See sel_updatefunc() for description of the parameters.
733 * \p data is not used.
735 * Returns the mass for each atom in \p out->u.r.
738 evaluate_mass(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
739 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
744 for (i = 0; i < g->isize; ++i)
746 out->u.r[i] = top->atoms.atom[g->index[i]].m;
751 * See sel_updatefunc() for description of the parameters.
752 * \p data is not used.
754 * Returns the charge for each atom in \p out->u.r.
757 evaluate_charge(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
758 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
763 for (i = 0; i < g->isize; ++i)
765 out->u.r[i] = top->atoms.atom[g->index[i]].q;
770 check_pdbinfo(t_topology *top, int /* npar */, gmx_ana_selparam_t * /* param */, void * /* data */)
774 bOk = (top != NULL && top->atoms.pdbinfo != NULL);
777 GMX_THROW(gmx::InconsistentInputError("PDB info not available in topology"));
782 * See sel_updatefunc() for description of the parameters.
783 * \p data is not used.
785 * Returns the alternate location identifier for each atom in \p out->u.s.
788 evaluate_altloc(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
789 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
794 for (i = 0; i < g->isize; ++i)
796 out->u.s[i][0] = top->atoms.pdbinfo[g->index[i]].altloc;
801 * See sel_updatefunc() for description of the parameters.
802 * \p data is not used.
804 * Returns the occupancy numbers for each atom in \p out->u.r.
805 * Segfaults if PDB info is not found in the topology.
808 evaluate_occupancy(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
809 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
814 for (i = 0; i < g->isize; ++i)
816 out->u.r[i] = top->atoms.pdbinfo[g->index[i]].occup;
821 * See sel_updatefunc() for description of the parameters.
822 * \p data is not used.
824 * Returns the B-factors for each atom in \p out->u.r.
825 * Segfaults if PDB info is not found in the topology.
828 evaluate_betafactor(t_topology *top, t_trxframe * /* fr */, t_pbc * /* pbc */,
829 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */)
834 for (i = 0; i < g->isize; ++i)
836 out->u.r[i] = top->atoms.pdbinfo[g->index[i]].bfac;
841 * Internal utility function for position keyword evaluation.
843 * \param[out] out Output array.
844 * \param[in] pos Position data to use instead of atomic coordinates
846 * \param[in] d Coordinate index to evaluate (\p XX, \p YY or \p ZZ).
848 * This function is used internally by evaluate_x(), evaluate_y() and
849 * evaluate_z() to do the actual evaluation.
852 evaluate_coord(real out[], gmx_ana_pos_t *pos, int d)
854 for (int b = 0; b < pos->count(); ++b)
856 const real v = pos->x[b][d];
857 for (int i = pos->m.mapb.index[b]; i < pos->m.mapb.index[b+1]; ++i)
862 // TODO: Make this more efficient by directly extracting the coordinates
863 // from the frame coordinates for atomic positions instead of going through
864 // a position calculation.
868 * See sel_updatefunc_pos() for description of the parameters.
869 * \p data is not used.
871 * Returns the \p x coordinate for each atom in \p out->u.r.
874 evaluate_x(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
875 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
877 out->nr = pos->m.mapb.nra;
878 evaluate_coord(out->u.r, pos, XX);
882 * See sel_updatefunc() for description of the parameters.
883 * \p data is not used.
885 * Returns the \p y coordinate for each atom in \p out->u.r.
888 evaluate_y(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
889 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
891 out->nr = pos->m.mapb.nra;
892 evaluate_coord(out->u.r, pos, YY);
896 * See sel_updatefunc() for description of the parameters.
897 * \p data is not used.
899 * Returns the \p z coordinate for each atom in \p out->u.r.
902 evaluate_z(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
903 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
905 out->nr = pos->m.mapb.nra;
906 evaluate_coord(out->u.r, pos, ZZ);