--- trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/12 16:08:05 681 +++ trunk/OOPSE/libmdtools/ZConstraint.cpp 2003/08/12 17:51:33 682 @@ -8,70 +8,26 @@ template ZConstraint::ZConstraint(SimIn //get properties from SimInfo GenericData* data; - IndexData* index; + ZConsParaData* zConsParaData; DoubleData* sampleTime; + DoubleData* tolerance; StringData* filename; - - //retrieve index of z-constraint molecules - 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(); - - //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 - 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(); - - } - - } - - } + //estimate the force constant of harmonical potential + double Kb = 1.986E-3 ; //in kcal/K + 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"); + data = info->getProperty(ZCONSTIME_ID); if(!data) { @@ -99,9 +55,8 @@ template ZConstraint::ZConstraint(SimIn } - //retrieve output filename of z force - data = info->getProperty("zconsfilename"); + data = info->getProperty(ZCONSFILENAME_ID); if(!data) { @@ -130,12 +85,148 @@ template ZConstraint::ZConstraint(SimIn } + + //retrieve tolerance for z-constraint molecuels + data = info->getProperty(ZCONSTOL_ID); + if(!data) { + + sprintf( painCave.errMsg, + "ZConstraint error: can not get tolerance \n"); + painCave.isFatal = 1; + simError(); + } + else{ + + tolerance = dynamic_cast(data); + + if(!tolerance){ + + sprintf( painCave.errMsg, + "ZConstraint error: Can not get property from SimInfo\n"); + painCave.isFatal = 1; + simError(); + + } + else{ + this->zconsTol = tolerance->getData(); + } + + } + + //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{ + + zConsParaData = dynamic_cast(data); + + if(!zConsParaData){ + + sprintf( painCave.errMsg, + "ZConstraint error: Can not get parameters of zconstraint method from SimInfo\n"); + painCave.isFatal = 1; + simError(); + + } + else{ + + parameters = zConsParaData->getData(); + + //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 + + 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()].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[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[i].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 i = 0; i < nMols; i++) + if (molecules[i].getGlobalIndex() == (*parameters)[i].zconsIndex){ + molecules[i].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; - double COM[3]; for(int i = 0; i < nMols; i++){ @@ -145,6 +236,9 @@ 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); } @@ -184,8 +278,9 @@ template ZConstraint::ZConstraint(SimIn MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); #endif + checkZConsState(); - + // fzOut = new ZConsWriter(zconsOutput.c_str()); if(!fzOut){ @@ -217,6 +312,8 @@ template void ZConstraint::update() zconsMols.clear(); massOfZConsMols.clear(); + zPos.clear(); + kz.clear(); unconsMols.clear(); massOfUnconsMols.clear(); @@ -230,6 +327,9 @@ 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); @@ -294,14 +394,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; @@ -530,6 +630,7 @@ template void ZConstraint::doZconstrain 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){ @@ -624,12 +725,12 @@ template bool ZConstraint::checkZConsSt for(int i =0; i < zconsMols.size(); i++){ zconsMols[i]->getCOM(COM); - diff = fabs(COM[whichDirection] - ZPos[i]); - if ( diff <= ztol && states[i] == zcsMoving){ + diff = fabs(COM[whichDirection] - zPos[i]); + if ( diff <= zconsTol && states[i] == zcsMoving){ states[i] = zcsFixed; changed = true; } - else if ( diff > ztol && states[i] == zcsFixed){ + else if ( diff > zconsTol && states[i] == zcsFixed){ states[i] = zcsMoving; changed = true; } @@ -658,4 +759,43 @@ template bool ZConstraint::haveMovingZM return false; -} \ No newline at end of file +} + +/** + * + * + */ + +template void ZConstraint::doHarmonic(){ + double force[3]; + double harmonicU; + double COM[3]; + double diff; + + force[0] = 0; + force[1] = 0; + force[2] = 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->ltPot += harmonicU; + + force[whichDirection] = - kz[i] * diff / zconsMols[i]->getNAtoms(); + + Atom** movingZAtoms = zconsMols[i]->getMyAtoms(); + + for(int j = 0; j < zconsMols[i]->getNAtoms(); j++) + movingZAtoms[j]->addFrc(force); + } + + } + +}