e6905b206795a6668150021acaa8ffa4f3ee6415
[alexxy/gromacs.git] / src / programs / mdrun / tests / trajectorycomparison.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2016,2017,2018, 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.
8  *
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.
13  *
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.
18  *
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.
23  *
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.
31  *
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.
34  */
35
36 /*! \internal \file
37  * \brief Implemention of functions for comparing trajectories
38  * produced by mdrun.
39  *
40  * \author Mark Abraham <mark.j.abraham@gmail.com>
41  * \ingroup module_mdrun_integration_tests
42  */
43 #include "gmxpre.h"
44
45 #include "trajectorycomparison.h"
46
47 #include <gmock/gmock.h>
48
49 #include "gromacs/pbcutil/pbc.h"
50 #include "gromacs/trajectory/trajectoryframe.h"
51
52 #include "testutils/testasserts.h"
53 #include "testutils/testmatchers.h"
54
55 namespace gmx
56 {
57 namespace test
58 {
59
60 using ::testing::Pointwise;
61
62 /*! \brief Compares the box from \c reference and \c test
63  * according to the \c matchSettings and \c tolerance.
64  *
65  * \todo This could be streamlined when we have a proper 3D matrix
66  * class and view. */
67 static void compareBox(const TrajectoryFrame              &reference,
68                        const TrajectoryFrame              &test,
69                        const TrajectoryFrameMatchSettings &matchSettings,
70                        const FloatingPointTolerance        tolerance)
71 {
72     if (!matchSettings.mustCompareBox)
73     {
74         return;
75     }
76     bool canCompareBox = true;
77     if (!reference.hasBox())
78     {
79         ADD_FAILURE() << "Comparing the box was required, "
80         "but the reference frame did not have one";
81         canCompareBox = false;
82     }
83     if (!test.hasBox())
84     {
85         ADD_FAILURE() << "Comparing the box was required, "
86         "but the test frame did not have one";
87         canCompareBox = false;
88     }
89     if (!canCompareBox)
90     {
91         return;
92     }
93
94     // Do the comparing.
95     for (int d = 0; d < DIM; ++d)
96     {
97         for (int dd = 0; dd < DIM; ++dd)
98         {
99             EXPECT_REAL_EQ_TOL(reference.box()[d][dd], test.box()[d][dd], tolerance);
100         }
101     }
102 }
103
104 /*! \brief Help put all atom positions in \c frame into its box.
105  *
106  * This can perhaps go away when frame->x is a container. */
107 static std::vector<RVec>
108 putAtomsInBox(const TrajectoryFrame &frame)
109 {
110     std::vector<RVec> x(frame.x().begin(), frame.x().end());
111     matrix            box;
112     for (int d = 0; d < DIM; ++d)
113     {
114         for (int dd = 0; dd < DIM; ++dd)
115         {
116             box[d][dd] = frame.box()[d][dd];
117         }
118     }
119     // Note we don't need to compare bPBC because put_atoms_in_box
120     // implements a fallback if nothing specific was set in the
121     // trajectory frame.
122     put_atoms_in_box(frame.pbc(), box, x);
123     return x;
124 }
125
126 /*! \brief Compares the positions from \c reference and \c test
127  * according to the \c matchSettings and \c tolerance. */
128 static void comparePositions(const TrajectoryFrame              &reference,
129                              const TrajectoryFrame              &test,
130                              const TrajectoryFrameMatchSettings &matchSettings,
131                              const FloatingPointTolerance        tolerance)
132 {
133     bool canHandlePbc = true;
134     if (!reference.hasBox())
135     {
136         if (matchSettings.mustComparePositions)
137         {
138             ADD_FAILURE() << "Comparing positions required PBC handling, "
139             "but the reference frame did not have a box";
140         }
141         canHandlePbc = false;
142     }
143     if (!test.hasBox())
144     {
145         if (matchSettings.mustComparePositions)
146         {
147             ADD_FAILURE() << "Comparing positions required PBC handling, "
148             "but the test frame did not have a box";
149         }
150         canHandlePbc = false;
151     }
152
153     if (matchSettings.requirePbcHandling && !canHandlePbc)
154     {
155         ADD_FAILURE() << "Cannot compare positions for the above reason(s)";
156         return;
157     }
158
159     if ((matchSettings.handlePbcIfPossible || matchSettings.requirePbcHandling) && canHandlePbc)
160     {
161         EXPECT_THAT(putAtomsInBox(test), Pointwise(RVecEq(tolerance), putAtomsInBox(reference)));
162     }
163     else
164     {
165         EXPECT_THAT(test.x(), Pointwise(RVecEq(tolerance), reference.x()));
166     }
167 }
168
169 /*! \brief Compares the velocities from \c reference and \c test
170  * according to the \c matchSettings and \c tolerance. */
171 static void compareVelocities(const TrajectoryFrame              &reference,
172                               const TrajectoryFrame              &test,
173                               const TrajectoryFrameMatchSettings &matchSettings,
174                               const FloatingPointTolerance        tolerance)
175 {
176     if (!matchSettings.mustCompareVelocities)
177     {
178         return;
179     }
180     EXPECT_THAT(test.v(), Pointwise(RVecEq(tolerance), reference.v()));
181 }
182
183 /*! \brief Compares the forces from \c reference and \c test
184  * according to the \c matchSettings and \c tolerance. */
185 static void compareForces(const TrajectoryFrame              &reference,
186                           const TrajectoryFrame              &test,
187                           const TrajectoryFrameMatchSettings &matchSettings,
188                           const FloatingPointTolerance        tolerance)
189 {
190     if (!matchSettings.mustCompareForces)
191     {
192         return;
193     }
194     EXPECT_THAT(test.f(), Pointwise(RVecEq(tolerance), reference.f()));
195 }
196
197
198 void compareTrajectoryFrames(const TrajectoryFrame              &reference,
199                              const TrajectoryFrame              &test,
200                              const TrajectoryFrameMatchSettings &matchSettings,
201                              const TrajectoryTolerances         &tolerances)
202 {
203     SCOPED_TRACE("Comparing reference frame " + reference.frameName() + " and test frame " + test.frameName());
204     EXPECT_EQ(reference.step(), test.step());
205     EXPECT_EQ(reference.time(), test.time());
206     compareBox(reference, test, matchSettings, tolerances.box);
207     comparePositions(reference, test, matchSettings, tolerances.positions);
208     compareVelocities(reference, test, matchSettings, tolerances.velocities);
209     compareForces(reference, test, matchSettings, tolerances.forces);
210 }
211
212 } // namespace
213 } // namespace