--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/26 20:12:51 723 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/09/02 14:30:12 738 @@ -2,7 +2,7 @@ template ZConstraint::ZConstraint(SimIn #include "simError.h" #include template ZConstraint::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) - : T(theInfo, the_ff), fz(NULL), curZPos(NULL), + : T(theInfo, the_ff), fz(NULL), curZPos(NULL), fzOut(NULL), indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0) { @@ -27,16 +27,16 @@ 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 + //creat force Subtraction 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"); + "ZConstraint Warning: User does not set force Subtraction policy, " + "PolicyByMass is used\n"); painCave.isFatal = 0; simError(); - forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this); } else{ policy = dynamic_cast(data); @@ -44,23 +44,24 @@ template ZConstraint::ZConstraint(SimIn if(!policy){ sprintf( painCave.errMsg, "ZConstraint Error: Convertion from GenericData to StringData failure, " - "average force substraction policy is used\n"); + "PolicyByMass is used\n"); painCave.isFatal = 0; simError(); - forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this); } else{ if(policy->getData() == "BYNUMBER") - forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByNumber(this); else if(policy->getData() == "BYMASS") - forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this); else{ sprintf( painCave.errMsg, - "ZConstraint Warning: unknown force substraction policy, " - "average force substraction policy is used\n"); + "ZConstraint Warning: unknown force Subtraction policy, " + "PolicyByMass is used\n"); painCave.isFatal = 0; simError(); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this); } } } @@ -249,7 +250,7 @@ template ZConstraint::ZConstraint(SimIn (*parameters)[i].zPos = COM[whichDirection]; sprintf( painCave.errMsg, - "ZConstraint warningr: Does not specify zpos for z-constraint molecule " + "ZConstraint warning: Does not specify zpos for z-constraint molecule " "initial z coornidate will be used \n"); painCave.isFatal = 0; simError(); @@ -329,7 +330,6 @@ template ZConstraint::ZConstraint(SimIn MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); #endif - //get total number of unconstrained atoms int nUnconsAtoms_local; nUnconsAtoms_local = 0; @@ -343,16 +343,6 @@ template ZConstraint::ZConstraint(SimIn MPI_INT,MPI_SUM, MPI_COMM_WORLD); #endif - // 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(); - } - forcePolicy->update(); } @@ -505,7 +495,17 @@ template void ZConstraint::integrate(){ } 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(); @@ -531,33 +531,8 @@ template void ZConstraint::calcForce(in T::calcForce(calcPot, calcStress); - if (checkZConsState()){ - -#ifdef IS_MPI - if(worldRank == 0){ -#endif -// std::cerr << "\n" -// << "*******************************************\n" -// << " about to call zeroOutVel()\n" -// << "*******************************************\n" -// << "\n"; -#ifdef IS_MPI - } -#endif - zeroOutVel(); - -#ifdef IS_MPI - if(worldRank == 0){ -#endif -// std::cerr << "\n" -// << "*******************************************\n" -// << " finished zeroOutVel()\n" -// << "*******************************************\n" -// << "\n"; -#ifdef IS_MPI - } -#endif - + if (checkZConsState()){ + zeroOutVel(); forcePolicy->update(); } @@ -566,10 +541,10 @@ template void ZConstraint::calcForce(in #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 <getTime() << endl; + //cout << "center of mass at z: " << zsys << endl; + //cout << "before calcForce, the COMVel of system is " << zSysCOMVel < void ZConstraint::calcForce(in //do zconstraint force; if (haveFixedZMols()) this->doZconstraintForce(); - + //use harmonical poteintial to move the molecules to the specified positions if (haveMovingZMols()) this->doHarmonic(); @@ -608,11 +583,10 @@ template void ZConstraint::calcForce(in #ifdef IS_MPI if(worldRank == 0){ #endif - // cout << "after calcForce, the COMVel of system is " << zSysCOMVel < void ZConstraint::zeroOutVel() if (states[i] == zcsFixed){ zconsMols[i]->getCOMvel(COMvel); - //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; + //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; fixedZAtoms = zconsMols[i]->getMyAtoms(); @@ -708,7 +682,7 @@ template void ZConstraint::zeroOutVel() #ifdef IS_MPI if(worldRank == 0){ #endif -// cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl; + //cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl; #ifdef IS_MPI } #endif @@ -782,7 +756,7 @@ template void ZConstraint::zeroOutVel() #ifdef IS_MPI if(worldRank == 0){ #endif -// cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl; + //cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl; #ifdef IS_MPI } #endif @@ -802,14 +776,15 @@ template void ZConstraint::doZconstrain double COM[3]; double force[3]; - - //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 << "before zero out z-constraint force on fixed z-constraint molecuels " + // << "total force is " << calcTotalForce() << endl; for(int i = 0; i < zconsMols.size(); i++){ @@ -825,10 +800,11 @@ template void ZConstraint::doZconstrain } totalFZ_local += fz[i]; -// cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] -// <<"\tcurrent zpos: " << COM[whichDirection] -// << "\tcurrent fz: " < void ZConstraint::doZconstrain zconsAtoms = zconsMols[i]->getMyAtoms(); for(int j =0; j < nAtomOfCurZConsMol; j++) { - force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; - //force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]); + //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]); zconsAtoms[j]->addFrc(force); } @@ -863,25 +839,10 @@ 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; - - 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(&nMovingZAtoms_local, &nMovingZAtoms, 1, - MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else - nMovingZAtoms = nMovingZAtoms_local; -#endif - force[0]= 0; force[1]= 0; force[2]= 0; @@ -892,8 +853,8 @@ template void ZConstraint::doZconstrain 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); + //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); + force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ); unconsAtoms[j]->addFrc(force); } @@ -906,15 +867,14 @@ template void ZConstraint::doZconstrain 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); + //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; } @@ -944,11 +904,11 @@ template void ZConstraint::doHarmonic() zconsMols[i]->getCOM(COM); // cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] // << "\tcurrent zpos: " << COM[whichDirection] << endl; + + diff = COM[whichDirection] -zPos[i]; - diff = COM[whichDirection] -zPos[i]; - harmonicU = 0.5 * kz[i] * diff * diff; - info->lrPot += harmonicU; + info->lrPot += harmonicU; harmonicF = - kz[i] * diff; totalFZ_local += harmonicF; @@ -958,8 +918,8 @@ template void ZConstraint::doHarmonic() 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); + //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); + force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF); movingZAtoms[j]->addFrc(force); } } @@ -972,6 +932,9 @@ template void ZConstraint::doHarmonic() 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; @@ -982,11 +945,14 @@ template void ZConstraint::doHarmonic() 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); + //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; } @@ -1023,8 +989,9 @@ template bool ZConstraint::checkZConsSt #else MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); #endif - + return (changed > 0); + } template bool ZConstraint::haveFixedZMols(){ @@ -1046,7 +1013,7 @@ template bool ZConstraint::haveFixedZMo MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); #endif - return havingFixed > 0 ? true : false; + return (havingFixed > 0); } @@ -1072,7 +1039,7 @@ template bool ZConstraint::haveMovingZM MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); #endif - return havingMoving > 0 ? true : false; + return (havingMoving > 0); } @@ -1203,18 +1170,6 @@ template void ZConstraint::PolicyByNumb nMovingZAtoms = nMovingZAtoms_local; #endif totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms; - -#ifdef IS_MPI - if(worldRank == 0){ -#endif - // std::cerr << "\n" -// << "*******************************************\n" -// << " fiished Policy by numbr()\n" -// << "*******************************************\n" -// << "\n"; -#ifdef IS_MPI - } -#endif } templatedouble ZConstraint::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){ @@ -1252,7 +1207,7 @@ template void ZConstraint::PolicyByMass #else massOfMovingZAtoms = massOfMovingZAtoms_local; #endif - totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons; + totMassOfMovingAtoms = massOfMovingZAtoms + zconsIntegrator->totalMassOfUncons; } template double ZConstraint::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){