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 // Raymond Clay III, j.k.rofling@gmail.com, Lawrence Livermore National Laboratory 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 /** @file 18 * @brief Declare QMCUpdateBase class 19 */ 20 #ifndef QMCPLUSPLUS_QMCUPDATE_BASE_H 21 #define QMCPLUSPLUS_QMCUPDATE_BASE_H 22 #include "Particle/MCWalkerConfiguration.h" 23 #include "QMCWaveFunctions/TrialWaveFunction.h" 24 #include "QMCHamiltonians/QMCHamiltonian.h" 25 #include "QMCHamiltonians/NonLocalTOperator.h" 26 #include "QMCDrivers/GreenFunctionModifiers/DriftModifierBase.h" 27 #include "QMCDrivers/SimpleFixedNodeBranch.h" 28 #include "Estimators/EstimatorManagerBase.h" 29 30 namespace qmcplusplus 31 { 32 class TraceManager; 33 /** @ingroup QMC 34 * @brief Base class for update methods for each step 35 * 36 * QMCUpdateBase provides the common functions to update all the walkers for each time step. 37 * Derived classes should implement advanceWalkers to complete a step. 38 */ 39 class QMCUpdateBase : public QMCTraits 40 { 41 public: 42 typedef MCWalkerConfiguration::Walker_t Walker_t; 43 typedef MCWalkerConfiguration::iterator WalkerIter_t; 44 typedef SimpleFixedNodeBranch BranchEngineType; 45 #ifdef MIXED_PRECISION 46 typedef TinyVector<OHMMS_PRECISION_FULL, DIM> mPosType; 47 typedef Tensor<OHMMS_PRECISION_FULL, DIM> mTensorType; 48 #else 49 typedef PosType mPosType; 50 typedef TensorType mTensorType; 51 #endif 52 53 ///If true, terminate the simulation, but it is never checked 54 bool BadState; 55 ///number of steps per measurement 56 int nSubSteps; 57 ///MaxAge>0 indicates branch is done 58 IndexType MaxAge; 59 ///counter for number of moves accepted 60 IndexType nAccept; 61 ///counter for number of moves rejected 62 IndexType nReject; 63 ///Total number of the steps when all the particle moves are rejected. 64 IndexType nAllRejected; 65 ///Total number of node crossings per block 66 IndexType nNodeCrossing; 67 ///Total numer of non-local moves accepted 68 IndexType NonLocalMoveAccepted; 69 ///timestep 70 RealType Tau; 71 ///spin mass 72 RealType spinMass; 73 ///use Drift 74 bool UseDrift; 75 76 /// Constructor. 77 QMCUpdateBase(MCWalkerConfiguration& w, TrialWaveFunction& psi, QMCHamiltonian& h, RandomGenerator_t& rg); 78 ///Alt Constructor. 79 QMCUpdateBase(MCWalkerConfiguration& w, 80 TrialWaveFunction& psi, 81 TrialWaveFunction& guide, 82 QMCHamiltonian& h, 83 RandomGenerator_t& rg); 84 ///destructor 85 virtual ~QMCUpdateBase(); 86 acceptRatio()87 inline RealType acceptRatio() const 88 { 89 return static_cast<RealType>(nAccept) / static_cast<RealType>(nAccept + nReject); 90 } 91 92 /** reset the QMCUpdateBase parameters 93 * @param brancher engine which handles branching 94 * 95 * Update time-step variables to move walkers 96 */ 97 void resetRun(BranchEngineType* brancher, 98 EstimatorManagerBase* est, 99 TraceManager* traces, 100 const DriftModifierBase* driftmodifer); 101 getTau()102 inline RealType getTau() 103 { 104 //SpeciesSet tspecies(W.getSpeciesSet()); 105 //int massind=tspecies.addAttribute("mass"); 106 //RealType mass = tspecies(massind,0); 107 //return m_tauovermass*mass; 108 return Tau; 109 } 110 setTau(RealType t)111 inline void setTau(RealType t) 112 { 113 //SpeciesSet tspecies(W.getSpeciesSet()); 114 //int massind=tspecies.addAttribute("mass"); 115 //RealType mass = tspecies(massind,0); 116 //RealType oneovermass = 1.0/mass; 117 //RealType oneoversqrtmass = std::sqrt(oneovermass); 118 // // Tau=brancher->getTau(); 119 // // assert (Tau==i); 120 //m_tauovermass = i/mass; 121 Tau = t; 122 m_tauovermass = t * MassInvS[0]; 123 m_oneover2tau = 0.5 / (m_tauovermass); 124 m_sqrttau = std::sqrt(m_tauovermass); 125 } 126 getSpinMass()127 inline RealType getSpinMass() { return spinMass; } 128 setSpinMass(RealType m)129 inline void setSpinMass(RealType m) { spinMass = m; } 130 getLogs(std::vector<RealType> & logs)131 inline void getLogs(std::vector<RealType>& logs) { Psi.getLogs(logs); } 132 set_step(int step)133 inline void set_step(int step) { W.current_step = step; } 134 135 136 ///** start a run */ 137 void startRun(int blocks, bool record); 138 /** stop a run */ 139 void stopRun(); 140 void stopRun2(); 141 142 /** prepare to start a block 143 * @param steps number of steps within the block 144 */ 145 void startBlock(int steps); 146 147 /** stop a block 148 */ 149 void stopBlock(bool collectall = true); 150 151 /** set the multiplicity of the walkers to branch */ 152 void setMultiplicity(WalkerIter_t it, WalkerIter_t it_end); 153 setMultiplicity(Walker_t & awalker)154 inline void setMultiplicity(Walker_t& awalker) const 155 { 156 constexpr RealType onehalf(0.5); 157 constexpr RealType cone(1); 158 RealType M = awalker.Weight; 159 if (awalker.Age > MaxAge) 160 M = std::min(onehalf, M); 161 else if (awalker.Age > 0) 162 M = std::min(cone, M); 163 awalker.Multiplicity = M + RandomGen(); 164 } 165 166 /** set the multiplicity of the walkers to branch */ 167 void setReleasedNodeMultiplicity(WalkerIter_t it, WalkerIter_t it_end); 168 169 /** initialize Walker buffers for PbyP update 170 */ 171 virtual void initWalkersForPbyP(WalkerIter_t it, WalkerIter_t it_end); 172 173 /** initialize Walker for walker update 174 */ 175 virtual void initWalkers(WalkerIter_t it, WalkerIter_t it_end); 176 177 /** process options 178 */ 179 virtual bool put(xmlNodePtr cur); 180 accumulate(WalkerIter_t it,WalkerIter_t it_end)181 inline void accumulate(WalkerIter_t it, WalkerIter_t it_end) { Estimators->accumulate(W, it, it_end); } 182 183 /** advance walkers executed at each step 184 * 185 * Derived classes implement how to move walkers and accept/reject 186 * moves. 187 */ 188 virtual void advanceWalkers(WalkerIter_t it, WalkerIter_t it_end, bool recompute); 189 190 ///move a walker 191 virtual void advanceWalker(Walker_t& thisWalker, bool recompute) = 0; 192 advanceWalkerForEE(Walker_t & w1,std::vector<PosType> & dR,std::vector<int> & iats,std::vector<int> & rs,std::vector<RealType> & ratios)193 virtual RealType advanceWalkerForEE(Walker_t& w1, 194 std::vector<PosType>& dR, 195 std::vector<int>& iats, 196 std::vector<int>& rs, 197 std::vector<RealType>& ratios) 198 { 199 return 0.0; 200 }; 201 202 ///normalization offset for cs type runs. 203 RealType csoffset; 204 205 // virtual void estimateNormWalkers(std::vector<TrialWaveFunction*>& pclone 206 // , std::vector<MCWalkerConfiguration*>& wclone 207 // , std::vector<QMCHamiltonian*>& hclone 208 // , std::vector<RandomGenerator_t*>& rng 209 // , std::vector<RealType>& ratio_i_0){}; RMC_checkIndex(int N,int NMax)210 int RMC_checkIndex(int N, int NMax) 211 { 212 if (N < 0) 213 return N + NMax; 214 else if (N >= NMax) 215 return N - NMax; 216 else 217 return N; 218 } 219 RMC_checkWalkerBounds(WalkerIter_t & it,WalkerIter_t first,WalkerIter_t last)220 void RMC_checkWalkerBounds(WalkerIter_t& it, WalkerIter_t first, WalkerIter_t last) 221 { 222 if (it >= last) 223 it -= (last - first); 224 else if (it < first) 225 it += (last - first); 226 } 227 logBackwardGF(const ParticleSet::ParticlePos_t & displ)228 inline RealType logBackwardGF(const ParticleSet::ParticlePos_t& displ) 229 { 230 RealType logGb = 0.0; 231 for (int iat = 0; iat < W.getTotalNum(); ++iat) 232 { 233 RealType mass_over_tau = 1.0 / (SqrtTauOverMass[iat] * SqrtTauOverMass[iat]); 234 logGb += 0.5 * dot(displ[iat], displ[iat]) * mass_over_tau; 235 } 236 return -logGb; 237 } 238 239 public: 240 ///traces 241 TraceManager* Traces; 242 243 protected: 244 ///update particle-by-particle 245 bool UpdatePbyP; 246 ///number of particles 247 IndexType NumPtcl; 248 ///Time-step factor \f$ 1/(2\tau)\f$ 249 RealType m_oneover2tau; 250 ///Time-step factor \f$ \sqrt{\tau}\f$ 251 RealType m_sqrttau; 252 ///tau/mass 253 RealType m_tauovermass; 254 ///maximum displacement^2 255 RealType m_r2max; 256 ///walker ensemble 257 MCWalkerConfiguration& W; 258 ///trial function 259 TrialWaveFunction& Psi; 260 ///guide function 261 TrialWaveFunction& Guide; 262 ///Hamiltonian 263 QMCHamiltonian& H; 264 ///random number generator 265 RandomGenerator_t& RandomGen; 266 ///branch engine, stateless reference to the one in QMCDriver 267 const BranchEngineType* branchEngine; 268 ///drift modifer, stateless reference to the one in QMCDriver 269 const DriftModifierBase* DriftModifier; 270 ///estimator 271 EstimatorManagerBase* Estimators; 272 ///parameters 273 ParameterSet myParams; 274 ///1/Mass per species 275 std::vector<RealType> MassInvS; 276 ///1/Mass per particle 277 std::vector<RealType> MassInvP; 278 ///sqrt(tau/Mass) per particle 279 std::vector<RealType> SqrtTauOverMass; 280 ///temporary storage for drift 281 ParticleSet::ParticlePos_t drift; 282 ///temporary storage for random displacement 283 ParticleSet::ParticlePos_t deltaR; 284 ///temporart storage for spin displacement 285 ParticleSet::ParticleScalar_t deltaS; 286 ///storage for differential gradients for PbyP update 287 ParticleSet::ParticleGradient_t G, dG; 288 ///storage for differential laplacians for PbyP update 289 ParticleSet::ParticleLaplacian_t L, dL; 290 291 /** evaluate the ratio of scaled velocity and velocity 292 * @param g gradient 293 * @param gscaled scaled gradient 294 * @return the ratio 295 */ 296 RealType getNodeCorrection(const ParticleSet::ParticleGradient_t& g, ParticleSet::ParticlePos_t& gscaled); 297 298 ///copy constructor (disabled) 299 QMCUpdateBase(const QMCUpdateBase&) = delete; 300 301 /// check logpsi and grad and lap against values computed from scratch 302 static bool checkLogAndGL(ParticleSet& pset, TrialWaveFunction& twf); 303 304 private: 305 ///set default parameters 306 void setDefaults(); 307 /// Copy operator (disabled). 308 QMCUpdateBase& operator=(const QMCUpdateBase&) { return *this; } 309 /// 310 NewTimer* InitWalkersTimer; 311 }; 312 } // namespace qmcplusplus 313 314 #endif 315