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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 // Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory
9 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
11 //
12 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
13 //////////////////////////////////////////////////////////////////////////////////////
14
15
16 /**@file ParticleSetPool.cpp
17 * @brief Implements ParticleSetPool operators.
18 */
19 #include "ParticleSetPool.h"
20 #include "ParticleBase/RandomSeqGenerator.h"
21 #include "ParticleIO/XMLParticleIO.h"
22 #include "ParticleIO/ParticleLayoutIO.h"
23 #if OHMMS_DIM == 3
24 #include "ParticleIO/ESHDFParticleParser.h"
25 #endif
26 #include "Utilities/ProgressReportEngine.h"
27 #include "OhmmsData/AttributeSet.h"
28 #include "OhmmsData/Libxml2Doc.h"
29 #include "Particle/InitMolecularSystem.h"
30 #include "LongRange/LRCoulombSingleton.h"
31
32 namespace qmcplusplus
33 {
ParticleSetPool(Communicate * c,const char * aname)34 ParticleSetPool::ParticleSetPool(Communicate* c, const char* aname) : MPIObjectBase(c), TileMatrix(0)
35 {
36 TileMatrix.diagonal(1);
37 ClassName = "ParticleSetPool";
38 myName = aname;
39 }
40
ParticleSetPool(ParticleSetPool && other)41 ParticleSetPool::ParticleSetPool(ParticleSetPool&& other)
42 : MPIObjectBase(other.myComm),
43 SimulationCell(std::move(other.SimulationCell)),
44 TileMatrix(other.TileMatrix),
45 myPool(std::move(other.myPool))
46 {
47 ClassName = other.ClassName;
48 myName = other.myName;
49 }
50
~ParticleSetPool()51 ParticleSetPool::~ParticleSetPool()
52 {
53 PoolType::const_iterator it(myPool.begin()), it_end(myPool.end());
54 while (it != it_end)
55 {
56 delete (*it).second;
57 it++;
58 }
59 }
60
getParticleSet(const std::string & pname)61 ParticleSet* ParticleSetPool::getParticleSet(const std::string& pname)
62 {
63 std::map<std::string, ParticleSet*>::iterator pit(myPool.find(pname));
64 if (pit == myPool.end())
65 {
66 return 0;
67 }
68 else
69 {
70 return (*pit).second;
71 }
72 }
73
getWalkerSet(const std::string & pname)74 MCWalkerConfiguration* ParticleSetPool::getWalkerSet(const std::string& pname)
75 {
76 ParticleSet* mc = 0;
77 if (myPool.size() == 1)
78 mc = (*myPool.begin()).second;
79 else
80 mc = getParticleSet(pname);
81 if (mc == 0)
82 {
83 APP_ABORT("ParticleSePool::getWalkerSet missing " + pname);
84 }
85 return dynamic_cast<MCWalkerConfiguration*>(mc);
86 }
87
addParticleSet(std::unique_ptr<ParticleSet> && p)88 void ParticleSetPool::addParticleSet(std::unique_ptr<ParticleSet>&& p)
89 {
90 PoolType::iterator pit(myPool.find(p->getName()));
91 if (pit == myPool.end())
92 {
93 auto& pname = p->getName();
94 LOGMSG(" Adding " << pname << " ParticleSet to the pool")
95 myPool[pname] = p.release();
96 }
97 else
98 {
99 WARNMSG(" " << p->getName() << " exists. Ignore addition")
100 }
101 }
102
putTileMatrix(xmlNodePtr cur)103 bool ParticleSetPool::putTileMatrix(xmlNodePtr cur)
104 {
105 TileMatrix = 0;
106 TileMatrix.diagonal(1);
107 OhmmsAttributeSet pAttrib;
108 pAttrib.add(TileMatrix, "tilematrix");
109 pAttrib.put(cur);
110 return true;
111 }
112
putLattice(xmlNodePtr cur)113 bool ParticleSetPool::putLattice(xmlNodePtr cur)
114 {
115 ReportEngine PRE("ParticleSetPool", "putLattice");
116 bool printcell = false;
117 if (!SimulationCell)
118 {
119 app_debug() << " Creating global supercell " << std::endl;
120 SimulationCell = std::make_unique<ParticleSet::ParticleLayout_t>();
121 printcell = true;
122 }
123 else
124 {
125 app_log() << " Overwriting global supercell " << std::endl;
126 }
127 LatticeParser a(*SimulationCell);
128 bool lattice_defined = a.put(cur);
129 if (printcell && lattice_defined)
130 {
131 if (outputManager.isHighActive())
132 {
133 SimulationCell->print(app_log(), 2);
134 }
135 else
136 {
137 SimulationCell->print(app_summary(), 1);
138 }
139 }
140 return lattice_defined;
141 }
142
143 /** process an xml element
144 * @param cur current xmlNodePtr
145 * @return true, if successful.
146 *
147 * Creating MCWalkerConfiguration for all the ParticleSet
148 * objects.
149 */
put(xmlNodePtr cur)150 bool ParticleSetPool::put(xmlNodePtr cur)
151 {
152 ReportEngine PRE("ParticleSetPool", "put");
153 //const ParticleSet::ParticleLayout_t* sc=DistanceTable::getSimulationCell();
154 //ParticleSet::ParticleLayout_t* sc=0;
155 std::string id("e");
156 std::string role("none");
157 std::string randomR("no");
158 std::string randomsrc;
159 std::string useGPU;
160 OhmmsAttributeSet pAttrib;
161 pAttrib.add(id, "id");
162 pAttrib.add(id, "name");
163 pAttrib.add(role, "role");
164 pAttrib.add(randomR, "random");
165 pAttrib.add(randomsrc, "randomsrc");
166 pAttrib.add(randomsrc, "random_source");
167 #if defined(ENABLE_OFFLOAD)
168 pAttrib.add(useGPU, "gpu", {"yes", "no"});
169 #endif
170 pAttrib.put(cur);
171 //backward compatibility
172 if (id == "e" && role == "none")
173 role = "MC";
174 ParticleSet* pTemp = getParticleSet(id);
175 if (pTemp == 0)
176 {
177 app_summary() << std::endl;
178 app_summary() << " Particle Set" << std::endl;
179 app_summary() << " ------------" << std::endl;
180 app_summary() << " Name: " << id << " Offload : " << useGPU << std::endl;
181 app_summary() << std::endl;
182
183 // select OpenMP offload implementation in ParticleSet.
184 if (useGPU == "yes")
185 pTemp = new MCWalkerConfiguration(DynamicCoordinateKind::DC_POS_OFFLOAD);
186 else
187 pTemp = new MCWalkerConfiguration(DynamicCoordinateKind::DC_POS);
188 //if(role == "MC")
189 // pTemp = new MCWalkerConfiguration;
190 //else
191 // pTemp = new ParticleSet;
192 if (SimulationCell)
193 {
194 app_log() << " Initializing the lattice by the global supercell" << std::endl;
195 pTemp->Lattice = *SimulationCell;
196 }
197 myPool[id] = pTemp;
198 XMLParticleParser pread(*pTemp, TileMatrix);
199 bool success = pread.put(cur);
200 //if random_source is given, create a node <init target="" soruce=""/>
201 if (randomR == "yes" && !randomsrc.empty())
202 {
203 xmlNodePtr anode = xmlNewNode(NULL, (const xmlChar*)"init");
204 xmlNewProp(anode, (const xmlChar*)"source", (const xmlChar*)randomsrc.c_str());
205 xmlNewProp(anode, (const xmlChar*)"target", (const xmlChar*)id.c_str());
206 randomize_nodes.push_back(anode);
207 }
208 pTemp->setName(id);
209 app_summary() << " Particle set size: " << pTemp->getTotalNum() << std::endl;
210 app_summary() << std::endl;
211 return success;
212 }
213 else
214 {
215 app_warning() << "Particle set " << id << " is already created. Ignoring this section." << std::endl;
216 }
217 app_summary() << std::endl;
218 return true;
219 }
220
randomize()221 void ParticleSetPool::randomize()
222 {
223 app_log() << "ParticleSetPool::randomize " << randomize_nodes.size() << " ParticleSet"
224 << (randomize_nodes.size() == 1 ? "" : "s") << "." << std::endl;
225 bool success = true;
226 for (int i = 0; i < randomize_nodes.size(); ++i)
227 {
228 InitMolecularSystem moinit(*this);
229 success &= moinit.put(randomize_nodes[i]);
230 xmlFreeNode(randomize_nodes[i]);
231 }
232 randomize_nodes.clear();
233 if (!success)
234 APP_ABORT("ParticleSePool::randomize failed to randomize some Particlesets!");
235 }
236
get(std::ostream & os) const237 bool ParticleSetPool::get(std::ostream& os) const
238 {
239 os << "ParticleSetPool has: " << std::endl << std::endl;
240 os.setf(std::ios::scientific, std::ios::floatfield);
241 os.precision(14);
242 PoolType::const_iterator it(myPool.begin()), it_end(myPool.end());
243 while (it != it_end)
244 {
245 (*it).second->get(os);
246 ++it;
247 }
248 return true;
249 }
250
output_particleset_info(Libxml2Document & doc,xmlNodePtr root)251 void ParticleSetPool::output_particleset_info(Libxml2Document& doc, xmlNodePtr root)
252 {
253 xmlNodePtr particles_info = doc.addChild(root, "particles");
254 PoolType::const_iterator it(myPool.begin()), it_end(myPool.end());
255 while (it != it_end)
256 {
257 xmlNodePtr particle = doc.addChild(particles_info, "particle");
258 doc.addChild(particle, "name", (*it).second->getName());
259 doc.addChild(particle, "size", (*it).second->getTotalNum());
260 ++it;
261 }
262 }
263
264 /** reset is used to initialize and evaluate the distance tables
265 */
reset()266 void ParticleSetPool::reset()
267 {
268 PoolType::iterator it(myPool.begin()), it_end(myPool.end());
269 while (it != it_end)
270 {
271 ParticleSet* pt((*it).second);
272 pt->update();
273 ++it;
274 }
275 }
276
277 /** Create particlesets from ES-HDF file
278 */
createESParticleSet(xmlNodePtr cur,const std::string & target,ParticleSet * qp)279 ParticleSet* ParticleSetPool::createESParticleSet(xmlNodePtr cur, const std::string& target, ParticleSet* qp)
280 {
281 //TinyVector<int,OHMMS_DIM> tilefactor;
282 Tensor<int, OHMMS_DIM> eshdf_tilematrix(0);
283 eshdf_tilematrix.diagonal(1);
284 double lr_cut = 15;
285 std::string h5name;
286 std::string source("i");
287 std::string bc("p p p");
288 std::string spotype("0");
289 std::string lr_handler("opt_breakup");
290 OhmmsAttributeSet attribs;
291 attribs.add(h5name, "href");
292 attribs.add(eshdf_tilematrix, "tilematrix");
293 attribs.add(source, "source");
294 attribs.add(bc, "bconds");
295 attribs.add(lr_cut, "LR_dim_cutoff");
296 attribs.add(lr_handler, "LR_handler");
297 attribs.add(spotype, "type");
298 attribs.put(cur);
299
300 if (h5name.empty())
301 return qp;
302
303 #if OHMMS_DIM == 3
304 ParticleSet* ions = getParticleSet(source);
305 if (ions == 0)
306 {
307 ions = new MCWalkerConfiguration;
308 ions->setName(source);
309 //set the boundary condition
310 ions->Lattice.LR_dim_cutoff = lr_cut;
311 std::istringstream is(bc);
312 char c;
313 int idim = 0;
314 while (!is.eof() && idim < OHMMS_DIM)
315 {
316 if (is >> c)
317 ions->Lattice.BoxBConds[idim++] = (c == 'p');
318 }
319
320 tolower(lr_handler);
321 if (lr_handler == "ewald")
322 {
323 LRCoulombSingleton::this_lr_type = LRCoulombSingleton::EWALD;
324 }
325 else if (lr_handler == "opt_breakup")
326 {
327 LRCoulombSingleton::this_lr_type = LRCoulombSingleton::ESLER;
328 }
329 else if (lr_handler == "opt_breakup_original")
330 {
331 LRCoulombSingleton::this_lr_type = LRCoulombSingleton::NATOLI;
332 }
333 else
334 {
335 APP_ABORT("Long range breakup handler not recognized\n");
336 }
337
338 //initialize ions from hdf5
339 hid_t h5 = -1;
340 if (myComm->rank() == 0)
341 {
342 //Rather than turn off all H5errors, we're going to
343 //temporarily disable it.
344 //
345 //old_func: function pointer to current function that
346 // displays when H5 encounters error.
347 H5E_auto2_t old_func;
348 //old_client_data: null pointer to associated error stream.
349 void* old_client_data;
350 //Grab the current handler info.
351 H5Eget_auto2(H5E_DEFAULT, &old_func, &old_client_data);
352 //Now kill error notifications.
353 H5Eset_auto2(H5E_DEFAULT, NULL, NULL);
354 h5 = H5Fopen(h5name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
355 //and restore to defaults.
356 H5Eset_auto2(H5E_DEFAULT, old_func, old_client_data);
357 if (h5 < 0)
358 {
359 app_error() << "Could not open HDF5 file \"" << h5name
360 << "\" in ParticleSetPool::createESParticleSet. Aborting.\n"
361 << "(Please ensure that your path is correct, the file exists, and that "
362 << "you have read permissions.)\n";
363 APP_ABORT("ParticleSetPool::createESParticleSet");
364 }
365 }
366 ESHDFIonsParser ap(*ions, h5, myComm);
367 ap.put(cur);
368 ap.expand(eshdf_tilematrix);
369 if (h5 > -1)
370 H5Fclose(h5);
371 //failed to initialize the ions
372 if (ions->getTotalNum() == 0)
373 return 0;
374 typedef ParticleSet::SingleParticleIndex_t SingleParticleIndex_t;
375 std::vector<SingleParticleIndex_t> grid(OHMMS_DIM, SingleParticleIndex_t(1));
376 ions->Lattice.reset();
377
378 myPool[source] = ions;
379 }
380
381 if (SimulationCell == 0)
382 {
383 SimulationCell = std::make_unique<ParticleSet::ParticleLayout_t>(ions->Lattice);
384 }
385
386 if (qp == 0)
387 {
388 //create the electrons
389 qp = new MCWalkerConfiguration;
390 qp->setName(target);
391 qp->Lattice = ions->Lattice;
392
393 app_log() << " Simulation cell radius = " << qp->Lattice.SimulationCellRadius << std::endl;
394 app_log() << " Wigner-Seitz cell radius = " << qp->Lattice.WignerSeitzRadius << std::endl;
395 SimulationCell->print(app_log());
396
397 // Goback to the // and OhmmsXPathObject handles it internally
398 OhmmsXPathObject det("//determinant", cur);
399
400 if (det.size() > 2)
401 APP_ABORT("Only two electron groups are supported.");
402
403 std::vector<int> num_spin(det.size(), 0);
404 for (int i = 0; i < det.size(); ++i)
405 {
406 OhmmsAttributeSet a;
407 a.add(num_spin[i], "size");
408 a.put(det[i]);
409 }
410
411 {
412 //create species
413 SpeciesSet& species = qp->getSpeciesSet();
414 //add up and down
415 species.addSpecies("u");
416 if (num_spin.size() > 1)
417 species.addSpecies("d");
418 int chid = species.addAttribute("charge");
419 for (int i = 0; i < num_spin.size(); ++i)
420 species(chid, i) = -1.0;
421 int mid = species.addAttribute("membersize");
422 for (int i = 0; i < num_spin.size(); ++i)
423 species(mid, i) = num_spin[i];
424 mid = species.addAttribute("mass");
425 for (int i = 0; i < num_spin.size(); ++i)
426 species(mid, i) = 1.0;
427 qp->create(num_spin);
428 }
429 //name it with the target
430 qp->setName(target);
431
432 //if(qp->getTotalNum() == 0 || ions->getTotalNum() == 0)
433 //{
434 // delete qp;
435 // delete ions;
436 // APP_ABORT("ParticleSetPool failed to create particlesets for the electron structure calculation");
437 // return 0;
438 //}
439 //for PPP, use uniform random
440 if (qp->Lattice.SuperCellEnum == SUPERCELL_BULK)
441 {
442 makeUniformRandom(qp->R);
443 qp->R.setUnit(PosUnit::Lattice);
444 qp->convert2Cart(qp->R);
445 }
446 else
447 {
448 //assign non-trivial positions for the quanmtum particles
449 InitMolecularSystem mole(*this);
450 mole.initMolecule(ions, qp);
451 qp->R.setUnit(PosUnit::Cartesian);
452 }
453 //for(int i=0; i<qp->getTotalNum(); ++i)
454 // std::cout << qp->GroupID[i] << " " << qp->R[i] << std::endl;
455
456 #if !defined(QMC_CUDA)
457 makeUniformRandom(qp->spins);
458 qp->spins *= 2 * M_PI;
459 #endif
460
461 if (qp->Lattice.SuperCellEnum)
462 qp->createSK();
463 qp->resetGroups();
464 myPool[target] = qp;
465 }
466
467 #endif
468 return qp;
469 }
470 } // namespace qmcplusplus
471