--- trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/15 19:24:13	699
+++ trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/20 14:34:04	701
@@ -3,7 +3,7 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 #include <cmath>
 template<typename T> ZConstraint<T>::ZConstraint(SimInfo* theInfo, ForceFields* the_ff) 
                                     : T(theInfo, the_ff), fz(NULL), curZPos(NULL),
-				                         indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
+                                 indexOfZConsMols(NULL), forcePolicy(NULL), curZconsTime(0)
 {
 
   //get properties from SimInfo
@@ -40,8 +40,8 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   }
   else{
     policy = dynamic_cast<StringData*>(data);
-		
-	 if(!policy){
+    
+    if(!policy){
       sprintf( painCave.errMsg,
                  "ZConstraint Error: Convertion from GenericData to StringData failure, "
                  "average force substraction policy is used\n");
@@ -49,24 +49,23 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
       simError();      
 
       forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
-	 }
-	 else{
-	 	if(policy->getData() == "BYNUMBER")
-         forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
-	 	else if(policy->getData() == "BYMASS")
-         forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
-		else{
+    }
+    else{
+      if(policy->getData() == "BYNUMBER")
+        forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
+      else if(policy->getData() == "BYMASS")
+        forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
+      else{
         sprintf( painCave.errMsg,
                   "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_ID);
   
@@ -110,7 +109,7 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   }
   else{
 
-     filename = dynamic_cast<StringData*>(data);
+    filename = dynamic_cast<StringData*>(data);
     
     if(!filename){
 
@@ -124,7 +123,6 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
       this->zconsOutput = filename->getData();
     }
     
-
   }
 
   //retrieve tolerance for z-constraint molecuels
@@ -154,7 +152,7 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
     }
 
   }
-	
+  
   //retrieve index of z-constraint molecules
   data = info->getProperty(ZCONSPARADATA_ID);
   if(!data) {
@@ -195,7 +193,7 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
                "ZConstraint error: index is out of range\n");
         painCave.isFatal = 1;
         simError(); 
-      	}
+        }
 
       maxIndex = (*parameters)[parameters->size() - 1].zconsIndex;
 
@@ -216,50 +214,49 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
       //its initial z coordinate  will be used as default
       for(int i = 0; i < parameters->size(); i++){
 
-	      if(!(*parameters)[i].havingZPos){
-
+        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; 
-              }
+          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 
+          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; 
-               }
-				 
-           }
+          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);			  
+          MPI_Bcast(&COM[whichDirection], 1, MPI_DOUBLE_PRECISION, whichNode, MPI_COMM_WORLD);        
 #endif
             
-	         (*parameters)[i].zPos = COM[whichDirection];
+          (*parameters)[i].zPos = COM[whichDirection];
 
-            sprintf( painCave.errMsg,
+          sprintf( painCave.errMsg,
                      "ZConstraint warningr: Does not specify zpos for z-constraint molecule "
                      "initial z coornidate will be used \n");
-            painCave.isFatal = 0;
-            simError();  
-	  
-	      }
-	    }
-			
+           painCave.isFatal = 0;
+           simError();  
+    
+  }
+      }
+      
     }//end if (!zConsParaData)
   }//end  if (!data)
              
@@ -279,8 +276,9 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
       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);
+      cout << "index: "<< (*parameters)[searchResult].zconsIndex 
+             <<"\tzPos = " << (*parameters)[searchResult].zPos << endl; 
+      kz.push_back((*parameters)[searchResult]. kRatio * zForceConst);
       
       molecules[i].getCOM(COM);
     }
@@ -308,11 +306,11 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   for(int i = 0; i < zconsMols.size(); i++){
     indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
 
-	 zconsMols[i]->getCOM(COM);
+    zconsMols[i]->getCOM(COM);
     if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
-		states.push_back(zcsFixed);
-	 else
-		states.push_back(zcsMoving);
+      states.push_back(zcsFixed);
+    else
+      states.push_back(zcsMoving);
   }
   
 #endif 
@@ -327,7 +325,8 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 #ifndef IS_MPI
   totalMassOfUncons = totalMassOfUncons_local;
 #else
-  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
+  MPI_Allreduce(&totalMassOfUncons_local, &totalMassOfUncons, 1,
+                      MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);  
 #endif
 
 
@@ -340,7 +339,8 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 #ifndef IS_MPI
   totNumOfUnconsAtoms = nUnconsAtoms_local;
 #else
-  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);  
+  MPI_Allreduce(&nUnconsAtoms_local, &totNumOfUnconsAtoms, 1, 
+                      MPI_INT,MPI_SUM, MPI_COMM_WORLD);  
 #endif  
 
   // creat zconsWriter  
@@ -369,11 +369,16 @@ template<typename T> ZConstraint<T>::~ZConstraint()
   
   if(fzOut)
     delete fzOut;
-	
+  
   if(forcePolicy)
     delete forcePolicy;
 }
 
+
+/**
+ * 
+ */
+
 #ifdef IS_MPI
 template<typename T> void ZConstraint<T>::update()
 {
@@ -398,8 +403,8 @@ template<typename T> void ZConstraint<T>::update()
     
       zconsMols.push_back(&molecules[i]);      
       zPos.push_back((*parameters)[index].zPos);
-    	kz.push_back((*parameters)[index].kRatio * zForceConst);
-			
+      kz.push_back((*parameters)[index].kRatio * zForceConst);
+      
       massOfZConsMols.push_back(molecules[i].getTotalMass());  
        
       molecules[i].getCOM(COM);
@@ -415,11 +420,11 @@ 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);
+     zconsMols[i]->getCOM(COM);
       if (fabs(zPos[i] - COM[whichDirection]) < zconsTol)
-		  states.push_back(zcsFixed);
-	   else
-		  states.push_back(zcsMoving);
+      states.push_back(zcsFixed);
+     else
+      states.push_back(zcsMoving);
   }
 
     
@@ -427,7 +432,7 @@ template<typename T> void ZConstraint<T>::update()
   // that we want to make the MPI communication simple
   if(fz)
     delete[] fz;
-	
+  
   if(curZPos)
     delete[] curZPos;
     
@@ -436,7 +441,7 @@ template<typename T> void ZConstraint<T>::update()
     
   if (zconsMols.size() > 0){
     fz = new double[zconsMols.size()];
-	 curZPos = new double[zconsMols.size()];
+   curZPos = new double[zconsMols.size()];
     indexOfZConsMols =  new int[zconsMols.size()];
     
     if(!fz || !curZPos || !indexOfZConsMols){
@@ -453,10 +458,10 @@ template<typename T> void ZConstraint<T>::update()
   }
   else{
     fz = NULL;
-	 curZPos = NULL;
+   curZPos = NULL;
     indexOfZConsMols = NULL;
   }
-	
+  
   //
   forcePolicy->update();
   
@@ -464,12 +469,13 @@ template<typename T> void ZConstraint<T>::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)
@@ -524,16 +530,16 @@ template<typename T> void ZConstraint<T>::calcForce(in
 
   if (checkZConsState()){
     zeroOutVel();
-	 forcePolicy->update();
+   forcePolicy->update();
   }  
   zsys = calcZSys();
   cout << "---------------------------------------------------------------------" <<endl;
   cout << "current time: " << info->getTime() << endl;
-  cout << "center of mass at z: " << zsys << endl;   	 
+  cout << "center of mass at z: " << zsys << endl;     
   //cout << "before calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;
   cout << "before calcForce, the COMVel of system is " << calcSysCOMVel() <<endl;
 
-  //cout <<	 "before doZConstraintForce, totalForce is " << calcTotalForce() << endl;
+  //cout <<  "before doZConstraintForce, totalForce is " << calcTotalForce() << endl;
 
   //do zconstraint force; 
   if (haveFixedZMols())
@@ -543,32 +549,37 @@ template<typename T> void ZConstraint<T>::calcForce(in
   if (haveMovingZMols())
     this->doHarmonic();
 
-  //cout <<	 "after doHarmonic, totalForce is " << calcTotalForce() << endl;
+  //cout <<  "after doHarmonic, totalForce is " << calcTotalForce() << endl;
 
   //write out forces and current positions of z-constraint molecules
-  if(info->getTime() >= curZconsTime){		
-	 for(int i = 0; i < zconsMols.size(); i++){
+  if(info->getTime() >= curZconsTime){    
+   for(int i = 0; i < zconsMols.size(); i++){
       zconsMols[i]->getCOM(COM);
-		curZPos[i] = COM[whichDirection];
+    curZPos[i] = COM[whichDirection];
 
-		//if the z-constraint molecule is still moving, just record its force
-		if(states[i] == zcsMoving){
+    //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++){
+      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;
+   curZconsTime += zconsTime;
   }
-	
+  
   //cout << "after calcForce, the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; 
   cout << "after calcForce, the COMVel of system is " << calcSysCOMVel() <<endl; 
 }
+
+
+/**
+ * 
+ */
  
 template<typename T> double ZConstraint<T>::calcZSys()
 {
@@ -579,7 +590,7 @@ template<typename T> double ZConstraint<T>::calcZSys()
   double totalMZ;
   double massOfCurMol;
   double COM[3];
-	
+  
   totalMass_local = 0;
   totalMZ_local = 0;
   
@@ -594,8 +605,10 @@ template<typename T> double ZConstraint<T>::calcZSys()
  
   
 #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(&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;
@@ -618,8 +631,6 @@ template<typename T> void ZConstraint<T>::thermalize( 
 
 /**
  *
- *
- *
  */
 
 template<typename T> void ZConstraint<T>::zeroOutVel(){
@@ -627,40 +638,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);      
-		//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
+     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;
+    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;	
+    //cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl; 
 
+  zSysCOMVel = calcSysCOMVel();
 #ifdef IS_MPI
-  if (worldRank == 0){
+  if(worldRank == 0){
 #endif
-    cout << "before resetting the COMVel of sytem is " << calcSysCOMVel() << endl;	
+    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;
@@ -679,9 +692,9 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
     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
@@ -715,26 +728,32 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
     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){
+  if(worldRank == 0){
 #endif
-	cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl;
+    cout << "after resetting the COMVel of moving molecules is " << zSysCOMVel << endl;  
 #ifdef IS_MPI
   }
 #endif
 
 }
 
+/**
+ * 
+ */
+
 template<typename T> void ZConstraint<T>::doZconstraintForce(){
 
   Atom** zconsAtoms;
@@ -746,25 +765,17 @@ template<typename T> void ZConstraint<T>::doZconstrain
 
 
 
-  //constrain the molecules which do not reach the specified positions	
+  //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
 
-#ifdef IS_MPI
-  if (worldRank == 0){
-#endif
-    cout << "Fixed Molecules" << endl;
-#ifdef IS_MPI
-  }
-#endif
-
   for(int i = 0; i < zconsMols.size(); i++){
-		
+    
     if (states[i] == zcsFixed){
-	  	
+      
       zconsMols[i]->getCOM(COM);
       zconsAtoms = zconsMols[i]->getMyAtoms();  
 
@@ -772,15 +783,15 @@ 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] 
-				<< "\tcurrent fz: " <<fz[i] << endl;
+      cout << "Fixed Molecule --\tindex: " << indexOfZConsMols[i] 
+             <<"\tcurrent zpos: " << COM[whichDirection] 
+             << "\tcurrent fz: " <<fz[i] << endl;
 
     }
-	  
+    
   } 
 
   //calculate total z-constraint force
@@ -790,7 +801,7 @@ template<typename T> void ZConstraint<T>::doZconstrain
   totalFZ = totalFZ_local;
 #endif
 
-	
+  
   // apply negative to fixed z-constrained molecues;
   force[0]= 0;
   force[1]= 0;
@@ -799,34 +810,35 @@ template<typename T> void ZConstraint<T>::doZconstrain
   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] = -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;
+  //       << "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();
+     nMovingZAtoms_local += zconsMols[i]->getNAtoms();
   
 #ifdef IS_MPI
-  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, 
+                      MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
 #else
   nMovingZAtoms = nMovingZAtoms_local;
 #endif
@@ -851,7 +863,7 @@ template<typename T> void ZConstraint<T>::doZconstrain
  //modify the forces of moving z-constrained molecules
   for(int i = 0; i < zconsMols.size(); i++) {
     if (states[i] == zcsMoving){
-	 	
+    
       Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     
 
       for(int j = 0; j < zconsMols[i]->getNAtoms(); j++){
@@ -863,17 +875,92 @@ template<typename T> void ZConstraint<T>::doZconstrain
   }
 
   //cout << "after substracting z-constraint force from moving molecuels "
-  //		  << "total force is " << calcTotalForce()  << endl;
+  //      << "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 < 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
+
+  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;
   
   int changed_local;
   int changed;
-	
+  
   changed_local = 0;
   
   for(int i =0; i < zconsMols.size(); i++){
@@ -882,11 +969,11 @@ template<typename T> bool ZConstraint<T>::checkZConsSt
     diff = fabs(COM[whichDirection] - zPos[i]);  
     if (  diff <= zconsTol && states[i] == zcsMoving){
       states[i] = zcsFixed;
-	   changed_local = 1;
+     changed_local = 1;
     }
     else if ( diff > zconsTol && states[i] == zcsFixed){
       states[i] = zcsMoving;
-	   changed_local = 1;	  
+     changed_local = 1;   
     }
   
   }
@@ -910,7 +997,7 @@ template<typename T> bool ZConstraint<T>::haveFixedZMo
   for(int i = 0; i < zconsMols.size(); i++)
     if (states[i] == zcsFixed){
       havingFixed_local = 1;
-		break;
+    break;
     }
 
 #ifndef IS_MPI
@@ -936,7 +1023,7 @@ template<typename T> bool ZConstraint<T>::haveMovingZM
   for(int i = 0; i < zconsMols.size(); i++)
     if (states[i] == zcsMoving){
       havingMoving_local = 1;
-		break;
+    break;
     }
 
 #ifndef IS_MPI
@@ -950,85 +1037,9 @@ template<typename T> bool ZConstraint<T>::haveMovingZM
 }
 
 /**
-  * 
-  *
-  */
-
-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;
-
-#ifdef IS_MPI
-  if (worldRank == 0){
-#endif
-    cout << "Moving Molecules" << endl;	
-#ifdef IS_MPI
-  }
-#endif
-
-
-  for(int i = 0; i < zconsMols.size(); i++) {
-
-    if (states[i] == zcsMoving){
-      zconsMols[i]->getCOM(COM);
-      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
-		
-		diff = COM[whichDirection] -zPos[i];
-		
-      harmonicU = 0.5 * kz[i] * diff * diff;  
-		info->lrPot += harmonicU;
-
-      harmonicF =  - kz[i] * diff;
-      totalFZ_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> double ZConstraint<T>::calcMovingMolsCOMVel()
 {
   double MVzOfMovingMols_local;
@@ -1050,9 +1061,9 @@ template<typename T> double ZConstraint<T>::calcMoving
     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
@@ -1069,6 +1080,9 @@ template<typename T> double ZConstraint<T>::calcMoving
   return vzOfMovingMols;
 }
 
+/**
+ *
+ */
 
 template<typename T> double ZConstraint<T>::calcSysCOMVel()
 {
@@ -1083,11 +1097,11 @@ template<typename T> double ZConstraint<T>::calcSysCOM
 
   for(int i =0 ; i < nMols; i++){
     molecules[i].getCOMvel(COMvel);
-	 tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
+   tempMVz_local += molecules[i].getTotalMass()*COMvel[whichDirection];
   }
 
   massOfZCons_local = 0;
-	
+  
   for(int i = 0; i < massOfZConsMols.size(); i++){
     massOfZCons_local += massOfZConsMols[i];
   }
@@ -1102,6 +1116,10 @@ template<typename T> double ZConstraint<T>::calcTotalF
   return tempMVz /(totalMassOfUncons + massOfZCons);
 }
 
+/**
+ *
+ */
+
 template<typename T> double ZConstraint<T>::calcTotalForce(){
 
   double force[3];  
@@ -1133,11 +1151,11 @@ template<typename T> void ZConstraint<T>::PolicyByNumb
   //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();
+     nMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getNAtoms();
   
 #ifdef IS_MPI
   MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
@@ -1147,7 +1165,7 @@ template<typename T> double ZConstraint<T>::PolicyByNu
   totNumOfMovingAtoms = nMovingZAtoms + zconsIntegrator->totNumOfUnconsAtoms;
 }
 
-template<typename T> double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
+template<typename T>double ZConstraint<T>::PolicyByNumber::getZFOfFixedZMols(Molecule* mol, Atom* atom, double totalForce){
   return totalForce / mol->getNAtoms();
 }
 
@@ -1171,11 +1189,11 @@ template<typename T> void ZConstraint<T>::PolicyByMass
   //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();
+     massOfMovingZAtoms_local += (zconsIntegrator->zconsMols)[i]->getTotalMass();
   
 #ifdef IS_MPI
   MPI_Allreduce(&massOfMovingZAtoms_local, &massOfMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);