C++ OhmmsAttributeSet::put方法代码示例

  /** Parses the xml input file for parameter definitions for the wavefunction optimization.
   * @param q current xmlNode 
   * @return true if successful
   */
  bool
  QMCOptimize::put(xmlNodePtr q) {

    string vmcMove("pbyp");
    string useGPU("no");
    OhmmsAttributeSet oAttrib;
    oAttrib.add(vmcMove,"move");    
    oAttrib.add(useGPU,"gpu");
    oAttrib.put(q);

    xmlNodePtr qsave=q;
    xmlNodePtr cur=qsave->children;
    int pid=OHMMS::Controller->rank();
    while(cur != NULL) {
      string cname((const char*)(cur->name));
      if(cname == "mcwalkerset") {
        mcwalkerNodePtr.push_back(cur);
      } else if(cname == "optimizer") {
        xmlChar* att= xmlGetProp(cur,(const xmlChar*)"method");
        if(att) { optmethod = (const char*)att; }
        optNode=cur;
      } else if(cname == "optimize") {
        xmlChar* att= xmlGetProp(cur,(const xmlChar*)"method");
        if(att) { optmethod = (const char*)att; }
      }
      cur=cur->next;
    }  

    //no walkers exist, add 10
    if(W.getActiveWalkers() == 0) addWalkers(omp_get_max_threads()); 

    NumOfVMCWalkers=W.getActiveWalkers();

    //create VMC engine
    if(vmcEngine ==0) {
#if defined (QMC_CUDA)
      if (useGPU == "yes")
	vmcEngine = new VMCcuda(W,Psi,H);
      else
#endif

#if defined(ENABLE_OPENMP)
      if(omp_get_max_threads()>1)
        vmcEngine = new VMCSingleOMP(W,Psi,H,hamPool);
      else
#endif
        vmcEngine = new VMCSingle(W,Psi,H);
      vmcEngine->setUpdateMode(vmcMove[0] == 'p');
      vmcEngine->initCommunicator(myComm);
    }

    vmcEngine->setStatus(RootName,h5FileRoot,AppendRun);
    vmcEngine->process(qsave);

    if(optSolver ==0)
    {
      if(optmethod == "anneal") 
      {
        app_log() << " Annealing optimization using DampedDynamics"<<endl;
        optSolver = new DampedDynamics<RealType>;
      } 
      else
      {
        app_log() << " Conjugate-gradient optimization using CGOptimization"<<endl;
        optSolver = new CGOptimization<RealType>;
      }
      //set the stream
      optSolver->setOstream(&app_log());
    }

    if(optNode == NULL) 
      optSolver->put(qmcNode);
    else
      optSolver->put(optNode);

    bool success=true;
    if(optTarget == 0) {
#if defined (QMC_CUDA)
      if (useGPU == "yes")
	optTarget = new QMCCostFunctionCUDA(W,Psi,H,hamPool);
      else
#endif
#if defined(ENABLE_OPENMP)
      if(omp_get_max_threads()>1)
      {
        optTarget = new QMCCostFunctionOMP(W,Psi,H,hamPool);
      }
      else
#endif
        optTarget = new QMCCostFunctionSingle(W,Psi,H);

      optTarget->setStream(&app_log());
      success=optTarget->put(q);
    }
    return success;
  }

/** Parses the xml input file for parameter definitions for the wavefunction optimization.
* @param q current xmlNode
* @return true if successful
*/
bool
QMCSHLinearOptimize::put(xmlNodePtr q)
{
  string useGPU("no");
  string vmcMove("pbyp");
  OhmmsAttributeSet oAttrib;
//     oAttrib.add(useGPU,"gpu");
  oAttrib.add(vmcMove,"move");
  oAttrib.put(q);
  optNode=q;
  xmlNodePtr cur=optNode->children;
  int pid=OHMMS::Controller->rank();
//     while (cur != NULL)
//     {
//         string cname((const char*)(cur->name));
//         if (cname == "mcwalkerset")
//         {
//             mcwalkerNodePtr.push_back(cur);
//         }
//         cur=cur->next;
//     }
  //no walkers exist, add 10
  if (W.getActiveWalkers() == 0)
    addWalkers(omp_get_max_threads());
  NumOfVMCWalkers=W.getActiveWalkers();
  //create VMC engine
  if ((vmcEngine ==0)||(dmcEngine ==0))
  {
#if defined (QMC_CUDA)
    if (useGPU == "yes")
      vmcEngine = new VMCcuda(W,Psi,H,psiPool);
    else
#endif
      vmcEngine = dmcEngine = new DMCOMPOPT(W,Psi,H,hamPool,psiPool);
    dmcEngine->setUpdateMode(vmcMove[0] == 'p');
    dmcEngine->initCommunicator(myComm);
  }
//     app_log()<<RootName<<"   "<<h5FileRoot<<endl;
  dmcEngine->setStatus(RootName,h5FileRoot,AppendRun);
  dmcEngine->process(optNode);
//     dmcEngine->setBranchEngine(branchEngine);
  bool success=true;
  if (optTarget == 0)
  {
#if defined (QMC_CUDA)
    if (useGPU == "yes")
      optTarget = new QMCCostFunctionCUDA(W,Psi,H,hamPool);
    else
#endif
// #if defined(ENABLE_OPENMP)
//             if (omp_get_max_threads()>1)
//             {
//                 optTarget = new QMCCostFunctionOMP(W,Psi,H,hamPool);
//             }
//             else
// #endif
//                 optTarget = new QMCCostFunctionSingle(W,Psi,H);
      optTarget = new QMCCostFunctionOMP(W,Psi,H,hamPool);
    optTarget->setneedGrads(false);
    optTarget->setStream(&app_log());
    optTarget->setDMC();
    success=optTarget->put(optNode);
  }
  return success;
}

bool BsplineJastrowBuilder::createOneBodyJastrow(xmlNodePtr cur)
{
    ReportEngine PRE(ClassName,"createOneBodyJastrow(xmlNodePtr)");
    string j1name("J1");
    {
        OhmmsAttributeSet a;
        a.add(j1name,"name");
        a.put(cur);
    }
    int taskid=targetPsi.getGroupID();//(targetPsi.is_manager())?targetPsi.getGroupID():-1;
    OBJT* J1 =new OBJT(*sourcePtcl,targetPtcl);
    DOBJT *dJ1 = new DOBJT(*sourcePtcl, targetPtcl);
    xmlNodePtr kids = cur->xmlChildrenNode;
    // Find the number of the source species
    SpeciesSet &sSet = sourcePtcl->getSpeciesSet();
    SpeciesSet &tSet = targetPtcl.getSpeciesSet();
    int numSpecies = sSet.getTotalNum();
    bool success=false;
    bool Opt(false);
    while (kids != NULL)
    {
        std::string kidsname = (char*)kids->name;
        if (kidsname == "correlation")
        {
            RealType cusp=0.0;
            string speciesA;
            string speciesB;
            OhmmsAttributeSet rAttrib;
            rAttrib.add(speciesA,"elementType");
            rAttrib.add(speciesA,"speciesA");
            rAttrib.add(speciesB,"speciesB");
            rAttrib.add(cusp,"cusp");
            rAttrib.put(kids);
            BsplineFunctor<RealType> *functor = new BsplineFunctor<RealType>(cusp);
            functor->elementType = speciesA;
            int ig = sSet.findSpecies (speciesA);
            functor->cutoff_radius = sourcePtcl->Lattice.WignerSeitzRadius;
            int jg=-1;
            if(speciesB.size())
                jg=tSet.findSpecies(speciesB);
            if (ig < numSpecies)
            {
                //ignore
                functor->put (kids);
                if (functor->cutoff_radius < 1.0e-6)
                {
                    app_log()  << "  BsplineFunction rcut is currently zero.\n"
                               << "  Setting to Wigner-Seitz radius = "
                               << sourcePtcl->Lattice.WignerSeitzRadius << endl;
                    functor->cutoff_radius = sourcePtcl->Lattice.WignerSeitzRadius;
                    functor->reset();
                }
                J1->addFunc (ig,functor,jg);
                success = true;
                dJ1->addFunc(ig,functor,jg);
                Opt=(!functor->notOpt or Opt);
                if(qmc_common.io_node)
                {
                    char fname[128];
                    if(qmc_common.mpi_groups>1)
                    {
                        if(speciesB.size())
                            sprintf(fname,"%s.%s%s.g%03d.dat",j1name.c_str(),speciesA.c_str(),speciesB.c_str(),taskid);
                        else
                            sprintf(fname,"%s.%s.g%03d.dat",j1name.c_str(),speciesA.c_str(),taskid);
                    }
                    else
                    {
                        if(speciesB.size())
                            sprintf(fname,"%s.%s%s.dat",j1name.c_str(),speciesA.c_str(),speciesB.c_str());
                        else
                            sprintf(fname,"%s.%s.dat",j1name.c_str(),speciesA.c_str());
                    }
                    functor->setReportLevel(ReportLevel,fname);
                    functor->print();
                }
            }
        }
        kids = kids->next;
    }
    if(success)
    {
        J1->dPsi=dJ1;
        targetPsi.addOrbital(J1,"J1_bspline");
        J1->setOptimizable(Opt);
        return true;
    }
    else
    {
        PRE.warning("BsplineJastrowBuilder failed to add an One-Body Jastrow.");
        delete J1;
        delete dJ1;
        return false;
    }
}

  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;
      }
//.........这里部分代码省略.........

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);
            int maxWsize=blockSeries+2;
            int pindx = P.addPropertyHistory(maxWsize);
//.........这里部分代码省略.........

SPOSetBase*
EinsplineSetBuilder::createSPOSet(xmlNodePtr cur)
{
  //use 2 bohr as the default when truncated orbitals are used based on the extend of the ions
  BufferLayer=2.0;
  OhmmsAttributeSet attribs;
  int numOrbs = 0;
  qafm=0;
  int sortBands(1);
  string sourceName;
  string spo_prec("double");
  string truncate("no");
#if defined(QMC_CUDA)
  string useGPU="yes";
#else
  string useGPU="no";
#endif
  attribs.add (H5FileName, "href");
  attribs.add (TileFactor, "tile");
  attribs.add (sortBands,  "sort");
  attribs.add (qafm,  "afmshift");
  attribs.add (TileMatrix, "tilematrix");
  attribs.add (TwistNum,   "twistnum");
  attribs.add (givenTwist,   "twist");
  attribs.add (sourceName, "source");
  attribs.add (MeshFactor, "meshfactor");
  attribs.add (useGPU,     "gpu");
  attribs.add (spo_prec,   "precision");
  attribs.add (truncate,   "truncate");
  attribs.add (BufferLayer, "buffer");
  attribs.put (XMLRoot);
  attribs.add (numOrbs,    "size");
  attribs.add (numOrbs,    "norbs");
  attribs.put (cur);
  ///////////////////////////////////////////////
  // Read occupation information from XML file //
  ///////////////////////////////////////////////
  cur = cur->children;
  int spinSet = -1;
  vector<int> Occ_Old(0,0);
  Occ.resize(0,0);
  bool NewOcc(false);
  while (cur != NULL)
  {
    string cname((const char*)(cur->name));
    if(cname == "occupation")
    {
      string occ_mode("ground");
      occ_format="energy";
      particle_hole_pairs=0;
      OhmmsAttributeSet oAttrib;
      oAttrib.add(occ_mode,"mode");
      oAttrib.add(spinSet,"spindataset");
      oAttrib.add(occ_format,"format");
      oAttrib.add(particle_hole_pairs,"pairs");
      oAttrib.put(cur);
      if(occ_mode == "excited")
      {
        putContent(Occ,cur);
      }
      else
        if(occ_mode != "ground")
        {
          app_error() << "Only ground state occupation currently supported "
                      << "in EinsplineSetBuilder.\n";
          APP_ABORT("EinsplineSetBuilder::createSPOSet");
        }
    }
    cur = cur->next;
  }
  if (Occ != Occ_Old)
  {
    NewOcc=true;
    Occ_Old = Occ;
  }
  else
    NewOcc=false;
#if defined(QMC_CUDA)
  app_log() << "\t  QMC_CUDA=1 Overwriting the precision of the einspline storage on the host.\n";
  spo_prec="double"; //overwrite
#endif
  H5OrbSet aset(H5FileName, spinSet, numOrbs);
  std::map<H5OrbSet,SPOSetBase*,H5OrbSet>::iterator iter;
  iter = SPOSetMap.find (aset);
  if ((iter != SPOSetMap.end() ) && (!NewOcc) && (qafm==0))
  {
    qafm=0;
    app_log() << "SPOSet parameters match in EinsplineSetBuilder:  "
              << "cloning EinsplineSet object.\n";
    return iter->second->makeClone();
  }
  // The tiling can be set by a simple vector, (e.g. 2x2x2), or by a
  // full 3x3 matrix of integers.  If the tilematrix was not set in
  // the input file...
  bool matrixNotSet = true;
  for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
      matrixNotSet = matrixNotSet && (TileMatrix(i,j) == 0);
  // then set the matrix to what may have been specified in the
  // tiling vector
//.........这里部分代码省略.........

  /** Create a two-body Jatrow function with a template 
   *@param cur the current xmlNode 
   *@param dummy null pointer used to identify FN
   *
   *The template class JeeType is a functor which handles the
   *evaluation of the function value, gradients and laplacians using
   *distance tables. This is a specialized builder function for
   *spin-dependent Jastrow function,e.g., for electrons, two functions
   *are created for uu(dd) and ud(du).
   */
  template<class FN> bool JAABuilder::createJAA(xmlNodePtr cur, const string& jname) 
  {

    string corr_tag("correlation");
    int ng = targetPtcl.groups();

    int ia=0, ib=0, iab=0;
    xmlNodePtr gridPtr=NULL;
    cur = cur->children;
    const SpeciesSet& species(targetPtcl.getSpeciesSet());
    typedef TwoBodyJastrowOrbital<FN> JeeType;
    JeeType *J2 = new JeeType(targetPtcl);
    typedef DiffTwoBodyJastrowOrbital<FN> dJ2Type;
    dJ2Type *dJ2 = new dJ2Type(targetPtcl);

    RealType rc=targetPtcl.Lattice.WignerSeitzRadius;
    int pairs=0;
    while(cur != NULL) {
      string cname((const char*)(cur->name));
      if(cname == corr_tag) 
      {
        string spA("u");
        string spB("u");
        OhmmsAttributeSet rAttrib;
        rAttrib.add(spA, "speciesA"); rAttrib.add(spA, "species1");
        rAttrib.add(spB, "speciesB"); rAttrib.add(spB, "species2");
        rAttrib.put(cur);
        if(spA==targetPsi.getName()) //could have used the particle name
        {
          spA=species.speciesName[0];
          spB=species.speciesName[0];
        }
        int ia = species.findSpecies(spA);
        int ib = species.findSpecies(spB);
        if(ia==species.size() || ia == species.size())
        {
          APP_ABORT("JAABuilder::createJAA is trying to use invalid species");
        }
        string pairID=spA+spB;
        FN *j= new FN;
        j->cutoff_radius=rc;
        j->put(cur);
        J2->addFunc(pairID,ia,ib,j);
	dJ2->addFunc(pairID,ia,ib,j);
        ++pairs;
      }
      cur = cur->next;
    } // while cur

    if(pairs)
    {
      J2->dPsi=dJ2;
      string j2name="J2_"+jname;
      targetPsi.addOrbital(J2,j2name);
      return true;
    }
    else
    {//clean up and delete the twobody orbitals
      delete J2;
      return false;
    }
  }

  bool PadeJastrowBuilder::put(xmlNodePtr cur) 
  {

    ReportEngine PRE(ClassName,"put()");

    string sourceOpt=targetPtcl.getName();
    string jname="PadeJastrow";
    string spin="no";
    string id_b="jee_b";
    RealType pade_b=1.0;
    OhmmsAttributeSet pattrib;
    pattrib.add(jname,"name");
    pattrib.add(spin,"spin");
    pattrib.add(sourceOpt,"source");
    pattrib.put(cur);

    cur=cur->children;
    while(cur != NULL)
    {
      {//just to hide this
        string pname="0";
        OhmmsAttributeSet aa;
        aa.add(pname,"name");
        aa.add(id_b,"id");
        aa.put(cur);
        if(pname[0]=='B') putContent(pade_b,cur);
      }

      xmlNodePtr cur1=cur->children;
      while(cur1!= NULL)
      {
        string pname="0";
        OhmmsAttributeSet aa;
        aa.add(pname,"name");
        aa.add(id_b,"id");
        aa.put(cur1);
        if(pname[0]=='B') putContent(pade_b,cur1);
        cur1=cur1->next;
      }
      cur=cur->next;
    }

    app_log() << "PadeJastrowBuilder " << id_b << " = " << pade_b << endl;

    typedef PadeFunctor<RealType> FuncType;

    typedef TwoBodyJastrowOrbital<FuncType> JeeType;
    JeeType *J2 = new JeeType(targetPtcl);

    SpeciesSet& species(targetPtcl.getSpeciesSet());
    RealType q=species(0,species.addAttribute("charge"));

    if(spin == "no") 
    {
      RealType cusp=-0.5*q*q;
      FuncType *func=new FuncType(cusp,pade_b);
      func->setIDs("jee_cusp",id_b);//set the ID's

      J2->addFunc("pade_uu",0,0,func);

      //DerivFuncType *dfunc=new DerivFuncType(cusp,B);
      //dJ2->addFunc("pade_uu",0,0,dfunc);
      //dFuncList.push_back(dfunc);
      app_log() << "    Adding Spin-independent Pade Two-Body Jastrow Cusp " << cusp<< "\n";
    } 
    else 
    {
      //build uu functor
      RealType cusp_uu=-0.25*q*q;
      FuncType *funcUU=new FuncType(cusp_uu,pade_b);
      funcUU->setIDs("pade_uu",id_b);//set the ID's

      //build ud functor
      RealType cusp_ud=-0.5*q*q;
      FuncType *funcUD=new FuncType(cusp_ud,pade_b);
      funcUD->setIDs("pade_ud",id_b);//set the ID's

      J2->addFunc("pade_uu",0,0,funcUU);

      //DerivFuncType *dfuncUU=new DerivFuncType(cusp_uu,B);
      //DerivFuncType *dfuncUD=new DerivFuncType(cusp_ud,B);
      //dJ2->addFunc("pade_uu",0,0,dfuncUU);
      //dJ2->addFunc("pade_ud",0,1,dfuncUD);
      app_log() << "    Adding Spin-dependent Pade Two-Body Jastrow " << "\n";
      app_log() << "      parallel spin     " << cusp_uu << "\n";
      app_log() << "      antiparallel spin " << cusp_ud << "\n";
    }

    targetPsi.addOrbital(J2,"J2_pade");

    if(sourceOpt != targetPtcl.getName())
    {
      map<string,ParticleSet*>::iterator pa_it(ptclPool.find(sourceOpt));
      if(pa_it == ptclPool.end()) 
      {
        PRE.warning("PadeJastrowBuilder::put failed. "+sourceOpt+" does not exist.");
        return true;
      }
      ParticleSet& sourcePtcl= (*(*pa_it).second);

//.........这里部分代码省略.........

bool QMCDriverFactory::setQMCDriver(int curSeries, xmlNodePtr cur)
{
  string curName((const char*)cur->name);
  string update_mode("pbyp");
  string qmc_mode("invalid");
  string multi_tag("no");
  string warp_tag("no");
  string append_tag("no");
#if defined(QMC_CUDA)
  string gpu_tag("yes");
#else
  string gpu_tag("no");
#endif
  string traces_tag("none");
  OhmmsAttributeSet aAttrib;
  aAttrib.add(qmc_mode,"method");
  aAttrib.add(update_mode,"move");
  aAttrib.add(multi_tag,"multiple");
  aAttrib.add(warp_tag,"warp");
  aAttrib.add(append_tag,"append");
  aAttrib.add(gpu_tag,"gpu");
  aAttrib.add(traces_tag,"trace");
  aAttrib.put(cur);
  bitset<QMC_MODE_MAX>  WhatToDo;
  bool append_run =(append_tag == "yes");
  WhatToDo[SPACEWARP_MODE]= (warp_tag == "yes");
  WhatToDo[MULTIPLE_MODE]= (multi_tag == "yes");
  WhatToDo[UPDATE_MODE]= (update_mode == "pbyp");
#if defined(QMC_CUDA)
  WhatToDo[GPU_MODE      ] = (gpu_tag     == "yes");
#endif
  OhmmsInfo::flush();
  QMCRunType newRunType = DUMMY_RUN;
  if(curName != "qmc")
    qmc_mode=curName;
  int nchars=qmc_mode.size();
  if(qmc_mode.find("linear") < nchars)
  {
    if (qmc_mode.find("cslinear") < nchars)
      newRunType=CS_LINEAR_OPTIMIZE_RUN;
    else
      newRunType=LINEAR_OPTIMIZE_RUN;
  }
  else
    if(qmc_mode.find("opt") < nchars)
    {
      newRunType=OPTIMIZE_RUN;
    }
    else
    {
      if(qmc_mode.find("ptcl")<nchars)
        WhatToDo[UPDATE_MODE]=1;
      if(qmc_mode.find("mul")<nchars)
        WhatToDo[MULTIPLE_MODE]=1;
      if(qmc_mode.find("warp")<nchars)
        WhatToDo[SPACEWARP_MODE]=1;
//       if (qmc_mode.find("rmcPbyP")<nchars)
//       {
//         newRunType=RMC_PBYP_RUN;
//       }
//       else
      if(qmc_mode.find("rmc")<nchars)
      {
        newRunType=RMC_RUN;
      }
      else
        if(qmc_mode.find("vmc")<nchars)
        {
          newRunType=VMC_RUN;
        }
        else
          if(qmc_mode.find("dmc")<nchars)
          {
            newRunType=DMC_RUN;
          }
    }
  unsigned long newQmcMode=WhatToDo.to_ulong();
  //initialize to 0
  QMCDriver::BranchEngineType* branchEngine=0;
  if(qmcDriver)
  {
    if( newRunType != curRunType || newQmcMode != curQmcMode)
    {
      if(curRunType == DUMMY_RUN)
      {
        APP_ABORT("QMCDriverFactory::setQMCDriver\n Other qmc sections cannot come after <qmc method=\"test\">.\n");
      }
      //pass to the new driver
      branchEngine=qmcDriver->getBranchEngine();
      delete qmcDriver;
      //set to 0 so that a new driver is created
      qmcDriver = 0;
      //if the current qmc method is different from the previous one, append_run is set to false
      append_run = false;
    }
    else
    {
      app_log() << "  Reusing " << qmcDriver->getEngineName() << endl;
      //         if(curRunType == DMC_RUN)
      qmcDriver->resetComponents(cur);
//.........这里部分代码省略.........

  bool SpinDensity::put(xmlNodePtr cur)
  {
    using std::ceil;
    using std::sqrt;
    reset();
    string write_report = "no";
    OhmmsAttributeSet attrib;
    attrib.add(myName,"name");
    attrib.add(write_report,"report");
    attrib.put(cur);

    bool have_dr = false;
    bool have_grid = false;
    bool have_center = false;
    bool have_corner = false;
    bool have_cell = false;

    PosType dr;
    PosType center;
    Tensor<RealType,DIM> axes;

    int test_moves = 0;

    xmlNodePtr element = cur->xmlChildrenNode;
    while(element!=NULL)
    {
      string ename((const char*)element->name);
      if(ename=="parameter")
      {
        string name((const char*)(xmlGetProp(element,(const xmlChar*)"name")));
        if(name=="dr")      
        {
          have_dr = true;
          putContent(dr,element);   
        }
        else if(name=="grid") 
        {
          have_grid = true;
          putContent(grid,element);        
        }
        else if(name=="corner") 
        {
          have_corner = true;
          putContent(corner,element);        
        }
        else if(name=="center") 
        {
          have_center = true;
          putContent(center,element);        
        }
        else if(name=="cell") 
        {
          have_cell = true;
          putContent(axes,element);        
        }
        else if(name=="test_moves") 
          putContent(test_moves,element);        
      }
      element = element->next;
    }

    if(have_dr && have_grid)
    {
      APP_ABORT("SpinDensity::put  dr and grid are provided, this is ambiguous");
    }
    else if(!have_dr && !have_grid)
      APP_ABORT("SpinDensity::put  must provide dr or grid");

    if(have_corner && have_center)
      APP_ABORT("SpinDensity::put  corner and center are provided, this is ambiguous");
    if(have_cell)
    {
      cell.set(axes);
      if(!have_corner && !have_center)
        APP_ABORT("SpinDensity::put  must provide corner or center");
    }
    else
      cell = Ptmp->Lattice;

    if(have_center)
      corner = center-cell.Center;

    if(have_dr)
      for(int d=0;d<DIM;++d)
        grid[d] = (int)ceil(sqrt(dot(cell.Rv[d],cell.Rv[d]))/dr[d]);

    npoints = 1;
    for(int d=0;d<DIM;++d)
      npoints *= grid[d];
    gdims[0] = npoints/grid[0];
    for(int d=1;d<DIM;++d)
      gdims[d] = gdims[d-1]/grid[d];

    if(write_report=="yes")
      report("  ");
    if(test_moves>0)
      test(test_moves,*Ptmp);

    return true;
  }

  bool 
    ECPComponentBuilder::parseCasino(const std::string& fname, xmlNodePtr cur)
  {

    app_log() << "   Start ECPComponentBuilder::parseCasino" << endl;
    RealType rmax=2.0;
    Llocal=-1;
    Lmax=-1;
    OhmmsAttributeSet aAttrib;
    aAttrib.add(rmax,"cutoff");
    aAttrib.add(Llocal,"l-local");
    aAttrib.add(Lmax,"lmax");
    aAttrib.add(Nrule,"nrule");
    aAttrib.put(cur);
    //const xmlChar* rptr=xmlGetProp(cur,(const xmlChar*)"cutoff");
    //if(rptr != NULL) rmax = atof((const char*)rptr);


    //app_log() << "   Creating a Linear Grid Rmax=" << rmax << endl;
    //const RealType d=5e-4;
    //LinearGrid<RealType>* agrid = new LinearGrid<RealType>;
    //int ng=static_cast<int>(rmax/d)+1;
    //agrid->set(0,rmax,ng);

    ifstream fin(fname.c_str(),ios_base::in);
    if(!fin)
    {
      app_error() << "Could not open file " << fname << endl;
      APP_ABORT("ECPComponentBuilder::parseCasino");
    }      

    if(pp_nonloc==0) pp_nonloc=new NonLocalECPComponent; 

    OhmmsAsciiParser aParser;
    int atomNumber=0;
    int npts=0, idummy;
    string eunits("rydberg");

    app_log() << "    ECPComponentBuilder::parseCasino" <<endl;
    aParser.skiplines(fin,1);//Header
    aParser.skiplines(fin,1);//Atomic number and pseudo-charge
    aParser.getValue(fin,atomNumber,Zeff);
    app_log() << "      Atomic number = " << atomNumber << "  Zeff = " << Zeff << endl;
    aParser.skiplines(fin,1);//Energy units (rydberg/hartree/ev):
    aParser.getValue(fin,eunits);
    app_log() << "      Unit of the potentials = " << eunits << endl;
    RealType Vprefactor = (eunits == "rydberg")?0.5:1.0;

    aParser.skiplines(fin,1);//Angular momentum of local component (0=s,1=p,2=d..)
    aParser.getValue(fin,idummy);
    if(Lmax<0) Lmax=idummy;
    aParser.skiplines(fin,1);//NLRULE override (1) VMC/DMC (2) config gen (0 ==> input/default value)
    aParser.skiplines(fin,1);//0 0, not sure what to do yet
    aParser.skiplines(fin,1);//Number of grid points
    aParser.getValue(fin,npts);
    app_log() << "      Input Grid size = " << npts << endl;
    vector<RealType> temp(npts);

    aParser.skiplines(fin,1);//R(i) in atomic units
    aParser.getValues(fin,temp.begin(),temp.end());

    //create a global grid of numerical type
    grid_global= new NumericalGrid<RealType>(temp);

    Matrix<RealType> vnn(Lmax+1,npts);
    for(int l=0; l<=Lmax; l++)
    {
      aParser.skiplines(fin,1);
      aParser.getValues(fin,vnn[l],vnn[l]+npts);
    }

    vector<int> angList(Lmax+1);
    for(int l=0; l<=Lmax; l++) angList[l]=l;

    // Now, check to see what maximum cutoff should be
    if(vnn.size()>1)
    {
      const double tolerance=1.0e-5;
      double rc_check = grid_global->r(npts-1);
      for (int j=npts-1; j>0; j++) {
        bool closeEnough = true;
        for (int i=0; i<vnn.rows(); i++)
          for (int k=i+1; k<vnn.rows(); k++)
            if (std::fabs(vnn[i][j] - vnn[k][j]) > tolerance)
              closeEnough = false;
        if (!closeEnough) {
          rc_check = grid_global->r(j);
          break;
        }
      }
      app_log() << "  Maxium cutoff for non-local pseudopotentials " << rc_check << endl;
    }

    doBreakUp(angList,vnn,rmax,Vprefactor);

    SetQuadratureRule(Nrule);
    app_log() << "    Non-local pseudopotential parameters" <<endl;
    pp_nonloc->print(app_log());

    return true;
//.........这里部分代码省略.........

  bool QMCDriverFactory::setQMCDriver(int curSeries, xmlNodePtr cur) 
  {

    string curName((const char*)cur->name);
    string update_mode("walker");
    string qmc_mode("invalid");
    string multi_tag("no");
    string warp_tag("no");
    string append_tag("no");

    OhmmsAttributeSet aAttrib;
    aAttrib.add(qmc_mode,"method");
    aAttrib.add(update_mode,"move");
    aAttrib.add(multi_tag,"multiple");
    aAttrib.add(warp_tag,"warp");
    aAttrib.add(append_tag,"append");
    aAttrib.put(cur);

    bool append_run =(append_tag == "yes");
    bitset<3>  WhatToDo;
    WhatToDo[SPACEWARP_MODE]= (warp_tag == "yes");
    WhatToDo[MULTIPLE_MODE]= (multi_tag == "yes");
    WhatToDo[UPDATE_MODE]= (update_mode == "pbyp");

    QMCRunType newRunType = DUMMY_RUN;
    if(curName != "qmc") qmc_mode=curName;
    int nchars=qmc_mode.size();
    if(qmc_mode.find("opt") < nchars)
    {
      newRunType=OPTIMIZE_RUN;
    }
    else
    {
      if(qmc_mode.find("vmc")<nchars)
      {
        newRunType=VMC_RUN;
      }
      else if(qmc_mode.find("dmc")<nchars)
      {
        newRunType=DMC_RUN;
      }
      else if(qmc_mode.find("rmc")<nchars)
      {
        newRunType=RMC_RUN;
      }
      if(qmc_mode.find("ptcl")<nchars) WhatToDo[UPDATE_MODE]=1;
      if(qmc_mode.find("mul")<nchars) WhatToDo[MULTIPLE_MODE]=1;
      if(qmc_mode.find("warp")<nchars) WhatToDo[SPACEWARP_MODE]=1;
    } 

    unsigned long newQmcMode=WhatToDo.to_ulong();
   
    //initialize to 0
    QMCDriver::BranchEngineType* branchEngine=0;

    if(qmcDriver) {
      if(newRunType != curRunType || newQmcMode != curQmcMode) {
        //copy the pointer of the BranchEngine 
        branchEngine=qmcDriver->getBranchEngine();
        //remove the qmcDriver
        delete qmcDriver;
        //set to 0 so that a new driver is created
        qmcDriver = 0;
        //if the current qmc method is different from the previous one, append_run is set to false
        append_run = false;
      } else { 
        app_log() << "  Reusing " << qmcDriver->getEngineName() << endl;
      }
    }

    if(curSeries == 0) append_run = false;

    //carryon with the existing qmcDriver
    if(qmcDriver) return append_run;

    //need to create a qmcDriver
    curRunType = newRunType;
    curQmcMode = newQmcMode;
    curQmcModeBits = WhatToDo;
    createQMCDriver(cur);

    if(qmcComm)
      qmcDriver->setCommunicator(qmcComm);
    else
      qmcDriver->setCommunicator(OHMMS::Controller);

    //branchEngine has to be transferred to a new QMCDriver
    if(branchEngine) qmcDriver->setBranchEngine(branchEngine);

    return append_run;
  }

  bool JastrowBuilder::addThreeBody(xmlNodePtr cur) 
  {
    if(sourceOpt == targetPtcl.getName()) 
    {
      app_warning() << "  Three-Body Jastrow Function needs a source different from " << targetPtcl.getName() << endl;
      app_warning() << "  Exit JastrowBuilder::addOneBody." << endl;
      return false;
    }

    PtclPoolType::iterator pit(ptclPool.find(sourceOpt));
    if(pit == ptclPool.end()) {
      app_error() << "     JastrowBuilder::addThreeBody requires a center. " << sourceOpt << " is invalid " << endl;
      return false;
    }
    app_log() << "  JastrowBuilder::addThreeBody source="<< sourceOpt <<  endl;
    xmlNodePtr basisPtr=NULL;
    xmlNodePtr coeffPtr=NULL;
    cur = cur->xmlChildrenNode;
    string diagOnly("no");
    while(cur != NULL) {
      string cname((const char*)(cur->name));
      if(cname == basisset_tag) 
      {
        basisPtr=cur;
      } 
      else if(cname.find("coeff")<cname.size())
      {
        coeffPtr=cur;
        OhmmsAttributeSet oAttrib;
        oAttrib.add(diagOnly,"diagonal");
        oAttrib.put(cur);
      }
      cur=cur->next;
    }

    if(basisPtr == NULL)
    {
      app_error() << "     JastrowBuilder::addThreeBody exit. Missing <basisset/>"<< endl;
      return false;
    }

    ParticleSet* sourcePtcl=(*pit).second;
    JastrowBasisBuilder* basisBuilder = 
      new JastrowBasisBuilder(targetPtcl,*sourcePtcl,funcOpt,transformOpt == "yes");
    basisBuilder->put(basisPtr);

    if(diagOnly == "yes")
    {
      app_log() << "\n  creating Three-Body Jastrow function using only diagnoal blocks." << endl;
      ThreeBodyBlockSparse* J3 = new ThreeBodyBlockSparse(*sourcePtcl, targetPtcl);
      J3->setBasisSet(basisBuilder->myBasisSet);
      J3->put(coeffPtr,targetPsi.VarList);
      J3->setBlocks(basisBuilder->SizeOfBasisPerCenter);
      targetPsi.addOrbital(J3);
    } 
    else
    {
      app_log() << "\n  creating Three-Body Jastrow function using a complete Geminal matrix." << endl;
      ThreeBodyGeminal* J3 = new ThreeBodyGeminal(*sourcePtcl, targetPtcl);
      J3->setBasisSet(basisBuilder->myBasisSet);
      J3->put(coeffPtr,targetPsi.VarList);
      targetPsi.addOrbital(J3);
    }

    //if(jbuilder)
    //{
    //  jbuilder->put(cur);
    //  Children.push_back(jbuilder);
    //  return true;
    //}
    //    } else if (jasttype == "Three-Body-Pade") {
    //      app_log() << "\n  creating Three-Body-Pade Jastrow function " << endl;
    //      string source_name("i");
    //      const xmlChar* iptr = xmlGetProp(cur, (const xmlChar *)"source");
    //      //if(iptr != NULL) source_name=(const char*)iptr;
    //      PtclPoolType::iterator pit(ptclPool.find(source_name));
    //      if(pit != ptclPool.end()) {
    //        jbuilder = new ThreeBodyPadeBuilder(*targetPtcl,*targetPsi,*((*pit).second));
    //      }
    //    }
    return true;
  }

  bool RPAJastrow::put(xmlNodePtr cur){
    ReportEngine PRE("RPAJastrow","put");
    xmlNodePtr myNode=xmlCopyNode(cur,1);

    //capture attribute jastrow/@name
    MyName="RPA_Jee";
    string useL="yes";
    string useS="yes";
    rpafunc="breakup";
    
    OhmmsAttributeSet a;
    a.add(MyName,"name");
    a.add(useL,"longrange");
    a.add(useS,"shortrange");
    a.add(rpafunc,"function");
    a.put(cur);

    Rs=-1.0;
    Kc=-1.0;
    string ID_Rs="RPA_rs";
    ParameterSet params;
    params.add(Rs,"rs","double");
    params.add(Kc,"kc","double");
    params.put(cur);

    app_log() <<endl<<"   LongRangeForm is "<<rpafunc<<endl;

    DropLongRange = (useL == "no");
    DropShortRange = (useS=="no");

    app_log() << "    Rs can be optimized using ID=" << ID_Rs << endl;
    RealType tlen = std::pow(3.0/4.0/M_PI*targetPtcl.Lattice.Volume/ static_cast<RealType>(targetPtcl.getTotalNum()) ,1.0/3.0);
    
    if(Rs<0) {
      if(targetPtcl.Lattice.SuperCellEnum) {
        Rs=tlen;
      } else {
        cout<<"  Error finding rs. Is this an open system?!"<<endl;
        Rs=100.0;
      }
    }
    
    //Add Rs to optimizable list
    myVars.insert(ID_Rs,Rs,true);
    
    int indx = targetPtcl.SK->KLists.ksq.size()-1;
    double Kc_max=std::pow(targetPtcl.SK->KLists.ksq[indx],0.5);

    if(Kc<0){ 
      Kc = 2.0*  std::pow(2.25*M_PI,1.0/3.0)/tlen ;
    }
    
    if(Kc>Kc_max){
      Kc=Kc_max;
      app_log() << "    Kc set too high. Resetting to the maximum value"<<endl;
    }

    app_log() << "    RPAJastrowBuilder::addTwoBodyPart Rs = " << Rs <<  "  Kc= " << Kc << endl;
    
    if (rpafunc=="Yukawa" || rpafunc=="breakup"){
      myHandler= new LRHandlerTemp<YukawaBreakup<RealType>,LPQHIBasis>(targetPtcl,Kc);
    } else if (rpafunc=="RPA"){
      myHandler= new LRRPAHandlerTemp<RPABreakup<RealType>,LPQHIBasis>(targetPtcl,Kc);
    } else if (rpafunc=="dYukawa"){
      myHandler= new LRHandlerTemp<DerivYukawaBreakup<RealType>,LPQHIBasis >(targetPtcl,Kc);
    } else if (rpafunc=="dRPA"){
      myHandler= new LRRPAHandlerTemp<DerivRPABreakup<RealType>,LPQHIBasis >(targetPtcl,Kc);
    }
    
    
    myHandler->Breakup(targetPtcl,Rs);
    
    app_log() << "  Maximum K shell " << myHandler->MaxKshell << endl;
    app_log() << "  Number of k vectors " << myHandler->Fk.size() << endl;

    if(!DropLongRange) makeLongRange();
    if(!DropShortRange) makeShortRange();

    return true;
  }

  bool ECPotentialBuilder::put(xmlNodePtr cur) {

    if(localPot.empty()) {
      int ng(IonConfig.getSpeciesSet().getTotalNum());
      localZeff.resize(ng,1);
      localPot.resize(ng,0);
      nonLocalPot.resize(ng,0);
    }

    string ecpFormat("table");
    string pbc("yes");
    string forces("no");
    OhmmsAttributeSet pAttrib;
    pAttrib.add(ecpFormat,"format");
    pAttrib.add(pbc,"pbc");
    pAttrib.add(forces,"forces");
    pAttrib.put(cur);
    bool doForces = (forces == "yes") || (forces == "true");

    //const xmlChar* t=xmlGetProp(cur,(const xmlChar*)"format");
    //if(t != NULL) {
    //  ecpFormat= (const char*)t;
    //} 

    if(ecpFormat == "xml")  
    {
      useXmlFormat(cur);
    } 
    else 
    {
      useSimpleTableFormat();
    } 

    ///create LocalECPotential
    bool usePBC = 
      !(IonConfig.Lattice.SuperCellEnum == SUPERCELL_OPEN || pbc =="no");
    if(hasLocalPot) {
      if(IonConfig.Lattice.SuperCellEnum == SUPERCELL_OPEN || pbc =="no") 
      {
#ifdef QMC_CUDA
        LocalECPotential_CUDA* apot = 
	  new LocalECPotential_CUDA(IonConfig,targetPtcl);
#else
        LocalECPotential* apot = new LocalECPotential(IonConfig,targetPtcl);
#endif
        for(int i=0; i<localPot.size(); i++) {
          if(localPot[i]) apot->add(i,localPot[i],localZeff[i]);
        }
        targetH.addOperator(apot,"LocalECP");
      }
      else
      {
	if (doForces) 
	  app_log() << "  Will compute forces in CoulombPBCABTemp.\n" << endl;
#ifdef QMC_CUDA
        CoulombPBCAB_CUDA* apot=
	  new CoulombPBCAB_CUDA(IonConfig,targetPtcl, doForces);
#else
	CoulombPBCABTemp* apot =
	  new CoulombPBCABTemp(IonConfig,targetPtcl, doForces);
#endif
        for(int i=0; i<localPot.size(); i++) {
          if(localPot[i]) apot->add(i,localPot[i]);
        }
        targetH.addOperator(apot,"LocalECP");
      }
      //if(IonConfig.Lattice.BoxBConds[0]) {
      //  CoulombPBCABTemp* apot=new CoulombPBCABTemp(IonConfig,targetPtcl);
      //  for(int i=0; i<localPot.size(); i++) {
      //    if(localPot[i]) apot->add(i,localPot[i]);
      //  }
      //  targetH.addOperator(apot,"LocalECP");
      //} else {
      //  LocalECPotential* apot = new LocalECPotential(IonConfig,targetPtcl);
      //  for(int i=0; i<localPot.size(); i++) {
      //    if(localPot[i]) apot->add(i,localPot[i],localZeff[i]);
      //  }
      //  targetH.addOperator(apot,"LocalECP");
      //}
    }

    if(hasNonLocalPot) {
      //resize the sphere
      targetPtcl.resizeSphere(IonConfig.getTotalNum());
      RealType rc2=0.0;
#ifdef QMC_CUDA   
      NonLocalECPotential_CUDA* apot = 
	new NonLocalECPotential_CUDA(IonConfig,targetPtcl,targetPsi,usePBC,doForces);
#else
      NonLocalECPotential* apot = 
	new NonLocalECPotential(IonConfig,targetPtcl,targetPsi, doForces);
#endif

      for(int i=0; i<nonLocalPot.size(); i++) 
      {
        if(nonLocalPot[i]) 
        {
          rc2=std::max(rc2,nonLocalPot[i]->Rmax);
          apot->add(i,nonLocalPot[i]);
        }
//.........这里部分代码省略.........