C++ Link::atom1方法代码示例

   /// Evaluate and output x-velocity profile as a function of z.
   void LinkMSD::sample(long iStep) 
   { 
     if (isAtInterval(iStep))  {
      
       int    nLink = system().linkMaster().nLink();
       int    idLink, dt, dr, dr2, le, le1, i, t;     
       Link*  linkPtr;
               
       for (idLink=0; idLink < nLink; ++idLink){
	 linkPtr = &(system().linkMaster().link(idLink));
	 le = (linkPtr->tag())*2;
	 le1 = le+1;
	 if (t0_[le] != -1) {
	   t = iStep;
	   i = linkPtr->atom0().indexInMolecule();
	   dt = (t - t0_[le])/interval();
	   dr = i - i0_[le];
	   dr2 = dr*dr;
	   //std::cout << "Tag: " << linkPtr->tag() << " end: " << le << " dt: "  << dt << " dr :" << dr << std::endl;
           //std::cout << "i0: " << i0_[le] << " i: " << i << " t0: "  << t0_[le] << " t :" << t << std::endl;	   
	   if (dt < capacity_) {
	     accumulator_[dt].sample(double(dr2), outputFile_);
	   } 
	   /*else {
	      UTIL_THROW("link MSD accumulator capacity exceded");
	   }*/	   
	 }
	 else {
	   t0_[le] = iStep;
	   i0_[le] = linkPtr->atom0().indexInMolecule();
	 }
	 if (t0_[le1] != -1) {
	   t = iStep;
	   i = linkPtr->atom1().indexInMolecule();
	   dt = (t - t0_[le1])/interval();
	   dr = i - i0_[le1];
	   dr2 = dr*dr;
	   //std::cout << "Tag: " << linkPtr->tag() << " end: " << le1 << " dt: "  << dt << " dr :" << dr << std::endl;
           //std::cout << "i0: " << i0_[le1] << " i: " << i << " t0: "  << t0_[le1] << " t :" << t << std::endl;	   
	   if (dt < capacity_) {
	     accumulator_[dt].sample(double(dr2), outputFile_);
	   } 
	   /*else {
	      UTIL_THROW("link MSD accumulator capacity exceded");
	   }*/
	 }
	 else {
	   t0_[le1] = iStep;
	   i0_[le1] = linkPtr->atom1().indexInMolecule();
	 }
       }

     } // if isAtInterval
   }

   /**
   * Remove a Link.
   */
   void LinkMaster::removeLink(int id)
   {
      Link* linkPtr = &link(id);
      int  atom0Id;
      int  atom1Id;

      if (!link(id).isActive()) {
        UTIL_THROW("Attempt to remove a nonactive link");
      }
      atom0Id = linkPtr->atom0().id();
      atom1Id = linkPtr->atom1().id();
      // Remove link from atom0 and atom1 link sets
      if (!atomLinkSets_[atom0Id].isElement(*linkPtr)) {
        UTIL_THROW("Link is not in atomLinkSets of atom0");
      }
      if (!atomLinkSets_[atom1Id].isElement(*linkPtr)) {
        UTIL_THROW("Link is not in atomLinkSets of atom1");
      }
      atomLinkSets_[atom0Id].remove(*linkPtr);
      atomLinkSets_[atom1Id].remove(*linkPtr);

      // Notify observers of removal of this Link.
      // Notify before clearing Link so atomPtrs and typeId are available.
      LinkRemoveEvent event(linkPtr);
      Notifier<LinkRemoveEvent>::notifyObservers(event);

      // Clear the link: nullify atomPtrs, set typeId = -1, isActive = false
      linkPtr->clear();

      // Return link to reservoir 
      linkSet_.remove(*linkPtr);
      reservoir_.push(*linkPtr);

   }

   /**
   * connection to LinkResetEvent
   */ 
   void HoomdMove::update(const LinkResetEvent& event)
   {
      if (!HoomdIsInitialized_) return;

      Link *linkPtr = event.get();

      // First remove the link from the HOOMD bond table
      int linkTag = linkPtr->tag();
      boost::unordered_map<int, unsigned>::iterator it;
      it = linkBondMap_.find(linkTag);
      assert(it != linkBondMap_.end());
      unsigned int oldBondTag = it ->second;
      bondDataSPtr_->removeBond(oldBondTag);

      // Add the modified link 
      // link type ids are mapped according to
      // HOOMD bond id =  nBondType + link type id
      ::Bond bond(simulation().nBondType()+linkPtr->typeId(),
         linkPtr->atom0().id(),
         linkPtr->atom1().id());
      bondDataSPtr_->addBond(bond);

      // Update link-bond mapping
      unsigned int newBondTag = bondDataSPtr_->getLastTag();
      it->second = newBondTag;
   }

  /*
   * Modify one atom attached to a link
   */
   void LinkMaster::reSetAtom(Link& link, Atom& atom, int endId)
   {     
      int   atomId = atom.id();
      int   oldAtomId;
      
      if (endId==0){
        oldAtomId = link.atom0().id();
      }
      else {
        oldAtomId = link.atom1().id();
      }
      
      // Check preconditions    
      if (!link.isActive()) {
        UTIL_THROW("Attempt to reset an atom for a nonactive link");
      } 
      if (atomId < 0 || atomId >= atomCapacity_) {
         Log::file() << "AtomId       = " << atomId << std::endl;
         Log::file() << "atomCapacity = " << atomCapacity_ << std::endl;
         UTIL_THROW("Invalid atom id");
      }
 
      // Remove link from oldAtom link sets
      if (!atomLinkSets_[oldAtomId].isElement(link)) {
        UTIL_THROW("Link is not in atomLinkSets of atom");
      }     
      atomLinkSets_[oldAtomId].remove(link);       

      // Change the atom      
      if (endId==0){      
        link.setAtoms(atom,link.atom1());
      }
      else {
	link.setAtoms(link.atom0(),atom);
      }
      atomLinkSets_[atomId].append(link);
      
      ReSetAtomEvent event(&link, endId);
      Notifier<ReSetAtomEvent>::notifyObservers(event);     
 
   }      

   /*
   * Modify the atoms attached to a link
   */
   void LinkMaster::reSetAtoms(Link& link, Atom& atom0, Atom& atom1)
   {
      int   atom0Id = atom0.id();
      int   atom1Id = atom1.id();
      int   oldAtom0Id = link.atom0().id();
      int   oldAtom1Id = link.atom1().id();
     
      // Check preconditions    
      if (!link.isActive()) {
        UTIL_THROW("Attempt to reset atoms for a nonactive link");
      } 
      if (atom0Id < 0 || atom0Id >= atomCapacity_) {
         Log::file() << "Atom0Id       = " << atom0Id << std::endl;
         Log::file() << "atomCapacity = " << atomCapacity_ << std::endl;
         UTIL_THROW("Invalid atom0 id");
      }
      if (atom1Id < 0 || atom1Id >= atomCapacity_) {
         Log::file() << "Atom1Id       = " << atom1Id << std::endl;
         Log::file() << "atomCapacity = " << atomCapacity_ << std::endl;
         UTIL_THROW("Invalid atom1 id");
      }

      // Remove link from oldAtom0 and oldAtom1 link sets
      if (!atomLinkSets_[oldAtom0Id].isElement(link)) {
        UTIL_THROW("Link is not in atomLinkSets of atom0");
      }
      if (!atomLinkSets_[oldAtom1Id].isElement(link)) {
        UTIL_THROW("Link is not in atomLinkSets of atom1");
      }
      atomLinkSets_[oldAtom0Id].remove(link);
      atomLinkSets_[oldAtom1Id].remove(link);       

      // Change the atoms      
      link.setAtoms(atom0,atom1); 
      atomLinkSets_[atom0Id].append(link);
      atomLinkSets_[atom1Id].append(link);
      
      LinkResetEvent event(&link);
      Notifier<LinkResetEvent>::notifyObservers(event);      
 
   }

   /* 
   * Write the configuration file.
   */
   void McMdConfigIo::write(std::ostream &out)
   {
      using std::endl;

      // Write Boundary dimensions
      out << "BOUNDARY" << endl << endl;
      out << boundary() << endl;

      // Write atomic positions
      System::ConstMoleculeIterator molIter;
      Molecule::ConstAtomIterator   atomIter;
      Species* speciesPtr;
      int iSpecies, iMolecule;
      out << endl << "MOLECULES" << endl;
      for (iSpecies = 0; iSpecies < simulation().nSpecies(); ++iSpecies) {
         out << endl;
         out << "species    " << iSpecies << endl;
         out << "nMolecule  " << system().nMolecule(iSpecies) << endl;
         speciesPtr = &simulation().species(iSpecies);
         iMolecule = 0;
         system().begin(iSpecies, molIter); 
         for ( ; molIter.notEnd(); ++molIter) {
            out << endl;
            out << "molecule   " << iMolecule << endl;
            if (speciesPtr->isMutable()) {
               speciesPtr->mutator().writeMoleculeState(out, *molIter);
            }
            for (molIter->begin(atomIter); atomIter.notEnd(); ++atomIter) {
               //out << atomIter->position() << endl;
               writeAtom(out, *atomIter);
            }
            ++iMolecule;
         }
      }

      #ifdef INTER_TETHER
      { // Scope for local variables

         // Write Tethers
         Tether*   tetherPtr;
         Atom*     atomPtr;
         Molecule* molPtr;
         int       iTether, nTether, iAtom;
         out << std::endl;
         out << "TETHERS" << endl << endl;
         nTether = system().tetherMaster().nTether();
         out << Label("nTether") << nTether << std::endl;
         for (iTether = 0; iTether < nTether; ++iTether) {
            tetherPtr = &(system().tetherMaster().tether(iTether));
            atomPtr   = &tetherPtr->atom();
            molPtr    = &atomPtr->molecule();
            iAtom     = atomPtr->indexInMolecule();
            iMolecule = system().moleculeId(*molPtr);
            iSpecies  = molPtr->species().id();
            out << Int(iSpecies,5) << Int(iMolecule,9) << Int(iAtom,6) 
                << tetherPtr->anchor() << std::endl;
         }

      }
      #endif

      #ifdef MCMD_LINK
      { // Scope for local variables

         // Write Links
         Link*      linkPtr;
         Atom*      atomPtr;
         Molecule*  molPtr;
         int        iLink, nLink, iAtom;
         out << std::endl;
         out << "LINKS" << endl << endl;
         nLink = system().linkMaster().nLink();
         out << Label("nLink") << nLink << std::endl;
         for (iLink = 0; iLink < nLink; ++iLink) {
            linkPtr  = &(system().linkMaster().link(iLink));

            // Output species, molecule, atom ids for atom 0
            atomPtr   = &(linkPtr->atom0());
            molPtr    = &atomPtr->molecule();
            iAtom     = atomPtr->indexInMolecule();
            iMolecule = system().moleculeId(*molPtr);
            iSpecies  = molPtr->species().id();
            out << Int(iSpecies,8) << Int(iMolecule,8) << Int(iAtom,8);
            out << "   ";

            // Output species, molecule, atom ids for atom 1
            atomPtr   = &(linkPtr->atom1());
            molPtr    = &atomPtr->molecule();
            iAtom     = atomPtr->indexInMolecule();
            iMolecule = system().moleculeId(*molPtr);
            iSpecies  = molPtr->species().id();
            out << Int(iSpecies,8) << Int(iMolecule,8) << Int(iAtom,8);
            out << "   ";

            out << Int(linkPtr->typeId(),8) << std::endl;
         }

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

   // Create or destroy the slip-springs.
   bool SliplinkerAll::move() 
   {
     System::MoleculeIterator molIter;
     Molecule::AtomIterator   atomIter;
     Atom                    *atom0Ptr, *atom1Ptr;
     Atom*                    atomPtr;
     Molecule                *mol0Ptr, *mol1Ptr;
     Molecule*                molIPtr;
     double                   prob, dRSq, mindRSq=cutoff_*cutoff_, rnd, norm;
     int                      i, ntrials, j, nNeighbor, idLink, iAtom, id1, id0;
     int                      iAtom0, iAtom1, iMolecule0, iMolecule1, n0;
     Link*                    linkPtr;
     static const int         maxNeighbor = CellList::MaxNeighbor;
     double                   cdf[maxNeighbor], energy, sum;
     int                      idneighbors[maxNeighbor];

     
     ntrials = 2 * system().simulation().atomCapacity();
     for (i=0; i < ntrials; ++i){   
       
       //Choose to create or destroy a link with prob 0.5
       if (random().uniform(0.0, 1.0) > 0.5){         
	  // Try to create a link. 
	  incrementNAttempt();		
	  // Choose a molecule and atom at random
	  molIPtr  = &(system().randomMolecule(speciesId_));
	  iMolecule0 = system().moleculeId(*molIPtr);
	  iAtom   = random().uniformInt(0, molIPtr->nAtom());
	  atomPtr = &molIPtr->atom(iAtom);
	  id0 = atomPtr->id();
	  
	  // Get array of neighbors
	  system().pairPotential().cellList().getNeighbors(atomPtr->position(), neighbors_);
	  nNeighbor = neighbors_.size();
	    
	  n0 = 0;
	  sum = 0;
	  // Loop over neighboring atoms
	  for (j = 0; j < nNeighbor; ++j) {
	    atom1Ptr = neighbors_[j];
	    mol1Ptr = &atom1Ptr->molecule();
	    iMolecule1 = system().moleculeId(*mol1Ptr);	
	    id1 = atom1Ptr->id();
	    
            // Check if atoms are the same
            if (id0 != id1){	    	    
	      // Exclude masked pairs	    
	      if (!atomPtr->mask().isMasked(*atom1Ptr)) {
		// Identify possible partners and calculate the cumulative distribution function
		dRSq = system().boundary().distanceSq(atomPtr->position(), atom1Ptr->position());
		if (dRSq <= mindRSq) {
		  energy = system().linkPotential().energy(dRSq, 0);
		  //energy = 0.5*dRSq;
		  sum = sum + boltzmann(energy);
		  cdf[n0] = sum;
		  idneighbors[n0] = j;
		  n0++;
		}    
	      }
	    }
	  }   
	    
	  // If at least 1 candidate has been found.  
	  if (n0 > 0) {
	    // Choose a partner with probability cdf[j]/cdf[n0-1]
	    j = 0;
	    rnd = random().uniform(0.0, 1.0);
	    norm = 1.0/cdf[n0-1];
	    while (rnd > cdf[j]*norm ){
	      j = j + 1;
	    }     
	    atom1Ptr = neighbors_[idneighbors[j]];	  	  
	    // Create a slip-link between the selected atoms with probability = prob
	    prob = 2.0 * (system().linkMaster().nLink() + 1.0);
	    prob = system().simulation().atomCapacity() * boltzmann(-mu_) * cdf[n0-1]/ prob ;
	    //prob = 2.0 * system().nMolecule(speciesId_) * boltzmann(-mu_) * cdf[n0-1]/ prob ;	      
	    if (system().simulation().random().uniform(0.0, 1.0) < prob) {        
	      system().linkMaster().addLink(*atomPtr, *atom1Ptr, 0); 
	      incrementNAccept();	    
	    }      
	  }    
       }
       else {
	  // Try to destroy a link
	  incrementNAttempt();
	  // Choose a link at random
	  if (system().linkMaster().nLink() > 0){
	    idLink = random().uniformInt(0, system().linkMaster().nLink());
	    // Indentify the atoms for this link.
	    linkPtr = &(system().linkMaster().link(idLink));
	    atom0Ptr = &(linkPtr->atom0());
	    mol0Ptr = &atom0Ptr->molecule();
	    iMolecule0 = system().moleculeId(*mol0Ptr);
	    iAtom0 = atom0Ptr->indexInMolecule();	    
	    atom1Ptr = &(linkPtr->atom1());
	    mol1Ptr = &atom1Ptr->molecule();
	    iMolecule1 = system().moleculeId(*mol1Ptr);
	    iAtom1 = atom1Ptr->indexInMolecule();
	    // try to delete the slip-spring	
//.........这里部分代码省略.........