--- trunk/src/constraints/ZconstraintForceManager.cpp 2009/11/25 20:02:06 1390 +++ trunk/src/constraints/ZconstraintForceManager.cpp 2012/09/10 18:38:44 1796 @@ -36,7 +36,8 @@ * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -45,6 +46,10 @@ #include "utils/simError.h" #include "utils/PhysicalConstants.hpp" #include "utils/StringUtils.hpp" +#ifdef IS_MPI +#include +#endif + namespace OpenMD { ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) { currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); @@ -156,8 +161,8 @@ namespace OpenMD { #ifndef IS_MPI totMassUnconsMols_ = totMassUnconsMols_local; #else - MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_REALTYPE, - MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, + MPI::REALTYPE, MPI::SUM); #endif // creat zconsWriter @@ -240,8 +245,8 @@ namespace OpenMD { currZconsTime_ = currSnapshot_->getTime(); } - void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){ - ForceManager::calcForces(needPotential, needStress); + void ZconstraintForceManager::calcForces(){ + ForceManager::calcForces(); if (usingZconsGap_){ updateZPos(); @@ -282,18 +287,21 @@ namespace OpenMD { Vector3d vel; std::list::iterator i; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; Molecule::IntegrableObjectIterator ii; //zero out the velocities of center of mass of fixed z-constrained molecules for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) { + mol = i->mol; comVel = mol->getComVel(); - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - vel = integrableObject->getVel(); + + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + vel = sd->getVel(); vel[whichDirection] -= comVel[whichDirection]; - integrableObject->setVel(vel); + sd->setVel(vel); } } @@ -318,33 +326,37 @@ namespace OpenMD { #ifndef IS_MPI pzMovingMols = pzMovingMols_local; #else - MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_REALTYPE, - MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&pzMovingMols_local, &pzMovingMols, 1, + MPI::REALTYPE, MPI::SUM); #endif RealType vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_); //modify the velocities of moving z-constrained molecuels for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) { + mol = i->mol; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - vel = integrableObject->getVel(); + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + vel = sd->getVel(); vel[whichDirection] -= vzMovingMols; - integrableObject->setVel(vel); + sd->setVel(vel); } } //modify the velocites of unconstrained molecules for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) { + mol =*j; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - vel = integrableObject->getVel(); + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + vel = sd->getVel(); vel[whichDirection] -= vzMovingMols; - integrableObject->setVel(vel); + sd->setVel(vel); } } @@ -366,16 +378,18 @@ namespace OpenMD { //calculate the total z-contrained force of fixed z-contrained molecules std::list::iterator i; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; Molecule::IntegrableObjectIterator ii; for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) { + mol = i->mol; i->fz = 0.0; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - force = integrableObject->getFrc(); + for( sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + force = sd->getFrc(); i->fz += force[whichDirection]; } totalFZ_local += i->fz; @@ -383,7 +397,8 @@ namespace OpenMD { //calculate total z-constraint force #ifdef IS_MPI - MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&totalFZ_local, &totalFZ, 1, + MPI::REALTYPE, MPI::SUM); #else totalFZ = totalFZ_local; #endif @@ -391,35 +406,41 @@ namespace OpenMD { // apply negative to fixed z-constrained molecues; for ( i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) { + mol = i->mol; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { - force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz); - integrableObject->addFrc(force); + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + + force[whichDirection] = -getZFOfFixedZMols(mol, sd, i->fz); + sd->addFrc(force); } } //modify the forces of moving z-constrained molecules for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) { + mol = i->mol; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + force[whichDirection] = -getZFOfMovingMols(mol,totalFZ); - integrableObject->addFrc(force); + sd->addFrc(force); } } //modify the forces of unconstrained molecules std::vector::iterator j; for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) { + mol =*j; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + force[whichDirection] = -getZFOfMovingMols(mol, totalFZ); - integrableObject->addFrc(force); + sd->addFrc(force); } } @@ -431,8 +452,9 @@ namespace OpenMD { Vector3d force(0.0); Vector3d com; RealType totalFZ_local = 0; + RealType lrPot; std::list::iterator i; - StuntDouble* integrableObject; + StuntDouble* sd; Molecule::IntegrableObjectIterator ii; Molecule* mol; for ( i = movingZMols_.begin(); i != movingZMols_.end(); ++i) { @@ -441,34 +463,39 @@ namespace OpenMD { RealType resPos = usingSMD_? i->cantPos : i->param.zTargetPos; RealType diff = com[whichDirection] - resPos; RealType harmonicU = 0.5 * i->param.kz * diff * diff; - currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU; + lrPot = currSnapshot_->getLongRangePotential(); + lrPot += harmonicU; + currSnapshot_->setLongRangePotential(lrPot); RealType harmonicF = -i->param.kz * diff; totalFZ_local += harmonicF; //adjust force - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { - force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF); - integrableObject->addFrc(force); + force[whichDirection] = getHFOfFixedZMols(mol, sd, harmonicF); + sd->addFrc(force); } } #ifndef IS_MPI totalFZ = totalFZ_local; #else - MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&totalFZ_local, &totalFZ, 1, MPI::REALTYPE, + MPI::SUM); #endif //modify the forces of unconstrained molecules std::vector::iterator j; for ( j = unzconsMols_.begin(); j != unzconsMols_.end(); ++j) { + mol = *j; - for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + for(sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { + force[whichDirection] = getHFOfUnconsMols(mol, totalFZ); - integrableObject->addFrc(force); + sd->addFrc(force); } } @@ -527,7 +554,7 @@ namespace OpenMD { #ifndef IS_MPI changed = changed_local; #else - MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&changed_local, &changed, 1, MPI::INT, MPI::SUM); #endif return (changed > 0); @@ -540,8 +567,8 @@ namespace OpenMD { #ifndef IS_MPI havingFixed = havingFixed_local; #else - MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM, - MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&havingFixed_local, &havingFixed, 1, + MPI::INT, MPI::SUM); #endif return havingFixed > 0; @@ -557,8 +584,8 @@ namespace OpenMD { #ifndef IS_MPI havingMoving = havingMoving_local; #else - MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM, - MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&havingMoving_local, &havingMoving, 1, + MPI::INT, MPI::SUM); #endif return havingMoving > 0; @@ -573,8 +600,8 @@ namespace OpenMD { } #ifdef IS_MPI - MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_REALTYPE, - MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, + 1, MPI::REALTYPE, MPI::SUM); #else totMassMovingZMols_ = totMassMovingZMols_local; #endif @@ -621,7 +648,7 @@ namespace OpenMD { zTargetPos = com[whichDirection]; #else int whicProc = info_->getMolToProc(index); - MPI_Bcast(&zTargetPos, 1, MPI_REALTYPE, whicProc, MPI_COMM_WORLD); + MPI::COMM_WORLD.Bcast(&zTargetPos, 1, MPI::REALTYPE, whicProc); #endif return zTargetPos; }