C++ MCWalkerConfiguration::getActiveWalkers方法代码示例

 //NOTE: weights are not handled nicely.
 //weight should be done carefully, not valid for DMC
 //will add a function to MCWalkerConfiguration to track total weight
 void EstimatorManager::accumulate(MCWalkerConfiguration& W)
 {
   BlockWeight += W.getActiveWalkers();
   RealType norm=1.0/W.getGlobalNumWalkers();
   for(int i=0; i< Estimators.size(); i++) 
     Estimators[i]->accumulate(W,W.begin(),W.end(),norm);
 }

/** accumulate Local energies and collectables
 * @param W ensemble
 */
void EstimatorManager::accumulate(MCWalkerConfiguration& W)
{
  BlockWeight += W.getActiveWalkers();
  RealType norm=1.0/W.getGlobalNumWalkers();
  for(int i=0; i< Estimators.size(); i++)
    Estimators[i]->accumulate(W,W.begin(),W.end(),norm);
  if(Collectables)//collectables are normalized by QMC drivers
    Collectables->accumulate_all(W.Collectables,1.0);
}

/**  Write the set of walker configurations to the HDF5 file.  
 *@param W set of walker configurations
 *@param ic the number of frames
 *
 * \if ic==-1
 *  use only the last frame for a restart
 * \else if ic>=0
 *  use ic frames from the file for opitimizations
 */
bool  
HDFWalkerInput0::put(MCWalkerConfiguration& W, int ic){

  if(Counter<0) return false;

  int selected = ic;
  if(ic<0) {
    XMLReport("Will use the last set from " << NumSets << " of configurations.")
    selected = NumSets-1;
  }

  typedef MCWalkerConfiguration::PosType PosType;
  typedef MCWalkerConfiguration::PropertyContainer_t ProtertyContainer_t;

  typedef Matrix<PosType>  PosContainer_t;

  int nwt = 0;
  int npt = 0;
  //2D array of PosTypes (x,y,z) indexed by (walker,particle)
  PosContainer_t Pos_temp;

  //open the group
  char GrpName[128];
  sprintf(GrpName,"config%04d",selected);
  hid_t group_id = H5Gopen(h_config,GrpName);
    
  HDFAttribIO<PosContainer_t> Pos_in(Pos_temp);
  //read the dataset
  Pos_in.read(group_id,"coord");
  //close the group
  H5Gclose(group_id);

  /*check to see if the number of walkers and particles is  consistent with W */
  int nptcl = Pos_temp.cols();
  nwt = Pos_temp.rows();

  int curWalker = W.getActiveWalkers();
  if(curWalker) {
    LOGMSG("Adding " << nwt << " walkers to " << curWalker)
    W.createWalkers(nwt);
  } else {
    W.resize(nwt,nptcl); 
  }

  //assign configurations to W
  int iw=0;
  MCWalkerConfiguration::iterator it = W.begin()+curWalker; 
  MCWalkerConfiguration::iterator it_end = W.end(); 
  while(it != it_end) {
    std::copy(Pos_temp[iw],Pos_temp[iw+1], (*it)->R.begin());
    ++it;++iw;
  }

  return true;
}

bool  
HDFWalkerInput0::append(MCWalkerConfiguration& W, int blocks){

  if(Counter<0) return false;

  //if(nwalkers<0) return put(W,-1);

  typedef MCWalkerConfiguration::PosType PosType;
  typedef Matrix<PosType>  PosContainer_t;
  PosContainer_t Pos_temp;

  int nw_in=0;
  int firstConf=std::max(0,NumSets-blocks);
  if(blocks<0) firstConf=0;

  for(int iconf=firstConf; iconf<NumSets; iconf++) {
    //open the group
    char GrpName[128];
    sprintf(GrpName,"config%04d",iconf);
    hid_t group_id = H5Gopen(h_config,GrpName);
    HDFAttribIO<PosContainer_t> Pos_in(Pos_temp);
    //read the dataset
    Pos_in.read(group_id,"coord");
    //close the group
    H5Gclose(group_id);
    /*check to see if the number of walkers and particles is  consistent with W */
    int nptcl = Pos_temp.cols();
    int nwt = Pos_temp.rows();
    int curWalker=0;
    if(nptcl != W.getParticleNum()) {
      W.resize(nwt,nptcl); 
    } else {
      curWalker=W.getActiveWalkers();
      W.createWalkers(nwt);
    }
    MCWalkerConfiguration::iterator it = W.begin()+curWalker; 
    for(int iw=0; iw<nwt; iw++) {
      //std::copy(Post_temp[iw],Post_temp[iw+1], (*it)->R.begin());
      for(int iat=0; iat < nptcl; iat++){
        (*it)->R(iat) = Pos_temp(iw,iat);
      }
      ++it;
    }
    nw_in += nwt; 
  }

  LOGMSG("Total " << nw_in << " walkers are loaded using " << NumSets-firstConf << " blocks.")
  return true;
}

 void 
 TrialWaveFunction::evaluateOptimizableLog (MCWalkerConfiguration &W,  
					    vector<RealType>& logpsi_opt,  
					    GradMatrix_t&  optG,
					    ValueMatrix_t& optL)
 {
   for (int iw=0; iw<W.getActiveWalkers(); iw++) 
     logpsi_opt[iw] = RealType();
   optG = GradType();
   optL = RealType();
   
   // Sum optimizable part of log Psi
   for (int i=0,ii=RECOMPUTE_TIMER; i<Z.size(); i++,ii+=TIMER_SKIP) {
     myTimers[ii]->start();
     if (Z[i]->Optimizable) {
       Z[i]->addLog(W, logpsi_opt);
       Z[i]->gradLapl(W, optG, optL);
     }
     myTimers[ii]->stop();
   }
 }

  int WalkerControlBase::copyWalkers(MCWalkerConfiguration& W) {
    //clear the WalkerList to populate them with the good walkers
    W.clear();
    W.insert(W.begin(), good_w.begin(), good_w.end());

    int cur_walker = good_w.size();
    for(int i=0; i<good_w.size(); i++) { //,ie+=ncols) {
      for(int j=0; j<ncopy_w[i]; j++, cur_walker++) 
      {
        Walker_t* awalker=new Walker_t(*(good_w[i]));
        awalker->ID=(++NumWalkersCreated)*NumContexts+MyContext;
        awalker->ParentID=good_w[i]->ParentID;
        W.push_back(awalker);
      }
    }

    //clear good_w and ncopy_w for the next branch
    good_w.clear();
    ncopy_w.clear();
    return W.getActiveWalkers();
  }

int WalkerControlBase::doNotBranch(int iter, MCWalkerConfiguration& W)
{
    MCWalkerConfiguration::iterator it(W.begin());
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    for(; it!=it_end; ++it)
    {
        r2_accepted+=(*it)->Properties(R2ACCEPTED);
        r2_proposed+=(*it)->Properties(R2PROPOSED);
        RealType e((*it)->Properties(LOCALENERGY));
        int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
        RealType wgt((*it)->Weight);
        esum += wgt*e;
        e2sum += wgt*e*e;
        wsum += wgt;
        w2sum += wgt*wgt;
        ecum += e;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;

    myComm->allreduce(curData);

    measureProperties(iter);
    trialEnergy=EnsembleProperty.Energy;
    W.EnsembleProperty=EnsembleProperty;

    //return the current data
    return W.getGlobalNumWalkers();
}

  int WalkerControlBase::doNotBranch(int iter, MCWalkerConfiguration& W)
  {
    MCWalkerConfiguration::iterator it(W.begin());
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0, besum=0.0, bwgtsum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    int nrn(0),ncr(0),nfn(0),ngoodfn(0);
    for(; it!=it_end;++it)
    {
      bool inFN=(((*it)->ReleasedNodeAge)==0);
      int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
      
      if ((*it)->ReleasedNodeAge==1) ncr+=1;
      else if ((*it)->ReleasedNodeAge==0) 
      {
        nfn+=1;
        ngoodfn+=nc;
      }
      
      r2_accepted+=(*it)->Properties(R2ACCEPTED);
      r2_proposed+=(*it)->Properties(R2PROPOSED);
      RealType e((*it)->Properties(LOCALENERGY));
      RealType bfe((*it)->Properties(ALTERNATEENERGY));
      RealType rnwgt(0.0);
      if (inFN)
        rnwgt=((*it)->Properties(SIGN));
      else
        nrn+=nc;
      //       RealType wgt((*it)->Weight);
      RealType wgt(0.0);
      if (inFN)
        wgt=((*it)->Weight); 
      
      esum += wgt*e;
      e2sum += wgt*e*e;
      wsum += wgt;
      w2sum += wgt*wgt;
      ecum += e;
      besum += bfe*rnwgt*wgt;
      bwgtsum += rnwgt*wgt;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;
    curData[FNSIZE_INDEX]=static_cast<RealType>(nfn);
    curData[RNONESIZE_INDEX]=static_cast<RealType>(ncr);
    curData[RNSIZE_INDEX]=nrn;
    curData[B_ENERGY_INDEX]=besum;
    curData[B_WGT_INDEX]=bwgtsum;

    myComm->allreduce(curData);

    measureProperties(iter);
    trialEnergy=EnsembleProperty.Energy;
    W.EnsembleProperty=EnsembleProperty;

    return W.getActiveWalkers();
  }

  void SimpleFixedNodeBranch::initWalkerController(MCWalkerConfiguration& walkers, RealType tau, bool fixW) 
  {

    vParam[B_TAU]=tau;
    if(!BranchMode[B_DMCSTAGE])
      vParam[B_TAUEFF]=tau*R2Accepted.result()/R2Proposed.result();
    
    if(WalkerController == 0) 
    {
      if(iParam[B_TARGETWALKERS]==0) 
      {
        Communicate* acomm=MyEstimator->getCommunicator();
        int ncontexts=acomm->size();
        vector<int> nw(ncontexts,0),nwoff(ncontexts+1,0);
        nw[acomm->rank()]=walkers.getActiveWalkers();
        acomm->allreduce(nw);
        for(int ip=0; ip<ncontexts; ++ip) nwoff[ip+1]=nwoff[ip]+nw[ip];
        walkers.setGlobalNumWalkers(nwoff[ncontexts]);
        walkers.setWalkerOffsets(nwoff);
        iParam[B_TARGETWALKERS]=nwoff[ncontexts];
      }


      BranchMode.set(B_DMC,1);//set DMC
      BranchMode.set(B_POPCONTROL,!fixW);//fixW -> 0 
      WalkerController = createWalkerController(iParam[B_TARGETWALKERS], MyEstimator->getCommunicator(), myNode);
      iParam[B_MAXWALKERS]=WalkerController->Nmax;
      iParam[B_MINWALKERS]=WalkerController->Nmin;

      if(!fixW && sParam[MIXDMCOPT]=="yes")
      {
        app_log() << "Warmup DMC is done with a fixed population " << iParam[B_TARGETWALKERS] << endl;
        BackupWalkerController=WalkerController; //save the main controller
        WalkerController=createWalkerController(iParam[B_TARGETWALKERS],MyEstimator->getCommunicator(), myNode,true);
        BranchMode.set(B_POPCONTROL,0);
      }
      WalkerController->setWalkerID(walkers);
      PopHist.clear();
      PopHist.reserve(std::max(iParam[B_ENERGYUPDATEINTERVAL],5));
    }

    //save the BranchMode in anticipating state changes in reset
    bitset<B_MODE_MAX> bmode(BranchMode);
    //reset Feedback pararmeter
    this->reset();

    MyEstimator->reset();
    //update the simulation parameters
    WalkerController->put(myNode);
    //assign current Eref and a large number for variance
    WalkerController->setEnergyAndVariance(vParam[B_EREF],vParam[B_SIGMA]);

    //determine the branch cutoff to limit wild weights based on the sigma and sigmaBound
    RealType sigma=std::max(std::sqrt(static_cast<RealType>(iParam[B_TARGETWALKERS]))*vParam[B_SIGMA]*WalkerController->targetSigma,
			    static_cast<RealType>(100.0));
    vParam[B_BRANCHCUTOFF]=std::min(sigma,5.0/tau);
    //vParam[B_BRANCHCUTOFF]=vParam[B_SIGMA]*WalkerController->targetSigma;
    vParam[B_BRANCHMAX]=vParam[B_BRANCHCUTOFF]*1.5;
    vParam[B_BRANCHFILTER]=1.0/(vParam[B_BRANCHMAX]-vParam[B_BRANCHCUTOFF]);

    //reset controller 
    WalkerController->reset();
    if(BackupWalkerController) BackupWalkerController->reset();
    BranchMode=bmode;

    app_log() << "  QMC counter      = " << iParam[B_COUNTER] << endl;
    app_log() << "  time step        = " << vParam[B_TAU] << endl;
    app_log() << "  effective time step = " << vParam[B_TAUEFF] << endl;
    app_log() << "  trial energy     = " << vParam[B_ETRIAL] << endl;
    app_log() << "  reference energy = " << vParam[B_EREF] << endl;
    app_log() << "  Feedback = " << Feedback <<  endl;
    app_log() << "  reference variance = " << vParam[B_SIGMA] << endl;
    app_log() << "  target walkers = " << iParam[B_TARGETWALKERS] << endl;
    app_log() << "  branch cutoff = " <<  vParam[B_BRANCHCUTOFF] << " " << vParam[B_BRANCHMAX] << endl;
    app_log() << "  Max and mimum walkers per node= " << iParam[B_MAXWALKERS] << " " << iParam[B_MINWALKERS] << endl;
    app_log() << "  QMC Status (BranchMode) = " << BranchMode << endl;
  }

//void SimpleFixedNodeBranch::initWalkerController(MCWalkerConfiguration& walkers, RealType tau, bool fixW, bool killwalker) 
  void SimpleFixedNodeBranch::initWalkerController(MCWalkerConfiguration& walkers, bool fixW, bool killwalker) 
  {
    BranchMode.set(B_DMC,1);//set DMC
    BranchMode.set(B_DMCSTAGE,iParam[B_WARMUPSTEPS]==0);//use warmup

    //this is not necessary
    //check if tau is different and set the initial values
    //vParam[B_TAU]=tau;

    bool fromscratch=false;
    RealType tau=vParam[B_TAU];

    //this is the first time DMC is used
    if(WalkerController == 0) 
    {
      if(iParam[B_TARGETWALKERS]==0) 
      {
        Communicate* acomm=MyEstimator->getCommunicator();
        int ncontexts=acomm->size();
        vector<int> nw(ncontexts,0),nwoff(ncontexts+1,0);
        nw[acomm->rank()]=walkers.getActiveWalkers();
        acomm->allreduce(nw);
        for(int ip=0; ip<ncontexts; ++ip) nwoff[ip+1]=nwoff[ip]+nw[ip];
        walkers.setGlobalNumWalkers(nwoff[ncontexts]);
        walkers.setWalkerOffsets(nwoff);
        iParam[B_TARGETWALKERS]=nwoff[ncontexts];
      }

      WalkerController = createWalkerController(iParam[B_TARGETWALKERS], MyEstimator->getCommunicator(), myNode);

      if(!BranchMode[B_RESTART])
      {
        fromscratch=true;
        app_log() << "  START ALL OVER " << endl;
        vParam[B_TAUEFF]=tau;
        BranchMode.set(B_POPCONTROL,!fixW);//fixW -> 0 
        BranchMode.set(B_KILLNODES,killwalker);
        iParam[B_MAXWALKERS]=WalkerController->Nmax;
        iParam[B_MINWALKERS]=WalkerController->Nmin;

        if(!fixW && sParam[MIXDMCOPT]=="yes")
        {
          app_log() << "Warmup DMC is done with a fixed population " << iParam[B_TARGETWALKERS] << endl;
          BackupWalkerController=WalkerController; //save the main controller
          WalkerController=createWalkerController(iParam[B_TARGETWALKERS],MyEstimator->getCommunicator(), myNode,true);
          BranchMode.set(B_POPCONTROL,0);
        }
        //PopHist.clear();
        //PopHist.reserve(std::max(iParam[B_ENERGYUPDATEINTERVAL],5));
      }
      WalkerController->setWalkerID(walkers);
    }
    //else
    //{
    //  BranchMode.set(B_DMCSTAGE,0);//always reset warmup
    //}

    MyEstimator->reset();
    //update the simulation parameters
    WalkerController->put(myNode);
    //assign current Eref and a large number for variance
    WalkerController->setEnergyAndVariance(vParam[B_EREF],vParam[B_SIGMA]);

    this->reset();

    if(fromscratch)
    {
      //determine the branch cutoff to limit wild weights based on the sigma and sigmaBound
      //RealType sigma=std::max(std::sqrt(static_cast<RealType>(iParam[B_TARGETWALKERS]))*vParam[B_SIGMA]*WalkerController->targetSigma,100.0);
      RealType sigma=std::max(std::sqrt(vParam[B_SIGMA])*WalkerController->targetSigma,50.0);
      vParam[B_BRANCHCUTOFF]=std::min(sigma,2.5/tau);
      //vParam[B_BRANCHCUTOFF]=vParam[B_SIGMA]*WalkerController->targetSigma;
      vParam[B_BRANCHMAX]=vParam[B_BRANCHCUTOFF]*1.5;
      vParam[B_BRANCHFILTER]=1.0/(vParam[B_BRANCHMAX]-vParam[B_BRANCHCUTOFF]);
      vParam[B_TAUEFF]=tau*R2Accepted.result()/R2Proposed.result();
    }

    //reset controller 
    WalkerController->reset();
    if(BackupWalkerController) BackupWalkerController->reset();

    app_log() << "  QMC counter      = " << iParam[B_COUNTER] << endl;
    app_log() << "  time step        = " << vParam[B_TAU] << endl;
    app_log() << "  effective time step = " << vParam[B_TAUEFF] << endl;
    app_log() << "  trial energy     = " << vParam[B_ETRIAL] << endl;
    app_log() << "  reference energy = " << vParam[B_EREF] << endl;
    app_log() << "  Feedback = " << vParam[B_FEEDBACK] <<  endl;
    app_log() << "  reference variance = " << vParam[B_SIGMA] << endl;
    app_log() << "  target walkers = " << iParam[B_TARGETWALKERS] << endl;
    app_log() << "  branch cutoff = " <<  vParam[B_BRANCHCUTOFF] << " " << vParam[B_BRANCHMAX] << endl;
    app_log() << "  Max and mimum walkers per node= " << iParam[B_MAXWALKERS] << " " << iParam[B_MINWALKERS] << endl;
    app_log() << "  QMC Status (BranchMode) = " << BranchMode << endl;
  }

int WalkerReconfigurationMPI::swapWalkers(MCWalkerConfiguration& W) {
  //ostringstream o;
  //o << "check." << MyContext << ".dat";
  //ofstream fout(o.str().c_str(),ios::app);

  int nw=W.getActiveWalkers();
  if(TotalWalkers ==0) {
    FirstWalker=nw*MyContext;
    LastWalker=FirstWalker+nw;
    TotalWalkers = nw*NumContexts;
    nwInv = 1.0/static_cast<RealType>(TotalWalkers);
    DeltaStep = UnitZeta*nwInv;

    ncopy_w.resize(nw);
    wConf.resize(nw);
    //wSum.resize(NumContexts);
    wOffset.resize(NumContexts+1);
    dN.resize(NumContexts+1);
  }

  //std::fill(wSum.begin(),wSum.end(),0.0);

  MCWalkerConfiguration::iterator it(W.begin()), it_end(W.end());
  int iw=0;
  RealType esum=0.0,e2sum=0.0,wtot=0.0,ecum=0.0;
  while(it != it_end) {
    RealType wgt((*it)->Weight);
    RealType e((*it)->Properties(LOCALENERGY));
    esum += wgt*e;
    e2sum += wgt*e*e;
    wtot += wgt;
    ecum += e;
    wConf[iw++]=wgt;
    ++it;
  }
  //wSum[MyContext]=wtot;
  curData[ENERGY_INDEX]=esum;
  curData[ENERGY_SQ_INDEX]=e2sum;
  curData[WALKERSIZE_INDEX]=nw;
  curData[WEIGHT_INDEX]=wtot;
  curData[EREF_INDEX]=ecum;
  std::fill(curData.begin()+LE_MAX,curData.end(),0.0);
  curData[LE_MAX+MyContext]=wtot;

  //collect everything
  myComm->allreduce(curData);

  //wOffset[ip] is the partial sum update to ip
  wOffset[0]=0;
  //for(int ip=0; ip<NumContexts; ip++) wOffset[ip+1]=wOffset[ip]+wSum[ip];
  for(int ip=0,jp=LE_MAX; ip<NumContexts; ip++,jp++) 
    wOffset[ip+1]=wOffset[ip]+curData[jp];

  wtot=wOffset[NumContexts]; //wtot is the total weight

  //find the lower and upper bound of index
  int minIndex=static_cast<int>((wOffset[MyContext]/wtot-DeltaStep)*static_cast<RealType>(TotalWalkers))-1;
  int maxIndex=static_cast<int>((wOffset[MyContext+1]/wtot-DeltaStep)*static_cast<RealType>(TotalWalkers))+1;
  int nb=maxIndex-minIndex+1;
  vector<RealType> Zeta(nb);

  for(int i=minIndex, ii=0; i<maxIndex; i++,ii++) {
    Zeta[ii]= wtot*(DeltaStep+static_cast<RealType>(i)*nwInv);
  }

  RealType wCur=wOffset[MyContext];
  int ind=0;
  while(Zeta[ind]<wCur) {ind++;} 

  //surviving walkers
  int icdiff=0;
  for(iw=0; iw<nw; iw++) {
    RealType tryp=wCur+abs(wConf[iw]);
    int ni=0;
    while(Zeta[ind]<tryp && Zeta[ind] >= wCur) {
      ind++;
      ni++;
    }
    wCur+=abs(wConf[iw]);
    if(ni) {
      icdiff++;
    } 
    ncopy_w[iw]=ni;
  }

  vector<int> plus, minus;
  for(iw=0;iw<nw; iw++) {
    int m=ncopy_w[iw];
    if(m>1) {// add the index of this walker to plus, duplicate m-1 times
      plus.insert(plus.end(),m-1,iw);
    } else if(m==0) { // add the walker index to be killed/overwritten
      minus.push_back(iw);
    }
  }

  //copy within the local node
  int lower=std::min(plus.size(),minus.size()); 
  while(lower>0) {
    --lower;
    W[minus[lower]]->assign(*(W[plus[lower]]));
//.........这里部分代码省略.........

  // old ones
  void 
    RQMCEstimator
    ::initialize(MCWalkerConfiguration& W, 
        vector<QMCHamiltonian*>& h, 
        vector<TrialWaveFunction*>& psi,
        RealType tau,vector<RealType>& Norm,
        bool require_register) {

      NumWalkers = W.getActiveWalkers();
      //allocate UmbrellaEnergy
      int numPtcls(W.getTotalNum());
      RatioIJ.resize(NumWalkers,NumCopies*(NumCopies-1)/2);

      MCWalkerConfiguration::iterator it(W.begin()); 
      MCWalkerConfiguration::iterator it_end(W.end()); 

      vector<RealType> sumratio(NumCopies), logpsi(NumCopies);
      int iw(0);
      int DataSetSize((*it)->DataSet.size());
      while(it != it_end) {

        Walker_t& thisWalker(**it);
        (*it)->DataSet.rewind();

        if(require_register) {
          W.registerData(thisWalker,(*it)->DataSet);
        } else {
          W.R = thisWalker.R;
          W.update();
          if(DataSetSize) W.updateBuffer((*it)->DataSet);
        }

        //evalaute the wavefunction and hamiltonian
        for(int ipsi=0; ipsi< NumCopies;ipsi++) {
          psi[ipsi]->G.resize(numPtcls);
          psi[ipsi]->L.resize(numPtcls);
          //Need to modify the return value of OrbitalBase::registerData
          if(require_register) {
            logpsi[ipsi]=psi[ipsi]->registerData(W,(*it)->DataSet);
          } else {
            if(DataSetSize)logpsi[ipsi]=psi[ipsi]->updateBuffer(W,(*it)->DataSet);
            else logpsi[ipsi]=psi[ipsi]->evaluateLog(W); 		 
          }
          psi[ipsi]->G=W.G;
          thisWalker.Properties(ipsi,LOGPSI)=logpsi[ipsi];
          thisWalker.Properties(ipsi,LOCALENERGY)=h[ipsi]->evaluate(W);
          h[ipsi]->saveProperty(thisWalker.getPropertyBase(ipsi));
          sumratio[ipsi]=1.0;
        } 							

        //Check SIMONE's note
        //Compute the sum over j of Psi^2[j]/Psi^2[i] for each i
        int indexij(0);
        RealType *rPtr=RatioIJ[iw];
        for(int ipsi=0; ipsi< NumCopies-1; ipsi++) {			  
          for(int jpsi=ipsi+1; jpsi< NumCopies; jpsi++){     		 
            RealType r= std::exp(2.0*(logpsi[jpsi]-logpsi[ipsi])); 
            rPtr[indexij++]=r*Norm[ipsi]/Norm[jpsi];
            sumratio[ipsi] += r;                            
            sumratio[jpsi] += 1.0/r;		
          }                                              
        }                                               

        //Re-use Multiplicity as the sumratio
        thisWalker.Multiplicity=sumratio[0];
        //DON't forget DRIFT!!!
        thisWalker.Drift=0.0;

        for(int ipsi=0; ipsi< NumCopies; ipsi++) {
          RealType wgt=1.0/sumratio[ipsi];
          thisWalker.Properties(ipsi,UMBRELLAWEIGHT)=wgt;

          //thisWalker.Drift += wgt*psi[ipsi]->G;
          PAOps<RealType,DIM>::axpy(wgt,psi[ipsi]->G,thisWalker.Drift);
        }
        thisWalker.Drift *= tau;
        ++it;++iw;
      }
    }

/** swap walkers using (low,high) ordered pairs
 *
 * The communication occur only between the (low,high) pairs.
 * This does not guarantee a perfect balance swapWalkersBlocked and swapWalkersAsync
 * try to achieve. However, the number of messages and their sizes are less than
 * other methods.
 */
void GlobalWalkerControl::swapWalkersMap(MCWalkerConfiguration& W) {

  NumSwaps++;
  multimap<int,int> nw_map;
  for(int i=0; i<NumContexts; i++) {
    nw_map.insert(pair<int,int>(NumPerNode[i],i));
  }
  // multimap key is sorted with ascending order 
  multimap<int,int>::iterator it(nw_map.begin());
  multimap<int,int>::reverse_iterator it_b(nw_map.end());
  bool notpaired=true;
  int target_context=-1;
  int half=NumContexts/2;
  int item=0;
  bool minorcontext;
  while(notpaired &&item<half) {
    int i=(*it).second;
    int j=(*it_b).second;
    if(i == MyContext) {
      target_context=j;
      notpaired=false;
      minorcontext=true;
    } else if(j == MyContext) {
      target_context= i;
      notpaired=false;
      minorcontext=false;
    } 
    ++it; ++it_b; ++item;
  }

  int nw_tot=NumPerNode[MyContext]+NumPerNode[target_context];
  int nw_L=nw_tot/2;
  int nw_R=nw_tot-nw_L;
  int dnw(0);
  if(minorcontext) {
    dnw=nw_L-NumPerNode[MyContext];//how many to get
  } else {
    dnw=NumPerNode[MyContext]-nw_R;//how many to send
  }

  if(dnw) {//something to swap
    if(minorcontext) {//open recv buffer
      Walker_t& wRef(*W[0]);
      OOMPI_Packed recvBuffer(dnw*wRef.byteSize(),OOMPI_COMM_WORLD);
      OOMPI_COMM_WORLD[target_context].Recv(recvBuffer);
      //create walkers
      int last = W.getActiveWalkers();
      while(dnw) {
        Walker_t *awalker= new Walker_t(wRef);
        awalker->getMessage(recvBuffer);
        W.push_back(awalker);
        --dnw; ++last;
      }
    } else {
      Walker_t& wRef(*W[0]);
      OOMPI_Packed sendBuffer(dnw*wRef.byteSize(),OOMPI_COMM_WORLD);
      int last=W.getActiveWalkers()-1;
      while(dnw) {
        W[last]->putMessage(sendBuffer);
        --dnw; --last;
      }
      OOMPI_COMM_WORLD[target_context].Send(sendBuffer);
      //last=WalkerList.size()-1;
      //while(dnw_save) {
      //  delete WalkerList[last]; 
      //  WalkerList.pop_back();
      //  --dnw_save; --last;
      //}
      //destroyWalkers(WalkerList.begin()+nsub[MyContext], WalkerList.end());
      W.destroyWalkers(W.begin()+nw_R, W.end());
    }
  }

  /* not used yet
  struct lessNode {
    inline bool operator()(const pair<int,int>& a,
        const pair<int,int>& b) const {
      return a.second < b.second;
    }
  };
  //typedef pair<int,int> mytype;
  //vector<mytype> id2n(NumContexts);
  //for(int i=0; i<NumContexts; i++) {
  //  id2n=pair<int,int>(i,NumPerNode[i]);
  //}
  //
  //std::sort(id2n.begin(),id2n.end(),lessNode);
  */
}

/** evaluate curData and mark the bad/good walkers
 */
void WalkerControlBase::sortWalkers(MCWalkerConfiguration& W) {

    MCWalkerConfiguration::iterator it(W.begin());

    vector<Walker_t*> bad;
    NumWalkers=0;
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    //RealType sigma=std::max(5.0*targetVar,targetEnergyBound);
    //RealType ebar= targetAvg;
    while(it != it_end)
    {
        r2_accepted+=(*it)->Properties(R2ACCEPTED);
        r2_proposed+=(*it)->Properties(R2PROPOSED);
        RealType e((*it)->Properties(LOCALENERGY));
        int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
        /// HACK HACK HACK
        //RealType wgt((*it)->Weight);
        RealType wgt = (RealType)nc;

        esum += wgt*e;
        e2sum += wgt*e*e;
        wsum += wgt;
        w2sum += wgt*wgt;
        ecum += e;

        if(nc)
        {
            NumWalkers += nc;
            good_w.push_back(*it);
            ncopy_w.push_back(nc-1);
        } else {
            bad.push_back(*it);
        }
        ++it;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    //evaluate variance of this block
    //curVar=(e2sum-esum*esum/wsum)/wsum;
    //THIS IS NOT USED ANYMORE:BELOW
    //if(curVar>sigma) {
    //  app_error() << "Unphysical block variance is detected. Stop simulations." << endl;
    //  Write2XYZ(W);
    //  //Does not work some reason
    //  //OHMMS::Controller->abort();
    //#if defined(HAVE_MPI)
    //      OOMPI_COMM_WORLD.Abort();
    //#else
    //      abort();
    //#endif
    //    }
    //THIS IS NOT USED ANYMORE:ABOVE

    //update curData
    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;

    ////this should be move
    //W.EnsembleProperty.NumSamples=curData[WALKERSIZE_INDEX];
    //W.EnsembleProperty.Weight=curData[WEIGHT_INDEX];
    //W.EnsembleProperty.Energy=(esum/=wsum);
    //W.EnsembleProperty.Variance=(e2sum/wsum-esum*esum);
    //W.EnsembleProperty.Variance=(e2sum*wsum-esum*esum)/(wsum*wsum-w2sum);

    //remove bad walkers empty the container
    for(int i=0; i<bad.size(); i++) delete bad[i];

    if(good_w.empty()) {
        app_error() << "All the walkers have died. Abort. " << endl;
        APP_ABORT("WalkerControlBase::sortWalkers");
    }

    int sizeofgood = good_w.size();

    //check if the projected number of walkers is too small or too
    //large
    if(NumWalkers>Nmax) {
        int nsub=0;
        int nsub_target=NumWalkers-static_cast<int>(0.9*Nmax);
        int i=0;
        while(i< sizeofgood && nsub<nsub_target) {
            if(ncopy_w[i]) {
                ncopy_w[i]--;
                nsub++;
            }
            ++i;
        }
        NumWalkers -= nsub;
    } else  if(NumWalkers < Nmin) {
//.........这里部分代码省略.........

  /** evaluate curData and mark the bad/good walkers
   */
  int WalkerControlBase::sortWalkers(MCWalkerConfiguration& W) {

    MCWalkerConfiguration::iterator it(W.begin());

    vector<Walker_t*> bad,good_rn;
    vector<int> ncopy_rn;
    NumWalkers=0;
    MCWalkerConfiguration::iterator it_end(W.end());
    RealType esum=0.0,e2sum=0.0,wsum=0.0,ecum=0.0, w2sum=0.0, besum=0.0, bwgtsum=0.0;
    RealType r2_accepted=0.0,r2_proposed=0.0;
    int nrn(0),ncr(0);
    while(it != it_end) 
    {
      bool inFN=(((*it)->ReleasedNodeAge)==0);
      if ((*it)->ReleasedNodeAge==1) ncr+=1;
      int nc= std::min(static_cast<int>((*it)->Multiplicity),MaxCopy);
      r2_accepted+=(*it)->Properties(R2ACCEPTED);
      r2_proposed+=(*it)->Properties(R2PROPOSED);
      RealType e((*it)->Properties(LOCALENERGY));
      RealType bfe((*it)->Properties(ALTERNATEENERGY));
      RealType rnwgt(0.0);
      if (inFN)
        rnwgt=((*it)->Properties(SIGN));
      
//       RealType wgt((*it)->Weight);
      RealType wgt(0.0);
      if (inFN)
        wgt=((*it)->Weight); 
      
      esum += wgt*e;
      e2sum += wgt*e*e;
      wsum += wgt;
      w2sum += wgt*wgt;
      ecum += e;
      besum += bfe*rnwgt*wgt;
      bwgtsum += rnwgt*wgt;

      if((nc) && (inFN))
      {
        NumWalkers += nc;
        good_w.push_back(*it);
        ncopy_w.push_back(nc-1);
      }
      else if (nc)
      {
        NumWalkers += nc;
        nrn+=nc;
        good_rn.push_back(*it);
        ncopy_rn.push_back(nc-1);
      }
      else
      {
        bad.push_back(*it);
      }
      ++it;
    }

    //temp is an array to perform reduction operations
    std::fill(curData.begin(),curData.end(),0);

    //update curData
    curData[ENERGY_INDEX]=esum;
    curData[ENERGY_SQ_INDEX]=e2sum;
    curData[WALKERSIZE_INDEX]=W.getActiveWalkers();
    curData[WEIGHT_INDEX]=wsum;
    curData[EREF_INDEX]=ecum;
    curData[R2ACCEPTED_INDEX]=r2_accepted;
    curData[R2PROPOSED_INDEX]=r2_proposed;
    curData[FNSIZE_INDEX]=static_cast<RealType>(good_w.size());
    curData[RNONESIZE_INDEX]=static_cast<RealType>(ncr);
    curData[RNSIZE_INDEX]=nrn;
    curData[B_ENERGY_INDEX]=besum;
    curData[B_WGT_INDEX]=bwgtsum;
    
    ////this should be move
    //W.EnsembleProperty.NumSamples=curData[WALKERSIZE_INDEX];
    //W.EnsembleProperty.Weight=curData[WEIGHT_INDEX];
    //W.EnsembleProperty.Energy=(esum/=wsum);
    //W.EnsembleProperty.Variance=(e2sum/wsum-esum*esum);
    //W.EnsembleProperty.Variance=(e2sum*wsum-esum*esum)/(wsum*wsum-w2sum);

    //remove bad walkers empty the container
    for(int i=0; i<bad.size(); i++) delete bad[i];
    if (!WriteRN)
    {
     if(good_w.empty()) {
      app_error() << "All the walkers have died. Abort. " << endl;
      APP_ABORT("WalkerControlBase::sortWalkers");
     }

     int sizeofgood = good_w.size();
     //check if the projected number of walkers is too small or too large
     if(NumWalkers>Nmax) {
      int nsub=0;
      int nsub_target=(NumWalkers-nrn)-static_cast<int>(0.9*Nmax);
      int i=0;
      while(i< sizeofgood && nsub<nsub_target) {
        if(ncopy_w[i]) {ncopy_w[i]--; nsub++;}
//.........这里部分代码省略.........