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