3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
31 /*! \page libtrajana Library for trajectory analysis
33 * This is a trajectory analysis library for Gromacs.
35 * The main features of the library are:
36 * - \subpage selengine "Flexible handling of textual selections" that can
37 * also be dynamic, i.e., depend of the trajectory frame through
38 * positions of the particles.
39 * Selections evaluate directly to positions, which can then be used in
40 * the analysis program.
41 * - \subpage selmethods "Custom selection keywords" can also be easily
42 * implemented, also on a tool-by-tool basis.
43 * - Efficient \subpage nbsearch "neighborhood searching"
44 * (currently not very efficient...).
45 * - \subpage displacements "On-the-fly calculation of displacement vectors"
46 * during a single pass over the trajectory.
47 * - \subpage histograms "Calculation of histograms" with error estimates
48 * through block averaging.
50 * The functions also automate several things common to most analysis programs
51 * such as making molecules whole if required, reading input files, and
52 * setting up index groups for analysis.
53 * The library unifies the structure of analysis programs (at least a bit)
54 * and makes it easier to add common functionality to all analysis programs.
57 * \section main_using Using the library
59 * There is a \ref share/template/template.c "template" for writing
60 * analysis programs, the documentation for it and links from there should
61 * help getting started.
65 * \section libtrajana_impl Implementation details
67 * Some internal implementation details of the library are documented on
69 * - \subpage poscalcengine
70 * - \subpage selparser
71 * - \subpage selcompiler
73 /*! \page selengine Text-based selections
75 * \section selection_basics Basics
77 * Selections are enabled automatically for an analysis program that uses
78 * the library. The selection syntax is described in an online help that is
79 * accessible from all tools that use the library.
80 * By default, dynamic selections are allowed, and the user can freely
81 * choose whether to analyze atoms or centers of mass/geometry of
83 * These defaults, as well as some others, can be changed by specifying
84 * flags for gmx_ana_traj_create().
86 * The analysis program can then access the selected positions for each frame
87 * through a \p gmx_ana_selection_t array that is passed to the frame
88 * analysis function (see gmx_analysisfunc()).
89 * As long as the analysis program is written such that it does not assume
90 * that the number of positions or the atoms in the groups groups remain
91 * constant, any kind of selection expression can be used.
93 * Some analysis programs may require a special structure for the index groups
94 * (e.g., \c g_angle requires the index group to be made of groups of three or
96 * For such programs, it is up to the user to provide a proper selection
97 * expression that always returns such positions.
98 * Such analysis program can define \ref ANA_REQUIRE_WHOLE to make the
99 * default behavior appropriate for the most common uses where the groups
100 * should consist of atoms within a single residue/molecule.
102 * \section selection_methods Implementing new keywords
104 * New selection keywords can be easily implemented, either directly into the
105 * library or as part of analysis programs (the latter may be useful for
106 * testing or methods very specific to some analysis).
107 * For both cases, you should first create a \c gmx_ana_selmethod_t structure
108 * and fill it with the necessary information.
109 * Details can be found on a separate page: \ref selmethods.
112 * \brief Implementation of functions in trajana.h.
114 /*! \internal \dir src/gmxlib/trajana
115 * \brief Source code for common trajectory analysis functions.
117 * Selection handling is found in \ref src/gmxlib/selection.
131 #include <statutil.h>
132 #include <typedefs.h>
137 #include <selection.h>
138 #include <selmethod.h>
142 * Data structure for trajectory analysis tools.
144 struct gmx_ana_traj_t
147 * Flags that alter the behavior of the analysis library.
149 * This variable stores the flags passed to gmx_ana_traj_create() for use
150 * in the other functions.
153 /** Number of input reference groups. */
156 * Number of input analysis groups.
158 * This is the number of groups in addition to the reference groups
160 * If -1, any number of groups (at least one) is acceptable.
164 * Flags for init_first_frame() to specify what to read from the
168 /** TRUE if molecules should be made whole for each frame. */
171 * TRUE if periodic boundary conditions should be used.
173 * If the flag is FALSE, the \p pbc pointer passed to gmx_analysisfunc()
178 /** Name of the trajectory file (NULL if not provided). */
180 /** Name of the topology file (NULL if no topology loaded). */
182 /** Non-NULL name of the topology file. */
183 char *topfile_notnull;
184 /** Name of the index file (NULL if no index file provided). */
186 /** Name of the selection file (NULL if no file provided). */
188 /** The selection string (NULL if not provided). */
191 /** The topology structure, or \p NULL if no topology loaded. */
193 /** TRUE if full tpx file was loaded, FALSE otherwise. */
195 /** Coordinates from the topology (see \p bTopX). */
197 /** The box loaded from the topology file. */
199 /** The ePBC field loaded from the topology file. */
202 /** The current frame, or \p NULL if no frame loaded yet. */
204 /** Used to store the status variable from read_first_frame(). */
206 /** The number of frames read. */
209 /** Number of elements in the \p sel array. */
212 * Array of selection information (one element for each index group).
214 * After the call to gmx_ana_init_selections(), this array contains
215 * information about the selections the user has provided.
216 * The array contains \p ngrps elements;
217 * if \p nanagrps was -1, the number may vary.
218 * See \c gmx_ana_selection_t for details of the contents.
220 * The largest possible index groups for dynamic selections can be found
221 * in \p sel[i]->g, i.e., the program can assume that any index group
222 * passed to gmx_analysisfunc() is a subset of the provided group.
223 * After gmx_ana_do(), the same groups can be used in post-processing.
225 gmx_ana_selection_t **sel;
226 /** Array of names of the selections for convenience. */
228 /** Position calculation data. */
229 gmx_ana_poscalc_coll_t *pcc;
230 /** Selection data. */
231 gmx_ana_selcollection_t *sc;
233 /** Data for statutil.c utilities. */
237 /** Loads the topology. */
238 static int load_topology(gmx_ana_traj_t *d, gmx_bool bReq);
239 /** Loads the first frame and does some checks. */
240 static int init_first_frame(gmx_ana_traj_t *d);
242 static int add_fnmarg(int nfile, t_filenm *fnm, t_filenm *fnm_add)
244 memcpy(&(fnm[nfile]), fnm_add, sizeof(*fnm_add));
248 /* Copied from src/gmxlib/statutil.c */
250 add_parg(int npargs, t_pargs *pa, t_pargs *pa_add)
252 memcpy(&(pa[npargs]), pa_add, sizeof(*pa_add));
257 * \param[out] data Trajectory analysis data structure poitner to initialize.
258 * \param[in] flags Combination of flags (see \ref analysis_flags).
259 * \returns 0 on success.
262 gmx_ana_traj_create(gmx_ana_traj_t **data, unsigned long flags)
271 d->frflags = TRX_NEED_X;
293 d->topfile_notnull = NULL;
294 rc = gmx_ana_poscalc_coll_create(&d->pcc);
301 rc = gmx_ana_selcollection_create(&d->sc, d->pcc);
304 gmx_ana_poscalc_coll_free(d->pcc);
317 * \param d Trajectory analysis data to free.
320 gmx_ana_traj_free(gmx_ana_traj_t *d)
326 sfree(d->topfile_notnull);
336 /* Gromacs does not seem to have a function for freeing frame data */
344 gmx_ana_selcollection_free(d->sc);
345 gmx_ana_poscalc_coll_free(d->pcc);
347 output_env_done(d->oenv);
352 * \param[in,out] d Trajectory analysis data structure.
353 * \param[in] flags Additional flags to set.
354 * \returns 0 on success, a non-zero error code on error.
356 * Currently there are no checks whether the flags make sense or not.
359 gmx_ana_add_flags(gmx_ana_traj_t *d, unsigned long flags)
366 * \param[in,out] d Trajectory analysis data structure.
367 * \param[in] bPBC TRUE if periodic boundary conditions should be used.
368 * \returns 0 on success.
370 * If this function is called before parse_trjana_args(), it sets the default
371 * for whether PBC are used in the analysis or not.
372 * If \ref ANA_NOUSER_PBC is not set, a command-line option is provided for the
373 * user to override the default value.
374 * If called after parse_trjana_args(), it overrides the setting provided by
375 * the user or an earlier call.
377 * If this function is not called, the default is to use PBC.
379 * If PBC are not used, the \p pbc pointer passed to gmx_analysisfunc()
382 * \see \ref ANA_NOUSER_PBC
385 gmx_ana_set_pbc(gmx_ana_traj_t *d, gmx_bool bPBC)
392 * \param[in] d Trajectory analysis data structure.
393 * \returns TRUE if periodic boundary conditions are set to be used.
396 gmx_ana_has_pbc(gmx_ana_traj_t *d)
402 * \param[in,out] d Trajectory analysis data structure.
403 * \param[in] bRmPBC TRUE if molecules should be made whole.
404 * \returns 0 on success.
406 * If this function is called before parse_trjana_args(), it sets the default
407 * for whether molecules are made whole.
408 * If \ref ANA_NOUSER_RMPBC is not set, a command-line option is provided for
409 * the user to override the default value.
410 * If called after parse_trjana_args(), it overrides the setting provided by
411 * the user or an earlier call.
413 * If this function is not called, the default is to make molecules whole.
415 * The main use of this function is to call it with FALSE if your analysis
416 * program does not require whole molecules as this can increase the
418 * In such a case, you can also specify \ref ANA_NOUSER_RMPBC to not to
419 * confuse the user with an option that would only slow the program
422 * \see \ref ANA_NOUSER_RMPBC
425 gmx_ana_set_rmpbc(gmx_ana_traj_t *d, gmx_bool bRmPBC)
432 * \param[in,out] d Trajectory analysis data structure.
433 * \param[in] frflags Flags for what to read from the trajectory file.
434 * \returns 0 on success, an error code on error.
436 * The TRX_NEED_X flag is always set.
437 * If the analysis tools needs some other information (velocities, forces),
438 * it can call this function to load additional information from the
442 gmx_ana_set_frflags(gmx_ana_traj_t *d, int frflags)
446 gmx_call("cannot set trajectory flags after initializing selections");
451 gmx_call("cannot set trajectory flags after the first frame has been read");
454 frflags |= TRX_NEED_X;
455 d->frflags = frflags;
460 * \param[in,out] d Trajectory analysis data structure.
461 * \param[in] nrefgrps Number of reference groups required.
462 * \returns 0 on success, a non-zero error code on error.
464 * \p nrefgrps should be a non-negative integer.
465 * If this function is not called (or \p nrefgrps is 0), all selections are
466 * treated as reference groups.
469 gmx_ana_set_nrefgrps(gmx_ana_traj_t *d, int nrefgrps)
474 gmx_incons("number of reference groups is negative");
477 d->nrefgrps = nrefgrps;
482 * \param[in,out] d Trajectory analysis data structure.
483 * \param[in] nanagrps Number of analysis groups required
484 * (-1 stands for any number of analysis groups).
485 * \returns 0 on success, a non-zero error code on error.
487 * \p nanagrps should be a positive integer or -1.
488 * In the latter case, any number of groups (but at least one) is acceptable.
489 * gmx_ana_get_nanagrps() can be used to access the actual value after
490 * gmx_ana_init_selections() has been called.
491 * If this function is not called, a single analysis group is expected.
494 gmx_ana_set_nanagrps(gmx_ana_traj_t *d, int nanagrps)
496 if (nanagrps <= 0 && nanagrps != -1)
499 gmx_incons("number of analysis groups is invalid");
502 d->nanagrps = nanagrps;
507 * \param[in] d Trajectory analysis data structure.
508 * \param[out] ngrps Total number of selections specified by the user.
509 * \returns 0 on success, a non-zero error code on error.
511 * The value stored in \p *ngrps is the sum of the number of reference groups
512 * and the number of analysis groups.
513 * If a specific number (not -1) of analysis groups has been set with
514 * gmx_ana_set_nanagrps(), the value is always the sum of the values provided
515 * to gmx_ana_set_nrefgrps() and gmx_ana_set_nanagrps().
517 * Should only be called after gmx_ana_init_selections().
520 gmx_ana_get_ngrps(gmx_ana_traj_t *d, int *ngrps)
522 if (d->nanagrps == -1)
525 gmx_call("gmx_ana_init_selections() not called");
528 *ngrps = d->nrefgrps + d->nanagrps;
533 * \param[in] d Trajectory analysis data structure.
534 * \param[out] nanagrps Number of analysis groups specified by the user.
535 * \returns 0 on success, a non-zero error code on error.
537 * If a specific number (not -1) of analysis groups has been set with
538 * gmx_ana_set_nanagrps(), the value is always the same value.
539 * Hence, you only need to call this function if gmx_ana_set_nanagrps() has
540 * been called with \p nanagrps set to -1.
542 * Should only be called after gmx_ana_init_selections().
545 gmx_ana_get_nanagrps(gmx_ana_traj_t *d, int *nanagrps)
547 if (d->nanagrps == -1)
550 gmx_call("gmx_ana_init_selections() not called");
553 *nanagrps = d->nanagrps;
558 * \param[in] d Trajectory analysis data structure.
559 * \param[in] i Ordinal number of the reference selection to get.
560 * \param[out] sel Selection object for the \p i'th reference group.
561 * \returns 0 on success, a non-zero error code on error.
563 * The pointer returned in \p *sel should not be freed.
564 * Should only be called after gmx_ana_init_selections().
567 gmx_ana_get_refsel(gmx_ana_traj_t *d, int i, gmx_ana_selection_t **sel)
569 if (i < 0 || i >= d->nrefgrps)
572 gmx_call("invalid reference group number");
575 *sel = gmx_ana_selcollection_get_selection(d->sc, i);
578 gmx_incons("gmx_ana_init_selections() not called");
585 * \param[in] d Trajectory analysis data structure.
586 * \param[out] sel Array of analysis selections.
587 * \returns 0 on success, a non-zero error code on error.
589 * The pointer returned in \p *sel should not be freed.
590 * Should only be called after gmx_ana_init_selections().
593 gmx_ana_get_anagrps(gmx_ana_traj_t *d, gmx_ana_selection_t ***sel)
598 gmx_incons("gmx_ana_init_selections() not called");
606 * \param[in] d Trajectory analysis data structure.
607 * \param[out] grpnames Array of selection names.
608 * \returns 0 on success, a non-zero error code on error.
610 * The pointer returned in \p *grpnames should not be freed.
611 * Should only be called after gmx_ana_init_selections().
614 gmx_ana_get_grpnames(gmx_ana_traj_t *d, char ***grpnames)
619 gmx_call("gmx_ana_init_selections() not called");
622 *grpnames = d->grpnames;
627 * \param[in] d Trajectory analysis data structure.
628 * \param[out] sc Selection collection object.
629 * \returns 0 on success.
631 * This function is mostly useful for debugging purposes.
632 * The information commonly required in analysis programs is accessible
633 * more conveniently through other means.
635 * The pointer returned in \p *sc should not be freed.
636 * Can be called at any point.
639 gmx_ana_get_selcollection(gmx_ana_traj_t *d, gmx_ana_selcollection_t **sc)
646 * \param[in,out] d Trajectory analysis data structure.
647 * \returns 0 on success, a non-zero error code on error.
649 * This function should be called first in the analysis program, much in
650 * the same way as parse_common_args() in traditional Gromacs analysis
651 * programs. It adds some command-line arguments of its own, and uses
652 * parse_common_args() to parse the command-line arguments.
653 * It also loads topology information if required or if a topology is
654 * provided on the command line.
655 * Selection handling is also initialized if it is enabled and
656 * the user has selected it on the command line.
658 * The rest of the parameters are passed on to the Gromacs routine
659 * parse_common_args(), and the use of this function should be identical
660 * to parse_common_args(), with the exception that for \p pca_flags,
661 * \p PCA_CAN_TIME and \p PCA_BE_NICE flags are automatically added.
676 parse_trjana_args(gmx_ana_traj_t *d,
677 int *argc, char *argv[], unsigned long pca_flags,
678 int nfile, t_filenm fnm[], int npargs, t_pargs *pa,
679 int ndesc, const char **desc,
680 int nbugs, const char **bugs,
683 t_filenm *all_fnm = NULL;
686 t_pargs *all_pa = NULL;
693 t_filenm def_fnm[] = {
694 {efTRX, NULL, NULL, ffOPTRD},
695 {efTPS, NULL, NULL, ffREAD},
696 {efDAT, "-sf", "selection", ffOPTRD},
697 {efNDX, NULL, NULL, ffOPTRD},
699 gmx_bool bPBC = TRUE;
701 {"-pbc", FALSE, etBOOL, {&bPBC},
702 "Use periodic boundary conditions for distance calculation"},
704 gmx_bool bRmPBC = TRUE;
705 t_pargs rmpbc_pa[] = {
706 {"-rmpbc", FALSE, etBOOL, {&bRmPBC},
707 "Make molecules whole for each frame"},
709 char *selection = NULL;
710 const char **rpost = NULL;
711 gmx_bool bSelDump = FALSE;
713 {"-select", FALSE, etSTR, {&selection},
714 "Selection string (use 'help' for help). Note that the "
715 "whole selection string will need to be quoted so that "
716 "your shell will pass it in as a string. Example: "
717 "[TT]g_select -select '\"Nearby water\" resname SOL "
718 "and within 0.25 of group Protein'[tt]"},
719 {"-seldebug", FALSE, etBOOL, {&bSelDump},
720 "HIDDENPrint out the parsed and compiled selection trees"},
722 t_pargs dsel_pa[] = {
723 {"-selrpos", FALSE, etENUM, {NULL},
724 "Selection reference position"},
726 const char **spost = NULL;
727 t_pargs selpt_pa[] = {
728 {"-seltype", FALSE, etENUM, {NULL},
729 "Default analysis positions"},
731 #define MAX_PA asize(sel_pa)+asize(dsel_pa)+5
735 gmx_incons("number of reference groups is negative");
739 if (d->flags & ANA_DEBUG_SELECTION)
744 rpost = gmx_ana_poscalc_create_type_enum(!(d->flags & ANA_REQUIRE_WHOLE));
749 spost = gmx_ana_poscalc_create_type_enum(TRUE);
756 /* Construct the file name argument array */
757 max_fnm = nfile + asize(def_fnm);
758 snew(all_fnm, max_fnm);
760 if (!(d->flags & ANA_REQUIRE_TOP))
762 def_fnm[1].flag |= ffOPT;
764 snew(fnm_map, nfile);
765 for (k = 0; k < nfile; ++k)
770 for (i = 0; i < asize(def_fnm); ++i)
772 for (k = 0; k < nfile; ++k)
774 if (fnm_map[k] == -1 && def_fnm[i].opt == NULL
775 && fnm[k].opt == NULL && fnm[k].ftp == def_fnm[i].ftp)
783 nfall = add_fnmarg(nfall, all_fnm, &(fnm[k]));
787 nfall = add_fnmarg(nfall, all_fnm, &(def_fnm[i]));
791 for (k = 0; k < nfile; ++k)
793 if (fnm_map[k] == -1)
796 nfall = add_fnmarg(nfall, all_fnm, &(fnm[k]));
800 /* Construct the argument array */
801 max_pa = npargs + MAX_PA;
802 snew(all_pa, max_pa);
805 if (!(d->flags & ANA_NOUSER_RMPBC))
807 for (i = 0; i < asize(rmpbc_pa); ++i)
809 npall = add_parg(npall, all_pa, &(rmpbc_pa[i]));
812 if (!(d->flags & ANA_NOUSER_PBC))
814 for (i = 0; i < asize(pbc_pa); ++i)
816 npall = add_parg(npall, all_pa, &(pbc_pa[i]));
820 for (i = 0; i < asize(sel_pa); ++i)
822 npall = add_parg(npall, all_pa, &(sel_pa[i]));
824 if (!(d->flags & ANA_NO_DYNSEL))
826 dsel_pa[0].u.c = rpost;
827 for (i = 0; i < asize(dsel_pa); ++i)
829 npall = add_parg(npall, all_pa, &(dsel_pa[i]));
833 if (!(d->flags & ANA_ONLY_ATOMPOS))
835 selpt_pa[0].u.c = spost;
836 for (i = 0; i < asize(selpt_pa); ++i)
838 npall = add_parg(npall, all_pa, &(selpt_pa[i]));
842 for (k = 0; k < npargs; ++k)
844 npall = add_parg(npall, all_pa, &(pa[k]));
847 pca_flags |= PCA_CAN_TIME | PCA_BE_NICE;
848 parse_common_args(argc, argv, pca_flags, nfall, all_fnm, npall, all_pa,
849 ndesc, desc, nbugs, bugs,oenv);
852 /* Process our own options.
853 * Make copies of file names for easier memory management. */
854 tmp_fnm = ftp2fn_null(efTRX, nfall, all_fnm);
855 d->trjfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
856 tmp_fnm = ftp2fn_null(efTPS, nfall, all_fnm);
857 d->topfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
858 d->topfile_notnull = strdup(ftp2fn(efTPS, nfall, all_fnm));
859 tmp_fnm = ftp2fn_null(efNDX, nfall, all_fnm);
860 d->ndxfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
861 if (!(d->flags & ANA_NOUSER_RMPBC))
865 if (!(d->flags & ANA_NOUSER_PBC))
869 d->selection = selection;
870 tmp_fnm = opt2fn_null("-sf", nfall, all_fnm);
871 d->selfile = tmp_fnm ? strdup(tmp_fnm) : NULL;
873 /* Copy the results back */
874 for (k = 0; k < nfile; ++k)
876 memcpy(&(fnm[k]), &(all_fnm[fnm_map[k]]), sizeof(fnm[k]));
877 /* Delegate responsibility of freeing the file names to caller. */
878 all_fnm[fnm_map[k]].nfiles = 0;
879 all_fnm[fnm_map[k]].fns = NULL;
881 for (i = 0, k = npall - npargs; i < (size_t)npargs; ++i, ++k)
883 memcpy(&(pa[i]), &(all_pa[k]), sizeof(pa[i]));
886 /* Free memory we have allocated. */
887 done_filenms(nfall, all_fnm);
892 if (!(d->flags & ANA_NO_DYNSEL))
894 gmx_ana_selcollection_set_refpostype(d->sc, rpost[0]);
898 gmx_ana_selcollection_set_refpostype(d->sc, rpost[1]);
903 d->flags |= ANA_DEBUG_SELECTION;
907 d->flags &= ~ANA_DEBUG_SELECTION;
910 if (!(d->flags & ANA_ONLY_ATOMPOS))
912 gmx_ana_selcollection_set_outpostype(d->sc, spost[0], d->flags & ANA_USE_POSMASK);
916 gmx_ana_selcollection_set_outpostype(d->sc, spost[1], d->flags & ANA_USE_POSMASK);
920 /* Check if the user requested help on selections.
921 * If so, call gmx_ana_init_selections() to print the help and exit. */
922 if (selection && strncmp(selection, "help", 4) == 0)
924 gmx_ana_init_selections(d);
927 /* If no trajectory file is given, we need to set some flags to be able
928 * to prepare a frame from the loaded topology information. Also, check
929 * that a topology is provided. */
934 gmx_input("No trajectory or topology provided, nothing to do!");
937 d->flags |= ANA_REQUIRE_TOP;
938 d->flags |= ANA_USE_TOPX;
941 /* Load the topology if so requested. */
942 rc = load_topology(d, (d->flags & ANA_REQUIRE_TOP));
948 /* Initialize the selections/index groups */
949 if (!(d->flags & ANA_USER_SELINIT))
951 rc = gmx_ana_init_selections(d);
958 * \param[in,out] d Trajectory analysis data structure.
959 * \param[in] bReq If TRUE, topology loading is forced.
960 * \returns 0 on success, a non-zero error code on error.
962 * Initializes the \c gmx_ana_traj_t::top, \c gmx_ana_traj_t::bTop,
963 * \c gmx_ana_traj_t::boxtop and \c gmx_ana_traj_t::ePBC fields of the
964 * analysis structure.
965 * If \p bReq is TRUE, the topology is loaded even if it is not given on
968 * The function can be called multiple times safely; extra calls are
971 static int load_topology(gmx_ana_traj_t *d, gmx_bool bReq)
975 /* Return if topology already loaded */
981 if (d->topfile || bReq)
984 d->bTop = read_tps_conf(d->topfile_notnull, title, d->top,
985 &d->ePBC, &d->xtop, NULL, d->boxtop, TRUE);
986 if (!(d->flags & ANA_USE_TOPX))
996 * \param[in] d Trajectory analysis data structure.
997 * \param[in] bReq If TRUE, topology loading is forced.
998 * \param[out] top Topology data pointer to initialize.
999 * \param[out] bTop TRUE if a full tpx file was loaded, FALSE otherwise
1000 * (can be NULL, in which case it is not used).
1001 * \returns 0 on success, a non-zero error code on error.
1003 * If \ref ANA_REQUIRE_TOP has not been specified and \p bReq is FALSE,
1004 * the pointer stored in \p *top may be NULL if no topology has been provided
1005 * on the command line.
1007 * The pointer returned in \p *top should not be freed.
1010 gmx_ana_get_topology(gmx_ana_traj_t *d, gmx_bool bReq, t_topology **top, gmx_bool *bTop)
1014 rc = load_topology(d, bReq);
1029 * \param[in] d Trajectory analysis data structure.
1030 * \param[out] x Topology data pointer to initialize.
1031 * (can be NULL, in which case it is not used).
1032 * \param[out] box Box size from the topology file
1033 * (can be NULL, in which case it is not used).
1034 * \param[out] ePBC The ePBC field from the topology
1035 * (can be NULL, in which case it is not used).
1036 * \returns 0 on success, a non-zero error code on error.
1038 * If \ref ANA_USE_TOPX has not been specified, the \p x parameter should be
1041 * The pointer returned in \p *x should not be freed.
1044 gmx_ana_get_topconf(gmx_ana_traj_t *d, rvec **x, matrix box, int *ePBC)
1048 copy_mat(d->boxtop, box);
1056 if (!(d->flags & ANA_USE_TOPX))
1058 gmx_incons("topology coordinates requested by ANA_USE_TOPX not set");
1068 * Loads default index groups from a selection file.
1070 * \param[in,out] d Trajectory analysis data structure.
1071 * \param[out] grps Pointer to receive the default groups.
1072 * \returns 0 on success, a non-zero error code on error.
1075 init_default_selections(gmx_ana_traj_t *d, gmx_ana_indexgrps_t **grps)
1077 gmx_ana_selcollection_t *sc;
1079 int nr, nr_notempty, i;
1082 /* If an index file is provided, just load it and exit. */
1085 gmx_ana_indexgrps_init(grps, d->top, d->ndxfile);
1088 /* Initialize groups to NULL if we return prematurely. */
1090 /* Return immediately if no topology provided. */
1096 /* Find the default selection file, return if none found. */
1097 fnm = low_gmxlibfn("defselection.dat", TRUE, FALSE);
1103 /* Create a temporary selection collection. */
1104 rc = gmx_ana_selcollection_create(&sc, d->pcc);
1110 rc = gmx_ana_selmethod_register_defaults(sc);
1113 gmx_ana_selcollection_free(sc);
1115 gmx_fatal(FARGS, "default selection method registration failed");
1118 /* FIXME: It would be better to not have the strings here hard-coded. */
1119 gmx_ana_selcollection_set_refpostype(sc, "atom");
1120 gmx_ana_selcollection_set_outpostype(sc, "atom", FALSE);
1122 /* Parse and compile the file with no external groups. */
1123 rc = gmx_ana_selcollection_parse_file(sc, fnm, NULL);
1127 gmx_ana_selcollection_free(sc);
1128 fprintf(stderr, "\nWARNING: default selection(s) could not be parsed\n");
1131 gmx_ana_selcollection_set_topology(sc, d->top, -1);
1132 rc = gmx_ana_selcollection_compile(sc);
1135 gmx_ana_selcollection_free(sc);
1136 fprintf(stderr, "\nWARNING: default selection(s) could not be compiled\n");
1140 /* Count the non-empty groups and check that there are no dynamic
1142 nr = gmx_ana_selcollection_get_count(sc);
1144 for (i = 0; i < nr; ++i)
1146 gmx_ana_selection_t *sel;
1148 sel = gmx_ana_selcollection_get_selection(sc, i);
1151 fprintf(stderr, "\nWARNING: dynamic default selection ignored\n");
1153 else if (sel->g->isize > 0)
1159 /* Copy the groups to the output structure */
1160 gmx_ana_indexgrps_alloc(grps, nr_notempty);
1162 for (i = 0; i < nr; ++i)
1164 gmx_ana_selection_t *sel;
1166 sel = gmx_ana_selcollection_get_selection(sc, i);
1167 if (!sel->bDynamic && sel->g->isize > 0)
1171 g = gmx_ana_indexgrps_get_grp(*grps, nr_notempty);
1172 gmx_ana_index_copy(g, sel->g, TRUE);
1173 g->name = strdup(sel->name);
1178 gmx_ana_selcollection_free(sc);
1183 * \param[in,out] d Trajectory analysis data structure.
1184 * \returns 0 on success, a non-zero error code on error.
1186 * Initializes the selection data in \c gmx_ana_traj_t based on
1187 * the selection options and/or index files provided on the command line.
1189 * This function is called automatically by parse_trjana_args() and should
1190 * not be called directly unless \ref ANA_USER_SELINIT is specified.
1192 * \see ANA_USER_SELINIT
1195 gmx_ana_init_selections(gmx_ana_traj_t *d)
1200 gmx_ana_indexgrps_t *grps;
1203 gmx_bool bInteractive;
1208 gmx_call("init_selections called more than once\n"
1209 "perhaps you forgot ANA_USER_SELINIT");
1213 gmx_ana_selcollection_set_veloutput(d->sc,
1214 d->frflags & (TRX_READ_V | TRX_NEED_V));
1215 gmx_ana_selcollection_set_forceoutput(d->sc,
1216 d->frflags & (TRX_READ_F | TRX_NEED_F));
1217 /* Check if we need some information from the topology */
1218 if (gmx_ana_selcollection_requires_top(d->sc))
1220 rc = load_topology(d, TRUE);
1226 /* Initialize the default selection methods */
1227 rc = gmx_ana_selmethod_register_defaults(d->sc);
1230 gmx_fatal(FARGS, "default selection method registration failed");
1233 /* Initialize index groups.
1234 * We ignore the return value to continue without the default groups if
1235 * there is an error there. */
1236 init_default_selections(d, &grps);
1237 /* Parse the selections */
1238 bStdIn = (d->selfile && d->selfile[0] == '-' && d->selfile[1] == 0)
1239 || (d->selection && d->selection[0] == 0)
1240 || (!d->selfile && !d->selection);
1241 /* Behavior is not very pretty if we cannot check for interactive input,
1242 * but at least it should compile and work in most cases. */
1243 #ifdef HAVE_UNISTD_H
1244 bInteractive = bStdIn && isatty(fileno(stdin));
1246 bInteractive = bStdIn;
1248 if (bStdIn && bInteractive)
1250 /* Parse from stdin */
1251 /* First we parse the reference groups if there are any */
1252 if (d->nrefgrps > 0)
1254 fprintf(stderr, "\nSpecify ");
1255 if (d->nrefgrps == 1)
1257 fprintf(stderr, "a reference selection");
1261 fprintf(stderr, "%d reference selections", d->nrefgrps);
1263 fprintf(stderr, ":\n");
1264 fprintf(stderr, "(one selection per line, 'help' for help)\n");
1265 rc = gmx_ana_selcollection_parse_stdin(d->sc, d->nrefgrps, grps, TRUE);
1266 nr = gmx_ana_selcollection_get_count(d->sc);
1267 if (rc != 0 || nr != d->nrefgrps)
1269 gmx_ana_traj_free(d);
1270 gmx_input("unrecoverable error in selection parsing");
1274 /* Then, we parse the analysis groups */
1275 fprintf(stderr, "\nSpecify ");
1276 if (d->nanagrps == 1)
1278 fprintf(stderr, "a selection");
1280 else if (d->nanagrps == -1)
1282 fprintf(stderr, "any number of selections");
1286 fprintf(stderr, "%d selections", d->nanagrps);
1288 fprintf(stderr, " for analysis:\n");
1289 fprintf(stderr, "(one selection per line, 'help' for help%s)\n",
1290 d->nanagrps == -1 ? ", Ctrl-D to end" : "");
1291 rc = gmx_ana_selcollection_parse_stdin(d->sc, d->nanagrps, grps, TRUE);
1292 fprintf(stderr, "\n");
1296 rc = gmx_ana_selcollection_parse_stdin(d->sc, -1, grps, FALSE);
1298 else if (d->selection)
1300 rc = gmx_ana_selcollection_parse_str(d->sc, d->selection, grps);
1304 rc = gmx_ana_selcollection_parse_file(d->sc, d->selfile, grps);
1308 gmx_ana_indexgrps_free(grps);
1312 /* Free memory for memory leak checking */
1313 gmx_ana_traj_free(d);
1314 gmx_input("selection(s) could not be parsed");
1318 /* Check the number of groups */
1319 nr = gmx_ana_selcollection_get_count(d->sc);
1322 /* TODO: Don't print this if the user has requested help */
1323 fprintf(stderr, "Nothing selected, finishing up.\n");
1324 gmx_ana_traj_free(d);
1325 /* TODO: It would be better to return some code that tells the caller
1326 * that one should exit. */
1329 if (nr <= d->nrefgrps)
1331 gmx_input("selection does not specify enough index groups");
1334 if (d->nanagrps <= 0)
1336 d->nanagrps = nr - d->nrefgrps;
1338 else if (nr != d->nrefgrps + d->nanagrps)
1340 gmx_input("selection does not specify the correct number of index groups");
1344 if (d->flags & ANA_DEBUG_SELECTION)
1346 gmx_ana_selcollection_print_tree(stderr, d->sc, FALSE);
1348 if (gmx_ana_selcollection_requires_top(d->sc))
1350 rc = load_topology(d, TRUE);
1362 rc = init_first_frame(d);
1367 natoms = d->fr->natoms;
1369 gmx_ana_selcollection_set_topology(d->sc, d->top, natoms);
1370 rc = gmx_ana_selcollection_compile(d->sc);
1373 /* Free memory for memory leak checking */
1374 gmx_ana_traj_free(d);
1375 gmx_input("selection could not be compiled");
1378 /* Create the selection array */
1379 d->ngrps = gmx_ana_selcollection_get_count(d->sc);
1380 if (!(d->flags & ANA_USE_FULLGRPS))
1382 d->ngrps -= d->nrefgrps;
1384 snew(d->sel, d->ngrps);
1385 for (i = 0; i < d->ngrps; ++i)
1387 if (d->flags & ANA_USE_FULLGRPS)
1389 d->sel[i] = gmx_ana_selcollection_get_selection(d->sc, i);
1393 d->sel[i] = gmx_ana_selcollection_get_selection(d->sc, i + d->nrefgrps);
1396 if (d->flags & ANA_DEBUG_SELECTION)
1398 fprintf(stderr, "\n");
1399 gmx_ana_selcollection_print_tree(stderr, d->sc, FALSE);
1400 fprintf(stderr, "\n");
1401 gmx_ana_poscalc_coll_print_tree(stderr, d->pcc);
1402 fprintf(stderr, "\n");
1405 /* Initialize the position evaluation */
1406 gmx_ana_poscalc_init_eval(d->pcc);
1407 if (d->flags & ANA_DEBUG_SELECTION)
1409 gmx_ana_poscalc_coll_print_tree(stderr, d->pcc);
1410 fprintf(stderr, "\n");
1413 /* Check that dynamic selections are not provided if not allowed */
1414 if (d->flags & ANA_NO_DYNSEL)
1416 for (i = 0; i < d->nrefgrps + d->nanagrps; ++i)
1418 gmx_ana_selection_t *sel;
1420 sel = gmx_ana_selcollection_get_selection(d->sc, i);
1423 gmx_fatal(FARGS, "%s does not support dynamic selections",
1429 /* Check that non-atom positions are not provided if not allowed.
1430 * TODO: It would be better to have these checks in the parser. */
1431 if (d->flags & ANA_ONLY_ATOMPOS)
1433 for (i = 0; i < d->nanagrps; ++i)
1435 gmx_ana_selection_t *sel;
1437 sel = gmx_ana_selcollection_get_selection(d->sc, i + d->nrefgrps);
1438 if (sel->p.m.type != INDEX_ATOM)
1440 gmx_fatal(FARGS, "%s does not support non-atom positions",
1446 /* Create the names array */
1447 snew(d->grpnames, d->ngrps);
1448 for (i = 0; i < d->ngrps; ++i)
1450 d->grpnames[i] = gmx_ana_selection_name(d->sel[i]);
1457 * \param[in,out] d Trajectory analysis data structure.
1458 * \param[in] type Type of covered fractions to calculate.
1459 * \returns 0 on success, a non-zero error code on error.
1461 * By default, covered fractions are not calculated.
1462 * If this function is called, the covered fraction calculation is
1463 * initialize to calculate the fractions of type \p type for each selection.
1464 * The function must be called after gmx_ana_init_selections() has been
1467 * For more fine-grained control of the calculation, you can use
1468 * gmx_ana_selection_init_coverfrac(): if you initialize some selections
1469 * this function to something else than CFRAC_NONE before calling
1470 * gmx_ana_init_coverfrac(), these settings are not overwritten.
1472 * You only need to call this function if your program needs to have
1473 * information about the covered fractions, e.g., for normalization.
1475 * \see gmx_ana_selection_init_coverfrac()
1478 gmx_ana_init_coverfrac(gmx_ana_traj_t *d, e_coverfrac_t type)
1482 if (type == CFRAC_NONE)
1487 for (g = 0; g < d->ngrps; ++g)
1489 if (d->sel[g]->cfractype == CFRAC_NONE)
1491 gmx_ana_selection_init_coverfrac(d->sel[g], type);
1498 * \param[in] out Output file.
1499 * \param[in] d Trajectory analysis data structure.
1500 * \returns 0 on success, a non-zero error code on error.
1502 int xvgr_selections(FILE *out, gmx_ana_traj_t *d)
1504 xvgr_selcollection(out, d->sc, d->oenv);
1509 * \param[in,out] d Trajectory analysis data structure.
1510 * \returns 0 on success, a non-zero error code on error.
1512 static int init_first_frame(gmx_ana_traj_t *d)
1516 /* Return if we have already initialized the trajectory */
1522 d->frflags |= TRX_NEED_X;
1528 if (!read_first_frame(d->oenv, &d->status, d->trjfile, d->fr, d->frflags))
1530 gmx_input("could not read coordinates from trajectory");
1534 if (d->top && d->fr->natoms > d->top->atoms.nr)
1536 gmx_fatal(FARGS, "Trajectory (%d atoms) does not match topology (%d atoms)",
1537 d->fr->natoms, d->top->atoms.nr);
1540 /* check index groups */
1541 for (i = 0; i < d->ngrps; ++i)
1543 gmx_ana_index_check(d->sel[i]->g, d->fr->natoms);
1548 /* Prepare a frame from topology information */
1549 /* TODO: Initialize more of the fields */
1550 if (d->frflags & (TRX_NEED_V))
1552 gmx_impl("Velocity reading from a topology not implemented");
1555 if (d->frflags & (TRX_NEED_F))
1557 gmx_input("Forces cannot be read from a topology");
1560 d->fr->flags = d->frflags;
1561 d->fr->natoms = d->top->atoms.nr;
1563 snew(d->fr->x, d->fr->natoms);
1564 memcpy(d->fr->x, d->xtop, sizeof(*d->fr->x)*d->fr->natoms);
1566 copy_mat(d->boxtop, d->fr->box);
1569 set_trxframe_ePBC(d->fr, d->ePBC);
1575 * \param[in,out] d Trajectory analysis data structure.
1576 * \param[out] fr First frame in the trajectory.
1577 * \returns 0 on success, a non-zero error code on error.
1579 * The pointer stored in \p *fr should not be freed by the caller.
1581 * You only need to call this function if your program needs to do some
1582 * initialization for which it requires data from the first frame.
1586 int gmx_ana_get_first_frame(gmx_ana_traj_t *d, t_trxframe **fr)
1590 rc = init_first_frame(d);
1601 * \param[in,out] d Trajectory analysis data structure.
1602 * \param[in] flags Combination of flags
1603 * (currently, there are no flags defined).
1604 * \param[in] analyze Pointer to frame analysis function.
1605 * \param data User data to be passed to \p analyze.
1606 * \returns 0 on success, a non-zero error code on error.
1608 * This function performs the actual analysis of the trajectory.
1609 * It loops through all the frames in the trajectory, and calls
1610 * \p analyze for each frame to perform the actual analysis.
1611 * Before calling \p analyze, selections are updated (if there are any).
1612 * See gmx_analysisfunc() for description of \p analyze parameters.
1614 * This function also calculates the number of frames during the run.
1616 int gmx_ana_do(gmx_ana_traj_t *d, int flags, gmx_analysisfunc analyze, void *data)
1621 gmx_rmpbc_t gpbc=NULL;
1623 rc = init_first_frame(d);
1629 ppbc = d->bPBC ? &pbc : 0;
1636 gpbc = gmx_rmpbc_init(&d->top->idef,d->ePBC,d->fr->natoms,d->fr->box);
1643 gmx_rmpbc_trxfr(gpbc,d->fr);
1647 set_pbc(&pbc, d->ePBC, d->fr->box);
1650 gmx_ana_poscalc_init_frame(d->pcc);
1651 rc = gmx_ana_selcollection_evaluate(d->sc, d->fr, ppbc);
1654 close_trj(d->status);
1655 gmx_fatal(FARGS, "selection evaluation failed");
1658 rc = analyze(d->top, d->fr, ppbc, d->ngrps, d->sel, data);
1661 close_trj(d->status);
1667 while (d->trjfile && read_next_frame(d->oenv, d->status, d->fr));
1670 gmx_rmpbc_done(gpbc);
1674 close_trj(d->status);
1675 fprintf(stderr, "Analyzed %d frames, last time %.3f\n",
1676 d->nframes, d->fr->time);
1680 fprintf(stderr, "Analyzed topology coordinates\n");
1683 /* Restore the maximal groups for dynamic selections */
1684 rc = gmx_ana_selcollection_evaluate_fin(d->sc, d->nframes);
1687 gmx_fatal(FARGS, "selection evaluation failed");
1694 * \param[in,out] d Trajectory analysis data structure.
1695 * \param[out] nframes Number of frames.
1696 * \returns 0 on success, a non-zero error code on error.
1698 * If called before gmx_ana_do(), the behavior is currently undefined.
1701 gmx_ana_get_nframes(gmx_ana_traj_t *d, int *nframes)
1703 *nframes = d->nframes;