022317375df13ea0027661d0987332eb1e1515d4
[alexxy/gromacs.git] / api / nblib / tests / nbkernelsystem.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2020,2021, 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 /*! \internal \file
36  * \brief
37  * This implements topology setup tests
38  *
39  * \author Victor Holanda <victor.holanda@cscs.ch>
40  * \author Joe Jordan <ejjordan@kth.se>
41  * \author Prashanth Kanduri <kanduri@cscs.ch>
42  * \author Sebastian Keller <keller@cscs.ch>
43  * \author Artem Zhmurov <zhmurov@gmail.com>
44  */
45 #include <gtest/gtest.h>
46
47 #include "gromacs/topology/exclusionblocks.h"
48 #include "nblib/forcecalculator.h"
49 #include "nblib/gmxsetup.h"
50 #include "nblib/integrator.h"
51 #include "nblib/tests/testhelpers.h"
52 #include "nblib/tests/testsystems.h"
53 #include "nblib/topology.h"
54 #include "nblib/util/setup.h"
55
56 namespace nblib
57 {
58 namespace test
59 {
60 namespace
61 {
62
63 TEST(NBlibTest, SpcMethanolForcesAreCorrect)
64 {
65     auto options        = NBKernelOptions();
66     options.nbnxmSimd   = SimdKernels::SimdNo;
67     options.coulombType = CoulombType::Cutoff;
68
69     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
70
71     auto simState        = spcMethanolSystemBuilder.setupSimulationState();
72     auto forceCalculator = ForceCalculator(simState, options);
73
74     gmx::ArrayRef<Vec3> forces(simState.forces());
75     ASSERT_NO_THROW(forceCalculator.compute(simState.coordinates(), forces));
76
77     RefDataChecker forcesOutputTest(5e-5);
78     forcesOutputTest.testArrays<Vec3>(forces, "SPC-methanol forces");
79 }
80
81 TEST(NBlibTest, ExpectedNumberOfForces)
82 {
83     auto options      = NBKernelOptions();
84     options.nbnxmSimd = SimdKernels::SimdNo;
85
86     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
87
88     auto simState        = spcMethanolSystemBuilder.setupSimulationState();
89     auto forceCalculator = ForceCalculator(simState, options);
90
91     gmx::ArrayRef<Vec3> forces(simState.forces());
92     forceCalculator.compute(simState.coordinates(), forces);
93     EXPECT_EQ(simState.topology().numParticles(), forces.size());
94 }
95
96 TEST(NBlibTest, CanIntegrateSystem)
97 {
98     auto options          = NBKernelOptions();
99     options.nbnxmSimd     = SimdKernels::SimdNo;
100     options.numIterations = 1;
101
102     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
103
104     auto simState        = spcMethanolSystemBuilder.setupSimulationState();
105     auto forceCalculator = ForceCalculator(simState, options);
106
107     LeapFrog integrator(simState.topology(), simState.box());
108
109     for (int iter = 0; iter < options.numIterations; iter++)
110     {
111         gmx::ArrayRef<Vec3> forces(simState.forces());
112         forceCalculator.compute(simState.coordinates(), simState.forces());
113         EXPECT_NO_THROW(integrator.integrate(
114                 1.0, simState.coordinates(), simState.velocities(), simState.forces()));
115     }
116 }
117
118 /*!
119  * Check if the following aspects of the ForceCalculator and
120  * LeapFrog (integrator) work as expected:
121  *
122  * 1. Calling the ForceCalculator::compute() function makes no change
123  *    to the internal representation of the system. Calling it repeatedly
124  *    without update should return the same vector of forces.
125  *
126  * 2. Once the LeapFrog objects integrates for the given time using the
127  *    forces, there the coordinates in SimulationState must change.
128  *    Calling the compute() function must now generate a new set of forces.
129  *
130  */
131 TEST(NBlibTest, UpdateChangesForces)
132 {
133     auto options          = NBKernelOptions();
134     options.nbnxmSimd     = SimdKernels::SimdNo;
135     options.numIterations = 1;
136
137     SpcMethanolSimulationStateBuilder spcMethanolSystemBuilder;
138
139     auto simState        = spcMethanolSystemBuilder.setupSimulationState();
140     auto forceCalculator = ForceCalculator(simState, options);
141
142     LeapFrog integrator(simState.topology(), simState.box());
143
144     // step 1
145     gmx::ArrayRef<Vec3> forces(simState.forces());
146     forceCalculator.compute(simState.coordinates(), simState.forces());
147
148     // copy computed forces to another array
149     std::vector<Vec3> forces_1(forces.size());
150     std::copy(forces.begin(), forces.end(), begin(forces_1));
151
152     // zero original force buffer
153     zeroCartesianArray(forces);
154
155     // check if forces change without update step
156     forceCalculator.compute(simState.coordinates(), forces);
157
158     // check if forces change without update
159     for (size_t i = 0; i < forces_1.size(); i++)
160     {
161         for (int j = 0; j < dimSize; j++)
162         {
163             EXPECT_EQ(forces[i][j], forces_1[i][j]);
164         }
165     }
166
167     // update
168     integrator.integrate(1.0, simState.coordinates(), simState.velocities(), simState.forces());
169
170     // zero original force buffer
171     zeroCartesianArray(forces);
172
173     // step 2
174     forceCalculator.compute(simState.coordinates(), forces);
175     std::vector<Vec3> forces_2(forces.size());
176     std::copy(forces.begin(), forces.end(), begin(forces_2));
177
178     // check if forces change after update
179     for (size_t i = 0; i < forces_1.size(); i++)
180     {
181         for (int j = 0; j < dimSize; j++)
182         {
183             EXPECT_NE(forces_1[i][j], forces_2[i][j]);
184         }
185     }
186 }
187
188 TEST(NBlibTest, ArgonOplsaForcesAreCorrect)
189 {
190     auto options        = NBKernelOptions();
191     options.nbnxmSimd   = SimdKernels::SimdNo;
192     options.coulombType = CoulombType::Cutoff;
193
194     ArgonSimulationStateBuilder argonSystemBuilder(fftypes::OPLSA);
195
196     auto simState        = argonSystemBuilder.setupSimulationState();
197     auto forceCalculator = ForceCalculator(simState, options);
198
199     gmx::ArrayRef<Vec3> testForces(simState.forces());
200     forceCalculator.compute(simState.coordinates(), simState.forces());
201
202     RefDataChecker forcesOutputTest(1e-7);
203     forcesOutputTest.testArrays<Vec3>(testForces, "Argon forces");
204 }
205
206 TEST(NBlibTest, ArgonGromos43A1ForcesAreCorrect)
207 {
208     auto options        = NBKernelOptions();
209     options.nbnxmSimd   = SimdKernels::SimdNo;
210     options.coulombType = CoulombType::Cutoff;
211
212     ArgonSimulationStateBuilder argonSystemBuilder(fftypes::GROMOS43A1);
213
214     auto simState        = argonSystemBuilder.setupSimulationState();
215     auto forceCalculator = ForceCalculator(simState, options);
216
217     gmx::ArrayRef<Vec3> testForces(simState.forces());
218     forceCalculator.compute(simState.coordinates(), simState.forces());
219
220     RefDataChecker forcesOutputTest(1e-7);
221     forcesOutputTest.testArrays<Vec3>(testForces, "Argon forces");
222 }
223
224 } // namespace
225 } // namespace test
226 } // namespace nblib