C++ StuntDouble类代码示例

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

}

  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;
  }

  void NVE::moveB(){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* sd;
    Vector3d vel;
    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();
	frc = sd->getFrc();
	mass = sd->getMass();
                
	// velocity half step
	vel += (dt2 /mass * PhysicalConstants::energyConvert) * frc;
                
	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;

	  sd->setJ(ji);
	}

            
      }
    }
  
    flucQ_->moveB();
    rattle_->constraintB();
  }

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

    velCom /= totalMass;

    return velCom;
}

  void RCorrFuncZ::computeFrame(int istep) {
    hmat_ = currentSnapshot_->getHmat();
    halfBoxZ_ = hmat_(2,2) / 2.0;      

    StuntDouble* sd;

    int isd1, isd2;
    unsigned int index;

    if (evaluator1_.isDynamic()) {
      seleMan1_.setSelectionSet(evaluator1_.evaluate());
    }
    
    if (uniqueSelections_ && evaluator2_.isDynamic()) {
      seleMan2_.setSelectionSet(evaluator2_.evaluate());
    }      
    
    for (sd = seleMan1_.beginSelected(isd1); sd != NULL;
         sd = seleMan1_.nextSelected(isd1)) {

      index = computeProperty1(istep, sd);        
      if (index == sele1ToIndex_[istep].size()) {
        sele1ToIndex_[istep].push_back(sd->getGlobalIndex());
      } else {
        sele1ToIndex_[istep].resize(index+1);
        sele1ToIndex_[istep][index] = sd->getGlobalIndex();
      }
    }
   
    if (uniqueSelections_) { 
      for (sd = seleMan2_.beginSelected(isd2); sd != NULL;
           sd = seleMan2_.nextSelected(isd2)) {        

        index = computeProperty1(istep, sd);

        if (index == sele2ToIndex_[istep].size()) {
          sele2ToIndex_[istep].push_back(sd->getGlobalIndex());
        } else {
          sele2ToIndex_[istep].resize(index+1);
          sele2ToIndex_[istep][index] = sd->getGlobalIndex();
        }
      }
    }
  }

  PotDiff::PotDiff(SimInfo* info, const std::string& filename, 
                   const std::string& sele)
    : StaticAnalyser(info, filename), selectionScript_(sele), 
      seleMan_(info), evaluator_(info) {
    
    StuntDouble* sd;
    int i;
    
    setOutputName(getPrefix(filename) + ".potDiff");

    // The PotDiff is computed by negating the charge on the atom type
    // using fluctuating charge values.  If we don't have any
    // fluctuating charges in the simulation, we need to expand
    // storage to hold them.
    int storageLayout = info_->getStorageLayout();
    storageLayout |= DataStorage::dslFlucQPosition;
    storageLayout |= DataStorage::dslFlucQVelocity;
    storageLayout |= DataStorage::dslFlucQForce;
    info_->setStorageLayout(storageLayout);
    info_->setSnapshotManager(new SimSnapshotManager(info_, storageLayout));

    // now we have to figure out which AtomTypes to convert to fluctuating
    // charges
    evaluator_.loadScriptString(sele);    
    seleMan_.setSelectionSet(evaluator_.evaluate());
    for (sd = seleMan_.beginSelected(i); sd != NULL;
         sd = seleMan_.nextSelected(i)) {      
      AtomType* at = static_cast<Atom*>(sd)->getAtomType();
      FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(at);
      if (fqa.isFluctuatingCharge()) {
        selectionWasFlucQ_.push_back(true);
      } else {
        selectionWasFlucQ_.push_back(false);                
        // make a fictitious fluctuating charge with an unphysical
        // charge mass and slaterN, but we need to zero out the
        // electronegativity and hardness to remove the self
        // contribution:
        fqa.makeFluctuatingCharge(1.0e9, 0.0, 0.0, 1);
        sd->setFlucQPos(0.0);
      }
    }
    info_->getSnapshotManager()->advance();
  }

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

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

    thermostat = snap->getThermostat();
    RealType oldChi  = thermostat.first;
    RealType prevChi;

    loadEta();
    
    //save velocity and angular momentum
    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
	   sd = mol->nextIntegrableObject(j)) {
                
	oldVel[index] = sd->getVel();

        if (sd->isDirectional())
	   oldJi[index] = sd->getJ();

	++index;
      }
    }

    // do the iteration:
    instaVol =thermo.getVolume();

    for(int k = 0; k < maxIterNum_; k++) {
      instaTemp =thermo.getTemperature();
      instaPress =thermo.getPressure();

      // evolve chi another half step using the temperature at t + dt/2
      prevChi = thermostat.first;
      thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;

      //evolve eta
      this->evolveEtaB();
      this->calcVelScale();

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

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

	  frc = sd->getFrc();
	  mass = sd->getMass();

	  getVelScaleB(sc, index);

	  // velocity half step
	  vel = oldVel[index] 
            + 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());

	    ji = oldJi[index] 
              + dt2*PhysicalConstants::energyConvert*Tb 
              - dt2*thermostat.first*oldJi[index];

	    sd->setJ(ji);
	  }

	  ++index;
	}
      }
        
      rattle_->constraintB();

      if ((fabs(prevChi - thermostat.first) <= chiTolerance) && 
          this->etaConverged())
	break;
    }

    //calculate integral of chidt
    thermostat.second += dt2 * thermostat.first;

    snap->setThermostat(thermostat);

    flucQ_->moveB();
    saveEta();
  }

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

  
  }

int main(int argc, char* argv[]){
  
  gengetopt_args_info args_info;
  string dumpFileName;
  string outFileName;
  
  //parse the command line option
  if (cmdline_parser (argc, argv, &args_info) != 0) {
    exit(1) ;
  }
  
  //get the dumpfile name and meta-data file name
  if (args_info.input_given){
    dumpFileName = args_info.input_arg;
  } else {
    strcpy( painCave.errMsg,
            "No input file name was specified.\n" );
    painCave.isFatal = 1;
    simError();
  }
  
  if (args_info.output_given){
    outFileName = args_info.output_arg;
  } else {
    strcpy( painCave.errMsg,
            "No output file name was specified.\n" );
    painCave.isFatal = 1;
    simError();
  }

  Vector3i repeat = Vector3i(args_info.repeatX_arg,
                             args_info.repeatY_arg,
                             args_info.repeatZ_arg);
  Mat3x3d repeatD = Mat3x3d(0.0);
  repeatD(0,0) = repeat.x();
  repeatD(1,1) = repeat.y();
  repeatD(2,2) = repeat.z();

  Vector3d translate = Vector3d(args_info.translateX_arg,
                                args_info.translateY_arg,
                                args_info.translateZ_arg);

  //parse md file and set up the system

  SimCreator oldCreator;
  SimInfo* oldInfo = oldCreator.createSim(dumpFileName, false);
  Globals* simParams = oldInfo->getSimParams();
  std::vector<Component*> components = simParams->getComponents();
  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() );
//.........这里部分代码省略.........

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

int main(int argc, char* argv[]){
  registerHydrodynamicsModels();
  
  gengetopt_args_info args_info;
  std::string dumpFileName;
  std::string mdFileName;
  std::string prefix;
  
  //parse the command line option
  if (cmdline_parser (argc, argv, &args_info) != 0) {
    exit(1) ;
  }
  
  //get the dumpfile name and meta-data file name
  if (args_info.input_given){
    dumpFileName = args_info.input_arg;
  } else {
    strcpy( painCave.errMsg,
            "No input file name was specified.\n" );
    painCave.isFatal = 1;
    simError();
  }
  
  if (args_info.output_given){
    prefix = args_info.output_arg;
  } else {
    prefix = getPrefix(dumpFileName);    
  }
  std::string outputFilename = prefix + ".diff";
    
  //parse md file and set up the system
  SimCreator creator;
  SimInfo* info = creator.createSim(dumpFileName, true);
    
  SimInfo::MoleculeIterator mi;
  Molecule* mol;
  Molecule::IntegrableObjectIterator  ii;
  StuntDouble* sd;
  Mat3x3d identMat;
  identMat(0,0) = 1.0;
  identMat(1,1) = 1.0;
  identMat(2,2) = 1.0;

  Globals* simParams = info->getSimParams();
  RealType temperature(0.0);
  RealType viscosity(0.0);

  if (simParams->haveViscosity()) {
    viscosity = simParams->getViscosity();
  } else {
    sprintf(painCave.errMsg, "viscosity must be set\n");
    painCave.isFatal = 1;
    simError();  
  }

  if (simParams->haveTargetTemp()) {
    temperature = simParams->getTargetTemp();
  } else {
    sprintf(painCave.errMsg, "target temperature must be set\n");
    painCave.isFatal = 1;
    simError();  
  }
 
  std::map<std::string, SDShape> uniqueStuntDoubles;
  
  for (mol = info->beginMolecule(mi); mol != NULL; 
       mol = info->nextMolecule(mi)) {
    
    for (sd = mol->beginIntegrableObject(ii); sd != NULL;
         sd = mol->nextIntegrableObject(ii)) {
      
      if (uniqueStuntDoubles.find(sd->getType()) ==  uniqueStuntDoubles.end()) {
        
        sd->setPos(V3Zero);
        sd->setA(identMat);
        if (sd->isRigidBody()) {
          RigidBody* rb = static_cast<RigidBody*>(sd);
          rb->updateAtoms();
        }
        
        SDShape tmp;
        tmp.shape = ShapeBuilder::createShape(sd);
        tmp.sd = sd;    
        uniqueStuntDoubles.insert(std::map<std::string, SDShape>::value_type(sd->getType(), tmp));
        
      }
    }
  }
  
  
  
  std::map<std::string, SDShape>::iterator si;
  for (si = uniqueStuntDoubles.begin(); si != uniqueStuntDoubles.end(); ++si) {
    HydrodynamicsModel* model;
    Shape* shape = si->second.shape;
    StuntDouble* sd = si->second.sd;;
    if (args_info.model_given) {  
      model = HydrodynamicsModelFactory::getInstance()->createHydrodynamicsModel(args_info.model_arg, sd, info);
    } else if (shape->hasAnalyticalSolution()) {
      model = new AnalyticalModel(sd, info);
//.........这里部分代码省略.........

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