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

 void Write2XYZ(MCWalkerConfiguration& W)
 {
   ofstream fout("bad.xyz");
   MCWalkerConfiguration::iterator it(W.begin());
   MCWalkerConfiguration::iterator it_end(W.end());
   int nptcls(W.getTotalNum());
   while(it != it_end) {
     fout << nptcls << endl 
       << "# E = " << (*it)->Properties(LOCALENERGY) 
       << " Wgt= " << (*it)->Weight << endl;
     for(int i=0; i<nptcls; i++)
       fout << "H " << (*it)->R[i] << endl;
     ++it;
   }
 }

bool
HDFWalkerInputCollect::read(MCWalkerConfiguration& W, int firstConf, int lastConf) {

    int myID = OHMMS::Controller->mycontext();
    hid_t mastercf = H5Gopen(fileID,"config_collection");

    char confName[128];
    char coordName[128];

#if H5_VERS_RELEASE < 4
    hssize_t offset[3];
#else
    hsize_t offset[3];
#endif
    hsize_t dimIn[3],dimTot[3];
    offset[0]=0;
    offset[1]=0;
    offset[2]=0;

    typedef MCWalkerConfiguration::PosType PosType;

    vector<PosType> pos;
    int nwRead=0;

    for(int iconf=firstConf; iconf<lastConf; iconf++) {

        sprintf(coordName,"config%04d/coord",iconf);
        hid_t dataset = H5Dopen(mastercf,coordName);
        hid_t dataspace = H5Dget_space(dataset);
        int rank = H5Sget_simple_extent_ndims(dataspace);
        int status_n = H5Sget_simple_extent_dims(dataspace, dimTot, NULL);

        if(CollectMode) {
            distribute(dimTot[0]);
        } else  {
            OffSet[myID]=0;
            OffSet[myID+1]=dimTot[0];
        }

        //get the input dimension
        dimIn[0]=OffSet[myID+1]-OffSet[myID];
        dimIn[1]=dimTot[1];
        dimIn[2]=dimTot[2];
        offset[0]=OffSet[myID];

        vector<PosType> posIn(dimIn[0]*dimIn[1]);

        hid_t memspace = H5Screate_simple(3, dimIn, NULL);
        herr_t status = H5Sselect_hyperslab(dataspace,H5S_SELECT_SET, offset,NULL,dimIn,NULL);
        status = H5Dread(dataset, H5T_NATIVE_DOUBLE, memspace, dataspace, H5P_DEFAULT, &(posIn[0][0]));

        H5Sclose(memspace);
        H5Dclose(dataset);
        H5Sclose(dataspace);

        pos.insert(pos.end(), posIn.begin(), posIn.end());
        nwRead += dimIn[0];
    }

    H5Gclose(mastercf);

    int curWalker = W.getActiveWalkers();
    int nptcl=W.getTotalNum();
    if(curWalker) {
        W.createWalkers(nwRead);
    } else {
        W.resize(nwRead,nptcl);
    }

    MCWalkerConfiguration::iterator it = W.begin()+curWalker;
    int ii=0;
    for(int iw=0; iw<nwRead; iw++) {
        //std::copy(Post_temp[iw],Post_temp[iw+1], (*it)->R.begin());
        for(int iat=0; iat < nptcl; iat++,ii++) {
            (*it)->R(iat) = pos[ii];
        }
        ++it;
    }

    return true;
}

//.........这里部分代码省略.........
    xmlNodePtr el_ptr=NULL, ion_ptr=NULL;
    for(int i=0; i<result->nodesetval->nodeNr; i++) {
        xmlNodePtr cur=result->nodesetval->nodeTab[i];
        xmlChar* aname= xmlGetProp(cur,(const xmlChar*)"name");
        if(aname) {
            char fc = aname[0];
            if(fc == 'e') {
                el_ptr=cur;
            }
            else if(fc == 'i') {
                ion_ptr=cur;
            }
        }
    }

    bool donotresize = false;
    if(init_els) {
        el.setName("e");
        XMLReport("The configuration for electrons is already determined by the wave function")
        donotresize = true;
    }
    if(el_ptr) {
        XMLParticleParser pread(el,donotresize);
        pread.put(el_ptr);
    }

    if(ion_ptr) {
        XMLParticleParser pread(ion);
        pread.put(ion_ptr);
    }

    xmlXPathFreeObject(result);

    if(!ion.getTotalNum()) {
        ion.setName("i");
        ion.create(1);
        ion.R[0] = 0.0;
    }

    //The ion-ion distance-table
    DistanceTableData* d_ii = DistanceTable::getTable(DistanceTable::add(ion));
    d_ii->create(1);
    d_ii->evaluate(ion);

    vector<double> Cut, Core;
    int Centers = ion.getTotalNum();
    //attribute id for cut
    int icut = ion.Species.addAttribute("cut");
    //store the max distance from atom
    Cut.resize(Centers);
    for(int iat=0; iat<Centers; iat++) {
        int id = ion.GroupID[iat];
        Cut[iat] = ion.Species(icut,id);
    }

    int icore = ion.Species.addAttribute("core");
    //store the max distance from atom
    Core.resize(Centers);
    for(int iat=0; iat<Centers; iat++) {
        Core[iat]=ion.Species(icore,ion.GroupID[iat]);
    }

    //3N-dimensional Gaussian
    ParticleSet::ParticlePos_t chi(el.getTotalNum());
    makeGaussRandom(chi);

  ParticleSet* ParticleSetPool::createESParticleSet(xmlNodePtr cur, const string& target)
  {
    TinyVector<int,OHMMS_DIM> tilefactor;
    Tensor<int,OHMMS_DIM> tilematrix(1,0,0,0,1,0,0,0,1);
    double lr_cut=10;
    string h5name;
    string source("i");
    string bc("p p p");

    OhmmsAttributeSet attribs;
    attribs.add(h5name, "href");
    attribs.add(tilefactor, "tile");
    attribs.add(tilematrix, "tilematrix");
    attribs.add(source, "source");
    attribs.add(bc, "bconds");
    attribs.add(lr_cut, "LR_dim_cutoff");
    attribs.put(cur);

    ParticleSet* ions=getParticleSet(source);
    if(ions==0)
    {
      ions=new MCWalkerConfiguration;
      ions->setName(source);
    }

    //set the boundary condition
    ions->Lattice.LR_dim_cutoff=lr_cut;
    std::istringstream  is(bc);
    char c;
    int idim=0;
    while(!is.eof() && idim<OHMMS_DIM) 
    { 
      if(is>>c) ions->Lattice.BoxBConds[idim++]=(c=='p');
    }
    
    //initialize ions from hdf5
    hid_t h5=-1;
    if(myComm->rank()==0)
      h5 = H5Fopen(h5name.c_str(),H5F_ACC_RDONLY,H5P_DEFAULT);
    ESHDFIonsParser ap(*ions,h5,myComm);
    ap.put(cur);
    ap.expand(tilematrix);
    if(h5>-1) H5Fclose(h5);

    //failed to initialize the ions
    if(ions->getTotalNum() == 0) return 0;

    typedef ParticleSet::SingleParticleIndex_t SingleParticleIndex_t;
    vector<SingleParticleIndex_t> grid(OHMMS_DIM,SingleParticleIndex_t(1));
    ions->Lattice.reset();
    ions->Lattice.makeGrid(grid);

    if(SimulationCell==0)
    {
      SimulationCell = new ParticleSet::ParticleLayout_t(ions->Lattice);
    }

    //create the electrons
    MCWalkerConfiguration* qp = new MCWalkerConfiguration;
    qp->setName(target);
    qp->Lattice.copy(ions->Lattice);

    //qp->Lattice.reset();
    //qp->Lattice.makeGrid(grid);

    app_log() << "  Simulation cell radius = " << qp->Lattice.SimulationCellRadius << endl;
    app_log() << "  Wigner-Seitz    radius = " << qp->Lattice.WignerSeitzRadius    << endl;
    SimulationCell->print(app_log());

    myPool[target]=qp;
    myPool[source]=ions;
    //addParticleSet(qp);
    //addParticleSet(ions);

    {//get the number of electrons per spin
      vector<int> num_spin;
      xmlNodePtr cur1=cur->children;
      while(cur1!=NULL)
      {
        string cname1((const char*)cur1->name);
        if(cname1 == OrbitalBuilderBase::sd_tag)
        {
          num_spin.clear();
          xmlNodePtr cur2=cur1->children;
          while(cur2!=NULL)
          {
            string cname2((const char*)cur2->name);
            if(cname2 == OrbitalBuilderBase::det_tag)
            {
              int n=0;
              OhmmsAttributeSet a;
              a.add(n,"size");
              a.put(cur2);
              if(num_spin.size()<2) num_spin.push_back(n);
            }
            cur2=cur2->next;
          }
        }
        cur1=cur1->next;
      }
//.........这里部分代码省略.........

  // 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;
      }
    }

int main(int argc, char **argv) {

  int nblocks = 10;
  int nsteps = 10000;
  int nup = 1;
  int ndown = 1;
  int nw = 1;
  RealType Tau = 0.01;

  int iargc = 0;
  while(iargc<argc) {
    if(!strcmp(argv[iargc],"--blocks")) {
      nblocks = atoi(argv[++iargc]);
    } else if(!strcmp(argv[iargc],"--steps")) {
      nsteps = atoi(argv[++iargc]);
    } else if(!strcmp(argv[iargc],"--walkers")) {
      nw = atoi(argv[++iargc]);
    } else if(!strcmp(argv[iargc],"--up")) {
      nup = atoi(argv[++iargc]);
    } else if(!strcmp(argv[iargc],"--down")) {
      ndown = atoi(argv[++iargc]);
    } else if(!strcmp(argv[iargc],"--tau")) {
      Tau = atof(argv[++iargc]);
    }
    iargc++;
  }

  Random.init(0,1,0);
  DOMProcessor reader(argv[1]);

  MCWalkerConfiguration el;
  initQuantumParticle(el,reader);

  ///create ions
  ParticleBase ion;
  ion.create(1);
  ion.R[0] = 0.0;

  IndexType iee = DistanceTable::add(el,"ee");
  IndexType iei = DistanceTable::add(ion,el,"ie");

  ///create a trial wave function
  TrialWaveFunction Psi;
  initTrialWaveFunction(Psi,reader);

  nup = el.last(0);
  QMCHamiltonian H;

  DistanceTableData* d_ee = DistanceTable::getTable(iee);
  DistanceTableData* d_ei = DistanceTable::getTable(iei);
  H.add(new CoulombPotential(el.getTotalNum()), d_ei);
  H.add(new HartreePotential, d_ee);
  H.add(new BareKineticEnergy, NULL);

  EstimatorManager Estimators;

  Estimators.reset("vmc");
  VMC vmc(el,Psi,H,Estimators);
  vmc.initialize(nw);
  vmc.run(nblocks,nsteps,Tau);

/*
  Estimators.reset("dmc");
  DMC dmc(el,Psi,H,Estimators);
  dmc.run(nblocks,nsteps,Tau);
*/
}