--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/12 17:51:33 682 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/13 19:21:53 693 @@ -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); @@ -342,6 +362,16 @@ template void ZConstraint::update() } } + + //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 @@ -434,22 +460,31 @@ template void ZConstraint::calcForce(in template void ZConstraint::calcForce(int calcPot, int calcStress){ + double zsys; T::calcForce(calcPot, calcStress); if (checkZConsState()) - zeroOutVel(); + zeroOutVel(); + + zsys = calcZSys(); + cout << "---------------------------------------------------------------------" <getTime() <<"\tcenter of mass at z: " << zsys << endl; + cout << "before calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <doZconstraintForce(); - + + + //use harmonical poteintial to move the molecules to the specified positions if (haveMovingZMols()) - //this->doHarmonic(); + this->doHarmonic(); fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); - + cout << "after calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() < double ZConstraint::calcZSys() @@ -521,26 +556,43 @@ template void ZConstraint::zeroOutVel() Atom** fixedZAtoms; double COMvel[3]; double vel[3]; - + + double tempMass = 0; + double tempMVz =0; + double tempVz = 0; + for(int i = 0; i < nMols; i++){ + molecules[i].getCOMvel(COMvel); + tempMass +=molecules[i].getTotalMass(); + tempMVz += molecules[i].getTotalMass() * COMvel[whichDirection]; + tempVz += COMvel[whichDirection]; + } + cout << "initial velocity of center of mass is " << tempMVz / tempMass << endl; + //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 unconstraint molecules is " << calcCOMVel() < 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 +642,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 unconstraint molecules is " << calcCOMVel() < void ZConstraint::doZconstraintForce(){ @@ -616,15 +669,11 @@ 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; //Zero Out the force of z-contrained molecules totalFZ_local = 0; @@ -653,10 +702,10 @@ template void ZConstraint::doZconstrain //calculate the number of atoms of moving z-constrained molecules nMovingZMols_local = 0; - for(int i = 0; zconsMols.size(); i++){ + for(int i = 0; i < zconsMols.size(); i++) if(states[i] == zcsMoving) - nMovingZMols_local += massOfZConsMols[i]; - } + nMovingZMols_local += massOfZConsMols[i]; + #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); @@ -670,7 +719,7 @@ template void ZConstraint::doZconstrain force[2]= 0; force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols); - //modify the velocites of unconstrained molecules + //modify the forces of unconstrained molecules for(int i = 0; i < unconsMols.size(); i++){ Atom** unconsAtoms = unconsMols[i]->getMyAtoms(); @@ -680,7 +729,7 @@ 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){ @@ -769,26 +818,32 @@ 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 / zconsMols[i]->getNAtoms(); + force[whichDirection] = harmonicF; + totalFZ += harmonicF; Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); @@ -798,4 +853,80 @@ 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::calcCOMVel() +{ + 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::calcCOMVel2() +{ + double COMvel[3]; + double tempMVz = 0; + int index; + + for(int i =0 ; i < nMols; i++){ + index = isZConstraintMol(&molecules[i]); + if( index == -1){ + molecules[i].getCOMvel(COMvel); + tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection]; + } + else if(states[i] == zcsMoving){ + zconsMols[index]->getCOMvel(COMvel); + tempMVz += massOfZConsMols[index]*COMvel[whichDirection]; + } + } + + return tempMVz /totalMassOfUncons; +}