From bb129bff4e2c2452d9553ef94559d18f2399e75f Mon Sep 17 00:00:00 2001 From: Jeongnim Kim Date: Tue, 7 Dec 2004 23:34:18 +0000 Subject: [PATCH] ratio/update for Green functions. git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@141 e5b18d87-469d-4833-9cc0-8cdfa06e9491 --- src/QMC/VMCParticleByParticle.cpp | 43 +++++++++--- src/QMCWaveFunctions/DiracDeterminant.h | 70 +++++++++++++++++-- src/QMCWaveFunctions/MultiSlaterDeterminant.h | 29 ++++++-- src/QMCWaveFunctions/OneBodyJastrowFunction.h | 19 +++++ src/QMCWaveFunctions/OrbitalBase.h | 18 +++++ src/QMCWaveFunctions/PolarizedJastrow.h | 15 ++++ src/QMCWaveFunctions/SlaterDeterminant.h | 20 ++++++ src/QMCWaveFunctions/TrialWaveFunction.cpp | 23 +++++- src/QMCWaveFunctions/TrialWaveFunction.h | 9 ++- .../TwoBodyJastrowFunction.Shared.h | 30 ++++++++ src/QMCWaveFunctions/TwoBodyJastrowFunction.h | 32 ++++++++- 11 files changed, 283 insertions(+), 25 deletions(-) diff --git a/src/QMC/VMCParticleByParticle.cpp b/src/QMC/VMCParticleByParticle.cpp index 861fd8cb7..60aac1845 100644 --- a/src/QMC/VMCParticleByParticle.cpp +++ b/src/QMC/VMCParticleByParticle.cpp @@ -50,7 +50,6 @@ namespace ohmmsqmc { //set the data members to start a new run getReady(); - //going to add routines to calculate how much we need bool require_register = W.createAuxDataSet(); @@ -80,8 +79,8 @@ namespace ohmmsqmc { MCWalkerConfiguration::PropertyContainer_t Properties; ParticleSet::ParticlePos_t deltaR(W.getTotalNum()); - ParticleSet::ParticleGradient_t G(W.getTotalNum()); - ParticleSet::ParticleLaplacian_t L(W.getTotalNum()); + ParticleSet::ParticleGradient_t G(W.getTotalNum()), dG(W.getTotalNum()); + ParticleSet::ParticleLaplacian_t L(W.getTotalNum()), dL(W.getTotalNum()); IndexType accstep=0; IndexType nAcceptTot = 0; @@ -108,19 +107,48 @@ namespace ohmmsqmc { //create a 3N-Dimensional Gaussian with variance=1 makeGaussRandom(deltaR); bool moved = false; + for(int iat=0; iatDrift[iat]; PosType dr = g*deltaR[iat]; W.makeMove(iat,dr); - RealType ratio = Psi.ratio(W,iat); + + //RealType ratio = Psi.ratio(W,iat); + RealType ratio = Psi.ratio(W,iat,dG,dL); + + G = W.G+dG; + + /* green function + deltaR = -scale*G; + dletaR[iat] -= dr; //subtract dr = (*it)->R(iat) - W.R(iat) + RealType backwardGF = exp(-oneover2tau*Dot(deltaR,deltaR)); + //a better way to do this + RealType forwardGF = exp(-0.5*dot(dr,dr)); + ValueType vsq = Dot(G,G); + ValueType scale = ((-1.0+sqrt(1.0+2.0*Tau*vsq))/vsq); + ValueType scale2 = scale*scale; + RealType dist = scale2*vsq+2.0*scale*dot(dr,G(iat))+dot(dr,dr); + RealType backwardGF = exp(-oneover2tau*dist); + */ + RealType ratio2 = pow(ratio,2); + + //alternatively + //if(Random() > backwardGF/forwardGF*ratio2) { if(Random() < ratio2) { moved = true; ++nAccept; W.acceptMove(iat); - Psi.update(W,iat); + //Psi.update(W,iat); + Psi.update2(W,iat); + W.G = G; + W.L += dL; + //Need to change the drift + //(*it)->Drift = scale*G; } else { ++nReject; + Psi.restore(iat); } } @@ -137,10 +165,7 @@ namespace ohmmsqmc { (*it)->Properties(PSI) = psi; (*it)->Properties(LOCALENERGY) = H.evaluate(W); H.copy((*it)->getEnergyBase()); - //GreenFunction is not used here - //ValueType vsq = Dot(W.G,W.G); - //ValueType scale = ((-1.0+sqrt(1.0+2.0*Tau*vsq))/vsq); - //(*it)->Drift = scale*W.G; + (*it)->Properties(LOCALPOTENTIAL) = H.getLocalPotential(); } //Keep until everything is tested: debug routine // if(moved) { diff --git a/src/QMCWaveFunctions/DiracDeterminant.h b/src/QMCWaveFunctions/DiracDeterminant.h index 1f012ced4..d054df4c4 100644 --- a/src/QMCWaveFunctions/DiracDeterminant.h +++ b/src/QMCWaveFunctions/DiracDeterminant.h @@ -160,6 +160,9 @@ namespace ohmmsqmc { psiM.resize(nel,norb); dpsiM.resize(nel,norb); d2psiM.resize(nel,norb); + psiM_temp.resize(nel,norb); + dpsiM_temp.resize(nel,norb); + d2psiM_temp.resize(nel,norb); psiMinv.resize(nel,norb); LastIndex = FirstIndex + nel; } @@ -176,6 +179,8 @@ namespace ohmmsqmc { NP=P.getTotalNum(); myG.resize(NP); myL.resize(NP); + myG_temp.resize(NP); + myL_temp.resize(NP); FirstAddressOfG = &myG[0][0]; LastAddressOfG = FirstAddressOfG + NP*DIM; FirstAddressOfdV = &(dpsiM(0,0)[0]); //(*dpsiM.begin())[0]); @@ -209,6 +214,61 @@ namespace ohmmsqmc { buf.get(myL.begin(), myL.end()); buf.get(FirstAddressOfG,LastAddressOfG); buf.get(CurrentDet); + + //need extra copy for gradient/laplacian calculations without updating it + psiM_temp = psiM; + dpsiM_temp = dpsiM; + d2psiM_temp = d2psiM; + } + + ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL) { + + Phi.evaluate(P, iat, psiV, dpsiV, d2psiV); + WorkingIndex = iat-FirstIndex; + + curRatio= DetRatio(psiM_temp, psiV.begin(),WorkingIndex); + + DetUpdate(psiM_temp,psiV,workV1,workV2,WorkingIndex,curRatio); + + for(int j=0; j psiV; @@ -326,8 +386,8 @@ namespace ohmmsqmc { ValueType *FirstAddressOfdV; ValueType *LastAddressOfdV; - ParticleSet::ParticleGradient_t myG; - ParticleSet::ParticleLaplacian_t myL; + ParticleSet::ParticleGradient_t myG, myG_temp; + ParticleSet::ParticleLaplacian_t myL, myL_temp; }; /** Calculate the value of the Dirac determinant for particles diff --git a/src/QMCWaveFunctions/MultiSlaterDeterminant.h b/src/QMCWaveFunctions/MultiSlaterDeterminant.h index a81683c25..39a17b8be 100644 --- a/src/QMCWaveFunctions/MultiSlaterDeterminant.h +++ b/src/QMCWaveFunctions/MultiSlaterDeterminant.h @@ -72,13 +72,6 @@ public: void initParameters() { } - inline ValueType - ratio(ParticleSet& P, int iat) { - std::cerr << "MultiSlaterDeterminant should not be used by Particle-By-Particle update" - << std::endl; - return 1.0; - } - inline ValueType evaluate(ParticleSet& P, //const DistanceTableData* dtable, ParticleSet::ParticleGradient_t& G, @@ -123,6 +116,28 @@ public: std::cerr << "MultiSlaterDeterminant::putData is empty" << std::endl; } + inline ValueType + ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& G, + ParticleSet::ParticleLaplacian_t& L) { + std::cerr << "MultiSlaterDeterminant should not be used by Particle-By-Particle update" + << std::endl; + return 1.0; + } + + inline void restore(int iat) { } + + void update(ParticleSet& P, int iat) { + std::cerr << "MultiSlaterDeterminant::update is empty" << std::endl; + } + + inline ValueType + ratio(ParticleSet& P, int iat) { + std::cerr << "MultiSlaterDeterminant should not be used by Particle-By-Particle update" + << std::endl; + return 1.0; + } + void update(ParticleSet& P, ParticleSet::ParticleGradient_t& G, ParticleSet::ParticleLaplacian_t& L, diff --git a/src/QMCWaveFunctions/OneBodyJastrowFunction.h b/src/QMCWaveFunctions/OneBodyJastrowFunction.h index fb9d0ef83..8fa8f860e 100644 --- a/src/QMCWaveFunctions/OneBodyJastrowFunction.h +++ b/src/QMCWaveFunctions/OneBodyJastrowFunction.h @@ -131,6 +131,25 @@ namespace ohmmsqmc { return exp(-sumu); } + + ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL) { + ValueType v=ratio(P,iat); + dG[iat] += curGrad-dU[iat]; + dL[iat] += curLap-d2U[iat]; + return v; + } + + inline void restore(int iat) {} + + void update(ParticleSet& P, int iat) { + U[iat] = curVal; + dU[iat]=curGrad; + d2U[iat]=curLap; + } + + ValueType ratio(ParticleSet& P, int iat) { int n=d_table->size(VisitorIndex); curVal=0.0; diff --git a/src/QMCWaveFunctions/OrbitalBase.h b/src/QMCWaveFunctions/OrbitalBase.h index f454f2780..effb6b9c2 100644 --- a/src/QMCWaveFunctions/OrbitalBase.h +++ b/src/QMCWaveFunctions/OrbitalBase.h @@ -77,6 +77,23 @@ namespace ohmmsqmc { WalkerSetRef::WalkerLaplacian_t& L) = 0; + /** evaluate the ratio of the new to old orbital value + *@param P the active ParticleSet + *@param iat the index of a particle + *@param dG the differential gradient + *@param dL the differential laplacian + *@return \f$\phi(\{\bf R}^{'}\})/\phi(\{\bf R}^{'}\})\f$. + * + *Paired with update(ParticleSet& P, int iat). + */ + virtual ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL) = 0; + + virtual void update(ParticleSet& P, int iat) =0; + + virtual void restore(int iat) = 0; + /** evalaute the ratio of the new to old orbital value *@param P the active ParticleSet *@param iat the index of a particle @@ -101,6 +118,7 @@ namespace ohmmsqmc { ParticleSet::ParticleLaplacian_t& dL, int iat) =0; + /** equivalent to evaluate(P,G,L) with write-back function */ virtual ValueType evaluate(ParticleSet& P,PooledData& buf)=0; diff --git a/src/QMCWaveFunctions/PolarizedJastrow.h b/src/QMCWaveFunctions/PolarizedJastrow.h index 060f43ff6..029e1d9c5 100644 --- a/src/QMCWaveFunctions/PolarizedJastrow.h +++ b/src/QMCWaveFunctions/PolarizedJastrow.h @@ -77,6 +77,20 @@ namespace ohmmsqmc { } #endif + ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& G, + ParticleSet::ParticleLaplacian_t& L) { + std::cerr << "PolarizedJastrow::ratio for particle-by-particle is empty " << std::endl; + return 1.0; + } + + inline void restore(int iat) { } + + //@todo implement the virutal functions for particle-by-particle move + void update(ParticleSet& P, int iat) { + std::cerr << "PolarizedJastrow::update for particle-by-particle is empty " << std::endl; + } + ValueType ratio(ParticleSet& P, int iat) { std::cerr << "PolarizedJastrow::ratio for particle-by-particle is empty " << std::endl; return 1.0; @@ -90,6 +104,7 @@ namespace ohmmsqmc { std::cerr << "PolarizedJastrow::update for particle-by-particle is empty " << std::endl; } + ValueType evaluate(ParticleSet& P,PooledData& buf) { std::cerr << "PolarizedJastrow::evaluate for particle-by-particle is empty " << std::endl; return 1.0; diff --git a/src/QMCWaveFunctions/SlaterDeterminant.h b/src/QMCWaveFunctions/SlaterDeterminant.h index ebc31ffac..8ce9b7bc1 100644 --- a/src/QMCWaveFunctions/SlaterDeterminant.h +++ b/src/QMCWaveFunctions/SlaterDeterminant.h @@ -119,6 +119,26 @@ namespace ohmmsqmc { return r; } + inline ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL) { + int i=1; + while(iat>=M[i]) {i++;} + return Dets[i-1]->ratio(P,iat, dG, dL); + } + + inline void restore(int iat) { + int i=1; + while(iat>=M[i]) {i++;} + Dets[i-1]->restore(iat); + } + + inline void update(ParticleSet& P, int iat) { + int i=1; + while(iat>=M[i]) {i++;} + Dets[i-1]->update(P,iat); + } + ValueType ratio(ParticleSet& P, int iat) { int i=1; diff --git a/src/QMCWaveFunctions/TrialWaveFunction.cpp b/src/QMCWaveFunctions/TrialWaveFunction.cpp index fc5197b4c..ff800ace8 100644 --- a/src/QMCWaveFunctions/TrialWaveFunction.cpp +++ b/src/QMCWaveFunctions/TrialWaveFunction.cpp @@ -82,7 +82,6 @@ namespace ohmmsqmc { for(int i=0; iratio(P,iat); return r; } - void TrialWaveFunction::update(ParticleSet& P,int iat) { //ready to collect "changes" in the gradients and laplacians by the move @@ -92,6 +91,28 @@ namespace ohmmsqmc { P.L += delta_L; } + + TrialWaveFunction::ValueType + TrialWaveFunction::ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL) { + dG = 0.0; + dL = 0.0; + RealType r=1.0; + for(int i=0; iratio(P,iat,dG,dL); + return r; + } + + void + TrialWaveFunction::restore(int iat) { + for(int i=0; irestore(iat); + } + + void + TrialWaveFunction::update2(ParticleSet& P,int iat) { + for(int i=0; iupdate(P,iat); + } + void TrialWaveFunction::resizeByWalkers(int nwalkers){ for(int i=0; iresizeByWalkers(nwalkers); } diff --git a/src/QMCWaveFunctions/TrialWaveFunction.h b/src/QMCWaveFunctions/TrialWaveFunction.h index 920dd9c63..64298d59f 100644 --- a/src/QMCWaveFunctions/TrialWaveFunction.h +++ b/src/QMCWaveFunctions/TrialWaveFunction.h @@ -68,9 +68,16 @@ namespace ohmmsqmc { /** functions to handle particle-by-particle update */ ValueType ratio(ParticleSet& P, int iat); + void update(ParticleSet& P, int iat); + + ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL); + void restore(int iat); + void update2(ParticleSet& P, int iat); + void registerData(ParticleSet& P, PooledData& buf); void copyFromBuffer(ParticleSet& P, PooledData& buf); - void update(ParticleSet& P, int iat); ValueType evaluate(ParticleSet& P, PooledData& buf); /** evalaute the values of the wavefunction, gradient and laplacian for all the walkers */ diff --git a/src/QMCWaveFunctions/TwoBodyJastrowFunction.Shared.h b/src/QMCWaveFunctions/TwoBodyJastrowFunction.Shared.h index 21a4e074a..6ee587a54 100644 --- a/src/QMCWaveFunctions/TwoBodyJastrowFunction.Shared.h +++ b/src/QMCWaveFunctions/TwoBodyJastrowFunction.Shared.h @@ -89,6 +89,36 @@ namespace ohmmsqmc { return exp(-sumu); } + /** later merge the loop */ + ValueType ratio(ParticleSet& P, int iat, + ParticleSet::ParticleGradient_t& dG, + ParticleSet::ParticleLaplacian_t& dL) { + ValueType v=ratio(P,iat); + GradType sumg,dg; + ValueType suml=0.0,dl; + for(int jat=0,ij=iat*N,ji=iat; jat