--- trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/11 19:40:06	676
+++ trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/20 14:50:32	702
@@ -2,76 +2,72 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 #include "simError.h"
 #include <cmath>
 template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) 
-                                    : T(theInfo, the_ff), fz(NULL), 
-				      indexOfZConsMols(NULL)
+                                    : T(theInfo, the_ff), fz(NULL), curZPos(NULL),
+                                 indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
 {
 
   //get properties from SimInfo
   GenericData* data;
-  IndexData* index;
+  ZConsParaData* zConsParaData;
   DoubleData* sampleTime;
+  DoubleData* tolerance;
+  StringData* policy;
   StringData* filename; 
+  double COM[3];
+
+  //by default, the direction of constraint is z
+  // 0 --> x
+  // 1 --> y
+  // 2 --> z
+  whichDirection = 2;
+
+  //estimate the force constant of harmonical potential
+  double Kb = 1.986E-3 ; //in kcal/K
   
-  //retrieve index of z-constraint molecules
-  data = info->getProperty("zconsindex");
-  if(!data) {
+  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /2;
+  zForceConst = Kb * info->target_temp /(halfOfLargestBox * halfOfLargestBox);
 
+  //creat force substraction policy
+  data = info->getProperty(ZCONSFORCEPOLICY_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();  
+               "ZConstraint Warning: User does not set force substraction policy, "
+               "average force substraction policy is used\n");
+    painCave.isFatal = 0;
+    simError();      
+
+    forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
   }
   else{
-    index = dynamic_cast<IndexData*>(data);
+    policy = dynamic_cast<StringData*>(data);
     
-    if(!index){
-
+    if(!policy){
       sprintf( painCave.errMsg,
-                 "ZConstraint error: Can not get property from SimInfo\n");
-      painCave.isFatal = 1;
-      simError();  
-    
+                 "ZConstraint Error: Convertion from GenericData to StringData failure, "
+                 "average force substraction policy is used\n");
+      painCave.isFatal = 0;
+      simError();      
+
+      forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
     }
     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){
+      if(policy->getData() == "BYNUMBER")
+        forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
+      else if(policy->getData() == "BYMASS")
+        forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
+      else{
         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(); 
-                 
-      }
-      
+                  "ZConstraint Warning: unknown force substraction policy, "
+                  "average force substraction policy is used\n");
+        painCave.isFatal = 0;
+        simError();      
+      }   
     }
-         
   }
   
+  
   //retrieve sample time of z-contraint 
-  data = info->getProperty("zconstime");
+  data = info->getProperty(ZCONSTIME_ID);
   
   if(!data) {
       
@@ -99,9 +95,8 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 
   }
   
-  
   //retrieve output filename of z force
-  data = info->getProperty("zconsfilename");
+  data = info->getProperty(ZCONSFILENAME_ID);
   if(!data) {
 
       
@@ -114,7 +109,7 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   }
   else{
 
-     filename = dynamic_cast<StringData*>(data);
+    filename = dynamic_cast<StringData*>(data);
     
     if(!filename){
 
@@ -128,14 +123,148 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
       this->zconsOutput = filename->getData();
     }
     
+  }
 
+  //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<DoubleData*>(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<ZConsParaData*>(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() - 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;
-  double COM[3];
       
   for(int i = 0; i < nMols; i++){
     
@@ -145,6 +274,11 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
     
       zconsMols.push_back(&molecules[i]);      
       massOfZConsMols.push_back(molecules[i].getTotalMass());  
+
+      zPos.push_back((*parameters)[searchResult].zPos);
+      cout << "index: "<< (*parameters)[searchResult].zconsIndex 
+             <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; 
+      kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
       
       molecules[i].getCOM(COM);
     }
@@ -158,20 +292,44 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   }
 
   fz = new double[zconsMols.size()];
+  curZPos = new double[zconsMols.size()]; 
   indexOfZConsMols = new int [zconsMols.size()];
 
-  if(!fz || !indexOfZConsMols){
+  if(!fz || !curZPos || !indexOfZConsMols){
     sprintf( painCave.errMsg,
              "Memory allocation failure in class Zconstraint\n");
     painCave.isFatal = 1;
     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;
@@ -181,34 +339,46 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 #ifndef IS_MPI
   totNumOfUnconsAtoms = nUnconsAtoms_local;
 #else
-  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
-#endif
+  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, 
+                      MPI_INT,MPI_SUM, MPI_COMM_WORLD);  
+#endif  
 
-
+  // creat zconsWriter  
+  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);   
   
-  fzOut = new ZConsWriter(zconsOutput.c_str());   
-  
   if(!fzOut){
     sprintf( painCave.errMsg,
              "Memory allocation failure in class Zconstraint\n");
     painCave.isFatal = 1;
     simError();
   }
-  
+
+  forcePolicy->update();
 }
 
 template<typename T> ZConstraint<T>::~ZConstraint()
 {
   if(fz)
     delete[] fz;
+
+  if(curZPos)
+    delete[] curZPos;
   
   if(indexOfZConsMols)
     delete[] indexOfZConsMols;
   
   if(fzOut)
     delete fzOut;
+  
+  if(forcePolicy)
+    delete forcePolicy;
 }
 
+
+/**
+ * 
+ */
+
 #ifdef IS_MPI
 template<typename T> void ZConstraint<T>::update()
 {
@@ -217,6 +387,8 @@ template<typename T> void ZConstraint<T>::update()
   
   zconsMols.clear();
   massOfZConsMols.clear();
+  zPos.clear();
+  kz.clear();
   
   unconsMols.clear();
   massOfUnconsMols.clear();
@@ -230,6 +402,9 @@ template<typename T> void ZConstraint<T>::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);
@@ -242,20 +417,34 @@ template<typename T> void ZConstraint<T>::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
   if(fz)
     delete[] fz;
+  
+  if(curZPos)
+    delete[] curZPos;
     
   if(indexOfZConsMols)
     delete[] indexOfZConsMols;
     
   if (zconsMols.size() > 0){
     fz = new double[zconsMols.size()];
+   curZPos = new double[zconsMols.size()];
     indexOfZConsMols =  new int[zconsMols.size()];
     
-    if(!fz || !indexOfZConsMols){
+    if(!fz || !curZPos || !indexOfZConsMols){
       sprintf( painCave.errMsg,
                "Memory allocation failure in class Zconstraint\n");
       painCave.isFatal = 1;
@@ -269,19 +458,24 @@ template<typename T> void ZConstraint<T>::update()
   }
   else{
     fz = NULL;
+   curZPos = NULL;
     indexOfZConsMols = NULL;
   }
   
+  //
+  forcePolicy->update();
+  
 }
 
 #endif
 
-/**  Function Name: isZConstraintMol
- **  Parameter
- **    Molecule* mol
- **  Return value:
- **    -1, if the molecule is not z-constraint molecule, 
- **    other non-negative values, its index in indexOfAllZConsMols vector
+/**
+ *  Function Name: isZConstraintMol
+ *  Parameter
+ *    Molecule* mol
+ *  Return value:
+ *    -1, if the molecule is not z-constraint molecule, 
+ *    other non-negative values, its index in indexOfAllZConsMols vector
  */
 
 template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
@@ -294,14 +488,14 @@ template<typename T> int ZConstraint<T>::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; 
@@ -310,10 +504,6 @@ template<typename T> int ZConstraint<T>::isZConstraint
   return -1;
 }
 
-/** 
- * Description:
- *  Reset the z coordinates
- */
 template<typename T> void ZConstraint<T>::integrate(){
   
   //zero out the velocities of center of mass of unconstrained molecules 
@@ -330,27 +520,78 @@ template<typename T> void ZConstraint<T>::integrate(){
  *
  *
  *
- */
-
- 
+ */ 
 template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
+  double zsys;
+  double COM[3];
+  double force[3];
+  double zSysCOMVel;
 
   T::calcForce(calcPot, calcStress);
 
-  if (checkZConsState())
+  if (checkZConsState()){
     zeroOutVel();
+   forcePolicy->update();
+  }  
 
+  zsys = calcZSys();
+  zSysCOMVel = calcSysCOMVel();
+#ifdef IS_MPI
+  if(worldRank == 0){
+#endif
+    cout << "---------------------------------------------------------------------" <<endl;
+    cout << "current time: " << info->getTime() << endl;
+    cout << "center of mass at z: " << zsys << endl;     
+    cout << "before calcForce, the COMVel of system is " << zSysCOMVel <<endl;
+
+#ifdef IS_MPI
+  }
+#endif
+
   //do zconstraint force; 
   if (haveFixedZMols())
     this->doZconstraintForce();
-  
+     
   //use harmonical poteintial to move the molecules to the specified positions
   if (haveMovingZMols())
-    //this->doHarmonic();
-   
-  fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
-      
+    this->doHarmonic();
+
+  //write out forces and current positions of z-constraint molecules
+  if(info->getTime() >= curZconsTime){    
+   for(int i = 0; i < zconsMols.size(); i++){
+      zconsMols[i]->getCOM(COM);
+    curZPos[i] = COM[whichDirection];
+
+    //if the z-constraint molecule is still moving, just record its force
+    if(states[i] == zcsMoving){
+         fz[i] = 0;
+      Atom** movingZAtoms;
+      movingZAtoms = zconsMols[i]->getMyAtoms();
+      for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
+           movingZAtoms[j]->getFrc(force);
+           fz[i] += force[whichDirection];
+      }
+     }
+   }
+    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz, curZPos);
+   curZconsTime += zconsTime;
+  }
+
+  zSysCOMVel = calcSysCOMVel();  
+#ifdef IS_MPI
+  if(worldRank == 0){
+#endif
+    cout << "after calcForce, the COMVel of system is " << zSysCOMVel <<endl;
+#ifdef IS_MPI
+  }
+#endif
+
 }
+
+
+/**
+ * 
+ */
  
 template<typename T> double ZConstraint<T>::calcZSys()
 {
@@ -359,16 +600,11 @@ template<typename T> double ZConstraint<T>::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(); 
@@ -376,24 +612,19 @@ template<typename T> double ZConstraint<T>::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
+  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;
-  totalMassOfUncons = massOfUncons_local;
-#endif  
+ #endif  
 
   double zsys;
   zsys = totalMZ / totalMass;
@@ -412,8 +643,6 @@ template<typename T> void ZConstraint<T>::thermalize( 
 
 /**
  *
- *
- *
  */
 
 template<typename T> void ZConstraint<T>::zeroOutVel(){
@@ -421,26 +650,42 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
   Atom** fixedZAtoms;  
   double COMvel[3];
   double vel[3];
-  
+  double zSysCOMVel;
+
   //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){	
+    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() <<endl; 
+
+  zSysCOMVel = calcSysCOMVel();
+#ifdef IS_MPI
+  if(worldRank == 0){
+#endif
+    cout << "before resetting the COMVel of sytem is " << zSysCOMVel << endl;  
+#ifdef IS_MPI
+  }
+#endif
+       
   // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
   double MVzOfMovingMols_local;
   double MVzOfMovingMols;
@@ -455,14 +700,13 @@ template<typename T> void ZConstraint<T>::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];   
-      totalMassOfMovingZMols_local += massOfZConsMols[i]; 	  	
+      totalMassOfMovingZMols_local += massOfZConsMols[i];       
     }
-		
+    
   }
 
 #ifndef IS_MPI
@@ -491,23 +735,37 @@ template<typename T> void ZConstraint<T>::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);
     }
+    
+   }
 
+ }
+
+
+  zSysCOMVel = calcSysCOMVel();
+#ifdef IS_MPI
+  if(worldRank == 0){
+#endif
+    cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
+#ifdef IS_MPI
   }
+#endif
 
 }
 
+/**
+ * 
+ */
+
 template<typename T> void ZConstraint<T>::doZconstraintForce(){
 
   Atom** zconsAtoms;
@@ -516,24 +774,20 @@ template<typename T> void ZConstraint<T>::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;
 
   //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();  
 
@@ -541,110 +795,230 @@ template<typename T> void ZConstraint<T>::doZconstrain
       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;
+      cout << "Fixed Molecule\tindex: " << indexOfZConsMols[i] 
+             <<"\tcurrent zpos: " << COM[whichDirection] 
+             << "\tcurrent fz: " <<fz[i] << 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];
-  }
+  //calculate total z-constraint force
 #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
 
+  
+  // apply negative to fixed z-constrained molecues;
   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 < 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;
+        //force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
+        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(&nMovingZAtoms_local, &nMovingZAtoms, 1, 
+                      MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+#else
+  nMovingZAtoms = nMovingZAtoms_local;
+#endif
+
+  force[0]= 0;
+  force[1]= 0;
+  force[2]= 0;
+
+  //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++)          
+     for(int j = 0; j < unconsMols[i]->getNAtoms(); j++){          
+       force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
+       //force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
        unconsAtoms[j]->addFrc(force);
+     }
      
   }      
 
- //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){
-	 	
-     Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     
+    if (states[i] == zcsMoving){
+    
+      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++){
+        force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
+        //force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
+        movingZAtoms[j]->addFrc(force);
+      }
+    }
   }
 
-  // apply negative to fixed z-constrained molecues;
-  force[0]= 0;
-  force[1]= 0;
-  force[2]= 0;
+  //cout << "after substracting z-constraint force from moving molecuels "
+  //      << "total force is " << calcTotalForce()  << endl;
 
-  for(int i = 0; i < zconsMols.size(); i++){
+}
 
-    if (states[i] == zcsFixed){  
-	
-      int nAtomOfCurZConsMol = zconsMols[i]->getNAtoms();
-      zconsAtoms = zconsMols[i]->getMyAtoms();  
+/**
+  * 
+  *
+  */
+
+template<typename T> void ZConstraint<T>::doHarmonic(){
+  double force[3];
+  double harmonicU;
+  double harmonicF;
+  double COM[3];
+  double diff;
+  double totalFZ_local;
+  double totalFZ;
+  
+  force[0] = 0;
+  force[1] = 0;
+  force[2] = 0;
+
+  totalFZ_local = 0;
+
+  for(int i = 0; i < zconsMols.size(); i++) {
+
+    if (states[i] == zcsMoving){
+      zconsMols[i]->getCOM(COM);
+      cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
     
-      for(int j =0; j < nAtomOfCurZConsMol; j++) {
-        force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
-        zconsAtoms[j]->addFrc(force);
-      }
-		
+    diff = COM[whichDirection] -zPos[i];
+    
+      harmonicU = 0.5 * kz[i] * diff * diff;  
+    info->lrPot += harmonicU;
+
+      harmonicF =  - kz[i] * diff;
+      totalFZ_local += harmonicF;
+
+       //adjust force
+    
+      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     
+
+       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){          
+        force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
+         //force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
+         movingZAtoms[j]->addFrc(force);
+       }
     }
-	
-  } 
 
+  }
+
+#ifndef IS_MPI
+  totalFZ = totalFZ_local;
+#else
+  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);  
+#endif
+
+  force[0]= 0;
+  force[1]= 0;
+  force[2]= 0;
+
+  //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++){          
+       force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
+       //force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
+       unconsAtoms[j]->addFrc(force);     
+     }
+  }   
+
 }
 
+/**
+ * 
+ */
+
 template<typename T> bool ZConstraint<T>::checkZConsState(){
   double COM[3];
   double diff;
   
-  bool changed;
+  int changed_local;
+  int changed;
   
-  changed = false;
+  changed_local = 0;
   
   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;
+     changed_local = 1;
     }
-    else if ( diff > ztol && states[i] == zcsFixed){
+    else if ( diff > zconsTol && states[i] == zcsFixed){
       states[i] = zcsMoving;
-	changed = true;	  
+     changed_local = 1;   
     }
   
   }
 
-  return changed;
+#ifndef IS_MPI
+  changed =changed_local; 
+#else
+  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+  return changed > 0 ? true : false;
 }
 
 template<typename T> bool ZConstraint<T>::haveFixedZMols(){
+
+  int havingFixed_local;
+  int havingFixed;
+
+  havingFixed_local = 0;
+
   for(int i = 0; i < zconsMols.size(); i++)
-    if (states[i] == zcsFixed)
-      return true;
+    if (states[i] == zcsFixed){
+      havingFixed_local = 1;
+    break;
+    }
 
-  return false;
+#ifndef IS_MPI
+  havingFixed = havingFixed_local;
+#else
+  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+  return havingFixed > 0 ? true : false;
 }
 
 
@@ -652,10 +1026,208 @@ template<typename T> bool ZConstraint<T>::haveMovingZM
  *
  */
 template<typename T> bool ZConstraint<T>::haveMovingZMols(){
+
+  int havingMoving_local;
+  int havingMoving;
+
+  havingMoving_local = 0;
+
   for(int i = 0; i < zconsMols.size(); i++)
-    if (states[i] == zcsMoving)
-      return true;
+    if (states[i] == zcsMoving){
+      havingMoving_local = 1;
+    break;
+    }
 
-  return false;
+#ifndef IS_MPI
+  havingMoving = havingMoving_local;
+#else
+  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+  return havingMoving > 0 ? true : false;
   
-}
\ No newline at end of file
+}
+
+/**
+ *
+ */
+
+template<typename T> double ZConstraint<T>::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<typename T> double ZConstraint<T>::calcSysCOMVel()
+{
+  double COMvel[3];
+  double tempMVz_local;
+  double tempMVz;
+  double massOfZCons_local;
+  double massOfZCons;
+
+
+ tempMVz_local = 0;
+
+  for(int i =0 ; i < nMols; i++){
+    molecules[i].getCOMvel(COMvel);
+   tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
+  }
+
+  massOfZCons_local = 0;
+  
+  for(int i = 0; i < massOfZConsMols.size(); i++){
+    massOfZCons_local += massOfZConsMols[i];
+  }
+#ifndef IS_MPI
+  massOfZCons = massOfZCons_local;
+  tempMVz = tempMVz_local;
+#else
+  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&tempMVz_local, &tempMVz, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+  return tempMVz /(totalMassOfUncons + massOfZCons);
+}
+
+/**
+ *
+ */
+
+template<typename T> double ZConstraint<T>::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;
+
+}
+
+/**
+ *
+ */
+
+template<typename T> void ZConstraint<T>::PolicyByNumber::update(){
+  //calculate the number of atoms of moving z-constrained molecules
+  int nMovingZAtoms_local;
+  int nMovingZAtoms;
+  
+  nMovingZAtoms_local = 0;
+  for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
+    if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
+     nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
+  
+#ifdef IS_MPI
+  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+#else
+  nMovingZAtoms = nMovingZAtoms_local;
+#endif
+  totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
+}
+
+template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
+  return totalForce / mol->getNAtoms();
+}
+
+template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfMovingMols(Atom* atom, double totalForce){
+  return totalForce / totNumOfMovingAtoms;
+}
+
+template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
+    return totalForce / mol->getNAtoms();
+}
+
+template<typename T> double ZConstraint<T>::PolicyByNumber::getHFOfUnconsMols(Atom* atom, double totalForce){
+  return totalForce / zconsIntegrator->totNumOfUnconsAtoms;
+}
+
+/**
+ *
+ */
+
+template<typename T> void ZConstraint<T>::PolicyByMass::update(){
+  //calculate the number of atoms of moving z-constrained molecules
+  double massOfMovingZAtoms_local;
+  double massOfMovingZAtoms;
+  
+  massOfMovingZAtoms_local = 0;
+  for(int i = 0; i < (zconsIntegrator->zconsMols).size(); i++)
+    if((zconsIntegrator->states)[i] == (zconsIntegrator->zcsMoving))
+     massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
+  
+#ifdef IS_MPI
+  MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+#else
+  massOfMovingZAtoms = massOfMovingZAtoms_local;
+#endif
+  totMassOfMovingAtoms = massOfMovingZAtoms_local + zconsIntegrator->totalMassOfUncons;
+}
+
+template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
+  return totalForce * atom->getMass() / mol->getTotalMass();
+}
+
+template<typename T> double ZConstraint<T>::PolicyByMass::getZFOfMovingMols( Atom* atom, double totalForce){
+    return totalForce * atom->getMass() / totMassOfMovingAtoms;
+}
+
+template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
+  return totalForce * atom->getMass() / mol->getTotalMass();
+}
+
+template<typename T> double ZConstraint<T>::PolicyByMass::getHFOfUnconsMols(Atom* atom, double totalForce){
+    return totalForce * atom->getMass() / zconsIntegrator->totalMassOfUncons;
+}
+