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: Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
9 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
10 //
11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
12 //////////////////////////////////////////////////////////////////////////////////////
13
14
15 #include "DMCUpdateAll.h"
16 #include "ParticleBase/ParticleUtility.h"
17 #include "ParticleBase/RandomSeqGenerator.h"
18 #include "QMCDrivers/DriftOperators.h"
19 #include "QMCDrivers/WalkerProperties.h"
20
21 namespace qmcplusplus
22 {
23 using WP = WalkerProperties::Indexes;
24
25 /// Constructor.
DMCUpdateAllWithRejection(MCWalkerConfiguration & w,TrialWaveFunction & psi,QMCHamiltonian & h,RandomGenerator_t & rg)26 DMCUpdateAllWithRejection::DMCUpdateAllWithRejection(MCWalkerConfiguration& w,
27 TrialWaveFunction& psi,
28 QMCHamiltonian& h,
29 RandomGenerator_t& rg)
30 : QMCUpdateBase(w, psi, h, rg)
31 {
32 UpdatePbyP = false;
33 }
34
35 /// destructor
~DMCUpdateAllWithRejection()36 DMCUpdateAllWithRejection::~DMCUpdateAllWithRejection() {}
37
38 /** advance all the walkers with killnode==no
39 * @param nat number of particles to move
40 *
41 * When killnode==no, any move resulting in node-crossing is treated
42 * as a normal rejection.
43 */
advanceWalker(Walker_t & thisWalker,bool recompute)44 void DMCUpdateAllWithRejection::advanceWalker(Walker_t& thisWalker, bool recompute)
45 {
46 W.loadWalker(thisWalker, false);
47 //create a 3N-Dimensional Gaussian with variance=1
48 RealType nodecorr = setScaledDriftPbyPandNodeCorr(Tau, MassInvP, W.G, drift);
49 //RealType nodecorr = setScaledDriftPbyPandNodeCorr(m_tauovermass,W.G,drift);
50 makeGaussRandomWithEngine(deltaR, RandomGen);
51 //save old local energy
52 RealType eold = thisWalker.Properties(WP::LOCALENERGY);
53 RealType enew = eold;
54 bool accepted = false;
55 RealType rr_accepted = 0.0;
56 RealType rr_proposed = 0.0;
57 RealType logpsi;
58
59 if (W.makeMoveAllParticlesWithDrift(thisWalker, drift, deltaR, SqrtTauOverMass))
60 {
61 //evaluate the new wave function
62 logpsi = Psi.evaluateLog(W);
63 //fixed node
64 if (!branchEngine->phaseChanged(Psi.getPhaseDiff()))
65 {
66 RealType logGf = -0.5 * Dot(deltaR, deltaR);
67 nodecorr = setScaledDriftPbyPandNodeCorr(Tau, MassInvP, W.G, drift);
68 deltaR = thisWalker.R - W.R - drift;
69 RealType logGb = logBackwardGF(deltaR);
70 //RealType logGb = -m_oneover2tau*Dot(deltaR,deltaR);
71 RealType prob = std::min(std::exp(logGb - logGf + 2.0 * (logpsi - thisWalker.Properties(WP::LOGPSI))), 1.0);
72 //calculate rr_proposed here
73 deltaR = W.R - thisWalker.R;
74 rr_proposed = Dot(deltaR, deltaR);
75 if (RandomGen() <= prob)
76 {
77 accepted = true;
78 rr_accepted = rr_proposed;
79 }
80 }
81 }
82
83 // recompute Psi if the move is rejected
84 if (!accepted)
85 {
86 W.loadWalker(thisWalker, false);
87 W.update();
88 logpsi = Psi.evaluateLog(W);
89 }
90
91 // evaluate Hamiltonian
92 enew = H.evaluateWithToperator(W);
93 H.auxHevaluate(W, thisWalker);
94 H.saveProperty(thisWalker.getPropertyBase());
95
96 // operate on thisWalker.
97 if (accepted)
98 {
99 W.saveWalker(thisWalker);
100 thisWalker.resetProperty(logpsi, Psi.getPhase(), enew, rr_accepted, rr_proposed, nodecorr);
101 }
102 else
103 {
104 thisWalker.Age++;
105 thisWalker.Properties(WP::R2ACCEPTED) = 0.0;
106 thisWalker.Properties(WP::R2PROPOSED) = rr_proposed;
107 }
108
109 const int NonLocalMoveAcceptedTemp = H.makeNonLocalMoves(W);
110 if (NonLocalMoveAcceptedTemp > 0)
111 {
112 W.saveWalker(thisWalker);
113 thisWalker.resetProperty(Psi.getLogPsi(), Psi.getPhase(), enew);
114 // debugging lines
115 //logpsi = Psi.getLogPsi();
116 //W.update(true);
117 //RealType logpsi2 = Psi.evaluateLog(W);
118 //if(logpsi!=logpsi2) std::cout << " logpsi " << logpsi << " logpsi2 " << logpsi2
119 // << " diff " << logpsi2-logpsi << std::endl;
120
121 NonLocalMoveAccepted += NonLocalMoveAcceptedTemp;
122 }
123
124 thisWalker.Weight *= branchEngine->branchWeight(enew, eold);
125 //branchEngine->accumulate(eold,1);
126 if (accepted)
127 ++nAccept;
128 else
129 ++nReject;
130
131 setMultiplicity(thisWalker);
132 }
133
134 /*
135 * DMCUpdateAllWithKill member functions
136 */
137 /// Constructor.
DMCUpdateAllWithKill(MCWalkerConfiguration & w,TrialWaveFunction & psi,QMCHamiltonian & h,RandomGenerator_t & rg)138 DMCUpdateAllWithKill::DMCUpdateAllWithKill(MCWalkerConfiguration& w,
139 TrialWaveFunction& psi,
140 QMCHamiltonian& h,
141 RandomGenerator_t& rg)
142 : QMCUpdateBase(w, psi, h, rg)
143 {
144 UpdatePbyP = false;
145 }
146
147 /// destructor
~DMCUpdateAllWithKill()148 DMCUpdateAllWithKill::~DMCUpdateAllWithKill() {}
149
150 /** advance all the walkers with killnode==yes
151 */
advanceWalker(Walker_t & thisWalker,bool recompute)152 void DMCUpdateAllWithKill::advanceWalker(Walker_t& thisWalker, bool recompute)
153 {
154 W.loadWalker(thisWalker, false);
155 //RealType nodecorr = setScaledDriftPbyPandNodeCorr(m_tauovermass,W.G,drift);
156 RealType nodecorr = setScaledDriftPbyPandNodeCorr(Tau, MassInvP, W.G, drift);
157 //create a 3N-Dimensional Gaussian with variance=1
158 makeGaussRandomWithEngine(deltaR, RandomGen);
159 //if(!W.makeMoveAllParticlesWithDrift(thisWalker,drift,deltaR, m_sqrttau))
160 if (!W.makeMoveAllParticlesWithDrift(thisWalker, drift, deltaR, SqrtTauOverMass))
161 {
162 H.rejectedMove(W, thisWalker);
163 return;
164 }
165 //save old local energy
166 RealType eold = thisWalker.Properties(WP::LOCALENERGY);
167 RealType enew = eold;
168 RealType logpsi(Psi.evaluateLog(W));
169 bool accepted = false;
170 RealType rr_accepted = 0.0;
171 nodecorr = 0.0;
172 if (branchEngine->phaseChanged(Psi.getPhaseDiff()))
173 {
174 thisWalker.Age++;
175 thisWalker.willDie();
176 }
177 else
178 {
179 enew = H.evaluate(W);
180 RealType logGf = -0.5 * Dot(deltaR, deltaR);
181 //nodecorr = setScaledDriftPbyPandNodeCorr(m_tauovermass,W.G,drift);
182 nodecorr = setScaledDriftPbyPandNodeCorr(Tau, MassInvP, W.G, drift);
183 deltaR = thisWalker.R - W.R - drift;
184 RealType logGb = logBackwardGF(deltaR);
185 //RealType logGb = -m_oneover2tau*Dot(deltaR,deltaR);
186 RealType prob = std::min(std::exp(logGb - logGf + 2.0 * (logpsi - thisWalker.Properties(WP::LOGPSI))), 1.0);
187 //calculate rr_proposed here
188 deltaR = W.R - thisWalker.R;
189 RealType rr_proposed = Dot(deltaR, deltaR);
190 if (RandomGen() > prob)
191 {
192 enew = eold;
193 thisWalker.Age++;
194 thisWalker.Properties(WP::R2ACCEPTED) = 0.0;
195 thisWalker.Properties(WP::R2PROPOSED) = rr_proposed;
196 H.rejectedMove(W, thisWalker);
197 }
198 else
199 {
200 thisWalker.Age = 0;
201 accepted = true;
202 W.saveWalker(thisWalker);
203 rr_accepted = rr_proposed;
204 thisWalker.resetProperty(logpsi, Psi.getPhase(), enew, rr_accepted, rr_proposed, nodecorr);
205 H.auxHevaluate(W, thisWalker);
206 H.saveProperty(thisWalker.getPropertyBase());
207 }
208 // std::cout <<logpsi<<" "<<Psi.getPhase()<<" "<<enew<<" "<<rr_accepted<<" "<<rr_proposed<<" "<<nodecorr<< std::endl;
209 thisWalker.Weight *= branchEngine->branchWeight(enew, eold);
210 }
211 if (accepted)
212 ++nAccept;
213 else
214 ++nReject;
215
216 setMultiplicity(thisWalker);
217 }
218 } // namespace qmcplusplus
219