--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/11 19:40:06 676 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/14 16:16:39 696 @@ -8,70 +8,26 @@ template ZConstraint::ZConstraint(SimIn //get properties from SimInfo GenericData* data; - IndexData* index; + ZConsParaData* zConsParaData; DoubleData* sampleTime; + DoubleData* tolerance; StringData* filename; - - //retrieve index of z-constraint molecules - data = info->getProperty("zconsindex"); - if(!data) { + 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(); - } - else{ - index = dynamic_cast(data); - - if(!index){ + //by default, the direction of constraint is z + // 0 --> x + // 1 --> y + // 2 --> z + whichDirection = 2; - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); - painCave.isFatal = 1; - simError(); - - } - 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(); - - } - - } - - } + //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); + //retrieve sample time of z-contraint - data = info->getProperty("zconstime"); + data = info->getProperty(ZCONSTIME_ID); if(!data) { @@ -99,9 +55,8 @@ template ZConstraint::ZConstraint(SimIn } - //retrieve output filename of z force - data = info->getProperty("zconsfilename"); + data = info->getProperty(ZCONSFILENAME_ID); if(!data) { @@ -130,12 +85,148 @@ template ZConstraint::ZConstraint(SimIn } + + //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); + + if(!tolerance){ + + sprintf( painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); + painCave.isFatal = 1; + simError(); + + } + else{ + this->zconsTol = tolerance->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(); + } + 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->getTotNmol(); +#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 < 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; + 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); +#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(); + + } + } + + }//end if (!zConsParaData) + }//end if (!data) + +// #ifdef IS_MPI update(); #else int searchResult; - double COM[3]; for(int i = 0; i < nMols; i++){ @@ -145,6 +236,9 @@ template ZConstraint::ZConstraint(SimIn 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); } @@ -167,11 +261,33 @@ template ZConstraint::ZConstraint(SimIn simError(); } - for(int i = 0; i < zconsMols.size(); i++) + //determine the states of z-constraint molecules + for(int 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); + } #endif + //get total masss of unconstraint molecules + double totalMassOfUncons_local; + totalMassOfUncons_local = 0; + + for(int 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; @@ -182,10 +298,9 @@ template ZConstraint::ZConstraint(SimIn totNumOfUnconsAtoms = nUnconsAtoms_local; #else MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#endif +#endif - - + // creat zconsWriter fzOut = new ZConsWriter(zconsOutput.c_str()); if(!fzOut){ @@ -217,6 +332,8 @@ template void ZConstraint::update() zconsMols.clear(); massOfZConsMols.clear(); + zPos.clear(); + kz.clear(); unconsMols.clear(); massOfUnconsMols.clear(); @@ -230,6 +347,9 @@ template void ZConstraint::update() 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); @@ -242,6 +362,16 @@ template void ZConstraint::update() } } + + //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); + } + //The reason to declare fz and indexOfZconsMols as pointer to array is // that we want to make the MPI communication simple @@ -294,14 +424,14 @@ template int ZConstraint::isZConstraint 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) + if ((*parameters)[mid].zconsIndex == index) return mid; - else if (indexOfAllZConsMols[mid] > index ) + else if ((*parameters)[mid].zconsIndex > index ) high = mid -1; else low = mid + 1; @@ -310,10 +440,6 @@ template int ZConstraint::isZConstraint return -1; } -/** - * Description: - * Reset the z coordinates - */ template void ZConstraint::integrate(){ //zero out the velocities of center of mass of unconstrained molecules @@ -330,26 +456,37 @@ template void ZConstraint::integrate(){ * * * - */ - - + */ template void ZConstraint::calcForce(int calcPot, int calcStress){ + double zsys; T::calcForce(calcPot, calcStress); if (checkZConsState()) zeroOutVel(); + + zsys = calcZSys(); + cout << "---------------------------------------------------------------------" <getTime() << endl; + cout << "center of mass at z: " << zsys << endl; + //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <doZconstraintForce(); - + //use harmonical poteintial to move the molecules to the specified positions if (haveMovingZMols()) - //this->doHarmonic(); - + this->doHarmonic(); + + //cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; + fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); - + //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() < double ZConstraint::calcZSys() @@ -359,16 +496,11 @@ template double ZConstraint::calcZSys() 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++){ massOfCurMol = molecules[i].getTotalMass(); @@ -376,24 +508,17 @@ template double ZConstraint::calcZSys() 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); -#else + #else totalMass = totalMass_local; totalMZ = totalMZ_local; - totalMassOfUncons = massOfUncons_local; -#endif + #endif double zsys; zsys = totalMZ / totalMass; @@ -421,26 +546,33 @@ template void ZConstraint::zeroOutVel() Atom** fixedZAtoms; double COMvel[3]; double vel[3]; - + //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){ - zconsMols[i]->getCOMvel(COMvel); + 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++){ 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() < void ZConstraint::zeroOutVel() MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; } - for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ - + for(int i = 0; i < zconsMols.size(); i++){ if (states[i] == zcsMoving){ zconsMols[i]->getCOMvel(COMvel); MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; @@ -491,21 +622,23 @@ template void ZConstraint::zeroOutVel() //modify the velocities of moving z-constrained molecuels Atom** movingZAtoms; - for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ + for(int i = 0; i < 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(){ @@ -516,20 +649,16 @@ template void ZConstraint::doZconstrain 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 - zsys = calcZSys(); - cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl; + //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++){ if (states[i] == zcsFixed){ @@ -544,32 +673,61 @@ template void ZConstraint::doZconstrain } totalFZ_local += fz[i]; - cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; + 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; + 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 - nMovingZMols_local = 0; - for(int i = 0; zconsMols.size(); i++){ + int nMovingZAtoms_local; + int nMovingZAtoms; + + nMovingZAtoms_local = 0; + for(int i = 0; i < zconsMols.size(); i++) if(states[i] == zcsMoving) - nMovingZMols_local += massOfZConsMols[i]; - } + nMovingZAtoms_local += zconsMols[i]->getNAtoms(); + #ifdef IS_MPI MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 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; - nMovingZMols = nMovingZMols_local; + nMovingZAtoms = nMovingZAtoms_local; #endif force[0]= 0; force[1]= 0; force[2]= 0; - force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols); + force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); - //modify the velocites of unconstrained molecules + //modify the forces of unconstrained molecules + int accessCount = 0; for(int i = 0; i < unconsMols.size(); i++){ Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); @@ -579,37 +737,19 @@ template void ZConstraint::doZconstrain } - //modify the velocities of moving z-constrained molecules + //modify the forces of moving z-constrained molecules for(int i = 0; i < zconsMols.size(); i++) { - if (states[i] == zcsMoving){ + if (states[i] == zcsMoving){ - Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + 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++) + 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){ - int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); - zconsAtoms = zconsMols[i]->getMyAtoms(); - - for(int j =0; j < nAtomOfCurZConsMol; j++) { - force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; - zconsAtoms[j]->addFrc(force); - } - - } - - } + cout << "after substracting z-constraint force from moving molecuels " + << "total force is " << calcTotalForce() << endl; } @@ -624,12 +764,12 @@ template bool ZConstraint::checkZConsSt for(int i =0; i < zconsMols.size(); i++){ zconsMols[i]->getCOM(COM); - diff = fabs(COM[whichDirection] - ZPos[i]); - if ( diff <= ztol && states[i] == zcsMoving){ + diff = fabs(COM[whichDirection] - zPos[i]); + if ( diff <= zconsTol && states[i] == zcsMoving){ states[i] = zcsFixed; changed = true; } - else if ( diff > ztol && states[i] == zcsFixed){ + else if ( diff > zconsTol && states[i] == zcsFixed){ states[i] = zcsMoving; changed = true; } @@ -658,4 +798,156 @@ template bool ZConstraint::haveMovingZM return false; -} \ No newline at end of file +} + +/** + * + * + */ + +template void ZConstraint::doHarmonic(){ + double force[3]; + double harmonicU; + double harmonicF; + double COM[3]; + double diff; + double totalFZ; + + force[0] = 0; + force[1] = 0; + force[2] = 0; + + totalFZ = 0; + + cout << "Moving Molecules" << endl; + for(int i = 0; i < 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]; + + harmonicU = 0.5 * kz[i] * diff * diff; + info->lrPot += harmonicU; + + harmonicF = - kz[i] * diff; + totalFZ += harmonicF; + + // + force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); + + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + + for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) + movingZAtoms[j]->addFrc(force); + } + + } + + force[0]= 0; + force[1]= 0; + force[2]= 0; + force[whichDirection] = -totalFZ /totNumOfUnconsAtoms; + + //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); + } + +} + +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 = 0; + + for(int i =0 ; i < nMols; i++){ + molecules[i].getCOMvel(COMvel); + tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection]; + } + + double massOfZCons_local; + double massOfZCons; + + massOfZCons_local = 0; + + for(int i = 0; i < massOfZConsMols.size(); i++){ + massOfZCons_local += massOfZConsMols[i]; + } +#ifndef IS_MPI + massOfZCons = massOfZCons_local; +#else + MPI_Allreduce(&massOfZCons_local, &massOfZCons, 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; + +}