--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/12 17:51:33 682 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/14 16:16:39 696 @@ -156,7 +156,7 @@ template ZConstraint::ZConstraint(SimIn simError(); } - maxIndex = (*parameters)[parameters->size()].zconsIndex; + maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; #ifndef IS_MPI totalNumMol = nMols; @@ -179,8 +179,8 @@ template ZConstraint::ZConstraint(SimIn #ifndef IS_MPI for(int j = 0; j < nMols; j++){ - if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ - molecules[i].getCOM(COM); + if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[j].getCOM(COM); break; } } @@ -197,9 +197,9 @@ template ZConstraint::ZConstraint(SimIn if (worldRank == whichNode ){ - for(int i = 0; i < nMols; i++) - if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ - molecules[i].getCOM(COM); + for(int j = 0; j < nMols; j++) + if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[j].getCOM(COM); break; } @@ -261,11 +261,33 @@ template ZConstraint::ZConstraint(SimIn simError(); } - for(int i = 0; i < zconsMols.size(); i++) + //determine the states of z-constraint molecules + for(int i = 0; i < zconsMols.size(); i++){ indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); + + zconsMols[i]->getCOM(COM); + if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) + states.push_back(zcsFixed); + else + states.push_back(zcsMoving); + } #endif + + //get total masss of unconstraint molecules + double totalMassOfUncons_local; + totalMassOfUncons_local = 0; + + for(int i = 0; i < unconsMols.size(); i++) + totalMassOfUncons_local += unconsMols[i]->getTotalMass(); + +#ifndef IS_MPI + totalMassOfUncons = totalMassOfUncons_local; +#else + MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + //get total number of unconstrained atoms int nUnconsAtoms_local; nUnconsAtoms_local = 0; @@ -276,11 +298,9 @@ template ZConstraint::ZConstraint(SimIn totNumOfUnconsAtoms = nUnconsAtoms_local; #else MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#endif +#endif - checkZConsState(); - - // + // creat zconsWriter fzOut = new ZConsWriter(zconsOutput.c_str()); if(!fzOut){ @@ -328,8 +348,8 @@ template void ZConstraint::update() zconsMols.push_back(&molecules[i]); zPos.push_back((*parameters)[index].zPos); - kz.push_back((*parameters)[index].kRatio * zForceConst); - + kz.push_back((*parameters)[index].kRatio * zForceConst); + massOfZConsMols.push_back(molecules[i].getTotalMass()); molecules[i].getCOM(COM); @@ -341,7 +361,17 @@ template void ZConstraint::update() massOfUnconsMols.push_back(molecules[i].getTotalMass()); } + } + + //determine the states of z-constraint molecules + for(int i = 0; i < zconsMols.size(); i++){ + zconsMols[i]->getCOM(COM); + if (fabs(zPos[i] - COM[whichDirection]) < zconsTol) + states.push_back(zcsFixed); + else + states.push_back(zcsMoving); } + //The reason to declare fz and indexOfZconsMols as pointer to array is // that we want to make the MPI communication simple @@ -410,10 +440,6 @@ template int ZConstraint::isZConstraint return -1; } -/** - * Description: - * Reset the z coordinates - */ template void ZConstraint::integrate(){ //zero out the velocities of center of mass of unconstrained molecules @@ -430,26 +456,37 @@ template void ZConstraint::integrate(){ * * * - */ - - + */ template void ZConstraint::calcForce(int calcPot, int calcStress){ + double zsys; T::calcForce(calcPot, calcStress); if (checkZConsState()) zeroOutVel(); + + zsys = calcZSys(); + cout << "---------------------------------------------------------------------" <getTime() << endl; + cout << "center of mass at z: " << zsys << endl; + //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <doZconstraintForce(); - + //use harmonical poteintial to move the molecules to the specified positions if (haveMovingZMols()) - //this->doHarmonic(); - + this->doHarmonic(); + + //cout << "after doHarmonic, totalForce is " << calcTotalForce() << endl; + fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); - + //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() < double ZConstraint::calcZSys() @@ -459,16 +496,11 @@ template double ZConstraint::calcZSys() double totalMass; double totalMZ_local; double totalMZ; - double massOfUncons_local; double massOfCurMol; double COM[3]; - + totalMass_local = 0; - totalMass = 0; totalMZ_local = 0; - totalMZ = 0; - massOfUncons_local = 0; - for(int i = 0; i < nMols; i++){ massOfCurMol = molecules[i].getTotalMass(); @@ -476,24 +508,17 @@ template double ZConstraint::calcZSys() totalMass_local += massOfCurMol; totalMZ_local += massOfCurMol * COM[whichDirection]; - - if(isZConstraintMol(&molecules[i]) == -1){ - - massOfUncons_local += massOfCurMol; - } - + } + - #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(&massOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else + #else totalMass = totalMass_local; totalMZ = totalMZ_local; - totalMassOfUncons = massOfUncons_local; -#endif + #endif double zsys; zsys = totalMZ / totalMass; @@ -521,26 +546,33 @@ template void ZConstraint::zeroOutVel() Atom** fixedZAtoms; double COMvel[3]; double vel[3]; - + //zero out the velocities of center of mass of fixed z-constrained molecules for(int i = 0; i < zconsMols.size(); i++){ if (states[i] == zcsFixed){ - zconsMols[i]->getCOMvel(COMvel); + zconsMols[i]->getCOMvel(COMvel); + //cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; + fixedZAtoms = zconsMols[i]->getMyAtoms(); for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ fixedZAtoms[j]->getVel(vel); - vel[whichDirection] -= COMvel[whichDirection]; - fixedZAtoms[j]->setVel(vel); + vel[whichDirection] -= COMvel[whichDirection]; + fixedZAtoms[j]->setVel(vel); } - + + zconsMols[i]->getCOMvel(COMvel); + //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; } } - + + //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() < void ZConstraint::zeroOutVel() MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; } - for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ - + for(int i = 0; i < zconsMols.size(); i++){ if (states[i] == zcsMoving){ zconsMols[i]->getCOMvel(COMvel); MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; @@ -591,21 +622,23 @@ template void ZConstraint::zeroOutVel() //modify the velocities of moving z-constrained molecuels Atom** movingZAtoms; - for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ + for(int i = 0; i < zconsMols.size(); i++){ if (states[i] ==zcsMoving){ movingZAtoms = zconsMols[i]->getMyAtoms(); - for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ + for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ movingZAtoms[j]->getVel(vel); vel[whichDirection] -= vzOfMovingMols; - movingZAtoms[j]->setVel(vel); - } + movingZAtoms[j]->setVel(vel); + } - } + } - } + } + cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() < void ZConstraint::doZconstraintForce(){ @@ -616,15 +649,10 @@ template void ZConstraint::doZconstrain double COMvel[3]; double COM[3]; double force[3]; - double zsys; - int nMovingZMols_local; - int nMovingZMols; - //constrain the molecules which do not reach the specified positions - zsys = calcZSys(); - cout <<"current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl; + //constrain the molecules which do not reach the specified positions //Zero Out the force of z-contrained molecules totalFZ_local = 0; @@ -645,32 +673,61 @@ template void ZConstraint::doZconstrain } totalFZ_local += fz[i]; - cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; + cout << "index: " << indexOfZConsMols[i] + <<"\tcurrent zpos: " << COM[whichDirection] + << "\tcurrent fz: " <getNAtoms(); + zconsAtoms = zconsMols[i]->getMyAtoms(); + + for(int j =0; j < nAtomOfCurZConsMol; j++) { + force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + zconsAtoms[j]->addFrc(force); + } + + } + + } + + 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 - nMovingZMols_local = 0; - for(int i = 0; zconsMols.size(); i++){ + int nMovingZAtoms_local; + int nMovingZAtoms; + + nMovingZAtoms_local = 0; + for(int i = 0; i < zconsMols.size(); i++) if(states[i] == zcsMoving) - nMovingZMols_local += massOfZConsMols[i]; - } + nMovingZAtoms_local += zconsMols[i]->getNAtoms(); + #ifdef IS_MPI MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); - MPI_Allreduce(&nMovingZMols_local, &nMovingZMols, 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; - nMovingZMols = nMovingZMols_local; + nMovingZAtoms = nMovingZAtoms_local; #endif force[0]= 0; force[1]= 0; force[2]= 0; - force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols); + force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); - //modify the velocites of unconstrained molecules + //modify the forces of unconstrained molecules + int accessCount = 0; for(int i = 0; i < unconsMols.size(); i++){ Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); @@ -680,37 +737,19 @@ template void ZConstraint::doZconstrain } - //modify the velocities of moving z-constrained molecules + //modify the forces of moving z-constrained molecules for(int i = 0; i < zconsMols.size(); i++) { - if (states[i] == zcsMoving){ + if (states[i] == zcsMoving){ - Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); - for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) - movingZAtoms[j]->addFrc(force); - } + for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) + movingZAtoms[j]->addFrc(force); + } } - - // apply negative to fixed z-constrained molecues; - force[0]= 0; - force[1]= 0; - force[2]= 0; - - for(int i = 0; i < zconsMols.size(); i++){ - - if (states[i] == zcsFixed){ - int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms(); - zconsAtoms = zconsMols[i]->getMyAtoms(); - - for(int j =0; j < nAtomOfCurZConsMol; j++) { - force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; - zconsAtoms[j]->addFrc(force); - } - - } - - } + cout << "after substracting z-constraint force from moving molecuels " + << "total force is " << calcTotalForce() << endl; } @@ -769,26 +808,34 @@ template void ZConstraint::doHarmonic() template void ZConstraint::doHarmonic(){ double force[3]; double harmonicU; + double harmonicF; double COM[3]; double diff; + double totalFZ; force[0] = 0; force[1] = 0; force[2] = 0; + totalFZ = 0; + cout << "Moving Molecules" << endl; for(int i = 0; i < zconsMols.size(); i++) { if (states[i] == zcsMoving){ - zconsMols[i]->getCOM(COM): + zconsMols[i]->getCOM(COM); cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; diff = COM[whichDirection] -zPos[i]; harmonicU = 0.5 * kz[i] * diff * diff; - info->ltPot += harmonicU; + info->lrPot += harmonicU; - force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms(); + harmonicF = - kz[i] * diff; + totalFZ += harmonicF; + + // + force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); @@ -797,5 +844,110 @@ template void ZConstraint::doHarmonic() } } + + 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++) + unconsAtoms[j]->addFrc(force); + } + } + +template double ZConstraint::calcMovingMolsCOMVel() +{ + double MVzOfMovingMols_local; + double MVzOfMovingMols; + double totalMassOfMovingZMols_local; + double totalMassOfMovingZMols; + double COMvel[3]; + + MVzOfMovingMols_local = 0; + totalMassOfMovingZMols_local = 0; + + for(int i =0; i < unconsMols.size(); i++){ + unconsMols[i]->getCOMvel(COMvel); + MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection]; + } + + for(int i = 0; i < zconsMols.size(); i++){ + + if (states[i] == zcsMoving){ + zconsMols[i]->getCOMvel(COMvel); + MVzOfMovingMols_local += massOfZConsMols[i] * COMvel[whichDirection]; + totalMassOfMovingZMols_local += massOfZConsMols[i]; + } + + } + +#ifndef IS_MPI + MVzOfMovingMols = MVzOfMovingMols_local; + totalMassOfMovingZMols = totalMassOfMovingZMols_local; +#else + MPI_Allreduce(&MVzOfMovingMols_local, &MVzOfMovingMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&totalMassOfMovingZMols_local, &totalMassOfMovingZMols, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + double vzOfMovingMols; + vzOfMovingMols = MVzOfMovingMols / (totalMassOfUncons + totalMassOfMovingZMols); + + return vzOfMovingMols; +} + + +template double ZConstraint::calcSysCOMVel() +{ + double COMvel[3]; + double tempMVz = 0; + + for(int i =0 ; i < nMols; i++){ + molecules[i].getCOMvel(COMvel); + tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection]; + } + + double massOfZCons_local; + double massOfZCons; + + massOfZCons_local = 0; + + for(int i = 0; i < massOfZConsMols.size(); i++){ + massOfZCons_local += massOfZConsMols[i]; + } +#ifndef IS_MPI + massOfZCons = massOfZCons_local; +#else + MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return tempMVz /(totalMassOfUncons + massOfZCons); +} + +template double ZConstraint::calcTotalForce(){ + + double force[3]; + double totalForce_local; + double totalForce; + + totalForce_local = 0; + + for(int i = 0; i < nAtoms; i++){ + atoms[i]->getFrc(force); + totalForce_local += force[whichDirection]; + } + +#ifndef IS_MPI + totalForce = totalForce_local; +#else + MPI_Allreduce(&totalForce_local, &totalForce, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return totalForce; + +}