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

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

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

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

//.........这里部分代码省略.........
      surfaceMesh_ = new ConvexHull();
      break;
    case hullAlphaShape :
      surfaceMesh_ = new AlphaHull(simParams->getAlpha());
      break;
    case hullUnknown :
    default :
      sprintf(painCave.errMsg, 
              "LangevinHallForceManager: Unknown Hull_Method was requested!\n");
      painCave.isFatal = 1;
      simError();      
      break;
    }

    doThermalCoupling_ = true;
    doPressureCoupling_ = true;

    /* Check that the simulation has target pressure and target
       temperature set */    
    if (!simParams->haveTargetTemp()) {
      sprintf(painCave.errMsg, 
              "LangevinHullForceManager: no targetTemp (K) was set.\n"
              "\tOpenMD is turning off the thermal coupling to the bath.\n");
      painCave.isFatal = 0;
      painCave.severity = OPENMD_INFO;
      simError();
      doThermalCoupling_ = false;
    } else {
      targetTemp_ = simParams->getTargetTemp();

      if (!simParams->haveViscosity()) {
        sprintf(painCave.errMsg, 
                "LangevinHullForceManager: no viscosity was set.\n"
                "\tOpenMD is turning off the thermal coupling to the bath.\n");
        painCave.isFatal = 0;
        painCave.severity = OPENMD_INFO;
        simError();      
        doThermalCoupling_ = false;
      }else{
        viscosity_ = simParams->getViscosity();
        if ( fabs(viscosity_) < 1e-6 ) {
          sprintf(painCave.errMsg, 
                  "LangevinHullDynamics: The bath viscosity was set lower\n"
                  "\tthan 1e-6 poise.  OpenMD is turning off the thermal\n"
                  "\tcoupling to the bath.\n");
          painCave.isFatal = 0;
          painCave.severity = OPENMD_INFO;
          simError();
          doThermalCoupling_ = false;
        }      
      }      
    }
    if (!simParams->haveTargetPressure()) {
      sprintf(painCave.errMsg, 
              "LangevinHullForceManager: no targetPressure (atm) was set.\n"
              "\tOpenMD is turning off the pressure coupling to the bath.\n");
      painCave.isFatal = 0;
      painCave.severity = OPENMD_INFO;
      simError();
      doPressureCoupling_ = false;
    } else {
      // Convert pressure from atm -> amu/(fs^2*Ang)
      targetPressure_ = simParams->getTargetPressure() / 
        PhysicalConstants::pressureConvert;
    }
    if (simParams->getUsePeriodicBoundaryConditions()) {
      sprintf(painCave.errMsg, 
              "LangevinHallForceManager: You can't use the Langevin Hull\n"
              "\tintegrator with periodic boundary conditions turned on!\n");
      painCave.isFatal = 1;
      simError();
    } 

    dt_ = simParams->getDt();

    if (doThermalCoupling_)
      variance_ = 2.0 * PhysicalConstants::kb * targetTemp_ / dt_;

    // Build a vector of integrable objects to determine if the are
    // surface atoms
    Molecule* mol;
    StuntDouble* sd;
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {          
      for (sd = mol->beginIntegrableObject(j); 
           sd != NULL;
           sd = mol->nextIntegrableObject(j)) {	
	localSites_.push_back(sd);
      }
    }
    
    // We need to make an initial guess at the bounding box in order
    // to compute long range forces in ForceMatrixDecomposition:

    // Compute surface Mesh
    surfaceMesh_->computeHull(localSites_);
  }  

  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 ThermoIntegrationForceManager::calcForces(){
    Snapshot* curSnapshot;
    SimInfo::MoleculeIterator mi;
    Molecule* mol;
    Molecule::IntegrableObjectIterator ii;
    StuntDouble* sd;
    Vector3d frc;
    Vector3d trq;
    Mat3x3d tempTau;
    
    // perform the standard calcForces first
    ForceManager::calcForces();
    
    curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();

    // now scale forces and torques of all the sds
      
    for (mol = info_->beginMolecule(mi); mol != NULL; 
         mol = info_->nextMolecule(mi)) {

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

        frc = sd->getFrc();
        frc *= factor_;
        sd->setFrc(frc);
        
        if (sd->isDirectional()){
          trq = sd->getTrq();
          trq *= factor_;
          sd->setTrq(trq);
        }
      }
    }
    
    // set rawPotential to be the unmodulated potential
    lrPot_ = curSnapshot->getLongRangePotential();
    curSnapshot->setRawPotential(lrPot_);
    
    // modulate the potential and update the snapshot
    lrPot_ *= factor_;
    curSnapshot->setLongRangePotential(lrPot_);
    
    // scale the pressure tensor
    tempTau = curSnapshot->getStressTensor();
    tempTau *= factor_;
    curSnapshot->setStressTensor(tempTau);

    // now, on to the applied restraining potentials (if needed):
    RealType restPot_local = 0.0;
    RealType vHarm_local = 0.0;
    
    if (simParam->getUseRestraints()) {
      // do restraints from RestraintForceManager:
      restPot_local = doRestraints(1.0 - factor_);      
      vHarm_local = getUnscaledPotential();
    }
      
#ifdef IS_MPI
    RealType restPot;
    MPI::COMM_WORLD.Allreduce(&restPot_local, &restPot, 1, 
                              MPI::REALTYPE, MPI::SUM);
    MPI::COMM_WORLD.Allreduce(&vHarm_local, &vHarm_, 1, 
                              MPI::REALTYPE, MPI::SUM);         
    lrPot_ += restPot;
#else
    lrPot_ += restPot_local;
    vHarm_ = vHarm_local;
#endif

    // give the final values to stats
    curSnapshot->setLongRangePotential(lrPot_);
    curSnapshot->setRestraintPotential(vHarm_);
  }  

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