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