--- trunk/OOPSE/libmdtools/ZConstraint.cpp	2004/03/16 19:22:56	1091
+++ trunk/OOPSE/libmdtools/ZConstraint.cpp	2004/03/17 14:22:59	1093
@@ -5,12 +5,10 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 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),
+                                                                       usingSMD(false),
                                                                        hasZConsGap(false){
   //get properties from SimInfo
   GenericData* data;
@@ -21,6 +19,7 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
   DoubleData* fixtime;
   StringData* policy;
   StringData* filename; 
+  IntData* smdFlag;
   double COM[3];
 
   //by default, the direction of constraint is z
@@ -150,33 +149,59 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
 
   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();
+    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();
+    }
   }
-  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 zconsFixTime from SimInfo\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    else{
+      this->zconsFixTime = fixtime->getData();
+    }
   }
-
-  if (!fixtime){
-    sprintf(painCave.errMsg,
-            "ZConstraint error: Can not get property from SimInfo\n");
-    painCave.isFatal = 1;
-    simError();
+  else if(hasZConsGap){
+      sprintf(painCave.errMsg,
+              "ZConstraint error: must set fixtime if already set zconsGap\n");
+      painCave.isFatal = 1;
+      simError();
   }
-  else{
-    this->zconsFixTime = fixtime->getData();
+
+
+
+  data = info->getProperty(ZCONSUSINGSMD_ID);
+
+  if (data){
+    smdFlag = dynamic_cast<IntData*>(data);
+
+    if (!smdFlag){
+      sprintf(painCave.errMsg,
+              "ZConstraint error: Can not get property from SimInfo\n");
+      painCave.isFatal = 1;
+      simError();
+    }
+    else{
+      this->usingSMD= smdFlag->getData() ? true : false;
+    }
+
   }
 
 
@@ -292,11 +317,11 @@ 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);
-      molecules[i].getCOM(COM);
+      
+      if(usingSMD)
+        cantVel.push_back((*parameters)[searchResult].cantVel);
+
     }
     else{
       unconsMols.push_back(&molecules[i]);
@@ -304,21 +329,16 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
     }
   }
 
-  fz = new double[zconsMols.size()];
-  curZPos = new double[zconsMols.size()]; 
-  indexOfZConsMols = new int [zconsMols.size()];
+  fz.resize(zconsMols.size());
+  curZPos.resize(zconsMols.size());
+  indexOfZConsMols.resize(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++){
+  for (size_t i = 0; i < zconsMols.size(); i++){
     indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
 
     zconsMols[i]->getCOM(COM);
+    
     if (fabs(zPos[i] - COM[whichDirection]) < zconsTol){
       states.push_back(zcsFixed);
 
@@ -331,15 +351,19 @@ template<typename T> ZConstraint<T>::ZConstraint(SimIn
       if (hasZConsGap)
         endFixTime.push_back(INFINITE_TIME);
     }
+
+    if(usingSMD)
+      cantPos.push_back(COM[whichDirection]);    
   }
 
 #endif 
 
+  
   //get total masss of unconstraint molecules
   double totalMassOfUncons_local;
   totalMassOfUncons_local = 0;
 
-  for (int i = 0; i < (int) (unconsMols.size()); i++)
+  for (size_t i = 0; i < unconsMols.size(); i++)
     totalMassOfUncons_local += unconsMols[i]->getTotalMass();
 
 #ifndef IS_MPI
@@ -366,18 +390,7 @@ template<typename T> ZConstraint<T>::~ZConstraint(){
 }
 
 template<typename T> ZConstraint<T>::~ZConstraint(){
-  if (fz){
-    delete[] fz;
-  }
 
-  if (curZPos){
-    delete[] curZPos;
-  }
-
-  if (indexOfZConsMols){
-    delete[] indexOfZConsMols;
-  }
-
   if (fzOut){
     delete fzOut;
   }
@@ -401,6 +414,8 @@ template<typename T> void ZConstraint<T>::update(){
   massOfZConsMols.clear();
   zPos.clear();
   kz.clear();
+  cantPos.clear();
+  cantVel.clear();
 
   unconsMols.clear();
   massOfUnconsMols.clear();
@@ -414,9 +429,11 @@ 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);
-      massOfZConsMols.push_back(molecules[i].getTotalMass());  
+      massOfZConsMols.push_back(molecules[i].getTotalMass());
+      
+      if(usingSMD)
+        cantVel.push_back((*parameters)[index].cantVel);
 
-      molecules[i].getCOM(COM);
     }
     else{
       unconsMols.push_back(&molecules[i]);
@@ -424,46 +441,19 @@ template<typename T> void ZConstraint<T>::update(){
     }
   }
 
-
-  //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;
+  fz.resize(zconsMols.size());
+  curZPos.resize(zconsMols.size());
+  indexOfZConsMols.resize(zconsMols.size());  
+ 
+  for (size_t i = 0; i < zconsMols.size(); i++){
+    indexOfZConsMols[i] = zconsMols[i]->getGlobalIndex();
   }
-
+    
   //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);
 
@@ -476,6 +466,9 @@ template<typename T> void ZConstraint<T>::update(){
       if (hasZConsGap)
         endFixTime.push_back(INFINITE_TIME);
     }
+
+    if(usingSMD)
+      cantPos.push_back(COM[whichDirection]);        
   } 
   //
   forcePolicy->update();
@@ -580,9 +573,12 @@ template<typename T> void ZConstraint<T>::calcForce(in
     this->doZconstraintForce();
   }
 
-  //use harmonical poteintial to move the molecules to the specified positions
+  //use external force to move the molecules to the specified positions
   if (haveMovingZMols()){
-    this->doHarmonic();
+    if (usingSMD)
+      this->doHarmonic(cantPos);
+    else
+      this->doHarmonic(zPos);      
   }
 
   //write out forces and current positions of z-constraint molecules
@@ -602,8 +598,8 @@ template<typename T> void ZConstraint<T>::calcForce(in
         }
       }
     }
-    fzOut->writeFZ(info->getTime(), zconsMols.size(), indexOfZConsMols, fz,
-                   curZPos, &zPos[0]);
+    fzOut->writeFZ(info->getTime(), zconsMols.size(), &indexOfZConsMols[0], &fz[0],
+                   &curZPos[0], &zPos[0]);
     curZconsTime += zconsTime;
   }
 
@@ -887,7 +883,7 @@ template<typename T> void ZConstraint<T>::doHarmonic()
   *
   */
 
-template<typename T> void ZConstraint<T>::doHarmonic(){
+template<typename T> void ZConstraint<T>::doHarmonic(vector<double>& resPos){
   double force[3];
   double harmonicU;
   double harmonicF;
@@ -908,7 +904,7 @@ template<typename T> void ZConstraint<T>::doHarmonic()
       //       cout << "Moving Molecule\tindex: " << indexOfZConsMols[i] 
       //     << "\tcurrent zpos: " << COM[whichDirection] << endl;
 
-      diff = COM[whichDirection] - zPos[i];
+      diff = COM[whichDirection] - resPos[i];
 
       harmonicU = 0.5 * kz[i] * diff * diff;  
       info->lrPot += harmonicU;
@@ -1249,12 +1245,37 @@ template<typename T> void ZConstraint<T>::updateZPos()
 
 template<typename T> void ZConstraint<T>::updateZPos(){
   double curTime;
-
+  double COM[3];
+  
   curTime = info->getTime();
 
   for (size_t i = 0; i < zconsMols.size(); i++){
+
     if (states[i] == zcsFixed && curTime >= endFixTime[i]){
       zPos[i] += zconsGap;
+
+      if (usingSMD){
+        zconsMols[i]->getCOM(COM);
+        cantPos[i] = COM[whichDirection];
+      }
+      
     }
+    
   }
+  
 }
+
+template<typename T> void ZConstraint<T>::updateCantPos(){
+  double curTime;
+  double dt;
+
+  curTime = info->getTime();
+  dt = info->dt;
+
+  for (size_t i = 0; i < zconsMols.size(); i++){
+    if (states[i] == zcsMoving){
+      cantPos[i] += cantVel[i] * dt;
+    }
+  }
+
+}