2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * 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);
70 /** Checks whether molecule information is present in the topology. */
72 check_molecules(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
73 /** Evaluates the \p molindex selection keyword. */
75 evaluate_molindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
76 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
77 /** Evaluates the \p atomname selection keyword. */
79 evaluate_atomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
80 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
81 /** Evaluates the \p pdbatomname selection keyword. */
83 evaluate_pdbatomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
84 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
85 /** Checks whether atom types are present in the topology. */
87 check_atomtype(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
88 /** Evaluates the \p atomtype selection keyword. */
90 evaluate_atomtype(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 insertcode selection keyword. */
94 evaluate_insertcode(t_topology *top, t_trxframe *fr, t_pbc *pbc,
95 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
96 /** Evaluates the \p chain selection keyword. */
98 evaluate_chain(t_topology *top, t_trxframe *fr, t_pbc *pbc,
99 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
100 /** Evaluates the \p mass selection keyword. */
102 evaluate_mass(t_topology *top, t_trxframe *fr, t_pbc *pbc,
103 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
104 /** Evaluates the \p charge selection keyword. */
106 evaluate_charge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
107 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
108 /** Checks whether PDB info is present in the topology. */
110 check_pdbinfo(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
111 /** Evaluates the \p altloc selection keyword. */
113 evaluate_altloc(t_topology *top, t_trxframe *fr, t_pbc *pbc,
114 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
115 /** Evaluates the \p occupancy selection keyword. */
117 evaluate_occupancy(t_topology *top, t_trxframe *fr, t_pbc *pbc,
118 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
119 /** Evaluates the \p betafactor selection keyword. */
121 evaluate_betafactor(t_topology *top, t_trxframe *fr, t_pbc *pbc,
122 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
123 /** Evaluates the \p resname selection keyword. */
125 evaluate_resname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
126 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
128 /** Evaluates the \p x selection keyword. */
130 evaluate_x(t_topology *top, t_trxframe *fr, t_pbc *pbc,
131 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
132 /** Evaluates the \p y selection keyword. */
134 evaluate_y(t_topology *top, t_trxframe *fr, t_pbc *pbc,
135 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
136 /** Evaluates the \p z selection keyword. */
138 evaluate_z(t_topology *top, t_trxframe *fr, t_pbc *pbc,
139 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void *data);
141 /** Help text for atom name selection keywords. */
142 static const char *help_atomname[] = {
143 "ATOM NAME SELECTION KEYWORDS[PAR]",
145 "[TT]name[tt] [TT]pdbname[tt] [TT]atomname[tt] [TT]pdbatomname[tt][PAR]",
147 "These keywords select atoms by name. [TT]name[tt] selects atoms using",
148 "the Gromacs atom naming convention.",
149 "For input formats other than PDB, the atom names are matched exactly",
150 "as they appear in the input file. For PDB files, 4 character atom names",
151 "that start with a digit are matched after moving the digit to the end",
152 "(e.g., to match 3HG2 from a PDB file, use [TT]name HG23[tt]).",
153 "[TT]pdbname[tt] can only be used with a PDB input file, and selects",
154 "atoms based on the exact name given in the input file, without the",
155 "transformation described above.[PAR]",
157 "[TT]atomname[tt] and [TT]pdbatomname[tt] are synonyms for the above two",
161 /** \internal Selection method data for \p all selection keyword. */
162 gmx_ana_selmethod_t sm_all = {
163 "all", GROUP_VALUE, 0,
175 /** \internal Selection method data for \p none selection keyword. */
176 gmx_ana_selmethod_t sm_none = {
177 "none", GROUP_VALUE, 0,
189 /** \internal Selection method data for \p atomnr selection keyword. */
190 gmx_ana_selmethod_t sm_atomnr = {
191 "atomnr", INT_VALUE, 0,
203 /** \internal Selection method data for \p resnr selection keyword. */
204 gmx_ana_selmethod_t sm_resnr = {
205 "resnr", INT_VALUE, SMETH_REQTOP,
217 /** \internal Selection method data for \p resindex selection keyword. */
218 gmx_ana_selmethod_t sm_resindex = {
219 "resindex", INT_VALUE, SMETH_REQTOP,
231 /** \internal Selection method data for \p molindex selection keyword. */
232 gmx_ana_selmethod_t sm_molindex = {
233 "molindex", INT_VALUE, SMETH_REQTOP,
245 /** \internal Selection method data for \p atomname selection keyword. */
246 gmx_ana_selmethod_t sm_atomname = {
247 "atomname", STR_VALUE, SMETH_REQTOP,
257 {NULL, asize(help_atomname), help_atomname}
260 /** \internal Selection method data for \p pdbatomname selection keyword. */
261 gmx_ana_selmethod_t sm_pdbatomname = {
262 "pdbatomname", STR_VALUE, SMETH_REQTOP,
270 &evaluate_pdbatomname,
272 {NULL, asize(help_atomname), help_atomname}
275 /** \internal Selection method data for \p atomtype selection keyword. */
276 gmx_ana_selmethod_t sm_atomtype = {
277 "atomtype", STR_VALUE, SMETH_REQTOP,
289 /** \internal Selection method data for \p resname selection keyword. */
290 gmx_ana_selmethod_t sm_resname = {
291 "resname", STR_VALUE, SMETH_REQTOP,
303 /** \internal Selection method data for \p chain selection keyword. */
304 gmx_ana_selmethod_t sm_insertcode = {
305 "insertcode", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
313 &evaluate_insertcode,
317 /** \internal Selection method data for \p chain selection keyword. */
318 gmx_ana_selmethod_t sm_chain = {
319 "chain", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
331 /** \internal Selection method data for \p mass selection keyword. */
332 gmx_ana_selmethod_t sm_mass = {
333 "mass", REAL_VALUE, SMETH_REQTOP,
345 /** \internal Selection method data for \p charge selection keyword. */
346 gmx_ana_selmethod_t sm_charge = {
347 "charge", REAL_VALUE, SMETH_REQTOP,
359 /** \internal Selection method data for \p chain selection keyword. */
360 gmx_ana_selmethod_t sm_altloc = {
361 "altloc", STR_VALUE, SMETH_REQTOP | SMETH_CHARVAL,
373 /** \internal Selection method data for \p occupancy selection keyword. */
374 gmx_ana_selmethod_t sm_occupancy = {
375 "occupancy", REAL_VALUE, SMETH_REQTOP,
387 /** \internal Selection method data for \p betafactor selection keyword. */
388 gmx_ana_selmethod_t sm_betafactor = {
389 "betafactor", REAL_VALUE, SMETH_REQTOP,
397 &evaluate_betafactor,
401 /** \internal Selection method data for \p x selection keyword. */
402 gmx_ana_selmethod_t sm_x = {
403 "x", REAL_VALUE, SMETH_DYNAMIC,
415 /** \internal Selection method data for \p y selection keyword. */
416 gmx_ana_selmethod_t sm_y = {
417 "y", REAL_VALUE, SMETH_DYNAMIC,
429 /** \internal Selection method data for \p z selection keyword. */
430 gmx_ana_selmethod_t sm_z = {
431 "z", REAL_VALUE, SMETH_DYNAMIC,
444 * See sel_updatefunc() for description of the parameters.
445 * \p data is not used.
447 * Copies \p g to \p out->u.g.
450 evaluate_all(t_topology *top, t_trxframe *fr, t_pbc *pbc,
451 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
453 gmx_ana_index_copy(out->u.g, g, false);
457 * See sel_updatefunc() for description of the parameters.
458 * \p data is not used.
460 * Returns an empty \p out->u.g.
463 evaluate_none(t_topology *top, t_trxframe *fr, t_pbc *pbc,
464 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
470 * See sel_updatefunc() for description of the parameters.
471 * \p data is not used.
473 * Returns the indices for each atom in \p out->u.i.
476 evaluate_atomnr(t_topology *top, t_trxframe *fr, t_pbc *pbc,
477 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
482 for (i = 0; i < g->isize; ++i)
484 out->u.i[i] = g->index[i] + 1;
489 * See sel_updatefunc() for description of the parameters.
490 * \p data is not used.
492 * Returns the residue numbers for each atom in \p out->u.i.
495 evaluate_resnr(t_topology *top, t_trxframe *fr, t_pbc *pbc,
496 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
502 for (i = 0; i < g->isize; ++i)
504 resind = top->atoms.atom[g->index[i]].resind;
505 out->u.i[i] = top->atoms.resinfo[resind].nr;
510 * See sel_updatefunc() for description of the parameters.
511 * \p data is not used.
513 * Returns the residue indices for each atom in \p out->u.i.
516 evaluate_resindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
517 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
522 for (i = 0; i < g->isize; ++i)
524 out->u.i[i] = top->atoms.atom[g->index[i]].resind + 1;
529 * \param[in] top Topology structure.
530 * \param npar Not used.
531 * \param param Not used.
532 * \param data Not used.
533 * \returns 0 if molecule info is present in the topology, -1 otherwise.
535 * If molecule information is not found, also prints an error message.
538 check_molecules(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
542 bOk = (top != NULL && top->mols.nr > 0);
545 GMX_THROW(gmx::InconsistentInputError("Molecule information not available in topology"));
550 * See sel_updatefunc() for description of the parameters.
551 * \p data is not used.
553 * Returns the molecule indices for each atom in \p out->u.i.
556 evaluate_molindex(t_topology *top, t_trxframe *fr, t_pbc *pbc,
557 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
562 for (i = j = 0; i < g->isize; ++i)
564 while (top->mols.index[j + 1] <= g->index[i])
573 * See sel_updatefunc() for description of the parameters.
574 * \p data is not used.
576 * Returns the atom name for each atom in \p out->u.s.
579 evaluate_atomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
580 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
585 for (i = 0; i < g->isize; ++i)
587 out->u.s[i] = *top->atoms.atomname[g->index[i]];
592 * See sel_updatefunc() for description of the parameters.
593 * \p data is not used.
595 * Returns the PDB atom name for each atom in \p out->u.s.
598 evaluate_pdbatomname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
599 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
604 for (i = 0; i < g->isize; ++i)
606 char *s = top->atoms.pdbinfo[g->index[i]].atomnm;
607 while (std::isspace(*s))
616 * \param[in] top Topology structure.
617 * \param npar Not used.
618 * \param param Not used.
619 * \param data Not used.
620 * \returns 0 if atom types are present in the topology, -1 otherwise.
622 * If the atom types are not found, also prints an error message.
625 check_atomtype(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
629 bOk = (top != NULL && top->atoms.atomtype != NULL);
632 GMX_THROW(gmx::InconsistentInputError("Atom types not available in topology"));
637 * See sel_updatefunc() for description of the parameters.
638 * \p data is not used.
640 * Returns the atom type for each atom in \p out->u.s.
641 * Segfaults if atom types are not found in the topology.
644 evaluate_atomtype(t_topology *top, t_trxframe *fr, t_pbc *pbc,
645 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
650 for (i = 0; i < g->isize; ++i)
652 out->u.s[i] = *top->atoms.atomtype[g->index[i]];
657 * See sel_updatefunc() for description of the parameters.
658 * \p data is not used.
660 * Returns the residue name for each atom in \p out->u.s.
663 evaluate_resname(t_topology *top, t_trxframe *fr, t_pbc *pbc,
664 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
670 for (i = 0; i < g->isize; ++i)
672 resind = top->atoms.atom[g->index[i]].resind;
673 out->u.s[i] = *top->atoms.resinfo[resind].name;
678 * See sel_updatefunc() for description of the parameters.
679 * \p data is not used.
681 * Returns the insertion code for each atom in \p out->u.s.
684 evaluate_insertcode(t_topology *top, t_trxframe *fr, t_pbc *pbc,
685 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
691 for (i = 0; i < g->isize; ++i)
693 resind = top->atoms.atom[g->index[i]].resind;
694 out->u.s[i][0] = top->atoms.resinfo[resind].ic;
699 * See sel_updatefunc() for description of the parameters.
700 * \p data is not used.
702 * Returns the chain for each atom in \p out->u.s.
705 evaluate_chain(t_topology *top, t_trxframe *fr, t_pbc *pbc,
706 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
712 for (i = 0; i < g->isize; ++i)
714 resind = top->atoms.atom[g->index[i]].resind;
715 out->u.s[i][0] = top->atoms.resinfo[resind].chainid;
720 * See sel_updatefunc() for description of the parameters.
721 * \p data is not used.
723 * Returns the mass for each atom in \p out->u.r.
726 evaluate_mass(t_topology *top, t_trxframe *fr, t_pbc *pbc,
727 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
732 for (i = 0; i < g->isize; ++i)
734 out->u.r[i] = top->atoms.atom[g->index[i]].m;
739 * See sel_updatefunc() for description of the parameters.
740 * \p data is not used.
742 * Returns the charge for each atom in \p out->u.r.
745 evaluate_charge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
746 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
751 for (i = 0; i < g->isize; ++i)
753 out->u.r[i] = top->atoms.atom[g->index[i]].q;
758 * \param[in] top Topology structure.
759 * \param npar Not used.
760 * \param param Not used.
761 * \param data Not used.
762 * \returns 0 if PDB info is present in the topology, -1 otherwise.
764 * If PDB info is not found, also prints an error message.
767 check_pdbinfo(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
771 bOk = (top != NULL && top->atoms.pdbinfo != NULL);
774 GMX_THROW(gmx::InconsistentInputError("PDB info not available in topology"));
779 * See sel_updatefunc() for description of the parameters.
780 * \p data is not used.
782 * Returns the alternate location identifier for each atom in \p out->u.s.
785 evaluate_altloc(t_topology *top, t_trxframe *fr, t_pbc *pbc,
786 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
791 for (i = 0; i < g->isize; ++i)
793 out->u.s[i][0] = top->atoms.pdbinfo[g->index[i]].altloc;
798 * See sel_updatefunc() for description of the parameters.
799 * \p data is not used.
801 * Returns the occupancy numbers for each atom in \p out->u.r.
802 * Segfaults if PDB info is not found in the topology.
805 evaluate_occupancy(t_topology *top, t_trxframe *fr, t_pbc *pbc,
806 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
811 for (i = 0; i < g->isize; ++i)
813 out->u.r[i] = top->atoms.pdbinfo[g->index[i]].occup;
818 * See sel_updatefunc() for description of the parameters.
819 * \p data is not used.
821 * Returns the B-factors for each atom in \p out->u.r.
822 * Segfaults if PDB info is not found in the topology.
825 evaluate_betafactor(t_topology *top, t_trxframe *fr, t_pbc *pbc,
826 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
831 for (i = 0; i < g->isize; ++i)
833 out->u.r[i] = top->atoms.pdbinfo[g->index[i]].bfac;
838 * Internal utility function for position keyword evaluation.
840 * \param[out] out Output array.
841 * \param[in] pos Position data to use instead of atomic coordinates
843 * \param[in] d Coordinate index to evaluate (\p XX, \p YY or \p ZZ).
845 * This function is used internally by evaluate_x(), evaluate_y() and
846 * evaluate_z() to do the actual evaluation.
849 evaluate_coord(real out[], gmx_ana_pos_t *pos, int d)
851 for (int b = 0; b < pos->nr; ++b)
853 const real v = pos->x[b][d];
854 for (int i = pos->m.mapb.index[b]; i < pos->m.mapb.index[b+1]; ++i)
859 // TODO: Make this more efficient by directly extracting the coordinates
860 // from the frame coordinates for atomic positions instead of going through
861 // a position calculation.
865 * See sel_updatefunc_pos() for description of the parameters.
866 * \p data is not used.
868 * Returns the \p x coordinate for each atom in \p out->u.r.
871 evaluate_x(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
872 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
874 out->nr = pos->m.mapb.nra;
875 evaluate_coord(out->u.r, pos, XX);
879 * See sel_updatefunc() for description of the parameters.
880 * \p data is not used.
882 * Returns the \p y coordinate for each atom in \p out->u.r.
885 evaluate_y(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
886 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
888 out->nr = pos->m.mapb.nra;
889 evaluate_coord(out->u.r, pos, YY);
893 * See sel_updatefunc() for description of the parameters.
894 * \p data is not used.
896 * Returns the \p z coordinate for each atom in \p out->u.r.
899 evaluate_z(t_topology * /*top*/, t_trxframe * /*fr*/, t_pbc * /*pbc*/,
900 gmx_ana_pos_t *pos, gmx_ana_selvalue_t *out, void * /*data*/)
902 out->nr = pos->m.mapb.nra;
903 evaluate_coord(out->u.r, pos, ZZ);