--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/14 16:16:39 696 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2004/05/27 18:59:17 1203 @@ -1,17 +1,25 @@ #include "Integrator.hpp" #include "simError.h" -#include -template ZConstraint::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) - : T(theInfo, the_ff), fz(NULL), - indexOfZConsMols(NULL) -{ +#include +const double INFINITE_TIME = 10e30; +template ZConstraint::ZConstraint(SimInfo* theInfo, + ForceFields* the_ff): T(theInfo, the_ff), + fzOut(NULL), + curZconsTime(0), + forcePolicy(NULL), + usingSMD(false), + hasZConsGap(false){ //get properties from SimInfo GenericData* data; ZConsParaData* zConsParaData; DoubleData* sampleTime; DoubleData* tolerance; + DoubleData* gap; + DoubleData* fixtime; + StringData* policy; StringData* filename; + IntData* smdFlag; double COM[3]; //by default, the direction of constraint is z @@ -22,122 +30,201 @@ template ZConstraint::ZConstraint(SimIn //estimate the force constant of harmonical potential double Kb = 1.986E-3 ; //in kcal/K - - double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2; - zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox); + double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) / + 2; + zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox); + + //creat force Subtraction policy + data = info->getProperty(ZCONSFORCEPOLICY_ID); + if (!data){ + sprintf(painCave.errMsg, + "ZConstraint Warning: User does not set force Subtraction policy, " + "PolicyByMass is used\n"); + painCave.isFatal = 0; + simError(); + + forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); + } + else{ + policy = dynamic_cast(data); + + if (!policy){ + sprintf(painCave.errMsg, + "ZConstraint Error: Convertion from GenericData to StringData failure, " + "PolicyByMass is used\n"); + painCave.isFatal = 0; + simError(); + + forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); + } + else{ + if (policy->getData() == "BYNUMBER") + forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this); + else if (policy->getData() == "BYMASS") + forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); + else{ + sprintf(painCave.errMsg, + "ZConstraint Warning: unknown force Subtraction policy, " + "PolicyByMass is used\n"); + painCave.isFatal = 0; + simError(); + forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this); + } + } + } + + //retrieve sample time of z-contraint data = info->getProperty(ZCONSTIME_ID); - - if(!data) { - - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set sample time.\n"); + + if (!data){ + sprintf(painCave.errMsg, + "ZConstraint error: If you use an ZConstraint\n" + " , you must set sample time.\n"); painCave.isFatal = 1; - simError(); + simError(); } else{ - sampleTime = dynamic_cast(data); - - if(!sampleTime){ - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); + if (!sampleTime){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); painCave.isFatal = 1; - simError(); - + simError(); } else{ this->zconsTime = sampleTime->getData(); } - } - + //retrieve output filename of z force data = info->getProperty(ZCONSFILENAME_ID); - if(!data) { - - - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set output filename of z-force.\n"); + if (!data){ + sprintf(painCave.errMsg, + "ZConstraint error: If you use an ZConstraint\n" + " , you must set output filename of z-force.\n"); painCave.isFatal = 1; - simError(); - + simError(); } else{ + filename = dynamic_cast(data); - filename = dynamic_cast(data); - - if(!filename){ - - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); + if (!filename){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); painCave.isFatal = 1; - simError(); - + simError(); } else{ this->zconsOutput = filename->getData(); } - - } //retrieve tolerance for z-constraint molecuels data = info->getProperty(ZCONSTOL_ID); - - if(!data) { - - sprintf( painCave.errMsg, - "ZConstraint error: can not get tolerance \n"); + + if (!data){ + sprintf(painCave.errMsg, "ZConstraint error: can not get tolerance \n"); painCave.isFatal = 1; - simError(); + simError(); } else{ - tolerance = dynamic_cast(data); - - if(!tolerance){ - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); + if (!tolerance){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); painCave.isFatal = 1; - simError(); - + simError(); } else{ this->zconsTol = tolerance->getData(); } + } + //quick hack here + data = info->getProperty(ZCONSGAP_ID); + + if (data){ + gap = dynamic_cast(data); + + if (!gap){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); + painCave.isFatal = 1; + simError(); + } + else{ + this->hasZConsGap = true; + this->zconsGap = gap->getData(); + } } - - //retrieve index of z-constraint molecules - data = info->getProperty(ZCONSPARADATA_ID); - if(!data) { - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set index of z-constraint molecules.\n"); - painCave.isFatal = 1; - simError(); + + + data = info->getProperty(ZCONSFIXTIME_ID); + + if (data){ + fixtime = dynamic_cast(data); + if (!fixtime){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get zconsFixTime from SimInfo\n"); + painCave.isFatal = 1; + simError(); + } + else{ + this->zconsFixTime = fixtime->getData(); + } } - else{ - - zConsParaData = dynamic_cast(data); - - if(!zConsParaData){ + else if(hasZConsGap){ + sprintf(painCave.errMsg, + "ZConstraint error: must set fixtime if already set zconsGap\n"); + painCave.isFatal = 1; + simError(); + } - sprintf( painCave.errMsg, - "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); + + + data = info->getProperty(ZCONSUSINGSMD_ID); + + if (data){ + smdFlag = dynamic_cast(data); + + if (!smdFlag){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); painCave.isFatal = 1; - simError(); - + simError(); + } + else{ + this->usingSMD= smdFlag->getData() ? true : false; } + + } + + + + //retrieve index of z-constraint molecules + data = info->getProperty(ZCONSPARADATA_ID); + if (!data){ + sprintf(painCave.errMsg, + "ZConstraint error: If you use an ZConstraint\n" + " , you must set index of z-constraint molecules.\n"); + painCave.isFatal = 1; + simError(); + } + else{ + zConsParaData = dynamic_cast(data); + + if (!zConsParaData){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); + painCave.isFatal = 1; + simError(); + } else{ - parameters = zConsParaData->getData(); //check the range of zconsIndex @@ -149,308 +236,306 @@ template ZConstraint::ZConstraint(SimIn int totalNumMol; minIndex = (*parameters)[0].zconsIndex; - if(minIndex < 0){ - sprintf( painCave.errMsg, - "ZConstraint error: index is out of range\n"); + if (minIndex < 0){ + sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); painCave.isFatal = 1; - simError(); - } + simError(); + } maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; #ifndef IS_MPI totalNumMol = nMols; #else - totalNumMol = mpiSim->getTotNmol(); + totalNumMol = mpiSim->getNMolGlobal(); #endif - - if(maxIndex > totalNumMol - 1){ - sprintf( painCave.errMsg, - "ZConstraint error: index is out of range\n"); + + if (maxIndex > totalNumMol - 1){ + sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); painCave.isFatal = 1; - simError(); + simError(); } //if user does not specify the zpos for the zconstraint molecule //its initial z coordinate will be used as default - for(int i = 0; i < parameters->size(); i++){ + for (int i = 0; i < (int) (parameters->size()); i++){ + if (!(*parameters)[i].havingZPos){ +#ifndef IS_MPI + for (int j = 0; j < nMols; j++){ + if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[j].getCOM(COM); + break; + } + } +#else + //query which processor current zconstraint molecule belongs to + int* MolToProcMap; + int whichNode; - if(!(*parameters)[i].havingZPos){ + MolToProcMap = mpiSim->getMolToProcMap(); + whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; -#ifndef IS_MPI - for(int j = 0; j < nMols; j++){ + //broadcast the zpos of current z-contraint molecule + //the node which contain this + + if (worldRank == whichNode){ + for (int j = 0; j < nMols; j++) if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ - molecules[j].getCOM(COM); - break; + molecules[j].getCOM(COM); + break; } - } -#else - //query which processor current zconstraint molecule belongs to - int *MolToProcMap; - int whichNode; - double initZPos; - MolToProcMap = mpiSim->getMolToProcMap(); - whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; - - //broadcast the zpos of current z-contraint molecule - //the node which contain this - - if (worldRank == whichNode ){ - - for(int j = 0; j < nMols; j++) - if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ - molecules[j].getCOM(COM); - break; - } - - } + } - MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD); + MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE, whichNode, + MPI_COMM_WORLD); #endif - - (*parameters)[i].zPos = COM[whichDirection]; - sprintf( painCave.errMsg, - "ZConstraint warningr: Does not specify zpos for z-constraint molecule " - "initial z coornidate will be used \n"); - painCave.isFatal = 0; - simError(); - - } - } - + (*parameters)[i].zPos = COM[whichDirection]; + + sprintf(painCave.errMsg, + "ZConstraint warning: Does not specify zpos for z-constraint molecule " + "initial z coornidate will be used \n"); + painCave.isFatal = 0; + simError(); + } + } }//end if (!zConsParaData) + }//end if (!data) - -// + + // #ifdef IS_MPI update(); #else int searchResult; - - for(int i = 0; i < nMols; i++){ - + + for (int i = 0; i < nMols; i++){ searchResult = isZConstraintMol(&molecules[i]); - - if(searchResult > -1){ - + + if (searchResult > -1){ zconsMols.push_back(&molecules[i]); massOfZConsMols.push_back(molecules[i].getTotalMass()); zPos.push_back((*parameters)[searchResult].zPos); - kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); + kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); - molecules[i].getCOM(COM); + if(usingSMD) + cantVel.push_back((*parameters)[searchResult].cantVel); + } - else - { - + else{ unconsMols.push_back(&molecules[i]); massOfUnconsMols.push_back(molecules[i].getTotalMass()); - } } - fz = new double[zconsMols.size()]; - indexOfZConsMols = new int [zconsMols.size()]; + fz.resize(zconsMols.size()); + curZPos.resize(zconsMols.size()); + indexOfZConsMols.resize(zconsMols.size()); - if(!fz || !indexOfZConsMols){ - sprintf( painCave.errMsg, - "Memory allocation failure in class Zconstraint\n"); - painCave.isFatal = 1; - simError(); - } - //determine the states of z-constraint molecules - for(int i = 0; i < zconsMols.size(); i++){ + for (size_t i = 0; i < zconsMols.size(); i++){ indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); - zconsMols[i]->getCOM(COM); - if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) - states.push_back(zcsFixed); - else - states.push_back(zcsMoving); + zconsMols[i]->getCOM(COM); + + if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ + states.push_back(zcsFixed); + + if (hasZConsGap) + endFixTime.push_back(info->getTime() + zconsFixTime); + } + else{ + states.push_back(zcsMoving); + + if (hasZConsGap) + endFixTime.push_back(INFINITE_TIME); + } + + if(usingSMD) + cantPos.push_back(COM[whichDirection]); } - + + if(usingSMD) + prevCantPos = cantPos; #endif + //get total masss of unconstraint molecules double totalMassOfUncons_local; totalMassOfUncons_local = 0; - - for(int i = 0; i < unconsMols.size(); i++) + + for (size_t i = 0; i < unconsMols.size(); i++) totalMassOfUncons_local += unconsMols[i]->getTotalMass(); - + #ifndef IS_MPI totalMassOfUncons = totalMassOfUncons_local; #else - MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); #endif - //get total number of unconstrained atoms int nUnconsAtoms_local; nUnconsAtoms_local = 0; - for(int i = 0; i < unconsMols.size(); i++) + for (int i = 0; i < (int) (unconsMols.size()); i++) nUnconsAtoms_local += unconsMols[i]->getNAtoms(); - + #ifndef IS_MPI totNumOfUnconsAtoms = nUnconsAtoms_local; #else - MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); #endif - // creat zconsWriter - fzOut = new ZConsWriter(zconsOutput.c_str()); - - if(!fzOut){ - sprintf( painCave.errMsg, - "Memory allocation failure in class Zconstraint\n"); - painCave.isFatal = 1; - simError(); - } - + forcePolicy->update(); } -template ZConstraint::~ZConstraint() -{ - if(fz) - delete[] fz; - - if(indexOfZConsMols) - delete[] indexOfZConsMols; - - if(fzOut) +template ZConstraint::~ZConstraint(){ + + if (fzOut){ delete fzOut; + } + + if (forcePolicy){ + delete forcePolicy; + } } + +/** + * + */ + #ifdef IS_MPI -template void ZConstraint::update() -{ +template void ZConstraint::update(){ double COM[3]; int index; - + zconsMols.clear(); massOfZConsMols.clear(); zPos.clear(); kz.clear(); - + cantPos.clear(); + cantVel.clear(); + unconsMols.clear(); massOfUnconsMols.clear(); - + //creat zconsMol and unconsMol lists - for(int i = 0; i < nMols; i++){ - + for (int i = 0; i < nMols; i++){ index = isZConstraintMol(&molecules[i]); - - if(index > -1){ - + + if (index > -1){ zconsMols.push_back(&molecules[i]); zPos.push_back((*parameters)[index].zPos); - kz.push_back((*parameters)[index].kRatio * zForceConst); - - massOfZConsMols.push_back(molecules[i].getTotalMass()); - - molecules[i].getCOM(COM); + kz.push_back((*parameters)[index].kRatio * zForceConst); + massOfZConsMols.push_back(molecules[i].getTotalMass()); + + if(usingSMD) + cantVel.push_back((*parameters)[index].cantVel); + } - else - { - + else{ unconsMols.push_back(&molecules[i]); massOfUnconsMols.push_back(molecules[i].getTotalMass()); - } } - //determine the states of z-constraint molecules - for(int i = 0; i < zconsMols.size(); i++){ - zconsMols[i]->getCOM(COM); - if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) - states.push_back(zcsFixed); - else - states.push_back(zcsMoving); + fz.resize(zconsMols.size()); + curZPos.resize(zconsMols.size()); + indexOfZConsMols.resize(zconsMols.size()); + + for (size_t i = 0; i < zconsMols.size(); i++){ + indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); } + + //determine the states of z-constraint molecules + for (int i = 0; i < (int) (zconsMols.size()); i++){ + zconsMols[i]->getCOM(COM); - //The reason to declare fz and indexOfZconsMols as pointer to array is - // that we want to make the MPI communication simple - if(fz) - delete[] fz; - - if(indexOfZConsMols) - delete[] indexOfZConsMols; - - if (zconsMols.size() > 0){ - fz = new double[zconsMols.size()]; - indexOfZConsMols = new int[zconsMols.size()]; - - if(!fz || !indexOfZConsMols){ - sprintf( painCave.errMsg, - "Memory allocation failure in class Zconstraint\n"); - painCave.isFatal = 1; - simError(); + if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){ + states.push_back(zcsFixed); + + if (hasZConsGap) + endFixTime.push_back(info->getTime() + zconsFixTime); } - - for(int i = 0; i < zconsMols.size(); i++){ - indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); + else{ + states.push_back(zcsMoving); + + if (hasZConsGap) + endFixTime.push_back(INFINITE_TIME); } - } - else{ - fz = NULL; - indexOfZConsMols = NULL; - } - + if(usingSMD) + cantPos.push_back(COM[whichDirection]); + } + + if(usingSMD) + prevCantPos = cantPos; + + // + forcePolicy->update(); } #endif -/** Function Name: isZConstraintMol - ** Parameter - ** Molecule* mol - ** Return value: - ** -1, if the molecule is not z-constraint molecule, - ** other non-negative values, its index in indexOfAllZConsMols vector +/** + * Function Name: isZConstraintMol + * Parameter + * Molecule* mol + * Return value: + * -1, if the molecule is not z-constraint molecule, + * other non-negative values, its index in indexOfAllZConsMols vector */ -template int ZConstraint::isZConstraintMol(Molecule* mol) -{ +template int ZConstraint::isZConstraintMol(Molecule* mol){ int index; int low; int high; int mid; index = mol->getGlobalIndex(); - + low = 0; high = parameters->size() - 1; - + //Binary Search (we have sorted the array) - while(low <= high){ - mid = (low + high) /2; + while (low <= high){ + mid = (low + high) / 2; if ((*parameters)[mid].zconsIndex == index) return mid; - else if ((*parameters)[mid].zconsIndex > index ) - high = mid -1; - else - low = mid + 1; + else if ((*parameters)[mid].zconsIndex > index) + high = mid - 1; + else + low = mid + 1; } - + return -1; } template void ZConstraint::integrate(){ - + // creat zconsWriter + fzOut = new ZConsWriter(zconsOutput.c_str(), parameters); + + if (!fzOut){ + sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); + painCave.isFatal = 1; + simError(); + } + //zero out the velocities of center of mass of unconstrained molecules //and the velocities of center of mass of every single z-constrained molecueles zeroOutVel(); - - T::integrate(); + curZconsTime = zconsTime + info->getTime(); + + T::integrate(); } - + /** * * @@ -459,38 +544,87 @@ template void ZConstraint::calcForce(in */ template void ZConstraint::calcForce(int calcPot, int calcStress){ double zsys; + double COM[3]; + double force[3]; + double zSysCOMVel; T::calcForce(calcPot, calcStress); - if (checkZConsState()) - zeroOutVel(); - + + if (hasZConsGap){ + updateZPos(); + } + + if (checkZConsState()){ + zeroOutVel(); + forcePolicy->update(); + } + zsys = calcZSys(); - cout << "---------------------------------------------------------------------" <getTime() << endl; - cout << "center of mass at z: " << zsys << endl; - //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <getTime() << endl; + //cout << "center of mass at z: " << zsys << endl; + //cout << "before calcForce, the COMVel of system is " << zSysCOMVel <doZconstraintForce(); - - //use harmonical poteintial to move the molecules to the specified positions - if (haveMovingZMols()) - this->doHarmonic(); + } - //cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; + //use external force to move the molecules to the specified positions + if (haveMovingZMols()){ + if (usingSMD) + this->doHarmonic(cantPos); + else + this->doHarmonic(zPos); + } - fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); - //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <getTime() >= curZconsTime){ + for (int i = 0; i < (int) (zconsMols.size()); i++){ + zconsMols[i]->getCOM(COM); + curZPos[i] = COM[whichDirection]; + + //if the z-constraint molecule is still moving, just record its force + if (states[i] == zcsMoving){ + fz[i] = 0; + Atom** movingZAtoms; + movingZAtoms = zconsMols[i]->getMyAtoms(); + for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ + movingZAtoms[j]->getFrc(force); + fz[i] += force[whichDirection]; + } + } + } + fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0], + &curZPos[0], &zPos[0]); + curZconsTime += zconsTime; + } + + zSysCOMVel = calcSysCOMVel(); +#ifdef IS_MPI + if (worldRank == 0){ +#endif + //cout << "after calcForce, the COMVel of system is " << zSysCOMVel < double ZConstraint::calcZSys() -{ + + +/** + * + */ + +template double ZConstraint::calcZSys(){ //calculate reference z coordinate for z-constraint molecules double totalMass_local; double totalMass; @@ -498,27 +632,27 @@ template double ZConstraint::calcZSys() double totalMZ; double massOfCurMol; double COM[3]; - + totalMass_local = 0; totalMZ_local = 0; - - for(int i = 0; i < nMols; i++){ + + for (int i = 0; i < nMols; i++){ massOfCurMol = molecules[i].getTotalMass(); molecules[i].getCOM(COM); - + totalMass_local += massOfCurMol; totalMZ_local += massOfCurMol * COM[whichDirection]; - } - - + + #ifdef IS_MPI - MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - #else + MPI_Allreduce(&totalMass_local, &totalMass, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(&totalMZ_local, &totalMZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#else totalMass = totalMass_local; totalMZ = totalMZ_local; - #endif +#endif double zsys; zsys = totalMZ / totalMass; @@ -529,415 +663,491 @@ template double ZConstraint::calcZSys() /** * */ -template void ZConstraint::thermalize( void ){ - +template void ZConstraint::thermalize(void){ T::thermalize(); zeroOutVel(); } /** * - * - * */ template void ZConstraint::zeroOutVel(){ - Atom** fixedZAtoms; double COMvel[3]; double vel[3]; + double zSysCOMVel; //zero out the velocities of center of mass of fixed z-constrained molecules - - for(int i = 0; i < zconsMols.size(); i++){ - if (states[i] == zcsFixed){ + for (int i = 0; i < (int) (zconsMols.size()); i++){ + if (states[i] == zcsFixed){ + zconsMols[i]->getCOMvel(COMvel); + //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; - zconsMols[i]->getCOMvel(COMvel); - //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; - fixedZAtoms = zconsMols[i]->getMyAtoms(); - - for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ + + for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ fixedZAtoms[j]->getVel(vel); - vel[whichDirection] -= COMvel[whichDirection]; - fixedZAtoms[j]->setVel(vel); + vel[whichDirection] -= COMvel[whichDirection]; + fixedZAtoms[j]->setVel(vel); } - zconsMols[i]->getCOMvel(COMvel); - //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; + zconsMols[i]->getCOMvel(COMvel); + //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; } - } - //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <getCOMvel(COMvel); - MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; + MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; } - for(int i = 0; i < zconsMols.size(); i++){ + for (int i = 0; i < (int) (zconsMols.size()); i++){ if (states[i] == zcsMoving){ zconsMols[i]->getCOMvel(COMvel); MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; - totalMassOfMovingZMols_local += massOfZConsMols[i]; + totalMassOfMovingZMols_local += massOfZConsMols[i]; } - } #ifndef IS_MPI MVzOfMovingMols = MVzOfMovingMols_local; totalMassOfMovingZMols = totalMassOfMovingZMols_local; #else - MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, + MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif double vzOfMovingMols; - vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); + vzOfMovingMols = MVzOfMovingMols / + (totalMassOfUncons + totalMassOfMovingZMols); //modify the velocites of unconstrained molecules Atom** unconsAtoms; - for(int i = 0; i < unconsMols.size(); i++){ - + for (int i = 0; i < (int) (unconsMols.size()); i++){ unconsAtoms = unconsMols[i]->getMyAtoms(); - for(int j = 0; j < unconsMols[i]->getNAtoms();j++){ + for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ unconsAtoms[j]->getVel(vel); vel[whichDirection] -= vzOfMovingMols; unconsAtoms[j]->setVel(vel); } - } //modify the velocities of moving z-constrained molecuels Atom** movingZAtoms; - for(int i = 0; i < zconsMols.size(); i++){ - - if (states[i] ==zcsMoving){ - + for (int i = 0; i < (int) (zconsMols.size()); i++){ + if (states[i] == zcsMoving){ movingZAtoms = zconsMols[i]->getMyAtoms(); - for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ + for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ movingZAtoms[j]->getVel(vel); vel[whichDirection] -= vzOfMovingMols; - movingZAtoms[j]->setVel(vel); - } - - } + movingZAtoms[j]->setVel(vel); + } + } + } - } - cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() < void ZConstraint::doZconstraintForce(){ +/** + * + */ +template void ZConstraint::doZconstraintForce(){ Atom** zconsAtoms; double totalFZ; double totalFZ_local; - double COMvel[3]; double COM[3]; double force[3]; + //constrain the molecules which do not reach the specified positions - - //constrain the molecules which do not reach the specified positions - //Zero Out the force of z-contrained molecules totalFZ_local = 0; //calculate the total z-contrained force of fixed z-contrained molecules - cout << "Fixed Molecules" << endl; - for(int i = 0; i < zconsMols.size(); i++){ - + + //cout << "before zero out z-constraint force on fixed z-constraint molecuels " + // << "total force is " << calcTotalForce() << endl; + + for (int i = 0; i < (int) (zconsMols.size()); i++){ if (states[i] == zcsFixed){ - zconsMols[i]->getCOM(COM); zconsAtoms = zconsMols[i]->getMyAtoms(); fz[i] = 0; - for(int j =0; j < zconsMols[i]->getNAtoms(); j++) { + for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ zconsAtoms[j]->getFrc(force); - fz[i] += force[whichDirection]; - } + fz[i] += force[whichDirection]; + } totalFZ_local += fz[i]; - cout << "index: " << indexOfZConsMols[i] - <<"\tcurrent zpos: " << COM[whichDirection] - << "\tcurrent fz: " <getNAtoms(); zconsAtoms = zconsMols[i]->getMyAtoms(); - - for(int j =0; j < nAtomOfCurZConsMol; j++) { - force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + + for (int j = 0; j < nAtomOfCurZConsMol; j++){ + //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + force[whichDirection] = -forcePolicy->getZFOfFixedZMols(zconsMols[i], + zconsAtoms[j], + fz[i]); zconsAtoms[j]->addFrc(force); } - } - } - cout << "after zero out z-constraint force on fixed z-constraint molecuels " - << "total force is " << calcTotalForce() << endl; - //calculate the number of atoms of moving z-constrained molecules - int nMovingZAtoms_local; - int nMovingZAtoms; - - nMovingZAtoms_local = 0; - for(int i = 0; i < zconsMols.size(); i++) - if(states[i] == zcsMoving) - nMovingZAtoms_local += zconsMols[i]->getNAtoms(); - -#ifdef IS_MPI - MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else - totalFZ = totalFZ_local; - nMovingZAtoms = nMovingZAtoms_local; -#endif + //cout << "after zero out z-constraint force on fixed z-constraint molecuels " + // << "total force is " << calcTotalForce() << endl; - force[0]= 0; - force[1]= 0; - force[2]= 0; - force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); + force[0] = 0; + force[1] = 0; + force[2] = 0; + //modify the forces of unconstrained molecules - int accessCount = 0; - for(int i = 0; i < unconsMols.size(); i++){ - - Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); - - for(int j = 0; j < unconsMols[i]->getNAtoms(); j++) - unconsAtoms[j]->addFrc(force); - + for (int i = 0; i < (int) (unconsMols.size()); i++){ + Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); + + for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ + //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); + force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j], + totalFZ); + unconsAtoms[j]->addFrc(force); + } } - //modify the forces of moving z-constrained molecules - for(int i = 0; i < zconsMols.size(); i++) { + //modify the forces of moving z-constrained molecules + for (int i = 0; i < (int) (zconsMols.size()); i++){ if (states[i] == zcsMoving){ - Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); - for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) + for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ + //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); + force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j], + totalFZ); movingZAtoms[j]->addFrc(force); + } } } - - cout << "after substracting z-constraint force from moving molecuels " - << "total force is " << calcTotalForce() << endl; - + // cout << "after substracting z-constraint force from moving molecuels " + // << "total force is " << calcTotalForce() << endl; } -template bool ZConstraint::checkZConsState(){ - double COM[3]; - double diff; - - bool changed; - - changed = false; - - for(int i =0; i < zconsMols.size(); i++){ - - zconsMols[i]->getCOM(COM); - diff = fabs(COM[whichDirection] - zPos[i]); - if ( diff <= zconsTol && states[i] == zcsMoving){ - states[i] = zcsFixed; - changed = true; - } - else if ( diff > zconsTol && states[i] == zcsFixed){ - states[i] = zcsMoving; - changed = true; - } - - } - - return changed; -} - -template bool ZConstraint::haveFixedZMols(){ - for(int i = 0; i < zconsMols.size(); i++) - if (states[i] == zcsFixed) - return true; - - return false; -} - - /** - * - */ -template bool ZConstraint::haveMovingZMols(){ - for(int i = 0; i < zconsMols.size(); i++) - if (states[i] == zcsMoving) - return true; - - return false; - -} - -/** * * */ -template void ZConstraint::doHarmonic(){ +template void ZConstraint::doHarmonic(vector& resPos){ double force[3]; double harmonicU; double harmonicF; double COM[3]; double diff; + double totalFZ_local; double totalFZ; - + force[0] = 0; force[1] = 0; force[2] = 0; - totalFZ = 0; + totalFZ_local = 0; - cout << "Moving Molecules" << endl; - for(int i = 0; i < zconsMols.size(); i++) { - + for (int i = 0; i < (int) (zconsMols.size()); i++){ if (states[i] == zcsMoving){ zconsMols[i]->getCOM(COM); - cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; - - diff = COM[whichDirection] -zPos[i]; - + // cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] + // << "\tcurrent zpos: " << COM[whichDirection] << endl; + + diff = COM[whichDirection] - resPos[i]; + harmonicU = 0.5 * kz[i] * diff * diff; - info->lrPot += harmonicU; + info->lrPot += harmonicU; - harmonicF = - kz[i] * diff; - totalFZ += harmonicF; + harmonicF = -kz[i] * diff; + totalFZ_local += harmonicF; - // - force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); - + //adjust force + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); - for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) - movingZAtoms[j]->addFrc(force); + for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){ + //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); + force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], + movingZAtoms[j], + harmonicF); + movingZAtoms[j]->addFrc(force); + } } - } - - force[0]= 0; - force[1]= 0; - force[2]= 0; - force[whichDirection] = -totalFZ /totNumOfUnconsAtoms; +#ifndef IS_MPI + totalFZ = totalFZ_local; +#else + MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + + //cout << "before substracting harmonic force from moving molecuels " + // << "total force is " << calcTotalForce() << endl; + + force[0] = 0; + force[1] = 0; + force[2] = 0; + //modify the forces of unconstrained molecules - for(int i = 0; i < unconsMols.size(); i++){ - - Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); - - for(int j = 0; j < unconsMols[i]->getNAtoms(); j++) - unconsAtoms[j]->addFrc(force); + for (int i = 0; i < (int) (unconsMols.size()); i++){ + Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); + + for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){ + //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms; + force[whichDirection] = -forcePolicy->getHFOfUnconsMols(unconsAtoms[j], + totalFZ); + unconsAtoms[j]->addFrc(force); + } } + //cout << "after substracting harmonic force from moving molecuels " + // << "total force is " << calcTotalForce() << endl; } -template double ZConstraint::calcMovingMolsCOMVel() -{ +/** + * + */ + +template bool ZConstraint::checkZConsState(){ + double COM[3]; + double diff; + + int changed_local; + int changed; + + changed_local = 0; + + for (int i = 0; i < (int) (zconsMols.size()); i++){ + zconsMols[i]->getCOM(COM); + diff = fabs(COM[whichDirection] - zPos[i]); + if (diff <= zconsTol && states[i] == zcsMoving){ + states[i] = zcsFixed; + changed_local = 1; + + if(usingSMD) + prevCantPos = cantPos; + + if (hasZConsGap) + endFixTime[i] = info->getTime() + zconsFixTime; + } + else if (diff > zconsTol && states[i] == zcsFixed){ + states[i] = zcsMoving; + changed_local = 1; + + if(usingSMD) + cantPos = prevCantPos; + + if (hasZConsGap) + endFixTime[i] = INFINITE_TIME; + } + } + +#ifndef IS_MPI + changed = changed_local; +#else + MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); +#endif + + return (changed > 0); +} + +template bool ZConstraint::haveFixedZMols(){ + int havingFixed_local; + int havingFixed; + + havingFixed_local = 0; + + for (int i = 0; i < (int) (zconsMols.size()); i++) + if (states[i] == zcsFixed){ + havingFixed_local = 1; + break; + } + +#ifndef IS_MPI + havingFixed = havingFixed_local; +#else + MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#endif + + return (havingFixed > 0); +} + + +/** + * + */ +template bool ZConstraint::haveMovingZMols(){ + int havingMoving_local; + int havingMoving; + + havingMoving_local = 0; + + for (int i = 0; i < (int) (zconsMols.size()); i++) + if (states[i] == zcsMoving){ + havingMoving_local = 1; + break; + } + +#ifndef IS_MPI + havingMoving = havingMoving_local; +#else + MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#endif + + return (havingMoving > 0); +} + +/** + * + */ + +template double ZConstraint::calcMovingMolsCOMVel(){ double MVzOfMovingMols_local; double MVzOfMovingMols; double totalMassOfMovingZMols_local; double totalMassOfMovingZMols; double COMvel[3]; - + MVzOfMovingMols_local = 0; totalMassOfMovingZMols_local = 0; - for(int i =0; i < unconsMols.size(); i++){ + for (int i = 0; i < unconsMols.size(); i++){ unconsMols[i]->getCOMvel(COMvel); - MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; + MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; } - for(int i = 0; i < zconsMols.size(); i++){ - + for (int i = 0; i < zconsMols.size(); i++){ if (states[i] == zcsMoving){ zconsMols[i]->getCOMvel(COMvel); MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; - totalMassOfMovingZMols_local += massOfZConsMols[i]; + totalMassOfMovingZMols_local += massOfZConsMols[i]; } - } #ifndef IS_MPI MVzOfMovingMols = MVzOfMovingMols_local; totalMassOfMovingZMols = totalMassOfMovingZMols_local; #else - MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, + MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif double vzOfMovingMols; - vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); + vzOfMovingMols = MVzOfMovingMols / + (totalMassOfUncons + totalMassOfMovingZMols); return vzOfMovingMols; } +/** + * + */ -template double ZConstraint::calcSysCOMVel() -{ +template double ZConstraint::calcSysCOMVel(){ double COMvel[3]; - double tempMVz = 0; - - for(int i =0 ; i < nMols; i++){ - molecules[i].getCOMvel(COMvel); - tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection]; - } - + double tempMVz_local; + double tempMVz; double massOfZCons_local; double massOfZCons; + + tempMVz_local = 0; + + for (int i = 0 ; i < nMols; i++){ + molecules[i].getCOMvel(COMvel); + tempMVz_local += molecules[i].getTotalMass() * COMvel[whichDirection]; + } + massOfZCons_local = 0; - - for(int i = 0; i < massOfZConsMols.size(); i++){ + + for (int i = 0; i < (int) (massOfZConsMols.size()); i++){ massOfZCons_local += massOfZConsMols[i]; } #ifndef IS_MPI massOfZCons = massOfZCons_local; + tempMVz = tempMVz_local; #else - MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); + MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif - return tempMVz /(totalMassOfUncons + massOfZCons); + return tempMVz / (totalMassOfUncons + massOfZCons); } -template double ZConstraint::calcTotalForce(){ +/** + * + */ +template double ZConstraint::calcTotalForce(){ double force[3]; double totalForce_local; double totalForce; totalForce_local = 0; - for(int i = 0; i < nAtoms; i++){ + for (int i = 0; i < nAtoms; i++){ atoms[i]->getFrc(force); totalForce_local += force[whichDirection]; } @@ -945,9 +1155,139 @@ template double ZConstraint::calcTotalF #ifndef IS_MPI totalForce = totalForce_local; #else - MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE, MPI_SUM, + MPI_COMM_WORLD); #endif return totalForce; +} +/** + * + */ + +template void ZConstraint::PolicyByNumber::update(){ + //calculate the number of atoms of moving z-constrained molecules + int nMovingZAtoms_local; + int nMovingZAtoms; + + nMovingZAtoms_local = 0; + for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) + if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ + nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms(); + } + +#ifdef IS_MPI + MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#else + nMovingZAtoms = nMovingZAtoms_local; +#endif + totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms; } + +template double ZConstraint::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, + Atom* atom, + double totalForce){ + return totalForce / mol->getNAtoms(); +} + +template double ZConstraint::PolicyByNumber::getZFOfMovingMols(Atom* atom, + double totalForce){ + return totalForce / totNumOfMovingAtoms; +} + +template double ZConstraint::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, + Atom* atom, + double totalForce){ + return totalForce / mol->getNAtoms(); +} + +template double ZConstraint::PolicyByNumber::getHFOfUnconsMols(Atom* atom, + double totalForce){ + return totalForce / zconsIntegrator->totNumOfUnconsAtoms; +} + +/** + * + */ + +template void ZConstraint::PolicyByMass::update(){ + //calculate the number of atoms of moving z-constrained molecules + double massOfMovingZAtoms_local; + double massOfMovingZAtoms; + + massOfMovingZAtoms_local = 0; + for (int i = 0; i < (int) ((zconsIntegrator->zconsMols).size()); i++) + if ((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)){ + massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass(); + } + +#ifdef IS_MPI + MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE, + MPI_SUM, MPI_COMM_WORLD); +#else + massOfMovingZAtoms = massOfMovingZAtoms_local; +#endif + totMassOfMovingAtoms = massOfMovingZAtoms + + zconsIntegrator->totalMassOfUncons; +} + +template double ZConstraint::PolicyByMass::getZFOfFixedZMols(Molecule* mol, + Atom* atom, + double totalForce){ + return totalForce * atom->getMass() / mol->getTotalMass(); +} + +template double ZConstraint::PolicyByMass::getZFOfMovingMols(Atom* atom, + double totalForce){ + return totalForce * atom->getMass() / totMassOfMovingAtoms; +} + +template double ZConstraint::PolicyByMass::getHFOfFixedZMols(Molecule* mol, + Atom* atom, + double totalForce){ + return totalForce * atom->getMass() / mol->getTotalMass(); +} + +template double ZConstraint::PolicyByMass::getHFOfUnconsMols(Atom* atom, + double totalForce){ + return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons; +} + +template void ZConstraint::updateZPos(){ + double curTime; + double COM[3]; + + curTime = info->getTime(); + + for (size_t i = 0; i < zconsMols.size(); i++){ + + if (states[i] == zcsFixed && curTime >= endFixTime[i]){ + zPos[i] += zconsGap; + + if (usingSMD){ + zconsMols[i]->getCOM(COM); + cantPos[i] = COM[whichDirection]; + } + + } + + } + +} + +template void ZConstraint::updateCantPos(){ + double curTime; + double dt; + + curTime = info->getTime(); + dt = info->dt; + + for (size_t i = 0; i < zconsMols.size(); i++){ + if (states[i] == zcsMoving){ + cantPos[i] += cantVel[i] * dt; + } + } + +}