2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2009, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /*! \page libtrajana Library for trajectory analysis
40 * This is a trajectory analysis library for Gromacs.
42 * The main features of the library are:
43 * - \subpage selengine "Flexible handling of textual selections" that can
44 * also be dynamic, i.e., depend of the trajectory frame through
45 * positions of the particles.
46 * Selections evaluate directly to positions, which can then be used in
47 * the analysis program.
48 * - \subpage selmethods "Custom selection keywords" can also be easily
49 * implemented, also on a tool-by-tool basis.
50 * - Efficient \subpage nbsearch "neighborhood searching"
51 * (currently not very efficient...).
52 * - \subpage displacements "On-the-fly calculation of displacement vectors"
53 * during a single pass over the trajectory.
54 * - \subpage histograms "Calculation of histograms" with error estimates
55 * through block averaging.
57 * The functions also automate several things common to most analysis programs
58 * such as making molecules whole if required, reading input files, and
59 * setting up index groups for analysis.
60 * The library unifies the structure of analysis programs (at least a bit)
61 * and makes it easier to add common functionality to all analysis programs.
64 * \section main_using Using the library
66 * There is a \ref share/template/template.c "template" for writing
67 * analysis programs, the documentation for it and links from there should
68 * help getting started.
72 * \section libtrajana_impl Implementation details
74 * Some internal implementation details of the library are documented on
76 * - \subpage poscalcengine
77 * - \subpage selparser
78 * - \subpage selcompiler
80 /*! \page selengine Text-based selections
82 * \section selection_basics Basics
84 * Selections are enabled automatically for an analysis program that uses
85 * the library. The selection syntax is described in an online help that is
86 * accessible from all tools that use the library.
87 * By default, dynamic selections are allowed, and the user can freely
88 * choose whether to analyze atoms or centers of mass/geometry of
90 * These defaults, as well as some others, can be changed by specifying
91 * flags for gmx_ana_traj_create().
93 * The analysis program can then access the selected positions for each frame
94 * through a \p gmx_ana_selection_t array that is passed to the frame
95 * analysis function (see gmx_analysisfunc()).
96 * As long as the analysis program is written such that it does not assume
97 * that the number of positions or the atoms in the groups groups remain
98 * constant, any kind of selection expression can be used.
100 * Some analysis programs may require a special structure for the index groups
101 * (e.g., \c g_angle requires the index group to be made of groups of three or
103 * For such programs, it is up to the user to provide a proper selection
104 * expression that always returns such positions.
105 * Such analysis program can define \ref ANA_REQUIRE_WHOLE to make the
106 * default behavior appropriate for the most common uses where the groups
107 * should consist of atoms within a single residue/molecule.
109 * \section selection_methods Implementing new keywords
111 * New selection keywords can be easily implemented, either directly into the
112 * library or as part of analysis programs (the latter may be useful for
113 * testing or methods very specific to some analysis).
114 * For both cases, you should first create a \c gmx_ana_selmethod_t structure
115 * and fill it with the necessary information.
116 * Details can be found on a separate page: \ref selmethods.
119 * \brief Implementation of functions in trajana.h.
121 /*! \internal \dir src/gmxlib/trajana
122 * \brief Source code for common trajectory analysis functions.
124 * Selection handling is found in \ref src/gmxlib/selection.
138 #include <statutil.h>
139 #include <typedefs.h>
144 #include <selection.h>
145 #include <selmethod.h>
149 * Data structure for trajectory analysis tools.
151 struct gmx_ana_traj_t
154 * Flags that alter the behavior of the analysis library.
156 * This variable stores the flags passed to gmx_ana_traj_create() for use
157 * in the other functions.
160 /** Number of input reference groups. */
163 * Number of input analysis groups.
165 * This is the number of groups in addition to the reference groups
167 * If -1, any number of groups (at least one) is acceptable.
171 * Flags for init_first_frame() to specify what to read from the
175 /** TRUE if molecules should be made whole for each frame. */
178 * TRUE if periodic boundary conditions should be used.
180 * If the flag is FALSE, the \p pbc pointer passed to gmx_analysisfunc()
185 /** Name of the trajectory file (NULL if not provided). */
187 /** Name of the topology file (NULL if no topology loaded). */
189 /** Non-NULL name of the topology file. */
190 char *topfile_notnull;
191 /** Name of the index file (NULL if no index file provided). */
193 /** Name of the selection file (NULL if no file provided). */
195 /** The selection string (NULL if not provided). */
198 /** The topology structure, or \p NULL if no topology loaded. */
200 /** TRUE if full tpx file was loaded, FALSE otherwise. */
202 /** Coordinates from the topology (see \p bTopX). */
204 /** The box loaded from the topology file. */
206 /** The ePBC field loaded from the topology file. */
209 /** The current frame, or \p NULL if no frame loaded yet. */
211 /** Used to store the status variable from read_first_frame(). */
213 /** The number of frames read. */
216 /** Number of elements in the \p sel array. */
219 * Array of selection information (one element for each index group).
221 * After the call to gmx_ana_init_selections(), this array contains
222 * information about the selections the user has provided.
223 * The array contains \p ngrps elements;
224 * if \p nanagrps was -1, the number may vary.
225 * See \c gmx_ana_selection_t for details of the contents.
227 * The largest possible index groups for dynamic selections can be found
228 * in \p sel[i]->g, i.e., the program can assume that any index group
229 * passed to gmx_analysisfunc() is a subset of the provided group.
230 * After gmx_ana_do(), the same groups can be used in post-processing.
232 gmx_ana_selection_t **sel;
233 /** Array of names of the selections for convenience. */
235 /** Position calculation data. */
236 gmx_ana_poscalc_coll_t *pcc;
237 /** Selection data. */
238 gmx_ana_selcollection_t *sc;
240 /** Data for statutil.c utilities. */
244 /** Loads the topology. */
245 static int load_topology(gmx_ana_traj_t *d, gmx_bool bReq);
246 /** Loads the first frame and does some checks. */
247 static int init_first_frame(gmx_ana_traj_t *d);
249 static int add_fnmarg(int nfile, t_filenm *fnm, t_filenm *fnm_add)
251 memcpy(&(fnm[nfile]), fnm_add, sizeof(*fnm_add));
255 /* Copied from src/gmxlib/statutil.c */
257 add_parg(int npargs, t_pargs *pa, t_pargs *pa_add)
259 memcpy(&(pa[npargs]), pa_add, sizeof(*pa_add));
264 * \param[out] data Trajectory analysis data structure poitner to initialize.
265 * \param[in] flags Combination of flags (see \ref analysis_flags).
266 * \returns 0 on success.
269 gmx_ana_traj_create(gmx_ana_traj_t **data, unsigned long flags)
278 d->frflags = TRX_NEED_X;
300 d->topfile_notnull = NULL;
301 rc = gmx_ana_poscalc_coll_create(&d->pcc);
308 rc = gmx_ana_selcollection_create(&d->sc, d->pcc);
311 gmx_ana_poscalc_coll_free(d->pcc);
324 * \param d Trajectory analysis data to free.
327 gmx_ana_traj_free(gmx_ana_traj_t *d)
333 sfree(d->topfile_notnull);
343 /* Gromacs does not seem to have a function for freeing frame data */
351 gmx_ana_selcollection_free(d->sc);
352 gmx_ana_poscalc_coll_free(d->pcc);
354 output_env_done(d->oenv);
359 * \param[in,out] d Trajectory analysis data structure.
360 * \param[in] flags Additional flags to set.
361 * \returns 0 on success, a non-zero error code on error.
363 * Currently there are no checks whether the flags make sense or not.
366 gmx_ana_add_flags(gmx_ana_traj_t *d, unsigned long flags)
373 * \param[in,out] d Trajectory analysis data structure.
374 * \param[in] bPBC TRUE if periodic boundary conditions should be used.
375 * \returns 0 on success.
377 * If this function is called before parse_trjana_args(), it sets the default
378 * for whether PBC are used in the analysis or not.
379 * If \ref ANA_NOUSER_PBC is not set, a command-line option is provided for the
380 * user to override the default value.
381 * If called after parse_trjana_args(), it overrides the setting provided by
382 * the user or an earlier call.
384 * If this function is not called, the default is to use PBC.
386 * If PBC are not used, the \p pbc pointer passed to gmx_analysisfunc()
389 * \see \ref ANA_NOUSER_PBC
392 gmx_ana_set_pbc(gmx_ana_traj_t *d, gmx_bool bPBC)
399 * \param[in] d Trajectory analysis data structure.
400 * \returns TRUE if periodic boundary conditions are set to be used.
403 gmx_ana_has_pbc(gmx_ana_traj_t *d)
409 * \param[in,out] d Trajectory analysis data structure.
410 * \param[in] bRmPBC TRUE if molecules should be made whole.
411 * \returns 0 on success.
413 * If this function is called before parse_trjana_args(), it sets the default
414 * for whether molecules are made whole.
415 * If \ref ANA_NOUSER_RMPBC is not set, a command-line option is provided for
416 * the user to override the default value.
417 * If called after parse_trjana_args(), it overrides the setting provided by
418 * the user or an earlier call.
420 * If this function is not called, the default is to make molecules whole.
422 * The main use of this function is to call it with FALSE if your analysis
423 * program does not require whole molecules as this can increase the
425 * In such a case, you can also specify \ref ANA_NOUSER_RMPBC to not to
426 * confuse the user with an option that would only slow the program
429 * \see \ref ANA_NOUSER_RMPBC
432 gmx_ana_set_rmpbc(gmx_ana_traj_t *d, gmx_bool bRmPBC)
439 * \param[in,out] d Trajectory analysis data structure.
440 * \param[in] frflags Flags for what to read from the trajectory file.
441 * \returns 0 on success, an error code on error.
443 * The TRX_NEED_X flag is always set.
444 * If the analysis tools needs some other information (velocities, forces),
445 * it can call this function to load additional information from the
449 gmx_ana_set_frflags(gmx_ana_traj_t *d, int frflags)
453 gmx_call("cannot set trajectory flags after initializing selections");
458 gmx_call("cannot set trajectory flags after the first frame has been read");
461 frflags |= TRX_NEED_X;
462 d->frflags = frflags;
467 * \param[in,out] d Trajectory analysis data structure.
468 * \param[in] nrefgrps Number of reference groups required.
469 * \returns 0 on success, a non-zero error code on error.
471 * \p nrefgrps should be a non-negative integer.
472 * If this function is not called (or \p nrefgrps is 0), all selections are
473 * treated as reference groups.
476 gmx_ana_set_nrefgrps(gmx_ana_traj_t *d, int nrefgrps)
481 gmx_incons("number of reference groups is negative");
484 d->nrefgrps = nrefgrps;
489 * \param[in,out] d Trajectory analysis data structure.
490 * \param[in] nanagrps Number of analysis groups required
491 * (-1 stands for any number of analysis groups).
492 * \returns 0 on success, a non-zero error code on error.
494 * \p nanagrps should be a positive integer or -1.
495 * In the latter case, any number of groups (but at least one) is acceptable.
496 * gmx_ana_get_nanagrps() can be used to access the actual value after
497 * gmx_ana_init_selections() has been called.
498 * If this function is not called, a single analysis group is expected.
501 gmx_ana_set_nanagrps(gmx_ana_traj_t *d, int nanagrps)
503 if (nanagrps <= 0 && nanagrps != -1)
506 gmx_incons("number of analysis groups is invalid");
509 d->nanagrps = nanagrps;
514 * \param[in] d Trajectory analysis data structure.
515 * \param[out] ngrps Total number of selections specified by the user.
516 * \returns 0 on success, a non-zero error code on error.
518 * The value stored in \p *ngrps is the sum of the number of reference groups
519 * and the number of analysis groups.
520 * If a specific number (not -1) of analysis groups has been set with
521 * gmx_ana_set_nanagrps(), the value is always the sum of the values provided
522 * to gmx_ana_set_nrefgrps() and gmx_ana_set_nanagrps().
524 * Should only be called after gmx_ana_init_selections().
527 gmx_ana_get_ngrps(gmx_ana_traj_t *d, int *ngrps)
529 if (d->nanagrps == -1)
532 gmx_call("gmx_ana_init_selections() not called");
535 *ngrps = d->nrefgrps + d->nanagrps;
540 * \param[in] d Trajectory analysis data structure.
541 * \param[out] nanagrps Number of analysis groups specified by the user.
542 * \returns 0 on success, a non-zero error code on error.
544 * If a specific number (not -1) of analysis groups has been set with
545 * gmx_ana_set_nanagrps(), the value is always the same value.
546 * Hence, you only need to call this function if gmx_ana_set_nanagrps() has
547 * been called with \p nanagrps set to -1.
549 * Should only be called after gmx_ana_init_selections().
552 gmx_ana_get_nanagrps(gmx_ana_traj_t *d, int *nanagrps)
554 if (d->nanagrps == -1)
557 gmx_call("gmx_ana_init_selections() not called");
560 *nanagrps = d->nanagrps;
565 * \param[in] d Trajectory analysis data structure.
566 * \param[in] i Ordinal number of the reference selection to get.
567 * \param[out] sel Selection object for the \p i'th reference group.
568 * \returns 0 on success, a non-zero error code on error.
570 * The pointer returned in \p *sel should not be freed.
571 * Should only be called after gmx_ana_init_selections().
574 gmx_ana_get_refsel(gmx_ana_traj_t *d, int i, gmx_ana_selection_t **sel)
576 if (i < 0 || i >= d->nrefgrps)
579 gmx_call("invalid reference group number");
582 *sel = gmx_ana_selcollection_get_selection(d->sc, i);
585 gmx_incons("gmx_ana_init_selections() not called");
592 * \param[in] d Trajectory analysis data structure.
593 * \param[out] sel Array of analysis selections.
594 * \returns 0 on success, a non-zero error code on error.
596 * The pointer returned in \p *sel should not be freed.
597 * Should only be called after gmx_ana_init_selections().
600 gmx_ana_get_anagrps(gmx_ana_traj_t *d, gmx_ana_selection_t ***sel)
605 gmx_incons("gmx_ana_init_selections() not called");
613 * \param[in] d Trajectory analysis data structure.
614 * \param[out] grpnames Array of selection names.
615 * \returns 0 on success, a non-zero error code on error.
617 * The pointer returned in \p *grpnames should not be freed.
618 * Should only be called after gmx_ana_init_selections().
621 gmx_ana_get_grpnames(gmx_ana_traj_t *d, char ***grpnames)
626 gmx_call("gmx_ana_init_selections() not called");
629 *grpnames = d->grpnames;
634 * \param[in] d Trajectory analysis data structure.
635 * \param[out] sc Selection collection object.
636 * \returns 0 on success.
638 * This function is mostly useful for debugging purposes.
639 * The information commonly required in analysis programs is accessible
640 * more conveniently through other means.
642 * The pointer returned in \p *sc should not be freed.
643 * Can be called at any point.
646 gmx_ana_get_selcollection(gmx_ana_traj_t *d, gmx_ana_selcollection_t **sc)
653 * \param[in,out] d Trajectory analysis data structure.
654 * \returns 0 on success, a non-zero error code on error.
656 * This function should be called first in the analysis program, much in
657 * the same way as parse_common_args() in traditional Gromacs analysis
658 * programs. It adds some command-line arguments of its own, and uses
659 * parse_common_args() to parse the command-line arguments.
660 * It also loads topology information if required or if a topology is
661 * provided on the command line.
662 * Selection handling is also initialized if it is enabled and
663 * the user has selected it on the command line.
665 * The rest of the parameters are passed on to the Gromacs routine
666 * parse_common_args(), and the use of this function should be identical
667 * to parse_common_args(), with the exception that for \p pca_flags,
668 * \p PCA_CAN_TIME and \p PCA_BE_NICE flags are automatically added.
683 parse_trjana_args(gmx_ana_traj_t *d,
684 int *argc, char *argv[], unsigned long pca_flags,
685 int nfile, t_filenm fnm[], int npargs, t_pargs *pa,
686 int ndesc, const char **desc,
687 int nbugs, const char **bugs,
690 t_filenm *all_fnm = NULL;
693 t_pargs *all_pa = NULL;
700 t_filenm def_fnm[] = {
701 {efTRX, NULL, NULL, ffOPTRD},
702 {efTPS, NULL, NULL, ffREAD},
703 {efDAT, "-sf", "selection", ffOPTRD},
704 {efNDX, NULL, NULL, ffOPTRD},
706 gmx_bool bPBC = TRUE;
708 {"-pbc", FALSE, etBOOL, {&bPBC},
709 "Use periodic boundary conditions for distance calculation"},
711 gmx_bool bRmPBC = TRUE;
712 t_pargs rmpbc_pa[] = {
713 {"-rmpbc", FALSE, etBOOL, {&bRmPBC},
714 "Make molecules whole for each frame"},
716 char *selection = NULL;
717 const char **rpost = NULL;
718 gmx_bool bSelDump = FALSE;
720 {"-select", FALSE, etSTR, {&selection},
721 "Selection string (use 'help' for help). Note that the "
722 "whole selection string will need to be quoted so that "
723 "your shell will pass it in as a string. Example: "
724 "[TT]g_select -select '\"Nearby water\" resname SOL "
725 "and within 0.25 of group Protein'[tt]"},
726 {"-seldebug", FALSE, etBOOL, {&bSelDump},
727 "HIDDENPrint out the parsed and compiled selection trees"},
729 t_pargs dsel_pa[] = {
730 {"-selrpos", FALSE, etENUM, {NULL},
731 "Selection reference position"},
733 const char **spost = NULL;
734 t_pargs selpt_pa[] = {
735 {"-seltype", FALSE, etENUM, {NULL},
736 "Default analysis positions"},
738 #define MAX_PA asize(sel_pa)+asize(dsel_pa)+5
742 gmx_incons("number of reference groups is negative");
746 if (d->flags & ANA_DEBUG_SELECTION)
751 rpost = gmx_ana_poscalc_create_type_enum(!(d->flags & ANA_REQUIRE_WHOLE));
756 spost = gmx_ana_poscalc_create_type_enum(TRUE);
763 /* Construct the file name argument array */
764 max_fnm = nfile + asize(def_fnm);
765 snew(all_fnm, max_fnm);
767 if (!(d->flags & ANA_REQUIRE_TOP))
769 def_fnm[1].flag |= ffOPT;
771 snew(fnm_map, nfile);
772 for (k = 0; k < nfile; ++k)
777 for (i = 0; i < asize(def_fnm); ++i)
779 for (k = 0; k < nfile; ++k)
781 if (fnm_map[k] == -1 && def_fnm[i].opt == NULL
782 && fnm[k].opt == NULL && fnm[k].ftp == def_fnm[i].ftp)
790 nfall = add_fnmarg(nfall, all_fnm, &(fnm[k]));
794 nfall = add_fnmarg(nfall, all_fnm, &(def_fnm[i]));
798 for (k = 0; k < nfile; ++k)
800 if (fnm_map[k] == -1)
803 nfall = add_fnmarg(nfall, all_fnm, &(fnm[k]));
807 /* Construct the argument array */
808 max_pa = npargs + MAX_PA;
809 snew(all_pa, max_pa);
812 if (!(d->flags & ANA_NOUSER_RMPBC))
814 for (i = 0; i < asize(rmpbc_pa); ++i)
816 npall = add_parg(npall, all_pa, &(rmpbc_pa[i]));
819 if (!(d->flags & ANA_NOUSER_PBC))
821 for (i = 0; i < asize(pbc_pa); ++i)
823 npall = add_parg(npall, all_pa, &(pbc_pa[i]));
827 for (i = 0; i < asize(sel_pa); ++i)
829 npall = add_parg(npall, all_pa, &(sel_pa[i]));
831 if (!(d->flags & ANA_NO_DYNSEL))
833 dsel_pa[0].u.c = rpost;
834 for (i = 0; i < asize(dsel_pa); ++i)
836 npall = add_parg(npall, all_pa, &(dsel_pa[i]));
840 if (!(d->flags & ANA_ONLY_ATOMPOS))
842 selpt_pa[0].u.c = spost;
843 for (i = 0; i < asize(selpt_pa); ++i)
845 npall = add_parg(npall, all_pa, &(selpt_pa[i]));
849 for (k = 0; k < npargs; ++k)
851 npall = add_parg(npall, all_pa, &(pa[k]));
854 pca_flags |= PCA_CAN_TIME | PCA_BE_NICE;
855 parse_common_args(argc, argv, pca_flags, nfall, all_fnm, npall, all_pa,
856 ndesc, desc, nbugs, bugs,oenv);
859 /* Process our own options.
860 * Make copies of file names for easier memory management. */
861 tmp_fnm = ftp2fn_null(efTRX, nfall, all_fnm);
862 d->trjfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
863 tmp_fnm = ftp2fn_null(efTPS, nfall, all_fnm);
864 d->topfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
865 d->topfile_notnull = strdup(ftp2fn(efTPS, nfall, all_fnm));
866 tmp_fnm = ftp2fn_null(efNDX, nfall, all_fnm);
867 d->ndxfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
868 if (!(d->flags & ANA_NOUSER_RMPBC))
872 if (!(d->flags & ANA_NOUSER_PBC))
876 d->selection = selection;
877 tmp_fnm = opt2fn_null("-sf", nfall, all_fnm);
878 d->selfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
880 /* Copy the results back */
881 for (k = 0; k < nfile; ++k)
883 memcpy(&(fnm[k]), &(all_fnm[fnm_map[k]]), sizeof(fnm[k]));
884 /* Delegate responsibility of freeing the file names to caller. */
885 all_fnm[fnm_map[k]].nfiles = 0;
886 all_fnm[fnm_map[k]].fns = NULL;
888 for (i = 0, k = npall - npargs; i < (size_t)npargs; ++i, ++k)
890 memcpy(&(pa[i]), &(all_pa[k]), sizeof(pa[i]));
893 /* Free memory we have allocated. */
894 done_filenms(nfall, all_fnm);
899 if (!(d->flags & ANA_NO_DYNSEL))
901 gmx_ana_selcollection_set_refpostype(d->sc, rpost[0]);
905 gmx_ana_selcollection_set_refpostype(d->sc, rpost[1]);
910 d->flags |= ANA_DEBUG_SELECTION;
914 d->flags &= ~ANA_DEBUG_SELECTION;
917 if (!(d->flags & ANA_ONLY_ATOMPOS))
919 gmx_ana_selcollection_set_outpostype(d->sc, spost[0], d->flags & ANA_USE_POSMASK);
923 gmx_ana_selcollection_set_outpostype(d->sc, spost[1], d->flags & ANA_USE_POSMASK);
927 /* Check if the user requested help on selections.
928 * If so, call gmx_ana_init_selections() to print the help and exit. */
929 if (selection && strncmp(selection, "help", 4) == 0)
931 gmx_ana_init_selections(d);
934 /* If no trajectory file is given, we need to set some flags to be able
935 * to prepare a frame from the loaded topology information. Also, check
936 * that a topology is provided. */
941 gmx_input("No trajectory or topology provided, nothing to do!");
944 d->flags |= ANA_REQUIRE_TOP;
945 d->flags |= ANA_USE_TOPX;
948 /* Load the topology if so requested. */
949 rc = load_topology(d, (d->flags & ANA_REQUIRE_TOP));
955 /* Initialize the selections/index groups */
956 if (!(d->flags & ANA_USER_SELINIT))
958 rc = gmx_ana_init_selections(d);
965 * \param[in,out] d Trajectory analysis data structure.
966 * \param[in] bReq If TRUE, topology loading is forced.
967 * \returns 0 on success, a non-zero error code on error.
969 * Initializes the \c gmx_ana_traj_t::top, \c gmx_ana_traj_t::bTop,
970 * \c gmx_ana_traj_t::boxtop and \c gmx_ana_traj_t::ePBC fields of the
971 * analysis structure.
972 * If \p bReq is TRUE, the topology is loaded even if it is not given on
975 * The function can be called multiple times safely; extra calls are
978 static int load_topology(gmx_ana_traj_t *d, gmx_bool bReq)
982 /* Return if topology already loaded */
988 if (d->topfile || bReq)
991 d->bTop = read_tps_conf(d->topfile_notnull, title, d->top,
992 &d->ePBC, &d->xtop, NULL, d->boxtop, TRUE);
993 if (!(d->flags & ANA_USE_TOPX))
1003 * \param[in] d Trajectory analysis data structure.
1004 * \param[in] bReq If TRUE, topology loading is forced.
1005 * \param[out] top Topology data pointer to initialize.
1006 * \param[out] bTop TRUE if a full tpx file was loaded, FALSE otherwise
1007 * (can be NULL, in which case it is not used).
1008 * \returns 0 on success, a non-zero error code on error.
1010 * If \ref ANA_REQUIRE_TOP has not been specified and \p bReq is FALSE,
1011 * the pointer stored in \p *top may be NULL if no topology has been provided
1012 * on the command line.
1014 * The pointer returned in \p *top should not be freed.
1017 gmx_ana_get_topology(gmx_ana_traj_t *d, gmx_bool bReq, t_topology **top, gmx_bool *bTop)
1021 rc = load_topology(d, bReq);
1036 * \param[in] d Trajectory analysis data structure.
1037 * \param[out] x Topology data pointer to initialize.
1038 * (can be NULL, in which case it is not used).
1039 * \param[out] box Box size from the topology file
1040 * (can be NULL, in which case it is not used).
1041 * \param[out] ePBC The ePBC field from the topology
1042 * (can be NULL, in which case it is not used).
1043 * \returns 0 on success, a non-zero error code on error.
1045 * If \ref ANA_USE_TOPX has not been specified, the \p x parameter should be
1048 * The pointer returned in \p *x should not be freed.
1051 gmx_ana_get_topconf(gmx_ana_traj_t *d, rvec **x, matrix box, int *ePBC)
1055 copy_mat(d->boxtop, box);
1063 if (!(d->flags & ANA_USE_TOPX))
1065 gmx_incons("topology coordinates requested by ANA_USE_TOPX not set");
1075 * Loads default index groups from a selection file.
1077 * \param[in,out] d Trajectory analysis data structure.
1078 * \param[out] grps Pointer to receive the default groups.
1079 * \returns 0 on success, a non-zero error code on error.
1082 init_default_selections(gmx_ana_traj_t *d, gmx_ana_indexgrps_t **grps)
1084 gmx_ana_selcollection_t *sc;
1086 int nr, nr_notempty, i;
1089 /* If an index file is provided, just load it and exit. */
1092 gmx_ana_indexgrps_init(grps, d->top, d->ndxfile);
1095 /* Initialize groups to NULL if we return prematurely. */
1097 /* Return immediately if no topology provided. */
1103 /* Find the default selection file, return if none found. */
1104 fnm = low_gmxlibfn("defselection.dat", TRUE, FALSE);
1110 /* Create a temporary selection collection. */
1111 rc = gmx_ana_selcollection_create(&sc, d->pcc);
1117 rc = gmx_ana_selmethod_register_defaults(sc);
1120 gmx_ana_selcollection_free(sc);
1122 gmx_fatal(FARGS, "default selection method registration failed");
1125 /* FIXME: It would be better to not have the strings here hard-coded. */
1126 gmx_ana_selcollection_set_refpostype(sc, "atom");
1127 gmx_ana_selcollection_set_outpostype(sc, "atom", FALSE);
1129 /* Parse and compile the file with no external groups. */
1130 rc = gmx_ana_selcollection_parse_file(sc, fnm, NULL);
1134 gmx_ana_selcollection_free(sc);
1135 fprintf(stderr, "\nWARNING: default selection(s) could not be parsed\n");
1138 gmx_ana_selcollection_set_topology(sc, d->top, -1);
1139 rc = gmx_ana_selcollection_compile(sc);
1142 gmx_ana_selcollection_free(sc);
1143 fprintf(stderr, "\nWARNING: default selection(s) could not be compiled\n");
1147 /* Count the non-empty groups and check that there are no dynamic
1149 nr = gmx_ana_selcollection_get_count(sc);
1151 for (i = 0; i < nr; ++i)
1153 gmx_ana_selection_t *sel;
1155 sel = gmx_ana_selcollection_get_selection(sc, i);
1158 fprintf(stderr, "\nWARNING: dynamic default selection ignored\n");
1160 else if (sel->g->isize > 0)
1166 /* Copy the groups to the output structure */
1167 gmx_ana_indexgrps_alloc(grps, nr_notempty);
1169 for (i = 0; i < nr; ++i)
1171 gmx_ana_selection_t *sel;
1173 sel = gmx_ana_selcollection_get_selection(sc, i);
1174 if (!sel->bDynamic && sel->g->isize > 0)
1178 g = gmx_ana_indexgrps_get_grp(*grps, nr_notempty);
1179 gmx_ana_index_copy(g, sel->g, TRUE);
1180 g->name = strdup(sel->name);
1185 gmx_ana_selcollection_free(sc);
1190 * \param[in,out] d Trajectory analysis data structure.
1191 * \returns 0 on success, a non-zero error code on error.
1193 * Initializes the selection data in \c gmx_ana_traj_t based on
1194 * the selection options and/or index files provided on the command line.
1196 * This function is called automatically by parse_trjana_args() and should
1197 * not be called directly unless \ref ANA_USER_SELINIT is specified.
1199 * \see ANA_USER_SELINIT
1202 gmx_ana_init_selections(gmx_ana_traj_t *d)
1207 gmx_ana_indexgrps_t *grps;
1210 gmx_bool bInteractive;
1215 gmx_call("init_selections called more than once\n"
1216 "perhaps you forgot ANA_USER_SELINIT");
1220 gmx_ana_selcollection_set_veloutput(d->sc,
1221 d->frflags & (TRX_READ_V | TRX_NEED_V));
1222 gmx_ana_selcollection_set_forceoutput(d->sc,
1223 d->frflags & (TRX_READ_F | TRX_NEED_F));
1224 /* Check if we need some information from the topology */
1225 if (gmx_ana_selcollection_requires_top(d->sc))
1227 rc = load_topology(d, TRUE);
1233 /* Initialize the default selection methods */
1234 rc = gmx_ana_selmethod_register_defaults(d->sc);
1237 gmx_fatal(FARGS, "default selection method registration failed");
1240 /* Initialize index groups.
1241 * We ignore the return value to continue without the default groups if
1242 * there is an error there. */
1243 init_default_selections(d, &grps);
1244 /* Parse the selections */
1245 bStdIn = (d->selfile && d->selfile[0] == '-' && d->selfile[1] == 0)
1246 || (d->selection && d->selection[0] == 0)
1247 || (!d->selfile && !d->selection);
1248 /* Behavior is not very pretty if we cannot check for interactive input,
1249 * but at least it should compile and work in most cases. */
1250 #ifdef HAVE_UNISTD_H
1251 bInteractive = bStdIn && isatty(fileno(stdin));
1253 bInteractive = bStdIn;
1255 if (bStdIn && bInteractive)
1257 /* Parse from stdin */
1258 /* First we parse the reference groups if there are any */
1259 if (d->nrefgrps > 0)
1261 fprintf(stderr, "\nSpecify ");
1262 if (d->nrefgrps == 1)
1264 fprintf(stderr, "a reference selection");
1268 fprintf(stderr, "%d reference selections", d->nrefgrps);
1270 fprintf(stderr, ":\n");
1271 fprintf(stderr, "(one selection per line, 'help' for help)\n");
1272 rc = gmx_ana_selcollection_parse_stdin(d->sc, d->nrefgrps, grps, TRUE);
1273 nr = gmx_ana_selcollection_get_count(d->sc);
1274 if (rc != 0 || nr != d->nrefgrps)
1276 gmx_ana_traj_free(d);
1277 gmx_input("unrecoverable error in selection parsing");
1281 /* Then, we parse the analysis groups */
1282 fprintf(stderr, "\nSpecify ");
1283 if (d->nanagrps == 1)
1285 fprintf(stderr, "a selection");
1287 else if (d->nanagrps == -1)
1289 fprintf(stderr, "any number of selections");
1293 fprintf(stderr, "%d selections", d->nanagrps);
1295 fprintf(stderr, " for analysis:\n");
1296 fprintf(stderr, "(one selection per line, 'help' for help%s)\n",
1297 d->nanagrps == -1 ? ", Ctrl-D to end" : "");
1298 rc = gmx_ana_selcollection_parse_stdin(d->sc, d->nanagrps, grps, TRUE);
1299 fprintf(stderr, "\n");
1303 rc = gmx_ana_selcollection_parse_stdin(d->sc, -1, grps, FALSE);
1305 else if (d->selection)
1307 rc = gmx_ana_selcollection_parse_str(d->sc, d->selection, grps);
1311 rc = gmx_ana_selcollection_parse_file(d->sc, d->selfile, grps);
1315 gmx_ana_indexgrps_free(grps);
1319 /* Free memory for memory leak checking */
1320 gmx_ana_traj_free(d);
1321 gmx_input("selection(s) could not be parsed");
1325 /* Check the number of groups */
1326 nr = gmx_ana_selcollection_get_count(d->sc);
1329 /* TODO: Don't print this if the user has requested help */
1330 fprintf(stderr, "Nothing selected, finishing up.\n");
1331 gmx_ana_traj_free(d);
1332 /* TODO: It would be better to return some code that tells the caller
1333 * that one should exit. */
1336 if (nr <= d->nrefgrps)
1338 gmx_input("selection does not specify enough index groups");
1341 if (d->nanagrps <= 0)
1343 d->nanagrps = nr - d->nrefgrps;
1345 else if (nr != d->nrefgrps + d->nanagrps)
1347 gmx_input("selection does not specify the correct number of index groups");
1351 if (d->flags & ANA_DEBUG_SELECTION)
1353 gmx_ana_selcollection_print_tree(stderr, d->sc, FALSE);
1355 if (gmx_ana_selcollection_requires_top(d->sc))
1357 rc = load_topology(d, TRUE);
1369 rc = init_first_frame(d);
1374 natoms = d->fr->natoms;
1376 gmx_ana_selcollection_set_topology(d->sc, d->top, natoms);
1377 rc = gmx_ana_selcollection_compile(d->sc);
1380 /* Free memory for memory leak checking */
1381 gmx_ana_traj_free(d);
1382 gmx_input("selection could not be compiled");
1385 /* Create the selection array */
1386 d->ngrps = gmx_ana_selcollection_get_count(d->sc);
1387 if (!(d->flags & ANA_USE_FULLGRPS))
1389 d->ngrps -= d->nrefgrps;
1391 snew(d->sel, d->ngrps);
1392 for (i = 0; i < d->ngrps; ++i)
1394 if (d->flags & ANA_USE_FULLGRPS)
1396 d->sel[i] = gmx_ana_selcollection_get_selection(d->sc, i);
1400 d->sel[i] = gmx_ana_selcollection_get_selection(d->sc, i + d->nrefgrps);
1403 if (d->flags & ANA_DEBUG_SELECTION)
1405 fprintf(stderr, "\n");
1406 gmx_ana_selcollection_print_tree(stderr, d->sc, FALSE);
1407 fprintf(stderr, "\n");
1408 gmx_ana_poscalc_coll_print_tree(stderr, d->pcc);
1409 fprintf(stderr, "\n");
1412 /* Initialize the position evaluation */
1413 gmx_ana_poscalc_init_eval(d->pcc);
1414 if (d->flags & ANA_DEBUG_SELECTION)
1416 gmx_ana_poscalc_coll_print_tree(stderr, d->pcc);
1417 fprintf(stderr, "\n");
1420 /* Check that dynamic selections are not provided if not allowed */
1421 if (d->flags & ANA_NO_DYNSEL)
1423 for (i = 0; i < d->nrefgrps + d->nanagrps; ++i)
1425 gmx_ana_selection_t *sel;
1427 sel = gmx_ana_selcollection_get_selection(d->sc, i);
1430 gmx_fatal(FARGS, "%s does not support dynamic selections",
1436 /* Check that non-atom positions are not provided if not allowed.
1437 * TODO: It would be better to have these checks in the parser. */
1438 if (d->flags & ANA_ONLY_ATOMPOS)
1440 for (i = 0; i < d->nanagrps; ++i)
1442 gmx_ana_selection_t *sel;
1444 sel = gmx_ana_selcollection_get_selection(d->sc, i + d->nrefgrps);
1445 if (sel->p.m.type != INDEX_ATOM)
1447 gmx_fatal(FARGS, "%s does not support non-atom positions",
1453 /* Create the names array */
1454 snew(d->grpnames, d->ngrps);
1455 for (i = 0; i < d->ngrps; ++i)
1457 d->grpnames[i] = gmx_ana_selection_name(d->sel[i]);
1464 * \param[in,out] d Trajectory analysis data structure.
1465 * \param[in] type Type of covered fractions to calculate.
1466 * \returns 0 on success, a non-zero error code on error.
1468 * By default, covered fractions are not calculated.
1469 * If this function is called, the covered fraction calculation is
1470 * initialize to calculate the fractions of type \p type for each selection.
1471 * The function must be called after gmx_ana_init_selections() has been
1474 * For more fine-grained control of the calculation, you can use
1475 * gmx_ana_selection_init_coverfrac(): if you initialize some selections
1476 * this function to something else than CFRAC_NONE before calling
1477 * gmx_ana_init_coverfrac(), these settings are not overwritten.
1479 * You only need to call this function if your program needs to have
1480 * information about the covered fractions, e.g., for normalization.
1482 * \see gmx_ana_selection_init_coverfrac()
1485 gmx_ana_init_coverfrac(gmx_ana_traj_t *d, e_coverfrac_t type)
1489 if (type == CFRAC_NONE)
1494 for (g = 0; g < d->ngrps; ++g)
1496 if (d->sel[g]->cfractype == CFRAC_NONE)
1498 gmx_ana_selection_init_coverfrac(d->sel[g], type);
1505 * \param[in] out Output file.
1506 * \param[in] d Trajectory analysis data structure.
1507 * \returns 0 on success, a non-zero error code on error.
1509 int xvgr_selections(FILE *out, gmx_ana_traj_t *d)
1511 xvgr_selcollection(out, d->sc, d->oenv);
1516 * \param[in,out] d Trajectory analysis data structure.
1517 * \returns 0 on success, a non-zero error code on error.
1519 static int init_first_frame(gmx_ana_traj_t *d)
1523 /* Return if we have already initialized the trajectory */
1529 d->frflags |= TRX_NEED_X;
1535 if (!read_first_frame(d->oenv, &d->status, d->trjfile, d->fr, d->frflags))
1537 gmx_input("could not read coordinates from trajectory");
1541 if (d->top && d->fr->natoms > d->top->atoms.nr)
1543 gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
1544 d->fr->natoms, d->top->atoms.nr);
1547 /* check index groups */
1548 for (i = 0; i < d->ngrps; ++i)
1550 gmx_ana_index_check(d->sel[i]->g, d->fr->natoms);
1555 /* Prepare a frame from topology information */
1556 /* TODO: Initialize more of the fields */
1557 if (d->frflags & (TRX_NEED_V))
1559 gmx_impl("Velocity reading from a topology not implemented");
1562 if (d->frflags & (TRX_NEED_F))
1564 gmx_input("Forces cannot be read from a topology");
1567 d->fr->flags = d->frflags;
1568 d->fr->natoms = d->top->atoms.nr;
1570 snew(d->fr->x, d->fr->natoms);
1571 memcpy(d->fr->x, d->xtop, sizeof(*d->fr->x)*d->fr->natoms);
1573 copy_mat(d->boxtop, d->fr->box);
1576 set_trxframe_ePBC(d->fr, d->ePBC);
1582 * \param[in,out] d Trajectory analysis data structure.
1583 * \param[out] fr First frame in the trajectory.
1584 * \returns 0 on success, a non-zero error code on error.
1586 * The pointer stored in \p *fr should not be freed by the caller.
1588 * You only need to call this function if your program needs to do some
1589 * initialization for which it requires data from the first frame.
1593 int gmx_ana_get_first_frame(gmx_ana_traj_t *d, t_trxframe **fr)
1597 rc = init_first_frame(d);
1608 * \param[in,out] d Trajectory analysis data structure.
1609 * \param[in] flags Combination of flags
1610 * (currently, there are no flags defined).
1611 * \param[in] analyze Pointer to frame analysis function.
1612 * \param data User data to be passed to \p analyze.
1613 * \returns 0 on success, a non-zero error code on error.
1615 * This function performs the actual analysis of the trajectory.
1616 * It loops through all the frames in the trajectory, and calls
1617 * \p analyze for each frame to perform the actual analysis.
1618 * Before calling \p analyze, selections are updated (if there are any).
1619 * See gmx_analysisfunc() for description of \p analyze parameters.
1621 * This function also calculates the number of frames during the run.
1623 int gmx_ana_do(gmx_ana_traj_t *d, int flags, gmx_analysisfunc analyze, void *data)
1628 gmx_rmpbc_t gpbc=NULL;
1630 rc = init_first_frame(d);
1636 ppbc = d->bPBC ? &pbc : 0;
1643 gpbc = gmx_rmpbc_init(&d->top->idef,d->ePBC,d->fr->natoms,d->fr->box);
1650 gmx_rmpbc_trxfr(gpbc,d->fr);
1654 set_pbc(&pbc, d->ePBC, d->fr->box);
1657 gmx_ana_poscalc_init_frame(d->pcc);
1658 rc = gmx_ana_selcollection_evaluate(d->sc, d->fr, ppbc);
1661 close_trj(d->status);
1662 gmx_fatal(FARGS, "selection evaluation failed");
1665 rc = analyze(d->top, d->fr, ppbc, d->ngrps, d->sel, data);
1668 close_trj(d->status);
1674 while (d->trjfile && read_next_frame(d->oenv, d->status, d->fr));
1677 gmx_rmpbc_done(gpbc);
1681 close_trj(d->status);
1682 fprintf(stderr, "Analyzed %d frames, last time %.3f\n",
1683 d->nframes, d->fr->time);
1687 fprintf(stderr, "Analyzed topology coordinates\n");
1690 /* Restore the maximal groups for dynamic selections */
1691 rc = gmx_ana_selcollection_evaluate_fin(d->sc, d->nframes);
1694 gmx_fatal(FARGS, "selection evaluation failed");
1701 * \param[in,out] d Trajectory analysis data structure.
1702 * \param[out] nframes Number of frames.
1703 * \returns 0 on success, a non-zero error code on error.
1705 * If called before gmx_ana_do(), the behavior is currently undefined.
1708 gmx_ana_get_nframes(gmx_ana_traj_t *d, int *nframes)
1710 *nframes = d->nframes;