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: Jordan E. Vincent, University of Illinois at Urbana-Champaign
8 //                    Bryan Clark, bclark@Princeton.edu, Princeton University
9 //                    Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
10 //                    Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
11 //                    Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
12 //                    Cynthia Gu, zg1@ornl.gov, Oak Ridge National Laboratory
13 //                    Ye Luo, yeluo@anl.gov, Argonne National Laboratory
14 //                    Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
15 //
16 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
17 //////////////////////////////////////////////////////////////////////////////////////
18 
19 
20 #include "MCWalkerConfiguration.h"
21 #include "Particle/DistanceTableData.h"
22 #include "ParticleBase/RandomSeqGenerator.h"
23 #include "Message/Communicate.h"
24 #include "Message/CommOperators.h"
25 #include "Utilities/IteratorUtility.h"
26 #include "LongRange/StructFact.h"
27 #include "Particle/HDFWalkerOutput.h"
28 #include "Particle/MCSample.h"
29 #include "hdf/hdf_hyperslab.h"
30 #include "hdf/HDFVersion.h"
31 #include <map>
32 
33 #ifdef QMC_CUDA
34 #include "Particle/accept_kernel.h"
35 #endif
36 
37 namespace qmcplusplus
38 {
MCWalkerConfiguration(const DynamicCoordinateKind kind)39 MCWalkerConfiguration::MCWalkerConfiguration(const DynamicCoordinateKind kind)
40     : ParticleSet(kind),
41 #ifdef QMC_CUDA
42       RList_GPU("MCWalkerConfiguration::RList_GPU"),
43       GradList_GPU("MCWalkerConfiguration::GradList_GPU"),
44       LapList_GPU("MCWalkerConfiguration::LapList_GPU"),
45       Rnew_GPU("MCWalkerConfiguration::Rnew_GPU"),
46       NLlist_GPU("MCWalkerConfiguration::NLlist_GPU"),
47       iatList_GPU("iatList_GPU"),
48       AcceptList_GPU("MCWalkerConfiguration::AcceptList_GPU"),
49 #endif
50       ReadyForPbyP(false),
51       UpdateMode(Update_Walker),
52       reptile(0),
53       Polymer(0)
54 {
55   //move to ParticleSet
56   //initPropertyList();
57 }
58 
MCWalkerConfiguration(const MCWalkerConfiguration & mcw)59 MCWalkerConfiguration::MCWalkerConfiguration(const MCWalkerConfiguration& mcw)
60     : ParticleSet(mcw),
61 #ifdef QMC_CUDA
62       RList_GPU("MCWalkerConfiguration::RList_GPU"),
63       GradList_GPU("MCWalkerConfiguration::GradList_GPU"),
64       LapList_GPU("MCWalkerConfiguration::LapList_GPU"),
65       Rnew_GPU("MCWalkerConfiguration::Rnew_GPU"),
66       NLlist_GPU("MCWalkerConfiguration::NLlist_GPU"),
67       iatList_GPU("iatList_GPU"),
68       AcceptList_GPU("MCWalkerConfiguration::AcceptList_GPU"),
69 #endif
70       ReadyForPbyP(false),
71       UpdateMode(Update_Walker),
72       Polymer(0)
73 {
74   samples.clearEnsemble();
75   samples.setMaxSamples(mcw.getMaxSamples());
76   GlobalNumWalkers = mcw.GlobalNumWalkers;
77   WalkerOffsets    = mcw.WalkerOffsets;
78   Properties       = mcw.Properties;
79   //initPropertyList();
80 }
81 
createWalkers(int n)82 void MCWalkerConfiguration::createWalkers(int n)
83 {
84   const int old_nw = getActiveWalkers();
85   WalkerConfigurations::createWalkers(n, TotalNum);
86   // no pre-existing walkers, need to initialized based on particleset.
87   if (old_nw == 0)
88     for (auto* awalker : WalkerList)
89     {
90       awalker->R     = R;
91       awalker->spins = spins;
92     }
93   resizeWalkerHistories();
94 }
95 
96 
resize(int numWalkers,int numPtcls)97 void MCWalkerConfiguration::resize(int numWalkers, int numPtcls)
98 {
99   if (TotalNum && WalkerList.size())
100     app_warning() << "MCWalkerConfiguration::resize cleans up the walker list." << std::endl;
101   const int old_nw = getActiveWalkers();
102   ParticleSet::resize(unsigned(numPtcls));
103   WalkerConfigurations::resize(numWalkers, TotalNum);
104   // no pre-existing walkers, need to initialized based on particleset.
105   if (old_nw == 0)
106     for (auto* awalker : WalkerList)
107     {
108       awalker->R     = R;
109       awalker->spins = spins;
110     }
111 }
112 
113 /** Make Metropolis move to the walkers and save in a temporary array.
114  * @param it the iterator of the first walker to work on
115  * @param tauinv  inverse of the time step
116  *
117  * R + D + X
118  */
sample(iterator it,RealType tauinv)119 void MCWalkerConfiguration::sample(iterator it, RealType tauinv)
120 {
121   APP_ABORT("MCWalkerConfiguration::sample obsolete");
122   //  makeGaussRandom(R);
123   //  R *= tauinv;
124   //  R += (*it)->R + (*it)->Drift;
125 }
126 
127 //void MCWalkerConfiguration::clearAuxDataSet() {
128 //  UpdateMode=Update_Particle;
129 //  int nbytes=128*TotalNum*sizeof(RealType);//could be pagesize
130 //  if(WalkerList.size())//check if capacity is bigger than the estimated one
131 //    nbytes = (WalkerList[0]->DataSet.capacity()>nbytes)?WalkerList[0]->DataSet.capacity():nbytes;
132 //  iterator it(WalkerList.begin());
133 //  iterator it_end(WalkerList.end());
134 //  while(it!=it_end) {
135 //    (*it)->DataSet.clear();
136 //    //CHECK THIS WITH INTEL 10.1
137 //    //(*it)->DataSet.reserve(nbytes);
138 //    ++it;
139 //  }
140 //  ReadyForPbyP = true;
141 //}
142 //
143 //bool MCWalkerConfiguration::createAuxDataSet(int nfield) {
144 //
145 //  if(ReadyForPbyP) return false;
146 //
147 //  ReadyForPbyP=true;
148 //  UpdateMode=Update_Particle;
149 //  iterator it(WalkerList.begin());
150 //  iterator it_end(WalkerList.end());
151 //  while(it!=it_end) {
152 //    (*it)->DataSet.reserve(nfield); ++it;
153 //  }
154 //
155 //  return true;
156 //}
157 
158 
159 /** reset the Property container of all the walkers
160  */
resetWalkerProperty(int ncopy)161 void MCWalkerConfiguration::resetWalkerProperty(int ncopy)
162 {
163   int m(PropertyList.size());
164   app_log() << "  Resetting Properties of the walkers " << ncopy << " x " << m << std::endl;
165   try
166   {
167     Properties.resize(ncopy, m);
168   }
169   catch (std::domain_error& de)
170   {
171     app_error() << de.what() << '\n'
172                 << "This is likely because some object has attempted to add walker properties\n"
173                 << " in excess of WALKER_MAX_PROPERTIES.\n"
174                 << "build with cmake ... -DWALKER_MAX_PROPERTIES=at_least_properties_required" << std::endl;
175     APP_ABORT("Fatal Exception");
176   }
177 
178   iterator it(WalkerList.begin()), it_end(WalkerList.end());
179   while (it != it_end)
180   {
181     (*it)->resizeProperty(ncopy, m);
182     (*it)->Weight = 1;
183     ++it;
184   }
185   resizeWalkerHistories();
186 }
187 
resizeWalkerHistories()188 void MCWalkerConfiguration::resizeWalkerHistories()
189 {
190   //using std::vector<std::vector<RealType> > is too costly.
191   int np = PropertyHistory.size();
192   if (np)
193     for (int iw = 0; iw < WalkerList.size(); ++iw)
194       WalkerList[iw]->PropertyHistory = PropertyHistory;
195   np = PHindex.size();
196   if (np)
197     for (int iw = 0; iw < WalkerList.size(); ++iw)
198       WalkerList[iw]->PHindex = PHindex;
199   ;
200 }
201 
202 /** allocate the SampleStack
203  * @param n number of samples per thread
204  */
setNumSamples(int n)205 void MCWalkerConfiguration::setNumSamples(int n)
206 {
207   samples.clearEnsemble();
208   samples.setMaxSamples(n);
209 }
210 
211 /** save the current walkers to SampleStack
212  */
saveEnsemble()213 void MCWalkerConfiguration::saveEnsemble() { saveEnsemble(WalkerList.begin(), WalkerList.end()); }
214 
215 /** save the [first,last) walkers to SampleStack
216  */
saveEnsemble(iterator first,iterator last)217 void MCWalkerConfiguration::saveEnsemble(iterator first, iterator last)
218 {
219   for (; first != last; first++)
220   {
221     samples.appendSample(MCSample(**first));
222   }
223 }
224 /** load a single sample from SampleStack
225  */
loadSample(ParticleSet & pset,size_t iw) const226 void MCWalkerConfiguration::loadSample(ParticleSet& pset, size_t iw) const { samples.loadSample(pset, iw); }
227 
228 /** load SampleStack to WalkerList
229  */
loadEnsemble()230 void MCWalkerConfiguration::loadEnsemble()
231 {
232   using WP     = WalkerProperties::Indexes;
233   int nsamples = std::min(samples.getMaxSamples(), samples.getNumSamples());
234   if (samples.empty() || nsamples == 0)
235     return;
236   Walker_t::PropertyContainer_t prop(1, PropertyList.size(), 1, WP::MAXPROPERTIES);
237   delete_iter(WalkerList.begin(), WalkerList.end());
238   WalkerList.resize(nsamples);
239   for (int i = 0; i < nsamples; ++i)
240   {
241     Walker_t* awalker = new Walker_t(TotalNum);
242     awalker->Properties.copy(prop);
243     samples.getSample(i).convertToWalker(*awalker);
244     WalkerList[i] = awalker;
245   }
246   resizeWalkerHistories();
247   samples.clearEnsemble();
248 }
249 
250 ///** load SampleStack to WalkerList
251 // */
252 //void MCWalkerConfiguration::loadEnsemble(const Walker_t& wcopy)
253 //{
254 //  int nsamples=std::min(MaxSamples,CurSampleCount);
255 //  if(SampleStack.empty() || nsamples==0) return;
256 //
257 //  Walker_t::PropertyContainer_t prop(1,PropertyList.size());
258 //
259 //  while(WalkerList.size()) pop_back();
260 //  WalkerList.resize(nsamples);
261 //
262 //  for(int i=0; i<nsamples; ++i)
263 //  {
264 //    Walker_t* awalker=new Walker_t(TotalNum);
265 //    awalker->Properties.copy(prop);
266 //    SampleStack[i]->get(*awalker);
267 //    WalkerList[i]=awalker;
268 //  }
269 //  resizeWalkerHistories();
270 //  clearEnsemble();
271 //}
272 //
273 //void MCWalkerConfiguration::loadEnsemble(MCWalkerConfiguration& other)
274 //{
275 //  if(SampleStack.empty()) return;
276 //
277 //  Walker_t twalker(*WalkerList[0]);
278 //  for(int i=0; i<MaxSamples; ++i)
279 //  {
280 //    Walker_t* awalker=new Walker_t(twalker);
281 //    SampleStack[i]->get(*awalker);
282 //    other.WalkerList.push_back(awalker);
283 //  }
284 //
285 //  clearEnsemble();
286 //}
287 
dumpEnsemble(std::vector<MCWalkerConfiguration * > & others,HDFWalkerOutput * out,int np,int nBlock)288 bool MCWalkerConfiguration::dumpEnsemble(std::vector<MCWalkerConfiguration*>& others,
289                                          HDFWalkerOutput* out,
290                                          int np,
291                                          int nBlock)
292 {
293   return samples.dumpEnsemble(others, out, np, nBlock);
294 }
295 
getMaxSamples() const296 int MCWalkerConfiguration::getMaxSamples() const { return samples.getMaxSamples(); }
297 
loadEnsemble(std::vector<MCWalkerConfiguration * > & others,bool doclean)298 void MCWalkerConfiguration::loadEnsemble(std::vector<MCWalkerConfiguration*>& others, bool doclean)
299 {
300   using WP = WalkerProperties::Indexes;
301   std::vector<int> off(others.size() + 1, 0);
302   for (int i = 0; i < others.size(); ++i)
303   {
304     off[i + 1] = off[i] + std::min(others[i]->getMaxSamples(), others[i]->numSamples());
305   }
306   int nw_tot = off.back();
307   if (nw_tot)
308   {
309     Walker_t::PropertyContainer_t prop(1, PropertyList.size(), 1, WP::MAXPROPERTIES);
310     while (WalkerList.size())
311       pop_back();
312     WalkerList.resize(nw_tot);
313     for (int i = 0; i < others.size(); ++i)
314     {
315       SampleStack& astack(others[i]->getSampleStack());
316       for (int j = 0, iw = off[i]; iw < off[i + 1]; ++j, ++iw)
317       {
318         Walker_t* awalker = new Walker_t(TotalNum);
319         awalker->Properties.copy(prop);
320         astack.getSample(j).convertToWalker(*awalker);
321         WalkerList[iw] = awalker;
322       }
323       if (doclean)
324         others[i]->clearEnsemble();
325     }
326   }
327   if (doclean)
328     resizeWalkerHistories();
329 }
330 
clearEnsemble()331 void MCWalkerConfiguration::clearEnsemble() { samples.clearEnsemble(); }
332 
333 #ifdef QMC_CUDA
updateLists_GPU()334 void MCWalkerConfiguration::updateLists_GPU()
335 {
336   int nw         = WalkerList.size();
337   int NumSpecies = getSpeciesSet().TotalNum;
338   if (Rnew_GPU.size() != nw * kblocksize)
339   {
340     Rnew_GPU.resize(nw * kblocksize);
341     RhokLists_GPU.resize(NumSpecies);
342     for (int isp = 0; isp < NumSpecies; isp++)
343       RhokLists_GPU[isp].resize(nw);
344     Rnew_host.resize(nw * kblocksize);
345     Rnew.resize(nw * kblocksize);
346     AcceptList_GPU.resize(nw);
347     AcceptList_host.resize(nw);
348     RList_GPU.resize(nw);
349     GradList_GPU.resize(nw);
350     LapList_GPU.resize(nw);
351     DataList_GPU.resize(nw);
352   }
353   hostlist.resize(nw);
354   hostlist_valueType.resize(nw);
355   hostlist_AA.resize(nw);
356 
357   for (int iw = 0; iw < nw; iw++)
358   {
359     if (WalkerList[iw]->R_GPU.size() != R.size())
360       std::cerr << "Error in R_GPU size for iw = " << iw << "!\n";
361     hostlist[iw] = (CTS::RealType*)WalkerList[iw]->R_GPU.data();
362   }
363   RList_GPU = hostlist;
364 
365   for (int iw = 0; iw < nw; iw++)
366   {
367     if (WalkerList[iw]->Grad_GPU.size() != R.size())
368       std::cerr << "Error in Grad_GPU size for iw = " << iw << "!\n";
369     hostlist_valueType[iw] = (CTS::ValueType*)WalkerList[iw]->Grad_GPU.data();
370   }
371   GradList_GPU = hostlist_valueType;
372 
373   for (int iw = 0; iw < nw; iw++)
374   {
375     if (WalkerList[iw]->Lap_GPU.size() != R.size())
376       std::cerr << "Error in Lap_GPU size for iw = " << iw << "!\n";
377     hostlist_valueType[iw] = (CTS::ValueType*)WalkerList[iw]->Lap_GPU.data();
378   }
379   LapList_GPU = hostlist_valueType;
380 
381   for (int iw = 0; iw < nw; iw++)
382     hostlist_valueType[iw] = WalkerList[iw]->cuda_DataSet.data();
383   DataList_GPU = hostlist_valueType;
384 
385   for (int isp = 0; isp < NumSpecies; isp++)
386   {
387     for (int iw = 0; iw < nw; iw++)
388       hostlist_AA[iw] = WalkerList[iw]->get_rhok_ptr(isp);
389     RhokLists_GPU[isp] = hostlist_AA;
390   }
391 }
392 
allocateGPU(size_t buffersize)393 void MCWalkerConfiguration::allocateGPU(size_t buffersize)
394 {
395   int N    = WalkerList[0]->R.size();
396   int Numk = 0;
397   if (SK)
398     Numk = SK->KLists.numk;
399   int NumSpecies = getSpeciesSet().TotalNum;
400   for (int iw = 0; iw < WalkerList.size(); iw++)
401   {
402     Walker_t& walker = *(WalkerList[iw]);
403     walker.resizeCuda(buffersize, NumSpecies, Numk);
404   }
405 }
406 
407 
copyWalkersToGPU(bool copyGrad)408 void MCWalkerConfiguration::copyWalkersToGPU(bool copyGrad)
409 {
410   R_host.resize(WalkerList[0]->R.size());
411   for (int iw = 0; iw < WalkerList.size(); iw++)
412   {
413     for (int i = 0; i < WalkerList[iw]->size(); i++)
414       for (int dim = 0; dim < OHMMS_DIM; dim++)
415         R_host[i][dim] = WalkerList[iw]->R[i][dim];
416     WalkerList[iw]->R_GPU = R_host;
417   }
418   if (copyGrad)
419     copyWalkerGradToGPU();
420 }
421 
copyWalkerGradToGPU()422 void MCWalkerConfiguration::copyWalkerGradToGPU()
423 {
424   Grad_host.resize(WalkerList[0]->G.size());
425   for (int iw = 0; iw < WalkerList.size(); iw++)
426   {
427     for (int i = 0; i < WalkerList[iw]->size(); i++)
428       for (int dim = 0; dim < OHMMS_DIM; dim++)
429         Grad_host[i][dim] = WalkerList[iw]->G[i][dim];
430     WalkerList[iw]->Grad_GPU = Grad_host;
431   }
432 }
433 
proposeMove_GPU(std::vector<PosType> & newPos,int iat)434 void MCWalkerConfiguration::proposeMove_GPU(std::vector<PosType>& newPos, int iat)
435 {
436   int nw = newPos.size();
437   if (Rnew_host.size() < nw * kblocksize)
438   {
439     Rnew.resize(nw * kblocksize);
440     Rnew_host.resize(nw * kblocksize);
441   }
442   // store things sequentially with k to make evaluation more straight-forward:
443   //           k=0     k=1     k=kblocksize
444   // Rnew = [0,..,nw|0,..,nw|...|0,..,nw]
445   int offset = kcurr * nw;
446   for (int i = 0; i < nw; i++)
447   {
448     for (int dim = 0; dim < OHMMS_DIM; dim++)
449     {
450       Rnew[i + offset][dim]      = newPos[i][dim];
451       Rnew_host[i + offset][dim] = newPos[i][dim];
452     }
453   }
454   if (kDelay)
455   {
456     kcurr  = (kcurr + 1) % kblocksize; // loop kcurr around every k blocks
457     kstart = kblock * kblocksize;
458     if (kcurr == 0)
459       kblock++; // keep increasing kblock (even beyond available matrix blocks) - the update check takes care of self-consistency
460     // only copy new position matrix when needed (when update is imminent)
461     if (klinear)
462     {
463       Rnew_GPU.asyncCopy(&(Rnew_host[offset]), nw * kblocksize, offset, nw);
464     }
465     else if (kcurr == 0 || (kcurr + kblock * kblocksize >= getnat(iat)))
466       Rnew_GPU.asyncCopy(Rnew_host);
467   }
468   else
469     Rnew_GPU.asyncCopy(Rnew_host);
470   CurrentParticle = iat;
471 }
472 
473 
acceptMove_GPU(std::vector<bool> & toAccept,int k)474 void MCWalkerConfiguration::acceptMove_GPU(std::vector<bool>& toAccept, int k)
475 {
476   if (AcceptList_host.size() < toAccept.size())
477     AcceptList_host.resize(toAccept.size());
478   for (int i = 0; i < toAccept.size(); i++)
479     AcceptList_host[i] = (int)toAccept[i];
480   AcceptList_GPU.asyncCopy(AcceptList_host);
481   //   app_log() << "toAccept.size()        = " << toAccept.size() << std::endl;
482   //   app_log() << "AcceptList_host.size() = " << AcceptList_host.size() << std::endl;
483   //   app_log() << "AcceptList_GPU.size()  = " << AcceptList_GPU.size() << std::endl;
484   //   app_log() << "WalkerList.size()      = " << WalkerList.size() << std::endl;
485   //   app_log() << "Rnew_GPU.size()        = " << Rnew_GPU.size() << std::endl;
486   //   app_log() << "RList_GPU.size()       = " << RList_GPU.size() << std::endl;
487   if (RList_GPU.size() != WalkerList.size())
488     std::cerr << "Error in RList_GPU size.\n";
489   if (Rnew_GPU.size() != WalkerList.size() * kblocksize)
490     std::cerr << "Error in Rnew_GPU size.\n";
491   if (AcceptList_GPU.size() != WalkerList.size())
492     std::cerr << "Error in AcceptList_GPU size.\n";
493   accept_move_GPU_cuda(RList_GPU.data(), (CUDA_PRECISION*)Rnew_GPU.data(), AcceptList_GPU.data(), CurrentParticle++,
494                        WalkerList.size(), k);
495 }
496 
497 
NLMove_GPU(std::vector<Walker_t * > & walkers,std::vector<PosType> & newpos,std::vector<int> & iat)498 void MCWalkerConfiguration::NLMove_GPU(std::vector<Walker_t*>& walkers,
499                                        std::vector<PosType>& newpos,
500                                        std::vector<int>& iat)
501 {
502   int N = walkers.size();
503   if (NLlist_GPU.size() < N)
504   {
505     NLlist_GPU.resize(N);
506     NLlist_host.resize(N);
507   }
508   if (Rnew_GPU.size() < N)
509   {
510     Rnew_host.resize(N);
511     Rnew_GPU.resize(N);
512   }
513   for (int iw = 0; iw < N; iw++)
514   {
515     Rnew_host[iw]   = newpos[iw];
516     NLlist_host[iw] = (CUDA_PRECISION*)(walkers[iw]->R_GPU.data()) + OHMMS_DIM * iat[iw];
517   }
518   Rnew_GPU   = Rnew_host;
519   NLlist_GPU = NLlist_host;
520   NL_move_cuda(NLlist_GPU.data(), (CUDA_PRECISION*)Rnew_GPU.data(), N);
521 }
522 #endif
523 
524 } // namespace qmcplusplus
525