--- trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/13 19:21:53	693
+++ trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/14 16:16:39	696
@@ -456,35 +456,37 @@ template<typename T> void ZConstraint<T>::integrate(){
  *
  *
  *
- */
-
- 
+ */ 
 template<typename T> void ZConstraint<T>::calcForce(int calcPot, int calcStress){
   double zsys;
 
   T::calcForce(calcPot, calcStress);
 
   if (checkZConsState())
-  zeroOutVel();
+    zeroOutVel();
   
   zsys = calcZSys();
   cout << "---------------------------------------------------------------------" <<endl;
-  cout << "current time: " << info->getTime() <<"\tcenter of mass at z: " << zsys << endl;   	 
-  cout << "before calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
-	
+  cout << "current time: " << info->getTime() << 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;
+
   //do zconstraint force; 
   if (haveFixedZMols())
     this->doZconstraintForce();
-
-
-      
+     
   //use harmonical poteintial to move the molecules to the specified positions
   if (haveMovingZMols())
     this->doHarmonic();
-   
+
+  //cout <<	 "after doHarmonic, totalForce is " << calcTotalForce() << endl;
+
   fzOut->writeFZ(info->getTime(), zconsMols.size(),indexOfZConsMols, fz);
-  cout << "after calcForce, the COMVel of unconstraint molecules is " << calcCOMVel() <<endl; 
+  //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()
@@ -494,16 +496,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(); 
@@ -511,24 +508,17 @@ 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
+ #else
   totalMass = totalMass_local;
   totalMZ = totalMZ_local;
-  totalMassOfUncons = massOfUncons_local;
-#endif  
+ #endif  
 
   double zsys;
   zsys = totalMZ / totalMass;
@@ -557,17 +547,6 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
   double COMvel[3];
   double vel[3];
 
-  double tempMass = 0;
-  double tempMVz =0;
-  double tempVz = 0;
-  for(int i = 0; i < nMols; i++){
-    molecules[i].getCOMvel(COMvel);
-    tempMass +=molecules[i].getTotalMass();  
- 	 tempMVz += molecules[i].getTotalMass() * COMvel[whichDirection];
-	 tempVz += COMvel[whichDirection];
-  }
-  cout << "initial velocity of center of mass is " << tempMVz / tempMass << endl;
-
   //zero out the velocities of center of mass of fixed z-constrained molecules
   
   for(int i = 0; i < zconsMols.size(); i++){
@@ -575,7 +554,7 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
     if (states[i] == zcsFixed){	
 
 	   zconsMols[i]->getCOMvel(COMvel);      
-		cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
+		//cout << "before resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
 
       fixedZAtoms = zconsMols[i]->getMyAtoms(); 
 	  
@@ -586,12 +565,13 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
       }
 
 		zconsMols[i]->getCOMvel(COMvel);
-		cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
+		//cout << "after resetting " << indexOfZConsMols[i] <<"'s vz is " << COMvel[whichDirection] << endl;
     }
 	
   }
 
-  	cout << "before resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;	
+  	//cout << "before resetting the COMVel of moving molecules is " << calcMovingMolsCOMVel() <<endl;	
+  	cout << "before resetting the COMVel of sytem is " << calcSysCOMVel() << endl;	
 		   
   // calculate the vz of center of mass of unconstrained molecules and moving z-constrained molecules
   double MVzOfMovingMols_local;
@@ -657,7 +637,7 @@ template<typename T> void ZConstraint<T>::zeroOutVel()
 
  }
 
-	cout << "after resetting the COMVel of unconstraint molecules is " << calcCOMVel() <<endl;
+	cout << "after resetting the COMVel of moving molecules is " << calcSysCOMVel() <<endl;
 
 }
 
@@ -670,9 +650,8 @@ template<typename T> void ZConstraint<T>::doZconstrain
   double COM[3];
   double force[3];
 
-  int nMovingZMols_local;
-  int nMovingZMols;
 
+
   //constrain the molecules which do not reach the specified positions	
     
   //Zero Out the force of z-contrained molecules    
@@ -694,32 +673,61 @@ template<typename T> void ZConstraint<T>::doZconstrain
       }	
       totalFZ_local += fz[i];
 
-      cout << "index: " << indexOfZConsMols[i] <<"\tcurrent zpos: " << COM[whichDirection] << endl;
+      cout << "index: " << indexOfZConsMols[i] 
+				<<"\tcurrent zpos: " << COM[whichDirection] 
+				<< "\tcurrent fz: " <<fz[i] << endl;
 
     }
 	  
   } 
+	
+  // apply negative to fixed z-constrained molecues;
+  force[0]= 0;
+  force[1]= 0;
+  force[2]= 0;
 
+  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;
+        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
-  nMovingZMols_local = 0;
+  int nMovingZAtoms_local;
+  int nMovingZAtoms;
+	
+  nMovingZAtoms_local = 0;
   for(int i = 0; i < zconsMols.size(); i++)
     if(states[i] == zcsMoving)
-	   nMovingZMols_local += massOfZConsMols[i];
+	   nMovingZAtoms_local += zconsMols[i]->getNAtoms();
   
 #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);
+  MPI_Allreduce(&nMovingZAtoms_local, &nMovingZAtoms, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
 #else
   totalFZ = totalFZ_local;
-  nMovingZMols = nMovingZMols_local;
+  nMovingZAtoms = nMovingZAtoms_local;
 #endif
 
   force[0]= 0;
   force[1]= 0;
   force[2]= 0;
-  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZMols);
+  force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
 
   //modify the forces of unconstrained molecules
+  int accessCount = 0;
   for(int i = 0; i < unconsMols.size(); i++){
      
      Atom** unconsAtoms = unconsMols[i]->getMyAtoms();
@@ -731,35 +739,17 @@ 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){
+    if (states[i] == zcsMoving){
 	 	
-     Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     
+      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++)
+        movingZAtoms[j]->addFrc(force);
+    }
   }
-
-  // apply negative to fixed z-constrained molecues;
-  force[0]= 0;
-  force[1]= 0;
-  force[2]= 0;
-
-  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;
-        zconsAtoms[j]->addFrc(force);
-      }
-		
-    }
-	
-  } 
+  cout << "after substracting z-constraint force from moving molecuels "
+		  << "total force is " << calcTotalForce()  << endl;
 
 }
 
@@ -841,9 +831,11 @@ template<typename T> void ZConstraint<T>::doHarmonic()
       harmonicU = 0.5 * kz[i] * diff * diff;  
 		info->lrPot += harmonicU;
 
-      harmonicF =  - kz[i] * diff / zconsMols[i]->getNAtoms();
-		force[whichDirection] = harmonicF;
+      harmonicF =  - kz[i] * diff;
       totalFZ += harmonicF;
+
+       //
+		force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
 		
       Atom** movingZAtoms = zconsMols[i]->getMyAtoms();     
 
@@ -852,7 +844,7 @@ template<typename T> void ZConstraint<T>::doHarmonic()
     }
 
   }
-
+	
   force[0]= 0;
   force[1]= 0;
   force[2]= 0;
@@ -869,7 +861,7 @@ template<typename T> void ZConstraint<T>::doHarmonic()
 
 }
 
-template<typename T> double ZConstraint<T>::calcCOMVel()
+template<typename T> double ZConstraint<T>::calcMovingMolsCOMVel()
 {
   double MVzOfMovingMols_local;
   double MVzOfMovingMols;
@@ -910,23 +902,52 @@ template<typename T> double ZConstraint<T>::calcCOMVel
 }
 
 
-template<typename T> double ZConstraint<T>::calcCOMVel2()
+template<typename T> double ZConstraint<T>::calcSysCOMVel()
 {
   double COMvel[3];
   double tempMVz = 0;
-  int index;
 		
   for(int i =0 ; i < nMols; i++){
-	 index = isZConstraintMol(&molecules[i]);
-    if( index == -1){
-       molecules[i].getCOMvel(COMvel);
-	    tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection];
-    }
-	 else if(states[i] == zcsMoving){
-       zconsMols[index]->getCOMvel(COMvel);
-	    tempMVz += massOfZConsMols[index]*COMvel[whichDirection];    
-	 }
+    molecules[i].getCOMvel(COMvel);
+	 tempMVz += molecules[i].getTotalMass()*COMvel[whichDirection];
   }
+
+  double massOfZCons_local;
+  double massOfZCons;
+
+  massOfZCons_local = 0;
 	
-  return tempMVz /totalMassOfUncons;
+  for(int i = 0; i < massOfZConsMols.size(); i++){
+    massOfZCons_local += massOfZConsMols[i];
+  }
+#ifndef IS_MPI
+  massOfZCons = massOfZCons_local;
+#else
+  MPI_Allreduce(&massOfZCons_local, &massOfZCons, 1, MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
+#endif
+
+  return tempMVz /(totalMassOfUncons + massOfZCons);
 }
+
+template<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;
+
+}