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

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

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

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

    //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();
	  pos = sd->getPos();

	  this->getPosScale(pos, COM, index, sc);

	  pos = oldPos[index] + dt * (vel + sc);
	  sd->setPos(pos);     

	  ++index;
	}
      }

      rattle_->constraintA();
    }

    // Scale the box after all the positions have been moved:

    this->scaleSimBox();

    snap->setThermostat(thermostat);

    saveEta();
  }

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