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

void
initialize(TrialWaveFunction& Psi, MCWalkerConfiguration& el,
           ParticleBase& ion) {

    enum {Up = 0, Down};

    ///add a distance table for ee IndexTypeeractions
    IndexType iee = DistanceTable::add(el,"ee");
    DistanceTableData* d_ee = DistanceTable::getTable(iee);
    ///add a distance table for ie interactions
    int iei = DistanceTable::add(ion,el,"ie");
    DistanceTableData* d_ei = DistanceTable::getTable(iei);

    ///create molecular orbital
    MolecularOrbitals<SlaterTypeOrbitals_t> *MO =
        new MolecularOrbitals<SlaterTypeOrbitals_t>;

    AtomicOrbitalBuilder<SlaterTypeOrbitals_t> orbitalbuilder;
    for(int iat=0; iat<ion.getTotalNum(); iat++) {
        orbitalbuilder.add("He",*MO);
    }
    MO->setTable(d_ei);

    typedef SlaterDeterminant<MolecularOrbitals<SlaterTypeOrbitals_t> > Det_t;
    Det_t *DetU = new Det_t(*MO,el.first(Up));
    DetU->set(el.first(Up),el.last(Up)-el.first(Up));

    Det_t* DetD = new Det_t(*MO,el.first(Down));
    DetD->set(el.first(Down),el.last(Down)-el.first(Down));

    LinearSlaterDeterminant<MolecularOrbitals<SlaterTypeOrbitals_t> >
    *asymmpsi
        = new LinearSlaterDeterminant<MolecularOrbitals<SlaterTypeOrbitals_t> >;

    asymmpsi->add(DetU);
    asymmpsi->add(DetD);

    Psi.add(asymmpsi,d_ei);

    /*
    OneBodyJastrow<PadeJastrow<ValueType>,FastWalkerIndex> *Jie
      = new OneBodyJastrow<PadeJastrow<ValueType>,FastWalkerIndex>;
    Jie->F.push_back(PadeJastrow<ValueType>(1.0,1.0));
    Psi.add(Jie,d_ei);
    */

    TwoBodyJastrow<PadeJastrow<ValueType>,FastWalkerIndex> *Jee
        = new TwoBodyJastrow<PadeJastrow<ValueType>,FastWalkerIndex>;

    Jee->F.push_back(PadeJastrow<double>(0.5,1.0));
    Jee->F.push_back(PadeJastrow<double>(1.0/6.0,1.0));
    Jee->F.push_back(PadeJastrow<double>(1.0/6.0,1.0));
    Jee->F.push_back(PadeJastrow<double>(0.5,1.0));
    Psi.add(Jee,d_ee);

}

//.........这里部分代码省略.........
    for(int iat=0; iat<Centers; iat++) {
        double sep=0.8*Cut[iat];
        for(int iel=0; iel<Core[iat]/2; iel++,ncore++) {
            el.R[ncore]=ion.R[iat]+sep*chi[ncore];
            el.R[ncore+ihalf]=ion.R[iat]+sep*chi[ncore+ihalf];
        }
    }

    int ipart = ncore;
    int isave_iat=0;

    for(int iat=0; iat<Centers; iat++) {
        for(int nn=d_ii->M[iat]; nn<d_ii->M[iat+1]; nn++) {
            double bondlength = d_ii->r(nn);
            int jat = d_ii->J[nn];
            //only assign if the half bond-length < cutoff
            if(bondlength < Cut[iat]+Cut[jat]) {
                if(ipart < ihalf) {
                    XMLReport("Assigning particles = " << ipart << " and " << ipart+ihalf)
                    /*place 2 electrons (an up and a down) at half
                      the bond-length plus a random number multiplied
                      by 10% of the bond-length*/
                    el.R[ipart] = ion.R[iat]+0.5*d_ii->dr(nn)+0.1*bondlength*chi[ipart];
                    el.R[ipart+ihalf] = ion.R[iat]+0.5*d_ii->dr(nn)+0.1*bondlength*chi[ipart+ihalf];
                    ipart++;
                    isave_iat = iat;
                }
            }
        }
    }

    //assign the last particle (if odd number of particles)
    int flag = 1;
    ipart = el.getTotalNum()-1;
    if(irem) {
        XMLReport("Assigning last particle.")
        for(int iat = isave_iat+1; iat<Centers; iat++) {
            for(int nn=d_ii->M[iat]; nn<d_ii->M[iat+1]; nn++) {
                double bondlength = d_ii->r(nn);
                if((0.5*bondlength < Cut[iat]) && flag) {
                    XMLReport("Assigning particle = " << ipart)
                    el.R[ipart] = ion.R[iat]+0.5*d_ii->dr(nn)+0.1*bondlength*chi[ipart];
                    flag = 0;
                }
            }
        }
    }

    cout << "Ionic configuration : " << ion.getName() << endl;
    ion.get(cout);

    cout << "Electronic configuration : " << el.getName() << endl;
    el.get(cout);

    string newxml(myProject.CurrentRoot());
    newxml.append(".ptcl.xml");
    ofstream ptcl_out(newxml.c_str());

    /*
      ofstream molmol("assign.xyz");

      molmol << Centers+el.getTotalNum() << endl;
      molmol << endl;

      for(int iat=0; iat<Centers; iat++)
      molmol << ion.Species.speciesName[ion.GroupID[iat]] << 0.5292*ion.R[iat] << endl;

      for(int ipart=0; ipart<el.getTotalNum(); ipart++)
      molmol << "He" << 0.5292*el.R[ipart] << endl;

      molmol.close();
    */

    xmlXPathFreeContext(m_context);
    xmlFreeDoc(m_doc);

    int nup = el.last(0);
    int ndown = el.last(1)-el.last(0);

    ptcl_out << "<?xml version=\"1.0\"?>" << endl;
    ptcl_out << "<particleset name=\"e\">" << endl;
    ptcl_out << "<group name=\"u\" size=\"" << nup << "\">" << endl;
    ptcl_out << "<parameter name=\"charge\">-1</parameter>" << endl;
    ptcl_out << "<attrib name=\"position\" datatype=\"posArray\">" << endl;
    for (int ipart=0; ipart<nup; ++ipart)
        ptcl_out << el.R[ipart] << endl;
    ptcl_out << "</attrib>" << endl;
    ptcl_out << "</group>" << endl;
    ptcl_out << "<group name=\"d\" size=\"" << ndown << "\">" << endl;
    ptcl_out << "<parameter name=\"charge\">-1</parameter>" << endl;
    ptcl_out << "<attrib name=\"position\" datatype=\"posArray\">" << endl;
    for (int ipart=nup; ipart<el.getTotalNum(); ++ipart)
        ptcl_out << el.R[ipart] << endl;
    ptcl_out << "</attrib>" << endl;
    ptcl_out << "</group>" << endl;
    ptcl_out << "</particleset>" << endl;

    OHMMS::Controller->finalize();
    return 0;
}

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