--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/11 19:40:06 676 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2004/05/27 18:59:17 1203 @@ -1,398 +1,657 @@ #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; - IndexData* index; + ZConsParaData* zConsParaData; DoubleData* sampleTime; + DoubleData* tolerance; + DoubleData* gap; + DoubleData* fixtime; + StringData* policy; StringData* filename; - - //retrieve index of z-constraint molecules - data = info->getProperty("zconsindex"); - if(!data) { + IntData* smdFlag; + double COM[3]; - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set index of z-constraint molecules.\n"); - painCave.isFatal = 1; - simError(); + //by default, the direction of constraint is z + // 0 --> x + // 1 --> y + // 2 --> z + whichDirection = 2; + + //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); + + //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{ - index = dynamic_cast(data); - - if(!index){ + policy = dynamic_cast(data); - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); - painCave.isFatal = 1; - simError(); - + 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{ - - indexOfAllZConsMols = index->getIndexData(); - - //the maximum value of index is the last one(we sorted the index data in SimSetup.cpp) - int maxIndex; - int minIndex; - int totalNumMol; - - minIndex = indexOfAllZConsMols[0]; - if(minIndex < 0){ - sprintf( painCave.errMsg, - "ZConstraint error: index is out of range\n"); - painCave.isFatal = 1; - simError(); - } - - maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1]; - -#ifndef IS_MPI - totalNumMol = nMols; -#else - totalNumMol = mpiSim->getTotNmol(); -#endif - - if(maxIndex > totalNumMol - 1){ - sprintf( painCave.errMsg, - "ZConstraint error: index is out of range\n"); - painCave.isFatal = 1; - simError(); - + 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"); - - if(!data) { - - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set sample time.\n"); + data = info->getProperty(ZCONSTIME_ID); + + 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"); - if(!data) { - - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set output filename of z-force.\n"); + //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"); painCave.isFatal = 1; - simError(); + simError(); + } + else{ + filename = dynamic_cast(data); + if (!filename){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); + painCave.isFatal = 1; + 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"); + painCave.isFatal = 1; + simError(); + } else{ + tolerance = dynamic_cast(data); - filename = dynamic_cast(data); - - if(!filename){ + if (!tolerance){ + sprintf(painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); + painCave.isFatal = 1; + simError(); + } + else{ + this->zconsTol = tolerance->getData(); + } + } - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); + //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(); - + simError(); } else{ - this->zconsOutput = filename->getData(); + this->hasZConsGap = true; + this->zconsGap = gap->getData(); } - + } + + + 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 if(hasZConsGap){ + sprintf(painCave.errMsg, + "ZConstraint error: must set fixtime if already set zconsGap\n"); + painCave.isFatal = 1; + simError(); + } + + + + 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(); + } + 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 + //and the minimum value of index is the first one (we already sorted the data) + //the maximum value of index is the last one + + int maxIndex; + int minIndex; + int totalNumMol; + + minIndex = (*parameters)[0].zconsIndex; + if (minIndex < 0){ + sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); + painCave.isFatal = 1; + simError(); + } + + maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; + +#ifndef IS_MPI + totalNumMol = nMols; +#else + totalNumMol = mpiSim->getNMolGlobal(); +#endif + + if (maxIndex > totalNumMol - 1){ + sprintf(painCave.errMsg, "ZConstraint error: index is out of range\n"); + painCave.isFatal = 1; + 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 < (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; + + 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, whichNode, + MPI_COMM_WORLD); +#endif + + (*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; - double COM[3]; - - 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); - 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 (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); + + 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]); } - for(int i = 0; i < zconsMols.size(); i++) - indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); - + if(usingSMD) + prevCantPos = cantPos; #endif + + //get total masss of unconstraint molecules + double totalMassOfUncons_local; + totalMassOfUncons_local = 0; + + 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); +#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); -#endif + MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#endif - - - 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]); - massOfZConsMols.push_back(molecules[i].getTotalMass()); - - molecules[i].getCOM(COM); + zPos.push_back((*parameters)[index].zPos); + 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()); - } } + + 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(); + } - //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; + //determine the states of z-constraint molecules + for (int i = 0; i < (int) (zconsMols.size()); i++){ + + zconsMols[i]->getCOM(COM); - 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 = indexOfAllZConsMols.size() - 1; - + high = parameters->size() - 1; + //Binary Search (we have sorted the array) - while(low <= high){ - mid = (low + high) /2; - if (indexOfAllZConsMols[mid] == index) + while (low <= high){ + mid = (low + high) / 2; + if ((*parameters)[mid].zconsIndex == index) return mid; - else if (indexOfAllZConsMols[mid] > index ) - high = mid -1; - else - low = mid + 1; + else if ((*parameters)[mid].zconsIndex > index) + high = mid - 1; + else + low = mid + 1; } - + return -1; } -/** - * Description: - * Reset the z coordinates - */ 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(); } - + /** * * * * - */ - - + */ 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(); + zSysCOMVel = calcSysCOMVel(); +#ifdef IS_MPI + if (worldRank == 0){ +#endif + //cout << "---------------------------------------------------------------------" <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(); - - fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); - + } + + //use external force to move the molecules to the specified positions + if (haveMovingZMols()){ + if (usingSMD) + this->doHarmonic(cantPos); + else + this->doHarmonic(zPos); + } + + //write out forces and current positions of z-constraint molecules + if (info->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; double totalMZ_local; double totalMZ; - double massOfUncons_local; double massOfCurMol; double COM[3]; - + totalMass_local = 0; - totalMass = 0; totalMZ_local = 0; - totalMZ = 0; - massOfUncons_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]; - - if(isZConstraintMol(&molecules[i]) == -1){ - - massOfUncons_local += massOfCurMol; - } - } - - + + #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); - MPI_Allreduce(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + 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; - totalMassOfUncons = massOfUncons_local; #endif double zsys; @@ -404,258 +663,631 @@ 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); 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; } - } - + + //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[i]->getNAtoms(); 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[i]->getNAtoms(); 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); + } } - } + + zSysCOMVel = calcSysCOMVel(); +#ifdef IS_MPI + if (worldRank == 0){ +#endif + //cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl; +#ifdef IS_MPI + } +#endif } -template void ZConstraint::doZconstraintForce(){ +/** + * + */ +template void ZConstraint::doZconstraintForce(){ Atom** zconsAtoms; double totalFZ; double totalFZ_local; - double COMvel[3]; double COM[3]; double force[3]; - double zsys; - int nMovingZMols_local; - int nMovingZMols; + //constrain the molecules which do not reach the specified positions - //constrain the molecules which do not reach the specified positions - - zsys = calcZSys(); - cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl; - //Zero Out the force of z-contrained molecules totalFZ_local = 0; //calculate the total z-contrained force of fixed z-contrained molecules - 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] << endl; - + //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] + // <<"\tcurrent zpos: " << COM[whichDirection] + // << "\tcurrent fz: " <getMyAtoms(); - - for(int j = 0; j < unconsMols[i]->getNAtoms(); j++) - unconsAtoms[j]->addFrc(force); - - } + // apply negative to fixed z-constrained molecues; + force[0] = 0; + force[1] = 0; + force[2] = 0; - //modify the velocities of moving z-constrained molecules - for(int i = 0; i < zconsMols.size(); i++) { - if (states[i] == zcsMoving){ - - Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); - - for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) - movingZAtoms[j]->addFrc(force); - } - } - - // apply negative to fixed z-constrained molecues; - force[0]= 0; - force[1]= 0; - force[2]= 0; - - 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){ int nAtomOfCurZConsMol = zconsMols[i]->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; + + + force[0] = 0; + force[1] = 0; + force[2] = 0; + + //modify the forces of unconstrained molecules + 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 < (int) (zconsMols.size()); i++){ + if (states[i] == zcsMoving){ + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + + 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; } -template bool ZConstraint::checkZConsState(){ +/** + * + * + */ + +template void ZConstraint::doHarmonic(vector& resPos){ + double force[3]; + double harmonicU; + double harmonicF; double COM[3]; double diff; - - bool changed; - - changed = false; - - for(int i =0; i < zconsMols.size(); i++){ + double totalFZ_local; + double totalFZ; - zconsMols[i]->getCOM(COM); - diff = fabs(COM[whichDirection] - ZPos[i]); - if ( diff <= ztol && states[i] == zcsMoving){ - states[i] = zcsFixed; - changed = true; + force[0] = 0; + force[1] = 0; + force[2] = 0; + + totalFZ_local = 0; + + for (int i = 0; i < (int) (zconsMols.size()); i++){ + if (states[i] == zcsMoving){ + zconsMols[i]->getCOM(COM); + // 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; + + harmonicF = -kz[i] * diff; + totalFZ_local += harmonicF; + + //adjust force + + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + + 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); + } } - else if ( diff > ztol && states[i] == zcsFixed){ - states[i] = zcsMoving; - changed = true; - } - } - return changed; -} +#ifndef IS_MPI + totalFZ = totalFZ_local; +#else + MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif -template bool ZConstraint::haveFixedZMols(){ - for(int i = 0; i < zconsMols.size(); i++) - if (states[i] == zcsFixed) - return true; + //cout << "before substracting harmonic force from moving molecuels " + // << "total force is " << calcTotalForce() << endl; - return false; + force[0] = 0; + force[1] = 0; + force[2] = 0; + + //modify the forces of unconstrained molecules + 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 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(){ - for(int i = 0; i < zconsMols.size(); i++) - if (states[i] == zcsMoving) - return true; + int havingMoving_local; + int havingMoving; - return false; + 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++){ + unconsMols[i]->getCOMvel(COMvel); + MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; + } + + 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]; + } + } + +#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); +#endif + + double vzOfMovingMols; + vzOfMovingMols = MVzOfMovingMols / + (totalMassOfUncons + totalMassOfMovingZMols); + + return vzOfMovingMols; +} + +/** + * + */ + +template double ZConstraint::calcSysCOMVel(){ + double COMvel[3]; + 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 < (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(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); +#endif + + return tempMVz / (totalMassOfUncons + massOfZCons); +} + +/** + * + */ + +template double ZConstraint::calcTotalForce(){ + double force[3]; + double totalForce_local; + double totalForce; + + totalForce_local = 0; + + for (int i = 0; i < nAtoms; i++){ + atoms[i]->getFrc(force); + totalForce_local += force[whichDirection]; + } + +#ifndef IS_MPI + totalForce = totalForce_local; +#else + 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]; -} \ No newline at end of file + 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; + } + } + +}