--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/28 21:09:47 736 +++ 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, " + "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); } } } @@ -213,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 @@ -227,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]; @@ -279,8 +279,8 @@ template ZConstraint::ZConstraint(SimIn zPos.push_back((*parameters)[searchResult].zPos); // cout << "index: "<< (*parameters)[searchResult].zconsIndex // <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; - kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); + kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); molecules[i].getCOM(COM); } else @@ -304,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); @@ -320,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 @@ -330,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 @@ -344,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(); } @@ -405,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); @@ -420,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); @@ -452,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(); } @@ -506,7 +494,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 +530,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 +540,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 //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]; @@ -609,7 +582,7 @@ template void ZConstraint::calcForce(in #ifdef IS_MPI if(worldRank == 0){ #endif - // cout << "after calcForce, the COMVel of system is " << zSysCOMVel < 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(); @@ -708,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 @@ -722,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]; @@ -749,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++){ @@ -762,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){ @@ -782,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 @@ -798,7 +771,6 @@ template void ZConstraint::doZconstrain Atom** zconsAtoms; double totalFZ; double totalFZ_local; - double COMvel[3]; double COM[3]; double force[3]; @@ -808,8 +780,11 @@ 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++){ + for(int i = 0; i < (int)(zconsMols.size()); i++){ if (states[i] == zcsFixed){ @@ -825,7 +800,7 @@ template void ZConstraint::doZconstrain //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){ @@ -862,15 +837,16 @@ 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; + 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(); @@ -883,7 +859,7 @@ template void ZConstraint::doZconstrain } //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(); @@ -895,10 +871,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; - } /** @@ -921,7 +896,7 @@ 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); @@ -955,12 +930,15 @@ 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(); @@ -971,6 +949,9 @@ template void ZConstraint::doHarmonic() } } + //cout << "after substracting harmonic force from moving molecuels " + // << "total force is " << calcTotalForce() << endl; + } /** @@ -986,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]); @@ -1018,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; @@ -1044,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; @@ -1126,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 @@ -1177,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(); @@ -1187,18 +1168,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){ @@ -1227,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(); @@ -1236,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){