C++ QMCHamiltonian类代码示例

  /** constructor
   * @param h QMCHamiltonian to define the components
   * @param hcopy number of copies of QMCHamiltonians
   */
  RQMCEstimator::RQMCEstimator(QMCHamiltonian& h, int hcopy) : CurrentWalker(0) {

    NumCopies=hcopy;
    NumOperators = h.sizeOfObservables();
    FirstHamiltonian=h.startIndex();

    esum.resize(NumCopies,LE_INDEX);
    esum_name.resize(LE_INDEX,NumCopies);

    elocal.resize(NumCopies,NumOperators);
    elocal_name.resize(NumCopies,NumOperators);

    char aname[32];
    //make the name tables
    //(localenergy, variance, weight)* + (each hamiltonian term)*
    for(int i=0; i<NumCopies; i++) {
      sprintf(aname,"LE%i",i);   esum_name(ENERGY_INDEX,i)=aname;
      sprintf(aname,"Var%i",i);  esum_name(ENERGY_SQ_INDEX,i)=aname;
      sprintf(aname,"WPsi%i",i); esum_name(WEIGHT_INDEX,i)=aname;

      for(int j=0; j<NumOperators; j++) {
        //sprintf(aname,"%s%i",h.getName(j).c_str(),i);
        sprintf(aname,"%s%i",h.getObservableName(j).c_str(),i);
        elocal_name(i,j)=aname;
      }
    }

    int ipair=0;
    for(int i=0; i<NumCopies-1; i++) {
      for(int j=i+1; j<NumCopies; j++) {
        sprintf(aname,"DiffS%iS%i",i,j); ediff_name.push_back(aname);
      }
    }
  }

  /** constructor
   * @param h QMCHamiltonian to define the components
   * @param hcopy number of copies of QMCHamiltonians
   */
  HFDHE2PolymerEstimator::HFDHE2PolymerEstimator(QMCHamiltonian& h, int hcopy, MultiChain* polymer)
  {
    Reptile=polymer;
    NumCopies=hcopy;
//     Hpointer = &h;

    SizeOfHamiltonians = h.sizeOfObservables();
    FirstHamiltonian = h.startIndex();
//     cout<<"size of Hamiltonian "<<SizeOfHamiltonians<<" First one "<<FirstHamiltonian<<endl;

    elocal_name.push_back("LocalEnergy");
    for(int i=0; i < SizeOfHamiltonians; i++)
    {
      elocal_name.push_back(h.getObservableName(i));
      //elocal_name.push_back(h.getName(i));
    };
    
    elocal_name.push_back("SumPot");
    elocal_name.push_back("ElSumPot");
    elocal_name.push_back("Virial");
    elocal_name.push_back("MaxAge");

    scalars.resize(SizeOfHamiltonians+5);
    scalars_saved=scalars;
    pNorm=0.0;
    
    Pindex = h.getHamiltonian("HePress")->myIndex;
    HDFHE2index = h.getHamiltonian("HFDHE2")->myIndex;
  }

 LocalEnergyEstimatorHDF::LocalEnergyEstimatorHDF(QMCHamiltonian& h)
   :refH(h)
 { 
   SizeOfHamiltonians = h.sizeOfObservables();
   FirstHamiltonian = h.startIndex();
   scalars.resize(SizeOfHamiltonians+LE_MAX);
   scalars_saved.resize(SizeOfHamiltonians+LE_MAX);
 }

 LocalEnergyEstimator::LocalEnergyEstimator(QMCHamiltonian& h):Href(h) { 
   int hterms(h.size());
   SizeOfHamiltonians = hterms;
   FirstHamiltonian = h.startIndex();
   elocal.resize(SizeOfHamiltonians+LE_MAX);
   elocal_name.resize(SizeOfHamiltonians+LE_MAX);
   elocal_name[ENERGY_INDEX] = "LocalEnergy";
   elocal_name[ENERGY_SQ_INDEX] = "Variance";
   elocal_name[POTENTIAL_INDEX] = "LocalPotential";
   int ii(LE_MAX);
   for(int i=0; i < SizeOfHamiltonians; i++) elocal_name[ii++] = h.getName(i);
 }

  void CloneManager::makeClones(MCWalkerConfiguration& w, 
      TrialWaveFunction& psi, QMCHamiltonian& ham)
  {

    if(wClones.size()) {
      app_log() << "  Cannot make clones again. Use existing " << NumThreads << " clones" << endl;
      return;
    }

    app_log() << "Number of threads = " << NumThreads << endl;
    wClones.resize(NumThreads,0);
    psiClones.resize(NumThreads,0);
    hClones.resize(NumThreads,0);

    wClones[0]=&w;
    psiClones[0]=&psi;
    hClones[0]=&ham;

#if defined(ENABLE_CLONE_PSI_AND_H)
    char pname[16];
    for(int ip=1; ip<NumThreads; ++ip) 
    {
      sprintf(pname,"%s.c%i",w.getName().c_str(),ip);
      wClones[ip]=new MCWalkerConfiguration(w);
      wClones[ip]->setName(pname);
      psiClones[ip]=psi.makeClone(*wClones[ip]);
      hClones[ip]=ham.makeClone(*wClones[ip],*psiClones[ip]);
    }
#else
    cloneEngine.clone(w,psi,ham,wClones,psiClones,hClones);
#endif
  }

  /** constructor
   * @param h QMCHamiltonian to define the components
   * @param hcopy number of copies of QMCHamiltonians
   */
  CSPolymerEstimator::CSPolymerEstimator(QMCHamiltonian& h, int hcopy, 
      MultiChain* polymer)
//   : Reptile(polymer)
  {
    Reptile = polymer;
    NumCopies=hcopy;
    NumObservables = h.size();
    scalars.resize(NumCopies+NumCopies*(NumCopies-1)/2);
    scalars_saved=scalars;
    //d_data.resize(NumCopies*3+NumCopies*(NumCopies-1)/2);
  }

/** constructor
 * @param h QMCHamiltonian to define the components
 * @param hcopy number of copies of QMCHamiltonians
 */
MJPolymerEstimator::MJPolymerEstimator(QMCHamiltonian& h, int hcopy, MultiChain* polymer)
{
    Reptile=polymer;
    NumCopies=hcopy;
//     Hpointer = &h;

    SizeOfHamiltonians = h.sizeOfObservables();
    FirstHamiltonian = h.startIndex();
//     cout<<"size of Hamiltonian "<<SizeOfHamiltonians<<" First one "<<FirstHamiltonian<<endl;

    elocal_name.push_back("LocalEnergy");
    for(int i=0; i < SizeOfHamiltonians; i++)
    {
        elocal_name.push_back(h.getObservableName(i));
        //elocal_name.push_back(h.getName(i));
    };

    elocal_name.push_back("SumPot");
    elocal_name.push_back("ElSumPot");
    elocal_name.push_back("CenterTruncSumPot");
    elocal_name.push_back("Null");
    elocal_name.push_back("CenterTruncElSumPot");
    elocal_name.push_back("Virial");
    elocal_name.push_back("MaxAge");
    elocal_name.push_back("MaxTouch");
    elocal_name.push_back("centerV");
    elocal_name.push_back("EcorrFun");
    elocal_name.push_back("Ehead");
    elocal_name.push_back("Etail");
    elocal_name.push_back("Edecorr");
    elocal_name.push_back("Vdecorr");
    elocal_name.push_back("RMC_HFCep_1_0");
    elocal_name.push_back("RMC_HFCep_1_1");

    scalars.resize(SizeOfHamiltonians+23);
    scalars_saved=scalars;
    pNorm=0.0;
    Findex = h.getObservable("HFCep_1_0");
//     app_log()<<"Force Index "<<Findex<<endl;
};

/** This should be moved to branch engine */
bool EstimatorManager::put(MCWalkerConfiguration& W, QMCHamiltonian& H, xmlNodePtr cur)
{
  vector<string> extra;
  cur = cur->children;
  while(cur != NULL)
  {
    string cname((const char*)(cur->name));
    if(cname == "estimator")
    {
      string est_name(MainEstimatorName);
      string use_hdf5("yes");
      OhmmsAttributeSet hAttrib;
      hAttrib.add(est_name, "name");
      hAttrib.add(use_hdf5, "hdf5");
      hAttrib.put(cur);
      if( (est_name == MainEstimatorName) || (est_name=="elocal") )
      {
        max4ascii=H.sizeOfObservables()+3;
        add(new LocalEnergyEstimator(H,use_hdf5=="yes"),MainEstimatorName);
      }
      else
        extra.push_back(est_name);
    }
    cur = cur->next;
  }
  if(Estimators.empty())
  {
    app_log() << "  Adding a default LocalEnergyEstimator for the MainEstimator " << endl;
    max4ascii=H.sizeOfObservables()+3;
    add(new LocalEnergyEstimator(H,true),MainEstimatorName);
  }
  //Collectables is special and should not be added to Estimators
  if(Collectables == 0 && H.sizeOfCollectables())
  {
    app_log() << "  Using CollectablesEstimator for collectables, e.g. sk, gofr, density " << endl;
    Collectables=new CollectablesEstimator(H);
  }
  return true;
}

  void CloneManager::makeClones(MCWalkerConfiguration& w, 
      TrialWaveFunction& psi, QMCHamiltonian& ham)
  {

    if(wClones.size()) 
    {
      app_log() << "  Cannot make clones again. Use existing " << NumThreads << " clones" << endl;
      return;
    }

    wClones.resize(NumThreads,0);
    psiClones.resize(NumThreads,0);
    hClones.resize(NumThreads,0);

    wClones[0]=&w;
    psiClones[0]=&psi;
    hClones[0]=&ham;

    if(NumThreads==1) return;

    app_log() << "  CloneManager::makeClones makes " << NumThreads << " clones for W/Psi/H." <<endl;
#if defined(ENABLE_CLONE_PSI_AND_H)
    app_log() << "  Cloning methods for both Psi and H are used" << endl;
    OhmmsInfo::Log->turnoff();
    OhmmsInfo::Warn->turnoff();
    char pname[16];
    for(int ip=1; ip<NumThreads; ++ip) 
    {
      wClones[ip]=new MCWalkerConfiguration(w);
      psiClones[ip]=psi.makeClone(*wClones[ip]);
      hClones[ip]=ham.makeClone(*wClones[ip],*psiClones[ip]);
    }
    OhmmsInfo::Log->reset();
    OhmmsInfo::Warn->reset();
#else
    app_log() << "Old parse method is used." << endl;
    cloneEngine.clone(w,psi,ham,wClones,psiClones,hClones);
#endif
  }

  void HamiltonianPool::clone(const MCWalkerConfiguration& qp, const TrialWaveFunction& psi, const QMCHamiltonian& h,
      vector<MCWalkerConfiguration*>& plist, 
      vector<TrialWaveFunction*>& olist, 
      vector<QMCHamiltonian*>& hlist) 
  {
    int np=omp_get_max_threads();
    {
      ReportEngine PRE("HamiltonianPool","clone");
      app_log() << "Number of threads = " << np << endl;
    }

    //temporarily turnoff stream buffer for the clones
    OhmmsInfo::Log->turnoff();
    OhmmsInfo::Warn->turnoff();

    //clone ParticleSet and TrialWaveFunction
    WaveFunctionFactory* psiFac=psiPool->getWaveFunctionFactory(psi.getName());
    psiFac->setCloneSize(np);

    //capture cloned WaveFunctionFactory*
    vector<WaveFunctionFactory*> otemp;
    otemp.resize(np,0);

    //allocate the data on each thread
//#pragma omp parallel 
//    {
//      int ip=omp_get_thread_num();
//#pragma omp critical 
//      {
//        if(ip) {
//          char pname[16],oname[16];
//          sprintf(pname,"%s.c%i",qp.getName().c_str(),ip);
//          plist[ip]=new MCWalkerConfiguration(qp);
//          plist[ip]->setName(pname);
//
//          sprintf(oname,"%s.c%i",psi.getName().c_str(),ip);
//          otemp[ip]= psiFac->clone(plist[ip],ip,oname);
//        }
//      }
//    }
//    
//    //add the Clones to the pools
//    for(int ip=1; ip<np; ip++) 
//    {
//      ptclPool->addParticleSet(plist[ip]);
//      psiPool->addFactory(otemp[ip]);
//      olist[ip]=otemp[ip]->targetPsi;
//      if(ip%2==1) olist[ip]->reverse();
//    }
//
    
    //add the Clones to the pools
    for(int ip=1; ip<np; ++ip) 
    {
      plist[ip]=new MCWalkerConfiguration(qp);

      char oname[16];
      sprintf(oname,"%s.c%i",psi.getName().c_str(),ip);
      // Don't recreate with a factory anymore
      //otemp[ip]= psiFac->clone(plist[ip],ip,oname);
      //ptclPool->addParticleSet(plist[ip]);
      //psiPool->addFactory(otemp[ip]);
      //olist[ip]=otemp[ip]->targetPsi;

      // Just clone the TrialWaveFunction.
      ptclPool->addParticleSet(plist[ip]);
      olist[ip]=psiFac->targetPsi->makeClone(*plist[ip]);
      olist[ip]->setName(oname);

      //need to add a WaveFunctionFactory so that Hamiltonian can use them
      otemp[ip]=new WaveFunctionFactory(*psiFac);//make a shallow copy
      otemp[ip]->setPsi(olist[ip]);
      psiPool->addFactory(otemp[ip]);
    }

    //find the HamiltonianFactory* to be cloned
    HamiltonianFactory* hFac=0;
    PoolType::iterator hit(myPool.find(h.getName()));
    if(hit == myPool.end()) {
      hFac=(*(myPool.begin())).second;
    } else {
      hFac=(*hit).second;
    }
    hFac->setCloneSize(np);

    vector<HamiltonianFactory*> htemp;
    htemp.resize(np,0);

//#pragma omp parallel 
//    {
//      int ip=omp_get_thread_num();
//#pragma omp critical 
//      {
//        if(ip) {
//          char hname[16];
//          sprintf(hname,"%s.c%i",h.getName().c_str(),ip);
//          htemp[ip]= hFac->clone(plist[ip],olist[ip],ip,hname);
//        }
//      }
//    }
//.........这里部分代码省略.........

void HFDBHE_smoothed_phy::add2Hamiltonian(ParticleSet& qp, TrialWaveFunction& psi, QMCHamiltonian& targetH)
{
  HFDBHE_smoothed_phy* myClone = new HFDBHE_smoothed_phy(qp);
  targetH.addOperator(myClone, myName, true);
  myClone->addCorrection(targetH);
}

void HFDBHE_smoothed_phy::addCorrection(QMCHamiltonian& targetH)
{
  HFDBHE_smoothed_aux* auxTerm = new HFDBHE_smoothed_aux(this);
  string auxName = myName+"_aux";
  targetH.addOperator(auxTerm, auxName, false);
}

  bool TrialDMCCorrection::putSpecial(xmlNodePtr cur, QMCHamiltonian& h, ParticleSet& P ) 
  {
    FirstHamiltonian = h.startIndex();
    nObservables=0;
    nValues=0;
    resum=100000;
    int blockSeries(0);
    int blockFreq(0);
    
    OhmmsAttributeSet attrib;
    attrib.add(resum,"resum");
    attrib.add(blockSeries,"max");
    attrib.add(blockFreq,"frequency");
    attrib.put(cur);
    //       app_log()<<"  Forward walking block size is "<< blockT<<"*Tau"<<endl;
    //       P.phLength=0;
    bool FIRST=true;
    CountIndex = P.addPropertyHistory(1);
    P.PropertyHistory[CountIndex][0]=0;
    xmlNodePtr tcur = cur->children;
    while(tcur != NULL) {
      string cname((const char*)tcur->name);
      //         app_log()<<cname<<endl;
      if(cname == "Observable") 
      {
        string tagName("none");
        int Hindex(-100);
//         int blockSeries(0);
//         int blockFreq(0);
        OhmmsAttributeSet Tattrib;
        Tattrib.add(tagName,"name");
//         Tattrib.add(blockSeries,"max");
//         Tattrib.add(blockFreq,"frequency");
        Tattrib.put(tcur);

        int numProps = P.PropertyList.Names.size();
        // 	  Hindex = P.PropertyList.add(tagName);
        Hindex = h.getObservable(tagName)+NUMPROPERTIES;
        if(tagName=="LocalPotential") {
          Hindex=LOCALPOTENTIAL ;
          tagName="LocPot";
        } else if (Hindex==(NUMPROPERTIES-1)){
          app_log()<<"Not a valid H element("<<Hindex<<") Valid names are:";
          for (int jk=0;jk<h.sizeOfObservables();jk++) app_log()<<"  "<<h.getObservableName(jk);
          app_log()<<endl;
          exit(-1);
        }
        /*
           if ((Hindex==-100)){
           app_log()<<" Hamiltonian Element "<<tagName<<" does not exist!! "<<Hindex<<endl;
           assert(Hindex>=0);
           }*/
        Names.push_back(tagName);
        Hindices.push_back( Hindex);
        app_log()<<" Hamiltonian Element "<<tagName<<" was found at "<< Hindex<<endl;

        int numT=blockSeries/blockFreq ;
        nObservables+=1;
        nValues+=numT;

        app_log()<<"   "<<numT<<" values will be calculated every "<<blockFreq<<"*tau H^-1"<<endl;
        vector<int> pms(3);
        pms[0]=blockFreq;
        pms[1]=numT;
        pms[2]=blockSeries+2;
        walkerLengths.push_back(pms);

        int maxWsize=blockSeries+2;
        int pindx = P.addPropertyHistory(maxWsize);
        // summed values.
        P.addPropertyHistory(numT);
        // number of times accumulated. For resum
        
        Pindices.push_back(pindx);
//         app_log()<<"pindex "<<pindx<<endl;

      }
      tcur = tcur->next;
    }
    app_log()<<"Total number of observables calculated:"<<nObservables<<endl;
    app_log()<<"Total number of values calculated:"<<nValues<<endl;
    Values.resize(nValues,0.0);
    EValues.resize(nValues,0.0);
    FWValues.resize(nValues,0.0);
    return true;
  }

  bool ForwardWalking::putSpecial(xmlNodePtr cur, QMCHamiltonian& h, ParticleSet& P)
  {
    FirstHamiltonian = h.startIndex();
    nObservables=0;
    nValues=0;
    blockT=1;
    //       OhmmsAttributeSet attrib;
    //       attrib.add(blockT,"blockSize");
    //       attrib.put(cur);
    xmlNodePtr tcur = cur->children;
    while(tcur != NULL) {
      string cname((const char*)tcur->name);
      if(cname == "Observable") 
      {
        string tagName("none");
        int Hindex(-100);
        int blockSeries(0);
        int blockFreq(0);
        OhmmsAttributeSet Tattrib;
        Tattrib.add(tagName,"name");
        Tattrib.add(blockSeries,"max");
        Tattrib.add(blockFreq,"frequency");
        Tattrib.put(tcur);

        if (tagName.find("*")==string::npos)
        {
          //Single Observable case
          int numProps = P.PropertyList.Names.size();
          Hindex = h.getObservable(tagName)+NUMPROPERTIES;
          if(tagName=="LocalPotential") {
            Hindex=LOCALPOTENTIAL ;
            tagName="LocPot";
          }
          else if(tagName=="LocalEnergy") {
            Hindex=LOCALENERGY ;
            tagName="LocEn";
          }	  
          else if (Hindex==(NUMPROPERTIES-1)){
            app_log()<<"Not a valid H element("<<Hindex<<") Valid names are:";
            for (int jk=0;jk<h.sizeOfObservables();jk++) app_log()<<"  "<<h.getObservableName(jk);
            app_log()<<endl;
            exit(-1);
          }

          Names.push_back(tagName);
          Hindices.push_back( Hindex);
          app_log()<<" Hamiltonian Element "<<tagName<<" was found at "<< Hindex<<endl;

          int numT=blockSeries/blockFreq ;
          nObservables+=1;
          nValues+=numT;

          app_log()<<"   "<<numT<<" values will be calculated every "<<blockFreq<<"*tau H^-1"<<endl;
          vector<int> pms(3);
          pms[0]=blockFreq;
          pms[1]=numT;
          pms[2]=blockSeries+2;
          walkerLengths.push_back(pms);

          int maxWsize=blockSeries+2;
          int pindx = P.addPropertyHistory(maxWsize);
          // summed values.
          //         P.addPropertyHistory(numT);

          Pindices.push_back(pindx);

        }
        else
        {
          bool FOUNDH(false);
          // 	    Multiple observables for this tag
          int found=tagName.rfind("*");
          tagName.replace (found,1,"");

          int numProps = P.PropertyList.Names.size();
          for(int j=0;j<h.sizeOfObservables();j++)
          {
            string Hname = h.getObservableName(j);
            if (Hname.find(tagName) != string::npos)
            {
              //               vector<int> Parameters;
              //               if(blockSeries==0) 
              //                 putContent(Parameters,tcur);
              //               else
              //                 for( int pl=blockFreq;pl<=blockSeries;pl+=blockFreq) Parameters.push_back(pl);
              FOUNDH=true;
              app_log()<<" Hamiltonian Element "<<Hname<<" was found at "<< j<<endl;
              Names.push_back(Hname);
              Hindex = j+NUMPROPERTIES;
              Hindices.push_back( Hindex);
              int numT=blockSeries/blockFreq ;
              nObservables+=1;
              nValues+=numT;

              app_log()<<"   "<<numT<<" values will be calculated every "<<blockFreq<<"*tau H^-1"<<endl;
              vector<int> pms(3);
              pms[0]=blockFreq;
              pms[1]=numT;
              pms[2]=blockSeries+2;
              walkerLengths.push_back(pms);
//.........这里部分代码省略.........

  /** This should be moved to branch engine */
  bool EstimatorManager::put(MCWalkerConfiguration& W, QMCHamiltonian& H, xmlNodePtr cur) 
  {
    vector<string> extra;
    cur = cur->children;
    while(cur != NULL) 
    {
      string cname((const char*)(cur->name));
      if(cname == "estimator") 
      {
        string est_name(MainEstimatorName);
        string use_hdf5("yes");
        OhmmsAttributeSet hAttrib;
        hAttrib.add(est_name, "name");
        hAttrib.add(use_hdf5, "hdf5");
        hAttrib.put(cur);
        if( (est_name == MainEstimatorName) || (est_name=="elocal") )
        {
          if(use_hdf5 == "yes")
          {
            max4ascii=H.size()+3;//write only physical energies
            add(new LocalEnergyEstimatorHDF(H),MainEstimatorName);
          }
          else
          {//fall back to the ascii file
            max4ascii=H.sizeOfObservables()+3;
            add(new LocalEnergyEstimator(H),MainEstimatorName);
          }
        }
        else if (est_name=="WFMConly")
        {
          max4ascii=H.sizeOfObservables()+10;
          app_log() << "  Using WFMConly for the MainEstimator " << endl;
          add(new WFMCOnlyEstimator(H),MainEstimatorName);
          est_name=MainEstimatorName;
        }
        else if (est_name=="releasednode")
        { 
          int Smax(100);
          int primary(1);
          OhmmsAttributeSet hAttrib;
          hAttrib.add(Smax, "Smax");
          hAttrib.add(primary, "primary");
          hAttrib.put(cur);
        
          max4ascii=H.sizeOfObservables()+ 4 + 3*(Smax+1);
          app_log() << "  Using ReleasedNode for the MainEstimator with Smax="<<Smax<<" and max4ascii="<<max4ascii << endl;
          if (primary==2) add(new ReleasedNodeEnergyEstimator(H,Smax),MainEstimatorName);
          else add(new AlternateReleasedNodeEnergyEstimator(H,Smax),MainEstimatorName);
          est_name=MainEstimatorName;
        }
        else if (est_name=="forwardwalking")
        {
          max4ascii=2*H.sizeOfObservables()+4;
          app_log() << "  Doing forwardwalking on hdf5 " << endl;
          add(new ForwardWalkingEstimator(H),MainEstimatorName);
          est_name=MainEstimatorName;
        }        else 
          extra.push_back(est_name);
      } 
      cur = cur->next;
    }

    if(Estimators.empty()) 
    {
      app_log() << "  Adding a default LocalEnergyEstimator for the MainEstimator " << endl;
      max4ascii=H.sizeOfObservables()+3;
      add(new LocalEnergyEstimator(H),MainEstimatorName);
      //add(new LocalEnergyOnlyEstimator(),MainEstimatorName);
    } 

    if(Collectables == 0 && H.sizeOfCollectables())
    {
      app_log() << "  Using CollectablesEstimator for collectables, e.g. sk, gofr, density " << endl;
      Collectables=new CollectablesEstimator(H);
      add(Collectables,"collectables");
    }
    return true;
  }