--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/20 14:34:04 701 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/10/28 16:03:37 829 @@ -1,9 +1,9 @@ #include "Integrator.hpp" #include "simError.h" -#include +#include template ZConstraint::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) - : T(theInfo, the_ff), fz(NULL), curZPos(NULL), - indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0) + : T(theInfo, the_ff), indexOfZConsMols(NULL), fz(NULL), curZPos(NULL), + fzOut(NULL), curZconsTime(0), forcePolicy(NULL) { //get properties from SimInfo @@ -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); } } } @@ -212,7 +213,7 @@ template ZConstraint::ZConstraint(SimIn //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++){ + for(int i = 0; i < (int)(parameters->size()); i++){ if(!(*parameters)[i].havingZPos){ #ifndef IS_MPI @@ -226,7 +227,7 @@ template ZConstraint::ZConstraint(SimIn //query which processor current zconstraint molecule belongs to int *MolToProcMap; int whichNode; - double initZPos; + MolToProcMap = mpiSim->getMolToProcMap(); whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; @@ -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(); @@ -276,10 +277,10 @@ 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); +// cout << "index: "<< (*parameters)[searchResult].zconsIndex +// <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; + kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); molecules[i].getCOM(COM); } else @@ -303,7 +304,7 @@ template ZConstraint::ZConstraint(SimIn } //determine the states of z-constraint molecules - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); zconsMols[i]->getCOM(COM); @@ -319,7 +320,7 @@ template ZConstraint::ZConstraint(SimIn double totalMassOfUncons_local; totalMassOfUncons_local = 0; - for(int i = 0; i < unconsMols.size(); i++) + for(int i = 0; i < (int)(unconsMols.size()); i++) totalMassOfUncons_local += unconsMols[i]->getTotalMass(); #ifndef IS_MPI @@ -329,11 +330,10 @@ template ZConstraint::ZConstraint(SimIn 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 @@ -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(); } @@ -404,7 +394,6 @@ template void ZConstraint::update() 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); @@ -419,7 +408,7 @@ template void ZConstraint::update() } //determine the states of z-constraint molecules - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ zconsMols[i]->getCOM(COM); if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) states.push_back(zcsFixed); @@ -451,7 +440,7 @@ template void ZConstraint::update() simError(); } - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); } @@ -505,10 +494,22 @@ 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(); + + curZconsTime = zconsTime + info->getTime(); T::integrate(); @@ -525,35 +526,40 @@ template void ZConstraint::calcForce(in double zsys; double COM[3]; double force[3]; + double zSysCOMVel; T::calcForce(calcPot, calcStress); - if (checkZConsState()){ - zeroOutVel(); - forcePolicy->update(); + if (checkZConsState()){ + zeroOutVel(); + forcePolicy->update(); } + zsys = calcZSys(); - cout << "---------------------------------------------------------------------" <getTime() << endl; - cout << "center of mass at z: " << zsys << endl; - //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <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(); - //cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; - //write out forces and current positions of z-constraint molecules if(info->getTime() >= curZconsTime){ - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ zconsMols[i]->getCOM(COM); curZPos[i] = COM[whichDirection]; @@ -571,9 +577,15 @@ template void ZConstraint::calcForce(in fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos); curZconsTime += zconsTime; } - - //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() < void ZConstraint::zeroOutVel() //zero out the velocities of center of mass of fixed z-constrained molecules - for(int i = 0; i < zconsMols.size(); i++){ + 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; + //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; fixedZAtoms = zconsMols[i]->getMyAtoms(); @@ -669,7 +681,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 @@ -683,12 +695,12 @@ template void ZConstraint::zeroOutVel() MVzOfMovingMols_local = 0; totalMassOfMovingZMols_local = 0; - for(int i =0; i < unconsMols.size(); i++){ + for(int i =0; i < (int)(unconsMols.size()); i++){ unconsMols[i]->getCOMvel(COMvel); MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; } - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ if (states[i] == zcsMoving){ zconsMols[i]->getCOMvel(COMvel); MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; @@ -710,7 +722,7 @@ template void ZConstraint::zeroOutVel() //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++){ @@ -723,7 +735,7 @@ template void ZConstraint::zeroOutVel() //modify the velocities of moving z-constrained molecuels Atom** movingZAtoms; - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ if (states[i] ==zcsMoving){ @@ -743,7 +755,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 @@ -759,20 +771,20 @@ template void ZConstraint::doZconstrain Atom** zconsAtoms; double totalFZ; double totalFZ_local; - double COMvel[3]; 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++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ if (states[i] == zcsFixed){ @@ -786,10 +798,11 @@ template void ZConstraint::doZconstrain } totalFZ_local += fz[i]; - cout << "Fixed Molecule --\tindex: " << indexOfZConsMols[i] - <<"\tcurrent zpos: " << COM[whichDirection] - << "\tcurrent fz: " < void ZConstraint::doZconstrain force[1]= 0; force[2]= 0; - for(int i = 0; i < zconsMols.size(); i++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ if (states[i] == zcsFixed){ @@ -815,8 +828,8 @@ template 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); } @@ -825,58 +838,42 @@ template void ZConstraint::doZconstrain } //cout << "after zero out z-constraint force on fixed z-constraint molecuels " - // << "total force is " << calcTotalForce() << endl; + // << "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; //modify the forces of unconstrained molecules - for(int i = 0; i < unconsMols.size(); i++){ + 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); + //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 < zconsMols.size(); i++) { + 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); + //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; - } /** @@ -899,16 +896,17 @@ template void ZConstraint::doHarmonic() totalFZ_local = 0; - for(int i = 0; i < zconsMols.size(); i++) { + 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; +// 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; @@ -918,8 +916,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); } } @@ -932,22 +930,28 @@ 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; //modify the forces of unconstrained molecules - for(int i = 0; i < unconsMols.size(); i++){ + 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); + //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; + } /** @@ -963,7 +967,7 @@ template bool ZConstraint::checkZConsSt changed_local = 0; - for(int i =0; i < zconsMols.size(); i++){ + for(int i =0; i < (int)(zconsMols.size()); i++){ zconsMols[i]->getCOM(COM); diff = fabs(COM[whichDirection] - zPos[i]); @@ -984,7 +988,8 @@ template bool ZConstraint::checkZConsSt MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD); #endif - return changed > 0 ? true : false; + return (changed > 0); + } template bool ZConstraint::haveFixedZMols(){ @@ -994,7 +999,7 @@ template bool ZConstraint::haveFixedZMo havingFixed_local = 0; - for(int i = 0; i < zconsMols.size(); i++) + for(int i = 0; i < (int)(zconsMols.size()); i++) if (states[i] == zcsFixed){ havingFixed_local = 1; break; @@ -1006,7 +1011,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); } @@ -1020,7 +1025,7 @@ template bool ZConstraint::haveMovingZM havingMoving_local = 0; - for(int i = 0; i < zconsMols.size(); i++) + for(int i = 0; i < (int)(zconsMols.size()); i++) if (states[i] == zcsMoving){ havingMoving_local = 1; break; @@ -1032,7 +1037,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); } @@ -1102,7 +1107,7 @@ template double ZConstraint::calcSysCOM massOfZCons_local = 0; - for(int i = 0; i < massOfZConsMols.size(); i++){ + for(int i = 0; i < (int)(massOfZConsMols.size()); i++){ massOfZCons_local += massOfZConsMols[i]; } #ifndef IS_MPI @@ -1153,7 +1158,7 @@ template void ZConstraint::PolicyByNumb int nMovingZAtoms; nMovingZAtoms_local = 0; - for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++) + for(int i = 0; i < (int)((zconsIntegrator->zconsMols).size()); i++) if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)) nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms(); @@ -1191,7 +1196,7 @@ template void ZConstraint::PolicyByMass double massOfMovingZAtoms; massOfMovingZAtoms_local = 0; - for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++) + for(int i = 0; i < (int)((zconsIntegrator->zconsMols).size()); i++) if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving)) massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass(); @@ -1200,7 +1205,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){