111 |
|
|
112 |
|
//estimate the force constant of harmonical potential |
113 |
|
Mat3x3d hmat = currSnapshot_->getHmat(); |
114 |
< |
double halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2; |
115 |
< |
double targetTemp; |
114 |
> |
RealType halfOfLargestBox = std::max(hmat(0, 0), std::max(hmat(1, 1), hmat(2, 2))) /2; |
115 |
> |
RealType targetTemp; |
116 |
|
if (simParam->haveTargetTemp()) { |
117 |
|
targetTemp = simParam->getTargetTemp(); |
118 |
|
} else { |
119 |
|
targetTemp = 298.0; |
120 |
|
} |
121 |
< |
double zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox); |
121 |
> |
RealType zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox); |
122 |
|
|
123 |
|
int nZconstraints = simParam->getNZconsStamps(); |
124 |
|
std::vector<ZConsStamp*> stamp = simParam->getZconsStamps(); |
148 |
|
update(); |
149 |
|
|
150 |
|
//calculate masss of unconstraint molecules in the whole system (never change during the simulation) |
151 |
< |
double totMassUnconsMols_local = 0.0; |
151 |
> |
RealType totMassUnconsMols_local = 0.0; |
152 |
|
std::vector<Molecule*>::iterator j; |
153 |
|
for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) { |
154 |
|
totMassUnconsMols_local += (*j)->getMass(); |
156 |
|
#ifndef IS_MPI |
157 |
|
totMassUnconsMols_ = totMassUnconsMols_local; |
158 |
|
#else |
159 |
< |
MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_DOUBLE, |
159 |
> |
MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_REALTYPE, |
160 |
|
MPI_SUM, MPI_COMM_WORLD); |
161 |
|
#endif |
162 |
|
|
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]); |
197 |
> |
RealType diff = fabs(zmol.param.zTargetPos - com[whichDirection]); |
198 |
|
if (diff < zconsTol_) { |
199 |
|
fixedZMols_.push_back(zmol); |
200 |
|
} else { |
299 |
|
|
300 |
|
// calculate the vz of center of mass of moving molecules(include unconstrained molecules |
301 |
|
// and moving z-constrained molecules) |
302 |
< |
double pzMovingMols_local = 0.0; |
303 |
< |
double pzMovingMols; |
302 |
> |
RealType pzMovingMols_local = 0.0; |
303 |
> |
RealType pzMovingMols; |
304 |
|
|
305 |
|
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) { |
306 |
|
mol = i->mol; |
318 |
|
#ifndef IS_MPI |
319 |
|
pzMovingMols = pzMovingMols_local; |
320 |
|
#else |
321 |
< |
MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_DOUBLE, |
321 |
> |
MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_REALTYPE, |
322 |
|
MPI_SUM, MPI_COMM_WORLD); |
323 |
|
#endif |
324 |
|
|
325 |
< |
double vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_); |
325 |
> |
RealType vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_); |
326 |
|
|
327 |
|
//modify the velocities of moving z-constrained molecuels |
328 |
|
for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) { |
352 |
|
|
353 |
|
|
354 |
|
void ZconstraintForceManager::doZconstraintForce(){ |
355 |
< |
double totalFZ; |
356 |
< |
double totalFZ_local; |
355 |
> |
RealType totalFZ; |
356 |
> |
RealType totalFZ_local; |
357 |
|
Vector3d com; |
358 |
|
Vector3d force(0.0); |
359 |
|
|
383 |
|
|
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_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
387 |
|
#else |
388 |
|
totalFZ = totalFZ_local; |
389 |
|
#endif |
427 |
|
|
428 |
|
|
429 |
|
void ZconstraintForceManager::doHarmonic(){ |
430 |
< |
double totalFZ; |
430 |
> |
RealType totalFZ; |
431 |
|
Vector3d force(0.0); |
432 |
|
Vector3d com; |
433 |
< |
double totalFZ_local = 0; |
433 |
> |
RealType totalFZ_local = 0; |
434 |
|
std::list<ZconstraintMol>::iterator i; |
435 |
|
StuntDouble* integrableObject; |
436 |
|
Molecule::IntegrableObjectIterator ii; |
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; |
441 |
> |
RealType resPos = usingSMD_? i->cantPos : i->param.zTargetPos; |
442 |
> |
RealType diff = com[whichDirection] - resPos; |
443 |
> |
RealType harmonicU = 0.5 * i->param.kz * diff * diff; |
444 |
|
currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU; |
445 |
< |
double harmonicF = -i->param.kz * diff; |
445 |
> |
RealType harmonicF = -i->param.kz * diff; |
446 |
|
totalFZ_local += harmonicF; |
447 |
|
|
448 |
|
//adjust force |
457 |
|
#ifndef IS_MPI |
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_REALTYPE, MPI_SUM, MPI_COMM_WORLD); |
461 |
|
#endif |
462 |
|
|
463 |
|
//modify the forces of unconstrained molecules |
476 |
|
|
477 |
|
bool ZconstraintForceManager::checkZConsState(){ |
478 |
|
Vector3d com; |
479 |
< |
double diff; |
479 |
> |
RealType diff; |
480 |
|
int changed_local = 0; |
481 |
|
|
482 |
|
std::list<ZconstraintMol>::iterator i; |
566 |
|
|
567 |
|
void ZconstraintForceManager::calcTotalMassMovingZMols(){ |
568 |
|
|
569 |
< |
double totMassMovingZMols_local = 0.0; |
569 |
> |
RealType 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(); |
573 |
|
} |
574 |
|
|
575 |
|
#ifdef IS_MPI |
576 |
< |
MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_DOUBLE, |
576 |
> |
MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_REALTYPE, |
577 |
|
MPI_SUM, MPI_COMM_WORLD); |
578 |
|
#else |
579 |
|
totMassMovingZMols_ = totMassMovingZMols_local; |
581 |
|
|
582 |
|
} |
583 |
|
|
584 |
< |
double ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, double totalForce){ |
584 |
> |
RealType ZconstraintForceManager::getZFOfFixedZMols(Molecule* mol, StuntDouble* sd, RealType totalForce){ |
585 |
|
return totalForce * sd->getMass() / mol->getMass(); |
586 |
|
} |
587 |
|
|
588 |
< |
double ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, double totalForce){ |
588 |
> |
RealType ZconstraintForceManager::getZFOfMovingMols(Molecule* mol, RealType totalForce){ |
589 |
|
return totalForce * mol->getMass() / (totMassUnconsMols_ + totMassMovingZMols_); |
590 |
|
} |
591 |
|
|
592 |
< |
double ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, double totalForce){ |
592 |
> |
RealType ZconstraintForceManager::getHFOfFixedZMols(Molecule* mol, StuntDouble*sd, RealType totalForce){ |
593 |
|
return totalForce * sd->getMass() / mol->getMass(); |
594 |
|
} |
595 |
|
|
596 |
< |
double ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, double totalForce){ |
596 |
> |
RealType ZconstraintForceManager::getHFOfUnconsMols(Molecule* mol, RealType totalForce){ |
597 |
|
return totalForce * mol->getMass() / totMassUnconsMols_; |
598 |
|
} |
599 |
|
|
600 |
|
void ZconstraintForceManager::updateZPos(){ |
601 |
< |
double curTime = currSnapshot_->getTime(); |
601 |
> |
RealType curTime = currSnapshot_->getTime(); |
602 |
|
std::list<ZconstraintMol>::iterator i; |
603 |
|
for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) { |
604 |
|
i->param.zTargetPos += zconsGap_; |
612 |
|
} |
613 |
|
} |
614 |
|
|
615 |
< |
double ZconstraintForceManager::getZTargetPos(int index){ |
616 |
< |
double zTargetPos; |
615 |
> |
RealType ZconstraintForceManager::getZTargetPos(int index){ |
616 |
> |
RealType zTargetPos; |
617 |
|
#ifndef IS_MPI |
618 |
|
Molecule* mol = info_->getMoleculeByGlobalIndex(index); |
619 |
|
assert(mol); |
621 |
|
zTargetPos = com[whichDirection]; |
622 |
|
#else |
623 |
|
int whicProc = info_->getMolToProc(index); |
624 |
< |
MPI_Bcast(&zTargetPos, 1, MPI_DOUBLE, whicProc, MPI_COMM_WORLD); |
624 |
> |
MPI_Bcast(&zTargetPos, 1, MPI_REALTYPE, whicProc, MPI_COMM_WORLD); |
625 |
|
#endif |
626 |
|
return zTargetPos; |
627 |
|
} |