1 |
#include <cmath>
|
2 |
#include "constraints/ZconstraintForceManager.hpp"
|
3 |
#include "integrators/Integrator.hpp"
|
4 |
#include "utils/simError.h"
|
5 |
#include "utils/OOPSEConstant.hpp"
|
6 |
#include "utils/StringUtils.hpp"
|
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 |
usingZconsGap_ = true;
|
49 |
zconsGap_ = simParam->getZconsGap();
|
50 |
}else {
|
51 |
usingZconsGap_ = false;
|
52 |
zconsGap_ = 0.0;
|
53 |
}
|
54 |
|
55 |
//set zcons fixtime
|
56 |
if (simParam->haveZConsFixTime()){
|
57 |
zconsFixingTime_ = simParam->getZconsFixtime();
|
58 |
} else {
|
59 |
zconsFixingTime_ = infiniteTime;
|
60 |
}
|
61 |
|
62 |
//set zconsUsingSMD
|
63 |
if (simParam->haveZConsUsingSMD()){
|
64 |
usingSMD_ = simParam->getZconsUsingSMD();
|
65 |
}else {
|
66 |
usingSMD_ =false;
|
67 |
}
|
68 |
|
69 |
zconsOutput_ = getPrefix(info_->getFinalConfigFileName()) + ".fz";
|
70 |
|
71 |
//estimate the force constant of harmonical potential
|
72 |
Mat3x3d hmat = currSnapshot_->getHmat();
|
73 |
double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2;
|
74 |
double targetTemp;
|
75 |
if (simParam->haveTargetTemp()) {
|
76 |
targetTemp = simParam->getTargetTemp();
|
77 |
} else {
|
78 |
targetTemp = 298.0;
|
79 |
}
|
80 |
double zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox);
|
81 |
|
82 |
int nZconstraints = simParam->getNzConstraints();
|
83 |
ZconStamp** stamp = simParam->getZconStamp();
|
84 |
//
|
85 |
for (int i = 0; i < nZconstraints; i++){
|
86 |
|
87 |
ZconstraintParam param;
|
88 |
int zmolIndex = stamp[i]->getMolIndex();
|
89 |
if (stamp[i]->haveZpos()) {
|
90 |
param.zTargetPos = stamp[i]->getZpos();
|
91 |
} else {
|
92 |
param.zTargetPos = getZTargetPos(zmolIndex);
|
93 |
}
|
94 |
|
95 |
param.kz = zforceConstant * stamp[i]->getKratio();
|
96 |
|
97 |
if (stamp[i]->haveCantVel()) {
|
98 |
param.cantVel = stamp[i]->getCantVel();
|
99 |
} else {
|
100 |
param.cantVel = 0.0;
|
101 |
}
|
102 |
|
103 |
allZMolIndices_.insert(std::make_pair(zmolIndex, param));
|
104 |
}
|
105 |
|
106 |
//create fixedMols_, movingMols_ and unconsMols lists
|
107 |
update();
|
108 |
|
109 |
//calculate masss of unconstraint molecules in the whole system (never change during the simulation)
|
110 |
double totMassUnconsMols_local = 0.0;
|
111 |
std::vector<Molecule*>::iterator j;
|
112 |
for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
|
113 |
totMassUnconsMols_local += (*j)->getMass();
|
114 |
}
|
115 |
#ifndef IS_MPI
|
116 |
totMassUnconsMols_ = totMassUnconsMols_local;
|
117 |
#else
|
118 |
MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE,
|
119 |
MPI_SUM, MPI_COMM_WORLD);
|
120 |
#endif
|
121 |
|
122 |
// creat zconsWriter
|
123 |
fzOut = new ZConsWriter(info_, zconsOutput_.c_str());
|
124 |
|
125 |
if (!fzOut){
|
126 |
sprintf(painCave.errMsg, "Fail to create ZConsWriter\n");
|
127 |
painCave.isFatal = 1;
|
128 |
simError();
|
129 |
}
|
130 |
|
131 |
}
|
132 |
|
133 |
ZconstraintForceManager::~ZconstraintForceManager(){
|
134 |
|
135 |
if (fzOut){
|
136 |
delete fzOut;
|
137 |
}
|
138 |
|
139 |
}
|
140 |
|
141 |
void ZconstraintForceManager::update(){
|
142 |
fixedZMols_.clear();
|
143 |
movingZMols_.clear();
|
144 |
unzconsMols_.clear();
|
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 = 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 {
|
160 |
movingZMols_.push_back(zmol);
|
161 |
}
|
162 |
|
163 |
#ifdef IS_MPI
|
164 |
}
|
165 |
#endif
|
166 |
}
|
167 |
|
168 |
calcTotalMassMovingZMols();
|
169 |
|
170 |
std::set<int> zmolSet;
|
171 |
for (std::list<ZconstraintMol>::iterator i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
|
172 |
zmolSet.insert(i->mol->getGlobalIndex());
|
173 |
}
|
174 |
|
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(mol);
|
184 |
}
|
185 |
}
|
186 |
|
187 |
}
|
188 |
|
189 |
bool ZconstraintForceManager::isZMol(Molecule* mol){
|
190 |
return allZMolIndices_.find(mol->getGlobalIndex()) == allZMolIndices_.end() ? false : true;
|
191 |
}
|
192 |
|
193 |
void ZconstraintForceManager::init(){
|
194 |
|
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 |
currZconsTime_ = currSnapshot_->getTime();
|
200 |
}
|
201 |
|
202 |
void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){
|
203 |
ForceManager::calcForces(needPotential, needStress);
|
204 |
|
205 |
if (usingZconsGap_){
|
206 |
updateZPos();
|
207 |
}
|
208 |
|
209 |
if (checkZConsState()){
|
210 |
zeroVelocity();
|
211 |
calcTotalMassMovingZMols();
|
212 |
}
|
213 |
|
214 |
//do zconstraint force;
|
215 |
if (haveFixedZMols()){
|
216 |
doZconstraintForce();
|
217 |
}
|
218 |
|
219 |
//use external force to move the molecules to the specified positions
|
220 |
if (haveMovingZMols()){
|
221 |
doHarmonic();
|
222 |
}
|
223 |
|
224 |
//write out forces and current positions of z-constraint molecules
|
225 |
if (currSnapshot_->getTime() >= currZconsTime_){
|
226 |
std::list<ZconstraintMol>::iterator i;
|
227 |
Vector3d com;
|
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 |
currZconsTime_ += zconsTime_;
|
235 |
}
|
236 |
}
|
237 |
|
238 |
void ZconstraintForceManager::zeroVelocity(){
|
239 |
|
240 |
Vector3d comVel;
|
241 |
Vector3d vel;
|
242 |
std::list<ZconstraintMol>::iterator i;
|
243 |
Molecule* mol;
|
244 |
StuntDouble* integrableObject;
|
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) {
|
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 |
}
|
257 |
}
|
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 = 0.0;
|
262 |
double pzMovingMols;
|
263 |
|
264 |
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
|
265 |
mol = i->mol;
|
266 |
comVel = mol->getComVel();
|
267 |
pzMovingMols_local += mol->getMass() * comVel[whichDirection];
|
268 |
}
|
269 |
|
270 |
std::vector<Molecule*>::iterator j;
|
271 |
for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
|
272 |
mol =*j;
|
273 |
comVel = mol->getComVel();
|
274 |
pzMovingMols_local += mol->getMass() * comVel[whichDirection];
|
275 |
}
|
276 |
|
277 |
#ifndef IS_MPI
|
278 |
pzMovingMols = pzMovingMols_local;
|
279 |
#else
|
280 |
MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE,
|
281 |
MPI_SUM, MPI_COMM_WORLD);
|
282 |
#endif
|
283 |
|
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) {
|
288 |
mol = i->mol;
|
289 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
290 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
291 |
|
292 |
vel = integrableObject->getVel();
|
293 |
vel[whichDirection] -= vzMovingMols;
|
294 |
integrableObject->setVel(vel);
|
295 |
}
|
296 |
}
|
297 |
|
298 |
//modify the velocites of unconstrained molecules
|
299 |
for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
|
300 |
mol =*j;
|
301 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
302 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
303 |
|
304 |
vel = integrableObject->getVel();
|
305 |
vel[whichDirection] -= vzMovingMols;
|
306 |
integrableObject->setVel(vel);
|
307 |
}
|
308 |
}
|
309 |
|
310 |
}
|
311 |
|
312 |
|
313 |
void ZconstraintForceManager::doZconstraintForce(){
|
314 |
Atom** zconsAtoms;
|
315 |
double totalFZ;
|
316 |
double totalFZ_local;
|
317 |
Vector3d com;
|
318 |
Vector3d force(0.0);
|
319 |
|
320 |
//constrain the molecules which do not reach the specified positions
|
321 |
|
322 |
//Zero Out the force of z-contrained molecules
|
323 |
totalFZ_local = 0;
|
324 |
|
325 |
|
326 |
//calculate the total z-contrained force of fixed z-contrained molecules
|
327 |
std::list<ZconstraintMol>::iterator i;
|
328 |
Molecule* mol;
|
329 |
StuntDouble* integrableObject;
|
330 |
Molecule::IntegrableObjectIterator ii;
|
331 |
|
332 |
for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
|
333 |
mol = i->mol;
|
334 |
i->fz = 0.0;
|
335 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
336 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
337 |
|
338 |
force = integrableObject->getFrc();
|
339 |
i->fz += force[whichDirection];
|
340 |
}
|
341 |
totalFZ_local += i->fz;
|
342 |
}
|
343 |
|
344 |
//calculate total z-constraint force
|
345 |
#ifdef IS_MPI
|
346 |
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
347 |
#else
|
348 |
totalFZ = totalFZ_local;
|
349 |
#endif
|
350 |
|
351 |
|
352 |
// apply negative to fixed z-constrained molecues;
|
353 |
for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
|
354 |
mol = i->mol;
|
355 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
356 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
357 |
|
358 |
force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz);
|
359 |
integrableObject->addFrc(force);
|
360 |
}
|
361 |
}
|
362 |
|
363 |
//modify the forces of moving z-constrained molecules
|
364 |
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
|
365 |
mol = i->mol;
|
366 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
367 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
368 |
|
369 |
force[whichDirection] = -getZFOfMovingMols(mol,totalFZ);
|
370 |
integrableObject->addFrc(force);
|
371 |
}
|
372 |
}
|
373 |
|
374 |
//modify the forces of unconstrained molecules
|
375 |
std::vector<Molecule*>::iterator j;
|
376 |
for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
|
377 |
mol =*j;
|
378 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
379 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
380 |
|
381 |
force[whichDirection] = -getZFOfMovingMols(mol, totalFZ);
|
382 |
integrableObject->addFrc(force);
|
383 |
}
|
384 |
}
|
385 |
|
386 |
}
|
387 |
|
388 |
|
389 |
void ZconstraintForceManager::doHarmonic(){
|
390 |
double totalFZ;
|
391 |
Vector3d force(0.0);
|
392 |
Vector3d com;
|
393 |
double totalFZ_local = 0;
|
394 |
std::list<ZconstraintMol>::iterator i;
|
395 |
StuntDouble* integrableObject;
|
396 |
Molecule::IntegrableObjectIterator ii;
|
397 |
Molecule* mol;
|
398 |
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
|
399 |
mol = i->mol;
|
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
|
409 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
410 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
411 |
|
412 |
force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF);
|
413 |
integrableObject->addFrc(force);
|
414 |
}
|
415 |
}
|
416 |
|
417 |
#ifndef IS_MPI
|
418 |
totalFZ = totalFZ_local;
|
419 |
#else
|
420 |
MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
|
421 |
#endif
|
422 |
|
423 |
//modify the forces of unconstrained molecules
|
424 |
std::vector<Molecule*>::iterator j;
|
425 |
for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) {
|
426 |
mol = *j;
|
427 |
for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;
|
428 |
integrableObject = mol->nextIntegrableObject(ii)) {
|
429 |
|
430 |
force[whichDirection] = getHFOfUnconsMols(mol, totalFZ);
|
431 |
integrableObject->addFrc(force);
|
432 |
}
|
433 |
}
|
434 |
|
435 |
}
|
436 |
|
437 |
bool ZconstraintForceManager::checkZConsState(){
|
438 |
Vector3d com;
|
439 |
double diff;
|
440 |
int changed_local = 0;
|
441 |
|
442 |
std::list<ZconstraintMol>::iterator i;
|
443 |
std::list<ZconstraintMol>::iterator j;
|
444 |
|
445 |
std::list<ZconstraintMol> newMovingZMols;
|
446 |
for ( i = fixedZMols_.begin(); i != fixedZMols_.end();) {
|
447 |
com = i->mol->getCom();
|
448 |
diff = fabs(com[whichDirection] - i->param.zTargetPos);
|
449 |
if (diff > zconsTol_) {
|
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;
|
460 |
}
|
461 |
}
|
462 |
|
463 |
std::list<ZconstraintMol> newFixedZMols;
|
464 |
for ( i = movingZMols_.begin(); i != movingZMols_.end();) {
|
465 |
com = i->mol->getCom();
|
466 |
diff = fabs(com[whichDirection] - i->param.zTargetPos);
|
467 |
if (diff <= zconsTol_) {
|
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;
|
479 |
}
|
480 |
}
|
481 |
|
482 |
//merge the lists
|
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
|
488 |
changed = changed_local;
|
489 |
#else
|
490 |
MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
|
491 |
#endif
|
492 |
|
493 |
return (changed > 0);
|
494 |
}
|
495 |
|
496 |
bool ZconstraintForceManager::haveFixedZMols(){
|
497 |
int havingFixed;
|
498 |
int havingFixed_local = fixedZMols_.empty() ? 0 : 1;
|
499 |
|
500 |
#ifndef IS_MPI
|
501 |
havingFixed = havingFixed_local;
|
502 |
#else
|
503 |
MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
|
504 |
MPI_COMM_WORLD);
|
505 |
#endif
|
506 |
|
507 |
return havingFixed > 0;
|
508 |
}
|
509 |
|
510 |
|
511 |
bool ZconstraintForceManager::haveMovingZMols(){
|
512 |
int havingMoving_local;
|
513 |
int havingMoving;
|
514 |
|
515 |
havingMoving_local = movingZMols_.empty()? 0 : 1;
|
516 |
|
517 |
#ifndef IS_MPI
|
518 |
havingMoving = havingMoving_local;
|
519 |
#else
|
520 |
MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
|
521 |
MPI_COMM_WORLD);
|
522 |
#endif
|
523 |
|
524 |
return havingMoving > 0;
|
525 |
}
|
526 |
|
527 |
void ZconstraintForceManager::calcTotalMassMovingZMols(){
|
528 |
|
529 |
double totMassMovingZMols_local = 0.0;
|
530 |
std::list<ZconstraintMol>::iterator i;
|
531 |
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
|
532 |
totMassMovingZMols_local += i->mol->getMass();
|
533 |
}
|
534 |
|
535 |
#ifdef IS_MPI
|
536 |
MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE,
|
537 |
MPI_SUM, MPI_COMM_WORLD);
|
538 |
#else
|
539 |
totMassMovingZMols_ = totMassMovingZMols_local;
|
540 |
#endif
|
541 |
|
542 |
}
|
543 |
|
544 |
double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){
|
545 |
return totalForce * sd->getMass() / mol->getMass();
|
546 |
}
|
547 |
|
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(Molecule* mol, double totalForce){
|
557 |
return totalForce * mol->getMass() / totMassUnconsMols_;
|
558 |
}
|
559 |
|
560 |
void ZconstraintForceManager::updateZPos(){
|
561 |
double curTime = currSnapshot_->getTime();
|
562 |
std::list<ZconstraintMol>::iterator i;
|
563 |
for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) {
|
564 |
i->param.zTargetPos += zconsGap_;
|
565 |
}
|
566 |
}
|
567 |
|
568 |
void ZconstraintForceManager::updateCantPos(){
|
569 |
std::list<ZconstraintMol>::iterator i;
|
570 |
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) {
|
571 |
i->cantPos += i->param.cantVel * dt_;
|
572 |
}
|
573 |
}
|
574 |
|
575 |
double ZconstraintForceManager::getZTargetPos(int index){
|
576 |
double zTargetPos;
|
577 |
#ifndef IS_MPI
|
578 |
Molecule* mol = info_->getMoleculeByGlobalIndex(index);
|
579 |
assert(mol);
|
580 |
Vector3d com = mol->getCom();
|
581 |
zTargetPos = com[whichDirection];
|
582 |
#else
|
583 |
int whicProc = info_->getMolToProc(index);
|
584 |
MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD);
|
585 |
#endif
|
586 |
return zTargetPos;
|
587 |
}
|
588 |
|
589 |
}
|