--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/14 16:16:39 696 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/15 19:24:13 699 @@ -2,8 +2,8 @@ template ZConstraint::ZConstraint(SimIn #include "simError.h" #include template ZConstraint::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) - : T(theInfo, the_ff), fz(NULL), - indexOfZConsMols(NULL) + : T(theInfo, the_ff), fz(NULL), curZPos(NULL), + indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0) { //get properties from SimInfo @@ -11,6 +11,7 @@ template ZConstraint::ZConstraint(SimIn ZConsParaData* zConsParaData; DoubleData* sampleTime; DoubleData* tolerance; + StringData* policy; StringData* filename; double COM[3]; @@ -25,7 +26,47 @@ template ZConstraint::ZConstraint(SimIn double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2; zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox); + + //creat force substraction policy + data = info->getProperty(ZCONSFORCEPOLICY_ID); + if(!data){ + sprintf( painCave.errMsg, + "ZConstraint Warning: User does not set force substraction policy, " + "average force substraction policy is used\n"); + painCave.isFatal = 0; + simError(); + + forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); + } + else{ + policy = dynamic_cast(data); + + if(!policy){ + sprintf( painCave.errMsg, + "ZConstraint Error: Convertion from GenericData to StringData failure, " + "average force substraction policy is used\n"); + painCave.isFatal = 0; + simError(); + forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); + } + else{ + if(policy->getData() == "BYNUMBER") + forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); + else if(policy->getData() == "BYMASS") + forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); + else{ + sprintf( painCave.errMsg, + "ZConstraint Warning: unknown force substraction policy, " + "average force substraction policy is used\n"); + painCave.isFatal = 0; + simError(); + } + } + } + + + //retrieve sample time of z-contraint data = info->getProperty(ZCONSTIME_ID); @@ -238,6 +279,7 @@ template ZConstraint::ZConstraint(SimIn massOfZConsMols.push_back(molecules[i].getTotalMass()); zPos.push_back((*parameters)[searchResult].zPos); + cout << "index: "<< (*parameters)[searchResult].zconsIndex <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); molecules[i].getCOM(COM); @@ -252,9 +294,10 @@ template ZConstraint::ZConstraint(SimIn } fz = new double[zconsMols.size()]; + curZPos = new double[zconsMols.size()]; indexOfZConsMols = new int [zconsMols.size()]; - if(!fz || !indexOfZConsMols){ + if(!fz || !curZPos || !indexOfZConsMols){ sprintf( painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); painCave.isFatal = 1; @@ -284,7 +327,7 @@ template ZConstraint::ZConstraint(SimIn #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 @@ -297,11 +340,11 @@ template ZConstraint::ZConstraint(SimIn #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()); + fzOut = new ZConsWriter(zconsOutput.c_str(), parameters); if(!fzOut){ sprintf( painCave.errMsg, @@ -309,19 +352,26 @@ template ZConstraint::ZConstraint(SimIn painCave.isFatal = 1; simError(); } - + + forcePolicy->update(); } template ZConstraint::~ZConstraint() { if(fz) delete[] fz; + + if(curZPos) + delete[] curZPos; if(indexOfZConsMols) delete[] indexOfZConsMols; if(fzOut) delete fzOut; + + if(forcePolicy) + delete forcePolicy; } #ifdef IS_MPI @@ -377,15 +427,19 @@ template void ZConstraint::update() // that we want to make the MPI communication simple if(fz) delete[] fz; + + if(curZPos) + delete[] curZPos; if(indexOfZConsMols) delete[] indexOfZConsMols; if (zconsMols.size() > 0){ fz = new double[zconsMols.size()]; + curZPos = new double[zconsMols.size()]; indexOfZConsMols = new int[zconsMols.size()]; - if(!fz || !indexOfZConsMols){ + if(!fz || !curZPos || !indexOfZConsMols){ sprintf( painCave.errMsg, "Memory allocation failure in class Zconstraint\n"); painCave.isFatal = 1; @@ -399,8 +453,12 @@ template void ZConstraint::update() } else{ fz = NULL; + curZPos = NULL; indexOfZConsMols = NULL; } + + // + forcePolicy->update(); } @@ -459,12 +517,15 @@ template void ZConstraint::calcForce(in */ template void ZConstraint::calcForce(int calcPot, int calcStress){ double zsys; + double COM[3]; + double force[3]; T::calcForce(calcPot, calcStress); - if (checkZConsState()) + if (checkZConsState()){ zeroOutVel(); - + forcePolicy->update(); + } zsys = calcZSys(); cout << "---------------------------------------------------------------------" <getTime() << endl; @@ -484,7 +545,27 @@ template void ZConstraint::calcForce(in //cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; - fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); + //write out forces and current positions of z-constraint molecules + if(info->getTime() >= curZconsTime){ + for(int i = 0; i < 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, fz, curZPos); + curZconsTime += zconsTime; + } + //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() < double ZConstraint::calcZSys() #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(&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; @@ -571,7 +652,14 @@ template void ZConstraint::zeroOutVel() } //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() < void ZConstraint::zeroOutVel() } +#ifdef IS_MPI + if (worldRank == 0){ +#endif cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() < void ZConstraint::doZconstrain totalFZ_local = 0; //calculate the total z-contrained force of fixed z-contrained molecules - cout << "Fixed Molecules" << endl; + +#ifdef IS_MPI + if (worldRank == 0){ +#endif + cout << "Fixed Molecules" << endl; +#ifdef IS_MPI + } +#endif + for(int i = 0; i < zconsMols.size(); i++){ if (states[i] == zcsFixed){ @@ -680,6 +782,14 @@ template void ZConstraint::doZconstrain } } + + //calculate total z-constraint force +#ifdef IS_MPI + MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#else + totalFZ = totalFZ_local; +#endif + // apply negative to fixed z-constrained molecues; force[0]= 0; @@ -694,7 +804,8 @@ template void ZConstraint::doZconstrain zconsAtoms = zconsMols[i]->getMyAtoms(); for(int j =0; j < nAtomOfCurZConsMol; j++) { - force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + //force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]); zconsAtoms[j]->addFrc(force); } @@ -702,8 +813,9 @@ template void ZConstraint::doZconstrain } - cout << "after zero out z-constraint force on fixed z-constraint molecuels " - << "total force is " << calcTotalForce() << endl; + //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; @@ -714,26 +826,25 @@ template void ZConstraint::doZconstrain 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 force[0]= 0; force[1]= 0; force[2]= 0; - force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); //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++) + 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); + } } @@ -743,48 +854,72 @@ template void ZConstraint::doZconstrain 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; + int changed_local; + int changed; + + changed_local = 0; - 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; + changed_local = 1; } else if ( diff > zconsTol && states[i] == zcsFixed){ states[i] = zcsMoving; - changed = true; + changed_local = 1; } } - return changed; +#ifndef IS_MPI + changed =changed_local; +#else + MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#endif + + return changed > 0 ? true : false; } template bool ZConstraint::haveFixedZMols(){ + + int havingFixed_local; + int havingFixed; + + havingFixed_local = 0; + for(int i = 0; i < zconsMols.size(); i++) - if (states[i] == zcsFixed) - return true; + if (states[i] == zcsFixed){ + havingFixed_local = 1; + break; + } - return false; +#ifndef IS_MPI + havingFixed = havingFixed_local; +#else + MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#endif + + return havingFixed > 0 ? true : false; } @@ -792,11 +927,25 @@ template bool ZConstraint::haveMovingZM * */ template bool ZConstraint::haveMovingZMols(){ + + int havingMoving_local; + int havingMoving; + + havingMoving_local = 0; + for(int i = 0; i < zconsMols.size(); i++) - if (states[i] == zcsMoving) - return true; + if (states[i] == zcsMoving){ + havingMoving_local = 1; + break; + } - return false; +#ifndef IS_MPI + havingMoving = havingMoving_local; +#else + MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#endif + + return havingMoving > 0 ? true : false; } @@ -811,15 +960,24 @@ template void ZConstraint::doHarmonic() 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; +#ifdef IS_MPI + if (worldRank == 0){ +#endif + cout << "Moving Molecules" << endl; +#ifdef IS_MPI + } +#endif + + for(int i = 0; i < zconsMols.size(); i++) { if (states[i] == zcsMoving){ @@ -832,31 +990,41 @@ template void ZConstraint::doHarmonic() info->lrPot += harmonicU; harmonicF = - kz[i] * diff; - totalFZ += harmonicF; + 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++) + 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); + } } } - + +#ifndef IS_MPI + totalFZ = totalFZ_local; +#else + MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + 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++) + for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){ + force[whichDirection] = - totalFZ /totNumOfUnconsAtoms; + //force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ); unconsAtoms[j]->addFrc(force); + } } } @@ -905,16 +1073,19 @@ template double ZConstraint::calcSysCOM template double ZConstraint::calcSysCOMVel() { double COMvel[3]; - double tempMVz = 0; - + 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 += molecules[i].getTotalMass()*COMvel[whichDirection]; + tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection]; } - double massOfZCons_local; - double massOfZCons; - massOfZCons_local = 0; for(int i = 0; i < massOfZConsMols.size(); i++){ @@ -922,8 +1093,10 @@ template double ZConstraint::calcSysCOM } #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); @@ -951,3 +1124,80 @@ template double ZConstraint::calcTotalF 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 < (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 < (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_local + 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; +} +