ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-4/src/constraints/ZconstraintForceManager.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-4/src/constraints/ZconstraintForceManager.cpp (file contents):
Revision 1911 by tim, Mon Jan 10 18:05:45 2005 UTC vs.
Revision 1912 by tim, Mon Jan 10 20:52:07 2005 UTC

# Line 3 | Line 3
3   #include "integrators/Integrator.hpp"
4   #include "utils/simError.h"
5   #include "utils/OOPSEConstant.hpp"
6 <
6 > #include "utils/StringUtils.hpp"
7   namespace oopse {
8   ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info) {
9      currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
# Line 45 | Line 45 | ZconstraintForceManager::ZconstraintForceManager(SimIn
45  
46      //set zcons gap
47      if (simParam->haveZConsGap()){
48 +        usingZconsGap_ = true;
49          zconsGap_ = simParam->getZconsGap();
50      }else {
51 <        zconsGap_ = false;
51 >        usingZconsGap_ = false;
52 >        zconsGap_ = 0.0;
53      }
54  
55      //set zcons fixtime
# Line 64 | Line 66 | ZconstraintForceManager::ZconstraintForceManager(SimIn
66          usingSMD_ =false;
67      }
68      
69 <    std::string zconsOutput = getPrefix(info_->getFinalConfigFileName()) + ".fz";
69 >    zconsOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
70  
71      //estimate the force constant of harmonical potential
72      Mat3x3d hmat = currSnapshot_->getHmat();
# Line 85 | Line 87 | ZconstraintForceManager::ZconstraintForceManager(SimIn
87          ZconstraintParam param;
88          int zmolIndex = stamp[i]->getMolIndex();
89          if (stamp[i]->haveZpos()) {
88 –            param.zTargetPos = getZTargetPos(zmolIndex);
89 –        } else {
90              param.zTargetPos = stamp[i]->getZpos();
91 +        } else {
92 +            param.zTargetPos = getZTargetPos(zmolIndex);
93          }
94  
95          param.kz = zforceConstant * stamp[i]->getKratio();
# Line 118 | Line 120 | ZconstraintForceManager::ZconstraintForceManager(SimIn
120   #endif
121  
122      // creat zconsWriter  
123 <    fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);  
123 >    fzOut = new ZConsWriter(info_, zconsOutput_.c_str());  
124  
125      if (!fzOut){
126 <      sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
126 >      sprintf(painCave.errMsg, "Fail to create ZConsWriter\n");
127        painCave.isFatal = 1;
128        simError();
129      }
130  
131   }
132  
133 < ZconstraintForceManager::~ZConstraint(){
133 > ZconstraintForceManager::~ZconstraintForceManager(){
134  
135    if (fzOut){
136      delete fzOut;
# Line 140 | Line 142 | void ZconstraintForceManager::update(){
142      fixedZMols_.clear();
143      movingZMols_.clear();
144      unzconsMols_.clear();
145 <    
146 <    std::map<int, ZconstraintParam>::iterator i;
145 <    for (i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
145 >
146 >    for (std::map<int, ZconstraintParam>::iterator i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
147   #ifdef IS_MPI
148          if (info_->getMolToProc(i->first) == worldRank) {
149   #endif
150              ZconstraintMol zmol;
151 <            zmol->mol = info_->getMoleculeByGlobalIndex(i->first);
152 <            assert(zmol->mol);
153 <            zmol->param = info_->getMoleculeByGlobalIndex(i->second);
154 <            zmol->cantPos = zmol->param->zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
155 <            Vector3d com = zmol->mol->getCom();
156 <            double diff = fabs(zmol->param->zTargetPos - com[whichDirection]);
151 >            zmol.mol = info_->getMoleculeByGlobalIndex(i->first);
152 >            assert(zmol.mol);
153 >            zmol.param = i->second;
154 >            zmol.cantPos = zmol.param.zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
155 >            Vector3d com = zmol.mol->getCom();
156 >            double diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
157              if (diff < zconsTol_) {
158                  fixedZMols_.push_back(zmol);
159              } else {
# Line 167 | Line 168 | void ZconstraintForceManager::update(){
168      calcTotalMassMovingZMols();
169  
170      std::set<int> zmolSet;
171 <    std::list<ZconstraintMol>::iterator i;
171 <    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
171 >    for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
172          zmolSet.insert(i->mol->getGlobalIndex());
173      }    
174  
175 <    for ( i = fixedZMols_.begin(); i !=  movingZMols_.end(); ++i) {
175 >    for (std::list<ZconstraintMol>::iterator i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
176          zmolSet.insert(i->mol->getGlobalIndex());
177      }
178  
179      SimInfo::MoleculeIterator mi;
180      Molecule* mol;
181      for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
182 <        if (zmolSet->find(mol->getGlobalIndex()) == zmolSet.end()) {
183 <            unzconsMols_.push_back(zmolSet);
182 >        if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
183 >            unzconsMols_.push_back(mol);
184          }
185      }
186  
187   }
188  
189   bool ZconstraintForceManager::isZMol(Molecule* mol){
190 <    return allZMolIndices.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
190 >    return allZMolIndices_.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
191   }
192  
193 < void ZconstraintForceManager::integrate(){
193 > void ZconstraintForceManager::init(){
194  
195 –
195    //zero out the velocities of center of mass of unconstrained molecules
196    //and the velocities of center of mass of every single z-constrained molecueles
197    zeroVelocity();
198  
199 <  curZconsTime = zconsTime + info->getTime();
201 <
202 <  T::integrate();
199 >  currZconsTime_ = currSnapshot_->getTime();
200   }
201  
202 < void ZconstraintForceManager::calcForce(bool needPotential, bool needStress){
203 <    ForceManager::calcForce(needPotential, needStress);
202 > void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
203 >    ForceManager::calcForces(needPotential, needStress);
204      
205      if (usingZconsGap_){
206          updateZPos();
# Line 211 | Line 208 | void ZconstraintForceManager::calcForce(bool needPoten
208  
209      if (checkZConsState()){
210          zeroVelocity();    
211 <        forcePolicy->update();
211 >        calcTotalMassMovingZMols();
212      }  
213  
214      //do zconstraint force;
# Line 225 | Line 222 | void ZconstraintForceManager::calcForce(bool needPoten
222      }
223  
224      //write out forces and current positions of z-constraint molecules    
225 <    if (currSnapshot_->getTime() >= curZconsTime_){
225 >    if (currSnapshot_->getTime() >= currZconsTime_){
226          std::list<ZconstraintMol>::iterator i;
227          Vector3d com;
228 <        for(i = fixedZMols_.begin(); i != fixedZMols_.end() ++i) {
228 >        for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
229              com = i->mol->getCom();
230              i->zpos = com[whichDirection];
231          }
232          
233          fzOut->writeFZ(fixedZMols_);
234 <        curZconsTime_ += zconsTime_;
234 >        currZconsTime_ += zconsTime_;
235      }
236   }
237  
# Line 245 | Line 242 | void ZconstraintForceManager::zeroVelocity(){
242      std::list<ZconstraintMol>::iterator i;
243      Molecule* mol;
244      StuntDouble* integrableObject;
245 <    Molecule::IntegrableObjectIteraot ii;
245 >    Molecule::IntegrableObjectIterator ii;
246  
247      //zero out the velocities of center of mass of fixed z-constrained molecules
248 <    for(i = fixedZMols_.begin(); i != fixedZMols_.end() ++i) {
248 >    for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
249          mol = i->mol;
250          comVel = mol->getComVel();
251          for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
252              integrableObject = mol->nextIntegrableObject(ii)) {
253 +            vel = integrableObject->getVel();  
254              vel[whichDirection] -= comVel[whichDirection];
255              integrableObject->setVel(vel);
256          }
# Line 260 | Line 258 | void ZconstraintForceManager::zeroVelocity(){
258  
259      // calculate the vz of center of mass of moving molecules(include unconstrained molecules
260      // and moving z-constrained molecules)  
261 <    double pzMovingMols_local;
261 >    double pzMovingMols_local = 0.0;
262      double pzMovingMols;
263      
264      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
# Line 272 | Line 270 | void ZconstraintForceManager::zeroVelocity(){
270      std::vector<Molecule*>::iterator j;
271      for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
272          mol =*j;
273 <        comVel = mol->getCoMVel();
273 >        comVel = mol->getComVel();
274          pzMovingMols_local += mol->getMass() * comVel[whichDirection];
275      }
276      
# Line 283 | Line 281 | void ZconstraintForceManager::zeroVelocity(){
281          MPI_SUM, MPI_COMM_WORLD);
282   #endif
283  
284 <    double vzMovingMols = pzMovingMols / totMassMovingMols_;
284 >    double vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_);
285  
286      //modify the velocities of moving z-constrained molecuels
287      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
# Line 292 | Line 290 | void ZconstraintForceManager::zeroVelocity(){
290              integrableObject = mol->nextIntegrableObject(ii)) {
291  
292              vel = integrableObject->getVel();
293 <            vel[whichDirection] -= vzOfMovingMols;
293 >            vel[whichDirection] -= vzMovingMols;
294              integrableObject->setVel(vel);
295          }
296      }
# Line 304 | Line 302 | void ZconstraintForceManager::zeroVelocity(){
302              integrableObject = mol->nextIntegrableObject(ii)) {
303  
304              vel = integrableObject->getVel();
305 <            vel[whichDirection] -= vzOfMovingMols;
305 >            vel[whichDirection] -= vzMovingMols;
306              integrableObject->setVel(vel);
307          }
308      }
# Line 329 | Line 327 | void ZconstraintForceManager::doZconstraintForce(){
327      std::list<ZconstraintMol>::iterator i;
328      Molecule* mol;
329      StuntDouble* integrableObject;
330 <    Molecule::IntegrableObjectIteraot ii;
330 >    Molecule::IntegrableObjectIterator ii;
331  
332      for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
333          mol = i->mol;
336 –        com = mol->getCom();        
337 –
334          i->fz = 0.0;
339 –
340 –        totalFZ_local += harmonicF;
341 –
342 –        //adjust force
335          for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
336              integrableObject = mol->nextIntegrableObject(ii)) {
337  
338              force = integrableObject->getFrc();    
339 <            i->fz += force[whichDirection]
339 >            i->fz += force[whichDirection];
340          }
341           totalFZ_local += i->fz;
342      }
# Line 374 | Line 366 | void ZconstraintForceManager::doZconstraintForce(){
366          for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
367              integrableObject = mol->nextIntegrableObject(ii)) {
368  
369 <            force[whichDirection] = -getZFOfMovingMols(mol, integrableObject, totalFZ);
369 >            force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
370              integrableObject->addFrc(force);
371          }
372      }
# Line 386 | Line 378 | void ZconstraintForceManager::doZconstraintForce(){
378          for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
379              integrableObject = mol->nextIntegrableObject(ii)) {
380  
381 <            force[whichDirection] = -getZFOfMovingMols(mol, integrableObject, totalFZ);
381 >            force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
382              integrableObject->addFrc(force);
383          }
384      }
# Line 395 | Line 387 | void ZconstraintForceManager::doHarmonic(){
387  
388  
389   void ZconstraintForceManager::doHarmonic(){
398 –    double harmonicU;
399 –    double harmonicF;
400 –    double diff;
401 –    double totalFZ_local;
390      double totalFZ;
403 –
391      Vector3d force(0.0);
392      Vector3d com;
406 –
393      double totalFZ_local = 0;
408 –
394      std::list<ZconstraintMol>::iterator i;
395 +    StuntDouble* integrableObject;
396 +    Molecule::IntegrableObjectIterator ii;
397      Molecule* mol;
411 –     StuntDouble* integrableObject;
412 –    Molecule::IntegrableObjectIteraot ii;
413 –    
398      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
399          mol = i->mol;
400 <        com = mol->getCom();        
401 <        diff = com[whichDirection] - resPos;
402 <        harmonicU = 0.5 * kz * diff * diff;
403 <        currSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
404 <        harmonicF = -kz * diff;
400 >        com = mol->getCom();  
401 >        double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
402 >        double diff = com[whichDirection] - resPos;
403 >        double harmonicU = 0.5 * i->param.kz * diff * diff;
404 >        currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
405 >        double harmonicF = -i->param.kz * diff;
406          totalFZ_local += harmonicF;
407  
408          //adjust force
# Line 442 | Line 427 | void ZconstraintForceManager::doHarmonic(){
427          for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
428              integrableObject = mol->nextIntegrableObject(ii)) {
429  
430 <            force[whichDirection] = getHFOfUnconsMols(mol, harmonicF);
430 >            force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
431              integrableObject->addFrc(force);            
432          }
433      }
# Line 462 | Line 447 | bool ZconstraintForceManager::checkZConsState(){
447          com = i->mol->getCom();
448          diff = fabs(com[whichDirection] - i->param.zTargetPos);
449          if (diff > zconsTol_) {
465 –            //this fixed zconstraint molecule is about to moving
466 –            //moved this molecule to
467 –            j = i++;
468 –            newMovingZMols.push_back(fixedZMols_.erase(j));
450              if (usingZconsGap_) {
451                  i->endFixingTime = infiniteTime;
452              }
453 +            j = i++;
454 +            newMovingZMols.push_back(*j);
455 +            fixedZMols_.erase(j);
456 +            
457              changed_local = 1;            
458          }else {
459              ++i;
# Line 480 | Line 465 | bool ZconstraintForceManager::checkZConsState(){
465          com = i->mol->getCom();
466          diff = fabs(com[whichDirection] - i->param.zTargetPos);
467          if (diff <= zconsTol_) {
483 –            //this moving zconstraint molecule is about to fixed
484 –            //moved this molecule to
485 –            j = i++;
486 –            newFixedZMols.push_back(movingZMols_.erase(j));
468              if (usingZconsGap_) {
469                  i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
470              }
471 +            //this moving zconstraint molecule is about to fixed
472 +            //moved this molecule to
473 +            j = i++;
474 +            newFixedZMols.push_back(*j);
475 +            movingZMols_.erase(j);
476              changed_local = 1;            
477          }else {
478              ++i;
# Line 494 | Line 480 | bool ZconstraintForceManager::checkZConsState(){
480      }    
481  
482      //merge the lists
483 <    fixedZMols_.merge(newFixedZMols);
484 <    movingZMols_.merge(newFixedZMols);
483 >    fixedZMols_.insert(fixedZMols_.end(), newFixedZMols.begin(), newFixedZMols.end());
484 >    movingZMols_.insert(movingZMols_.end(), newMovingZMols.begin(), newMovingZMols.end());
485  
486      int changed;
487   #ifndef IS_MPI
# Line 550 | Line 536 | void ZconstraintForceManager::calcTotalMassMovingZMols
536      MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
537                  MPI_SUM, MPI_COMM_WORLD);
538   #else
539 <    totMassMovingZMols_ = massOfMovingZAtoms_local;
539 >    totMassMovingZMols_ = totMassMovingZMols_local;
540   #endif
541  
542   }
# Line 559 | Line 545 | double ZconstraintForceManager::getZFOfMovingMols(Stun
545    return totalForce * sd->getMass() / mol->getMass();
546   }
547  
548 < double ZconstraintForceManager::getZFOfMovingMols(StuntDouble*sd, double totalForce){
549 <  return totalForce * sd->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
548 > double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
549 >  return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
550   }
551  
552   double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
553    return totalForce * sd->getMass() / mol->getMass();
554   }
555  
556 < double ZconstraintForceManager::getHFOfUnconsMols(StuntDouble*sd, double totalForce){
557 <  return totalForce * sd->getMass() / totMassUnconsMols_;
556 > double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
557 >  return totalForce * mol->getMass() / totMassUnconsMols_;
558   }
559  
560   void ZconstraintForceManager::updateZPos(){
# Line 590 | Line 576 | double ZconstraintForceManager::getZTargetPos(int inde
576      double zTargetPos;
577   #ifndef IS_MPI    
578      Molecule* mol = info_->getMoleculeByGlobalIndex(index);
579 <    assert(!mol);
579 >    assert(mol);
580      Vector3d com = mol->getCom();
581      zTargetPos = com[whichDirection];
582   #else

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines