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: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
8 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
9 //
10 // File created by: Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
11 //////////////////////////////////////////////////////////////////////////////////////
12
13
14 #include "DensityMatrices1B.h"
15 #include "OhmmsData/AttributeSet.h"
16 #include "QMCWaveFunctions/TrialWaveFunction.h"
17 #include "Numerics/MatrixOperators.h"
18 #include "Utilities/IteratorUtility.h"
19 #include "Utilities/string_utils.h"
20 #include "QMCWaveFunctions/WaveFunctionFactory.h"
21
22
23 namespace qmcplusplus
24 {
25 using MatrixOperators::diag_product;
26 using MatrixOperators::product;
27 using MatrixOperators::product_AtB;
28
29
DensityMatrices1B(ParticleSet & P,TrialWaveFunction & psi,ParticleSet * Pcl,const WaveFunctionFactory & factory)30 DensityMatrices1B::DensityMatrices1B(ParticleSet& P,
31 TrialWaveFunction& psi,
32 ParticleSet* Pcl,
33 const WaveFunctionFactory& factory)
34 : Lattice(P.Lattice), Psi(psi), Pq(P), Pc(Pcl), wf_factory_(factory)
35 {
36 reset();
37 }
38
39
DensityMatrices1B(DensityMatrices1B & master,ParticleSet & P,TrialWaveFunction & psi)40 DensityMatrices1B::DensityMatrices1B(DensityMatrices1B& master, ParticleSet& P, TrialWaveFunction& psi)
41 : OperatorBase(master), Lattice(P.Lattice), Psi(psi), Pq(P), Pc(master.Pc), wf_factory_(master.wf_factory_)
42 {
43 reset();
44 set_state(master);
45 basis_functions.clone_from(master.basis_functions);
46 initialize();
47 for (int i = 0; i < basis_size; ++i)
48 basis_norms[i] = master.basis_norms[i];
49 }
50
51
~DensityMatrices1B()52 DensityMatrices1B::~DensityMatrices1B()
53 {
54 if (initialized)
55 finalize();
56 }
57
58
makeClone(ParticleSet & P,TrialWaveFunction & psi)59 OperatorBase* DensityMatrices1B::makeClone(ParticleSet& P, TrialWaveFunction& psi)
60 {
61 return new DensityMatrices1B(*this, P, psi);
62 }
63
64
reset()65 void DensityMatrices1B::reset()
66 {
67 // uninitialized data
68 w_trace = NULL;
69 T_trace = NULL;
70 Vq_trace = NULL;
71 Vc_trace = NULL;
72 Vqq_trace = NULL;
73 Vqc_trace = NULL;
74 Vcc_trace = NULL;
75 basis_size = -1;
76 nindex = -1;
77 eindex = -1;
78 uniform_random = NULL;
79 // basic HamiltonianBase info
80 UpdateMode.set(COLLECTABLE, 1);
81 // default values
82 energy_mat = false;
83 integrator = uniform_grid;
84 evaluator = loop;
85 sampling = volume_based;
86 scale = 1.0;
87 center = 0.0;
88 points = 10;
89 samples = 10;
90 warmup = 30;
91 timestep = .5;
92 use_drift = false;
93 warmed_up = false;
94 metric = 1.0;
95 write_acceptance_ratio = false;
96 write_rstats = false;
97 normalized = true;
98 volume_normed = true;
99 check_overlap = false;
100 check_derivatives = false;
101 // trace data is required
102 request.request_scalar("weight");
103 request.request_array("Kinetic_complex");
104 request.request_array("Vq");
105 request.request_array("Vc");
106 request.request_array("Vqq");
107 request.request_array("Vqc");
108 request.request_array("Vcc");
109 // has not been initialized
110 initialized = false;
111 }
112
113
put(xmlNodePtr cur)114 bool DensityMatrices1B::put(xmlNodePtr cur)
115 {
116 // read in parameters from the input xml
117 set_state(cur);
118
119 // resize local data and perform warmup sampling, if necessary
120 initialize();
121
122 // write a report to the log
123 report();
124
125 if (check_overlap)
126 test_overlap();
127
128 return true;
129 }
130
131
set_state(xmlNodePtr cur)132 void DensityMatrices1B::set_state(xmlNodePtr cur)
133 {
134 bool center_defined = false;
135 std::string emstr = "no";
136 std::string igstr = "uniform_grid";
137 std::string evstr = "loop";
138 std::string costr = "no";
139 std::string cdstr = "no";
140 std::string arstr = "no";
141 std::string udstr = "no";
142 std::string wrstr = "no";
143 std::string nmstr = "yes";
144 std::string vnstr = "yes";
145 std::vector<std::string> sposets;
146
147 xmlNodePtr element = cur->xmlChildrenNode;
148 while (element != NULL)
149 {
150 std::string ename((const char*)element->name);
151 if (ename == "parameter")
152 {
153 const XMLAttrString name(element, "name");
154 if (name == "basis")
155 putContent(sposets, element);
156 else if (name == "energy_matrix")
157 putContent(emstr, element);
158 else if (name == "integrator")
159 putContent(igstr, element);
160 else if (name == "evaluator")
161 putContent(evstr, element);
162 else if (name == "scale")
163 putContent(scale, element);
164 else if (name == "center")
165 {
166 putContent(center, element);
167 center_defined = true;
168 }
169 else if (name == "points")
170 putContent(points, element);
171 else if (name == "samples")
172 putContent(samples, element);
173 else if (name == "warmup")
174 putContent(warmup, element);
175 else if (name == "timestep")
176 putContent(timestep, element);
177 else if (name == "use_drift")
178 putContent(udstr, element);
179 else if (name == "check_overlap")
180 putContent(costr, element);
181 else if (name == "check_derivatives")
182 putContent(cdstr, element);
183 else if (name == "acceptance_ratio")
184 putContent(arstr, element);
185 else if (name == "rstats")
186 putContent(wrstr, element);
187 else if (name == "normalized")
188 putContent(nmstr, element);
189 else if (name == "volume_normed")
190 putContent(vnstr, element);
191 }
192 element = element->next;
193 }
194
195 if (scale > 1.0 + 1e-10)
196 {
197 APP_ABORT("DensityMatrices1B::put scale must be less than one");
198 }
199 else if (scale < 0.0 - 1e-10)
200 APP_ABORT("DensityMatrices1B::put scale must be greater than zero");
201
202 // get volume and cell information
203 Lattice.reset();
204 if (!center_defined)
205 center = Lattice.Center;
206 volume = Lattice.Volume * std::exp(DIM * std::log(scale));
207 periodic = Lattice.SuperCellEnum != SUPERCELL_OPEN;
208 rcorner = center - scale * Lattice.Center;
209
210 energy_mat = emstr == "yes";
211 if (igstr == "uniform_grid")
212 {
213 integrator = uniform_grid;
214 sampling = volume_based;
215 samples = pow(points, DIM);
216 metric = volume / samples;
217 ind_dims[0] = pow(points, DIM - 1);
218 for (int d = 1; d < DIM; ++d)
219 ind_dims[d] = ind_dims[d - 1] / points;
220 }
221 else if (igstr == "uniform")
222 {
223 integrator = uniform;
224 sampling = volume_based;
225 metric = volume / samples;
226 }
227 else if (igstr == "density")
228 {
229 integrator = density;
230 sampling = metropolis;
231 metric = 1.0 / samples;
232 }
233 else
234 APP_ABORT("DensityMatrices1B::set_state invalid integrator\n valid options are: uniform_grid, uniform, density");
235
236 if (evstr == "loop")
237 evaluator = loop;
238 else if (evstr == "matrix")
239 evaluator = matrix;
240 else
241 APP_ABORT("DensityMatrices1B::set_state invalid evaluator\n valid options are: loop, matrix");
242
243 normalized = nmstr == "yes";
244 volume_normed = vnstr == "yes";
245 use_drift = udstr == "yes";
246 check_overlap = costr == "yes";
247 check_derivatives = cdstr == "yes";
248 write_rstats = wrstr == "yes";
249 write_acceptance_ratio = arstr == "yes";
250
251
252 // get the sposets that form the basis
253 if (sposets.size() == 0)
254 APP_ABORT("DensityMatrices1B::put basis must have at least one sposet");
255
256 for (int i = 0; i < sposets.size(); ++i)
257 {
258 basis_functions.add(wf_factory_.getSPOSet(sposets[i]));
259 }
260 basis_size = basis_functions.size();
261
262 if (basis_size < 1)
263 APP_ABORT("DensityMatrices1B::put basis_size must be greater than one");
264 }
265
266
set_state(DensityMatrices1B & master)267 void DensityMatrices1B::set_state(DensityMatrices1B& master)
268 {
269 basis_size = master.basis_size;
270 energy_mat = master.energy_mat;
271 integrator = master.integrator;
272 evaluator = master.evaluator;
273 sampling = master.sampling;
274 scale = master.scale;
275 points = master.points;
276 samples = master.samples;
277 warmup = master.warmup;
278 timestep = master.timestep;
279 use_drift = master.use_drift;
280 volume = master.volume;
281 periodic = master.periodic;
282 metric = master.metric;
283 rcorner = master.rcorner;
284 normalized = master.normalized;
285 volume_normed = master.volume_normed;
286 for (int d = 0; d < DIM; ++d)
287 ind_dims[d] = master.ind_dims[d];
288 }
289
290
initialize()291 void DensityMatrices1B::initialize()
292 {
293 // get particle information
294 SpeciesSet& species = Pq.getSpeciesSet();
295 nparticles = Pq.getTotalNum();
296 nspecies = species.size();
297 int natt = species.numAttributes();
298 int isize = species.addAttribute("membersize");
299 if (isize == natt)
300 APP_ABORT("DensityMatrices1B::set_state Species set does not have the required attribute 'membersize'");
301 for (int s = 0; s < nspecies; ++s)
302 species_size.push_back(species(isize, s));
303 for (int s = 0; s < nspecies; ++s)
304 species_name.push_back(species.speciesName[s]);
305
306 // allocate space
307 basis_values.resize(basis_size);
308 integrated_values.resize(basis_size);
309 basis_norms.resize(basis_size);
310 RealType bn_standard = 1.0;
311 if (volume_normed)
312 bn_standard = 1.0 / std::sqrt(volume);
313 for (int i = 0; i < basis_size; ++i)
314 basis_norms[i] = bn_standard;
315
316 rsamples.resize(samples);
317 sample_weights.resize(samples);
318 psi_ratios.resize(nparticles);
319
320 if (evaluator == matrix)
321 {
322 Phi_MB.resize(samples, basis_size);
323 for (int s = 0; s < nspecies; ++s)
324 {
325 int specs_size = species_size[s];
326 Phi_NB.push_back(new Matrix_t(specs_size, basis_size));
327 Psi_NM.push_back(new Matrix_t(specs_size, samples));
328 Phi_Psi_NB.push_back(new Matrix_t(specs_size, basis_size));
329 N_BB.push_back(new Matrix_t(basis_size, basis_size));
330 if (energy_mat)
331 {
332 E_N.push_back(new Vector_t(specs_size));
333 E_BB.push_back(new Matrix_t(basis_size, basis_size));
334 }
335 }
336
337 #ifdef DMCHECK
338 Phi_MBtmp.resize(samples, basis_size);
339 for (int s = 0; s < nspecies; ++s)
340 {
341 int specs_size = species_size[s];
342 Phi_NBtmp.push_back(new Matrix_t(specs_size, basis_size));
343 Psi_NMtmp.push_back(new Matrix_t(specs_size, samples));
344 Phi_Psi_NBtmp.push_back(new Matrix_t(specs_size, basis_size));
345 N_BBtmp.push_back(new Matrix_t(basis_size, basis_size));
346 if (energy_mat)
347 {
348 E_Ntmp.push_back(new Vector_t(specs_size));
349 E_BBtmp.push_back(new Matrix_t(basis_size, basis_size));
350 }
351 }
352 #endif
353 }
354
355 if (sampling == metropolis)
356 {
357 basis_gradients.resize(basis_size);
358 basis_laplacians.resize(basis_size);
359 }
360
361 if (!normalized)
362 {
363 normalize();
364 }
365
366 const TimerNameList_t<DMTimers> DMTimerNames = {{DM_eval, "DensityMatrices1B::evaluate"},
367 {DM_gen_samples, "DensityMatrices1B::generate_samples"},
368 {DM_gen_sample_basis, "DensityMatrices1B::generate_sample_basis"},
369 {DM_gen_sample_ratios, "DensityMatrices1B::generate_sample_ratios"},
370 {DM_gen_particle_basis, "DensityMatrices1B::generate_particle_basis"},
371 {DM_matrix_products, "DensityMatrices1B::evaluate_matrix_products"},
372 {DM_accumulate, "DensityMatrices1B::evaluate_matrix_accum"}};
373 setup_timers(timers, DMTimerNames, timer_level_fine);
374
375 initialized = true;
376 }
377
378
finalize()379 void DensityMatrices1B::finalize()
380 {
381 delete_iter(Phi_NB.begin(), Phi_NB.end());
382 delete_iter(Psi_NM.begin(), Psi_NM.end());
383 delete_iter(Phi_Psi_NB.begin(), Phi_Psi_NB.end());
384 delete_iter(N_BB.begin(), N_BB.end());
385 if (energy_mat)
386 {
387 delete_iter(E_N.begin(), E_N.end());
388 delete_iter(E_BB.begin(), E_BB.end());
389 }
390
391 #ifdef DMCHECK
392 delete_iter(Phi_NBtmp.begin(), Phi_NBtmp.end());
393 delete_iter(Psi_NMtmp.begin(), Psi_NMtmp.end());
394 delete_iter(Phi_Psi_NBtmp.begin(), Phi_Psi_NBtmp.end());
395 delete_iter(N_BBtmp.begin(), N_BBtmp.end());
396 if (energy_mat)
397 {
398 delete_iter(E_Ntmp.begin(), E_Ntmp.end());
399 delete_iter(E_BBtmp.begin(), E_BBtmp.end());
400 }
401 #endif
402 }
403
404
report(const std::string & pad)405 void DensityMatrices1B::report(const std::string& pad)
406 {
407 std::vector<std::string> integrator_list;
408 integrator_list.push_back("uniform_grid");
409 integrator_list.push_back("uniform");
410 integrator_list.push_back("density");
411 integrator_list.push_back("no_integrator");
412 std::vector<std::string> evaluator_list;
413 evaluator_list.push_back("loop");
414 evaluator_list.push_back("matrix");
415 evaluator_list.push_back("no_evaluator");
416 std::vector<std::string> sampling_list;
417 sampling_list.push_back("volume_based");
418 sampling_list.push_back("metropolis");
419 sampling_list.push_back("no_sampling");
420
421 std::ostream& out = app_log();
422 //std::ostream& out = std::cout;
423
424 out << pad << "DensityMatrices1B" << std::endl;
425
426 out << pad << " integrator = " << integrator_list[(int)integrator] << std::endl;
427 out << pad << " sampling = " << sampling_list[(int)sampling] << std::endl;
428 out << pad << " evaluator = " << evaluator_list[(int)evaluator] << std::endl;
429 out << pad << " periodic = " << periodic << std::endl;
430 if (sampling == volume_based)
431 {
432 PosType rmax = rcorner + 2 * scale * Lattice.Center;
433 out << pad << " points = " << points << std::endl;
434 out << pad << " scale = " << scale << std::endl;
435 out << pad << " center = " << center << std::endl;
436 out << pad << " rmin = " << rcorner << std::endl;
437 out << pad << " rmax = " << rmax << std::endl;
438 out << pad << " volume = " << volume << std::endl;
439 }
440 else if (sampling == metropolis)
441 {
442 out << pad << " warmup = " << warmup << std::endl;
443 out << pad << " timestep = " << timestep << std::endl;
444 }
445 out << pad << " metric = " << metric << std::endl;
446 out << pad << " nparticles = " << nparticles << std::endl;
447 out << pad << " nspecies = " << nspecies << std::endl;
448 for (int s = 0; s < nspecies; ++s)
449 out << pad << " species " << s << " = " << species_size[s] << std::endl;
450 out << pad << " basis_size = " << basis_size << std::endl;
451 out << pad << " samples = " << samples << std::endl;
452 out << pad << " energy_mat = " << energy_mat << std::endl;
453 out << pad << " initialized = " << initialized << std::endl;
454 out << pad << " normalized = " << normalized << std::endl;
455 out << pad << " volume_normed = " << volume_normed << std::endl;
456 out << pad << " rsamples : " << rsamples.size() << std::endl;
457 if (evaluator == matrix)
458 {
459 out << pad << " Phi_MB : " << Phi_MB.rows() << " " << Phi_MB.cols() << std::endl;
460 for (int s = 0; s < nspecies; ++s)
461 {
462 out << pad << " matrices/vectors for species " << s << std::endl;
463 if (energy_mat)
464 if (energy_mat)
465 out << pad << " E_N : " << E_N[s]->size() << std::endl;
466 out << pad << " Phi_NB : " << Phi_NB[s]->rows() << " " << Phi_NB[s]->cols() << std::endl;
467 out << pad << " Psi_NM : " << Psi_NM[s]->rows() << " " << Psi_NM[s]->cols() << std::endl;
468 out << pad << " Phi_Psi_NB : " << Phi_Psi_NB[s]->rows() << " " << Phi_Psi_NB[s]->cols() << std::endl;
469 out << pad << " N_BB : " << N_BB[s]->rows() << " " << N_BB[s]->cols() << std::endl;
470 if (energy_mat)
471 if (energy_mat)
472 out << pad << " E_BB : " << E_BB[s]->rows() << " " << E_BB[s]->cols() << std::endl;
473 }
474 }
475 out << pad << " basis_norms" << std::endl;
476 for (int i = 0; i < basis_size; ++i)
477 out << pad << " " << i << " " << basis_norms[i] << std::endl;
478 out << pad << "end DensityMatrices1B" << std::endl;
479 }
480
481
get_required_traces(TraceManager & tm)482 void DensityMatrices1B::get_required_traces(TraceManager& tm)
483 {
484 w_trace = tm.get_real_trace("weight");
485 if (energy_mat)
486 {
487 T_trace = tm.get_complex_trace(Pq, "Kinetic_complex");
488 Vq_trace = tm.get_real_combined_trace(Pq, "Vq");
489 Vqq_trace = tm.get_real_combined_trace(Pq, "Vqq");
490 Vqc_trace = tm.get_real_combined_trace(Pq, "Vqc");
491 if (Pc)
492 {
493 Vc_trace = tm.get_real_combined_trace(*Pc, "Vc");
494 Vcc_trace = tm.get_real_combined_trace(*Pc, "Vcc");
495 }
496
497 E_samp.resize(nparticles);
498 }
499 have_required_traces = true;
500 }
501
502
setRandomGenerator(RandomGenerator_t * rng)503 void DensityMatrices1B::setRandomGenerator(RandomGenerator_t* rng) { uniform_random = rng; }
504
505
addObservables(PropertySetType & plist,BufferType & collectables)506 void DensityMatrices1B::addObservables(PropertySetType& plist, BufferType& collectables)
507 {
508 #if defined(QMC_COMPLEX)
509 int nentries = 2 * basis_size * basis_size * nspecies;
510 #else
511 int nentries = basis_size * basis_size * nspecies;
512 #endif
513 myIndex = collectables.current();
514 nindex = myIndex;
515 std::vector<RealType> ntmp(nentries);
516 collectables.add(ntmp.begin(), ntmp.end());
517 if (energy_mat)
518 {
519 eindex = nindex + nentries;
520 std::vector<RealType> etmp(nentries);
521 collectables.add(etmp.begin(), etmp.end());
522 }
523 }
524
525
registerCollectables(std::vector<observable_helper * > & h5desc,hid_t gid) const526 void DensityMatrices1B::registerCollectables(std::vector<observable_helper*>& h5desc, hid_t gid) const
527 {
528 #if defined(QMC_COMPLEX)
529 std::vector<int> ng(3);
530 ng[0] = basis_size;
531 ng[1] = basis_size;
532 ng[2] = 2;
533 int nentries = ng[0] * ng[1] * ng[2];
534 #else
535 std::vector<int> ng(2);
536 ng[0] = basis_size;
537 ng[1] = basis_size;
538 int nentries = ng[0] * ng[1];
539 #endif
540
541 std::string dname = myName;
542 hid_t dgid = H5Gcreate(gid, dname.c_str(), 0);
543
544 std::string nname = "number_matrix";
545 hid_t ngid = H5Gcreate(dgid, nname.c_str(), 0);
546 for (int s = 0; s < nspecies; ++s)
547 {
548 observable_helper* oh;
549 oh = new observable_helper(species_name[s]);
550 oh->set_dimensions(ng, nindex + s * nentries);
551 oh->open(ngid);
552 h5desc.push_back(oh);
553 }
554
555 if (energy_mat)
556 {
557 std::string ename = "energy_matrix";
558 hid_t egid = H5Gcreate(dgid, ename.c_str(), 0);
559 for (int s = 0; s < nspecies; ++s)
560 {
561 observable_helper* oh;
562 oh = new observable_helper(species_name[s]);
563 oh->set_dimensions(ng, eindex + s * nentries);
564 oh->open(egid);
565 h5desc.push_back(oh);
566 }
567 }
568 }
569
570
warmup_sampling()571 void DensityMatrices1B::warmup_sampling()
572 {
573 if (sampling == metropolis)
574 {
575 if (!warmed_up)
576 {
577 diffusion(std::sqrt(timestep), rpcur);
578 rpcur += center;
579 if (integrator == density)
580 density_drift(rpcur, rhocur, dpcur);
581 else
582 APP_ABORT("DensityMatrices1B::warmup_sampling invalid integrator");
583 }
584 generate_samples(1.0, warmup);
585 warmed_up = true;
586 }
587 }
588
589
evaluate(ParticleSet & P)590 DensityMatrices1B::Return_t DensityMatrices1B::evaluate(ParticleSet& P)
591 {
592 ScopedTimer t(timers[DM_eval]);
593 if (have_required_traces || !energy_mat)
594 {
595 if (check_derivatives)
596 test_derivatives();
597 if (evaluator == loop)
598 evaluate_loop(P);
599 else if (evaluator == matrix)
600 evaluate_matrix(P);
601 else
602 APP_ABORT("DensityMatrices1B::evaluate invalid evaluator");
603 }
604 return 0.0;
605 }
606
607
evaluate_matrix(ParticleSet & P)608 DensityMatrices1B::Return_t DensityMatrices1B::evaluate_matrix(ParticleSet& P)
609 {
610 //perform warmup sampling the first time
611 if (!warmed_up)
612 warmup_sampling();
613 // get weight and single particle energy trace data
614 RealType weight;
615 if (energy_mat)
616 weight = w_trace->sample[0] * metric;
617 else
618 weight = tWalker->Weight * metric;
619
620 if (energy_mat)
621 get_energies(E_N); // energies : particles x 1
622 // compute sample positions (monte carlo or deterministic)
623 generate_samples(weight);
624 // compute basis and wavefunction ratio values in matrix form
625 generate_sample_basis(Phi_MB); // basis : samples x basis_size
626 generate_sample_ratios(Psi_NM); // conj(Psi ratio) : particles x samples
627 generate_particle_basis(P, Phi_NB); // conj(basis) : particles x basis_size
628 // perform integration via matrix products
629 {
630 ScopedTimer local_timer(timers[DM_matrix_products]);
631 for (int s = 0; s < nspecies; ++s)
632 {
633 Matrix_t& Psi_nm = *Psi_NM[s];
634 Matrix_t& Phi_Psi_nb = *Phi_Psi_NB[s];
635 Matrix_t& Phi_nb = *Phi_NB[s];
636 diag_product(Psi_nm, sample_weights, Psi_nm);
637 product(Psi_nm, Phi_MB, Phi_Psi_nb); // ratio*basis : particles x basis_size
638 product_AtB(Phi_nb, Phi_Psi_nb, *N_BB[s]); // conj(basis)^T*ratio*basis : basis_size^2
639 if (energy_mat)
640 {
641 Vector_t& E = *E_N[s];
642 diag_product(E, Phi_nb, Phi_nb); // diag(energies)*qmcplusplus::conj(basis)
643 product_AtB(Phi_nb, Phi_Psi_nb, *E_BB[s]); // (energies*conj(basis))^T*ratio*basis
644 }
645 }
646 }
647 // accumulate data into collectables
648 {
649 ScopedTimer local_timer(timers[DM_accumulate]);
650 const int basis_size2 = basis_size * basis_size;
651 int ij = nindex;
652 for (int s = 0; s < nspecies; ++s)
653 {
654 //int ij=nindex; // for testing
655 const Matrix_t& NDM = *N_BB[s];
656 for (int n = 0; n < basis_size2; ++n)
657 {
658 Value_t val = NDM(n);
659 P.Collectables[ij] += real(val);
660 ij++;
661 #if defined(QMC_COMPLEX)
662 P.Collectables[ij] += imag(val);
663 ij++;
664 #endif
665 }
666 }
667 if (energy_mat)
668 {
669 int ij = eindex;
670 for (int s = 0; s < nspecies; ++s)
671 {
672 //int ij=eindex; // for testing
673 const Matrix_t& EDM = *E_BB[s];
674 for (int n = 0; n < basis_size2; ++n)
675 {
676 Value_t val = EDM(n);
677 P.Collectables[ij] += real(val);
678 ij++;
679 #if defined(QMC_COMPLEX)
680 P.Collectables[ij] += imag(val);
681 ij++;
682 #endif
683 }
684 }
685 }
686 }
687
688
689 // jtk come back to this
690 // why are matrices so similar across species?
691 // O atom, einspline, points=20, blocks=1, steps=1
692 //
693 //app_log()<<" ntraces = ";
694 //for(int s=0;s<nspecies;++s)
695 //{
696 // const Matrix_t& NDM = *N_BB[s];
697 // Value_t ntrace = 0.0;
698 // for(int i=0;i<basis_size;++i)
699 // ntrace+=NDM(i,i);
700 // app_log()<<ntrace<<" ";
701 //}
702 //app_log()<< std::endl;
703 //
704 //app_log()<<"nmatrices"<< std::endl;
705 //for(int s=0;s<nspecies;++s)
706 //{
707 // app_log()<<" species "<<s<< std::endl;
708 // const Matrix_t& NDM = *N_BB[s];
709 // for(int i=0;i<basis_size;++i)
710 // {
711 // for(int j=0;j<basis_size;++j)
712 // app_log()<<" "<<real(NDM(i,j));
713 // app_log()<< std::endl;
714 // }
715 //}
716 //
717 //app_log()<<"positions"<< std::endl;
718 //int p=0;
719 //for(int s=0;s<nspecies;++s)
720 //{
721 // app_log()<<" species "<<s<< std::endl;
722 // for(int ps=0;ps<species_size[s];++ps,++p)
723 // app_log()<<" "<<p<<" "<<P.R[p]-P.Lattice.Center<< std::endl;
724 //}
725 //
726 ////app_log()<<"basis_values"<< std::endl;
727 ////int p=0;
728 ////for(int s=0;s<nspecies;++s)
729 ////{
730 //// app_log()<<" species "<<s<< std::endl;
731 //// for(int ps=0;ps<species_size[s];++ps,++p)
732 //// app_log()<<" "<<p<<" "<<P.R[p]<< std::endl;
733 ////}
734
735
736 #ifdef DMCHECK
737 report();
738 app_log() << "DM Check" << std::endl;
739 evaluate_check(P);
740 compare(" Phi_MB", Phi_MB, Phi_MBtmp);
741 for (int s = 0; s < nspecies; ++s)
742 {
743 app_log() << " species " << s << std::endl;
744 compare(" E_N ", *E_N[s], *E_Ntmp[s]);
745 compare(" Phi_NB ", *Phi_NB[s], *Phi_NBtmp[s]);
746 compare(" Psi_NM ", *Psi_NM[s], *Psi_NMtmp[s]);
747 compare(" Phi_Psi_NB", *Phi_Psi_NB[s], *Phi_Psi_NBtmp[s], true);
748 compare(" N_BB ", *N_BB[s], *N_BBtmp[s], true);
749 if (energy_mat)
750 compare(" E_BB ", *E_BB[s], *E_BBtmp[s], true);
751 }
752 app_log() << "end DM Check" << std::endl;
753 APP_ABORT("DM Check");
754 #endif
755
756
757 return 0.0;
758 }
759
760
evaluate_check(ParticleSet & P)761 DensityMatrices1B::Return_t DensityMatrices1B::evaluate_check(ParticleSet& P)
762 {
763 #ifdef DMCHECK
764 APP_ABORT("DensityMatrices1B::evaluate_check use of E_trace in this function needs to be replaces with "
765 "get_energies() and E_samp");
766 int n = 0;
767 for (int s = 0; s < nspecies; ++s)
768 {
769 Matrix_t& Phi_mb = Phi_MBtmp;
770 Matrix_t& Psi_nm = *Psi_NMtmp[s];
771 Matrix_t& Phi_Psi_nb = *Phi_Psi_NBtmp[s];
772 Matrix_t& Phi_nb = *Phi_NBtmp[s];
773 Vector_t& E_n = *E_Ntmp[s];
774 Matrix_t& N_bb = *N_BBtmp[s];
775 Matrix_t& E_bb = *E_BBtmp[s];
776
777 for (int ij = 0; ij < basis_size * basis_size; ++ij)
778 N_bb(ij) = 0.0;
779 for (int ij = 0; ij < basis_size * basis_size; ++ij)
780 E_bb(ij) = 0.0;
781
782 for (int ns = 0; ns < species_size[s]; ++ns, ++n)
783 {
784 std::fill(integrated_values.begin(), integrated_values.end(), 0.0);
785 for (int m = 0; m < samples; ++m)
786 {
787 PosType& rsamp = rsamples[m];
788 update_basis(rsamp);
789 PosType dr = rsamp - P.R[n];
790 P.makeMove(n, dr);
791 Value_t ratio = sample_weights[m] * qmcplusplus::conj(Psi.calcRatio(P, n));
792 P.rejectMove(n);
793 for (int i = 0; i < basis_size; ++i)
794 {
795 integrated_values[i] += ratio * basis_values[i];
796 Phi_mb(m, i) = basis_values[i];
797 }
798 Psi_nm(ns, m) = ratio;
799 }
800 update_basis(P.R[n]);
801 for (int i = 0; i < basis_size; ++i)
802 Phi_Psi_nb(ns, i) = integrated_values[i];
803 for (int i = 0; i < basis_size; ++i)
804 Phi_nb(ns, i) = qmcplusplus::conj(basis_values[i]);
805 for (int i = 0; i < basis_size; ++i)
806 {
807 Value_t phi_i = qmcplusplus::conj(basis_values[i]);
808 for (int j = 0; j < basis_size; ++j)
809 {
810 Value_t val = phi_i * integrated_values[j];
811 N_bb(i, j) += val;
812 }
813 }
814 if (energy_mat)
815 {
816 RealType e_n = E_trace->sample[n]; //replace this with traces access later
817 E_n[ns] = e_n;
818 for (int i = 0; i < basis_size; ++i)
819 {
820 Value_t ephi_i = e_n * qmcplusplus::conj(basis_values[i]);
821 Phi_nb(ns, i) = ephi_i;
822 for (int j = 0; j < basis_size; ++j)
823 {
824 Value_t val = ephi_i * integrated_values[j];
825 E_bb(i, j) += val;
826 }
827 }
828 }
829 }
830 }
831 #endif
832 return 0.0;
833 }
834
835
evaluate_loop(ParticleSet & P)836 DensityMatrices1B::Return_t DensityMatrices1B::evaluate_loop(ParticleSet& P)
837 {
838 const int basis_size2 = basis_size * basis_size;
839 if (!warmed_up)
840 warmup_sampling();
841 RealType weight;
842 if (energy_mat)
843 weight = w_trace->sample[0] * metric;
844 else
845 weight = tWalker->Weight * metric;
846 int nparticles = P.getTotalNum();
847 generate_samples(weight);
848 int n = 0;
849 for (int s = 0; s < nspecies; ++s)
850 {
851 for (int ns = 0; ns < species_size[s]; ++ns, ++n)
852 {
853 integrate(P, n);
854 update_basis(P.R[n]);
855 int ij = nindex + s * basis_size2;
856 for (int i = 0; i < basis_size; ++i)
857 {
858 Value_t phi_i = qmcplusplus::conj(basis_values[i]);
859 for (int j = 0; j < basis_size; ++j)
860 {
861 Value_t val = phi_i * integrated_values[j];
862 P.Collectables[ij] += real(val);
863 ij++;
864 #if defined(QMC_COMPLEX)
865 P.Collectables[ij] += imag(val);
866 ij++;
867 #endif
868 }
869 }
870 if (energy_mat)
871 {
872 RealType e_n = E_trace->sample[n]; //replace this with traces access later
873 int ij = eindex + s * basis_size2;
874 for (int i = 0; i < basis_size; ++i)
875 {
876 Value_t ephi_i = e_n * qmcplusplus::conj(basis_values[i]);
877 for (int j = 0; j < basis_size; ++j)
878 {
879 Value_t val = ephi_i * integrated_values[j];
880 P.Collectables[ij] += real(val);
881 ij++;
882 #if defined(QMC_COMPLEX)
883 P.Collectables[ij] += imag(val);
884 ij++;
885 #endif
886 }
887 }
888 }
889 }
890 }
891 return 0.0;
892 }
893
894
generate_samples(RealType weight,int steps)895 inline void DensityMatrices1B::generate_samples(RealType weight, int steps)
896 {
897 ScopedTimer t(timers[DM_gen_samples]);
898 RandomGenerator_t& rng = *uniform_random;
899 bool save = false;
900 if (steps == 0)
901 {
902 save = true;
903 steps = samples;
904 }
905 if (integrator == uniform_grid)
906 generate_uniform_grid(rng);
907 else if (integrator == uniform)
908 generate_uniform_samples(rng);
909 else if (integrator == density)
910 generate_density_samples(save, steps, rng);
911
912 if (save)
913 {
914 if (sampling == metropolis)
915 {
916 sample_weights *= weight;
917 }
918 else
919 {
920 std::fill(sample_weights.begin(), sample_weights.end(), weight);
921 }
922 }
923
924
925 // temporary check
926 if (write_rstats && omp_get_thread_num() == 0)
927 {
928 PosType rmin = std::numeric_limits<RealType>::max();
929 PosType rmax = -std::numeric_limits<RealType>::max();
930 PosType rmean = 0.0;
931 PosType rstd = 0.0;
932 for (int s = 0; s < rsamples.size(); ++s)
933 for (int d = 0; d < DIM; ++d)
934 {
935 RealType rd = rsamples[s][d];
936 rmin[d] = std::min(rmin[d], rd);
937 rmax[d] = std::max(rmax[d], rd);
938 rmean[d] += rd;
939 rstd[d] += rd * rd;
940 }
941 rmean /= rsamples.size();
942 rstd /= rsamples.size();
943 for (int d = 0; d < DIM; ++d)
944 rstd[d] = std::sqrt(rstd[d] - rmean[d] * rmean[d]);
945 app_log() << "\nrsamples properties:" << std::endl;
946 app_log() << " rmin = " << rmin << std::endl;
947 app_log() << " rmax = " << rmax << std::endl;
948 app_log() << " rmean = " << rmean << std::endl;
949 app_log() << " rstd = " << rstd << std::endl;
950 }
951 }
952
953
generate_uniform_grid(RandomGenerator_t & rng)954 inline void DensityMatrices1B::generate_uniform_grid(RandomGenerator_t& rng)
955 {
956 PosType rp;
957 PosType ushift = 0.0;
958 RealType du = scale / points;
959 for (int d = 0; d < DIM; ++d)
960 ushift[d] += rng() * du;
961 for (int s = 0; s < samples; ++s)
962 {
963 int nrem = s;
964 for (int d = 0; d < DIM - 1; ++d)
965 {
966 int ind = nrem / ind_dims[d];
967 rp[d] = ind * du + ushift[d];
968 nrem -= ind * ind_dims[d];
969 }
970 rp[DIM - 1] = nrem * du + ushift[DIM - 1];
971 rsamples[s] = Lattice.toCart(rp) + rcorner;
972 }
973 }
974
975
generate_uniform_samples(RandomGenerator_t & rng)976 inline void DensityMatrices1B::generate_uniform_samples(RandomGenerator_t& rng)
977 {
978 PosType rp;
979 for (int s = 0; s < samples; ++s)
980 {
981 for (int d = 0; d < DIM; ++d)
982 rp[d] = scale * rng();
983 rsamples[s] = Lattice.toCart(rp) + rcorner;
984 }
985 }
986
987
generate_density_samples(bool save,int steps,RandomGenerator_t & rng)988 inline void DensityMatrices1B::generate_density_samples(bool save, int steps, RandomGenerator_t& rng)
989 {
990 RealType sqt = std::sqrt(timestep);
991 RealType ot = 1.0 / timestep;
992 PosType r = rpcur; //current position
993 PosType d = dpcur; //current drift
994 RealType rho = rhocur; //current density
995 for (int s = 0; s < steps; ++s)
996 {
997 nmoves++;
998 PosType n, rp, dp, ds; //diffusion, trial pos/drift, drift sum
999 RealType rhop, ratio, Pacc; //trial density, dens ratio, acc prob
1000 diffusion(sqt, n); //get diffusion
1001 if (use_drift)
1002 {
1003 rp = r + n + d; //update trial position
1004 density_drift(rp, rhop, dp); //get trial drift and density
1005 ratio = rhop / rho; //density ratio
1006 ds = dp + d; //drift sum
1007 Pacc = ratio * std::exp(-ot * (dot(n, ds) + .5 * dot(ds, ds))); //acceptance probability
1008 }
1009 else
1010 {
1011 rp = r + n; //update trial position
1012 density_only(rp, rhop); //get trial density
1013 ratio = rhop / rho; //density ratio
1014 Pacc = ratio; //acceptance probability
1015 }
1016 if (rng() < Pacc)
1017 { //accept move
1018 r = rp;
1019 d = dp;
1020 rho = rhop;
1021 naccepted++;
1022 }
1023 if (save)
1024 {
1025 rsamples[s] = r;
1026 sample_weights[s] = 1.0 / rho;
1027 }
1028 }
1029 acceptance_ratio = RealType(naccepted) / nmoves;
1030
1031 if (write_acceptance_ratio && omp_get_thread_num() == 0)
1032 app_log() << "dm1b acceptance_ratio = " << acceptance_ratio << std::endl;
1033
1034 rpcur = r;
1035 dpcur = d;
1036 rhocur = rho;
1037 }
1038
1039
diffusion(RealType sqt,PosType & diff)1040 inline void DensityMatrices1B::diffusion(RealType sqt, PosType& diff)
1041 {
1042 assignGaussRand(&diff[0], DIM, *uniform_random);
1043 diff *= sqt;
1044 }
1045
1046
density_only(const PosType & r,RealType & dens)1047 inline void DensityMatrices1B::density_only(const PosType& r, RealType& dens)
1048 {
1049 update_basis(r);
1050 dens = 0.0;
1051 for (int i = 0; i < basis_size; ++i)
1052 {
1053 Value_t b = basis_values[i];
1054 dens += std::abs(qmcplusplus::conj(b) * b);
1055 }
1056 dens /= basis_size;
1057 }
1058
1059
density_drift(const PosType & r,RealType & dens,PosType & drift)1060 inline void DensityMatrices1B::density_drift(const PosType& r, RealType& dens, PosType& drift)
1061 {
1062 update_basis_d012(r);
1063 dens = 0.0;
1064 drift = 0.0;
1065 for (int i = 0; i < basis_size; ++i)
1066 {
1067 const Grad_t& bg = basis_gradients[i];
1068 Value_t b = basis_values[i];
1069 Value_t bc = qmcplusplus::conj(b);
1070 dens += std::abs(bc * b);
1071 for (int d = 0; d < DIM; ++d)
1072 drift[d] += std::real(bc * bg[d]);
1073 }
1074 drift *= timestep / dens;
1075 dens /= basis_size;
1076 }
1077
1078
1079 typedef DensityMatrices1B::RealType RealType;
1080 typedef DensityMatrices1B::Value_t Value_t;
1081
1082
accum_constant(CombinedTraceSample<TraceReal> * etrace,RealType weight=1.0)1083 inline RealType accum_constant(CombinedTraceSample<TraceReal>* etrace, RealType weight = 1.0)
1084 {
1085 RealType E = 0.0;
1086 if (etrace)
1087 {
1088 etrace->combine();
1089 for (int p = 0; p < etrace->sample.size(); ++p)
1090 E += etrace->sample[p];
1091 }
1092 return weight * E;
1093 }
1094
1095 template<typename T>
accum_sample(std::vector<Value_t> & E_samp,CombinedTraceSample<T> * etrace,RealType weight=1.0)1096 inline void accum_sample(std::vector<Value_t>& E_samp, CombinedTraceSample<T>* etrace, RealType weight = 1.0)
1097 {
1098 if (etrace)
1099 {
1100 etrace->combine();
1101 for (int p = 0; p < etrace->sample.size(); ++p)
1102 E_samp[p] += weight * etrace->sample[p];
1103 }
1104 }
1105
1106 template<typename T>
accum_sample(std::vector<Value_t> & E_samp,TraceSample<T> * etrace,RealType weight=1.0)1107 inline void accum_sample(std::vector<Value_t>& E_samp, TraceSample<T>* etrace, RealType weight = 1.0)
1108 {
1109 #ifdef QMC_COMPLEX
1110 if (etrace)
1111 for (int p = 0; p < etrace->sample.size(); ++p)
1112 E_samp[p] += weight * etrace->sample[p];
1113 #else
1114 if (etrace)
1115 for (int p = 0; p < etrace->sample.size(); ++p)
1116 E_samp[p] += weight * real(etrace->sample[p]);
1117 #endif
1118 }
1119
1120
get_energies(std::vector<Vector_t * > & E_n)1121 void DensityMatrices1B::get_energies(std::vector<Vector_t*>& E_n)
1122 {
1123 Value_t Vc = 0;
1124 Vc += accum_constant(Vc_trace);
1125 Vc += accum_constant(Vcc_trace);
1126 Vc /= nparticles;
1127 std::fill(E_samp.begin(), E_samp.end(), Vc);
1128 accum_sample(E_samp, T_trace);
1129 accum_sample(E_samp, Vq_trace);
1130 accum_sample(E_samp, Vqq_trace);
1131 accum_sample(E_samp, Vqc_trace, 2.0);
1132
1133 int p = 0;
1134 for (int s = 0; s < nspecies; ++s)
1135 {
1136 Vector_t& E = *E_n[s];
1137 for (int ps = 0; ps < E.size(); ++ps, ++p)
1138 E[ps] = E_samp[p];
1139 }
1140
1141 //E_trace->combine();
1142 //RealType E = 0.0;
1143 //for(int p=0;p<nparticles;++p)
1144 // E += E_samp[p];
1145 //app_log()<<" E = "<<E<<" "<<E_trace->sample[0]<< std::endl;
1146 //APP_ABORT("dm1b::get_energies check sp traces");
1147 }
1148
1149
generate_sample_basis(Matrix_t & Phi_mb)1150 void DensityMatrices1B::generate_sample_basis(Matrix_t& Phi_mb)
1151 {
1152 ScopedTimer t(timers[DM_gen_sample_basis]);
1153 int mb = 0;
1154 for (int m = 0; m < samples; ++m)
1155 {
1156 update_basis(rsamples[m]);
1157 for (int b = 0; b < basis_size; ++b, ++mb)
1158 Phi_mb(mb) = basis_values[b];
1159 }
1160 }
1161
1162
generate_sample_ratios(std::vector<Matrix_t * > Psi_nm)1163 void DensityMatrices1B::generate_sample_ratios(std::vector<Matrix_t*> Psi_nm)
1164 {
1165 ScopedTimer t(timers[DM_gen_sample_ratios]);
1166 for (int m = 0; m < samples; ++m)
1167 {
1168 // get N ratios for the current sample point
1169 Pq.makeVirtualMoves(rsamples[m]);
1170 Psi.evaluateRatiosAlltoOne(Pq, psi_ratios);
1171
1172 // collect ratios into per-species matrices
1173 int p = 0;
1174 for (int s = 0; s < nspecies; ++s)
1175 {
1176 Matrix_t& P_nm = *Psi_nm[s];
1177 for (int n = 0; n < species_size[s]; ++n, ++p)
1178 {
1179 P_nm(n, m) = qmcplusplus::conj(psi_ratios[p]);
1180 }
1181 }
1182 }
1183 }
1184
1185
generate_particle_basis(ParticleSet & P,std::vector<Matrix_t * > & Phi_nb)1186 void DensityMatrices1B::generate_particle_basis(ParticleSet& P, std::vector<Matrix_t*>& Phi_nb)
1187 {
1188 ScopedTimer t(timers[DM_gen_particle_basis]);
1189 int p = 0;
1190 for (int s = 0; s < nspecies; ++s)
1191 {
1192 int nb = 0;
1193 Matrix_t& P_nb = *Phi_nb[s];
1194 for (int n = 0; n < species_size[s]; ++n, ++p)
1195 {
1196 update_basis(P.R[p]);
1197 for (int b = 0; b < basis_size; ++b, ++nb)
1198 P_nb(nb) = qmcplusplus::conj(basis_values[b]);
1199 }
1200 }
1201 }
1202
1203
integrate(ParticleSet & P,int n)1204 inline void DensityMatrices1B::integrate(ParticleSet& P, int n)
1205 {
1206 std::fill(integrated_values.begin(), integrated_values.end(), 0.0);
1207 for (int s = 0; s < samples; ++s)
1208 {
1209 PosType& rsamp = rsamples[s];
1210 update_basis(rsamp);
1211 P.makeMove(n, rsamp - P.R[n]);
1212 Value_t ratio = sample_weights[s] * qmcplusplus::conj(Psi.calcRatio(P, n));
1213 P.rejectMove(n);
1214 for (int i = 0; i < basis_size; ++i)
1215 integrated_values[i] += ratio * basis_values[i];
1216 }
1217 }
1218
1219
update_basis(const PosType & r)1220 inline void DensityMatrices1B::update_basis(const PosType& r)
1221 {
1222 Pq.makeMove(0, r - Pq.R[0]);
1223 basis_functions.evaluateValue(Pq, 0, basis_values);
1224 Pq.rejectMove(0);
1225 for (int i = 0; i < basis_size; ++i)
1226 basis_values[i] *= basis_norms[i];
1227 }
1228
1229
update_basis_d012(const PosType & r)1230 inline void DensityMatrices1B::update_basis_d012(const PosType& r)
1231 {
1232 Pq.makeMove(0, r - Pq.R[0]);
1233 basis_functions.evaluateVGL(Pq, 0, basis_values, basis_gradients, basis_laplacians);
1234 Pq.rejectMove(0);
1235 for (int i = 0; i < basis_size; ++i)
1236 basis_values[i] *= basis_norms[i];
1237 for (int i = 0; i < basis_size; ++i)
1238 basis_gradients[i] *= basis_norms[i];
1239 for (int i = 0; i < basis_size; ++i)
1240 basis_laplacians[i] *= basis_norms[i];
1241 }
1242
1243
normalize()1244 inline void DensityMatrices1B::normalize()
1245 {
1246 int ngrid = std::max(200, points);
1247 int ngtot = pow(ngrid, DIM);
1248 RealType du = scale / ngrid;
1249 RealType dV = volume / ngtot;
1250 PosType rp;
1251 ValueVector_t bnorms;
1252 int gdims[DIM];
1253 gdims[0] = pow(ngrid, DIM - 1);
1254 for (int d = 1; d < DIM; ++d)
1255 gdims[d] = gdims[d - 1] / ngrid;
1256 bnorms.resize(basis_size);
1257 for (int i = 0; i < basis_size; ++i)
1258 bnorms[i] = 0.0;
1259 std::fill(basis_norms.begin(), basis_norms.end(), 1.0);
1260 for (int p = 0; p < ngtot; ++p)
1261 {
1262 int nrem = p;
1263 for (int d = 0; d < DIM - 1; ++d)
1264 {
1265 int ind = nrem / gdims[d];
1266 rp[d] = ind * du + du / 2;
1267 nrem -= ind * gdims[d];
1268 }
1269 rp[DIM - 1] = nrem * du + du / 2;
1270 rp = Lattice.toCart(rp) + rcorner;
1271 update_basis(rp);
1272 for (int i = 0; i < basis_size; ++i)
1273 bnorms[i] += qmcplusplus::conj(basis_values[i]) * basis_values[i] * dV;
1274 }
1275 for (int i = 0; i < basis_size; ++i)
1276 basis_norms[i] = 1.0 / std::sqrt(real(bnorms[i]));
1277 normalized = true;
1278 }
1279
1280
test_overlap()1281 inline void DensityMatrices1B::test_overlap()
1282 {
1283 //int ngrid = std::max(200,points);
1284 int ngrid = std::max(50, points);
1285 int ngtot = pow(ngrid, DIM);
1286
1287 PosType rp;
1288 RealType du = scale / ngrid;
1289 RealType dV = volume / ngtot;
1290
1291 PosType rmin = std::numeric_limits<RealType>::max();
1292 PosType rmax = -std::numeric_limits<RealType>::max();
1293 int gdims[DIM];
1294 gdims[0] = pow(ngrid, DIM - 1);
1295 for (int d = 1; d < DIM; ++d)
1296 gdims[d] = gdims[d - 1] / ngrid;
1297
1298 Array<Value_t, 2> omat;
1299 omat.resize(basis_size, basis_size);
1300 for (int i = 0; i < basis_size; ++i)
1301 for (int j = 0; j < basis_size; ++j)
1302 omat(i, j) = 0.0;
1303
1304 for (int p = 0; p < ngtot; ++p)
1305 {
1306 int nrem = p;
1307 for (int d = 0; d < DIM - 1; ++d)
1308 {
1309 int ind = nrem / gdims[d];
1310 rp[d] = ind * du + du / 2;
1311 nrem -= ind * gdims[d];
1312 }
1313 rp[DIM - 1] = nrem * du + du / 2;
1314 rp = Lattice.toCart(rp) + rcorner;
1315 update_basis(rp);
1316 for (int i = 0; i < basis_size; ++i)
1317 for (int j = 0; j < basis_size; ++j)
1318 omat(i, j) += qmcplusplus::conj(basis_values[i]) * basis_values[j] * dV;
1319 for (int d = 0; d < DIM; ++d)
1320 {
1321 rmin[d] = std::min(rmin[d], rp[d]);
1322 rmax[d] = std::max(rmax[d], rp[d]);
1323 }
1324 }
1325
1326 app_log() << "DensityMatrices1B::test_overlap checking overlap matrix" << std::endl;
1327 app_log() << " rmin = " << rmin << std::endl;
1328 app_log() << " rmax = " << rmax << std::endl;
1329 app_log() << " overlap scale " << std::abs(omat(0, 0)) << std::endl;
1330 app_log() << " overlap matrix:" << std::endl;
1331 for (int i = 0; i < basis_size; ++i)
1332 {
1333 app_log() << std::endl;
1334 for (int j = 0; j < basis_size; ++j)
1335 app_log() << std::abs(omat(i, j)) / std::abs(omat(0, 0)) << " ";
1336 }
1337 app_log() << std::endl;
1338 APP_ABORT("DensityMatrices1B::test_overlap");
1339 }
1340
1341
test_derivatives()1342 void DensityMatrices1B::test_derivatives()
1343 {
1344 app_log() << "DensityMatrices1B::test_derivatives checking drift" << std::endl;
1345
1346 PosType r, rtmp;
1347
1348 RealType delta = 1e-5;
1349
1350 RealType dens, densp, densm;
1351 PosType drift, driftn, drifttmp;
1352
1353 app_log() << " warming up" << std::endl;
1354 warmup_sampling();
1355 app_log() << " generating samples" << std::endl;
1356 generate_samples(1.0);
1357
1358 app_log() << " testing derivatives at sample points" << std::endl;
1359 for (int s = 0; s < rsamples.size(); ++s)
1360 {
1361 r = rsamples[s];
1362
1363 density_drift(r, dens, drift);
1364
1365 for (int d = 0; d < DIM; ++d)
1366 {
1367 rtmp = r;
1368
1369 rtmp[d] = r[d] + delta;
1370 density_drift(rtmp, densp, drifttmp);
1371
1372 rtmp[d] = r[d] - delta;
1373 density_drift(rtmp, densm, drifttmp);
1374
1375 driftn[d] = (densp - densm) / (2 * delta);
1376 }
1377 driftn *= .5 * timestep / dens;
1378
1379 app_log() << s << std::endl;
1380 app_log() << " " << driftn << std::endl;
1381 app_log() << " " << drift << std::endl;
1382 }
1383
1384 APP_ABORT("DensityMatrices1B::test_derivatives");
1385 }
1386
1387
match(Value_t e1,Value_t e2,RealType tol)1388 bool DensityMatrices1B::match(Value_t e1, Value_t e2, RealType tol)
1389 {
1390 return std::abs(e1 - e2) < tol;
1391 //return std::abs(2*(e1-e2)/(e1+e2)) < tol;
1392 }
1393
1394
same(Vector_t & v1,Vector_t & v2,RealType tol)1395 bool DensityMatrices1B::same(Vector_t& v1, Vector_t& v2, RealType tol)
1396 {
1397 if (v1.size() != v2.size())
1398 APP_ABORT("DensityMatrices1B::same(vector) vectors differ in size");
1399 bool sm = true;
1400 for (int i = 0; i < v1.size(); ++i)
1401 sm &= match(v1[i], v2[i]);
1402 return sm;
1403 }
1404
same(Matrix_t & m1,Matrix_t & m2,RealType tol)1405 bool DensityMatrices1B::same(Matrix_t& m1, Matrix_t& m2, RealType tol)
1406 {
1407 if (m1.rows() != m2.rows() || m1.cols() != m2.cols())
1408 APP_ABORT("DensityMatrices1B::same(matrix) matrices differ in size");
1409 bool sm = true;
1410 int n = m1.rows() * m1.cols();
1411 for (int i = 0; i < n; ++i)
1412 sm &= match(m1(i), m2(i));
1413 return sm;
1414 }
1415
compare(const std::string & name,Vector_t & v1,Vector_t & v2,bool write,bool diff_only)1416 void DensityMatrices1B::compare(const std::string& name, Vector_t& v1, Vector_t& v2, bool write, bool diff_only)
1417 {
1418 bool sm = same(v1, v2);
1419 std::string result = "differ";
1420 if (sm)
1421 result = "agree";
1422 app_log() << name << " " << result << std::endl;
1423 if (write && !sm)
1424 for (int i = 0; i < v1.size(); ++i)
1425 app_log() << " " << i << " " << real(v1[i]) << " " << real(v2[i]) << " " << real(v1[i] / v2[i]) << " "
1426 << real(v2[i] / v1[i]) << std::endl;
1427 }
1428
compare(const std::string & name,Matrix_t & m1,Matrix_t & m2,bool write,bool diff_only)1429 void DensityMatrices1B::compare(const std::string& name, Matrix_t& m1, Matrix_t& m2, bool write, bool diff_only)
1430 {
1431 bool sm = same(m1, m2);
1432 std::string result = "differ";
1433 if (sm)
1434 result = "agree";
1435 app_log() << name << " " << result << std::endl;
1436 if (write && !sm)
1437 for (int i = 0; i < m1.rows(); ++i)
1438 for (int j = 0; j < m1.cols(); ++j)
1439 if (!diff_only || !match(m1(i, j), m2(i, j)))
1440 app_log() << " " << i << " " << j << " " << real(m1(i, j)) << " " << real(m2(i, j)) << " "
1441 << real(m1(i, j) / m2(i, j)) << std::endl;
1442 }
1443
1444
1445 } // namespace qmcplusplus
1446