diff --git a/src/QMC/DMCParticleByParticle.cpp b/src/QMC/DMCParticleByParticle.cpp index 2653b6399..93819ba1e 100644 --- a/src/QMC/DMCParticleByParticle.cpp +++ b/src/QMC/DMCParticleByParticle.cpp @@ -60,16 +60,16 @@ namespace ohmmsqmc { brancher.put(qmc_node,LogOut); if(BranchInfo != "default") brancher.read(BranchInfo); - //else { - // /*if VMC/DMC directly preceded DMC (Counter > 0) then - // use the average value of the energy estimator for - // the reference energy of the brancher*/ - // if(Counter) { - // RealType e_ref = W.getLocalEnergy(); - // LOGMSG("Overwriting the reference energy by the local energy " << e_ref) - // brancher.setEguess(e_ref); - // } - //} + else { + /*if VMC/DMC directly preceded DMC (Counter > 0) then + use the average value of the energy estimator for + the reference energy of the brancher*/ + if(Counter) { + RealType e_ref = W.getLocalEnergy(); + LOGMSG("Overwriting the reference energy by the local energy " << e_ref) + brancher.setEguess(e_ref); + } + } int PopIndex, E_TIndex; Estimators.resetReportSettings(RootName); @@ -85,8 +85,6 @@ namespace ohmmsqmc { (*it)->Properties(MULTIPLICITY) = 1.0; ++it; } - - //going to add routines to calculate how much we need bool require_register = W.createAuxDataSet(); @@ -110,6 +108,7 @@ namespace ohmmsqmc { int Population = W.getActiveWalkers(); int tPopulation = W.getActiveWalkers(); RealType Eest = brancher.E_T; + RealType E_T = Eest; RealType oneovertau = 1.0/Tau; RealType oneover2tau = 0.5*oneovertau; RealType g = sqrt(Tau); @@ -239,7 +238,8 @@ namespace ohmmsqmc { ++step;++accstep; Estimators.accumulate(W); - Eest = brancher.update(Population,Eest); //E_T = brancher.update(Population,Eest); + Eest = brancher.update(Population,Eest); + //E_T = brancher.update(Population,Eest); brancher.branch(accstep,W); } while(step(mlimit)-Feed*log(static_cast(pop_now))+logN; + //cout << "update trial energy = " << Eg.front()<< " " << Eg.back() << " " << E_T << endl; return E_T; } @@ -184,16 +185,17 @@ namespace ohmmsqmc { XMLReport("MaxCopy for branching = " << MaxCopy) XMLReport("Branching: Number of generations = " << InFeed) XMLReport("Branching: Feedback parameter = " << Feed) + //for(int i=0; i(EgBufferSize); - //vector etrial(Eg.begin(),Eg.end()); + vector etrial(Eg.begin(),Eg.end()); hid_t dataspace = H5Screate_simple(1, &dim, NULL); hid_t dataset = H5Dcreate(grp, "TrialEnergies", H5T_NATIVE_DOUBLE, dataspace, H5P_DEFAULT); hid_t ret = - H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,&Eg[0]); + H5Dwrite(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT,&etrial[0]); H5Dclose(dataset); H5Sclose(dataspace); //make a char array which contains the species names separated by space @@ -213,11 +215,15 @@ namespace ohmmsqmc { h5file.append(".config.h5"); hid_t h_file = H5Fopen(h5file.c_str(),H5F_ACC_RDWR,H5P_DEFAULT); hid_t h_config = H5Gopen(h_file,"config_collection"); - hsize_t dataset = H5Dopen(h_config, "TrialEnergies"); - if(dataset) { - hsize_t ret = H5Dread(dataset, H5T_NATIVE_DOUBLE, - H5S_ALL, H5S_ALL, H5P_DEFAULT, &Eg[0]); + hid_t dataset = H5Dopen(h_config, "TrialEnergies"); + if(dataset<0) { + Counter=0; + reset(); + } else { + vector etrial(EgBufferSize); + hsize_t ret = H5Dread(dataset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &etrial[0]); H5Dclose(dataset); + std::copy(etrial.begin(),etrial.end(),Eg.begin()); char temp[256]; dataset = H5Dopen(h_config, "BranchParameters"); ret = H5Dread(dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL, H5P_DEFAULT, temp); @@ -226,10 +232,9 @@ namespace ohmmsqmc { s >> Nideal>>InFeed >>Counter >> E_T; reset(); Counter = std::min(Counter,EgBufferSize); - } + } H5Gclose(h_config); H5Fclose(h_file); - } private: