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: Ken Esler, kpesler@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 // Jaron T. Krogel, krogeljt@ornl.gov, Oak Ridge National Laboratory
11 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
12 //
13 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
14 //////////////////////////////////////////////////////////////////////////////////////
15
16
17 #include "VMCUpdatePbyP.h"
18 #include "QMCDrivers/DriftOperators.h"
19 #include "Message/OpenMP.h"
20 #if !defined(REMOVE_TRACEMANAGER)
21 #include "Estimators/TraceManager.h"
22 #else
23 typedef int TraceManager;
24 #endif
25
26
27 namespace qmcplusplus
28 {
29 /// Constructor
VMCUpdatePbyP(MCWalkerConfiguration & w,TrialWaveFunction & psi,QMCHamiltonian & h,RandomGenerator_t & rg)30 VMCUpdatePbyP::VMCUpdatePbyP(MCWalkerConfiguration& w, TrialWaveFunction& psi, QMCHamiltonian& h, RandomGenerator_t& rg)
31 : QMCUpdateBase(w, psi, h, rg),
32 buffer_timer_(*timer_manager.createTimer("VMCUpdatePbyP::Buffer", timer_level_medium)),
33 movepbyp_timer_(*timer_manager.createTimer("VMCUpdatePbyP::MovePbyP", timer_level_medium)),
34 hamiltonian_timer_(*timer_manager.createTimer("VMCUpdatePbyP::Hamiltonian", timer_level_medium)),
35 collectables_timer_(*timer_manager.createTimer("VMCUpdatePbyP::Collectables", timer_level_medium))
36 {}
37
~VMCUpdatePbyP()38 VMCUpdatePbyP::~VMCUpdatePbyP() {}
39
advanceWalker(Walker_t & thisWalker,bool recompute)40 void VMCUpdatePbyP::advanceWalker(Walker_t& thisWalker, bool recompute)
41 {
42 buffer_timer_.start();
43 W.loadWalker(thisWalker, true);
44 Walker_t::WFBuffer_t& w_buffer(thisWalker.DataSet);
45 Psi.copyFromBuffer(W, w_buffer);
46 buffer_timer_.stop();
47
48 // start PbyP moves
49 movepbyp_timer_.start();
50 bool moved = false;
51 constexpr RealType mhalf(-0.5);
52 for (int iter = 0; iter < nSubSteps; ++iter)
53 {
54 //create a 3N-Dimensional Gaussian with variance=1
55 makeGaussRandomWithEngine(deltaR, RandomGen);
56 moved = false;
57 for (int ig = 0; ig < W.groups(); ++ig) //loop over species
58 {
59 RealType tauovermass = Tau * MassInvS[ig];
60 RealType oneover2tau = 0.5 / (tauovermass);
61 RealType sqrttau = std::sqrt(tauovermass);
62 Psi.prepareGroup(W, ig);
63 for (int iat = W.first(ig); iat < W.last(ig); ++iat)
64 {
65 PosType dr;
66 if (UseDrift)
67 {
68 GradType grad_now = Psi.evalGrad(W, iat);
69 DriftModifier->getDrift(tauovermass, grad_now, dr);
70 dr += sqrttau * deltaR[iat];
71 }
72 else
73 dr = sqrttau * deltaR[iat];
74
75 if (!W.makeMoveAndCheck(iat, dr))
76 {
77 ++nReject;
78 W.accept_rejectMove(iat, false);
79 continue;
80 }
81
82 RealType prob(0);
83 if (UseDrift)
84 {
85 GradType grad_new;
86 prob = std::norm(Psi.calcRatioGrad(W, iat, grad_new));
87 DriftModifier->getDrift(tauovermass, grad_new, dr);
88 dr = W.R[iat] - W.activePos - dr;
89 RealType logGb = -oneover2tau * dot(dr, dr);
90 RealType logGf = mhalf * dot(deltaR[iat], deltaR[iat]);
91 prob *= std::exp(logGb - logGf);
92 }
93 else
94 prob = std::norm(Psi.calcRatio(W, iat));
95
96 bool is_accepted = false;
97 if (prob >= std::numeric_limits<RealType>::epsilon() && RandomGen() < prob)
98 {
99 is_accepted = true;
100 moved = true;
101 ++nAccept;
102 Psi.acceptMove(W, iat, true);
103 }
104 else
105 {
106 ++nReject;
107 Psi.rejectMove(iat);
108 }
109 W.accept_rejectMove(iat, is_accepted);
110 }
111 }
112 Psi.completeUpdates();
113 }
114 W.donePbyP();
115 movepbyp_timer_.stop();
116 buffer_timer_.start();
117 RealType logpsi = Psi.updateBuffer(W, w_buffer, recompute);
118 assert(checkLogAndGL(W, Psi));
119 W.saveWalker(thisWalker);
120 buffer_timer_.stop();
121 // end PbyP moves
122 hamiltonian_timer_.start();
123 FullPrecRealType eloc = H.evaluate(W);
124 thisWalker.resetProperty(logpsi, Psi.getPhase(), eloc);
125 hamiltonian_timer_.stop();
126 collectables_timer_.start();
127 H.auxHevaluate(W, thisWalker);
128 H.saveProperty(thisWalker.getPropertyBase());
129 collectables_timer_.stop();
130 #if !defined(REMOVE_TRACEMANAGER)
131 Traces->buffer_sample(W.current_step);
132 #endif
133 if (!moved)
134 ++nAllRejected;
135 }
136
137 } // namespace qmcplusplus
138