--- trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/08 17:48:44	671
+++ trunk/OOPSE/libmdtools/ZConstraint.cpp	2004/03/16 19:22:56	1091
@@ -1,550 +1,1260 @@
-#include "Integrator.hpp"
-#include "simError.h"
-
-template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) 
-                                    : T(theInfo, the_ff), fz(NULL), indexOfZConsMols(NULL)
-{
-
-  //get properties from SimInfo
-  GenericData* data;
-  IndexData* index;
-  DoubleData* sampleTime;
-  StringData* filename; 
-  
-  
-  data = info->getProperty("zconsindex");
-  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{
-    index = dynamic_cast<IndexData*>(data);
-    
-    if(!index){
-
-      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 totalNumMol;
-      
-      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(); 
-                 
-      }
-      
-    }
-         
-  }
-  
-  //retrive sample time of z-contraint 
-  data = info->getProperty("zconstime");
-  
-  if(!data) {
-      
-    sprintf( painCave.errMsg,
-               "ZConstraint error: If you use an ZConstraint\n"
-               " , you must set sample time.\n");
-    painCave.isFatal = 1;
-    simError();      
-  }
-  else{
-  
-    sampleTime = dynamic_cast<DoubleData*>(data);
-     
-    if(!sampleTime){
-
-      sprintf( painCave.errMsg,
-                 "ZConstraint error: Can not get property from SimInfo\n");
-      painCave.isFatal = 1;
-      simError();  
-       
-    }
-    else{
-      this->zconsTime = sampleTime->getData();
-    }
-
-  }
-  
-  
-  //retrive output filename of z force
-  data = info->getProperty("zconsfilename");
-  if(!data) {
-
-      
-    sprintf( painCave.errMsg,
-               "ZConstraint error: If you use an ZConstraint\n"
-               " , you must set output filename of z-force.\n");
-    painCave.isFatal = 1;
-    simError();  
-
-  }
-  else{
-
-     filename = dynamic_cast<StringData*>(data);
-    
-    if(!filename){
-
-      sprintf( painCave.errMsg,
-                 "ZConstraint error: Can not get property from SimInfo\n");
-      painCave.isFatal = 1;
-      simError();  
-        
-    }
-    else{
-      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;
-
-  for(int i = 0; i < nMols; i++){
-    
-    searchResult = isZConstraintMol(&molecules[i]); 
-    
-    if(searchResult > -1){
-    
-      zconsMols.push_back(&molecules[i]);      
-      massOfZConsMols.push_back(molecules[i].getTotalMass());  
-       
-      molecules[i].getCOM(COM);
-    }
-    else
-    {
-    
-      unconsMols.push_back(&molecules[i]);
-      massOfUnconsMols.push_back(molecules[i].getTotalMass());
-
-    }
-  }
-
-  fz = new double[zconsMols.size()];
-  indexOfZConsMols = new int [zconsMols.size()];
-
-  if(!fz || !indexOfZConsMols){
-    sprintf( painCave.errMsg,
-             "Memory allocation failure in class Zconstraint\n");
-    painCave.isFatal = 1;
-    simError();
-  }
-
-  for(int i = 0; i < zconsMols.size(); i++)
-    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
-  
-#endif 
-  
-  fzOut = new ZConsWriter(zconsOutput.c_str());   
-  
-  if(!fzOut){
-    sprintf( painCave.errMsg,
-             "Memory allocation failure in class Zconstraint\n");
-    painCave.isFatal = 1;
-    simError();
-  }
-  
-  fzOut->writeRefZ(indexOfAllZConsMols, allRefZ);
-}
-
-template<typename T> ZConstraint<T>::~ZConstraint()
-{
-  if(fz)
-    delete[] fz;
-  
-  if(indexOfZConsMols)
-    delete[] indexOfZConsMols;
-  
-  if(fzOut)
-    delete fzOut;
-}
-
-#ifdef IS_MPI
-template<typename T> void ZConstraint<T>::update()
-{
-  double COM[3];
-  int index;
-  
-  zconsMols.clear();
-  massOfZConsMols.clear();
-  refZ.clear();
-  
-  unconsMols.clear();
-  massOfUnconsMols.clear();
- 
-
-  //creat zconsMol and unconsMol lists
-  for(int i = 0; i < nMols; i++){
-    
-    index = isZConstraintMol(&molecules[i]); 
-    
-    if(index > -1){
-    
-      zconsMols.push_back(&molecules[i]);      
-      massOfZConsMols.push_back(molecules[i].getTotalMass());  
-       
-      molecules[i].getCOM(COM);
-      refZ.push_back(allRefZ[index]);      
-    }
-    else
-    {
-    
-      unconsMols.push_back(&molecules[i]);
-      massOfUnconsMols.push_back(molecules[i].getTotalMass());
-
-    }
-  }
-    
-  //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(indexOfZConsMols)
-    delete[] indexOfZConsMols;
-    
-  if (zconsMols.size() > 0){
-    fz = new double[zconsMols.size()];
-    indexOfZConsMols =  new int[zconsMols.size()];
-    
-    if(!fz || !indexOfZConsMols){
-      sprintf( painCave.errMsg,
-               "Memory allocation failure in class Zconstraint\n");
-      painCave.isFatal = 1;
-      simError();
-    }
-        
-    for(int i = 0; i < zconsMols.size(); i++){
-      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
-    }
-
-  }
-  else{
-    fz = NULL;
-    indexOfZConsMols = NULL;
-  }
-  
-}
-
-#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
- */
-
-template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol)
-{
-  int index;
-  int low;
-  int high;
-  int mid;
-
-  index = mol->getGlobalIndex();
-  
-  low = 0;
-  high = indexOfAllZConsMols.size() - 1;
-  
-  //Binary Search (we have sorted the array)  
-  while(low <= high){
-    mid = (low + high) /2;
-    if (indexOfAllZConsMols[mid] == index)
-      return mid;
-    else if (indexOfAllZConsMols[mid] > index )
-       high = mid -1;
-    else     
-      low = mid + 1; 
-  }
-  
-  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<typename T> void ZConstraint<T>::integrateStep( int calcPot, int calcStress )
-{
-  T::integrateStep( calcPot, calcStress );
-  resetZ();
-  
-  double currZConsTime = 0;
-   
-  //write out forces of z constraint
-  if( info->getTime() >= currZConsTime){   
-      fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz); 
-  }     
-}
-
-/** Function Name: resetZ
- ** Description:
- **  Reset the z coordinates
- */
-
-template<typename T> void ZConstraint<T>::resetZ()
-{
-
-  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;
-
-  mzOfZCons = 0;
-  mzOfUncons  = 0;
-  
-  for(int i = 0; i < zconsMols.size(); i++){
-    mzOfZCons += massOfZConsMols[i] * refZ[i];    
-  }
-
-#ifdef IS_MPI
-  MPI_Allreduce(&mzOfZCons, &totalMZOfZCons, 1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
-#else
-  totalMZOfZCons = mzOfZCons;
-#endif
-
-  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  
-  
-  zsys = (totalMZOfZCons + totalMZOfUncons) /totalMassOfUncons; 
- 
-  for(int i = 0; i < zconsMols.size(); i++){   
-   
-    zconsMols[i]->getCOM(COM);
-    
-    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);  
-    }    
-    
-    //calculate z constrain force
-    fz[i] = massOfZConsMols[i]* deltaZ / dt2;
-    
-  }
-
-      
-}
+#include "Integrator.hpp"
+#include "simError.h"
+#include <math.h>
+
+const double INFINITE_TIME = 10e30;
+template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo,
+                                                 ForceFields* the_ff): T(theInfo, the_ff),
+                                                                       indexOfZConsMols(NULL),
+                                                                       fz(NULL),
+                                                                       curZPos(NULL),
+                                                                       fzOut(NULL),
+                                                                       curZconsTime(0),
+                                                                       forcePolicy(NULL),
+                                                                       hasZConsGap(false){
+  //get properties from SimInfo
+  GenericData* data;
+  ZConsParaData* zConsParaData;
+  DoubleData* sampleTime;
+  DoubleData* tolerance;
+  DoubleData* gap;
+  DoubleData* fixtime;
+  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
+
+  double halfOfLargestBox = max(info->boxL[0], max(info->boxL[1], info->boxL[2])) /
+                            2;
+  zForceConst = Kb * info->target_temp / (halfOfLargestBox * halfOfLargestBox);
+
+  //creat force Subtraction policy
+  data = info->getProperty(ZCONSFORCEPOLICY_ID);
+  if (!data){
+    sprintf(painCave.errMsg,
+            "ZConstraint Warning: User does not set force Subtraction policy, "
+            "PolicyByMass is used\n");
+    painCave.isFatal = 0;
+    simError();      
+
+    forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
+  }
+  else{
+    policy = dynamic_cast<StringData*>(data);
+
+    if (!policy){
+      sprintf(painCave.errMsg,
+              "ZConstraint Error: Convertion from GenericData to StringData failure, "
+              "PolicyByMass is used\n");
+      painCave.isFatal = 0;
+      simError();      
+
+      forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
+    }
+    else{
+      if (policy->getData() == "BYNUMBER")
+        forcePolicy = (ForceSubtractionPolicy *) new PolicyByNumber(this);
+      else if (policy->getData() == "BYMASS")
+        forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
+      else{
+        sprintf(painCave.errMsg,
+                "ZConstraint Warning: unknown force Subtraction policy, "
+                "PolicyByMass is used\n");
+        painCave.isFatal = 0;
+        simError();      
+        forcePolicy = (ForceSubtractionPolicy *) new PolicyByMass(this);
+      }
+    }
+  }
+
+
+  //retrieve sample time of z-contraint 
+  data = info->getProperty(ZCONSTIME_ID);
+
+  if (!data){
+    sprintf(painCave.errMsg,
+            "ZConstraint error: If you use an ZConstraint\n"
+            " , you must set sample time.\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+  else{
+    sampleTime = dynamic_cast<DoubleData*>(data);
+
+    if (!sampleTime){
+      sprintf(painCave.errMsg,
+              "ZConstraint error: Can not get property from SimInfo\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    else{
+      this->zconsTime = sampleTime->getData();
+    }
+  }
+
+  //retrieve output filename of z force
+  data = info->getProperty(ZCONSFILENAME_ID);
+  if (!data){
+    sprintf(painCave.errMsg,
+            "ZConstraint error: If you use an ZConstraint\n"
+            " , you must set output filename of z-force.\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+  else{
+    filename = dynamic_cast<StringData*>(data);
+
+    if (!filename){
+      sprintf(painCave.errMsg,
+              "ZConstraint error: Can not get property from SimInfo\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    else{
+      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();
+    }
+  }
+
+  //quick hack here
+  data = info->getProperty(ZCONSGAP_ID);
+
+  if (data){
+    gap = dynamic_cast<DoubleData*>(data);
+  }
+
+  if (!gap){
+    sprintf(painCave.errMsg,
+            "ZConstraint error: Can not get property from SimInfo\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+  else{
+    this->hasZConsGap = true;
+    this->zconsGap = gap->getData();
+  }
+
+  data = info->getProperty(ZCONSFIXTIME_ID);
+
+  if (data){
+    fixtime = dynamic_cast<DoubleData*>(data);
+  }
+
+  if (!fixtime){
+    sprintf(painCave.errMsg,
+            "ZConstraint error: Can not get property from SimInfo\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+  else{
+    this->zconsFixTime = fixtime->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 < (int) (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;
+
+          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, whichNode,
+                    MPI_COMM_WORLD);        
+#endif
+
+          (*parameters)[i].zPos = COM[whichDirection];
+
+          sprintf(painCave.errMsg,
+                  "ZConstraint warning: 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;
+
+  for (int i = 0; i < nMols; i++){
+    searchResult = isZConstraintMol(&molecules[i]); 
+
+    if (searchResult > -1){
+      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);
+    }
+    else{
+      unconsMols.push_back(&molecules[i]);
+      massOfUnconsMols.push_back(molecules[i].getTotalMass());
+    }
+  }
+
+  fz = new double[zconsMols.size()];
+  curZPos = new double[zconsMols.size()]; 
+  indexOfZConsMols = new int [zconsMols.size()];
+
+  if (!fz || !curZPos || !indexOfZConsMols){
+    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+
+  //determine the states of z-constraint molecules
+  for (int i = 0; i < (int) (zconsMols.size()); i++){
+    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
+
+    zconsMols[i]->getCOM(COM);
+    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
+      states.push_back(zcsFixed);
+
+      if (hasZConsGap)
+        endFixTime.push_back(info->getTime() + zconsFixTime);
+    }
+    else{
+      states.push_back(zcsMoving);
+
+      if (hasZConsGap)
+        endFixTime.push_back(INFINITE_TIME);
+    }
+  }
+
+#endif 
+
+  //get total masss of unconstraint molecules
+  double totalMassOfUncons_local;
+  totalMassOfUncons_local = 0;
+
+  for (int i = 0; i < (int) (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 < (int) (unconsMols.size()); i++)
+    nUnconsAtoms_local += unconsMols[i]->getNAtoms();
+
+#ifndef IS_MPI
+  totNumOfUnconsAtoms = nUnconsAtoms_local;
+#else
+  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT, MPI_SUM,
+                MPI_COMM_WORLD);  
+#endif  
+
+  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(){
+  double COM[3];
+  int index;
+
+  zconsMols.clear();
+  massOfZConsMols.clear();
+  zPos.clear();
+  kz.clear();
+
+  unconsMols.clear();
+  massOfUnconsMols.clear();
+
+
+  //creat zconsMol and unconsMol lists
+  for (int i = 0; i < nMols; i++){
+    index = isZConstraintMol(&molecules[i]); 
+
+    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);
+    }
+    else{
+      unconsMols.push_back(&molecules[i]);
+      massOfUnconsMols.push_back(molecules[i].getTotalMass());
+    }
+  }
+
+
+  //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 || !curZPos || !indexOfZConsMols){
+      sprintf(painCave.errMsg,
+              "Memory allocation failure in class Zconstraint\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+
+    for (int i = 0; i < (int) (zconsMols.size()); i++){
+      indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
+    }
+  }
+  else{
+    fz = NULL;
+    curZPos = NULL;
+    indexOfZConsMols = NULL;
+  }
+
+  //determine the states of z-constraint molecules
+  for (int i = 0; i < (int) (zconsMols.size()); i++){
+    zconsMols[i]->getCOM(COM);
+    if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
+      states.push_back(zcsFixed);
+
+      if (hasZConsGap)
+        endFixTime.push_back(info->getTime() + zconsFixTime);
+    }
+    else{
+      states.push_back(zcsMoving);
+
+      if (hasZConsGap)
+        endFixTime.push_back(INFINITE_TIME);
+    }
+  } 
+  //
+  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
+ */
+
+template<typename T> int ZConstraint<T>::isZConstraintMol(Molecule* mol){
+  int index;
+  int low;
+  int high;
+  int mid;
+
+  index = mol->getGlobalIndex();
+
+  low = 0;
+  high = parameters->size() - 1;
+
+  //Binary Search (we have sorted the array)  
+  while (low <= high){
+    mid = (low + high) / 2;
+    if ((*parameters)[mid].zconsIndex == index)
+      return mid;
+    else if ((*parameters)[mid].zconsIndex > index)
+      high = mid - 1;
+    else
+      low = mid + 1;
+  }
+
+  return -1;
+}
+
+template<typename T> void ZConstraint<T>::integrate(){
+  // creat zconsWriter  
+  fzOut = new ZConsWriter(zconsOutput.c_str(), parameters);   
+
+  if (!fzOut){
+    sprintf(painCave.errMsg, "Memory allocation failure in class Zconstraint\n");
+    painCave.isFatal = 1;
+    simError();
+  }
+
+  //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();
+
+  curZconsTime = zconsTime + info->getTime();
+
+  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 (hasZConsGap){
+    updateZPos();
+  }
+
+  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();
+  }
+
+  //write out forces and current positions of z-constraint molecules
+  if (info->getTime() >= curZconsTime){
+    for (int i = 0; i < (int) (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, &zPos[0]);
+    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(){
+  //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;
+  zsys = totalMZ / totalMass;
+
+  return zsys;
+}
+
+/**
+ *
+ */
+template<typename T> void ZConstraint<T>::thermalize(void){
+  T::thermalize();
+  zeroOutVel();
+}
+
+/**
+ *
+ */
+
+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 < (int) (zconsMols.size()); 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;
+    }
+  }
+
+  //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;
+  double totalMassOfMovingZMols_local;
+  double totalMassOfMovingZMols;
+
+  MVzOfMovingMols_local = 0;
+  totalMassOfMovingZMols_local = 0;
+
+  for (int i = 0; i < (int) (unconsMols.size()); i++){
+    unconsMols[i]->getCOMvel(COMvel);
+    MVzOfMovingMols_local += massOfUnconsMols[i] * COMvel[whichDirection];
+  } 
+
+  for (int i = 0; i < (int) (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);
+
+  //modify the velocites of unconstrained molecules  
+  Atom** unconsAtoms;
+  for (int i = 0; i < (int) (unconsMols.size()); i++){
+    unconsAtoms = unconsMols[i]->getMyAtoms();
+    for (int j = 0; j < unconsMols[i]->getNAtoms(); j++){
+      unconsAtoms[j]->getVel(vel);
+      vel[whichDirection] -= vzOfMovingMols;
+      unconsAtoms[j]->setVel(vel);
+    }
+  }  
+
+  //modify the velocities of moving z-constrained molecuels
+  Atom** movingZAtoms;
+  for (int i = 0; i < (int) (zconsMols.size()); i++){
+    if (states[i] == zcsMoving){
+      movingZAtoms = zconsMols[i]->getMyAtoms();
+      for (int j = 0; j < zconsMols[i]->getNAtoms(); j++){
+        movingZAtoms[j]->getVel(vel);
+        vel[whichDirection] -= vzOfMovingMols;
+        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;
+  double totalFZ; 
+  double totalFZ_local;
+  double COM[3];
+  double force[3];
+
+  //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
+
+  //cout << "before zero out z-constraint force on fixed z-constraint molecuels "
+  //       << "total force is " << calcTotalForce() << endl;
+
+  for (int i = 0; i < (int) (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 << "Fixed Molecule\tindex: " << indexOfZConsMols[i] 
+      //      <<"\tcurrent zpos: " << COM[whichDirection] 
+      //      << "\tcurrent fz: " <<fz[i] << endl;
+    }
+  } 
+
+  //calculate total z-constraint force
+#ifdef IS_MPI
+  MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
+#else
+  totalFZ = totalFZ_local;
+#endif
+
+
+  // apply negative to fixed z-constrained molecues;
+  force[0] = 0;
+  force[1] = 0;
+  force[2] = 0;
+
+  for (int i = 0; i < (int) (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;
+
+
+  force[0] = 0;
+  force[1] = 0;
+  force[2] = 0;
+
+  //modify the forces of unconstrained molecules
+  for (int i = 0; i < (int) (unconsMols.size()); i++){
+    Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
+
+    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 forces of moving z-constrained molecules
+  for (int i = 0; i < (int) (zconsMols.size()); i++){
+    if (states[i] == zcsMoving){
+      Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     
+
+      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);
+      }
+    }
+  }
+  //  cout << "after substracting z-constraint force from moving molecuels "
+  //        << "total force is " << calcTotalForce()  << endl;
+}
+
+/**
+  * 
+  *
+  */
+
+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 < (int) (zconsMols.size()); i++){
+    if (states[i] == zcsMoving){
+      zconsMols[i]->getCOM(COM);
+      //       cout << "Moving Molecule\tindex: " << 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_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
+
+  //cout << "before substracting harmonic force from moving molecuels "
+  //      << "total force is " << calcTotalForce()  << endl;
+
+  force[0] = 0;
+  force[1] = 0;
+  force[2] = 0;
+
+  //modify the forces of unconstrained molecules
+  for (int i = 0; i < (int) (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);
+    }
+  }   
+
+  //cout << "after substracting harmonic force from moving molecuels "
+  //      << "total force is " << calcTotalForce()  << endl;
+}
+
+/**
+ * 
+ */
+
+template<typename T> bool ZConstraint<T>::checkZConsState(){
+  double COM[3];
+  double diff;
+
+  int changed_local;
+  int changed;
+
+  changed_local = 0;
+
+  for (int i = 0; i < (int) (zconsMols.size()); i++){
+    zconsMols[i]->getCOM(COM);
+    diff = fabs(COM[whichDirection] - zPos[i]);  
+    if (diff <= zconsTol && states[i] == zcsMoving){
+      states[i] = zcsFixed;
+      changed_local = 1;
+
+      if (hasZConsGap)
+        endFixTime[i] = info->getTime() + zconsFixTime;
+    }
+    else if (diff > zconsTol && states[i] == zcsFixed){
+      states[i] = zcsMoving;
+      changed_local = 1;   
+
+      if (hasZConsGap)
+        endFixTime[i] = INFINITE_TIME;
+    }
+  }
+
+#ifndef IS_MPI
+  changed = changed_local; 
+#else
+  MPI_Allreduce(&changed_local, &changed, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+  return (changed > 0);
+}
+
+template<typename T> bool ZConstraint<T>::haveFixedZMols(){
+  int havingFixed_local;
+  int havingFixed;
+
+  havingFixed_local = 0;
+
+  for (int i = 0; i < (int) (zconsMols.size()); i++)
+    if (states[i] == zcsFixed){
+      havingFixed_local = 1;
+      break;
+    }
+
+#ifndef IS_MPI
+  havingFixed = havingFixed_local;
+#else
+  MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM,
+                MPI_COMM_WORLD);
+#endif
+
+  return (havingFixed > 0);
+}
+
+
+/**
+ *
+ */
+template<typename T> bool ZConstraint<T>::haveMovingZMols(){
+  int havingMoving_local;
+  int havingMoving;
+
+  havingMoving_local = 0;
+
+  for (int i = 0; i < (int) (zconsMols.size()); i++)
+    if (states[i] == zcsMoving){
+      havingMoving_local = 1;
+      break;
+    }
+
+#ifndef IS_MPI
+  havingMoving = havingMoving_local;
+#else
+  MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM,
+                MPI_COMM_WORLD);
+#endif
+
+  return (havingMoving > 0);
+}
+
+/**
+ *
+ */
+
+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 < (int) (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 < (int) ((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 < (int) ((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 +
+                         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;
+}
+
+template<typename T> void ZConstraint<T>::updateZPos(){
+  double curTime;
+
+  curTime = info->getTime();
+
+  for (size_t i = 0; i < zconsMols.size(); i++){
+    if (states[i] == zcsFixed && curTime >= endFixTime[i]){
+      zPos[i] += zconsGap;
+    }
+  }
+}