--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/08 17:48:44 671 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/11 19:40:06 676 @@ -1,8 +1,9 @@ #include "Integrator.hpp" #include "simError.h" - +#include template ZConstraint::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) - : T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL) + : T(theInfo, the_ff), fz(NULL), + indexOfZConsMols(NULL) { //get properties from SimInfo \ No newline at end of file @@ -11,7 +12,7 @@ template ZConstraint::ZConstraint(SimIn DoubleData* sampleTime; StringData* filename; - + //retrieve index of z-constraint molecules data = info->getProperty("zconsindex"); if(!data) { \ No newline at end of file @@ -38,8 +39,17 @@ template ZConstraint::ZConstraint(SimIn //the maximum value of index is the last one(we sorted the index data in SimSetup.cpp) int maxIndex; + int minIndex; int totalNumMol; - + + minIndex = indexOfAllZConsMols[0]; + if(minIndex < 0){ + sprintf( painCave.errMsg, + "ZConstraint error: index is out of range\n"); + painCave.isFatal = 1; + simError(); + } + maxIndex = indexOfAllZConsMols[indexOfAllZConsMols.size() - 1]; #ifndef IS_MPI \ No newline at end of file @@ -60,7 +70,7 @@ template ZConstraint::ZConstraint(SimIn } - //retrive sample time of z-contraint + //retrieve sample time of z-contraint data = info->getProperty("zconstime"); if(!data) { \ No newline at end of file @@ -90,7 +100,7 @@ template ZConstraint::ZConstraint(SimIn } - //retrive output filename of z force + //retrieve output filename of z force data = info->getProperty("zconsfilename"); if(!data) { \ No newline at end of file @@ -118,183 +128,15 @@ template ZConstraint::ZConstraint(SimIn this->zconsOutput = filename->getData(); } - - } - - - //calculate reference z coordinate for z-constraint molecules - double totalMass_local; - 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(); - molecules[i].getCOM(COM); - - totalMass_local += massOfCurMol; - totalMZ_local += massOfCurMol * COM[2]; - - 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 - totalMass = totalMass_local; - totalMZ = totalMZ_local; - totalMassOfUncons = massOfUncons_local; -#endif - - double zsys; - zsys = totalMZ / totalMass; - -#ifndef IS_MPI - for(int i = 0; i < nMols; i++){ - - if(isZConstraintMol(&molecules[i]) > -1 ){ - molecules[i].getCOM(COM); - allRefZ.push_back(COM[2] - zsys); - } - - } -#else - - int whichNode; - enum CommType { RequestMolZPos, EndOfRequest} status; - //int status; - double zpos; - int localIndex; - MPI_Status ierr; - int tag = 0; - - if(worldRank == 0){ - - int globalIndexOfCurMol; - int *MolToProcMap; - MolToProcMap = mpiSim->getMolToProcMap(); - - for(int i = 0; i < indexOfAllZConsMols.size(); i++){ - - whichNode = MolToProcMap[indexOfAllZConsMols[i]]; - globalIndexOfCurMol = indexOfAllZConsMols[i]; - - if(whichNode == 0){ - - for(int j = 0; j < nMols; j++) - if(molecules[j].getGlobalIndex() == globalIndexOfCurMol){ - localIndex = j; - break; - } - - molecules[localIndex].getCOM(COM); - allRefZ.push_back(COM[2] - zsys); - - } - else{ - status = RequestMolZPos; - MPI_Send(&status, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD); - MPI_Send(&globalIndexOfCurMol, 1, MPI_INT, whichNode, tag, MPI_COMM_WORLD); - MPI_Recv(&zpos, 1, MPI_DOUBLE_PRECISION, whichNode, tag, MPI_COMM_WORLD, &ierr); - - allRefZ.push_back(zpos - zsys); - - } - - } //End of Request Loop - - //Send ending request message to slave nodes - status = EndOfRequest; - for(int i =1; i < mpiSim->getNumberProcessors(); i++) - MPI_Send(&status, 1, MPI_INT, i, tag, MPI_COMM_WORLD); - - } - else{ - - int whichMol; - bool done = false; - - while (!done){ - - MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr); - - switch (status){ - - case RequestMolZPos : - - MPI_Recv(&whichMol, 1, MPI_INT, 0, tag, MPI_COMM_WORLD,&ierr); - - for(int i = 0; i < nMols; i++) - if(molecules[i].getGlobalIndex() == whichMol){ - localIndex = i; - break; - } - - molecules[localIndex].getCOM(COM); - zpos = COM[2]; - MPI_Send(&zpos, 1, MPI_DOUBLE_PRECISION, 0, tag, MPI_COMM_WORLD); - - break; - - case EndOfRequest : - - done = true; - break; - } - - } - - } - //Brocast the allRefZ to slave nodes; - double* allRefZBuf; - int nZConsMols; - nZConsMols = indexOfAllZConsMols.size(); - - allRefZBuf = new double[nZConsMols]; - - if(worldRank == 0){ - - for(int i = 0; i < nZConsMols; i++) - allRefZBuf[i] = allRefZ[i]; - } - - MPI_Bcast(allRefZBuf, nZConsMols, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD); - - if(worldRank != 0){ - - for(int i = 0; i < nZConsMols; i++) - allRefZ.push_back(allRefZBuf[i]); } - delete[] allRefZBuf; -#endif - - #ifdef IS_MPI update(); #else int searchResult; - - refZ = allRefZ; - + double COM[3]; + for(int i = 0; i < nMols; i++){ searchResult = isZConstraintMol(&molecules[i]); \ No newline at end of file @@ -303,7 +145,7 @@ template ZConstraint::ZConstraint(SimIn zconsMols.push_back(&molecules[i]); massOfZConsMols.push_back(molecules[i].getTotalMass()); - + molecules[i].getCOM(COM); } else \ No newline at end of file @@ -329,6 +171,20 @@ template ZConstraint::ZConstraint(SimIn indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex(); #endif + + //get total number of unconstrained atoms + int nUnconsAtoms_local; + nUnconsAtoms_local = 0; + for(int i = 0; i < unconsMols.size(); i++) + nUnconsAtoms_local += unconsMols[i]->getNAtoms(); + +#ifndef IS_MPI + totNumOfUnconsAtoms = nUnconsAtoms_local; +#else + MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + fzOut = new ZConsWriter(zconsOutput.c_str()); \ No newline at end of file @@ -339,7 +195,6 @@ template ZConstraint::ZConstraint(SimIn simError(); } - fzOut->writeRefZ(indexOfAllZConsMols, allRefZ); } template ZConstraint::~ZConstraint() \ No newline at end of file @@ -362,7 +217,6 @@ template void ZConstraint::update() zconsMols.clear(); massOfZConsMols.clear(); - refZ.clear(); unconsMols.clear(); massOfUnconsMols.clear(); \ No newline at end of file @@ -379,7 +233,6 @@ template void ZConstraint::update() massOfZConsMols.push_back(molecules[i].getTotalMass()); molecules[i].getCOM(COM); - refZ.push_back(allRefZ[index]); } else { \ No newline at end of file @@ -457,94 +310,352 @@ template int ZConstraint::isZConstraint return -1; } -/** Function Name: integrateStep - ** Parameter: - ** int calcPot; - ** int calcStress; - ** Description: - ** Advance One Step. - ** Memo: - ** The best way to implement z-constraint is to override integrateStep - ** Overriding constrainB is not a good choice, since in integrateStep, - ** constrainB is invoked by below line, - ** if(nConstrained) constrainB(); - ** For instance, we would like to apply z-constraint without bond contrain, - ** In that case, if we override constrainB, Z-constrain method will never be executed; +/** + * Description: + * Reset the z coordinates */ -template void ZConstraint::integrateStep( int calcPot, int calcStress ) -{ - T::integrateStep( calcPot, calcStress ); - resetZ(); +template void ZConstraint::integrate(){ - double currZConsTime = 0; - - //write out forces of z constraint - if( info->getTime() >= currZConsTime){ - fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); - } + //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(); + + T::integrate(); + } + -/** Function Name: resetZ - ** Description: - ** Reset the z coordinates +/** + * + * + * + * */ -template void ZConstraint::resetZ() -{ + +template void ZConstraint::calcForce(int calcPot, int calcStress){ - double pos[3]; - double deltaZ; - double mzOfZCons; //total sum of m*z of z-constrain molecules - double mzOfUncons; //total sum of m*z of unconstrain molecuels; - double totalMZOfZCons; - double totalMZOfUncons; - double COM[3]; - double zsys; - Atom** zconsAtoms; + T::calcForce(calcPot, calcStress); - mzOfZCons = 0; - mzOfUncons = 0; + if (checkZConsState()) + zeroOutVel(); + + //do zconstraint force; + if (haveFixedZMols()) + this->doZconstraintForce(); - for(int i = 0; i < zconsMols.size(); i++){ - mzOfZCons += massOfZConsMols[i] * refZ[i]; + //use harmonical poteintial to move the molecules to the specified positions + if (haveMovingZMols()) + //this->doHarmonic(); + + fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); + +} + +template double ZConstraint::calcZSys() +{ + //calculate reference z coordinate for z-constraint molecules + double totalMass_local; + 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(); + molecules[i].getCOM(COM); + + 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 + totalMass = totalMass_local; + totalMZ = totalMZ_local; + totalMassOfUncons = massOfUncons_local; +#endif -#ifdef IS_MPI - MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + double zsys; + zsys = totalMZ / totalMass; + + return zsys; +} + +/** + * + */ +template void ZConstraint::thermalize( void ){ + + T::thermalize(); + zeroOutVel(); +} + +/** + * + * + * + */ + +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); + 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); + } + + } + + } + + // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules + double MVzOfMovingMols_local; + double MVzOfMovingMols; + double totalMassOfMovingZMols_local; + double totalMassOfMovingZMols; + + 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[i]->getNAtoms(); 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 - totalMZOfZCons = mzOfZCons; + 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); + + //modify the velocites of unconstrained molecules + Atom** unconsAtoms; for(int i = 0; i < unconsMols.size(); i++){ - unconsMols[i]->getCOM(COM); - mzOfUncons += massOfUnconsMols[i] * COM[2]; - } -#ifdef IS_MPI - MPI_Allreduce(&mzOfUncons, &totalMZOfUncons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); -#else - totalMZOfUncons = mzOfUncons; -#endif + unconsAtoms = unconsMols[i]->getMyAtoms(); + for(int j = 0; j < unconsMols[i]->getNAtoms();j++){ + unconsAtoms[j]->getVel(vel); + vel[whichDirection] -= vzOfMovingMols; + unconsAtoms[j]->setVel(vel); + } - zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons; - - for(int i = 0; i < zconsMols.size(); i++){ + } + + //modify the velocities of moving z-constrained molecuels + Atom** movingZAtoms; + for(int i = 0; i < zconsMols[i]->getNAtoms(); i++){ + + if (states[i] ==zcsMoving){ - zconsMols[i]->getCOM(COM); + movingZAtoms = zconsMols[i]->getMyAtoms(); + for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){ + movingZAtoms[j]->getVel(vel); + vel[whichDirection] -= vzOfMovingMols; + movingZAtoms[j]->setVel(vel); + } + + } + + } + +} + +template void ZConstraint::doZconstraintForce(){ + + Atom** zconsAtoms; + double totalFZ; + double totalFZ_local; + 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; - deltaZ = zsys + refZ[i] - COM[2]; - //update z coordinate - zconsAtoms = zconsMols[i]->getMyAtoms(); - for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ - zconsAtoms[j]->getPos(pos); - pos[2] += deltaZ; - zconsAtoms[j]->setPos(pos); - } + //Zero Out the force of z-contrained molecules + totalFZ_local = 0; + + //calculate the total z-contrained force of fixed z-contrained molecules + for(int i = 0; i < zconsMols.size(); i++){ + + if (states[i] == zcsFixed){ + + zconsMols[i]->getCOM(COM); + zconsAtoms = zconsMols[i]->getMyAtoms(); + + fz[i] = 0; + for(int j =0; j < zconsMols[i]->getNAtoms(); j++) { + zconsAtoms[j]->getFrc(force); + fz[i] += force[whichDirection]; + } + totalFZ_local += fz[i]; + + cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; + + } + + } + + //calculate the number of atoms of moving z-constrained molecules + nMovingZMols_local = 0; + for(int i = 0; zconsMols.size(); i++){ + if(states[i] == zcsMoving) + 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); +#else + totalFZ = totalFZ_local; + nMovingZMols = nMovingZMols_local; +#endif + + force[0]= 0; + force[1]= 0; + force[2]= 0; + force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols); + + //modify the velocites 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); + + } + + //modify the velocities of moving z-constrained molecules + for(int i = 0; i < zconsMols.size(); i++) { + if (states[i] == zcsMoving){ + + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + + 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(); - //calculate z constrain force - fz[i] = massOfZConsMols[i]* deltaZ / dt2; - + for(int j =0; j < nAtomOfCurZConsMol; j++) { + force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol; + zconsAtoms[j]->addFrc(force); + } + + } + + } + +} + +template bool ZConstraint::checkZConsState(){ + double COM[3]; + double diff; + + bool changed; + + changed = false; + + for(int i =0; i < zconsMols.size(); i++){ + + zconsMols[i]->getCOM(COM); + diff = fabs(COM[whichDirection] - ZPos[i]); + if ( diff <= ztol && states[i] == zcsMoving){ + states[i] = zcsFixed; + changed = true; + } + else if ( diff > ztol && states[i] == zcsFixed){ + states[i] = zcsMoving; + changed = true; + } + } - + return changed; } + +template bool ZConstraint::haveFixedZMols(){ + for(int i = 0; i < zconsMols.size(); i++) + if (states[i] == zcsFixed) + return true; + + return false; +} + + +/** + * + */ +template bool ZConstraint::haveMovingZMols(){ + for(int i = 0; i < zconsMols.size(); i++) + if (states[i] == zcsMoving) + return true; + + return false; + +} \ No newline at end of file