C++ Boundary::distanceSq方法代码示例

   /**
   * Generate random molecules
   */
   void Linear::generateMolecules(int nMolecule, 
      DArray<double> exclusionRadius, System& system,
      BondPotential *bondPotentialPtr, const Boundary &boundary)
   {
      int iMol;

      // Set up a cell list with twice the maxium exclusion radius as the
      // cell size
      double maxExclusionRadius = 0.0;
      for (int iType = 0; iType < system.simulation().nAtomType(); iType++) {
         if (exclusionRadius[iType] > maxExclusionRadius)
            maxExclusionRadius = exclusionRadius[iType];
      }

      // the minimum cell size is twice the maxExclusionRadius,
      // but to save memory, we take 2 times that value 
      CellList cellList;
      cellList.allocate(system.simulation().atomCapacity(),
         boundary, 2.0*2.0*maxExclusionRadius);

      if (nMolecule > capacity())
         UTIL_THROW("nMolecule > Species.capacity()!"); 

      Simulation& sim = system.simulation();
      for (iMol = 0; iMol < nMolecule; ++iMol) {
         // Add a new molecule to the system
         Molecule &newMolecule= sim.getMolecule(id());
         system.addMolecule(newMolecule);

         // Try placing atoms
         bool moleculeHasBeenPlaced = false;
         for (int iAttempt = 0; iAttempt< maxPlacementAttempts_; iAttempt++) {
            // Place first atom
            Vector pos;
            system.boundary().randomPosition(system.simulation().random(),pos);
            Atom &thisAtom = newMolecule.atom(0);
 
            // check if the first atom can be placed at the new position
            CellList::NeighborArray neighbors;
            cellList.getNeighbors(pos, neighbors);
            int nNeighbor = neighbors.size();
            bool canBePlaced = true;
            for (int j = 0; j < nNeighbor; ++j) {
               Atom *jAtomPtr = neighbors[j];
         
               double r = sqrt(system.boundary().distanceSq(
                                        jAtomPtr->position(), pos));
               if (r < (exclusionRadius[thisAtom.typeId()] +
                  exclusionRadius[jAtomPtr->typeId()])) {
                  canBePlaced = false;
                  break;
               }
            } 
            if (canBePlaced)  {
               thisAtom.position() = pos;
               cellList.addAtom(thisAtom);

               // Try to recursively place other atoms
               if (tryPlaceAtom(newMolecule, 0, exclusionRadius, system,
                  cellList, bondPotentialPtr, system.boundary())) {
                  moleculeHasBeenPlaced = true;
                  break;
              } else {
                 cellList.deleteAtom(thisAtom); 
              }
            }
         }
         if (! moleculeHasBeenPlaced) {
           std::ostringstream oss;
           oss <<  "Failed to place molecule " << newMolecule.id();
           UTIL_THROW(oss.str().c_str());
         }

      }

      #if 0
      // Check
      for (int iMol =0; iMol < nMolecule; ++iMol) {
         Molecule::AtomIterator atomIter;
         system.molecule(id(),iMol).begin(atomIter);
         for (; atomIter.notEnd(); ++atomIter) {
            for (int jMol =0; jMol < nMolecule; ++jMol) {
               Molecule::AtomIterator atomIter2;
               system.molecule(id(),jMol).begin(atomIter2);
               for (; atomIter2.notEnd(); ++atomIter2 ) {
                  if (atomIter2->id() != atomIter->id()) {
                     double r = sqrt(boundary.distanceSq(
                        atomIter->position(),atomIter2->position()));
                     if (r < (exclusionRadius[atomIter->typeId()]+
                        exclusionRadius[atomIter2->typeId()])) {
                        std::cout << r << std::endl;
                        UTIL_THROW("ERROR");
                     }
                  }
               }
            }
         }
//.........这里部分代码省略.........

   /// Add particle pairs to RDF histogram.
   void RDF::sample(long iStep) 
   {
      if (isAtInterval(iStep))  {

         accumulator_.beginSnapshot();

         System::ConstMoleculeIterator molIter1, molIter2;
         Molecule::ConstAtomIterator   atomIter1, atomIter2;
         Vector    r1, r2;
         double    dRsq, dR;
         Boundary* boundaryPtr;
         int       iSpecies1, iSpecies2, nSpecies, i;
   
         for (i = 0; i < nAtomType_; ++i) {
            typeNumbers_[i] = 0;
         }

         boundaryPtr = &system().boundary();
         nSpecies    = system().simulation().nSpecies();

         // Loop over atom 1
         for (iSpecies1 = 0; iSpecies1 < nSpecies; ++iSpecies1) {
            system().begin(iSpecies1, molIter1); 
            for ( ; molIter1.notEnd(); ++molIter1) {
               molIter1->begin(atomIter1); 
               for ( ; atomIter1.notEnd(); ++atomIter1) {
                  r1 = atomIter1->position();
 
                  ++typeNumbers_[atomIter1->typeId()];
  
                  // Loop over atom 2 
                  for (iSpecies2 = 0; iSpecies2 < nSpecies; ++iSpecies2) {
                     system().begin(iSpecies2, molIter2); 
                     for ( ; molIter2.notEnd(); ++molIter2) {

                        //Check if molecules are the same  
                        //if ( &(*molIter2) != &(*molIter1)) {

                           molIter2->begin(atomIter2);
                           for ( ; atomIter2.notEnd(); ++atomIter2) {

                              if (selector_.match(*atomIter1, *atomIter2)) {
                                 r2 = atomIter2->position();
   
                                 dRsq = boundaryPtr->distanceSq(r1, r2);
                                 dR   = sqrt(dRsq);

                                 accumulator_.sample(dR);

                              }
               
                           }

                        //}

                     }
                  } // for iSpecies2
   
               }
            }
         } // for iSpecies1

         // Increment normSum_
         double number = 0;
         for (i = 0; i < nAtomType_; ++i) {
            number  += typeNumbers_[i];
         }
         normSum_ += number*number/boundaryPtr->volume();

      } // if isAtInterval

   }

   /**
   * Recursive function to try to place an atom.
   */
   bool Linear::tryPlaceAtom(Molecule& molecule, int atomId,
      DArray<double> exclusionRadius, System& system, CellList &cellList,
      BondPotential *bondPotentialPtr, const Boundary &boundary)
   {
      Atom& lastAtom = molecule.atom(atomId);
      Atom& thisAtom = molecule.atom(++atomId);
      Random& random = system.simulation().random();

      bool hasBeenPlaced = false;

      for (int iAttempt = 0; iAttempt < maxPlacementAttempts_; iAttempt++)
      {
         // draw a random bond vector
         int beta = 1;
         Vector v;
         random.unitVector(v);
         v *= bondPotentialPtr->randomBondLength(&random, beta,
            calculateBondTypeId(lastAtom.indexInMolecule()));

         Vector newPos;
         newPos = lastAtom.position();
         newPos += v;
         // shift into simulation cell
         boundary.shift(newPos);

         // check if the atom can be placed at the new position
         CellList::NeighborArray neighbors;
         cellList.getNeighbors(newPos, neighbors);
         int nNeighbor = neighbors.size();
         bool canBePlaced = true;
         for (int j = 0; j < nNeighbor; ++j) {
            Atom *jAtomPtr = neighbors[j];
         
            double r = sqrt(boundary.distanceSq(
               jAtomPtr->position(), newPos));
            if (r < (exclusionRadius[thisAtom.typeId()] +
               exclusionRadius[jAtomPtr->typeId()])) {
               canBePlaced = false;
               break;
            }
         } 
         if (canBePlaced) {
            // place the particle
            thisAtom.position() = newPos;

            // add to cell list
            cellList.addAtom(thisAtom);

            // are we add the end of the chain?
            if (atomId == molecule.nAtom()-1) 
              return true;
            
            // recursion step
            if (! tryPlaceAtom(molecule, atomId, exclusionRadius, system,
               cellList, bondPotentialPtr, boundary) ) {
               // If the next monomer cannot be inserted, delete this monomer
               // again
               cellList.deleteAtom(thisAtom);
            } else {
               hasBeenPlaced = true;
               break;
            }
         } 
      }

      return hasBeenPlaced;
   }