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.
39 * \page poscalcengine Position calculation engine
41 * The header file \ref poscalc.h defines an API for calculating positions
42 * in an automated way. This is useful mostly in the selection engine, in
43 * particular with dynamic selections, because the same COM/COG positions
44 * may be needed in several contexts. The API makes it possible to
45 * optimize the evaluation such that any heavy calculation is only done once,
46 * and the results just copied if needed more than once.
47 * The functions also provide a convenient interface for keeping the whole
48 * \c gmx_ana_pos_t structure up-to-date.
50 * A new collection of position calculations is allocated with
51 * gmx_ana_poscalc_coll_create().
52 * Calculations within one collection should share the same topology, and
53 * they are optimized. Calculations in different collections do not interact.
54 * The topology for a collection can be set with
55 * gmx_ana_poscalc_coll_set_topology().
56 * This needs to be done before calling gmx_ana_poscalc_set_maxindex() for
57 * any calculation in the collection, unless that calculation does not
58 * require topology information.
59 * All memory allocated for a collection and the calculations in it can be
60 * freed with gmx_ana_poscalc_coll_free().
62 * A new calculation is created with gmx_ana_poscalc_create().
63 * If flags need to be adjusted later, gmx_ana_poscalc_set_flags() can be
65 * After the flags are final, the largest possible index group for which the
66 * positions are needed has to be set with gmx_ana_poscalc_set_maxindex().
67 * gmx_ana_poscalc_coll_set_topology() should have been called before this
69 * After the above calls, gmx_ana_poscalc_init_pos() can be used to initialize
70 * output to a \c gmx_ana_pos_t structure. Several different structures can be
71 * initialized for the same calculation; the only requirement is that the
72 * structure passed later to gmx_ana_poscalc_update() has been initialized
74 * The memory allocated for a calculation can be freed with
75 * gmx_ana_poscalc_free().
77 * The position evaluation is simple: gmx_ana_poscalc_init_frame() should be
78 * called once for each frame, and gmx_ana_poscalc_update() can then be called
79 * for each calculation that is needed for that frame.
81 * It is also possible to initialize the calculations based on a type provided
83 * The possible strings are returned by gmx_ana_poscalc_create_type_enum(),
84 * and the string can be converted to the parameters for
85 * gmx_ana_poscalc_create() using gmx_ana_poscalc_type_from_enum().
86 * gmx_ana_poscalc_create_enum() is also provided for convenience.
89 * \brief Implementation of functions in poscalc.h.
92 * There is probably some room for optimization in the calculation of
93 * positions with bases.
94 * In particular, the current implementation may do a lot of unnecessary
96 * The interface would need to be changed to make it possible to use the
97 * same output positions for several calculations.
100 * The current algorithm for setting up base calculations could be improved
101 * in cases when there are calculations that cannot use a common base but
102 * still overlap partially (e.g., with three calculations A, B, and C
103 * such that A could use both B and C as a base, but B and C cannot use the
105 * Setting up the bases in an optimal manner in every possible situation can
106 * be quite difficult unless several bases are allowed for one calculation,
107 * but better heuristics could probably be implemented.
108 * For best results, the setup should probably be postponed (at least
109 * partially) to gmx_ana_poscalc_init_eval().
119 #include <typedefs.h>
123 #include <centerofmass.h>
124 #include <indexutil.h>
126 #include <position.h>
129 * Collection of \c gmx_ana_poscalc_t structures for the same topology.
131 * Calculations within the same structure are optimized to eliminate duplicate
134 struct gmx_ana_poscalc_coll_t
139 * Can be NULL if none of the calculations require topology data or if
140 * gmx_ana_poscalc_coll_set_topology() has not been called.
143 /** Pointer to the first data structure. */
144 gmx_ana_poscalc_t *first;
145 /** Pointer to the last data structure. */
146 gmx_ana_poscalc_t *last;
147 /** Whether the collection has been initialized for evaluation. */
152 * Data structure for position calculation.
154 struct gmx_ana_poscalc_t
157 * Type of calculation.
159 * This field may differ from the type requested by the user, because
160 * it is changed internally to the most effective calculation.
161 * For example, if the user requests a COM calculation for residues
162 * consisting of single atoms, it is simply set to POS_ATOM.
163 * To provide a consistent interface to the user, the field \p itype
164 * should be used when information should be given out.
168 * Flags for calculation options.
170 * See \ref poscalc_flags "documentation of the flags".
175 * Type for the created indices.
177 * This field always agrees with the type that the user requested, but
178 * may differ from \p type.
182 * Block data for the calculation.
186 * Mapping from the blocks to the blocks of \p sbase.
188 * If \p sbase is NULL, this field is also.
192 * Maximum evaluation group.
194 gmx_ana_index_t gmax;
196 /** Position storage for calculations that are used as a base. */
199 /** TRUE if the positions have been evaluated for the current frame. */
202 * Base position data for this calculation.
204 * If not NULL, the centers required by this calculation have
205 * already been calculated in \p sbase.
206 * The structure pointed by \p sbase is always a static calculation.
208 struct gmx_ana_poscalc_t *sbase;
209 /** Next structure in the linked list of calculations. */
210 struct gmx_ana_poscalc_t *next;
211 /** Previous structure in the linked list of calculations. */
212 struct gmx_ana_poscalc_t *prev;
213 /** Number of references to this structure. */
215 /** Collection this calculation belongs to. */
216 gmx_ana_poscalc_coll_t *coll;
219 static const char *const poscalc_enum_strings[] = {
221 "res_com", "res_cog",
222 "mol_com", "mol_cog",
223 "whole_res_com", "whole_res_cog",
224 "whole_mol_com", "whole_mol_cog",
225 "part_res_com", "part_res_cog",
226 "part_mol_com", "part_mol_cog",
227 "dyn_res_com", "dyn_res_cog",
228 "dyn_mol_com", "dyn_mol_cog",
231 #define NENUM asize(poscalc_enum_strings)
234 * Returns the partition type for a given position type.
236 * \param [in] type \c e_poscalc_t value to convert.
237 * \returns Corresponding \c e_indet_t.
240 index_type_for_poscalc(e_poscalc_t type)
244 case POS_ATOM: return INDEX_ATOM;
245 case POS_RES: return INDEX_RES;
246 case POS_MOL: return INDEX_MOL;
247 case POS_ALL: return INDEX_ALL;
248 case POS_ALL_PBC: return INDEX_ALL;
250 return INDEX_UNKNOWN;
254 * \param[in] post String (typically an enum command-line argument).
255 * Allowed values: 'atom', 'res_com', 'res_cog', 'mol_com', 'mol_cog',
256 * or one of the last four prepended by 'whole_', 'part_', or 'dyn_'.
257 * \param[out] type \c e_poscalc_t corresponding to \p post.
258 * \param[in,out] flags Flags corresponding to \p post.
259 * On input, the flags should contain the default flags.
260 * On exit, the flags \ref POS_MASS, \ref POS_COMPLMAX and
261 * \ref POS_COMPLWHOLE have been set according to \p post
262 * (the completion flags are left at the default values if no completion
264 * \returns 0 if \p post is one of the valid strings, EINVAL otherwise.
267 * Checking is not complete, and other values than those listed above
268 * may be accepted for \p post, but the results are undefined.
271 gmx_ana_poscalc_type_from_enum(const char *post, e_poscalc_t *type, int *flags)
278 *flags &= ~(POS_MASS | POS_COMPLMAX | POS_COMPLWHOLE);
282 /* Process the prefix */
286 *flags &= ~POS_COMPLMAX;
287 *flags |= POS_COMPLWHOLE;
290 else if (post[0] == 'p')
292 *flags &= ~POS_COMPLWHOLE;
293 *flags |= POS_COMPLMAX;
296 else if (post[0] == 'd')
298 *flags &= ~(POS_COMPLMAX | POS_COMPLWHOLE);
306 else if (ptr[0] == 'm')
312 gmx_incons("unknown position calculation type");
319 else if (ptr[6] == 'g')
325 gmx_incons("unknown position calculation type");
332 * \param[in] bAtom If TRUE, the "atom" value is included.
333 * \returns NULL-terminated array of strings that contains the string
334 * values acceptable for gmx_ana_poscalc_type_from_enum().
336 * The first string in the returned list is always NULL to allow the list to
337 * be used with Gromacs command-line parsing.
340 gmx_ana_poscalc_create_type_enum(gmx_bool bAtom)
347 snew(pcenum, NENUM+1);
348 for (i = 0; i < NENUM; ++i)
350 pcenum[i+1] = poscalc_enum_strings[i];
355 snew(pcenum, NENUM+1-1);
356 for (i = 1; i < NENUM; ++i)
358 pcenum[i] = poscalc_enum_strings[i];
366 * \param[out] pccp Allocated position calculation collection.
367 * \returns 0 for success.
370 gmx_ana_poscalc_coll_create(gmx_ana_poscalc_coll_t **pccp)
372 gmx_ana_poscalc_coll_t *pcc;
384 * \param[in,out] pcc Position calculation collection data structure.
385 * \param[in] top Topology data structure.
387 * This function should be called to set the topology before using
388 * gmx_ana_poscalc_set_maxindex() for any calculation that requires
389 * topology information.
392 gmx_ana_poscalc_coll_set_topology(gmx_ana_poscalc_coll_t *pcc, t_topology *top)
398 * \param[in] pcc Position calculation collection to free.
400 * The pointer \p pcc is invalid after the call.
401 * Any calculations in the collection are also freed, no matter how many
402 * references to them are left.
405 gmx_ana_poscalc_coll_free(gmx_ana_poscalc_coll_t *pcc)
409 gmx_ana_poscalc_free(pcc->first);
415 * \param[in] fp File handle to receive the output.
416 * \param[in] pcc Position calculation collection to print.
418 * The output is very technical, making this function mainly useful for
419 * debugging purposes.
422 gmx_ana_poscalc_coll_print_tree(FILE *fp, gmx_ana_poscalc_coll_t *pcc)
424 gmx_ana_poscalc_t *pc;
427 fprintf(fp, "Position calculations:\n");
432 fprintf(fp, "%2d ", i);
435 case POS_ATOM: fprintf(fp, "ATOM"); break;
436 case POS_RES: fprintf(fp, "RES"); break;
437 case POS_MOL: fprintf(fp, "MOL"); break;
438 case POS_ALL: fprintf(fp, "ALL"); break;
439 case POS_ALL_PBC: fprintf(fp, "ALL_PBC"); break;
441 if (pc->itype != index_type_for_poscalc(pc->type))
446 case INDEX_UNKNOWN: fprintf(fp, "???"); break;
447 case INDEX_ATOM: fprintf(fp, "ATOM"); break;
448 case INDEX_RES: fprintf(fp, "RES"); break;
449 case INDEX_MOL: fprintf(fp, "MOL"); break;
450 case INDEX_ALL: fprintf(fp, "ALL"); break;
454 fprintf(fp, " flg=");
455 if (pc->flags & POS_MASS)
459 if (pc->flags & POS_DYNAMIC)
463 if (pc->flags & POS_MASKONLY)
467 if (pc->flags & POS_COMPLMAX)
471 if (pc->flags & POS_COMPLWHOLE)
479 fprintf(fp, " nr=%d nra=%d", pc->b.nr, pc->b.nra);
480 fprintf(fp, " refc=%d", pc->refcount);
482 if (pc->gmax.nalloc_index > 0)
484 fprintf(fp, " Group: ");
485 if (pc->gmax.isize > 20)
487 fprintf(fp, " %d atoms", pc->gmax.isize);
491 for (j = 0; j < pc->gmax.isize; ++j)
493 fprintf(fp, " %d", pc->gmax.index[j] + 1);
498 if (pc->b.nalloc_a > 0)
500 fprintf(fp, " Atoms: ");
503 fprintf(fp, " %d atoms", pc->b.nra);
507 for (j = 0; j < pc->b.nra; ++j)
509 fprintf(fp, " %d", pc->b.a[j] + 1);
514 if (pc->b.nalloc_index > 0)
516 fprintf(fp, " Blocks:");
519 fprintf(fp, " %d pcs", pc->b.nr);
523 for (j = 0; j <= pc->b.nr; ++j)
525 fprintf(fp, " %d", pc->b.index[j]);
532 gmx_ana_poscalc_t *base;
534 fprintf(fp, " Base: ");
537 while (base && base != pc->sbase)
542 fprintf(fp, "%d", j);
543 if (pc->baseid && pc->b.nr <= 20)
546 for (j = 0; j < pc->b.nr; ++j)
548 fprintf(fp, " %d", pc->baseid[j]+1);
559 * Inserts a position calculation structure into its collection.
561 * \param pc Data structure to insert.
562 * \param before Data structure before which to insert
563 * (NULL = insert at end).
565 * Inserts \p pc to its collection before \p before.
566 * If \p before is NULL, \p pc is appended to the list.
569 insert_poscalc(gmx_ana_poscalc_t *pc, gmx_ana_poscalc_t *before)
574 pc->prev = pc->coll->last;
577 pc->coll->last->next = pc;
583 pc->prev = before->prev;
587 before->prev->next = pc;
593 pc->coll->first = pc;
598 * Removes a position calculation structure from its collection.
600 * \param pc Data structure to remove.
602 * Removes \p pc from its collection.
605 remove_poscalc(gmx_ana_poscalc_t *pc)
609 pc->prev->next = pc->next;
611 else if (pc == pc->coll->first)
613 pc->coll->first = pc->next;
617 pc->next->prev = pc->prev;
619 else if (pc == pc->coll->last)
621 pc->coll->last = pc->prev;
623 pc->prev = pc->next = NULL;
627 * Initializes position calculation using the maximum possible input index.
629 * \param[in,out] pc Position calculation data structure.
630 * \param[in] g Maximum index group for the calculation.
631 * \param[in] bBase Whether \p pc will be used as a base or not.
633 * \p bBase affects on how the \p pc->gmax field is initialized.
636 set_poscalc_maxindex(gmx_ana_poscalc_t *pc, gmx_ana_index_t *g, gmx_bool bBase)
638 gmx_ana_index_make_block(&pc->b, pc->coll->top, g, pc->itype, pc->flags & POS_COMPLWHOLE);
639 /* Set the type to POS_ATOM if the calculation in fact is such. */
640 if (pc->b.nr == pc->b.nra)
643 pc->flags &= ~(POS_MASS | POS_COMPLMAX | POS_COMPLWHOLE);
645 /* Set the POS_COMPLWHOLE flag if the calculation in fact always uses
646 * complete residues and molecules. */
647 if (!(pc->flags & POS_COMPLWHOLE)
648 && (!(pc->flags & POS_DYNAMIC) || (pc->flags & POS_COMPLMAX))
649 && (pc->type == POS_RES || pc->type == POS_MOL)
650 && gmx_ana_index_has_complete_elems(g, pc->itype, pc->coll->top))
652 pc->flags &= ~POS_COMPLMAX;
653 pc->flags |= POS_COMPLWHOLE;
655 /* Setup the gmax field */
656 if ((pc->flags & POS_COMPLWHOLE) && !bBase && pc->b.nra > g->isize)
658 gmx_ana_index_copy(&pc->gmax, g, TRUE);
659 sfree(pc->gmax.name);
660 pc->gmax.name = NULL;
664 gmx_ana_index_set(&pc->gmax, pc->b.nra, pc->b.a, NULL, 0);
669 * Checks whether a position calculation should use a base at all.
671 * \param[in] pc Position calculation data to check.
672 * \returns TRUE if \p pc can use a base and gets some benefit out of it,
676 can_use_base(gmx_ana_poscalc_t *pc)
678 /* For atoms, it should be faster to do a simple copy, so don't use a
680 if (pc->type == POS_ATOM)
684 /* For dynamic selections that do not use completion, it is not possible
686 if ((pc->type == POS_RES || pc->type == POS_MOL)
687 && (pc->flags & POS_DYNAMIC) && !(pc->flags & (POS_COMPLMAX | POS_COMPLWHOLE)))
691 /* Dynamic calculations for a single position cannot use a base. */
692 if ((pc->type == POS_ALL || pc->type == POS_ALL_PBC)
693 && (pc->flags & POS_DYNAMIC))
701 * Checks whether two position calculations should use a common base.
703 * \param[in] pc1 Calculation 1 to check for.
704 * \param[in] pc2 Calculation 2 to check for.
705 * \param[in] g1 Index group structure that contains the atoms from
707 * \param[in,out] g Working space, should have enough allocated memory to
708 * contain the intersection of the atoms in \p pc1 and \p pc2.
709 * \returns TRUE if the two calculations should be merged to use a common
710 * base, FALSE otherwise.
713 should_merge(gmx_ana_poscalc_t *pc1, gmx_ana_poscalc_t *pc2,
714 gmx_ana_index_t *g1, gmx_ana_index_t *g)
718 /* Do not merge calculations with different mass weighting. */
719 if ((pc1->flags & POS_MASS) != (pc2->flags & POS_MASS))
723 /* Avoid messing up complete calculations. */
724 if ((pc1->flags & POS_COMPLWHOLE) != (pc2->flags & POS_COMPLWHOLE))
728 /* Find the overlap between the calculations. */
729 gmx_ana_index_set(&g2, pc2->b.nra, pc2->b.a, NULL, 0);
730 gmx_ana_index_intersection(g, g1, &g2);
731 /* Do not merge if there is no overlap. */
736 /* Full completion calculations always match if the type is correct. */
737 if ((pc1->flags & POS_COMPLWHOLE) && (pc2->flags & POS_COMPLWHOLE)
738 && pc1->type == pc2->type)
742 /* The calculations also match if the intersection consists of full
744 if (gmx_ana_index_has_full_ablocks(g, &pc1->b)
745 && gmx_ana_index_has_full_ablocks(g, &pc2->b))
753 * Creates a static base for position calculation.
755 * \param pc Data structure to copy.
756 * \returns Pointer to a newly allocated base for \p pc.
758 * Creates and returns a deep copy of \p pc, but clears the
759 * \ref POS_DYNAMIC and \ref POS_MASKONLY flags.
760 * The newly created structure is set as the base (\c gmx_ana_poscalc_t::sbase)
761 * of \p pc and inserted in the collection before \p pc.
763 static gmx_ana_poscalc_t *
764 create_simple_base(gmx_ana_poscalc_t *pc)
766 gmx_ana_poscalc_t *base;
770 flags = pc->flags & ~(POS_DYNAMIC | POS_MASKONLY);
771 rc = gmx_ana_poscalc_create(&base, pc->coll, pc->type, flags);
774 gmx_fatal(FARGS, "position calculation base creation failed");
776 set_poscalc_maxindex(base, &pc->gmax, TRUE);
781 remove_poscalc(base);
782 insert_poscalc(base, pc);
788 * Merges a calculation into another calculation such that the new calculation
789 * can be used as a base.
791 * \param[in,out] base Base calculation to merge to.
792 * \param[in,out] pc Position calculation to merge to \p base.
794 * After the call, \p base can be used as a base for \p pc (or any calculation
795 * that used it as a base).
796 * It is assumed that any overlap between \p base and \p pc is in complete
797 * blocks, i.e., that the merge is possible.
800 merge_to_base(gmx_ana_poscalc_t *base, gmx_ana_poscalc_t *pc)
802 gmx_ana_index_t gp, gb, g;
804 int i, j, bi, bj, bo;
806 base->flags |= pc->flags & (POS_VELOCITIES | POS_FORCES);
807 gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, NULL, 0);
808 gmx_ana_index_set(&gb, base->b.nra, base->b.a, NULL, 0);
809 isize = gmx_ana_index_difference_size(&gp, &gb);
812 gmx_ana_index_clear(&g);
813 gmx_ana_index_reserve(&g, base->b.nra + isize);
814 /* Find the new blocks */
815 gmx_ana_index_difference(&g, &gp, &gb);
816 /* Count the blocks in g */
820 while (pc->b.a[pc->b.index[bi]] != g.index[i])
824 i += pc->b.index[bi+1] - pc->b.index[bi];
828 /* Merge the atoms into a temporary structure */
829 gmx_ana_index_merge(&g, &gb, &g);
830 /* Merge the blocks */
831 srenew(base->b.index, base->b.nr + bnr + 1);
835 bo = base->b.nr + bnr - 1;
836 base->b.index[bo+1] = i + 1;
839 if (bi < 0 || base->b.a[base->b.index[bi+1]-1] != g.index[i])
841 i -= pc->b.index[bj+1] - pc->b.index[bj];
846 if (bj >= 0 && pc->b.a[pc->b.index[bj+1]-1] == g.index[i])
850 i -= base->b.index[bi+1] - base->b.index[bi];
853 base->b.index[bo] = i + 1;
857 base->b.nalloc_index += bnr;
859 base->b.nra = g.isize;
861 base->b.nalloc_a = g.isize;
862 /* Refresh the gmax field */
863 gmx_ana_index_set(&base->gmax, base->b.nra, base->b.a, NULL, 0);
868 * Merges two bases into one.
870 * \param[in,out] tbase Base calculation to merge to.
871 * \param[in] mbase Base calculation to merge to \p tbase.
873 * After the call, \p mbase has been freed and \p tbase is used as the base
874 * for all calculations that previously had \p mbase as their base.
875 * It is assumed that any overlap between \p tbase and \p mbase is in complete
876 * blocks, i.e., that the merge is possible.
879 merge_bases(gmx_ana_poscalc_t *tbase, gmx_ana_poscalc_t *mbase)
881 gmx_ana_poscalc_t *pc;
883 merge_to_base(tbase, mbase);
884 remove_poscalc(mbase);
885 /* Set tbase as the base for all calculations that had mbase */
886 pc = tbase->coll->first;
889 if (pc->sbase == mbase)
898 gmx_ana_poscalc_free(mbase);
902 * Setups the static base calculation for a position calculation.
904 * \param[in,out] pc Position calculation to setup the base for.
907 setup_base(gmx_ana_poscalc_t *pc)
909 gmx_ana_poscalc_t *base, *pbase, *next;
910 gmx_ana_index_t gp, g;
912 /* Exit immediately if pc should not have a base. */
913 if (!can_use_base(pc))
918 gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, NULL, 0);
919 gmx_ana_index_clear(&g);
920 gmx_ana_index_reserve(&g, pc->b.nra);
922 base = pc->coll->first;
925 /* Save the next calculation so that we can safely delete base */
927 /* Skip pc, calculations that already have a base (we should match the
928 * base instead), as well as calculations that should not have a base.
929 * If the above conditions are met, check whether we should do a
932 if (base != pc && !base->sbase && can_use_base(base)
933 && should_merge(pbase, base, &gp, &g))
935 /* Check whether this is the first base found */
938 /* Create a real base if one is not present */
941 pbase = create_simple_base(base);
947 /* Make it a base for pc as well */
948 merge_to_base(pbase, pc);
952 else /* This was not the first base */
956 /* If it is not a real base, just make the new base as
957 * the base for it as well. */
958 merge_to_base(pbase, base);
964 /* If base is a real base, merge it with the new base
966 merge_bases(pbase, base);
969 gmx_ana_index_set(&gp, pbase->b.nra, pbase->b.a, NULL, 0);
970 gmx_ana_index_reserve(&g, pc->b.nra);
972 /* Proceed to the next unchecked calculation */
976 gmx_ana_index_deinit(&g);
978 /* If no base was found, create one if one is required */
979 if (!pc->sbase && (pc->flags & POS_DYNAMIC)
980 && (pc->flags & (POS_COMPLMAX | POS_COMPLWHOLE)))
982 create_simple_base(pc);
987 * \param[out] pcp Position calculation data structure pointer to initialize.
988 * \param[in,out] pcc Position calculation collection.
989 * \param[in] type Type of calculation.
990 * \param[in] flags Flags for setting calculation options
991 * (see \ref poscalc_flags "documentation of the flags").
992 * \returns 0 on success.
995 gmx_ana_poscalc_create(gmx_ana_poscalc_t **pcp, gmx_ana_poscalc_coll_t *pcc,
996 e_poscalc_t type, int flags)
998 gmx_ana_poscalc_t *pc;
1002 pc->itype = index_type_for_poscalc(type);
1003 gmx_ana_poscalc_set_flags(pc, flags);
1006 insert_poscalc(pc, NULL);
1012 * \param[out] pcp Position calculation data structure pointer to initialize.
1013 * \param[in,out] pcc Position calculation collection.
1014 * \param[in] post One of the strings acceptable for
1015 * gmx_ana_poscalc_type_from_enum().
1016 * \param[in] flags Flags for setting calculation options
1017 * (see \ref poscalc_flags "documentation of the flags").
1018 * \returns 0 on success, a non-zero error value on error.
1020 * This is a convenience wrapper for gmx_ana_poscalc_create().
1021 * \p flags sets the default calculation options if not overridden by \p post;
1022 * see gmx_ana_poscalc_type_from_enum().
1024 * \see gmx_ana_poscalc_create(), gmx_ana_poscalc_type_from_enum()
1027 gmx_ana_poscalc_create_enum(gmx_ana_poscalc_t **pcp, gmx_ana_poscalc_coll_t *pcc,
1028 const char *post, int flags)
1035 rc = gmx_ana_poscalc_type_from_enum(post, &type, &cflags);
1041 return gmx_ana_poscalc_create(pcp, pcc, type, cflags);
1045 * \param[in,out] pc Position calculation data structure.
1046 * \param[in] flags New flags.
1048 * \p flags are added to the old flags.
1049 * If calculation type is \ref POS_ATOM, \ref POS_MASS is automatically
1051 * If both \ref POS_DYNAMIC and \ref POS_MASKONLY are provided,
1052 * \ref POS_DYNAMIC is cleared.
1053 * If calculation type is not \ref POS_RES or \ref POS_MOL,
1054 * \ref POS_COMPLMAX and \ref POS_COMPLWHOLE are automatically cleared.
1057 gmx_ana_poscalc_set_flags(gmx_ana_poscalc_t *pc, int flags)
1059 if (pc->type == POS_ATOM)
1063 if (flags & POS_MASKONLY)
1065 flags &= ~POS_DYNAMIC;
1067 if (pc->type != POS_RES && pc->type != POS_MOL)
1069 flags &= ~(POS_COMPLMAX | POS_COMPLWHOLE);
1075 * \param[in,out] pc Position calculation data structure.
1076 * \param[in] g Maximum index group for the calculation.
1078 * Subsequent calls to gmx_ana_poscalc_update() should use only subsets of
1079 * \p g for evaluation.
1081 * The topology should have been set for the collection of which \p pc is
1085 gmx_ana_poscalc_set_maxindex(gmx_ana_poscalc_t *pc, gmx_ana_index_t *g)
1087 set_poscalc_maxindex(pc, g, FALSE);
1092 * \param[in] pc Position calculation data structure.
1093 * \param[out] p Output positions.
1095 * Calls to gmx_ana_poscalc_update() using \p pc should use only positions
1096 * initialized with this function.
1097 * The \c p->g pointer is initialized to point to an internal group that
1098 * contains the maximum index group set with gmx_ana_poscalc_set_maxindex().
1101 gmx_ana_poscalc_init_pos(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p)
1103 gmx_ana_indexmap_init(&p->m, &pc->gmax, pc->coll->top, pc->itype);
1104 /* Only do the static optimization when there is no completion */
1105 if (!(pc->flags & POS_DYNAMIC) && pc->b.nra == pc->gmax.isize)
1107 gmx_ana_indexmap_set_static(&p->m, &pc->b);
1109 gmx_ana_pos_reserve(p, p->m.nr, 0);
1110 if (pc->flags & POS_VELOCITIES)
1112 gmx_ana_pos_reserve_velocities(p);
1114 if (pc->flags & POS_FORCES)
1116 gmx_ana_pos_reserve_forces(p);
1118 gmx_ana_pos_set_nr(p, p->m.nr);
1119 gmx_ana_pos_set_evalgrp(p, &pc->gmax);
1123 * \param pc Position calculation data to be freed.
1125 * The \p pc pointer is invalid after the call.
1128 gmx_ana_poscalc_free(gmx_ana_poscalc_t *pc)
1136 if (pc->refcount > 0)
1142 if (pc->b.nalloc_index > 0)
1146 if (pc->b.nalloc_a > 0)
1150 if (pc->flags & POS_COMPLWHOLE)
1152 gmx_ana_index_deinit(&pc->gmax);
1156 gmx_ana_pos_free(pc->p);
1160 gmx_ana_poscalc_free(pc->sbase);
1167 * \param[in] pc Position calculation data to query.
1168 * \returns TRUE if \p pc requires topology for initialization and/or
1169 * evaluation, FALSE otherwise.
1172 gmx_ana_poscalc_requires_top(gmx_ana_poscalc_t *pc)
1174 if ((pc->flags & POS_MASS) || pc->type == POS_RES || pc->type == POS_MOL)
1182 * \param[in,out] pcc Position calculation collection to initialize.
1184 * This function does some final initialization of the data structures in the
1185 * collection to prepare them for evaluation.
1186 * After this function has been called, it is no longer possible to add new
1187 * calculations to the collection.
1189 * This function is automatically called by gmx_ana_poscalc_init_frame()
1190 * if not called by the user earlier.
1191 * Multiple calls to the function are ignored.
1194 gmx_ana_poscalc_init_eval(gmx_ana_poscalc_coll_t *pcc)
1196 gmx_ana_poscalc_t *pc;
1206 /* Initialize position storage for base calculations */
1209 gmx_ana_poscalc_init_pos(pc, pc->p);
1211 /* Construct the mapping of the base positions */
1214 snew(pc->baseid, pc->b.nr);
1215 for (bi = bj = 0; bi < pc->b.nr; ++bi, ++bj)
1217 while (pc->sbase->b.a[pc->sbase->b.index[bj]] != pc->b.a[pc->b.index[bi]])
1221 pc->baseid[bi] = bj;
1224 /* Free the block data for dynamic calculations */
1225 if (pc->flags & POS_DYNAMIC)
1227 if (pc->b.nalloc_index > 0)
1230 pc->b.nalloc_index = 0;
1232 if (pc->b.nalloc_a > 0)
1244 * \param[in,out] pcc Position calculation collection to initialize.
1246 * Clears the evaluation flag for all calculations.
1247 * Should be called for each frame before calling gmx_ana_poscalc_update().
1249 * This function is automatically called by gmx_ana_do() for each
1250 * frame, and should not be called by the user unless gmx_ana_do() is
1253 * This function calls gmx_ana_poscalc_init_eval() automatically if it has
1254 * not been called earlier.
1257 gmx_ana_poscalc_init_frame(gmx_ana_poscalc_coll_t *pcc)
1259 gmx_ana_poscalc_t *pc;
1263 gmx_ana_poscalc_init_eval(pcc);
1265 /* Clear the evaluation flags */
1275 * \param[in] pc Position calculation data.
1276 * \param[in,out] p Output positions, initialized previously with
1277 * gmx_ana_poscalc_init_pos() using \p pc.
1278 * \param[in] g Index group to use for the update.
1279 * \param[in] fr Current frame.
1280 * \param[in] pbc PBC data, or NULL if no PBC should be used.
1282 * gmx_ana_poscalc_init_frame() should be called for each frame before calling
1286 gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
1287 gmx_ana_index_t *g, t_trxframe *fr, t_pbc *pbc)
1291 if (pc->bEval == TRUE && !(pc->flags & POS_MASKONLY))
1297 gmx_ana_poscalc_update(pc->sbase, NULL, NULL, fr, pbc);
1307 gmx_ana_pos_set_evalgrp(p, g);
1309 /* Update the index map */
1310 if (pc->flags & POS_DYNAMIC)
1312 gmx_ana_indexmap_update(&p->m, g, FALSE);
1315 else if (pc->flags & POS_MASKONLY)
1317 gmx_ana_indexmap_update(&p->m, g, TRUE);
1321 if (!(pc->flags & POS_DYNAMIC))
1326 /* Evaluate the positions */
1329 /* TODO: It might be faster to evaluate the positions within this
1330 * loop instead of in the beginning. */
1331 if (pc->flags & POS_DYNAMIC)
1333 for (bi = 0; bi < p->nr; ++bi)
1335 bj = pc->baseid[p->m.refid[bi]];
1336 copy_rvec(pc->sbase->p->x[bj], p->x[bi]);
1340 for (bi = 0; bi < p->nr; ++bi)
1342 bj = pc->baseid[p->m.refid[bi]];
1343 copy_rvec(pc->sbase->p->v[bj], p->v[bi]);
1348 for (bi = 0; bi < p->nr; ++bi)
1350 bj = pc->baseid[p->m.refid[bi]];
1351 copy_rvec(pc->sbase->p->f[bj], p->f[bi]);
1357 for (bi = 0; bi < p->nr; ++bi)
1359 bj = pc->baseid[bi];
1360 copy_rvec(pc->sbase->p->x[bj], p->x[bi]);
1364 for (bi = 0; bi < p->nr; ++bi)
1366 bj = pc->baseid[bi];
1367 copy_rvec(pc->sbase->p->v[bj], p->v[bi]);
1372 for (bi = 0; bi < p->nr; ++bi)
1374 bj = pc->baseid[bi];
1375 copy_rvec(pc->sbase->p->f[bj], p->f[bi]);
1380 else /* pc->sbase is NULL */
1382 if (pc->flags & POS_DYNAMIC)
1384 pc->b.nr = p->m.mapb.nr;
1385 pc->b.index = p->m.mapb.index;
1386 pc->b.nra = g->isize;
1389 if (p->v && !fr->bV)
1391 for (i = 0; i < pc->b.nra; ++i)
1393 clear_rvec(p->v[i]);
1396 if (p->f && !fr->bF)
1398 for (i = 0; i < pc->b.nra; ++i)
1400 clear_rvec(p->f[i]);
1403 /* Here, we assume that the topology has been properly initialized,
1404 * and do not check the return values of gmx_calc_comg*(). */
1408 for (i = 0; i < pc->b.nra; ++i)
1410 copy_rvec(fr->x[pc->b.a[i]], p->x[i]);
1414 for (i = 0; i < pc->b.nra; ++i)
1416 copy_rvec(fr->v[pc->b.a[i]], p->v[i]);
1421 for (i = 0; i < pc->b.nra; ++i)
1423 copy_rvec(fr->f[pc->b.a[i]], p->f[i]);
1428 gmx_calc_comg(pc->coll->top, fr->x, pc->b.nra, pc->b.a,
1429 pc->flags & POS_MASS, p->x[0]);
1432 gmx_calc_comg(pc->coll->top, fr->v, pc->b.nra, pc->b.a,
1433 pc->flags & POS_MASS, p->v[0]);
1437 gmx_calc_comg_f(pc->coll->top, fr->f, pc->b.nra, pc->b.a,
1438 pc->flags & POS_MASS, p->f[0]);
1442 gmx_calc_comg_pbc(pc->coll->top, fr->x, pbc, pc->b.nra, pc->b.a,
1443 pc->flags & POS_MASS, p->x[0]);
1446 gmx_calc_comg(pc->coll->top, fr->v, pc->b.nra, pc->b.a,
1447 pc->flags & POS_MASS, p->v[0]);
1451 gmx_calc_comg_f(pc->coll->top, fr->f, pc->b.nra, pc->b.a,
1452 pc->flags & POS_MASS, p->f[0]);
1456 gmx_calc_comg_blocka(pc->coll->top, fr->x, &pc->b,
1457 pc->flags & POS_MASS, p->x);
1460 gmx_calc_comg_blocka(pc->coll->top, fr->v, &pc->b,
1461 pc->flags & POS_MASS, p->v);
1465 gmx_calc_comg_blocka(pc->coll->top, fr->f, &pc->b,
1466 pc->flags & POS_MASS, p->f);