C++ ParticleSet类代码示例

  TrialWaveFunction::RealType
  TrialWaveFunction::registerData(ParticleSet& P, PooledData<RealType>& buf) {
    delta_G.resize(P.getTotalNum());
    delta_L.resize(P.getTotalNum());

    P.G = 0.0;
    P.L = 0.0;

    ValueType logpsi(0.0);
    PhaseValue=0.0;
    vector<OrbitalBase*>::iterator it(Z.begin());
    vector<OrbitalBase*>::iterator it_end(Z.end());
    for(;it!=it_end; ++it)
    {
      logpsi += (*it)->registerData(P,buf);
      PhaseValue += (*it)->PhaseValue;
    }

    LogValue=real(logpsi);

    //append current gradients and laplacians to the buffer
    NumPtcls = P.getTotalNum();
    TotalDim = PosType::Size*NumPtcls;

    buf.add(&(P.G[0][0]), &(P.G[0][0])+TotalDim);
    buf.add(&(P.L[0]), &(P.L[P.getTotalNum()]));

    return LogValue;
//     cout << "Registering gradients and laplacians " << endl;
//     for(int i=0; i<P.getLocalNum(); i++) {
//       cout << P.G[i] << " " << P.L[i] << endl;
//     }
  }

/*!\fn int DistanceTable::add(const ParticleSet& s,const ParticleSet& t, const char* aname) 
 *\param s source particle set
 *\param s target particle set
 *\param aname of a new DistanceTableData
 *\return index of the distance table with the name
 *\brief Adding AsymmetricDTD to the list, e.g., el-nuclei distance table
 */
int
DistanceTable::add(const ParticleSet& s, const ParticleSet& t,  
		   const char* aname) {

  string newname;
  if(aname) {
    newname = aname;
  } else {
    newname = s.getName();
    newname.append(t.getName());
  }

  LOGMSG("Creating a distance table with " << newname)

  map<string,int>::iterator it = TableMap.find(newname);

  ///the named pair does not exist, add a new asymmetric metrics
  if(it == TableMap.end()) {
    int n = TableList.size();
    TableList.push_back(new AsymmetricDTD<NoBConds<double,3> > (s,t));
    TableMap[newname] = n;
    VisitorID.push_back(t.tag());
    return n;
  } else {
    return (*it).second;
  }
}

        /**
         * Updates the particle set according to the motion.
         *
         * @TODO Currently assumes odometry is how FakeOdometryModule creates it,
         *       ie (delta X, delta Y, etc...
         *       Verify when motion module is pulled
         * @return the updated ParticleSet.
         */
        void MotionSystem::update(ParticleSet& particles,
                                  const messages::RobotLocation& deltaMotionInfo)
        {
            ParticleIt iter;
            for(iter = particles.begin(); iter != particles.end(); iter++)
            {
                Particle* particle = &(*iter);

                /** Should be used if odometry gives global **/
                /** Should also be TESTED extensively       **/
                // float sinh, cosh;
                // sincosf(deltaMotionInfo.h() - particle->getLocation().h(),
                //         &sinh, &cosh);
                // float changeX = cosh * deltaMotionInfo.x() + sinh * deltaMotionInfo.y();
                // float changeY = cosh * deltaMotionInfo.y() - sinh * deltaMotionInfo.x();

                float changeX = deltaMotionInfo.x();
                float changeY = deltaMotionInfo.y();
                float changeH = deltaMotionInfo.h();

                particle->shift(changeX, changeY, changeH);
                randomlyShiftParticle(particle);
            }
//            std::cout << "\n\n Updated Particles w/ Motion \n";
        }

TrialWaveFunction::RealType TrialWaveFunction::registerData(ParticleSet& P, PooledData<RealType>& buf)
{
  delta_G.resize(P.getTotalNum());
  delta_L.resize(P.getTotalNum());
  P.G = 0.0;
  P.L = 0.0;
  //save the current position
  BufferCursor=buf.current();
  ValueType logpsi(0.0);
  PhaseValue=0.0;
  vector<OrbitalBase*>::iterator it(Z.begin());
  vector<OrbitalBase*>::iterator it_end(Z.end());
  for (; it!=it_end; ++it)
  {
    logpsi += (*it)->registerData(P,buf);
    PhaseValue += (*it)->PhaseValue;
  }
  convert(logpsi,LogValue);
  //LogValue=real(logpsi);
//append current gradients and laplacians to the buffer
  NumPtcls = P.getTotalNum();
  TotalDim = PosType::Size*NumPtcls;
  buf.add(PhaseValue);
  buf.add(LogValue);
  buf.add(&(P.G[0][0]), &(P.G[0][0])+TotalDim);
  buf.add(&(P.L[0]), &(P.L[P.getTotalNum()]));
  return LogValue;
}

ParticleSet::ParticleSet(const ParticleSet& p)
  : UseBoundBox(p.UseBoundBox), UseSphereUpdate(p.UseSphereUpdate),IsGrouped(p.IsGrouped)
  , ThreadID(0), mySpecies(p.getSpeciesSet()),SK(0), ParentTag(p.tag())
{
  initBase();
  initParticleSet();
  assign(p); //obly the base is copied, assumes that other properties are not assignable
  //need explicit copy:
  Mass=p.Mass;
  Z=p.Z;
  ostringstream o;
  o<<p.getName()<<ObjectTag;
  this->setName(o.str());
  app_log() << "  Copying a particle set " << p.getName() << " to " << this->getName() << " groups=" << groups() << endl;
  PropertyList.Names=p.PropertyList.Names;
  PropertyList.Values=p.PropertyList.Values;
  PropertyHistory=p.PropertyHistory;
  Collectables=p.Collectables;
  //construct the distance tables with the same order
  //first is always for this-this paier
  for (int i=1; i<p.DistTables.size(); ++i)
    addTable(p.DistTables[i]->origin());
  if(p.SK)
  {
    R.InUnit=p.R.InUnit;
    createSK();
    SK->DoUpdate=p.SK->DoUpdate;
  }
  if (p.Sphere.size())
    resizeSphere(p.Sphere.size());
  add_p_timer(myTimers);
  myTwist=p.myTwist;
}

LocalECPotential_CUDA::LocalECPotential_CUDA
(ParticleSet& ions, ParticleSet& elns) :
  LocalECPotential(ions,elns),
  ElecRef(elns), IonRef(ions),
  SumGPU("LocalECPotential::SumGPU"),
  IGPU("LocalECPotential::IGPU")
{
  SpeciesSet &sSet = ions.getSpeciesSet();
  NumIonSpecies = sSet.getTotalNum();
  NumIons  = ions.getTotalNum();
  NumElecs = elns.getTotalNum();
#ifdef QMC_CUDA
  // Copy center positions to GPU, sorting by GroupID
  gpu::host_vector<CUDA_PRECISION> I_host(OHMMS_DIM*NumIons);
  int index=0;
  for (int cgroup=0; cgroup<NumIonSpecies; cgroup++)
  {
    IonFirst.push_back(index);
    for (int i=0; i<NumIons; i++)
    {
      if (ions.GroupID[i] == cgroup)
      {
        for (int dim=0; dim<OHMMS_DIM; dim++)
          I_host[OHMMS_DIM*index+dim] = ions.R[i][dim];
        SortedIons.push_back(ions.R[i]);
        index++;
      }
    }
    IonLast.push_back(index-1);
  }
  IGPU = I_host;
  SRSplines.resize(NumIonSpecies,0);
#endif
}

  void ThreeDimMomDist::updateDistribution(ParticleSet& PtclSet, TrialWaveFunction& Psi, 
					   IndexType NumCycles) {
    TinyVector<IndexType, 3> Indexes;
    for (IndexType cycle = 0; cycle < NumCycles; cycle++) {

      pcp->NewWalker();
      PosType dr;
      for (Indexes[0] = 0; Indexes[0] < NumPts[0]; Indexes[0]++) {
	dr[1] = 0.0;
	for (Indexes[1] = 0; Indexes[1] < NumPts[1]; Indexes[1]++) {
	  dr[2] = 0.0;
	  for (Indexes[2] = 0; Indexes[2] < NumPts[2]; Indexes[2]++) {
	    IndexType partToDisplace = (*pcp)();

	    PtclSet.makeMove(partToDisplace, dr);

	    if (Indexes[0] == 0 && Indexes[1] == 0 && Indexes[2] == 0) {
	      placeIntsInBin(Indexes, 1.0);
	    } else {
	      placeIntsInBin(Indexes, Psi.ratio(PtclSet, partToDisplace));
	    }
	    totalNumSamples++;
	    PtclSet.rejectMove(partToDisplace);

	    dr[2] += Spacing[2];
	  }      
	  dr[1] += Spacing[1];
	}
	dr[0] += Spacing[0];
      }	

    } 
  }

/*void StressPBCAA::acceptMove(int active)
{
  if(is_active)
  {
    Return_t* restrict sr_ptr=SR2[active];
    Return_t* restrict pr_ptr=SR2.data()+active;
    for(int iat=0; iat<NumCenters; ++iat, ++sr_ptr,pr_ptr+=NumCenters)
      *pr_ptr = *sr_ptr += dSR[iat];
    Value=NewValue;
  }
}
*/
void StressPBCAA::initBreakup(ParticleSet& P)
{
  //SpeciesSet& tspecies(PtclRef->getSpeciesSet());
  SpeciesSet& tspecies(P.getSpeciesSet());
  //Things that don't change with lattice are done here instead of InitBreakup()
  ChargeAttribIndx = tspecies.addAttribute("charge");
  MemberAttribIndx = tspecies.addAttribute("membersize");
  NumCenters = P.getTotalNum();
  NumSpecies = tspecies.TotalNum;
 // V_const.resize(NumCenters);
  Zat.resize(NumCenters);
  Zspec.resize(NumSpecies);
  NofSpecies.resize(NumSpecies);
  for(int spec=0; spec<NumSpecies; spec++)
  {
    Zspec[spec] = tspecies(ChargeAttribIndx,spec);
    NofSpecies[spec] = static_cast<int>(tspecies(MemberAttribIndx,spec));
  }
  SpeciesID.resize(NumCenters);
  for(int iat=0; iat<NumCenters; iat++)
  {
    SpeciesID[iat]=P.GroupID[iat];
    Zat[iat] = Zspec[P.GroupID[iat]];
  }
  AA = LRCoulombSingleton::getDerivHandler(P);
  //AA->initBreakup(*PtclRef);
  myConst=evalConsts();
  myRcut=AA->get_rc();//Basis.get_rc();
  if(rVs==0)
  {
    rVs = LRCoulombSingleton::createSpline4RbyVs(AA,myRcut,myGrid);
  }
}

MomentumEstimator::MomentumEstimator(ParticleSet& elns, TrialWaveFunction& psi)
  :M(1), refPsi(psi), kgrid(4), Lattice(elns.Lattice), norm_nofK(1), hdf5_out(false)
{
  UpdateMode.set(COLLECTABLE,1);
  psi_ratios.resize(elns.getTotalNum());
  kdotp.resize(elns.getTotalNum());
  phases.resize(elns.getTotalNum());
  twist=elns.getTwist();
}

  NonLocalECPComponent::RealType 
  NonLocalECPComponent::evaluate(ParticleSet& W, TrialWaveFunction& psi,int iat, vector<NonLocalData>& Txy) {
    RealType esum=0.0;

    //int iel=0;
    for(int nn=myTable->M[iat],iel=0; nn<myTable->M[iat+1]; nn++,iel++){

      register RealType r(myTable->r(nn));
      if(r>Rmax) continue;

      register RealType rinv(myTable->rinv(nn));
      register PosType  dr(myTable->dr(nn));

      int txyCounter=Txy.size();
      // Compute ratio of wave functions
      for (int j=0; j < nknot ; j++){ 
        PosType deltar(r*rrotsgrid_m[j]-dr);
        PosType newpos(W.makeMove(iel,deltar)); 
        psiratio[j]=psi.ratio(W,iel)*sgridweight_m[j];
        W.rejectMove(iel);
        //psi.rejectMove(iel);
        //first, add a new NonLocalData with ratio
        Txy.push_back(NonLocalData(iel,psiratio[j],deltar));
      }
      // Compute radial potential
      for(int ip=0;ip< nchannel; ip++){
        vrad[ip]=nlpp_m[ip]->splint(r)*wgt_angpp_m[ip];
      }

      // Compute spherical harmonics on grid
      for (int j=0, jl=0; j<nknot ; j++){ 
        RealType zz=dot(dr,rrotsgrid_m[j])*rinv;
        // Forming the Legendre polynomials
        lpol[0]=1.0;
        RealType lpolprev=0.0;
        for (int l=0 ; l< lmax ; l++){
          //Not a big difference
          //lpol[l+1]=(2*l+1)*zz*lpol[l]-l*lpolprev;
          //lpol[l+1]/=(l+1);
          lpol[l+1]=Lfactor1[l]*zz*lpol[l]-l*lpolprev; 
          lpol[l+1]*=Lfactor2[l]; 
          lpolprev=lpol[l];
        }

        //for(int l=0; l <nchannel; l++,jl++) Amat[jl]=lpol[ angpp_m[l] ]; 
        RealType lsum=0;
        for(int l=0; l <nchannel; l++) lsum += vrad[l]*lpol[ angpp_m[l] ]; 
        esum += Txy[txyCounter++].Weight *= lsum;
      } 
     //BLAS::gemv(nknot, nchannel, &Amat[0], &psiratio[0], &wvec[0]);
     //esum += BLAS::dot(nchannel, &vrad[0], &wvec[0]);
    }   /* end loop over electron */
    return esum;

  }

 LocalECPotential::LocalECPotential(const ParticleSet& ions, ParticleSet& els):
   IonConfig(ions)
 { 
   NumIons=ions.getTotalNum();
   myTableIndex=els.addTable(ions);
   //allocate null
   PPset.resize(ions.getSpeciesSet().getTotalNum(),0);
   PP.resize(NumIons,0);
   Zeff.resize(NumIons,0.0);
   gZeff.resize(ions.getSpeciesSet().getTotalNum(),0);
 } 

RNDiracDeterminantBase::ValueType RNDiracDeterminantBase::alternateRatio(ParticleSet& P)
{
  //returns psi_T/psi_G
  for (int i=0, iat=FirstIndex; i<NumPtcls; i++, iat++)
  {
    P.G(iat) += myG_alternate(iat) - myG(iat);
    P.L(iat) += myL_alternate(iat) - myL(iat);
  }
  RealType sgn = std::cos(alternatePhaseValue);
  return sgn*std::exp(alternateLogValue-LogValue);
}

CoulombPBCAB_CUDA::CoulombPBCAB_CUDA
(ParticleSet& ions, ParticleSet& elns, bool cloning) :
  CoulombPBCAB(ions,elns,cloning),
  ElecRef(elns), IonRef(ions),
  SumGPU("CoulombPBCABTemp::SumGPU"),
  IGPU("CoulombPBCABTemp::IGPU"),
  L("CoulombPBCABTemp::L"),
  Linv("CoulombPBCABTemp::Linv"),
  kpointsGPU("CoulombPBCABTemp::kpointsGPU"),
  kshellGPU("CoulombPBCABTemp::kshellGPU"),
  FkGPU("CoulombPBCABTemp::FkGPU"),
  RhoklistGPU("CoulombPBCABTemp::RhoklistGPU"),
  RhokElecGPU("CoulombPBCABTemp::RhokElecGPU")
{
  MaxGridPoints = 8191;
  SpeciesSet &sSet = ions.getSpeciesSet();
  NumIonSpecies = sSet.getTotalNum();
  NumIons  = ions.getTotalNum();
  NumElecs = elns.getTotalNum();
#ifdef QMC_CUDA
  gpu::host_vector<CUDA_PRECISION> LHost(9), LinvHost(9);
  for (int i=0; i<3; i++)
    for (int j=0; j<3; j++)
    {
      LHost[3*i+j]    = elns.Lattice.a(j)[i];
      LinvHost[3*i+j] = elns.Lattice.b(i)[j];
    }
  L    = LHost;
  Linv = LinvHost;
  // Copy center positions to GPU, sorting by GroupID
  gpu::host_vector<CUDA_PRECISION> I_host(OHMMS_DIM*NumIons);
  int index=0;
  for (int cgroup=0; cgroup<NumIonSpecies; cgroup++)
  {
    IonFirst.push_back(index);
    for (int i=0; i<NumIons; i++)
    {
      if (ions.GroupID[i] == cgroup)
      {
        for (int dim=0; dim<OHMMS_DIM; dim++)
          I_host[OHMMS_DIM*index+dim] = ions.R[i][dim];
        SortedIons.push_back(ions.R[i]);
        index++;
      }
    }
    IonLast.push_back(index-1);
  }
  IGPU = I_host;
  SRSplines.resize(NumIonSpecies,0);
  setupLongRangeGPU();
#endif
}

StressPBCAA::StressPBCAA(ParticleSet& ref, bool active) :
  AA(0), myGrid(0), rVs(0), FirstTime(true), myConst(0.0), ForceBase(ref,ref), Ps(ref), is_active(active)
{
  ReportEngine PRE("StressPBCAA","StressPBCAA");
  //save source tag
  SourceID=ref.tag();
  //create a distance table: just to get the table name
  DistanceTableData *d_aa = DistanceTable::add(ref);
  PtclRefName=d_aa->Name;
  initBreakup(ref);
  prefix="S_"+PtclRefName;
  app_log() << "  Maximum K shell " << AA->MaxKshell << endl;
  app_log() << "  Number of k vectors " << AA->Fk.size() << endl;
  if(!is_active)
  {
    d_aa->evaluate(ref);
    update_source(ref);
   app_log()<<"Evaluating Stress SymTensor::Long Range\n"; 
    sLR=evalLR(ref);
   app_log()<<"Short Range...\n";
    sSR=evalSR(ref);
    stress=sLR+sSR+myConst;
    
    //RealType eL(0.0), eS(0.0);
    //if (computeForces)
    //{
    //  forces = 0.0;
    //  eS=evalSRwithForces(ref);
    //  // 1.3978248322
    //  eL=evalLRwithForces(ref);
    //  // 2.130267378
    //}
    //else
    //{
    //  eL=evalLR(ref);
    //  eS=evalSR(ref);
    //}
    //NewValue=Value = eL+eS+myConst;
    //app_log() << "  Fixed Coulomb potential for " << ref.getName();
    //app_log() << "\n    e-e Madelung Const. =" << MC0
    //          << "\n    Vtot     =" << Value << endl;
    
  }
  app_log() << "  Stress SymTensor components for  " << ref.getName();
  app_log() << "\n    e-e Madelung Const. =\n" << MC0
            << "\n    Stot     =\n" << stress 
            << "\n    S_SR     =\n" << sSR   
            << "\n    S_LR     =\n" << sLR
            << "\n    S_Const  =\n" << myConst<<endl;
}

MomentumEstimator::Return_t MomentumEstimator::evaluate(ParticleSet& P)
{
    const int np=P.getTotalNum();
    nofK=0.0;
    compQ=0.0;

    //will use temp[i].r1 for the Compton profile
    const vector<DistanceTableData::TempDistType>& temp(P.DistTables[0]->Temp);
    Vector<RealType> tmpn_k(nofK);
    for (int s=0; s<M; ++s)
    {
        PosType newpos;

        for (int i=0; i<OHMMS_DIM; ++i) newpos[i]=myRNG();
        //make it cartesian
        newpos=Lattice.toCart(newpos);
        P.makeVirtualMoves(newpos); //updated: temp[i].r1=|newpos-P.R[i]|, temp[i].dr1=newpos-P.R[i]
        refPsi.get_ratios(P,psi_ratios);
//         for (int i=0; i<np; ++i) app_log()<<i<<" "<<psi_ratios[i].real()<<" "<<psi_ratios[i].imag()<<endl;
        P.rejectMove(0); //restore P.R[0] to the orginal position

        for (int ik=0; ik < kPoints.size(); ++ik)
        {
            for (int i=0; i<np; ++i) kdotp[i]=dot(kPoints[ik],temp[i].dr1_nobox);
            eval_e2iphi(np,kdotp.data(),phases.data());
            RealType nofk_here(std::real(BLAS::dot(np,phases.data(),&psi_ratios[0])));//psi_ratios.data())));
            nofK[ik]+= nofk_here;
            tmpn_k[ik]=nofk_here;
        }

        for (int iq=0; iq < compQ.size(); ++iq)
            for (int i=0; i<mappedQtonofK[iq].size(); ++i)
                compQ[iq] += tmpn_k[mappedQtonofK[iq][i]];




    }
    for (int ik=0; ik<nofK.size(); ++ik) nofK[ik] *= norm_nofK;
    for (int iq=0; iq<compQ.size(); ++iq) compQ[iq] *= mappedQnorms[iq];

    if (hdf5_out)
    {
        int j=myIndex;
        for (int ik=0; ik<nofK.size(); ++ik,++j) P.Collectables[j]+= nofK[ik];
        for (int iq=0; iq<compQ.size(); ++iq,++j) P.Collectables[j]+= compQ[iq];
    }
    return 0.0;
}