QMCPACK
test_ReferencePoints.cpp
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2023 QMCPACK developers.
6 //
7 // File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab
8 //////////////////////////////////////////////////////////////////////////////////////
9 
10 #include "catch.hpp"
11 
12 #include "NEReferencePoints.h"
13 #include "ReferencePointsInput.h"
14 #include "ValidReferencePointsInput.h"
15 #include "OhmmsData/Libxml2Doc.h"
16 #include "EstimatorTesting.h"
17 #include "Particle/tests/MinimalParticlePool.h"
18 
19 /** \file
20  * This is a postfacto unit testing written for reference points during porting of EnergyDensity
21  * to the batched version of the estimators.
22  */
23 
24 namespace qmcplusplus
25 {
26 
27 namespace testing
28 {
29 class TestableNEReferencePoints : public NEReferencePoints
30 {
31 public:
32  TestableNEReferencePoints(const NEReferencePoints& nerp) : NEReferencePoints(nerp) {}
33  void write_testable_description(std::ostream& os) const;
34 };
35 
36 void TestableNEReferencePoints::write_testable_description(std::ostream& os) const
37 {
38  os << "{" << '\n';
39  std::map<std::string, Point>::const_iterator it, end = points_.end();
40  for (it = points_.begin(); it != end; ++it)
41  {
42  os << " {\"" << it->first << "\", {" << std::setw(16) << std::setprecision(16) << it->second[0] << ","
43  << it->second[1] << "," << it->second[2] << "}}," << '\n';
44  }
45  os << "};" << '\n';
46  return;
47 }
48 } // namespace testing
49 
50 std::ostream& operator<<(std::ostream& out, const testing::TestableNEReferencePoints& rhs)
51 {
52  rhs.write_testable_description(out);
53  return out;
54 }
55 
56 
57 std::ostream& operator<<(std::ostream& out, const testing::TestableNEReferencePoints& rhs);
58 
59 constexpr bool generate_test_data = false;
60 
61 template<typename T1, typename T2, unsigned D>
62 bool approxEquality(const TinyVector<T1, D>& val_a, const TinyVector<T2, D>& val_b)
63 {
64  for (int i = 0; i < D; ++i)
65  if (val_a[i] != Approx(val_b[i]))
66  return false;
67  return true;
68 }
69 
70 ReferencePointsInput makeTestRPI()
71 {
72  using Input = testing::ValidReferencePointsInputs;
73  Libxml2Document doc;
74  bool okay = doc.parseFromString(Input::xml[Input::valid::CELL]);
75  xmlNodePtr node = doc.getRoot();
76  return {node};
77 }
78 
79 struct PSetsAndRefList
80 {
81  ParticleSetPool ppool;
82  ParticleSet pset;
83  ParticleSet pset_ions;
84  RefVector<ParticleSet> ref_psets;
85 };
86 
87 PSetsAndRefList makePsets()
88 {
89  auto lattice = testing::makeTestLattice();
90  Communicate* comm;
91  comm = OHMMS::Controller;
92  auto particle_pool = MinimalParticlePool::make_diamondC_1x1x1(comm);
93  auto& pset = *(particle_pool.getParticleSet("e"));
94  auto& pset_ions = *(particle_pool.getParticleSet("ion"));
95 
96  // Setup particleset
97  pset.R = ParticleSet::ParticlePos{{1.751870349, 4.381521229, 2.865202269}, {3.244515371, 4.382273176, 4.21105285},
98  {3.000459944, 3.329603408, 4.265030556}, {3.748660329, 3.63420622, 5.393637791},
99  {3.033228526, 3.391869137, 4.654413566}, {3.114198787, 2.654334594, 5.231075822},
100  {3.657151589, 4.883870516, 4.201243939}, {2.97317591, 4.245644974, 4.284564732}};
101 
102  RefVector<ParticleSet> ref_psets;
103  ref_psets.push_back(pset_ions);
104  return {std::move(particle_pool), pset, pset_ions, ref_psets};
105 }
106 
107 auto expectedReferencePoints()
108 {
109  typename NEReferencePoints::Points expected_reference_points;
110  if constexpr (std::is_same_v<NEReferencePoints::Real, double>)
111  expected_reference_points = {
112  {"a1", {3.37316107749939, 3.37316107749939, 0}},
113  {"a2", {0, 3.37316107749939, 3.37316107749939}},
114  {"a3", {3.37316107749939, 0, 3.37316107749939}},
115  {"cmmm", {-3.37316107749939, -3.37316107749939, -3.37316107749939}},
116  {"cmmp", {0, -3.37316107749939, 0}},
117  {"cmpm", {-3.37316107749939, 0, 0}},
118  {"cmpp", {0, 0, 3.37316107749939}},
119  {"cpmm", {0, 0, -3.37316107749939}},
120  {"cpmp", {3.37316107749939, 0, 0}},
121  {"cppm", {0, 3.37316107749939, 0}},
122  {"cppp", {3.37316107749939, 3.37316107749939, 3.37316107749939}},
123  {"f1m", {-1.686580538749695, -1.686580538749695, 0}},
124  {"f1p", {1.686580538749695, 1.686580538749695, 0}},
125  {"f2m", {0, -1.686580538749695, -1.686580538749695}},
126  {"f2p", {0, 1.686580538749695, 1.686580538749695}},
127  {"f3m", {-1.686580538749695, 0, -1.686580538749695}},
128  {"f3p", {1.686580538749695, 0, 1.686580538749695}},
129  {"ion1", {0, 0, 0}},
130  {"ion2", {1.686580538749695, 1.686580538749695, 1.686580538749695}},
131  {"r1", {3.37316107749939, 3.37316107749939, 0}},
132  {"r2", {0, 3.37316107749939, 3.37316107749939}},
133  {"r3", {3.37316107749939, 0, 3.37316107749939}},
134  {"zero", {0, 0, 0}},
135  };
136  else
137  expected_reference_points = {
138  {"a1", {3.37316115, 3.37316115, 0}},
139  {"a2", {0, 3.37316115, 3.37316115}},
140  {"a3", {3.37316115, 0, 3.37316115}},
141  {"cmmm", {-3.37316115, -3.37316115, -3.37316115}},
142  {"cmmp", {0, -3.37316115, 0}},
143  {"cmpm", {-3.37316115, 0, 0}},
144  {"cmpp", {0, 0, 3.37316115}},
145  {"cpmm", {0, 0, -3.37316115}},
146  {"cpmp", {3.37316115, 0, 0}},
147  {"cppm", {0, 3.37316115, 0}},
148  {"cppp", {3.37316115, 3.37316115, 3.37316115}},
149  {"f1m", {-1.686580575, -1.686580575, 0}},
150  {"f1p", {1.686580575, 1.686580575, 0}},
151  {"f2m", {0, -1.686580575, -1.686580575}},
152  {"f2p", {0, 1.686580575, 1.686580575}},
153  {"f3m", {-1.686580575, 0, -1.686580575}},
154  {"f3p", {1.686580575, 0, 1.686580575}},
155  {"ion1", {0, 0, 0}},
156  {"ion2", {1.68658058, 1.68658058, 1.68658058}},
157  {"r1", {3.37316115, 3.37316115, 0}},
158  {"r2", {0, 3.37316115, 3.37316115}},
159  {"r3", {3.37316115, 0, 3.37316115}},
160  {"zero", {0, 0, 0}},
161  };
162  return expected_reference_points;
163 }
164 
165 TEST_CASE("ReferencePoints::DefaultConstruction", "[estimators]")
166 {
167  auto rpi = makeTestRPI();
168  auto ppr = makePsets();
169  NEReferencePoints ref_points(rpi, ppr.pset, ppr.ref_psets);
170 
171  if (generate_test_data)
172  {
173  testing::TestableNEReferencePoints tref_points(ref_points);
174  std::cout << "expected_reference_points" << tref_points;
175  }
176 
177  typename NEReferencePoints::Points expected_reference_points;
178  expected_reference_points = expectedReferencePoints();
179 
180  for (auto& [key, value] : ref_points.get_points())
181  {
182  bool coords_match = approxEquality(expected_reference_points[key], value);
183  CHECK(coords_match);
184  }
185 }
186 
187 TEST_CASE("ReferencePoints::Construction", "[estimators]")
188 {
189  auto rpi = makeTestRPI();
190  auto ppr = makePsets();
191  NEReferencePoints ref_points(rpi, ppr.pset, ppr.ref_psets);
192 
193  if (generate_test_data)
194  {
195  testing::TestableNEReferencePoints tref_points(ref_points);
196  std::cout << "expected_reference_points" << tref_points;
197  }
198 
199  typename NEReferencePoints::Points expected_reference_points;
200  expected_reference_points = expectedReferencePoints();
201 
202  for (auto& [key, value] : ref_points.get_points())
203  {
204  bool coords_match = approxEquality(expected_reference_points[key], value);
205  CHECK(coords_match);
206  }
207 }
208 
209 TEST_CASE("ReferencePoints::Description", "[estimators]")
210 {
211  auto rpi = makeTestRPI();
212  auto ppr = makePsets();
213  NEReferencePoints ref_points(rpi, ppr.pset, ppr.ref_psets);
214 
215  std::ostringstream ostr_stream;
216  ref_points.write_description(ostr_stream, " ");
217 
218  std::string expected_description;
219  // ParticleSet still lowers the coordinate precision in mixed
220  if constexpr (std::is_same_v<QMCTraits::RealType, double>)
221  expected_description = R"( reference_points
222  a1: 3.37316115 3.37316115 0
223  a2: 0 3.37316115 3.37316115
224  a3: 3.37316115 0 3.37316115
225  cmmm: -3.37316115 -3.37316115 -3.37316115
226  cmmp: 0 -3.37316115 0
227  cmpm: -3.37316115 0 0
228  cmpp: 0 0 3.37316115
229  cpmm: 0 0 -3.37316115
230  cpmp: 3.37316115 0 0
231  cppm: 0 3.37316115 0
232  cppp: 3.37316115 3.37316115 3.37316115
233  f1m: -1.686580575 -1.686580575 0
234  f1p: 1.686580575 1.686580575 0
235  f2m: 0 -1.686580575 -1.686580575
236  f2p: 0 1.686580575 1.686580575
237  f3m: -1.686580575 0 -1.686580575
238  f3p: 1.686580575 0 1.686580575
239  ion1: 0 0 0
240  ion2: 1.68658058 1.68658058 1.68658058
241  r1: 3.37316115 3.37316115 0
242  r2: 0 3.37316115 3.37316115
243  r3: 3.37316115 0 3.37316115
244  zero: 0 0 0
245  end reference_points
246 )";
247  else
248  expected_description = R"( reference_points
249  a1: 3.373161077 3.373161077 0
250  a2: 0 3.373161077 3.373161077
251  a3: 3.373161077 0 3.373161077
252  cmmm: -3.373161077 -3.373161077 -3.373161077
253  cmmp: 0 -3.373161077 0
254  cmpm: -3.373161077 0 0
255  cmpp: 0 0 3.373161077
256  cpmm: 0 0 -3.373161077
257  cpmp: 3.373161077 0 0
258  cppm: 0 3.373161077 0
259  cppp: 3.373161077 3.373161077 3.373161077
260  f1m: -1.686580539 -1.686580539 0
261  f1p: 1.686580539 1.686580539 0
262  f2m: 0 -1.686580539 -1.686580539
263  f2p: 0 1.686580539 1.686580539
264  f3m: -1.686580539 0 -1.686580539
265  f3p: 1.686580539 0 1.686580539
266  ion1: 0 0 0
267  ion2: 1.686580539 1.686580539 1.686580539
268  r1: 3.373161077 3.373161077 0
269  r2: 0 3.373161077 3.373161077
270  r3: 3.373161077 0 3.373161077
271  zero: 0 0 0
272  end reference_points
273 )";
274  CHECK(ostr_stream.str() == expected_description);
275  // This subclass and function are used to generate the test data and may be useful for further test cases in future.
276  testing::TestableNEReferencePoints test_ref_points(ref_points);
277  std::ostringstream ostr_testing_stream;
278  ostr_testing_stream << test_ref_points;
279  std::string expected_testable_description;
280  if constexpr (std::is_same_v<QMCTraits::RealType, double>)
281  {
282  std::cout << "ParticleSet coords were double\n";
283  expected_testable_description = R"({
284  {"a1", { 3.37316115,3.37316115,0}},
285  {"a2", { 0,3.37316115,3.37316115}},
286  {"a3", { 3.37316115,0,3.37316115}},
287  {"cmmm", { -3.37316115,-3.37316115,-3.37316115}},
288  {"cmmp", { 0,-3.37316115,0}},
289  {"cmpm", { -3.37316115,0,0}},
290  {"cmpp", { 0,0,3.37316115}},
291  {"cpmm", { 0,0,-3.37316115}},
292  {"cpmp", { 3.37316115,0,0}},
293  {"cppm", { 0,3.37316115,0}},
294  {"cppp", { 3.37316115,3.37316115,3.37316115}},
295  {"f1m", { -1.686580575,-1.686580575,0}},
296  {"f1p", { 1.686580575,1.686580575,0}},
297  {"f2m", { 0,-1.686580575,-1.686580575}},
298  {"f2p", { 0,1.686580575,1.686580575}},
299  {"f3m", { -1.686580575,0,-1.686580575}},
300  {"f3p", { 1.686580575,0,1.686580575}},
301  {"ion1", { 0,0,0}},
302  {"ion2", { 1.68658058,1.68658058,1.68658058}},
303  {"r1", { 3.37316115,3.37316115,0}},
304  {"r2", { 0,3.37316115,3.37316115}},
305  {"r3", { 3.37316115,0,3.37316115}},
306  {"zero", { 0,0,0}},
307 };
308 )";
309  }
310  else
311  {
312  std::cout << "ParticleSet coords were in float.\n";
313  expected_testable_description = R"({
314  {"a1", {3.37316107749939,3.37316107749939,0}},
315  {"a2", { 0,3.37316107749939,3.37316107749939}},
316  {"a3", {3.37316107749939,0,3.37316107749939}},
317  {"cmmm", {-3.37316107749939,-3.37316107749939,-3.37316107749939}},
318  {"cmmp", { 0,-3.37316107749939,0}},
319  {"cmpm", {-3.37316107749939,0,0}},
320  {"cmpp", { 0,0,3.37316107749939}},
321  {"cpmm", { 0,0,-3.37316107749939}},
322  {"cpmp", {3.37316107749939,0,0}},
323  {"cppm", { 0,3.37316107749939,0}},
324  {"cppp", {3.37316107749939,3.37316107749939,3.37316107749939}},
325  {"f1m", {-1.686580538749695,-1.686580538749695,0}},
326  {"f1p", {1.686580538749695,1.686580538749695,0}},
327  {"f2m", { 0,-1.686580538749695,-1.686580538749695}},
328  {"f2p", { 0,1.686580538749695,1.686580538749695}},
329  {"f3m", {-1.686580538749695,0,-1.686580538749695}},
330  {"f3p", {1.686580538749695,0,1.686580538749695}},
331  {"ion1", { 0,0,0}},
332  {"ion2", {1.686580538749695,1.686580538749695,1.686580538749695}},
333  {"r1", {3.37316107749939,3.37316107749939,0}},
334  {"r2", { 0,3.37316107749939,3.37316107749939}},
335  {"r3", {3.37316107749939,0,3.37316107749939}},
336  {"zero", { 0,0,0}},
337 };
338 )";
339  }
340  CHECK(ostr_testing_stream.str() == expected_testable_description);
341 }
342 
343 TEST_CASE("ReferencePoints::HDF5", "[estimators]")
344 {
345  auto rpi = makeTestRPI();
346  auto ppr = makePsets();
347  NEReferencePoints ref_points(rpi, ppr.pset, ppr.ref_psets);
348 
349  hdf_archive hd;
350  std::string test_file{"reference_points_test.hdf"};
351  bool okay = hd.create(test_file);
352  REQUIRE(okay);
353 
354  ref_points.write(hd);
355 
356  hd.close();
357 
358  hdf_archive hd_read;
359  bool okay_read = hd.open(test_file);
360 
361  hd.push("reference_points");
362 
363  typename NEReferencePoints::Points expected_reference_points;
364  expected_reference_points = expectedReferencePoints();
365 
366  for (auto& map_entry : expected_reference_points)
367  {
368  std::string key{map_entry.first};
369  NEReferencePoints::Point point;
370  hd.readEntry(point, key);
371  CHECK(approxEquality(point, map_entry.second));
372  }
373  hd.close();
374 }
375 
376 } // namespace qmcplusplus
qmcplusplus::TinyVector
Fixed-size array.
Definition: OhmmsTinyMeta.h:30
qmcplusplus::NEReferencePoints
This class creates, contains, and writes both user and machine readable referencepoints.
Definition: NEReferencePoints.h:37
Libxml2Document
class that handles xmlDoc
Definition: Libxml2Doc.h:76
qmcplusplus::hdf_archive::open
bool open(const std::filesystem::path &fname, unsigned flags=H5F_ACC_RDWR)
open a file
Definition: hdf_archive.cpp:155
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
qmcplusplus::NEReferencePoints::write
void write(hdf_archive &file) const
machine readable output
Definition: NEReferencePoints.cpp:102
qmcplusplus::node
if(!okay) throw std xmlNodePtr node
Definition: test_OneBodyDensityMatrices.cpp:372
qmcplusplus::makePsets
PSetsAndRefList makePsets()
Definition: test_ReferencePoints.cpp:87
qmcplusplus::NEReferencePoints::Points
std::map< std::string, Point > Points
Definition: NEReferencePoints.h:43
qmcplusplus::MinimalParticlePool::make_diamondC_1x1x1
static ParticleSetPool make_diamondC_1x1x1(Communicate *c)
Definition: MinimalParticlePool.h:190
qmcplusplus::TEST_CASE
TEST_CASE("complex_helper", "[type_traits]")
Definition: test_complex_helper.cpp:38
qmcplusplus::hdf_archive::close
void close()
close all the open groups and file
Definition: hdf_archive.cpp:38
Libxml2Document::getRoot
xmlNodePtr getRoot()
Definition: Libxml2Doc.h:88
qmcplusplus::generate_test_data
constexpr bool generate_test_data
set to true to generate new random R for particles.
Definition: test_QMCHamiltonian.cpp:32
qmcplusplus::hdf_archive
class to handle hdf file
Definition: hdf_archive.h:51
qmcplusplus::ParticleAttrib
Attaches a unit to a Vector for IO.
Definition: ParticleAttrib.h:33
OHMMS::Controller
Communicate * Controller
Global Communicator for a process.
Definition: Communicate.cpp:35
Communicate
Wrapping information on parallelism.
Definition: Communicate.h:68
qmcplusplus::lattice
CrystalLattice< OHMMS_PRECISION, OHMMS_DIM > lattice
Definition: test_SpinDensityNew.cpp:142
qmcplusplus::ParticleSet
Specialized paritlce class for atomistic simulations.
Definition: ParticleSet.h:55
qmcplusplus::PSetsAndRefList::ref_psets
RefVector< ParticleSet > ref_psets
Definition: test_ReferencePoints.cpp:84
qmcplusplus::REQUIRE
REQUIRE(std::filesystem::exists(filename))
qmcplusplus::ParticleSetPool
Manage a collection of ParticleSet objects.
Definition: ParticleSetPool.h:34
qmcplusplus::makeTestRPI
ReferencePointsInput makeTestRPI()
Definition: test_ReferencePoints.cpp:70
qmcplusplus::operator<<
std::ostream & operator<<(std::ostream &out, const AntiSymTensor< T, D > &rhs)
Definition: AntiSymTensor.h:615
qmcplusplus::approxEquality
bool approxEquality(T val_a, T val_b, std::optional< double > eps)
Definition: checkMatrix.hpp:26
qmcplusplus::hdf_archive::push
void push(const std::string &gname, bool createit=true)
push a group to the group stack
Definition: hdf_archive.cpp:198
qmcplusplus::NEReferencePoints::write_description
void write_description(std::ostream &os, const std::string &indent) const
writes a human readable representation of the reference points.
Definition: NEReferencePoints.cpp:90
qmcplusplus::RefVector
std::vector< std::reference_wrapper< T > > RefVector
Definition: template_types.hpp:32
qmcplusplus::hdf_archive::create
bool create(const std::filesystem::path &fname, unsigned flags=H5F_ACC_TRUNC)
create a file
Definition: hdf_archive.cpp:143
Libxml2Document::parseFromString
bool parseFromString(const std::string_view data)
Definition: Libxml2Doc.cpp:204
qmcplusplus::CHECK
CHECK(log_values[0]==ComplexApprox(std::complex< double >{ 5.603777579195571, -6.1586603331188225 }))
qmcplusplus::NEReferencePoints::get_points
const Points & get_points() const
return const ref to map of reference points.
Definition: NEReferencePoints.h:70
qmcplusplus::hdf_archive::readEntry
bool readEntry(T &data, const std::string &aname)
read the data from the group aname and return status use read() for inbuilt error checking ...
Definition: hdf_archive.h:293
qmcplusplus::testing::TestableNEReferencePoints::TestableNEReferencePoints
TestableNEReferencePoints(const NEReferencePoints &nerp)
Definition: test_ReferencePoints.cpp:32