--- trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/26 20:12:51	723
+++ trunk/OOPSE/libmdtools/ZConstraint.cpp	2003/08/27 19:23:29	733
@@ -32,11 +32,11 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   if(!data){
     sprintf( painCave.errMsg,
                "ZConstraint Warning: User does not set force substraction policy, "
-               "average force substraction policy is used\n");
+               "PolicyByMass is used\n");
     painCave.isFatal = 0;
     simError();      
 
-    forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
+    forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
   }
   else{
     policy = dynamic_cast<StringData*>(data);
@@ -44,11 +44,11 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
     if(!policy){
       sprintf( painCave.errMsg,
                  "ZConstraint Error: Convertion from GenericData to StringData failure, "
-                 "average force substraction policy is used\n");
+                 "PolicyByMass is used\n");
       painCave.isFatal = 0;
       simError();      
 
-      forcePolicy = (ForceSubstractionPolicy*) new PolicyByNumber(this);
+      forcePolicy = (ForceSubstractionPolicy*) new PolicyByMass(this);
     }
     else{
       if(policy->getData() == "BYNUMBER")
@@ -578,7 +578,7 @@ template<typename T> void ZConstraint<T>::calcForce(in
   //do zconstraint force; 
   if (haveFixedZMols())
     this->doZconstraintForce();
-     
+
   //use harmonical poteintial to move the molecules to the specified positions
   if (haveMovingZMols())
     this->doHarmonic();
@@ -612,7 +612,6 @@ template<typename T> void ZConstraint<T>::calcForce(in
 #ifdef IS_MPI
   }
 #endif
-
 }
 
 
@@ -802,8 +801,6 @@ template<typename T> void ZConstraint<T>::doZconstrain
   double COM[3];
   double force[3];
 
-
-
   //constrain the molecules which do not reach the specified positions  
     
   //Zero Out the force of z-contrained molecules    
@@ -825,10 +822,11 @@ template<typename T> void ZConstraint<T>::doZconstrain
       } 
       totalFZ_local += fz[i];
 
-//       cout << "Fixed Molecule\tindex: " << 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;
 
+
     }
     
   } 
@@ -854,8 +852,8 @@ template<typename T> void ZConstraint<T>::doZconstrain
       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]);
+        //force[whichDirection] = -fz[i]/ nAtomOfCurZConsMol;
+        force[whichDirection] = - forcePolicy->getZFOfFixedZMols(zconsMols[i], zconsAtoms[j], fz[i]);
         zconsAtoms[j]->addFrc(force);
       }
     
@@ -866,22 +864,6 @@ template<typename T> void ZConstraint<T>::doZconstrain
 //   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;
@@ -892,8 +874,8 @@ template<typename T> void ZConstraint<T>::doZconstrain
      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);
+       //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
+       force[whichDirection] = forcePolicy->getZFOfMovingMols(unconsAtoms[j],totalFZ);
        unconsAtoms[j]->addFrc(force);
      }
      
@@ -906,8 +888,8 @@ template<typename T> void ZConstraint<T>::doZconstrain
       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);
+        //force[whichDirection] = totalFZ / (totNumOfUnconsAtoms + nMovingZAtoms);
+        force[whichDirection] = forcePolicy->getZFOfMovingMols(movingZAtoms[j],totalFZ);
         movingZAtoms[j]->addFrc(force);
       }
     }
@@ -944,11 +926,11 @@ template<typename T> void ZConstraint<T>::doHarmonic()
       zconsMols[i]->getCOM(COM);
 //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] 
 // 	   << "\tcurrent zpos: " << COM[whichDirection] << endl;
+
+      diff = COM[whichDirection] -zPos[i];
     
-    diff = COM[whichDirection] -zPos[i];
-    
       harmonicU = 0.5 * kz[i] * diff * diff;  
-    info->lrPot += harmonicU;
+      info->lrPot += harmonicU;
 
       harmonicF =  - kz[i] * diff;
       totalFZ_local += harmonicF;
@@ -958,8 +940,8 @@ template<typename T> void ZConstraint<T>::doHarmonic()
       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);
+        //force[whichDirection] = harmonicF / zconsMols[i]->getNAtoms();
+        force[whichDirection] = forcePolicy->getHFOfFixedZMols(zconsMols[i], movingZAtoms[j], harmonicF);
          movingZAtoms[j]->addFrc(force);
        }
     }
@@ -982,8 +964,8 @@ template<typename T> void ZConstraint<T>::doHarmonic()
      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);
+       //force[whichDirection] = - totalFZ /totNumOfUnconsAtoms;
+       force[whichDirection] = - forcePolicy->getHFOfUnconsMols(unconsAtoms[j], totalFZ);
        unconsAtoms[j]->addFrc(force);     
      }
   }   
@@ -1023,8 +1005,9 @@ template<typename T> bool ZConstraint<T>::checkZConsSt
 #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(){
@@ -1046,7 +1029,7 @@ template<typename T> bool ZConstraint<T>::haveFixedZMo
   MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
 #endif
 
-  return havingFixed > 0 ? true : false;
+  return (havingFixed > 0);
 }
 
 
@@ -1072,7 +1055,7 @@ template<typename T> bool ZConstraint<T>::haveMovingZM
   MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT,MPI_SUM, MPI_COMM_WORLD);
 #endif
 
-  return havingMoving > 0 ? true : false;
+  return (havingMoving > 0);
   
 }