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