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

Comparing trunk/OOPSE-2.0/src/constraints/ZconstraintForceManager.cpp (file contents):
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2204 by gezelter, Fri Apr 15 22:04:00 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 46 | Line 46 | ZconstraintForceManager::ZconstraintForceManager(SimIn
46   #include "utils/OOPSEConstant.hpp"
47   #include "utils/StringUtils.hpp"
48   namespace oopse {
49 < ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) {
49 >  ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) {
50      currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot();
51      Globals* simParam = info_->getSimParams();
52  
53      if (simParam->haveDt()){
54 <        dt_ = simParam->getDt();
54 >      dt_ = simParam->getDt();
55      } else {
56 <        sprintf(painCave.errMsg,
57 <                "Integrator Error: dt is not set\n");
58 <        painCave.isFatal = 1;
59 <        simError();    
56 >      sprintf(painCave.errMsg,
57 >              "Integrator Error: dt is not set\n");
58 >      painCave.isFatal = 1;
59 >      simError();    
60      }
61  
62      if (simParam->haveZconstraintTime()){
63 <        zconsTime_ = simParam->getZconsTime();
63 >      zconsTime_ = simParam->getZconsTime();
64      }
65      else{
66 <        sprintf(painCave.errMsg,
67 <        "ZConstraint error: If you use a ZConstraint,\n"
68 <        "\tyou must set zconsTime.\n");
69 <        painCave.isFatal = 1;
70 <        simError();
66 >      sprintf(painCave.errMsg,
67 >              "ZConstraint error: If you use a ZConstraint,\n"
68 >              "\tyou must set zconsTime.\n");
69 >      painCave.isFatal = 1;
70 >      simError();
71      }
72  
73      if (simParam->haveZconsTol()){
74 <        zconsTol_ = simParam->getZconsTol();
74 >      zconsTol_ = simParam->getZconsTol();
75      }
76      else{
77 <        zconsTol_ = 0.01;
78 <        sprintf(painCave.errMsg,
79 <            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
80 <            "\tOOPSE will use a default value of %f.\n"
81 <            "\tTo set the tolerance, use the zconsTol variable.\n",
82 <            zconsTol_);
83 <        painCave.isFatal = 0;
84 <        simError();      
77 >      zconsTol_ = 0.01;
78 >      sprintf(painCave.errMsg,
79 >              "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
80 >              "\tOOPSE will use a default value of %f.\n"
81 >              "\tTo set the tolerance, use the zconsTol variable.\n",
82 >              zconsTol_);
83 >      painCave.isFatal = 0;
84 >      simError();      
85      }
86  
87      //set zcons gap
88      if (simParam->haveZConsGap()){
89 <        usingZconsGap_ = true;
90 <        zconsGap_ = simParam->getZconsGap();
89 >      usingZconsGap_ = true;
90 >      zconsGap_ = simParam->getZconsGap();
91      }else {
92 <        usingZconsGap_ = false;
93 <        zconsGap_ = 0.0;
92 >      usingZconsGap_ = false;
93 >      zconsGap_ = 0.0;
94      }
95  
96      //set zcons fixtime
97      if (simParam->haveZConsFixTime()){
98 <        zconsFixingTime_ = simParam->getZconsFixtime();
98 >      zconsFixingTime_ = simParam->getZconsFixtime();
99      } else {
100 <        zconsFixingTime_ = infiniteTime;
100 >      zconsFixingTime_ = infiniteTime;
101      }
102  
103      //set zconsUsingSMD
104      if (simParam->haveZConsUsingSMD()){
105 <        usingSMD_ = simParam->getZconsUsingSMD();
105 >      usingSMD_ = simParam->getZconsUsingSMD();
106      }else {
107 <        usingSMD_ =false;
107 >      usingSMD_ =false;
108      }
109      
110      zconsOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
# Line 114 | Line 114 | ZconstraintForceManager::ZconstraintForceManager(SimIn
114      double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;        
115      double targetTemp;
116      if (simParam->haveTargetTemp()) {
117 <        targetTemp = simParam->getTargetTemp();
117 >      targetTemp = simParam->getTargetTemp();
118      } else {
119 <        targetTemp = 298.0;
119 >      targetTemp = 298.0;
120      }
121      double zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox);
122          
# Line 125 | Line 125 | ZconstraintForceManager::ZconstraintForceManager(SimIn
125      //
126      for (int i = 0; i < nZconstraints; i++){
127  
128 <        ZconstraintParam param;
129 <        int zmolIndex = stamp[i]->getMolIndex();
130 <        if (stamp[i]->haveZpos()) {
131 <            param.zTargetPos = stamp[i]->getZpos();
132 <        } else {
133 <            param.zTargetPos = getZTargetPos(zmolIndex);
134 <        }
128 >      ZconstraintParam param;
129 >      int zmolIndex = stamp[i]->getMolIndex();
130 >      if (stamp[i]->haveZpos()) {
131 >        param.zTargetPos = stamp[i]->getZpos();
132 >      } else {
133 >        param.zTargetPos = getZTargetPos(zmolIndex);
134 >      }
135  
136 <        param.kz = zforceConstant * stamp[i]->getKratio();
136 >      param.kz = zforceConstant * stamp[i]->getKratio();
137  
138 <        if (stamp[i]->haveCantVel()) {
139 <            param.cantVel = stamp[i]->getCantVel();
140 <        } else {
141 <            param.cantVel = 0.0;
142 <        }
138 >      if (stamp[i]->haveCantVel()) {
139 >        param.cantVel = stamp[i]->getCantVel();
140 >      } else {
141 >        param.cantVel = 0.0;
142 >      }
143  
144 <        allZMolIndices_.insert(std::make_pair(zmolIndex, param));
144 >      allZMolIndices_.insert(std::make_pair(zmolIndex, param));
145      }
146  
147      //create fixedMols_, movingMols_ and unconsMols lists
# Line 151 | Line 151 | ZconstraintForceManager::ZconstraintForceManager(SimIn
151      double totMassUnconsMols_local = 0.0;    
152      std::vector<Molecule*>::iterator j;
153      for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
154 <        totMassUnconsMols_local += (*j)->getMass();
154 >      totMassUnconsMols_local += (*j)->getMass();
155      }    
156   #ifndef IS_MPI
157      totMassUnconsMols_ = totMassUnconsMols_local;
158   #else
159      MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
160 <        MPI_SUM, MPI_COMM_WORLD);  
160 >                  MPI_SUM, MPI_COMM_WORLD);  
161   #endif
162  
163      // creat zconsWriter  
# Line 169 | Line 169 | ZconstraintForceManager::ZconstraintForceManager(SimIn
169        simError();
170      }
171  
172 < }
172 >  }
173  
174 < ZconstraintForceManager::~ZconstraintForceManager(){
174 >  ZconstraintForceManager::~ZconstraintForceManager(){
175  
176 <  if (fzOut){
177 <    delete fzOut;
176 >    if (fzOut){
177 >      delete fzOut;
178 >    }
179 >
180    }
181  
182 < }
181 <
182 < void ZconstraintForceManager::update(){
182 >  void ZconstraintForceManager::update(){
183      fixedZMols_.clear();
184      movingZMols_.clear();
185      unzconsMols_.clear();
186  
187      for (std::map<int, ZconstraintParam>::iterator i = allZMolIndices_.begin(); i != allZMolIndices_.end(); ++i) {
188   #ifdef IS_MPI
189 <        if (info_->getMolToProc(i->first) == worldRank) {
189 >      if (info_->getMolToProc(i->first) == worldRank) {
190   #endif
191 <            ZconstraintMol zmol;
192 <            zmol.mol = info_->getMoleculeByGlobalIndex(i->first);
193 <            assert(zmol.mol);
194 <            zmol.param = i->second;
195 <            zmol.cantPos = zmol.param.zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
196 <            Vector3d com = zmol.mol->getCom();
197 <            double diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
198 <            if (diff < zconsTol_) {
199 <                fixedZMols_.push_back(zmol);
200 <            } else {
201 <                movingZMols_.push_back(zmol);            
202 <            }
191 >        ZconstraintMol zmol;
192 >        zmol.mol = info_->getMoleculeByGlobalIndex(i->first);
193 >        assert(zmol.mol);
194 >        zmol.param = i->second;
195 >        zmol.cantPos = zmol.param.zTargetPos; /**@todo fixed me when zmol migrate, it is incorrect*/
196 >        Vector3d com = zmol.mol->getCom();
197 >        double diff = fabs(zmol.param.zTargetPos - com[whichDirection]);
198 >        if (diff < zconsTol_) {
199 >          fixedZMols_.push_back(zmol);
200 >        } else {
201 >          movingZMols_.push_back(zmol);            
202 >        }
203  
204   #ifdef IS_MPI
205 <        }
205 >      }
206   #endif
207      }
208  
# Line 210 | Line 210 | void ZconstraintForceManager::update(){
210  
211      std::set<int> zmolSet;
212      for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
213 <        zmolSet.insert(i->mol->getGlobalIndex());
213 >      zmolSet.insert(i->mol->getGlobalIndex());
214      }    
215  
216      for (std::list<ZconstraintMol>::iterator i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
217 <        zmolSet.insert(i->mol->getGlobalIndex());
217 >      zmolSet.insert(i->mol->getGlobalIndex());
218      }
219  
220      SimInfo::MoleculeIterator mi;
221      Molecule* mol;
222      for(mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
223 <        if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
224 <            unzconsMols_.push_back(mol);
225 <        }
223 >      if (zmolSet.find(mol->getGlobalIndex()) == zmolSet.end()) {
224 >        unzconsMols_.push_back(mol);
225 >      }
226      }
227  
228 < }
228 >  }
229  
230 < bool ZconstraintForceManager::isZMol(Molecule* mol){
230 >  bool ZconstraintForceManager::isZMol(Molecule* mol){
231      return allZMolIndices_.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
232 < }
232 >  }
233  
234 < void ZconstraintForceManager::init(){
234 >  void ZconstraintForceManager::init(){
235  
236 <  //zero out the velocities of center of mass of unconstrained molecules
237 <  //and the velocities of center of mass of every single z-constrained molecueles
238 <  zeroVelocity();
236 >    //zero out the velocities of center of mass of unconstrained molecules
237 >    //and the velocities of center of mass of every single z-constrained molecueles
238 >    zeroVelocity();
239  
240 <  currZconsTime_ = currSnapshot_->getTime();
241 < }
240 >    currZconsTime_ = currSnapshot_->getTime();
241 >  }
242  
243 < void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
243 >  void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
244      ForceManager::calcForces(needPotential, needStress);
245      
246      if (usingZconsGap_){
247 <        updateZPos();
247 >      updateZPos();
248      }
249  
250      if (checkZConsState()){
251 <        zeroVelocity();    
252 <        calcTotalMassMovingZMols();
251 >      zeroVelocity();    
252 >      calcTotalMassMovingZMols();
253      }  
254  
255      //do zconstraint force;
256      if (haveFixedZMols()){
257 <        doZconstraintForce();
257 >      doZconstraintForce();
258      }
259  
260      //use external force to move the molecules to the specified positions
261      if (haveMovingZMols()){
262 <        doHarmonic();
262 >      doHarmonic();
263      }
264  
265      //write out forces and current positions of z-constraint molecules    
266      if (currSnapshot_->getTime() >= currZconsTime_){
267 <        std::list<ZconstraintMol>::iterator i;
268 <        Vector3d com;
269 <        for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
270 <            com = i->mol->getCom();
271 <            i->zpos = com[whichDirection];
272 <        }
267 >      std::list<ZconstraintMol>::iterator i;
268 >      Vector3d com;
269 >      for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
270 >        com = i->mol->getCom();
271 >        i->zpos = com[whichDirection];
272 >      }
273          
274 <        fzOut->writeFZ(fixedZMols_);
275 <        currZconsTime_ += zconsTime_;
274 >      fzOut->writeFZ(fixedZMols_);
275 >      currZconsTime_ += zconsTime_;
276      }
277 < }
277 >  }
278  
279 < void ZconstraintForceManager::zeroVelocity(){
279 >  void ZconstraintForceManager::zeroVelocity(){
280  
281      Vector3d comVel;
282      Vector3d vel;
# Line 287 | Line 287 | void ZconstraintForceManager::zeroVelocity(){
287  
288      //zero out the velocities of center of mass of fixed z-constrained molecules
289      for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
290 <        mol = i->mol;
291 <        comVel = mol->getComVel();
292 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
293 <            integrableObject = mol->nextIntegrableObject(ii)) {
294 <            vel = integrableObject->getVel();  
295 <            vel[whichDirection] -= comVel[whichDirection];
296 <            integrableObject->setVel(vel);
297 <        }
290 >      mol = i->mol;
291 >      comVel = mol->getComVel();
292 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
293 >          integrableObject = mol->nextIntegrableObject(ii)) {
294 >        vel = integrableObject->getVel();  
295 >        vel[whichDirection] -= comVel[whichDirection];
296 >        integrableObject->setVel(vel);
297 >      }
298      }
299  
300      // calculate the vz of center of mass of moving molecules(include unconstrained molecules
# Line 303 | Line 303 | void ZconstraintForceManager::zeroVelocity(){
303      double pzMovingMols;
304      
305      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
306 <        mol = i->mol;        
307 <        comVel = mol->getComVel();
308 <        pzMovingMols_local +=  mol->getMass() * comVel[whichDirection];  
306 >      mol = i->mol;        
307 >      comVel = mol->getComVel();
308 >      pzMovingMols_local +=  mol->getMass() * comVel[whichDirection];  
309      }
310  
311      std::vector<Molecule*>::iterator j;
312      for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
313 <        mol =*j;
314 <        comVel = mol->getComVel();
315 <        pzMovingMols_local += mol->getMass() * comVel[whichDirection];
313 >      mol =*j;
314 >      comVel = mol->getComVel();
315 >      pzMovingMols_local += mol->getMass() * comVel[whichDirection];
316      }
317      
318   #ifndef IS_MPI
319      pzMovingMols = pzMovingMols_local;
320   #else
321      MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
322 <        MPI_SUM, MPI_COMM_WORLD);
322 >                  MPI_SUM, MPI_COMM_WORLD);
323   #endif
324  
325      double vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_);
326  
327      //modify the velocities of moving z-constrained molecuels
328      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
329 <        mol = i->mol;
330 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
331 <            integrableObject = mol->nextIntegrableObject(ii)) {
329 >      mol = i->mol;
330 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
331 >          integrableObject = mol->nextIntegrableObject(ii)) {
332  
333 <            vel = integrableObject->getVel();
334 <            vel[whichDirection] -= vzMovingMols;
335 <            integrableObject->setVel(vel);
336 <        }
333 >        vel = integrableObject->getVel();
334 >        vel[whichDirection] -= vzMovingMols;
335 >        integrableObject->setVel(vel);
336 >      }
337      }
338  
339      //modify the velocites of unconstrained molecules  
340      for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
341 <        mol =*j;
342 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
343 <            integrableObject = mol->nextIntegrableObject(ii)) {
341 >      mol =*j;
342 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
343 >          integrableObject = mol->nextIntegrableObject(ii)) {
344  
345 <            vel = integrableObject->getVel();
346 <            vel[whichDirection] -= vzMovingMols;
347 <            integrableObject->setVel(vel);
348 <        }
345 >        vel = integrableObject->getVel();
346 >        vel[whichDirection] -= vzMovingMols;
347 >        integrableObject->setVel(vel);
348 >      }
349      }
350      
351 < }
351 >  }
352  
353  
354 < void ZconstraintForceManager::doZconstraintForce(){
355 <  double totalFZ;
356 <  double totalFZ_local;
357 <  Vector3d com;
358 <  Vector3d force(0.0);
354 >  void ZconstraintForceManager::doZconstraintForce(){
355 >    double totalFZ;
356 >    double totalFZ_local;
357 >    Vector3d com;
358 >    Vector3d force(0.0);
359  
360 <  //constrain the molecules which do not reach the specified positions  
360 >    //constrain the molecules which do not reach the specified positions  
361  
362 <  //Zero Out the force of z-contrained molecules    
363 <  totalFZ_local = 0;
362 >    //Zero Out the force of z-contrained molecules    
363 >    totalFZ_local = 0;
364  
365  
366      //calculate the total z-contrained force of fixed z-contrained molecules
# Line 370 | Line 370 | void ZconstraintForceManager::doZconstraintForce(){
370      Molecule::IntegrableObjectIterator ii;
371  
372      for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
373 <        mol = i->mol;
374 <        i->fz = 0.0;
375 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
376 <            integrableObject = mol->nextIntegrableObject(ii)) {
373 >      mol = i->mol;
374 >      i->fz = 0.0;
375 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
376 >          integrableObject = mol->nextIntegrableObject(ii)) {
377  
378 <            force = integrableObject->getFrc();    
379 <            i->fz += force[whichDirection];
380 <        }
381 <         totalFZ_local += i->fz;
378 >        force = integrableObject->getFrc();    
379 >        i->fz += force[whichDirection];
380 >      }
381 >      totalFZ_local += i->fz;
382      }
383  
384 <  //calculate total z-constraint force
384 >    //calculate total z-constraint force
385   #ifdef IS_MPI
386 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
386 >    MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
387   #else
388 <  totalFZ = totalFZ_local;
388 >    totalFZ = totalFZ_local;
389   #endif
390  
391  
392 <     // apply negative to fixed z-constrained molecues;
392 >    // apply negative to fixed z-constrained molecues;
393      for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
394 <        mol = i->mol;
395 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
396 <            integrableObject = mol->nextIntegrableObject(ii)) {
394 >      mol = i->mol;
395 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
396 >          integrableObject = mol->nextIntegrableObject(ii)) {
397  
398 <            force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
399 <            integrableObject->addFrc(force);
400 <        }
398 >        force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
399 >        integrableObject->addFrc(force);
400 >      }
401      }
402  
403 <  //modify the forces of moving z-constrained molecules
403 >    //modify the forces of moving z-constrained molecules
404      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
405 <        mol = i->mol;
406 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
407 <            integrableObject = mol->nextIntegrableObject(ii)) {
405 >      mol = i->mol;
406 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
407 >          integrableObject = mol->nextIntegrableObject(ii)) {
408  
409 <            force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
410 <            integrableObject->addFrc(force);
411 <        }
409 >        force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
410 >        integrableObject->addFrc(force);
411 >      }
412      }
413  
414 <  //modify the forces of unconstrained molecules
414 >    //modify the forces of unconstrained molecules
415      std::vector<Molecule*>::iterator j;
416      for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
417 <        mol =*j;
418 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419 <            integrableObject = mol->nextIntegrableObject(ii)) {
417 >      mol =*j;
418 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
419 >          integrableObject = mol->nextIntegrableObject(ii)) {
420  
421 <            force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
422 <            integrableObject->addFrc(force);
423 <        }
421 >        force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
422 >        integrableObject->addFrc(force);
423 >      }
424      }
425  
426 < }
426 >  }
427  
428  
429 < void ZconstraintForceManager::doHarmonic(){
429 >  void ZconstraintForceManager::doHarmonic(){
430      double totalFZ;
431      Vector3d force(0.0);
432      Vector3d com;
# Line 436 | Line 436 | void ZconstraintForceManager::doHarmonic(){
436      Molecule::IntegrableObjectIterator ii;
437      Molecule* mol;
438      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
439 <        mol = i->mol;
440 <        com = mol->getCom();  
441 <        double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
442 <        double diff = com[whichDirection] - resPos;
443 <        double harmonicU = 0.5 * i->param.kz * diff * diff;
444 <        currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
445 <        double harmonicF = -i->param.kz * diff;
446 <        totalFZ_local += harmonicF;
439 >      mol = i->mol;
440 >      com = mol->getCom();  
441 >      double resPos = usingSMD_? i->cantPos : i->param.zTargetPos;
442 >      double diff = com[whichDirection] - resPos;
443 >      double harmonicU = 0.5 * i->param.kz * diff * diff;
444 >      currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU;
445 >      double harmonicF = -i->param.kz * diff;
446 >      totalFZ_local += harmonicF;
447  
448 <        //adjust force
449 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
450 <            integrableObject = mol->nextIntegrableObject(ii)) {
448 >      //adjust force
449 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
450 >          integrableObject = mol->nextIntegrableObject(ii)) {
451  
452 <            force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
453 <            integrableObject->addFrc(force);            
454 <        }
452 >        force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
453 >        integrableObject->addFrc(force);            
454 >      }
455      }
456  
457   #ifndef IS_MPI
458 <  totalFZ = totalFZ_local;
458 >    totalFZ = totalFZ_local;
459   #else
460 <  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
460 >    MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
461   #endif
462  
463      //modify the forces of unconstrained molecules
464      std::vector<Molecule*>::iterator j;
465      for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) {
466 <        mol = *j;
467 <        for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
468 <            integrableObject = mol->nextIntegrableObject(ii)) {
466 >      mol = *j;
467 >      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
468 >          integrableObject = mol->nextIntegrableObject(ii)) {
469  
470 <            force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
471 <            integrableObject->addFrc(force);            
472 <        }
470 >        force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
471 >        integrableObject->addFrc(force);            
472 >      }
473      }
474  
475 < }
475 >  }
476  
477 < bool ZconstraintForceManager::checkZConsState(){
477 >  bool ZconstraintForceManager::checkZConsState(){
478      Vector3d com;
479      double diff;
480      int changed_local = 0;
# Line 484 | Line 484 | bool ZconstraintForceManager::checkZConsState(){
484      
485      std::list<ZconstraintMol> newMovingZMols;
486      for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end();) {
487 <        com = i->mol->getCom();
488 <        diff = fabs(com[whichDirection] - i->param.zTargetPos);
489 <        if (diff > zconsTol_) {
490 <            if (usingZconsGap_) {
491 <                i->endFixingTime = infiniteTime;
492 <            }
493 <            j = i++;
494 <            newMovingZMols.push_back(*j);
495 <            fixedZMols_.erase(j);
487 >      com = i->mol->getCom();
488 >      diff = fabs(com[whichDirection] - i->param.zTargetPos);
489 >      if (diff > zconsTol_) {
490 >        if (usingZconsGap_) {
491 >          i->endFixingTime = infiniteTime;
492 >        }
493 >        j = i++;
494 >        newMovingZMols.push_back(*j);
495 >        fixedZMols_.erase(j);
496              
497 <            changed_local = 1;            
498 <        }else {
499 <            ++i;
500 <        }
497 >        changed_local = 1;            
498 >      }else {
499 >        ++i;
500 >      }
501      }  
502  
503      std::list<ZconstraintMol> newFixedZMols;
504      for ( i = movingZMols_.begin(); i !=  movingZMols_.end();) {
505 <        com = i->mol->getCom();
506 <        diff = fabs(com[whichDirection] - i->param.zTargetPos);
507 <        if (diff <= zconsTol_) {
508 <            if (usingZconsGap_) {
509 <                i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
510 <            }
511 <            //this moving zconstraint molecule is about to fixed
512 <            //moved this molecule to
513 <            j = i++;
514 <            newFixedZMols.push_back(*j);
515 <            movingZMols_.erase(j);
516 <            changed_local = 1;            
517 <        }else {
518 <            ++i;
519 <        }
505 >      com = i->mol->getCom();
506 >      diff = fabs(com[whichDirection] - i->param.zTargetPos);
507 >      if (diff <= zconsTol_) {
508 >        if (usingZconsGap_) {
509 >          i->endFixingTime = currSnapshot_->getTime() + zconsFixingTime_;
510 >        }
511 >        //this moving zconstraint molecule is about to fixed
512 >        //moved this molecule to
513 >        j = i++;
514 >        newFixedZMols.push_back(*j);
515 >        movingZMols_.erase(j);
516 >        changed_local = 1;            
517 >      }else {
518 >        ++i;
519 >      }
520      }    
521  
522      //merge the lists
# Line 531 | Line 531 | bool ZconstraintForceManager::checkZConsState(){
531   #endif
532  
533      return (changed > 0);
534 < }
534 >  }
535  
536 < bool ZconstraintForceManager::haveFixedZMols(){
537 <  int havingFixed;
538 <  int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
536 >  bool ZconstraintForceManager::haveFixedZMols(){
537 >    int havingFixed;
538 >    int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
539  
540   #ifndef IS_MPI
541 <  havingFixed = havingFixed_local;
541 >    havingFixed = havingFixed_local;
542   #else
543 <  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
544 <                MPI_COMM_WORLD);
543 >    MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
544 >                  MPI_COMM_WORLD);
545   #endif
546  
547 <  return havingFixed > 0;
548 < }
547 >    return havingFixed > 0;
548 >  }
549  
550  
551 < bool ZconstraintForceManager::haveMovingZMols(){
552 <  int havingMoving_local;
553 <  int havingMoving;
551 >  bool ZconstraintForceManager::haveMovingZMols(){
552 >    int havingMoving_local;
553 >    int havingMoving;
554  
555 <  havingMoving_local = movingZMols_.empty()? 0 : 1;
555 >    havingMoving_local = movingZMols_.empty()? 0 : 1;
556  
557   #ifndef IS_MPI
558 <  havingMoving = havingMoving_local;
558 >    havingMoving = havingMoving_local;
559   #else
560 <  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
561 <                MPI_COMM_WORLD);
560 >    MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
561 >                  MPI_COMM_WORLD);
562   #endif
563  
564 <  return havingMoving > 0;
565 < }
564 >    return havingMoving > 0;
565 >  }
566  
567 < void ZconstraintForceManager::calcTotalMassMovingZMols(){
567 >  void ZconstraintForceManager::calcTotalMassMovingZMols(){
568  
569      double totMassMovingZMols_local = 0.0;
570      std::list<ZconstraintMol>::iterator i;
571      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
572 <        totMassMovingZMols_local += i->mol->getMass();
572 >      totMassMovingZMols_local += i->mol->getMass();
573      }
574      
575   #ifdef IS_MPI
576      MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
577 <                MPI_SUM, MPI_COMM_WORLD);
577 >                  MPI_SUM, MPI_COMM_WORLD);
578   #else
579      totMassMovingZMols_ = totMassMovingZMols_local;
580   #endif
581  
582 < }
582 >  }
583  
584 < double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
585 <  return totalForce * sd->getMass() / mol->getMass();
586 < }
584 >  double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
585 >    return totalForce * sd->getMass() / mol->getMass();
586 >  }
587  
588 < double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
589 <  return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
590 < }
588 >  double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){
589 >    return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_);
590 >  }
591  
592 < double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
593 <  return totalForce * sd->getMass() / mol->getMass();
594 < }
592 >  double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){
593 >    return totalForce * sd->getMass() / mol->getMass();
594 >  }
595  
596 < double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
597 <  return totalForce * mol->getMass() / totMassUnconsMols_;
598 < }
596 >  double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){
597 >    return totalForce * mol->getMass() / totMassUnconsMols_;
598 >  }
599  
600 < void ZconstraintForceManager::updateZPos(){
600 >  void ZconstraintForceManager::updateZPos(){
601      double curTime = currSnapshot_->getTime();
602      std::list<ZconstraintMol>::iterator i;
603      for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) {
604 <        i->param.zTargetPos += zconsGap_;    
604 >      i->param.zTargetPos += zconsGap_;    
605      }  
606 < }
606 >  }
607  
608 < void ZconstraintForceManager::updateCantPos(){
608 >  void ZconstraintForceManager::updateCantPos(){
609      std::list<ZconstraintMol>::iterator i;
610      for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) {
611 <        i->cantPos += i->param.cantVel * dt_;
611 >      i->cantPos += i->param.cantVel * dt_;
612      }
613 < }
613 >  }
614  
615 < double ZconstraintForceManager::getZTargetPos(int index){
615 >  double ZconstraintForceManager::getZTargetPos(int index){
616      double zTargetPos;
617   #ifndef IS_MPI    
618      Molecule* mol = info_->getMoleculeByGlobalIndex(index);
# Line 624 | Line 624 | double ZconstraintForceManager::getZTargetPos(int inde
624      MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
625   #endif
626      return zTargetPos;
627 < }
627 >  }
628  
629   }

Diff Legend

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