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) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
8 //
9 // File created by: Fionn Malone, malone14@llnl.gov, Lawrence Livermore National Laboratory
10 //////////////////////////////////////////////////////////////////////////////////////
11
12 #include "catch.hpp"
13
14 #include "Configuration.h"
15
16 #include "OhmmsData/Libxml2Doc.h"
17 #include "OhmmsApp/ProjectData.h"
18 #include "Utilities/TimerManager.h"
19 #include "hdf/hdf_archive.h"
20
21 #undef APP_ABORT
22 #define APP_ABORT(x) \
23 { \
24 std::cout << x << std::endl; \
25 throw; \
26 }
27
28 #include <stdio.h>
29 #include <string>
30
31 #include "AFQMC/config.h"
32 #include "AFQMC/Matrix/tests/matrix_helpers.h"
33 #include "AFQMC/Hamiltonians/HamiltonianFactory.h"
34 #include "AFQMC/Hamiltonians/Hamiltonian.hpp"
35 #include "AFQMC/Wavefunctions/WavefunctionFactory.h"
36 #include "AFQMC/Wavefunctions/Wavefunction.hpp"
37 #include "AFQMC/Walkers/WalkerSet.hpp"
38 #include "AFQMC/Estimators/EstimatorBase.h"
39 #include "AFQMC/Propagators/PropagatorFactory.h"
40 #include "AFQMC/Propagators/Propagator.hpp"
41 #include "AFQMC/Estimators/BackPropagatedEstimator.hpp"
42 #include "AFQMC/Utilities/test_utils.hpp"
43 #include "AFQMC/Memory/buffer_managers.h"
44
45 using std::cerr;
46 using std::complex;
47 using std::cout;
48 using std::endl;
49 using std::ifstream;
50 using std::setprecision;
51 using std::string;
52
53 extern std::string UTEST_HAMIL, UTEST_WFN;
54
55 namespace qmcplusplus
56 {
57 using namespace afqmc;
58
59 template<class Allocator>
reduced_density_matrix(boost::mpi3::communicator & world)60 void reduced_density_matrix(boost::mpi3::communicator& world)
61 {
62 using pointer = typename Allocator::pointer;
63
64 if (check_hamil_wfn_for_utest("reduced_density_matrix", UTEST_WFN, UTEST_HAMIL))
65 {
66 timer_manager.set_timer_threshold(timer_level_coarse);
67 setup_timers(AFQMCTimers, AFQMCTimerNames, timer_level_coarse);
68
69 // Global Task Group
70 afqmc::GlobalTaskGroup gTG(world);
71
72 int NMO, NAEA, NAEB;
73 std::tie(NMO, NAEA, NAEB) = read_info_from_hdf(UTEST_HAMIL);
74
75 std::map<std::string, AFQMCInfo> InfoMap;
76 InfoMap.insert(std::pair<std::string, AFQMCInfo>("info0", AFQMCInfo{"info0", NMO, NAEA, NAEB}));
77 HamiltonianFactory HamFac(InfoMap);
78 std::string hamil_xml = "<Hamiltonian name=\"ham0\" info=\"info0\"> \
79 <parameter name=\"filetype\">hdf5</parameter> \
80 <parameter name=\"filename\">" +
81 UTEST_HAMIL + "</parameter> \
82 <parameter name=\"cutoff_decomposition\">1e-5</parameter> \
83 </Hamiltonian> \
84 ";
85 const char* ham_xml_block = hamil_xml.c_str();
86 Libxml2Document doc;
87 bool okay = doc.parseFromString(ham_xml_block);
88 REQUIRE(okay);
89 std::string ham_name("ham0");
90 HamFac.push(ham_name, doc.getRoot());
91 Hamiltonian& ham = HamFac.getHamiltonian(gTG, ham_name);
92
93 WALKER_TYPES type = afqmc::getWalkerType(UTEST_WFN);
94 const char* wlk_xml_block_closed = "<WalkerSet name=\"wset0\"> \
95 <parameter name=\"walker_type\">closed</parameter> \
96 <parameter name=\"back_propagation_steps\">10</parameter> \
97 </WalkerSet> \
98 ";
99 const char* wlk_xml_block_coll = "<WalkerSet name=\"wset0\"> \
100 <parameter name=\"walker_type\">collinear</parameter> \
101 <parameter name=\"back_propagation_steps\">10</parameter> \
102 </WalkerSet> \
103 ";
104 const char* wlk_xml_block_noncol = "<WalkerSet name=\"wset0\"> \
105 <parameter name=\"walker_type\">noncollinear</parameter> \
106 <parameter name=\"back_propagation_steps\">10</parameter> \
107 </WalkerSet> \
108 ";
109
110 const char* wlk_xml_block =
111 ((type == CLOSED) ? (wlk_xml_block_closed) : (type == COLLINEAR ? wlk_xml_block_coll : wlk_xml_block_noncol));
112 Libxml2Document doc3;
113 okay = doc3.parseFromString(wlk_xml_block);
114 REQUIRE(okay);
115
116 std::string wfn_xml = "<Wavefunction name=\"wfn0\" info=\"info0\"> \
117 <parameter name=\"filetype\">ascii</parameter> \
118 <parameter name=\"filename\">" +
119 UTEST_WFN + "</parameter> \
120 <parameter name=\"cutoff\">1e-6</parameter> \
121 </Wavefunction> \
122 ";
123 const char* wfn_xml_block = wfn_xml.c_str();
124 auto TG = TaskGroup_(gTG, std::string("WfnTG"), 1, gTG.getTotalCores());
125 Allocator alloc_(make_localTG_allocator<ComplexType>(TG));
126 int nwalk = 1; // choose prime number to force non-trivial splits in shared routines
127 RandomGenerator_t rng;
128 Libxml2Document doc2;
129 okay = doc2.parseFromString(wfn_xml_block);
130 REQUIRE(okay);
131 std::string wfn_name("wfn0");
132 WavefunctionFactory WfnFac(InfoMap);
133 WfnFac.push(wfn_name, doc2.getRoot());
134 Wavefunction& wfn = WfnFac.getWavefunction(TG, TG, wfn_name, type, &ham, 1e-6, nwalk);
135
136 const char* propg_xml_block = "<Propagator name=\"prop0\"> \
137 </Propagator> \
138 ";
139 Libxml2Document doc5;
140 okay = doc5.parseFromString(propg_xml_block);
141 REQUIRE(okay);
142 std::string prop_name("prop0");
143 PropagatorFactory PropgFac(InfoMap);
144 PropgFac.push(prop_name, doc5.getRoot());
145 Propagator& prop = PropgFac.getPropagator(TG, prop_name, wfn, &rng);
146
147 WalkerSet wset(TG, doc3.getRoot(), InfoMap["info0"], &rng);
148 auto initial_guess = WfnFac.getInitialGuess(wfn_name);
149 REQUIRE(initial_guess.size(0) == 2);
150 REQUIRE(initial_guess.size(1) == NMO);
151 REQUIRE(initial_guess.size(2) == NAEA);
152 wset.resize(nwalk, initial_guess[0], initial_guess[0]);
153 using EstimPtr = std::shared_ptr<EstimatorBase>;
154 std::vector<EstimPtr> estimators;
155 const char* est_xml_block = "<Estimator name=\"back_propagation\"> \
156 <parameter name=\"nsteps\">1</parameter> \
157 <parameter name=\"block_size\">2</parameter> \
158 <OneRDM> </OneRDM> \
159 <parameter name=\"path_restoration\">false</parameter> \
160 </Estimator> \
161 ";
162 Libxml2Document doc4;
163 okay = doc4.parseFromString(est_xml_block);
164 REQUIRE(okay);
165 bool impsamp = true;
166 estimators.emplace_back(
167 static_cast<EstimPtr>(std::make_shared<BackPropagatedEstimator>(TG, InfoMap["info0"], "none", doc4.getRoot(),
168 type, wset, wfn, prop, impsamp)));
169
170 // generate P1 with dt=0
171 prop.generateP1(0.0, wset.getWalkerType());
172
173 std::string file = create_test_hdf(UTEST_WFN, UTEST_HAMIL);
174 hdf_archive dump;
175 std::ofstream out;
176 dump.create(file);
177 dump.open(file);
178 for (int iblock = 0; iblock < 10; iblock++)
179 {
180 wset.advanceBPPos();
181 estimators[0]->accumulate_block(wset);
182 estimators[0]->print(out, dump, wset);
183 }
184 dump.close();
185 boost::multi::array<ComplexType, 1> read_data(boost::multi::iextensions<1u>{2 * NMO * NMO});
186
187 ComplexType denom;
188 hdf_archive reader;
189 // Read from a particular block.
190 if (!reader.open(file, H5F_ACC_RDONLY))
191 {
192 app_error() << " Error opening estimates.h5. \n";
193 APP_ABORT("");
194 }
195 reader.read(read_data, "Observables/BackPropagated/FullOneRDM/Average_0/one_rdm_000000004");
196 reader.read(denom, "Observables/BackPropagated/FullOneRDM/Average_0/denominator_000000004");
197 // Test EstimatorHandler eventually.
198 //int NAEA_READ, NAEB_READ, NMO_READ, WALKER_TYPE_READ;
199 //reader.read(NAEA_READ, "Metadata/NAEA");
200 //REQUIRE(NAEA_READ==NAEA);
201 //reader.read(NAEB_READ, "Metadata/NAEB");
202 //REQUIRE(NAEB_READ==NAEB);
203 //reader.read(NMO_READ, "Metadata/NMO");
204 //REQUIRE(NMO_READ==NMO);
205 //reader.read(WALKER_TYPE_READ, "Metadata/WALKER_TYPE");
206 //REQUIRE(WALKER_TYPE_READ==type);
207 reader.close();
208 // Test the RDM. Since no back propagation has been performed the RDM should be
209 // identical to the mixed estimate.
210 if (type == CLOSED)
211 {
212 REQUIRE(read_data.num_elements() >= NMO * NMO);
213 boost::multi::array_ref<ComplexType, 2> BPRDM(read_data.origin(), {NMO, NMO});
214 ma::scal(1.0 / denom, BPRDM);
215 ComplexType trace = ComplexType(0.0);
216 for (int i = 0; i < NMO; i++)
217 trace += BPRDM[i][i];
218 REQUIRE(trace.real() == Approx(NAEA));
219 boost::multi::array<ComplexType, 2, Allocator> Gw({1, NMO * NMO}, alloc_);
220 wfn.MixedDensityMatrix(wset, Gw, false, true);
221 boost::multi::array_ref<ComplexType, 2, pointer> G(Gw.origin(), {NMO, NMO});
222 verify_approx(G, BPRDM);
223 }
224 else if (type == COLLINEAR)
225 {
226 REQUIRE(read_data.num_elements() >= 2 * NMO * NMO);
227 boost::multi::array_ref<ComplexType, 3> BPRDM(read_data.origin(), {2, NMO, NMO});
228 ma::scal(1.0 / denom, BPRDM[0]);
229 ma::scal(1.0 / denom, BPRDM[1]);
230 ComplexType trace = ComplexType(0.0);
231 for (int i = 0; i < NMO; i++)
232 trace += BPRDM[0][i][i] + BPRDM[1][i][i];
233 REQUIRE(trace.real() == Approx(NAEA + NAEB));
234 boost::multi::array<ComplexType, 2, Allocator> Gw({1, 2 * NMO * NMO}, alloc_);
235 wfn.MixedDensityMatrix(wset, Gw, false, true);
236 boost::multi::array_ref<ComplexType, 3, pointer> G(Gw.origin(), {2, NMO, NMO});
237 verify_approx(G, BPRDM);
238 }
239 else
240 {
241 APP_ABORT(" NONCOLLINEAR Wavefunction found.\n");
242 }
243 }
244 }
245
246 TEST_CASE("reduced_density_matrix", "[estimators]")
247 {
248 auto world = boost::mpi3::environment::get_world_instance();
249 if (not world.root())
250 infoLog.pause();
251 auto node = world.split_shared(world.rank());
252
253 #if defined(ENABLE_CUDA) || defined(ENABLE_HIP)
254 arch::INIT(node);
255 using Alloc = device::device_allocator<ComplexType>;
256 #else
257 using Alloc = shared_allocator<ComplexType>;
258 #endif
259 setup_memory_managers(node, 10uL * 1024uL * 1024uL);
260 reduced_density_matrix<Alloc>(world);
261 release_memory_managers();
262 }
263
264 } // namespace qmcplusplus
265