scnc
/
qmcpack
mirror of https://github.com/QMCPACK/qmcpack.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

ratio/update for Green functions. 

git-svn-id: https://subversion.assembla.com/svn/qmcdev/trunk@141 e5b18d87-469d-4833-9cc0-8cdfa06e9491
This commit is contained in:
Jeongnim Kim 2004-12-07 23:34:18 +00:00
parent 76cacd4c34
commit bb129bff4e
11 changed files with 283 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
src/QMC/VMCParticleByParticle.cpp
View File

@ -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; iat<W.getTotalNum(); iat++) {
//PosType dr = g*deltaR[iat]+(*it)->Drift[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) {

70
src/QMCWaveFunctions/DiracDeterminant.h
View File

@ -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<cols(); j++) {
dpsiM_temp(WorkingIndex,j)=dpsiV[j];
d2psiM_temp(WorkingIndex,j)=d2psiV[j];
}
int kat=FirstIndex;
for(int i=0; i<rows(); i++,kat++) {
PosType rv =psiM_temp(i,0)*dpsiM_temp(i,0);
ValueType lap=psiM_temp(i,0)*d2psiM_temp(i,0);
for(int j=1; j<cols(); j++) {
rv += psiM_temp(i,j)*dpsiM_temp(i,j);
lap += psiM_temp(i,j)*d2psiM_temp(i,j);
}
lap -= dot(rv,rv);
dG[kat] += rv - myG[kat]; myG_temp[kat]=rv;
dL[kat] += lap -myL[kat]; myL_temp[kat]=lap;
}
return curRatio;
}
void update(ParticleSet& P, int iat) {
CurrentDet *= curRatio;
myG = myG_temp;
myL = myL_temp;
psiM = psiM_temp;
for(int j=0; j<cols(); j++) {
dpsiM(WorkingIndex,j)=dpsiV[j];
d2psiM(WorkingIndex,j)=d2psiV[j];
}
}
void restore(int iat) {
psiM_temp = psiM;
for(int j=0; j<cols(); j++) {
dpsiM_temp(WorkingIndex,j)=dpsiM(WorkingIndex,j);
d2psiM_temp(WorkingIndex,j)=d2psiM(WorkingIndex,j);
}
}
ValueType ratio(ParticleSet& P, int iat) {
@ -298,16 +358,16 @@ namespace ohmmsqmc {
ValueType CurrentDet;
/// psiM(j,i) \f$= \psi_j({\bf r}_i)\f$
Determinant_t psiM;
Determinant_t psiM, psiM_temp;
/// temporary container for testing
Determinant_t psiMinv;
/// dpsiM(i,j) \f$= \nabla_i \psi_j({\bf r}_i)\f$
Gradient_t dpsiM;
Gradient_t dpsiM, dpsiM_temp;
/// d2psiM(i,j) \f$= \nabla_i^2 \psi_j({\bf r}_i)\f$
Laplacian_t d2psiM;
Laplacian_t d2psiM, d2psiM_temp;
/// value of single-particle orbital for particle-by-particle update
std::vector<ValueType> 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

29
src/QMCWaveFunctions/MultiSlaterDeterminant.h
View File

@ -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,

19
src/QMCWaveFunctions/OneBodyJastrowFunction.h
View File

@ -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;

18
src/QMCWaveFunctions/OrbitalBase.h
View File

@ -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<RealType>& buf)=0;

15
src/QMCWaveFunctions/PolarizedJastrow.h
View File

@ -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<RealType>& buf) {
std::cerr << "PolarizedJastrow::evaluate for particle-by-particle is empty " << std::endl;
return 1.0;

20
src/QMCWaveFunctions/SlaterDeterminant.h
View File

@ -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;

23
src/QMCWaveFunctions/TrialWaveFunction.cpp
View File

@ -82,7 +82,6 @@ namespace ohmmsqmc {
for(int i=0; i<Z.size(); i++) r *= Z[i]->ratio(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; i<Z.size(); i++) r *= Z[i]->ratio(P,iat,dG,dL);
return r;
}
void
TrialWaveFunction::restore(int iat) {
for(int i=0; i<Z.size(); i++) Z[i]->restore(iat);
}
void
TrialWaveFunction::update2(ParticleSet& P,int iat) {
for(int i=0; i<Z.size(); i++) Z[i]->update(P,iat);
}
void TrialWaveFunction::resizeByWalkers(int nwalkers){
for(int i=0; i<Z.size(); i++) Z[i]->resizeByWalkers(nwalkers);
}

9
src/QMCWaveFunctions/TrialWaveFunction.h
View File

@ -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<RealType>& buf);
void copyFromBuffer(ParticleSet& P, PooledData<RealType>& buf);
void update(ParticleSet& P, int iat);
ValueType evaluate(ParticleSet& P, PooledData<RealType>& buf);
/** evalaute the values of the wavefunction, gradient and laplacian for all the walkers */

30
src/QMCWaveFunctions/TwoBodyJastrowFunction.Shared.h
View File

@ -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<N; jat++,ij++,ji+=N) {
sumg += (dg=curGrad[jat]-dU[ij]);
suml += (dl=curLap[jat]-d2U[ij]);
dG[jat] -= dg;
dL[jat] += dl;
}
dG[iat] += sumg;
dL[iat] += suml;
return v;
}
inline void restore(int iat) {}
void update(ParticleSet& P, int iat) {
DiffValSum += DiffVal;
for(int jat=0,ij=iat*N,ji=iat; jat<N; jat++,ij++,ji+=N) {
dU[ij]=curGrad[jat];
dU[ji]=-1.0*curGrad[jat];
d2U[ij]=d2U[ji] = curLap[jat];
U[ij] = U[ji] = curVal[jat];
}
}
/** evalaute ratio
*@param P the active particle set
*@param iat the index of the particle that is moved

32
src/QMCWaveFunctions/TwoBodyJastrowFunction.h
View File

@ -173,6 +173,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<N; jat++,ij++,ji+=N) {
sumg += (dg=curGrad[jat]-dU[ij]);
suml += (dl=curLap[jat]-d2U[ij]);
dG[jat] -= dg;
dL[jat] += dl;
}
dG[iat] += sumg;
dL[iat] += suml;
return v;
}
inline void restore(int iat) {}
void update(ParticleSet& P, int iat) {
DiffValSum += DiffVal;
for(int jat=0,ij=iat*N,ji=iat; jat<N; jat++,ij++,ji+=N) {
dU[ij]=curGrad[jat];
dU[ji]=-1.0*curGrad[jat];
d2U[ij]=d2U[ji] = curLap[jat];
U[ij] = U[ji] = curVal[jat];
}
}
ValueType ratio(ParticleSet& P, int iat) {
register ValueType dudr, d2udr2,u;
register PosType gr;
@ -190,8 +220,6 @@ namespace ohmmsqmc {
}
return exp(DiffVal);
}
inline void update(ParticleSet& P,
ParticleSet::ParticleGradient_t& G,
ParticleSet::ParticleLaplacian_t& L,