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