--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/07/31 15:35:07 658 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/14 16:16:39 696 @@ -1,45 +1,33 @@ #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 GenericData* data; - IndexData* index; + ZConsParaData* zConsParaData; DoubleData* sampleTime; + DoubleData* tolerance; StringData* filename; - - - data = info->getProperty("zconsindex"); - if(!data) { + double COM[3]; - sprintf( painCave.errMsg, - "ZConstraint error: If you use an ZConstraint\n" - " , you must set index of z-constraint molecules.\n"); - painCave.isFatal = 1; - simError(); - } - else{ - index = dynamic_cast(data); - - if(!index){ + //by default, the direction of constraint is z + // 0 --> x + // 1 --> y + // 2 --> z + whichDirection = 2; - sprintf( painCave.errMsg, - "ZConstraint error: Can not get property from SimInfo\n"); - painCave.isFatal = 1; - simError(); - - } - else{ - indexOfAllZConsMols = index->getIndexData(); - } - - } + //estimate the force constant of harmonical potential + double Kb = 1.986E-3 ; //in kcal/K - //retrive sample time of z-contraint - data = info->getProperty("zconstime"); + double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2; + zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox); + + //retrieve sample time of z-contraint + data = info->getProperty(ZCONSTIME_ID); if(!data) { @@ -67,9 +55,8 @@ template ZConstraint::ZConstraint(SimIn } - - //retrive output filename of z force - data = info->getProperty("zconsfilename"); + //retrieve output filename of z force + data = info->getProperty(ZCONSFILENAME_ID); if(!data) { @@ -99,180 +86,148 @@ template ZConstraint::ZConstraint(SimIn } - - //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]; + //retrieve tolerance for z-constraint molecuels + data = info->getProperty(ZCONSTOL_ID); - 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; - } - + if(!data) { + + sprintf( painCave.errMsg, + "ZConstraint error: can not get tolerance \n"); + painCave.isFatal = 1; + simError(); } + else{ - -#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 + tolerance = dynamic_cast(data); + + if(!tolerance){ - 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); + sprintf( painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); + painCave.isFatal = 1; + simError(); + } - + else{ + this->zconsTol = tolerance->getData(); + } + } -#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); - + + //retrieve index of z-constraint molecules + data = info->getProperty(ZCONSPARADATA_ID); + if(!data) { + + sprintf( painCave.errMsg, + "ZConstraint error: If you use an ZConstraint\n" + " , you must set index of z-constraint molecules.\n"); + painCave.isFatal = 1; + simError(); } else{ - int whichMol; - bool done = false; + zConsParaData = dynamic_cast(data); + + if(!zConsParaData){ - while (!done){ - - MPI_Recv(&status, 1, MPI_INT, 0, tag, MPI_COMM_WORLD, &ierr); + sprintf( painCave.errMsg, + "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); + painCave.isFatal = 1; + simError(); - 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; - } - } - - } + else{ + + parameters = zConsParaData->getData(); - //Brocast the allRefZ to slave nodes; - double* allRefZBuf; - int nZConsMols; - nZConsMols = indexOfAllZConsMols.size(); - - allRefZBuf = new double[nZConsMols]; - - if(worldRank == 0){ + //check the range of zconsIndex + //and the minimum value of index is the first one (we already sorted the data) + //the maximum value of index is the last one - 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; + int maxIndex; + int minIndex; + int totalNumMol; + + minIndex = (*parameters)[0].zconsIndex; + if(minIndex < 0){ + sprintf( painCave.errMsg, + "ZConstraint error: index is out of range\n"); + painCave.isFatal = 1; + simError(); + } + + maxIndex = (*parameters)[parameters->size() - 1].zconsIndex; + +#ifndef IS_MPI + totalNumMol = nMols; +#else + totalNumMol = mpiSim->getTotNmol(); +#endif + + if(maxIndex > totalNumMol - 1){ + sprintf( painCave.errMsg, + "ZConstraint error: index is out of range\n"); + painCave.isFatal = 1; + simError(); + } + + //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++){ + + if(!(*parameters)[i].havingZPos){ + +#ifndef IS_MPI + for(int j = 0; j < nMols; j++){ + if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[j].getCOM(COM); + break; + } + } +#else + //query which processor current zconstraint molecule belongs to + int *MolToProcMap; + int whichNode; + double initZPos; + MolToProcMap = mpiSim->getMolToProcMap(); + whichNode = MolToProcMap[(*parameters)[i].zconsIndex]; + + //broadcast the zpos of current z-contraint molecule + //the node which contain this + + if (worldRank == whichNode ){ + + for(int j = 0; j < nMols; j++) + if (molecules[j].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[j].getCOM(COM); + break; + } + + } + + MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD); #endif + + (*parameters)[i].zPos = COM[whichDirection]; - + sprintf( painCave.errMsg, + "ZConstraint warningr: Does not specify zpos for z-constraint molecule " + "initial z coornidate will be used \n"); + painCave.isFatal = 0; + simError(); + + } + } + + }//end if (!zConsParaData) + }//end if (!data) + +// #ifdef IS_MPI update(); #else int searchResult; - - refZ = allRefZ; - + for(int i = 0; i < nMols; i++){ searchResult = isZConstraintMol(&molecules[i]); @@ -281,7 +236,10 @@ template ZConstraint::ZConstraint(SimIn zconsMols.push_back(&molecules[i]); massOfZConsMols.push_back(molecules[i].getTotalMass()); - + + zPos.push_back((*parameters)[searchResult].zPos); + kz.push_back((*parameters)[searchResult]. kRatio * zForceConst); + molecules[i].getCOM(COM); } else @@ -303,11 +261,46 @@ 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; + 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 + + // creat zconsWriter fzOut = new ZConsWriter(zconsOutput.c_str()); if(!fzOut){ @@ -317,7 +310,6 @@ template ZConstraint::ZConstraint(SimIn simError(); } - fzOut->writeRefZ(indexOfAllZConsMols, allRefZ); } template ZConstraint::~ZConstraint() @@ -340,7 +332,8 @@ template void ZConstraint::update() zconsMols.clear(); massOfZConsMols.clear(); - refZ.clear(); + zPos.clear(); + kz.clear(); unconsMols.clear(); massOfUnconsMols.clear(); @@ -354,10 +347,12 @@ template void ZConstraint::update() if(index > -1){ 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); - refZ.push_back(allRefZ[index]); } else { @@ -367,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 @@ -419,14 +424,14 @@ template int ZConstraint::isZConstraint index = mol->getGlobalIndex(); low = 0; - high = indexOfAllZConsMols.size() - 1; + high = parameters->size() - 1; //Binary Search (we have sorted the array) while(low <= high){ mid = (low + high) /2; - if (indexOfAllZConsMols[mid] == index) + if ((*parameters)[mid].zconsIndex == index) return mid; - else if (indexOfAllZConsMols[mid] > index ) + else if ((*parameters)[mid].zconsIndex > index ) high = mid -1; else low = mid + 1; @@ -435,94 +440,514 @@ 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; - */ -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::calcForce(int calcPot, int calcStress){ + double zsys; -template void ZConstraint::resetZ() + 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(); + + //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() { - 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; + //calculate reference z coordinate for z-constraint molecules + double totalMass_local; + double totalMass; + double totalMZ_local; + double totalMZ; + double massOfCurMol; double COM[3]; + + totalMass_local = 0; + totalMZ_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]; + + } + + +#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); + #else + totalMass = totalMass_local; + totalMZ = totalMZ_local; + #endif + double zsys; - Atom** zconsAtoms; + zsys = totalMZ / totalMass; - mzOfZCons = 0; - mzOfUncons = 0; + 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++){ - mzOfZCons += massOfZConsMols[i] * refZ[i]; + + if (states[i] == zcsFixed){ + + 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); + } + + zconsMols[i]->getCOMvel(COMvel); + //cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl; + } + } -#ifdef IS_MPI - MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <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 - 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; + } - cout << "current time: " << info->getTime() <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); + } + + } + + } + + cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() < void ZConstraint::doZconstraintForce(){ + + 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 - cout << "global index: " << zconsMols[i]->getGlobalIndex() << "\tZ: " << COM[2] << "\t"; - deltaZ = zsys + refZ[i] - COM[2]; - cout << "\tdistance: " << COM[2] +deltaZ - zsys; - //update z coordinate - zconsAtoms = zconsMols[i]->getMyAtoms(); - for(int j =0; j < zconsMols[i]->getNAtoms(); j++){ - zconsAtoms[j]->setZ(zconsAtoms[j]->getZ() + deltaZ); - } + //Zero Out the force of z-contrained molecules + totalFZ_local = 0; + + //calculate the total z-contrained force of fixed z-contrained molecules + cout << "Fixed Molecules" << endl; + 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] + << "\tcurrent fz: " <getNAtoms(); + zconsAtoms = zconsMols[i]->getMyAtoms(); - //calculate z constrain force - fz[i] = massOfZConsMols[i]* deltaZ / dt2; - - cout << "\tforce: " << fz[i] << endl; + 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 + 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(&totalFZ_local, &totalFZ, 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; + nMovingZAtoms = nMovingZAtoms_local; +#endif + + force[0]= 0; + force[1]= 0; + force[2]= 0; + force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms); + + //modify the forces of unconstrained molecules + int accessCount = 0; + 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 forces 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); + } } + + cout << "after substracting z-constraint force from moving molecuels " + << "total force is " << calcTotalForce() << endl; +} + +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 <= zconsTol && states[i] == zcsMoving){ + states[i] = zcsFixed; + changed = true; + } + else if ( diff > zconsTol && 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; + +} + +/** + * + * + */ + +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); + cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl; + + diff = COM[whichDirection] -zPos[i]; + + harmonicU = 0.5 * kz[i] * diff * diff; + info->lrPot += harmonicU; + + harmonicF = - kz[i] * diff; + totalFZ += harmonicF; + + // + force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms(); + + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + + for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) + movingZAtoms[j]->addFrc(force); + } + + } + + 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; + +}