2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2011,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements gmx::analysismodules::Angle.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_trajectoryanalysis
50 #include "gromacs/analysisdata/analysisdata.h"
51 #include "gromacs/analysisdata/modules/average.h"
52 #include "gromacs/analysisdata/modules/histogram.h"
53 #include "gromacs/analysisdata/modules/plot.h"
54 #include "gromacs/math/units.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/options/basicoptions.h"
57 #include "gromacs/options/filenameoption.h"
58 #include "gromacs/options/ioptionscontainer.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/selection/selection.h"
61 #include "gromacs/selection/selectionoption.h"
62 #include "gromacs/trajectory/trajectoryframe.h"
63 #include "gromacs/trajectoryanalysis/analysissettings.h"
64 #include "gromacs/utility/arrayref.h"
65 #include "gromacs/utility/exceptions.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/stringutil.h"
72 namespace analysismodules
78 /********************************************************************
83 * Helper to encapsulate logic for looping over input selections.
85 * This class provides two-dimensional iteration:
86 * - Over _angle groups_, corresponding to an input selection. If the input
87 * selection list contains a single selection, that selection gets used
88 * for all angle groups.
89 * - Within a group, over _values_, each consisting of a fixed number of
90 * selection positions. If there is only a single value within a selection,
91 * that value is returned over and over again.
92 * This transparently provides the semantics of using a single selection/vector
93 * to compute angles against multiple selections/vectors as described in the
96 * This class isn't perferctly self-contained and requires the caller to know
97 * some of the internals to use it properly, but it serves its purpose for this
98 * single analysis tool by simplifying the loops.
99 * Some methods have also been tailored to allow the caller to use it a bit
102 * \ingroup module_trajectoryanalysis
104 class AnglePositionIterator
108 * Creates an iterator to loop over input selection positions.
110 * \param[in] selections List of selections.
111 * \param[in] posCountPerValue Number of selection positions that
112 * constitute a single value for the iteration.
114 * If \p selections is empty, and/or \p posCountPerValue is zero, the
115 * iterator can still be advanced and hasValue()/hasSingleValue()
116 * called, but values cannot be accessed.
118 AnglePositionIterator(const SelectionList
&selections
,
119 int posCountPerValue
)
120 : selections_(selections
), posCountPerValue_(posCountPerValue
),
121 currentSelection_(0), nextPosition_(0)
125 //! Advances the iterator to the next group of angles.
128 if (selections_
.size() > 1)
134 //! Advances the iterator to the next angle in the current group.
137 if (!hasSingleValue())
139 nextPosition_
+= posCountPerValue_
;
144 * Returns whether this iterator represents any values.
146 * If the return value is `false`, only nextGroup(), nextValue() and
147 * hasSingleValue() are allowed to be called.
149 bool hasValue() const
151 return !selections_
.empty();
154 * Returns whether the current selection only contains a single value.
156 * Returns `false` if hasValue() returns false, which allows cutting
157 * some corners in consistency checks.
159 bool hasSingleValue() const
161 return hasValue() && currentSelection().posCount() == posCountPerValue_
;
163 //! Returns whether the current selection is dynamic.
164 bool isDynamic() const
166 return currentSelection().isDynamic();
169 * Returns whether positions in the current value are either all
170 * selected or all unselected.
172 bool allValuesConsistentlySelected() const
174 if (posCountPerValue_
<= 1)
178 const bool bSelected
= currentPosition(0).selected();
179 for (int i
= 1; i
< posCountPerValue_
; ++i
)
181 if (currentPosition(i
).selected() != bSelected
)
189 * Returns whether positions in the current value are selected.
191 * Only works reliably if allValuesConsistentlySelected() returns
194 bool currentValuesSelected() const
196 return selections_
.empty() || currentPosition(0).selected();
199 //! Returns the currently active selection.
200 const Selection
¤tSelection() const
202 GMX_ASSERT(currentSelection_
< static_cast<int>(selections_
.size()),
203 "Accessing an invalid selection");
204 return selections_
[currentSelection_
];
206 //! Returns the `i`th position for the current value.
207 SelectionPosition
currentPosition(int i
) const
209 return currentSelection().position(nextPosition_
+ i
);
212 * Extracts all coordinates corresponding to the current value.
214 * \param[out] x Array to which the positions are extracted.
216 * \p x should contain at minimum the number of positions per value
217 * passed to the constructor.
219 void getCurrentPositions(rvec x
[]) const
221 GMX_ASSERT(posCountPerValue_
> 0,
222 "Accessing positions for an invalid angle type");
223 GMX_ASSERT(nextPosition_
+ posCountPerValue_
<= currentSelection().posCount(),
224 "Accessing an invalid position");
225 for (int i
= 0; i
< posCountPerValue_
; ++i
)
227 copy_rvec(currentPosition(i
).x(), x
[i
]);
232 const SelectionList
&selections_
;
233 const int posCountPerValue_
;
234 int currentSelection_
;
237 GMX_DISALLOW_COPY_AND_ASSIGN(AnglePositionIterator
);
240 /********************************************************************
241 * Actual analysis module
244 //! How to interpret the selections in -group1.
252 //! How to interpret the selections in -group2.
260 Group2Type_SphereNormal
262 //! String values corresponding to Group1Type.
263 const char *const cGroup1TypeEnum
[] =
264 { "angle", "dihedral", "vector", "plane" };
265 //! String values corresponding to Group2Type.
266 const char *const cGroup2TypeEnum
[] =
267 { "none", "vector", "plane", "t0", "z", "sphnorm" };
269 class Angle
: public TrajectoryAnalysisModule
274 virtual void initOptions(IOptionsContainer
*options
,
275 TrajectoryAnalysisSettings
*settings
);
276 virtual void optionsFinished(TrajectoryAnalysisSettings
*settings
);
277 virtual void initAnalysis(const TrajectoryAnalysisSettings
&settings
,
278 const TopologyInformation
&top
);
280 virtual void analyzeFrame(int frnr
, const t_trxframe
&fr
, t_pbc
*pbc
,
281 TrajectoryAnalysisModuleData
*pdata
);
283 virtual void finishAnalysis(int nframes
);
284 virtual void writeOutput();
287 void initFromSelections(const SelectionList
&sel1
,
288 const SelectionList
&sel2
);
289 void checkSelections(const SelectionList
&sel1
,
290 const SelectionList
&sel2
) const;
294 SelectionOptionInfo
*sel1info_
;
295 SelectionOptionInfo
*sel2info_
;
296 std::string fnAverage_
;
298 std::string fnHistogram_
;
304 AnalysisData angles_
;
305 AnalysisDataFrameAverageModulePointer averageModule_
;
306 AnalysisDataSimpleHistogramModulePointer histogramModule_
;
308 std::vector
<int> angleCount_
;
311 std::vector
<std::vector
<RVec
> > vt0_
;
313 // Copy and assign disallowed by base.
317 : sel1info_(nullptr), sel2info_(nullptr),
318 g1type_(Group1Type_Angle
), g2type_(Group2Type_None
),
319 binWidth_(1.0), natoms1_(0), natoms2_(0)
321 averageModule_
.reset(new AnalysisDataFrameAverageModule());
322 angles_
.addModule(averageModule_
);
323 histogramModule_
.reset(new AnalysisDataSimpleHistogramModule());
324 angles_
.addModule(histogramModule_
);
326 registerAnalysisDataset(&angles_
, "angle");
327 registerBasicDataset(averageModule_
.get(), "average");
328 registerBasicDataset(&histogramModule_
->averager(), "histogram");
333 Angle::initOptions(IOptionsContainer
*options
, TrajectoryAnalysisSettings
*settings
)
335 static const char *const desc
[] = {
336 "[THISMODULE] computes different types of angles between vectors.",
337 "It supports both vectors defined by two positions and normals of",
338 "planes defined by three positions.",
339 "The z axis or the local normal of a sphere can also be used as",
340 "one of the vectors.",
341 "There are also convenience options 'angle' and 'dihedral' for",
342 "calculating bond angles and dihedrals defined by three/four",
344 "The type of the angle is specified with [TT]-g1[tt] and [TT]-g2[tt].",
345 "If [TT]-g1[tt] is [TT]angle[tt] or [TT]dihedral[tt], [TT]-g2[tt]",
346 "should not be specified.",
347 "In this case, [TT]-group1[tt] should specify one or more selections,",
348 "and each should contain triplets or quartets of positions that define",
349 "the angles to be calculated.[PAR]",
350 "If [TT]-g1[tt] is [TT]vector[tt] or [TT]plane[tt], [TT]-group1[tt]",
351 "should specify selections that contain either pairs ([TT]vector[tt])",
352 "or triplets ([TT]plane[tt]) of positions. For vectors, the positions",
353 "set the endpoints of the vector, and for planes, the three positions",
354 "are used to calculate the normal of the plane. In both cases,",
355 "[TT]-g2[tt] specifies the other vector to use (see below).[PAR]",
356 "With [TT]-g2 vector[tt] or [TT]-g2 plane[tt], [TT]-group2[tt] should",
357 "specify another set of vectors. [TT]-group1[tt] and [TT]-group2[tt]",
358 "should specify the same number of selections. It is also allowed to",
359 "only have a single selection for one of the options, in which case",
360 "the same selection is used with each selection in the other group.",
361 "Similarly, for each selection in [TT]-group1[tt], the corresponding",
362 "selection in [TT]-group2[tt] should specify the same number of",
363 "vectors or a single vector. In the latter case, the angle is",
364 "calculated between that single vector and each vector from the other",
366 "With [TT]-g2 sphnorm[tt], each selection in [TT]-group2[tt] should",
367 "specify a single position that is the center of the sphere.",
368 "The second vector is calculated as the vector from the center to the",
369 "midpoint of the positions specified by [TT]-group1[tt].[PAR]",
370 "With [TT]-g2 z[tt], [TT]-group2[tt] is not necessary, and angles",
371 "between the first vectors and the positive Z axis are calculated.[PAR]",
372 "With [TT]-g2 t0[tt], [TT]-group2[tt] is not necessary, and angles",
373 "are calculated from the vectors as they are in the first frame.[PAR]",
374 "There are three options for output:",
375 "[TT]-oav[tt] writes an xvg file with the time and the average angle",
377 "[TT]-oall[tt] writes all the individual angles.",
378 "[TT]-oh[tt] writes a histogram of the angles. The bin width can be",
379 "set with [TT]-binw[tt].",
380 "For [TT]-oav[tt] and [TT]-oh[tt], separate average/histogram is",
381 "computed for each selection in [TT]-group1[tt]."
384 settings
->setHelpText(desc
);
386 options
->addOption(FileNameOption("oav").filetype(eftPlot
).outputFile()
387 .store(&fnAverage_
).defaultBasename("angaver")
388 .description("Average angles as a function of time"));
389 options
->addOption(FileNameOption("oall").filetype(eftPlot
).outputFile()
390 .store(&fnAll_
).defaultBasename("angles")
391 .description("All angles as a function of time"));
392 options
->addOption(FileNameOption("oh").filetype(eftPlot
).outputFile()
393 .store(&fnHistogram_
).defaultBasename("anghist")
394 .description("Histogram of the angles"));
396 options
->addOption(EnumOption
<Group1Type
>("g1").enumValue(cGroup1TypeEnum
)
398 .description("Type of analysis/first vector group"));
399 options
->addOption(EnumOption
<Group2Type
>("g2").enumValue(cGroup2TypeEnum
)
401 .description("Type of second vector group"));
402 options
->addOption(DoubleOption("binw").store(&binWidth_
)
403 .description("Binwidth for -oh in degrees"));
405 sel1info_
= options
->addOption(SelectionOption("group1")
406 .required().dynamicMask().storeVector(&sel1_
)
408 .description("First analysis/vector selection"));
409 sel2info_
= options
->addOption(SelectionOption("group2")
410 .dynamicMask().storeVector(&sel2_
)
412 .description("Second analysis/vector selection"));
417 Angle::optionsFinished(TrajectoryAnalysisSettings
* /* settings */)
419 const bool bSingle
= (g1type_
== Group1Type_Angle
|| g1type_
== Group1Type_Dihedral
);
421 if (bSingle
&& g2type_
!= Group2Type_None
)
423 GMX_THROW(InconsistentInputError("Cannot use a second group (-g2) with "
424 "-g1 angle or dihedral"));
426 if (bSingle
&& sel2info_
->isSet())
428 GMX_THROW(InconsistentInputError("Cannot provide a second selection "
429 "(-group2) with -g1 angle or dihedral"));
431 if (!bSingle
&& g2type_
== Group2Type_None
)
433 GMX_THROW(InconsistentInputError("Should specify a second group (-g2) "
434 "if the first group is not an angle or a dihedral"));
437 // Set up the number of positions per angle.
440 case Group1Type_Angle
: natoms1_
= 3; break;
441 case Group1Type_Dihedral
: natoms1_
= 4; break;
442 case Group1Type_Vector
: natoms1_
= 2; break;
443 case Group1Type_Plane
: natoms1_
= 3; break;
445 GMX_THROW(InternalError("invalid -g1 value"));
449 case Group2Type_None
: natoms2_
= 0; break;
450 case Group2Type_Vector
: natoms2_
= 2; break;
451 case Group2Type_Plane
: natoms2_
= 3; break;
452 case Group2Type_TimeZero
: natoms2_
= 0; break;
453 case Group2Type_Z
: natoms2_
= 0; break;
454 case Group2Type_SphereNormal
: natoms2_
= 1; break;
456 GMX_THROW(InternalError("invalid -g2 value"));
458 if (natoms2_
== 0 && sel2info_
->isSet())
460 GMX_THROW(InconsistentInputError("Cannot provide a second selection (-group2) with -g2 t0 or z"));
462 // TODO: If bSingle is not set, the second selection option should be
468 Angle::initFromSelections(const SelectionList
&sel1
,
469 const SelectionList
&sel2
)
471 const int angleGroups
= std::max(sel1
.size(), sel2
.size());
472 const bool bHasSecondSelection
= natoms2_
> 0;
474 if (bHasSecondSelection
&& sel1
.size() != sel2
.size()
475 && std::min(sel1
.size(), sel2
.size()) != 1)
477 GMX_THROW(InconsistentInputError(
478 "-group1 and -group2 should specify the same number of selections"));
481 AnglePositionIterator
iter1(sel1
, natoms1_
);
482 AnglePositionIterator
iter2(sel2
, natoms2_
);
483 for (int g
= 0; g
< angleGroups
; ++g
, iter1
.nextGroup(), iter2
.nextGroup())
485 const int posCount1
= iter1
.currentSelection().posCount();
486 if (natoms1_
> 1 && posCount1
% natoms1_
!= 0)
488 GMX_THROW(InconsistentInputError(formatString(
489 "Number of positions in selection %d in the first group not divisible by %d",
490 static_cast<int>(g
+ 1), natoms1_
)));
492 const int angleCount1
= posCount1
/ natoms1_
;
493 int angleCount
= angleCount1
;
495 if (bHasSecondSelection
)
497 const int posCount2
= iter2
.currentSelection().posCount();
498 if (natoms2_
> 1 && posCount2
% natoms2_
!= 0)
500 GMX_THROW(InconsistentInputError(formatString(
501 "Number of positions in selection %d in the second group not divisible by %d",
502 static_cast<int>(g
+ 1), natoms2_
)));
504 if (g2type_
== Group2Type_SphereNormal
&& posCount2
!= 1)
506 GMX_THROW(InconsistentInputError(
507 "The second group should contain a single position with -g2 sphnorm"));
510 const int angleCount2
= posCount2
/ natoms2_
;
511 angleCount
= std::max(angleCount1
, angleCount2
);
512 if (angleCount1
!= angleCount2
513 && std::min(angleCount1
, angleCount2
) != 1)
515 GMX_THROW(InconsistentInputError(
516 "Number of vectors defined by the two groups are not the same"));
519 angleCount_
.push_back(angleCount
);
525 Angle::checkSelections(const SelectionList
&sel1
,
526 const SelectionList
&sel2
) const
528 AnglePositionIterator
iter1(sel1
, natoms1_
);
529 AnglePositionIterator
iter2(sel2
, natoms2_
);
530 for (size_t g
= 0; g
< angleCount_
.size(); ++g
, iter1
.nextGroup(), iter2
.nextGroup())
532 if (iter1
.isDynamic() || (iter2
.hasValue() && iter2
.isDynamic()))
534 for (int n
= 0; n
< angleCount_
[g
]; ++n
, iter1
.nextValue(), iter2
.nextValue())
537 if (!iter1
.allValuesConsistentlySelected())
541 if (!iter2
.allValuesConsistentlySelected())
545 if (angleCount_
[g
] > 1)
547 if (iter1
.hasSingleValue() && !iter1
.currentValuesSelected())
551 if (iter2
.hasSingleValue() && !iter2
.currentValuesSelected())
557 && (angleCount_
[g
] == 1
558 || (!iter1
.hasSingleValue() && !iter2
.hasSingleValue()))
559 && iter1
.currentValuesSelected() != iter2
.currentValuesSelected())
565 std::string message
=
566 formatString("Dynamic selection %d does not select "
567 "a consistent set of angles over the frames",
568 static_cast<int>(g
+ 1));
569 GMX_THROW(InconsistentInputError(message
));
578 Angle::initAnalysis(const TrajectoryAnalysisSettings
&settings
,
579 const TopologyInformation
& /* top */)
581 initFromSelections(sel1_
, sel2_
);
583 // checkSelections() ensures that both selection lists have the same size.
584 angles_
.setDataSetCount(angleCount_
.size());
585 for (size_t i
= 0; i
< angleCount_
.size(); ++i
)
587 angles_
.setColumnCount(i
, angleCount_
[i
]);
589 double histogramMin
= (g1type_
== Group1Type_Dihedral
? -180.0 : 0);
590 histogramModule_
->init(histogramFromRange(histogramMin
, 180.0)
591 .binWidth(binWidth_
).includeAll());
593 if (g2type_
== Group2Type_TimeZero
)
595 vt0_
.resize(sel1_
.size());
596 for (size_t g
= 0; g
< sel1_
.size(); ++g
)
598 vt0_
[g
].resize(sel1_
[g
].posCount() / natoms1_
);
602 if (!fnAverage_
.empty())
604 AnalysisDataPlotModulePointer
plotm(
605 new AnalysisDataPlotModule(settings
.plotSettings()));
606 plotm
->setFileName(fnAverage_
);
607 plotm
->setTitle("Average angle");
608 plotm
->setXAxisIsTime();
609 plotm
->setYLabel("Angle (degrees)");
610 // TODO: Consider adding information about the second selection,
611 // and/or a subtitle describing what kind of angle this is.
612 for (size_t g
= 0; g
< sel1_
.size(); ++g
)
614 plotm
->appendLegend(sel1_
[g
].name());
616 averageModule_
->addModule(plotm
);
621 AnalysisDataPlotModulePointer
plotm(
622 new AnalysisDataPlotModule(settings
.plotSettings()));
623 plotm
->setFileName(fnAll_
);
624 plotm
->setTitle("Angle");
625 plotm
->setXAxisIsTime();
626 plotm
->setYLabel("Angle (degrees)");
627 // TODO: Add legends? (there can be a massive amount of columns)
628 angles_
.addModule(plotm
);
631 if (!fnHistogram_
.empty())
633 AnalysisDataPlotModulePointer
plotm(
634 new AnalysisDataPlotModule(settings
.plotSettings()));
635 plotm
->setFileName(fnHistogram_
);
636 plotm
->setTitle("Angle histogram");
637 plotm
->setXLabel("Angle (degrees)");
638 plotm
->setYLabel("Probability");
639 // TODO: Consider adding information about the second selection,
640 // and/or a subtitle describing what kind of angle this is.
641 for (size_t g
= 0; g
< sel1_
.size(); ++g
)
643 plotm
->appendLegend(sel1_
[g
].name());
645 histogramModule_
->averager().addModule(plotm
);
650 //! Helper method to calculate a vector from two or three positions.
652 calc_vec(int natoms
, rvec x
[], t_pbc
*pbc
, rvec xout
, rvec cout
)
659 pbc_dx(pbc
, x
[1], x
[0], xout
);
663 rvec_sub(x
[1], x
[0], xout
);
665 svmul(0.5, xout
, cout
);
666 rvec_add(x
[0], cout
, cout
);
673 pbc_dx(pbc
, x
[1], x
[0], v1
);
674 pbc_dx(pbc
, x
[2], x
[0], v2
);
678 rvec_sub(x
[1], x
[0], v1
);
679 rvec_sub(x
[2], x
[0], v2
);
682 rvec_add(x
[0], x
[1], cout
);
683 rvec_add(cout
, x
[2], cout
);
684 svmul(1.0/3.0, cout
, cout
);
688 GMX_RELEASE_ASSERT(false, "Incorrectly initialized number of atoms");
694 Angle::analyzeFrame(int frnr
, const t_trxframe
&fr
, t_pbc
*pbc
,
695 TrajectoryAnalysisModuleData
*pdata
)
697 AnalysisDataHandle dh
= pdata
->dataHandle(angles_
);
698 const SelectionList
&sel1
= pdata
->parallelSelections(sel1_
);
699 const SelectionList
&sel2
= pdata
->parallelSelections(sel2_
);
701 checkSelections(sel1
, sel2
);
703 dh
.startFrame(frnr
, fr
.time
);
705 AnglePositionIterator
iter1(sel1
, natoms1_
);
706 AnglePositionIterator
iter2(sel2
, natoms2_
);
707 for (size_t g
= 0; g
< angleCount_
.size(); ++g
, iter1
.nextGroup(), iter2
.nextGroup())
712 // v2 & c2 are conditionally set in the switch statement below, and conditionally
713 // used in a different switch statement later. Apparently the clang static analyzer
714 // thinks there are cases where they can be used uninitialzed (which I cannot find),
715 // but to avoid trouble if we ever change just one of the switch statements it
716 // makes sense to clear them outside the first switch.
726 case Group2Type_SphereNormal
:
727 copy_rvec(sel2_
[g
].position(0).x(), c2
);
735 for (int n
= 0; n
< angleCount_
[g
]; ++n
, iter1
.nextValue(), iter2
.nextValue())
738 // x[] will be assigned below based on the number of atoms used to initialize iter1,
739 // which in turn should correspond perfectly to g1type_ (which determines how many we read),
740 // but unsurprisingly the static analyzer chokes a bit on that.
744 // checkSelections() ensures that this reflects all the involved
746 const bool bPresent
=
747 iter1
.currentValuesSelected() && iter2
.currentValuesSelected();
748 iter1
.getCurrentPositions(x
);
751 case Group1Type_Angle
:
754 pbc_dx(pbc
, x
[0], x
[1], v1
);
755 pbc_dx(pbc
, x
[2], x
[1], v2
);
759 rvec_sub(x
[0], x
[1], v1
);
760 rvec_sub(x
[2], x
[1], v2
);
762 angle
= gmx_angle(v1
, v2
);
764 case Group1Type_Dihedral
:
769 // cppcheck-suppress uninitvar
770 pbc_dx(pbc
, x
[0], x
[1], dx
[0]);
771 // cppcheck-suppress uninitvar
772 pbc_dx(pbc
, x
[2], x
[1], dx
[1]);
773 // cppcheck-suppress uninitvar
774 pbc_dx(pbc
, x
[2], x
[3], dx
[2]);
778 // cppcheck-suppress uninitvar
779 rvec_sub(x
[0], x
[1], dx
[0]);
780 // cppcheck-suppress uninitvar
781 rvec_sub(x
[2], x
[1], dx
[1]);
782 // cppcheck-suppress uninitvar
783 rvec_sub(x
[2], x
[3], dx
[2]);
785 cprod(dx
[0], dx
[1], v1
);
786 cprod(dx
[1], dx
[2], v2
);
787 angle
= gmx_angle(v1
, v2
);
788 real ipr
= iprod(dx
[0], v2
);
795 case Group1Type_Vector
:
796 case Group1Type_Plane
:
797 calc_vec(natoms1_
, x
, pbc
, v1
, c1
);
800 case Group2Type_Vector
:
801 case Group2Type_Plane
:
802 iter2
.getCurrentPositions(x
);
803 calc_vec(natoms2_
, x
, pbc
, v2
, c2
);
805 case Group2Type_TimeZero
:
806 // FIXME: This is not parallelizable.
809 copy_rvec(v1
, vt0_
[g
][n
]);
811 copy_rvec(vt0_
[g
][n
], v2
);
814 c1
[XX
] = c1
[YY
] = 0.0;
816 case Group2Type_SphereNormal
:
819 pbc_dx(pbc
, c1
, c2
, v2
);
823 rvec_sub(c1
, c2
, v2
);
827 GMX_THROW(InternalError("invalid -g2 value"));
829 angle
= gmx_angle(v1
, v2
);
832 GMX_THROW(InternalError("invalid -g1 value"));
834 dh
.setPoint(n
, angle
* RAD2DEG
, bPresent
);
842 Angle::finishAnalysis(int /*nframes*/)
844 AbstractAverageHistogram
&averageHistogram
= histogramModule_
->averager();
845 averageHistogram
.normalizeProbability();
846 averageHistogram
.done();
857 const char AngleInfo::name
[] = "gangle";
858 const char AngleInfo::shortDescription
[] =
861 TrajectoryAnalysisModulePointer
AngleInfo::create()
863 return TrajectoryAnalysisModulePointer(new Angle
);
866 } // namespace analysismodules