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) 2020 QMCPACK developers.
6 //
7 // File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab
8 //
9 // File created by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab
10 //////////////////////////////////////////////////////////////////////////////////////
11
12
13 #include "catch.hpp"
14
15 #include "SpinDensityInput.h"
16 #include "ValidSpinDensityInput.h"
17 #include "SpinDensityNew.h"
18 #include "RandomForTest.h"
19 #include "Utilities/SpeciesSet.h"
20 #include "ParticleSet.h"
21 #include "SpinDensityTesting.h"
22
23 #include "OhmmsData/Libxml2Doc.h"
24
25 #include <stdio.h>
26 #include <sstream>
27
28 namespace qmcplusplus
29 {
30
31 using QMCT = QMCTraits;
32
accumulateFromPsets(int ncrowds,SpinDensityNew & sdn,UPtrVector<OperatorEstBase> & crowd_sdns)33 void accumulateFromPsets(int ncrowds, SpinDensityNew& sdn, UPtrVector<OperatorEstBase>& crowd_sdns)
34 {
35 for (int iops = 0; iops < ncrowds; ++iops)
36 {
37 std::vector<OperatorEstBase::MCPWalker> walkers;
38 int nwalkers = 4;
39 for (int iw = 0; iw < nwalkers; ++iw)
40 walkers.emplace_back(2);
41
42 std::vector<ParticleSet> psets;
43
44 crowd_sdns.emplace_back(std::make_unique<SpinDensityNew>(sdn));
45 SpinDensityNew& crowd_sdn = dynamic_cast<SpinDensityNew&>(*(crowd_sdns.back()));
46
47 for (int iw = 0; iw < nwalkers; ++iw)
48 {
49 psets.emplace_back();
50 ParticleSet& pset = psets.back();
51 pset.create(2);
52 pset.R[0] = ParticleSet::PosType(0.00000000, 0.00000000, 0.00000000);
53 pset.R[1] = ParticleSet::PosType(0.68658058, 0.68658058, 0.68658058);
54 }
55
56 auto ref_walkers = makeRefVector<OperatorEstBase::MCPWalker>(walkers);
57 auto ref_psets = makeRefVector<ParticleSet>(psets);
58
59 crowd_sdn.accumulate(ref_walkers, ref_psets);
60 }
61 }
62
randomUpdateAccumulate(testing::RandomForTest<QMCT::RealType> & rft,UPtrVector<OperatorEstBase> & crowd_sdns)63 void randomUpdateAccumulate(testing::RandomForTest<QMCT::RealType>& rft, UPtrVector<OperatorEstBase>& crowd_sdns)
64 {
65 for (auto& uptr_crowd_sdn : crowd_sdns)
66 {
67 std::vector<OperatorEstBase::MCPWalker> walkers;
68 int nwalkers = 4;
69 for (int iw = 0; iw < nwalkers; ++iw)
70 walkers.emplace_back(2);
71
72 std::vector<ParticleSet> psets;
73
74 SpinDensityNew& crowd_sdn = dynamic_cast<SpinDensityNew&>(*(uptr_crowd_sdn));
75
76 std::vector<QMCT::RealType> rng_reals(nwalkers * QMCT::DIM * 2);
77 rft.makeRngReals(rng_reals);
78 auto it_rng_reals = rng_reals.begin();
79 for (int iw = 0; iw < nwalkers; ++iw)
80 {
81 psets.emplace_back();
82 ParticleSet& pset = psets.back();
83 pset.create(2);
84 pset.R[0] = ParticleSet::PosType(*it_rng_reals++, *it_rng_reals++, *it_rng_reals++);
85 pset.R[1] = ParticleSet::PosType(*it_rng_reals++, *it_rng_reals++, *it_rng_reals++);
86 }
87
88 auto ref_walkers = makeRefVector<OperatorEstBase::MCPWalker>(walkers);
89 auto ref_psets = makeRefVector<ParticleSet>(psets);
90
91 crowd_sdn.accumulate(ref_walkers, ref_psets);
92 }
93 }
94
95 TEST_CASE("SpinDensityNew::SpinDensityNew(SPInput, SpeciesSet)", "[estimators]")
96 {
97 Libxml2Document doc;
98 bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[testing::valid_spindensity_input_grid]);
99 REQUIRE(okay);
100 xmlNodePtr node = doc.getRoot();
101 SpinDensityInput sdi;
102 sdi.readXML(node);
103 SpeciesSet species_set;
104 int ispecies = species_set.addSpecies("C");
105 int iattribute = species_set.addAttribute("membersize");
106 species_set(iattribute, ispecies) = 2;
107 SpinDensityInput sdi_copy = sdi;
108 SpinDensityNew(std::move(sdi), species_set);
109 CrystalLattice<OHMMS_PRECISION, OHMMS_DIM> lattice;
110 CHECK_THROWS(SpinDensityNew(std::move(sdi_copy), lattice, species_set));
111 }
112
113 TEST_CASE("SpinDensityNew::SpinDensityNew(SPInput, Lattice, SpeciesSet)", "[estimators]")
114 {
115 Libxml2Document doc;
116 bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[testing::valid_spindensity_input_no_cell]);
117 REQUIRE(okay);
118 xmlNodePtr node = doc.getRoot();
119 SpinDensityInput sdi;
120 sdi.readXML(node);
121 SpeciesSet species_set;
122 int ispecies = species_set.addSpecies("C");
123 int iattribute = species_set.addAttribute("membersize");
124 species_set(iattribute, ispecies) = 2;
125 auto lattice = testing::makeTestLattice();
126 SpinDensityNew(std::move(sdi), lattice, species_set);
127 }
128
129 TEST_CASE("SpinDensityNew::accumulate", "[estimators]")
130 {
131 using MCPWalker = OperatorEstBase::MCPWalker;
132 using QMCT = QMCTraits;
133
134 Libxml2Document doc;
135 bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]);
136 REQUIRE(okay);
137 xmlNodePtr node = doc.getRoot();
138 SpinDensityInput sdi;
139 sdi.readXML(node);
140 SpeciesSet species_set;
141 int ispecies = species_set.addSpecies("u");
142 ispecies = species_set.addSpecies("d");
143 CHECK(ispecies == 1);
144 int iattribute = species_set.addAttribute("membersize");
145 species_set(iattribute, 0) = 1;
146 species_set(iattribute, 1) = 1;
147
148 SpinDensityNew sdn(std::move(sdi), species_set);
149 std::vector<MCPWalker> walkers;
150 int nwalkers = 4;
151 for (int iw = 0; iw < nwalkers; ++iw)
152 walkers.emplace_back(2);
153
154 std::vector<ParticleSet> psets;
155
156 for (int iw = 0; iw < nwalkers; ++iw)
157 {
158 psets.emplace_back();
159 ParticleSet& pset = psets.back();
160 pset.create(2);
161 pset.R[0] = ParticleSet::PosType(0.00000000, 0.00000000, 0.00000000);
162 pset.R[1] = ParticleSet::PosType(1.68658058, 1.68658058, 1.68658058);
163 }
164
165 auto ref_walkers = makeRefVector<MCPWalker>(walkers);
166 auto ref_psets = makeRefVector<ParticleSet>(psets);
167
168 sdn.accumulate(ref_walkers, ref_psets);
169
170 std::vector<QMCT::RealType>& data_ref = sdn.get_data_ref();
171 // There should be a check that the discretization of particle locations expressed in lattice coords
172 // is correct. This just checks it hasn't changed from how it was in SpinDensity which lacked testing.
173 CHECK(data_ref[555] == 4);
174 CHECK(data_ref[1777] == 4);
175 }
176
177 TEST_CASE("SpinDensityNew::collect(DataLocality::crowd)", "[estimators]")
178 {
179 {
180 using MCPWalker = OperatorEstBase::MCPWalker;
181 using QMCT = QMCTraits;
182
183 Libxml2Document doc;
184 bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]);
185 REQUIRE(okay);
186 xmlNodePtr node = doc.getRoot();
187 SpinDensityInput sdi;
188 sdi.readXML(node);
189 SpeciesSet species_set;
190 int ispecies = species_set.addSpecies("u");
191 ispecies = species_set.addSpecies("d");
192 CHECK(ispecies == 1);
193 int iattribute = species_set.addAttribute("membersize");
194 species_set(iattribute, 0) = 1;
195 species_set(iattribute, 1) = 1;
196
197 SpinDensityNew sdn(std::move(sdi), species_set);
198
199 UPtrVector<OperatorEstBase> crowd_sdns;
200 int ncrowds = 2;
201
202 accumulateFromPsets(ncrowds, sdn, crowd_sdns);
203
204 RefVector<OperatorEstBase> crowd_oeb_refs = convertUPtrToRefVector(crowd_sdns);
205 sdn.collect(crowd_oeb_refs);
206
207 std::vector<QMCT::RealType>& data_ref = sdn.get_data_ref();
208 // There should be a check that the discretization of particle locations expressed in lattice coords
209 // is correct. This just checks it hasn't changed from how it was in SpinDensity which lacked testing.
210 CHECK(data_ref[555] == 4 * ncrowds);
211 CHECK(data_ref[1666] == 4 * ncrowds);
212 }
213 }
214
215 TEST_CASE("SpinDensityNew::collect(DataLocality::rank)", "[estimators]")
216 {
217 {
218 using MCPWalker = OperatorEstBase::MCPWalker;
219 using QMCT = QMCTraits;
220
221 Libxml2Document doc;
222 bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]);
223 REQUIRE(okay);
224 xmlNodePtr node = doc.getRoot();
225 SpinDensityInput sdi;
226 sdi.readXML(node);
227 SpeciesSet species_set;
228 int ispecies = species_set.addSpecies("u");
229 ispecies = species_set.addSpecies("d");
230 CHECK(ispecies == 1);
231 int iattribute = species_set.addAttribute("membersize");
232 species_set(iattribute, 0) = 1;
233 species_set(iattribute, 1) = 1;
234
235 SpinDensityNew sdn(std::move(sdi), species_set, DataLocality::rank);
236
237 auto lattice = testing::makeTestLattice();
238
239 UPtrVector<OperatorEstBase> crowd_sdns;
240 int ncrowds = 2;
241
242 accumulateFromPsets(ncrowds, sdn, crowd_sdns);
243
244 RefVector<OperatorEstBase> crowd_oeb_refs = convertUPtrToRefVector(crowd_sdns);
245 sdn.collect(crowd_oeb_refs);
246
247 std::vector<QMCT::RealType>& data_ref = sdn.get_data_ref();
248 // There should be a check that the discretization of particle locations expressed in lattice coords
249 // is correct. This just checks it hasn't changed from how it was in SpinDensity which lacked testing.
250 CHECK(data_ref[555] == 4 * ncrowds);
251 CHECK(data_ref[1666] == 4 * ncrowds);
252 }
253 }
254
255 TEST_CASE("SpinDensityNew algorithm comparison", "[estimators]")
256 {
257 using MCPWalker = OperatorEstBase::MCPWalker;
258 using QMCT = QMCTraits;
259
260 Libxml2Document doc;
261 bool okay = doc.parseFromString(testing::valid_spin_density_input_sections[0]);
262 REQUIRE(okay);
263 xmlNodePtr node = doc.getRoot();
264 SpinDensityInput sdi;
265 sdi.readXML(node);
266 SpeciesSet species_set;
267 int ispecies = species_set.addSpecies("u");
268 ispecies = species_set.addSpecies("d");
269 CHECK(ispecies == 1);
270 int iattribute = species_set.addAttribute("membersize");
271 species_set(iattribute, 0) = 1;
272 species_set(iattribute, 1) = 1;
273
274 SpinDensityInput sdi_copy = sdi;
275
276 int ncrowds = 3;
277 int nsteps = 4;
278
279 SpinDensityNew sdn_rank(std::move(sdi), species_set, DataLocality::rank);
280 UPtrVector<OperatorEstBase> crowd_sdns_rank;
281 accumulateFromPsets(ncrowds, sdn_rank, crowd_sdns_rank);
282 testing::RandomForTest<QMCT::RealType> rng_for_test_rank;
283 for (int i = 0; i < nsteps; ++i)
284 randomUpdateAccumulate(rng_for_test_rank, crowd_sdns_rank);
285 RefVector<OperatorEstBase> crowd_oeb_refs_rank = convertUPtrToRefVector(crowd_sdns_rank);
286 sdn_rank.collect(crowd_oeb_refs_rank);
287 std::vector<QMCT::RealType>& data_ref_rank = sdn_rank.get_data_ref();
288
289 SpinDensityNew sdn_crowd(std::move(sdi), species_set, DataLocality::crowd);
290 UPtrVector<OperatorEstBase> crowd_sdns_crowd;
291 accumulateFromPsets(ncrowds, sdn_crowd, crowd_sdns_crowd);
292 testing::RandomForTest<QMCT::RealType> rng_for_test_crowd;
293 for (int i = 0; i < nsteps; ++i)
294 randomUpdateAccumulate(rng_for_test_crowd, crowd_sdns_crowd);
295 RefVector<OperatorEstBase> crowd_oeb_refs_crowd = convertUPtrToRefVector(crowd_sdns_crowd);
296 sdn_crowd.collect(crowd_oeb_refs_crowd);
297 std::vector<QMCT::RealType>& data_ref_crowd = sdn_crowd.get_data_ref();
298
299 for (size_t i = 0; i < data_ref_rank.size(); ++i)
300 {
301 if (data_ref_crowd[i] != data_ref_rank[i])
302 FAIL_CHECK("crowd local " << data_ref_crowd[i] << " != rank local " << data_ref_rank[i] << " at index " << i);
303 break;
304 }
305 }
306
307
308 } // namespace qmcplusplus
309