C++ Molecule::beginRigidBody方法代码示例

  DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) {
    nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies());
    nObjects_.push_back(info_->getNGlobalBonds());
    nObjects_.push_back(info_->getNGlobalBends());
    nObjects_.push_back(info_->getNGlobalTorsions());
    nObjects_.push_back(info_->getNGlobalInversions());

    stuntdoubles_.resize(nObjects_[STUNTDOUBLE]);
    bonds_.resize(nObjects_[BOND]);
    bends_.resize(nObjects_[BEND]);
    torsions_.resize(nObjects_[TORSION]);
    inversions_.resize(nObjects_[INVERSION]);
    
    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator rbIter;
    Molecule::BondIterator bondIter;
    Molecule::BendIterator bendIter;
    Molecule::TorsionIterator torsionIter;
    Molecule::InversionIterator inversionIter;

    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;    
    
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
        
      for(atom = mol->beginAtom(ai); atom != NULL; 
          atom = mol->nextAtom(ai)) {
	stuntdoubles_[atom->getGlobalIndex()] = atom;
      }
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
           rb = mol->nextRigidBody(rbIter)) {
	stuntdoubles_[rb->getGlobalIndex()] = rb;
      }
      for (bond = mol->beginBond(bondIter); bond != NULL; 
           bond = mol->nextBond(bondIter)) {
        bonds_[bond->getGlobalIndex()] = bond;
      }   
      for (bend = mol->beginBend(bendIter); bend != NULL; 
           bend = mol->nextBend(bendIter)) {
        bends_[bend->getGlobalIndex()] = bend;
      }   
      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; 
           torsion = mol->nextTorsion(torsionIter)) {
        torsions_[torsion->getGlobalIndex()] = torsion;
      }   
      for (inversion = mol->beginInversion(inversionIter); inversion != NULL; 
           inversion = mol->nextInversion(inversionIter)) {
        inversions_[inversion->getGlobalIndex()] = inversion;
      }   

    }
  }

  void IndexFinder::init() {

    SimInfo::MoleculeIterator mi;
    Molecule::AtomIterator ai;
    Molecule::RigidBodyIterator rbIter;
    Molecule::BondIterator bondIter;
    Molecule::BendIterator bendIter;
    Molecule::TorsionIterator torsionIter;
    Molecule::InversionIterator inversionIter;

    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    Bond* bond;
    Bend* bend;
    Torsion* torsion;
    Inversion* inversion;    
    
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
           
      SelectionSet ss(nObjects_);
      ss.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());

      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
	ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
      }
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
           rb = mol->nextRigidBody(rbIter)) {
        ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
      }
      for (bond = mol->beginBond(bondIter); bond != NULL; 
           bond = mol->nextBond(bondIter)) {
        ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
      }   
      for (bend = mol->beginBend(bendIter); bend != NULL; 
           bend = mol->nextBend(bendIter)) {
        ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
      }   
      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; 
           torsion = mol->nextTorsion(torsionIter)) {
        ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
      }   
      for (inversion = mol->beginInversion(inversionIter); inversion != NULL; 
           inversion = mol->nextInversion(inversionIter)) {
        ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
      }   

      selectionSets_[mol->getGlobalIndex()] = ss;
    }
  }

  void ObjectCount::process() {
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
  
    counts_.clear();
    counts_.resize(10, 0);
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    unsigned long int nsum = 0;
    unsigned long int n2sum = 0;

    for (int i = 0; i < nFrames; i += step_) {
      reader.readFrame(i);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    
      for (mol = info_->beginMolecule(mi); mol != NULL; 
	   mol = info_->nextMolecule(mi)) {
	//change the positions of atoms which belong to the rigidbodies
	for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
	     rb = mol->nextRigidBody(rbIter)) {
	  rb->updateAtoms();        
	}      
      }
    
      if (evaluator_.isDynamic()) {
	seleMan_.setSelectionSet(evaluator_.evaluate());
      }
        
      unsigned int count = seleMan_.getSelectionCount();

      if (counts_.size() <= count)  {
	counts_.resize(count, 0);
      }

      counts_[count]++;

      nsum += count;
      n2sum += count * count;
    }
   
    int nProcessed = nFrames /step_;

    nAvg = nsum / nProcessed;
    n2Avg = n2sum / nProcessed;
    sDev = sqrt(n2Avg - nAvg*nAvg);
    writeCounts();   
  }

void P2OrderParameter::process() {
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    StuntDouble* sd1;
    StuntDouble* sd2;
    int ii;
    int jj;
    int vecCount;
    bool usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();

    DumpReader reader(info_, dumpFilename_);
    int nFrames = reader.getNFrames();

    for (int i = 0; i < nFrames; i += step_) {
        reader.readFrame(i);
        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();

        for (mol = info_->beginMolecule(mi); mol != NULL;
                mol = info_->nextMolecule(mi)) {
            //change the positions of atoms which belong to the rigidbodies
            for (rb = mol->beginRigidBody(rbIter); rb != NULL;
                    rb = mol->nextRigidBody(rbIter)) {
                rb->updateAtoms();
            }
        }

        Mat3x3d orderTensor(0.0);
        vecCount = 0;

        seleMan1_.setSelectionSet(evaluator1_.evaluate());

        if (doVect_) {

            for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
                    sd1 = seleMan1_.nextSelected(ii)) {
                if (sd1->isDirectional()) {
                    Vector3d vec = sd1->getA().transpose()*V3Z;

                    vec.normalize();
                    orderTensor += outProduct(vec, vec);
                    vecCount++;
                }
            }

            orderTensor /= vecCount;

        } else {

            if (doOffset_) {

                for (sd1 = seleMan1_.beginSelected(ii); sd1 != NULL;
                        sd1 = seleMan1_.nextSelected(ii)) {

                    // This will require careful rewriting if StaticProps is
                    // ever parallelized.  For an example, see
                    // Thermo::getTaggedAtomPairDistance

                    int sd2Index = sd1->getGlobalIndex() + seleOffset_;
                    sd2 = info_->getIOIndexToIntegrableObject(sd2Index);

                    Vector3d vec = sd1->getPos() - sd2->getPos();

                    if (usePeriodicBoundaryConditions_)
                        currentSnapshot_->wrapVector(vec);

                    vec.normalize();

                    orderTensor +=outProduct(vec, vec);
                    vecCount++;
                }

                orderTensor /= vecCount;
            } else {

                seleMan2_.setSelectionSet(evaluator2_.evaluate());

                if (seleMan1_.getSelectionCount() != seleMan2_.getSelectionCount() ) {
                    sprintf( painCave.errMsg,
                             "In frame %d, the number of selected StuntDoubles are\n"
                             "\tnot the same in --sele1 and sele2\n", i);
                    painCave.severity = OPENMD_INFO;
                    painCave.isFatal = 0;
                    simError();
                }

                for (sd1 = seleMan1_.beginSelected(ii),
                        sd2 = seleMan2_.beginSelected(jj);
                        sd1 != NULL && sd2 != NULL;
                        sd1 = seleMan1_.nextSelected(ii),
                        sd2 = seleMan2_.nextSelected(jj)) {

                    Vector3d vec = sd1->getPos() - sd2->getPos();

                    if (usePeriodicBoundaryConditions_)
                        currentSnapshot_->wrapVector(vec);

                    vec.normalize();

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

  void DensityPlot::process() {
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    for (int i = 0; i < nFrames; i += step_) {
      reader.readFrame(i);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();

      for (mol = info_->beginMolecule(mi); mol != NULL; 
	   mol = info_->nextMolecule(mi)) {
        //change the positions of atoms which belong to the rigidbodies
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
	     rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }
        
      }
    
      if (evaluator_.isDynamic()) {
	seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      if (cmEvaluator_.isDynamic()) {
	cmSeleMan_.setSelectionSet(cmEvaluator_.evaluate());
      }

      Vector3d origin = calcNewOrigin();

      Mat3x3d hmat = currentSnapshot_->getHmat();
      RealType slabVolume = deltaR_ * hmat(0, 0) * hmat(1, 1);
      int k; 
      for (StuntDouble* sd = seleMan_.beginSelected(k); sd != NULL; 
	   sd = seleMan_.nextSelected(k)) {


        if (!sd->isAtom()) {
          sprintf( painCave.errMsg, 
		   "Can not calculate electron density if it is not atom\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError(); 
        }
            
        Atom* atom = static_cast<Atom*>(sd);
        GenericData* data = atom->getAtomType()->getPropertyByName("nelectron");
        if (data == NULL) {
          sprintf( painCave.errMsg, "Can not find Parameters for nelectron\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError(); 
        }
            
        DoubleGenericData* doubleData = dynamic_cast<DoubleGenericData*>(data);
        if (doubleData == NULL) {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to DoubleGenericData\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();   
        }
            
        RealType nelectron = doubleData->getData();
        LennardJonesAdapter lja = LennardJonesAdapter(atom->getAtomType());
        RealType sigma = lja.getSigma() * 0.5;
        RealType sigma2 = sigma * sigma;
            
        Vector3d pos = sd->getPos() - origin;
        for (int j =0; j < nRBins_; ++j) {
          Vector3d tmp(pos);
          RealType zdist =j * deltaR_ - halfLen_;
          tmp[2] += zdist;
          if (usePeriodicBoundaryConditions_) 
            currentSnapshot_->wrapVector(tmp);
              
          RealType wrappedZdist = tmp.z() + halfLen_;
          if (wrappedZdist < 0.0 || wrappedZdist > len_) {
            continue;
          }
              
          int which = int(wrappedZdist / deltaR_);
          density_[which] += nelectron * exp(-zdist*zdist/(sigma2*2.0)) /(slabVolume* sqrt(2*NumericConstant::PI*sigma*sigma));
              
        }            
      }        
    }
  
    int nProcessed = nFrames /step_;
    std::transform(density_.begin(), density_.end(), density_.begin(), 
		   std::bind2nd(std::divides<RealType>(), nProcessed));  
    writeDensity();
        

  
  }

  void BondAngleDistribution::process() {
    Molecule* mol;
    Atom* atom;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::AtomIterator ai;
    StuntDouble* sd;
    Vector3d vec;
    std::vector<Vector3d> bondvec;
    RealType r;    
    int nBonds;    
    int i;

    bool usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();
    
    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;
    
    nTotBonds_ = 0;
    
    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      // update the positions of atoms which belong to the rigidbodies
      
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }        
      }           
            
      // outer loop is over the selected StuntDoubles:

      for (sd = seleMan_.beginSelected(i); sd != NULL; 
           sd = seleMan_.nextSelected(i)) {

        myIndex = sd->getGlobalIndex();
        nBonds = 0;
        bondvec.clear();
        
        // inner loop is over all other atoms in the system:
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {
          for (atom = mol->beginAtom(ai); atom != NULL; 
               atom = mol->nextAtom(ai)) {

            if (atom->getGlobalIndex() != myIndex) {

              vec = sd->getPos() - atom->getPos();       

              if (usePeriodicBoundaryConditions_) 
                currentSnapshot_->wrapVector(vec);
              
              // Calculate "bonds" and make a pair list 
              
              r = vec.length();
              
              // Check to see if neighbor is in bond cutoff 
              
              if (r < rCut_) { 
                // Add neighbor to bond list's
                bondvec.push_back(vec);
                nBonds++;
                nTotBonds_++;
              }  
            }
          }
          
          
          for (int i = 0; i < nBonds-1; i++ ){
            Vector3d vec1 = bondvec[i];
            vec1.normalize();
            for(int j = i+1; j < nBonds; j++){
              Vector3d vec2 = bondvec[j];
              
              vec2.normalize();
	      
              RealType theta = acos(dot(vec1,vec2))*180.0/NumericConstant::PI;
              
              
              if (theta > 180.0){
                theta = 360.0 - theta;
              }
              int whichBin = int(theta/deltaTheta_);
              
              histogram_[whichBin] += 2;
            }
          }           
//.........这里部分代码省略.........

void TetrahedralityParamXYZ::process() {
    Molecule* mol;
    StuntDouble* sd;
    StuntDouble* sd2;
    StuntDouble* sdi;
    StuntDouble* sdj;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Vector3d vec;
    Vector3d ri, rj, rk, rik, rkj;
    RealType r;
    RealType cospsi;
    RealType Qk;
    std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
    //std::vector<std::pair<Vector3d, RealType> > qvals;
    //std::vector<std::pair<Vector3d, RealType> >::iterator qiter;
    int isd1;
    int isd2;


    int kMax = int(5.0 * gaussWidth_ / voxelSize_);
    int kSqLim = kMax*kMax;
    cerr << "gw = " << gaussWidth_ << " vS = " << voxelSize_ << " kMax = "
         << kMax << " kSqLim = " << kSqLim << "\n";

    DumpReader reader(info_, dumpFilename_);
    int nFrames = reader.getNFrames();

    for (int istep = 0; istep < nFrames; istep += step_) {
        reader.readFrame(istep);

        currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
        Mat3x3d hmat = currentSnapshot_->getHmat();
        Vector3d halfBox = Vector3d(hmat(0,0), hmat(1,1), hmat(2,2)) / 2.0;

        if (evaluator1_.isDynamic()) {
            seleMan1_.setSelectionSet(evaluator1_.evaluate());
        }

        if (evaluator2_.isDynamic()) {
            seleMan2_.setSelectionSet(evaluator2_.evaluate());
        }

        // update the positions of atoms which belong to the rigidbodies
        for (mol = info_->beginMolecule(mi); mol != NULL;
                mol = info_->nextMolecule(mi)) {
            for (rb = mol->beginRigidBody(rbIter); rb != NULL;
                    rb = mol->nextRigidBody(rbIter)) {
                rb->updateAtoms();
            }
        }

        //qvals.clear();

        // outer loop is over the selected StuntDoubles:
        for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
                sd = seleMan1_.nextSelected(isd1)) {

            myIndex = sd->getGlobalIndex();

            Qk = 1.0;
            myNeighbors.clear();

            for (sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
                    sd2 = seleMan2_.nextSelected(isd2)) {

                if (sd2->getGlobalIndex() != myIndex) {

                    vec = sd->getPos() - sd2->getPos();

                    if (usePeriodicBoundaryConditions_)
                        currentSnapshot_->wrapVector(vec);

                    r = vec.length();

                    // Check to see if neighbor is in bond cutoff

                    if (r < rCut_) {
                        myNeighbors.push_back(std::make_pair(r,sd2));
                    }
                }
            }

            // Sort the vector using predicate and std::sort
            std::sort(myNeighbors.begin(), myNeighbors.end());

            // Use only the 4 closest neighbors to do the rest of the work:

            int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
            int nang = int (0.5 * (nbors * (nbors - 1)));

            rk = sd->getPos();

            for (int i = 0; i < nbors-1; i++) {

                sdi = myNeighbors[i].second;
                ri = sdi->getPos();
                rik = rk - ri;
//.........这里部分代码省略.........

void NanoVolume::process() {
#if defined(HAVE_QHULL)
  Molecule* mol;
  RigidBody* rb;
  SimInfo::MoleculeIterator mi;
  Molecule::RigidBodyIterator rbIter;
  StuntDouble* sd;
  Vector3d vec;
  int i;

  AlphaHull* thishull = new AlphaHull(2.0);
  //ConvexHull* thishull = new ConvexHull();
  
  DumpReader reader(info_, dumpFilename_);
  int nFrames = reader.getNFrames();
  frameCounter_ = 0;

  theAtoms_.reserve(info_->getNGlobalAtoms());

  for (int istep = 0; istep < nFrames; istep += step_) {
    reader.readFrame(istep);
    frameCounter_++;
    currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    RealType time = currentSnapshot_->getTime();
    
    // Clear pos vector between each frame.
    theAtoms_.clear();
    
    if (evaluator_.isDynamic()) {
      seleMan_.setSelectionSet(evaluator_.evaluate());
    }
    
    // update the positions of atoms which belong to the rigid bodies
    
    for (mol = info_->beginMolecule(mi); mol != NULL;
	 mol = info_->nextMolecule(mi)) {
      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
	   rb = mol->nextRigidBody(rbIter)) {	
	rb->updateAtoms();
      }
    }
    
    // outer loop is over the selected StuntDoubles:
    
    for (sd = seleMan_.beginSelected(i); sd != NULL;
         sd = seleMan_.nextSelected(i)) {      
      theAtoms_.push_back(sd);      
    }
    
    /* variant below for signle atoms, not StuntDoubles:
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {
      for (atom = mol->beginAtom(ai); atom != NULL; 
           atom = mol->nextAtom(ai)) {
        theAtoms_.push_back(atom);
      }
    }
    */

    // Generate convex hull for this frame.
    thishull->computeHull(theAtoms_);
    RealType volume = thishull->getVolume();

    osq.precision(7);
    if (osq.is_open()){
      osq << time << "\t" << volume << std::endl;      
    }
  }
  osq.close();

#else
  sprintf(painCave.errMsg, "NanoVolume: qhull support was not compiled in!\n");
  painCave.isFatal = 1;
  simError();  
#endif

}

  void TetrahedralityParam::process() {
    Molecule* mol;
    StuntDouble* sd;
    StuntDouble* sd2;
    StuntDouble* sdi;
    StuntDouble* sdj;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Molecule::IntegrableObjectIterator ioi;
    Vector3d vec;
    Vector3d ri, rj, rk, rik, rkj, dposition, tposition;
    RealType r;
    RealType cospsi;
    RealType Qk;
    std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
    int isd;
    bool usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();
    frameCounter_ = 0;

    Distorted_.clear();
    Tetrahedral_.clear();

    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      frameCounter_++;
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      // update the positions of atoms which belong to the rigidbodies
      
      for (mol = info_->beginMolecule(mi); mol != NULL; 
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }        
      }           
            

      // outer loop is over the selected StuntDoubles:

      for (sd = seleMan_.beginSelected(isd); sd != NULL; 
           sd = seleMan_.nextSelected(isd)) {
	
        myIndex = sd->getGlobalIndex();
	Qk = 1.0;

	myNeighbors.clear();
                
        // inner loop is over all StuntDoubles in the system:
        
        for (mol = info_->beginMolecule(mi); mol != NULL; 
             mol = info_->nextMolecule(mi)) {

          for (sd2 = mol->beginIntegrableObject(ioi); sd2 != NULL; 
               sd2 = mol->nextIntegrableObject(ioi)) {
	    
            if (sd2->getGlobalIndex() != myIndex) {
	      
              vec = sd->getPos() - sd2->getPos();       
	      
              if (usePeriodicBoundaryConditions_) 
                currentSnapshot_->wrapVector(vec);
	      
              r = vec.length();             

              // Check to see if neighbor is in bond cutoff 
              
              if (r < rCut_) { 
		
		myNeighbors.push_back(std::make_pair(r,sd2));
	      }
	    }
	  }
	}

	// Sort the vector using predicate and std::sort
	std::sort(myNeighbors.begin(), myNeighbors.end());

	//std::cerr << myNeighbors.size() <<  " neighbors within " 
	//          << rCut_  << " A" << " \n";
	
	// Use only the 4 closest neighbors to do the rest of the work:
	
	int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
	int nang = int (0.5 * (nbors * (nbors - 1)));

	rk = sd->getPos();
	//std::cerr<<nbors<<endl;
	for (int i = 0; i < nbors-1; i++) {	  

	  sdi = myNeighbors[i].second;
//.........这里部分代码省略.........

  void TetrahedralityParamZ::process() {
    Molecule* mol;
    StuntDouble* sd;
    StuntDouble* sd2;
    StuntDouble* sdi;
    StuntDouble* sdj;
    RigidBody* rb;
    int myIndex;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    Vector3d vec;
    Vector3d ri, rj, rk, rik, rkj;
    RealType r;
    RealType cospsi;
    RealType Qk;
    std::vector<std::pair<RealType,StuntDouble*> > myNeighbors;
    int isd1;
    int isd2;

    DumpReader reader(info_, dumpFilename_);    
    int nFrames = reader.getNFrames();

    for (int istep = 0; istep < nFrames; istep += step_) {
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      Mat3x3d hmat = currentSnapshot_->getHmat();
      zBox_.push_back(hmat(2,2));
      
      RealType halfBoxZ_ = hmat(2,2) / 2.0;      

      if (evaluator1_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator1_.evaluate());
      }
      
      if (evaluator2_.isDynamic()) {
        seleMan2_.setSelectionSet(evaluator2_.evaluate());
      }
      
      // update the positions of atoms which belong to the rigidbodies
      for (mol = info_->beginMolecule(mi); mol != NULL;
           mol = info_->nextMolecule(mi)) {
        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
             rb = mol->nextRigidBody(rbIter)) {
          rb->updateAtoms();
        }
      }
      
      // outer loop is over the selected StuntDoubles:
      for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
           sd = seleMan1_.nextSelected(isd1)) {
        
        myIndex = sd->getGlobalIndex();
        
        Qk = 1.0;	  
        myNeighbors.clear();       

        for (sd2 = seleMan2_.beginSelected(isd2); sd2 != NULL;
             sd2 = seleMan2_.nextSelected(isd2)) {
          
          if (sd2->getGlobalIndex() != myIndex) {
            
            vec = sd->getPos() - sd2->getPos();       
            
            if (usePeriodicBoundaryConditions_) 
              currentSnapshot_->wrapVector(vec);
            
            r = vec.length();             
            
            // Check to see if neighbor is in bond cutoff 
            
            if (r < rCut_) {                
              myNeighbors.push_back(std::make_pair(r,sd2));
            }
          }
        }
        
        // Sort the vector using predicate and std::sort
        std::sort(myNeighbors.begin(), myNeighbors.end());
        
        // Use only the 4 closest neighbors to do the rest of the work:
        
        int nbors =  myNeighbors.size()> 4 ? 4 : myNeighbors.size();
        int nang = int (0.5 * (nbors * (nbors - 1)));
        
        rk = sd->getPos();
        
        for (int i = 0; i < nbors-1; i++) {       
          
          sdi = myNeighbors[i].second;
          ri = sdi->getPos();
          rik = rk - ri;
          if (usePeriodicBoundaryConditions_) 
            currentSnapshot_->wrapVector(rik);
          
          rik.normalize();
          
          for (int j = i+1; j < nbors; j++) {       
            
            sdj = myNeighbors[j].second;
//.........这里部分代码省略.........

  void GCN::process() {
    SelectionManager common(info_);
    
    std::vector<std::vector<int> > listNN;
    std::vector<int> globalToLocal;

    Molecule* mol;
    RigidBody* rb;
    StuntDouble* sd1;
    StuntDouble* sd2;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;

    Snapshot* currentSnapshot_;
    bool usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions();

    int iterator1;
    int iterator2;
    unsigned int mapIndex1(0);
    unsigned int mapIndex2(0);
    unsigned int tempIndex(0);
    unsigned int whichBin(0);
    RealType gcn(0.0);
    Vector3d pos1;
    Vector3d pos2;
    Vector3d diff;
    RealType distance;

    histogram_.clear();
    histogram_.resize(bins_, 0.0);

    DumpReader reader(info_, dumpFilename_);
    int nFrames = reader.getNFrames();

    //First have to calculate lists of nearest neighbors (listNN_):
    
    for(int istep = 0; istep < nFrames; istep += step_){
      reader.readFrame(istep);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
      
      for(mol = info_->beginMolecule(mi); mol != NULL;
          mol = info_->nextMolecule(mi)){
	for(rb = mol->beginRigidBody(rbIter); rb != NULL;
            rb = mol->nextRigidBody(rbIter)){
	  rb->updateAtoms();
	}
      }

      if (evaluator1_.isDynamic()) {
        seleMan1_.setSelectionSet(evaluator1_.evaluate());
	selectionCount1_ = seleMan1_.getSelectionCount();
      }
      if (evaluator2_.isDynamic()) {
        seleMan2_.setSelectionSet(evaluator2_.evaluate());
	selectionCount2_ = seleMan2_.getSelectionCount();
      }

      // We need a common selection set:
      common = seleMan1_ | seleMan2_;
      int commonCount = common.getSelectionCount();

      globalToLocal.clear();
      globalToLocal.resize(info_->getNGlobalAtoms() +
                           info_->getNGlobalRigidBodies(), -1);
      for (unsigned int i = 0; i < listNN.size(); i++)         
        listNN.at(i).clear();
      listNN.clear();
      listNN.resize(commonCount);

      mapIndex1 = 0;
      for(sd1 = common.beginSelected(iterator1); sd1 != NULL;
          sd1 = common.nextSelected(iterator1)) {
        
	globalToLocal.at(sd1->getGlobalIndex()) = mapIndex1;

        pos1 = sd1->getPos();

	mapIndex2 = 0;
 	for(sd2 = common.beginSelected(iterator2); sd2 != NULL;
            sd2 = common.nextSelected(iterator2)) {
          
	  if (mapIndex1 < mapIndex2) {
            pos2 = sd2->getPos();
            diff = pos2 - pos1;
	    if (usePeriodicBoundaryConditions_)
              currentSnapshot_->wrapVector(diff);
            distance = diff.length();
	    if (distance < rCut_) {
              listNN.at(mapIndex1).push_back(mapIndex2);
	      listNN.at(mapIndex2).push_back(mapIndex1);
	    }
	  }
	  mapIndex2++;
	}
	mapIndex1++;
      }
      
      // Fill up the histogram with gcn values
      for(sd1 = seleMan1_.beginSelected(iterator1); sd1 != NULL;
          sd1 = seleMan1_.nextSelected(iterator1)){
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  std::vector<int> nMol;
  for (vector<Component*>::iterator i = components.begin(); 
       i !=components.end(); ++i) {
    int nMolOld = (*i)->getNMol();
    int nMolNew = nMolOld * repeat.x() * repeat.y() * repeat.z();    
    nMol.push_back(nMolNew);
  }
  
  createMdFile(dumpFileName, outFileName, nMol);

  SimCreator newCreator;
  SimInfo* newInfo = newCreator.createSim(outFileName, false);

  DumpReader* dumpReader = new DumpReader(oldInfo, dumpFileName);
  int nframes = dumpReader->getNFrames();
  
  DumpWriter* writer = new DumpWriter(newInfo, outFileName);
  if (writer == NULL) {
    sprintf(painCave.errMsg, "error in creating DumpWriter");
    painCave.isFatal = 1;
    simError();
  }

  SimInfo::MoleculeIterator miter;
  Molecule::IntegrableObjectIterator  iiter;
  Molecule::RigidBodyIterator rbIter;
  Molecule* mol;
  StuntDouble* sd;
  StuntDouble* sdNew;
  RigidBody* rb;
  Mat3x3d oldHmat;
  Mat3x3d newHmat;
  Snapshot* oldSnap;
  Snapshot* newSnap;
  Vector3d oldPos;
  Vector3d newPos;
  
  for (int i = 0; i < nframes; i++){
    cerr << "frame = " << i << "\n";
    dumpReader->readFrame(i);        
    oldSnap = oldInfo->getSnapshotManager()->getCurrentSnapshot();
    newSnap = newInfo->getSnapshotManager()->getCurrentSnapshot();

    newSnap->setID( oldSnap->getID() );
    newSnap->setTime( oldSnap->getTime() );
    
    oldHmat = oldSnap->getHmat();
    newHmat = repeatD*oldHmat;
    newSnap->setHmat(newHmat);

    newSnap->setThermostat( oldSnap->getThermostat() );
    newSnap->setBarostat( oldSnap->getBarostat() );

    int newIndex = 0;
    for (mol = oldInfo->beginMolecule(miter); mol != NULL; 
         mol = oldInfo->nextMolecule(miter)) {
      
      for (int ii = 0; ii < repeat.x(); ii++) {
        for (int jj = 0; jj < repeat.y(); jj++) {
          for (int kk = 0; kk < repeat.z(); kk++) {

            Vector3d trans = Vector3d(ii, jj, kk);
            for (sd = mol->beginIntegrableObject(iiter); sd != NULL;
                 sd = mol->nextIntegrableObject(iiter)) {
              oldPos = sd->getPos() + translate;
              oldSnap->wrapVector(oldPos);
              newPos = oldPos + trans * oldHmat;
              sdNew = newInfo->getIOIndexToIntegrableObject(newIndex);
              sdNew->setPos( newPos );
              sdNew->setVel( sd->getVel() );
              if (sd->isDirectional()) {
                sdNew->setA( sd->getA() );
                sdNew->setJ( sd->getJ() );
              }
              newIndex++;
            }
          }
        }
      }      
    }
  
    //update atoms of rigidbody
    for (mol = newInfo->beginMolecule(miter); mol != NULL; 
         mol = newInfo->nextMolecule(miter)) {
      
      //change the positions of atoms which belong to the rigidbodies
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
           rb = mol->nextRigidBody(rbIter)) {
        
        rb->updateAtoms();
        rb->updateAtomVel();
      }
    }

    writer->writeDump();    
  }
  // deleting the writer will put the closing at the end of the dump file.
  delete writer;
  delete oldInfo;
}

  void PotDiff::process() {
    Molecule* mol;
    RigidBody* rb;
    SimInfo::MoleculeIterator mi;
    Molecule::RigidBodyIterator rbIter;
    StuntDouble* sd;
    int j;
  
    diff_.clear();
    DumpReader reader(info_, dumpFilename_);
    int nFrames = reader.getNFrames();

    // We'll need the force manager to compute the potential
    
    ForceManager* forceMan = new ForceManager(info_);

    // We'll need thermo to report the potential
    
    Thermo* thermo =  new Thermo(info_);

    for (int i = 0; i < nFrames; i += step_) {
      reader.readFrame(i);
      currentSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
    
      for (mol = info_->beginMolecule(mi); mol != NULL; 
	   mol = info_->nextMolecule(mi)) {
	//change the positions of atoms which belong to the rigidbodies
	for (rb = mol->beginRigidBody(rbIter); rb != NULL; 
	     rb = mol->nextRigidBody(rbIter)) {
	  rb->updateAtoms();        
	}      
      }

      for (sd = seleMan_.beginSelected(j); sd != NULL;
           sd = seleMan_.nextSelected(j)) {
        if (!selectionWasFlucQ_[j])  {
          sd->setFlucQPos(0.0);
        }
      }
            
      forceMan->calcForces();
      RealType pot1 = thermo->getPotential();    

      if (evaluator_.isDynamic()) {
        seleMan_.setSelectionSet(evaluator_.evaluate());
      }

      for (sd = seleMan_.beginSelected(j); sd != NULL;
           sd = seleMan_.nextSelected(j)) {

        AtomType* at = static_cast<Atom*>(sd)->getAtomType();

        FixedChargeAdapter fca = FixedChargeAdapter(at);
        FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(at);

        RealType charge = 0.0;
        
        if (fca.isFixedCharge()) charge += fca.getCharge();
        if (fqa.isFluctuatingCharge()) charge += sd->getFlucQPos();

        sd->setFlucQPos(-charge);
      }

      currentSnapshot_->clearDerivedProperties();
      forceMan->calcForces();
      RealType pot2 = thermo->getPotential();
      RealType diff = pot2-pot1;
      
      data_.add(diff);
      diff_.push_back(diff);
      times_.push_back(currentSnapshot_->getTime());

      info_->getSnapshotManager()->advance();
    }
   
    writeDiff();   
  }