C++ StuntDouble::getPos方法代码示例

void Molecule::moveCom(const Vetor3d& delta) {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
        s->setPos(sd->getPos() + delta);
    }

}

  void NVE::moveA(){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* sd;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d Tb;
    Vector3d ji;
    RealType mass;
    
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
	   sd = mol->nextIntegrableObject(j)) {

	vel = sd->getVel();
	pos = sd->getPos();
	frc = sd->getFrc();
	mass = sd->getMass();
                
	// velocity half step
	vel += (dt2 /mass * PhysicalConstants::energyConvert) * frc;

	// position whole step
	pos += dt * vel;

	sd->setVel(vel);
	sd->setPos(pos);

	if (sd->isDirectional()){

	  // get and convert the torque to body frame

	  Tb = sd->lab2Body(sd->getTrq());

	  // get the angular momentum, and propagate a half step

	  ji = sd->getJ();

	  ji += (dt2  * PhysicalConstants::energyConvert) * Tb;

	  rotAlgo_->rotate(sd, ji, dt);

	  sd->setJ(ji);
	}

            
      }
    }
    flucQ_->moveA();
    rattle_->constraintA();    
  }    

  Vector3d DensityPlot::calcNewOrigin() {

    int i;
    Vector3d newOrigin(0.0);
    RealType totalMass = 0.0;
    for (StuntDouble* sd = seleMan_.beginSelected(i); sd != NULL; sd = seleMan_.nextSelected(i)) {
      RealType mass = sd->getMass();
      totalMass += mass;
      newOrigin += sd->getPos() * mass;        
    }
    newOrigin /= totalMass;
    return newOrigin;
  }

Vector3d Molecule::getCom() {
    StuntDouble* sd;
    std::vector<StuntDouble*>::iterator i;
    Vector3d com;
    double totalMass = 0;
    double mass;
    
    for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){
        mass = sd->getMass();
        totalMass += mass;
        com += sd->getPos() * mass;    
    }

    com /= totalMass;

    return com;
}

  void NPT::moveA() {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* sd;
    Vector3d Tb, ji;
    RealType mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d sc;
    int index;

    thermostat = snap->getThermostat();
    loadEta();
    
    instaTemp =thermo.getTemperature();
    press = thermo.getPressureTensor();
    instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
    instaVol =thermo.getVolume();

    Vector3d  COM = thermo.getCom();

    //evolve velocity half step

    calcVelScale();

    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
	   sd = mol->nextIntegrableObject(j)) {
                
	vel = sd->getVel();
	frc = sd->getFrc();

	mass = sd->getMass();

	getVelScaleA(sc, vel);

	// velocity half step  (use chi from previous step here):

	vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
	sd->setVel(vel);

	if (sd->isDirectional()) {

	  // get and convert the torque to body frame

	  Tb = sd->lab2Body(sd->getTrq());

	  // get the angular momentum, and propagate a half step

	  ji = sd->getJ();

	  ji += dt2*PhysicalConstants::energyConvert * Tb 
            - dt2*thermostat.first* ji;
                
	  rotAlgo_->rotate(sd, ji, dt);

	  sd->setJ(ji);
	}
            
      }
    }
    // evolve chi and eta  half step

    thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
    
    evolveEtaA();

    //calculate the integral of chidt
    thermostat.second += dt2 * thermostat.first;
    
    flucQ_->moveA();


    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
	   sd = mol->nextIntegrableObject(j)) {

	oldPos[index++] = sd->getPos();            

      }
    }
    
    //the first estimation of r(t+dt) is equal to  r(t)

    for(int k = 0; k < maxIterNum_; k++) {
      index = 0;
      for (mol = info_->beginMolecule(i); mol != NULL; 
           mol = info_->nextMolecule(i)) {

	for (sd = mol->beginIntegrableObject(j); sd != NULL;
	     sd = mol->nextIntegrableObject(j)) {

	  vel = sd->getVel();
//.........这里部分代码省略.........

  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();
        

  
  }

  SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance, int frame ) {
    StuntDouble * center;
    Vector3d centerPos;
    Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame);
    SelectionSet bsResult(nObjects_);   
    assert(bsResult.size() == bs.size());

#ifdef IS_MPI
    int mol;
    int proc;
    RealType data[3];
    int worldRank = MPI::COMM_WORLD.Get_rank();
#endif
 
    for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
      if (stuntdoubles_[j]->isRigidBody()) {
        RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
        rb->updateAtoms(frame);
      }
    }
       
    SelectionSet bsTemp(nObjects_);
    bsTemp = bs;
    bsTemp.parallelReduce();

    for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1; 
         i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) {

      // Now, if we own stuntdouble i, we can use the position, but in
      // parallel, we'll need to let everyone else know what that
      // position is!

#ifdef IS_MPI
      mol = info_->getGlobalMolMembership(i);
      proc = info_->getMolToProc(mol);
     
      if (proc == worldRank) {
        center = stuntdoubles_[i];
        centerPos = center->getPos(frame);
        data[0] = centerPos.x();
        data[1] = centerPos.y();
        data[2] = centerPos.z();          
        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
      } else {
        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
        centerPos = Vector3d(data);
      }
#else
      center = stuntdoubles_[i];
      centerPos = center->getPos(frame);
#endif
      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
	Vector3d r =centerPos - stuntdoubles_[j]->getPos(frame);
	currSnapshot->wrapVector(r);
	if (r.length() <= distance) {
	  bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
	}
      }
      for (unsigned int j = 0; j < bonds_.size(); ++j) {       
        Vector3d loc = bonds_[j]->getAtomA()->getPos(frame);
        loc += bonds_[j]->getAtomB()->getPos(frame);
        loc = loc / 2.0;
	Vector3d r = centerPos - loc;
	currSnapshot->wrapVector(r);
	if (r.length() <= distance) {
	  bsResult.bitsets_[BOND].setBitOn(j);
	}
      }
      for (unsigned int j = 0; j < bends_.size(); ++j) {
        Vector3d loc = bends_[j]->getAtomA()->getPos(frame);
        loc += bends_[j]->getAtomB()->getPos(frame);
        loc += bends_[j]->getAtomC()->getPos(frame);
        loc = loc / 3.0;
	Vector3d r = centerPos - loc;
	currSnapshot->wrapVector(r);
	if (r.length() <= distance) {
	  bsResult.bitsets_[BEND].setBitOn(j);
	}
      }
      for (unsigned int j = 0; j < torsions_.size(); ++j) {
        Vector3d loc = torsions_[j]->getAtomA()->getPos(frame);
        loc += torsions_[j]->getAtomB()->getPos(frame);
        loc += torsions_[j]->getAtomC()->getPos(frame);
        loc += torsions_[j]->getAtomD()->getPos(frame);
        loc = loc / 4.0;
	Vector3d r = centerPos - loc;
	currSnapshot->wrapVector(r);
	if (r.length() <= distance) {
	  bsResult.bitsets_[TORSION].setBitOn(j);
	}
      }
      for (unsigned int j = 0; j < inversions_.size(); ++j) {
        Vector3d loc = inversions_[j]->getAtomA()->getPos(frame);
        loc += inversions_[j]->getAtomB()->getPos(frame);
        loc += inversions_[j]->getAtomC()->getPos(frame);
        loc += inversions_[j]->getAtomD()->getPos(frame);
        loc = loc / 4.0;
	Vector3d r = centerPos - loc;
	currSnapshot->wrapVector(r);
	if (r.length() <= distance) {
//.........这里部分代码省略.........

  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 ContactAngle1::doFrame(int frame) {
    StuntDouble* sd;
    int i;
    
    if (evaluator1_.isDynamic()) {
      seleMan1_.setSelectionSet(evaluator1_.evaluate());
    }
    
    
    RealType mtot = 0.0;
    Vector3d com(V3Zero);
    RealType mass;
    
    for (sd = seleMan1_.beginSelected(i); sd != NULL;
         sd = seleMan1_.nextSelected(i)) {      
      mass = sd->getMass();
      mtot += mass;
      com += sd->getPos() * mass;
    }
    
    com /= mtot;

    RealType dz = com.z() - solidZ_;

    if (dz < 0.0) {
      sprintf(painCave.errMsg, 
              "ContactAngle1: Z-center of mass of selection, %lf, was\n"
              "\tlocated below the solid reference plane, %lf\n",
              com.z(), solidZ_);
      painCave.isFatal = 1;
      painCave.severity = OPENMD_ERROR;
      simError();
    }

    if (dz > dropletRadius_) {
      values_.push_back(180.0);
    } else {
    
      RealType k = pow(2.0, -4.0/3.0) * dropletRadius_;
      
      RealType z2 = dz*dz;
      RealType z3 = z2 * dz;
      RealType k2 = k*k;
      RealType k3 = k2*k;
      
      Polynomial<RealType> poly;
      poly.setCoefficient(4,      z3 +      k3);
      poly.setCoefficient(3,  8.0*z3 +  8.0*k3);
      poly.setCoefficient(2, 24.0*z3 + 18.0*k3);
      poly.setCoefficient(1, 32.0*z3          );
      poly.setCoefficient(0, 16.0*z3 - 27.0*k3);
      vector<RealType> realRoots = poly.FindRealRoots();

      RealType ct;
      
      vector<RealType>::iterator ri;


      RealType maxct = -1.0;
      for (ri = realRoots.begin(); ri !=realRoots.end(); ++ri) {
        ct = *ri;
        if (ct > 1.0)  ct = 1.0;
        if (ct < -1.0) ct = -1.0;

        // use the largest magnitude of ct that it finds:
        if (ct > maxct) {
          maxct = ct;
        }                  
      }
      
      values_.push_back( acos(maxct)*(180.0/M_PI) );
    }
  }    

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 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 ContactAngle2::doFrame(int frame) {
    StuntDouble* sd;
    int i;

    // set up the bins for density analysis

    Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
    RealType len = std::min(hmat(0, 0), hmat(1, 1));
    RealType zLen = hmat(2,2);

    RealType dr = len / (RealType) nRBins_;
    RealType dz = zLen / (RealType) nZBins_;

    std::vector<std::vector<RealType> > histo;
    histo.resize(nRBins_);
    for (unsigned int i = 0; i < histo.size(); ++i){
      histo[i].resize(nZBins_);
      std::fill(histo[i].begin(), histo[i].end(), 0.0);
    }      
        
    if (evaluator1_.isDynamic()) {
      seleMan1_.setSelectionSet(evaluator1_.evaluate());
    }
    
    Vector3d com(centroidX_, centroidY_, solidZ_);

    // now that we have the centroid, we can make cylindrical density maps
    Vector3d pos;
    RealType r;
    RealType z;
    
    for (sd = seleMan1_.beginSelected(i); sd != NULL;
         sd = seleMan1_.nextSelected(i)) {      
      pos = sd->getPos() - com;

      // r goes from zero upwards
      r = sqrt(pow(pos.x(), 2) + pow(pos.y(), 2));
      // z is possibly symmetric around 0
      z = pos.z();
          
      int whichRBin = int(r / dr);
      int whichZBin = int( (zLen/2.0 + z) / dz);
      
      if ((whichRBin < int(nRBins_)) && (whichZBin >= 0)
          && (whichZBin < int(nZBins_))) {
        histo[whichRBin][whichZBin] += sd->getMass();
      }
      
    }
    
    for(unsigned int i = 0 ; i < histo.size(); ++i){

      RealType rL = i * dr;
      RealType rU = rL + dr;
      RealType volSlice = NumericConstant::PI * dz * (( rU*rU ) - ( rL*rL ));

      for (unsigned int j = 0; j < histo[i].size(); ++j) {
        histo[i][j] *= PhysicalConstants::densityConvert / volSlice;
      }
    }

    std::vector<Vector<RealType, 2> > points;
    points.clear();
    
    for (unsigned int j = 0; j < nZBins_;  ++j) {

      // The z coordinates were measured relative to the selection
      // center of mass.  However, we're interested in the elevation
      // above the solid surface.  Also, the binning was done around
      // zero with enough bins to cover the zLength of the box:
      
      RealType thez =  com.z() - solidZ_  - zLen/2.0 + dz * (j + 0.5);
      bool aboveThresh = false;
      bool foundThresh = false;
      int rloc = 0;
      
      for (std::size_t i = 0; i < nRBins_;  ++i) {

        if (histo[i][j] >= threshDens_) aboveThresh = true;

        if (aboveThresh && (histo[i][j] <= threshDens_)) {
          rloc = i;
          foundThresh = true;
          aboveThresh = false;
        }

      }
      if (foundThresh) {
        Vector<RealType,2> point;
        point[0] = dr*(rloc+0.5);
        point[1] = thez;

        if (thez > bufferLength_) {
          points.push_back( point );
        }
      }      
    }

    int numPoints = points.size();

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

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

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