--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/28 21:09:47 736 +++ 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, " + "ZConstraint Warning: User does not set force Subtraction policy, " "PolicyByMass is used\n"); painCave.isFatal = 0; simError(); - forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this); } else{ policy = dynamic_cast(data); @@ -48,20 +48,20 @@ template ZConstraint::ZConstraint(SimIn painCave.isFatal = 0; simError(); - forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(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, " + "ZConstraint Warning: unknown force Subtraction policy, " "PolicyByMass is used\n"); painCave.isFatal = 0; simError(); - forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this); + forcePolicy = (ForceSubtractionPolicy*) new PolicyByMass(this); } } } @@ -330,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 +342,6 @@ template ZConstraint::ZConstraint(SimIn MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, 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(); } @@ -506,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(); @@ -532,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(); } @@ -567,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 #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 @@ -808,6 +782,9 @@ template void ZConstraint::doZconstrain 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,7 +802,7 @@ template void ZConstraint::doZconstrain //cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] // <<"\tcurrent zpos: " << COM[whichDirection] - // << "\tcurrent fz: " < 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; + force[0]= 0; force[1]= 0; @@ -895,10 +873,9 @@ template void ZConstraint::doZconstrain } } } +// 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; - } /** @@ -955,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; @@ -971,6 +951,9 @@ template void ZConstraint::doHarmonic() } } + cout << "after substracting harmonic force from moving molecuels " + << "total force is " << calcTotalForce() << endl; + } /** @@ -1187,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){ @@ -1236,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){