Apply re-formatting to C++ in src/ tree.
[alexxy/gromacs.git] / api / nblib / tests / topology.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2020, 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  */
44 #include <gmock/gmock.h>
45 #include <gtest/gtest.h>
46
47 #include "gromacs/topology/exclusionblocks.h"
48 #include "nblib/exception.h"
49 #include "nblib/particletype.h"
50 #include "nblib/tests/testsystems.h"
51 #include "nblib/topology.h"
52
53 namespace nblib
54 {
55 namespace test
56 {
57 namespace
58 {
59
60 using ::testing::Eq;
61 using ::testing::Pointwise;
62
63 //! Compares all element between two lists of lists
64 //! Todo: unify this with the identical function in nbkernelsystem test make this a method
65 //!       of ListOfLists<>
66 template<typename T>
67 void compareLists(const gmx::ListOfLists<T>& list, const std::vector<std::vector<T>>& v)
68 {
69     ASSERT_EQ(list.size(), v.size());
70     for (std::size_t i = 0; i < list.size(); i++)
71     {
72         ASSERT_EQ(list[i].size(), v[i].size());
73         EXPECT_THAT(list[i], Pointwise(Eq(), v[i]));
74     }
75 }
76
77 // This is defined in src/gromacs/mdtypes/forcerec.h but there is also a
78 // legacy C6 macro defined there that conflicts with the nblib C6 type.
79 // Todo: Once that C6 has been refactored into a regular function, this
80 //       file can just include forcerec.h
81 //! Macro to marks particles to have Van der Waals interactions
82 #define SET_CGINFO_HAS_VDW(cgi) (cgi) = ((cgi) | (1 << 23))
83
84 TEST(NBlibTest, TopologyHasNumParticles)
85 {
86     WaterTopologyBuilder waters;
87     Topology             watersTopology = waters.buildTopology(2);
88     const int            test           = watersTopology.numParticles();
89     const int            ref            = 6;
90     EXPECT_EQ(ref, test);
91 }
92
93 TEST(NBlibTest, TopologyHasCharges)
94 {
95     WaterTopologyBuilder     waters;
96     Topology                 watersTopology = waters.buildTopology(2);
97     const std::vector<real>& test           = watersTopology.getCharges();
98     const std::vector<real>& ref = { Charges.at("Ow"), Charges.at("Hw"), Charges.at("Hw"),
99                                      Charges.at("Ow"), Charges.at("Hw"), Charges.at("Hw") };
100     EXPECT_EQ(ref, test);
101 }
102
103 TEST(NBlibTest, TopologyHasMasses)
104 {
105     WaterTopologyBuilder waters;
106     Topology             watersTopology = waters.buildTopology(2);
107
108     const Mass              refOwMass = waters.water().at("Ow").mass();
109     const Mass              refHwMass = waters.water().at("H").mass();
110     const std::vector<Mass> ref = { refOwMass, refHwMass, refHwMass, refOwMass, refHwMass, refHwMass };
111     const std::vector<Mass> test = expandQuantity(watersTopology, &ParticleType::mass);
112     EXPECT_EQ(ref, test);
113 }
114
115 TEST(NBlibTest, TopologyHasParticleTypes)
116 {
117     WaterTopologyBuilder             waters;
118     Topology                         watersTopology = waters.buildTopology(2);
119     const std::vector<ParticleType>& test           = watersTopology.getParticleTypes();
120     const ParticleType               refOw          = waters.water().at("Ow");
121     const ParticleType               refHw          = waters.water().at("H");
122     const std::vector<ParticleType>& ref            = { refOw, refHw };
123     const std::vector<ParticleType>& ref2           = { refHw, refOw };
124     EXPECT_TRUE(ref == test || ref2 == test);
125 }
126
127 TEST(NBlibTest, TopologyHasParticleTypeIds)
128 {
129     WaterTopologyBuilder waters;
130     Topology             watersTopology = waters.buildTopology(2);
131
132     const std::vector<int>&          testIds   = watersTopology.getParticleTypeIdOfAllParticles();
133     const std::vector<ParticleType>& testTypes = watersTopology.getParticleTypes();
134
135     std::vector<ParticleType> testTypesExpanded;
136     testTypesExpanded.reserve(testTypes.size());
137     for (int i : testIds)
138     {
139         testTypesExpanded.push_back(testTypes[i]);
140     }
141
142     const ParticleType              refOw = waters.water().at("Ow");
143     const ParticleType              refHw = waters.water().at("H");
144     const std::vector<ParticleType> ref   = { refOw, refHw, refHw, refOw, refHw, refHw };
145
146     EXPECT_TRUE(ref == testTypesExpanded);
147 }
148
149 TEST(NBlibTest, TopologyThrowsIdenticalParticleType)
150 {
151     //! User error: Two different ParticleTypes with the same name
152     ParticleType U235(ParticleTypeName("Uranium"), Mass(235));
153     ParticleType U238(ParticleTypeName("Uranium"), Mass(238));
154
155     Molecule ud235(MoleculeName("UraniumDimer235"));
156     ud235.addParticle(ParticleName("U1"), U235);
157     ud235.addParticle(ParticleName("U2"), U235);
158
159     Molecule ud238(MoleculeName("UraniumDimer238"));
160     ud238.addParticle(ParticleName("U1"), U238);
161     ud238.addParticle(ParticleName("U2"), U238);
162
163     TopologyBuilder topologyBuilder;
164     topologyBuilder.addMolecule(ud235, 1);
165     EXPECT_THROW(topologyBuilder.addMolecule(ud238, 1), InputException);
166 }
167
168 TEST(NBlibTest, TopologyHasExclusions)
169 {
170     WaterTopologyBuilder         waters;
171     Topology                     watersTopology = waters.buildTopology(2);
172     const gmx::ListOfLists<int>& testExclusions = watersTopology.getGmxExclusions();
173
174     const std::vector<std::vector<int>>& refExclusions = { { 0, 1, 2 }, { 0, 1, 2 }, { 0, 1, 2 },
175                                                            { 3, 4, 5 }, { 3, 4, 5 }, { 3, 4, 5 } };
176
177     compareLists(testExclusions, refExclusions);
178 }
179
180 TEST(NBlibTest, TopologyHasSequencing)
181 {
182     WaterTopologyBuilder waters;
183     Topology             watersTopology = waters.buildTopology(2);
184
185     EXPECT_EQ(0, watersTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen")));
186     EXPECT_EQ(1, watersTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1")));
187     EXPECT_EQ(2, watersTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2")));
188     EXPECT_EQ(3, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("Oxygen")));
189     EXPECT_EQ(4, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H1")));
190     EXPECT_EQ(5, watersTopology.sequenceID(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H2")));
191 }
192
193 TEST(NBlibTest, TopologyCanAggregateBonds)
194 {
195     Molecule water    = WaterMoleculeBuilder{}.waterMolecule();
196     Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
197
198     std::vector<std::tuple<Molecule, int>> molecules{ std::make_tuple(water, 2),
199                                                       std::make_tuple(methanol, 1) };
200     std::vector<HarmonicBondType>          bonds;
201     std::vector<size_t>                    bondsExpansion;
202     std::tie(bondsExpansion, bonds) = detail::collectInteractions<HarmonicBondType>(molecules);
203
204     std::vector<HarmonicBondType> bondsTest;
205     // use the expansionArray (bondsExpansion) to expand to the full list if bonds
206     std::transform(begin(bondsExpansion),
207                    end(bondsExpansion),
208                    std::back_inserter(bondsTest),
209                    [&bonds](size_t i) { return bonds[i]; });
210
211     std::vector<HarmonicBondType> waterBonds =
212             pickType<HarmonicBondType>(water.interactionData()).interactionTypes_;
213     std::vector<HarmonicBondType> methanolBonds =
214             pickType<HarmonicBondType>(methanol.interactionData()).interactionTypes_;
215
216     std::vector<HarmonicBondType> bondsReference;
217     std::copy(begin(waterBonds), end(waterBonds), std::back_inserter(bondsReference));
218     std::copy(begin(waterBonds), end(waterBonds), std::back_inserter(bondsReference));
219     std::copy(begin(methanolBonds), end(methanolBonds), std::back_inserter(bondsReference));
220
221     EXPECT_EQ(bondsTest, bondsReference);
222 }
223
224 TEST(NBlibTest, TopologyCanSequencePairIDs)
225 {
226     Molecule water    = WaterMoleculeBuilder{}.waterMolecule();
227     Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
228
229     std::vector<std::tuple<Molecule, int>> molecules{ std::make_tuple(water, 2),
230                                                       std::make_tuple(methanol, 1) };
231     detail::ParticleSequencer              particleSequencer;
232     particleSequencer.build(molecules);
233     auto pairs = detail::sequenceIDs<HarmonicBondType>(molecules, particleSequencer);
234
235     int Ow1 = particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen"));
236     int H11 = particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1"));
237     int H12 = particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2"));
238     int Ow2 = particleSequencer(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("Oxygen"));
239     int H21 = particleSequencer(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H1"));
240     int H22 = particleSequencer(MoleculeName("SOL"), 1, ResidueName("SOL"), ParticleName("H2"));
241
242     int Me  = particleSequencer(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("Me1"));
243     int MeO = particleSequencer(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("O2"));
244     int MeH = particleSequencer(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("H3"));
245
246     /// \cond DO_NOT_DOCUMENT
247 #define SORT(i, j) (i < j) ? i : j, (i < j) ? j : i
248
249     std::vector<CoordinateIndex<HarmonicBondType>> pairs_reference{
250         { SORT(Ow1, H11) }, { SORT(Ow1, H12) }, { SORT(Ow2, H21) },
251         { SORT(Ow2, H22) }, { SORT(MeO, MeH) }, { SORT(MeO, Me) }
252     };
253 #undef SORT
254     /// \endcond
255
256     EXPECT_EQ(pairs, pairs_reference);
257 }
258
259 TEST(NBlibTest, TopologySequenceIdThrows)
260 {
261     Molecule water    = WaterMoleculeBuilder{}.waterMolecule();
262     Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
263
264     std::vector<std::tuple<Molecule, int>> molecules{ std::make_tuple(water, 2),
265                                                       std::make_tuple(methanol, 1) };
266     detail::ParticleSequencer              particleSequencer;
267     particleSequencer.build(molecules);
268     auto pairs = detail::sequenceIDs<HarmonicBondType>(molecules, particleSequencer);
269
270     // Input error: no particle called O-Atom in molecule "water"
271     EXPECT_THROW(particleSequencer(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("O-Atom")),
272                  InputException);
273 }
274
275 TEST(NBlibTest, TopologyCanEliminateDuplicateBonds)
276 {
277     HarmonicBondType b1(1.0, 2.0);
278     HarmonicBondType b2(1.1, 2.1);
279     HarmonicBondType b3(1.2, 2.2);
280
281     // can be compressed to {b1,b2,b3} + {1,1,2,0,1,0,2,2}
282     std::vector<HarmonicBondType> bonds{ b2, b2, b3, b1, b2, b1, b3, b3 };
283
284     // expected output
285     std::vector<HarmonicBondType> uniqueBondsReference{ b1, b2, b3 };
286     std::vector<size_t>           indicesReference{ 1, 1, 2, 0, 1, 0, 2, 2 };
287
288     std::tuple<std::vector<size_t>, std::vector<HarmonicBondType>> bondData =
289             detail::eliminateDuplicateInteractions(bonds);
290
291     auto indices     = std::get<0>(bondData);
292     auto uniqueBonds = std::get<1>(bondData);
293
294     EXPECT_EQ(uniqueBondsReference, uniqueBonds);
295     EXPECT_EQ(indicesReference, indices);
296 }
297
298 TEST(NBlibTest, TopologyListedInteractions)
299 {
300     // Todo: add angles to SpcMethanolTopologyBuilder and extend interactions_reference below to
301     // Todo: include angles
302
303     Topology spcTopology = SpcMethanolTopologyBuilder{}.buildTopology(1, 2);
304
305     auto  interactionData = spcTopology.getInteractionData();
306     auto& harmonicBonds   = pickType<HarmonicBondType>(interactionData);
307
308     auto& indices = harmonicBonds.indices;
309     auto& bonds   = harmonicBonds.parameters;
310
311     std::map<std::tuple<int, int>, HarmonicBondType> interactions_test;
312     for (auto& ituple : indices)
313     {
314         interactions_test[std::make_tuple(std::get<0>(ituple), std::get<1>(ituple))] =
315                 bonds[std::get<2>(ituple)];
316     }
317
318     // there should be 3 unique HarmonicBondType instances
319     EXPECT_EQ(bonds.size(), 3);
320     // and 6 interaction pairs (bonds)
321     EXPECT_EQ(indices.size(), 6);
322
323     HarmonicBondType ohBond(1., 1.);
324     HarmonicBondType ohBondMethanol(1.01, 1.02);
325     HarmonicBondType ometBond(1.1, 1.2);
326
327     std::map<std::tuple<int, int>, HarmonicBondType> interactions_reference;
328
329     int Ow = spcTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen"));
330     int H1 = spcTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1"));
331     int H2 = spcTopology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2"));
332
333     int Me1 = spcTopology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("Me1"));
334     int MeO1 = spcTopology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("O2"));
335     int MeH1 = spcTopology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("H3"));
336
337     int Me2 = spcTopology.sequenceID(MoleculeName("MeOH"), 1, ResidueName("MeOH"), ParticleName("Me1"));
338     int MeO2 = spcTopology.sequenceID(MoleculeName("MeOH"), 1, ResidueName("MeOH"), ParticleName("O2"));
339     int MeH2 = spcTopology.sequenceID(MoleculeName("MeOH"), 1, ResidueName("MeOH"), ParticleName("H3"));
340
341     /// \cond DO_NOT_DOCUMENT
342 #define SORT(i, j) (i < j) ? i : j, (i < j) ? j : i
343     interactions_reference[std::make_tuple(SORT(Ow, H1))]     = ohBond;
344     interactions_reference[std::make_tuple(SORT(Ow, H2))]     = ohBond;
345     interactions_reference[std::make_tuple(SORT(MeO1, MeH1))] = ohBondMethanol;
346     interactions_reference[std::make_tuple(SORT(MeO1, Me1))]  = ometBond;
347     interactions_reference[std::make_tuple(SORT(MeO2, MeH2))] = ohBondMethanol;
348     interactions_reference[std::make_tuple(SORT(MeO2, Me2))]  = ometBond;
349 #undef SORT
350     /// \endcond
351
352     EXPECT_TRUE(std::equal(
353             begin(interactions_reference), end(interactions_reference), begin(interactions_test)));
354 }
355
356 TEST(NBlibTest, TopologyListedInteractionsMultipleTypes)
357 {
358     // Todo: add an angle type here
359
360     Molecule water    = WaterMoleculeBuilder{}.waterMolecule();
361     Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
362
363     CubicBondType testBond(1., 1., 1.);
364     DefaultAngle  testAngle(Degrees(1), 1);
365
366     water.addInteraction(ParticleName("H1"), ParticleName("H2"), testBond);
367     water.addInteraction(ParticleName("H1"), ParticleName("Oxygen"), ParticleName("H2"), testAngle);
368
369     ParticleTypesInteractions nbInteractions;
370     std::vector<std::string>  particleTypeNames = { "Ow", "H", "OMet", "CMet" };
371     for (const auto& name : particleTypeNames)
372     {
373         nbInteractions.add(ParticleTypeName(name), C6(0), C12(0));
374     }
375
376     TopologyBuilder topologyBuilder;
377     topologyBuilder.addMolecule(water, 1);
378     topologyBuilder.addMolecule(methanol, 1);
379     topologyBuilder.addParticleTypesInteractions(nbInteractions);
380
381     Topology topology = topologyBuilder.buildTopology();
382
383     auto  interactionData   = topology.getInteractionData();
384     auto& harmonicBonds     = pickType<HarmonicBondType>(interactionData);
385     auto& cubicBonds        = pickType<CubicBondType>(interactionData);
386     auto& angleInteractions = pickType<DefaultAngle>(interactionData);
387
388     HarmonicBondType              ohBond(1., 1.);
389     HarmonicBondType              ohBondMethanol(1.01, 1.02);
390     HarmonicBondType              ometBond(1.1, 1.2);
391     std::vector<HarmonicBondType> harmonicBondsReference{ ohBond, ohBondMethanol, ometBond };
392
393     EXPECT_EQ(harmonicBonds.parameters, harmonicBondsReference);
394
395     int H1 = topology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H1"));
396     int H2 = topology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("H2"));
397     int Ow = topology.sequenceID(MoleculeName("SOL"), 0, ResidueName("SOL"), ParticleName("Oxygen"));
398
399     int Me1 = topology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("Me1"));
400     int MeO1 = topology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("O2"));
401     int MeH1 = topology.sequenceID(MoleculeName("MeOH"), 0, ResidueName("MeOH"), ParticleName("H3"));
402
403     std::vector<CubicBondType>                   cubicBondsReference{ testBond };
404     std::vector<InteractionIndex<CubicBondType>> cubicIndicesReference{
405         { std::min(H1, H2), std::max(H1, H2), 0 }
406     };
407     EXPECT_EQ(cubicBondsReference, cubicBonds.parameters);
408     EXPECT_EQ(cubicIndicesReference, cubicBonds.indices);
409
410     DefaultAngle                                methanolAngle(Degrees(108.52), 397.5);
411     std::vector<DefaultAngle>                   angleReference{ testAngle, methanolAngle };
412     std::vector<InteractionIndex<DefaultAngle>> angleIndicesReference{
413         { std::min(H1, H2), Ow, std::max(H1, H2), 0 }, { std::min(MeH1, MeO1), Me1, std::max(MeO1, MeH1), 1 }
414     };
415     EXPECT_EQ(angleReference, angleInteractions.parameters);
416     EXPECT_EQ(angleIndicesReference, angleInteractions.indices);
417 }
418
419 TEST(NBlibTest, TopologyInvalidParticleInInteractionThrows)
420 {
421     Molecule water    = WaterMoleculeBuilder{}.waterMolecule();
422     Molecule methanol = MethanolMoleculeBuilder{}.methanolMolecule();
423
424     HarmonicBondType testBond(1., 1.);
425
426     // Invalid input: no particle named "Iron" in molecule water
427     water.addInteraction(ParticleName("H1"), ParticleName("Iron"), testBond);
428
429     ParticleTypesInteractions nbInteractions;
430     std::vector<std::string>  particleTypeNames = { "Ow", "H", "OMet", "CMet" };
431     for (const auto& name : particleTypeNames)
432     {
433         nbInteractions.add(ParticleTypeName(name), C6(0), C12(0));
434     }
435
436     TopologyBuilder topologyBuilder;
437     topologyBuilder.addMolecule(water, 1);
438     topologyBuilder.addMolecule(methanol, 1);
439     topologyBuilder.addParticleTypesInteractions(nbInteractions);
440
441     EXPECT_THROW(topologyBuilder.buildTopology(), InputException);
442 }
443
444
445 TEST(NBlibTest, toGmxExclusionBlockWorks)
446 {
447     std::vector<std::tuple<int, int>> testInput{ { 0, 0 }, { 0, 1 }, { 0, 2 }, { 1, 0 }, { 1, 1 },
448                                                  { 1, 2 }, { 2, 0 }, { 2, 1 }, { 2, 2 } };
449
450     std::vector<gmx::ExclusionBlock> reference;
451
452     gmx::ExclusionBlock localBlock;
453     localBlock.atomNumber.push_back(0);
454     localBlock.atomNumber.push_back(1);
455     localBlock.atomNumber.push_back(2);
456
457     reference.push_back(localBlock);
458     reference.push_back(localBlock);
459     reference.push_back(localBlock);
460
461     std::vector<gmx::ExclusionBlock> probe = detail::toGmxExclusionBlock(testInput);
462
463     ASSERT_EQ(reference.size(), probe.size());
464     for (size_t i = 0; i < reference.size(); ++i)
465     {
466         ASSERT_EQ(reference[i].atomNumber.size(), probe[i].atomNumber.size());
467         for (size_t j = 0; j < reference[i].atomNumber.size(); ++j)
468         {
469             EXPECT_EQ(reference[i].atomNumber[j], probe[i].atomNumber[j]);
470         }
471     }
472 }
473
474 } // namespace
475 } // namespace test
476 } // namespace nblib