--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/03/16 19:22:56	1091
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/04/28 22:06:29	1139
@@ -9,6 +9,7 @@
 #include "parse_me.h"
 #include "Integrator.hpp"
 #include "simError.h"
+#include "RigidBody.hpp"
 //#include "ConjugateMinimizer.hpp"
 #include "OOPSEMinimizer.hpp"
 
@@ -146,13 +147,6 @@ void SimSetup::createSim(void){
   // make the output filenames
 
   makeOutNames();
-
-  if (globals->haveMinimizer())
-    // make minimizer
-    makeMinimizer();
-  else
-    // make the integrator
-    makeIntegrator();
   
 #ifdef IS_MPI
   mpiSim->mpiRefresh();
@@ -161,65 +155,94 @@ void SimSetup::createSim(void){
   // initialize the Fortran
 
   initFortran();
+
+  if (globals->haveMinimizer())
+    // make minimizer
+    makeMinimizer();
+  else
+    // make the integrator
+    makeIntegrator();
+
 }
 
 
 void SimSetup::makeMolecules(void){
-  int k;
-  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
+  int i, j, k;
+  int exI, exJ, exK, exL, slI, slJ;
+  int tempI, tempJ, tempK, tempL;
+  int molI;
+  int stampID, atomOffset, rbOffset;
   molInit molInfo;
   DirectionalAtom* dAtom;
+  RigidBody* myRB;
+  StuntDouble* mySD;
   LinkedAssign* extras;
   LinkedAssign* current_extra;
   AtomStamp* currentAtom;
   BondStamp* currentBond;
   BendStamp* currentBend;
   TorsionStamp* currentTorsion;
+  RigidBodyStamp* currentRigidBody;
 
   bond_pair* theBonds;
   bend_set* theBends;
   torsion_set* theTorsions;
 
+  set<int> skipList;
+
+  double phi, theta, psi;
+  char* molName;
+  char rbName[100];
+
   //init the forceField paramters
 
   the_ff->readParams();
 
-
   // init the atoms
 
-  double phi, theta, psi;
-  double sux, suy, suz;
-  double Axx, Axy, Axz, Ayx, Ayy, Ayz, Azx, Azy, Azz;
-  double ux, uy, uz, u, uSqr;
+  int nMembers, nNew, rb1, rb2;
 
   for (k = 0; k < nInfo; k++){
     the_ff->setSimInfo(&(info[k]));
 
     atomOffset = 0;
-    excludeOffset = 0;
+
     for (i = 0; i < info[k].n_mol; i++){
       stampID = info[k].molecules[i].getStampID();
+      molName = comp_stamps[stampID]->getID();
 
       molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
       molInfo.nBonds = comp_stamps[stampID]->getNBonds();
       molInfo.nBends = comp_stamps[stampID]->getNBends();
       molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
-      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
-
+      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
+      
       molInfo.myAtoms = &(info[k].atoms[atomOffset]);
-      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
-      molInfo.myBonds = new Bond * [molInfo.nBonds];
-      molInfo.myBends = new Bend * [molInfo.nBends];
-      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
 
+      if (molInfo.nBonds > 0) 
+        molInfo.myBonds = new (Bond *) [molInfo.nBonds];
+      else 
+        molInfo.myBonds = NULL;
+
+      if (molInfo.nBends > 0) 
+        molInfo.myBends = new (Bend *) [molInfo.nBends];
+      else 
+        molInfo.myBends = NULL;
+
+      if (molInfo.nTorsions > 0) 
+        molInfo.myTorsions = new (Torsion *) [molInfo.nTorsions];
+      else 
+        molInfo.myTorsions = NULL;
+
       theBonds = new bond_pair[molInfo.nBonds];
       theBends = new bend_set[molInfo.nBends];
       theTorsions = new torsion_set[molInfo.nTorsions];
-
+      
       // make the Atoms
 
       for (j = 0; j < molInfo.nAtoms; j++){
         currentAtom = comp_stamps[stampID]->getAtom(j);
+
         if (currentAtom->haveOrientation()){
           dAtom = new DirectionalAtom((j + atomOffset),
                                       info[k].getConfiguration());
@@ -233,47 +256,20 @@ void SimSetup::makeMolecules(void){
           phi = currentAtom->getEulerPhi() * M_PI / 180.0;
           theta = currentAtom->getEulerTheta() * M_PI / 180.0;
           psi = currentAtom->getEulerPsi()* M_PI / 180.0;
+
+          dAtom->setUnitFrameFromEuler(phi, theta, psi);
             
-          Axx = (cos(phi) * cos(psi)) - (sin(phi) * cos(theta) * sin(psi));
-          Axy = (sin(phi) * cos(psi)) + (cos(phi) * cos(theta) * sin(psi));
-          Axz = sin(theta) * sin(psi);
-          
-          Ayx = -(cos(phi) * sin(psi)) - (sin(phi) * cos(theta) * cos(psi));
-          Ayy = -(sin(phi) * sin(psi)) + (cos(phi) * cos(theta) * cos(psi));
-          Ayz = sin(theta) * cos(psi);
-          
-          Azx = sin(phi) * sin(theta);
-          Azy = -cos(phi) * sin(theta);
-          Azz = cos(theta);
+        }
+        else{
 
-          sux = 0.0;
-          suy = 0.0;
-          suz = 1.0;
+          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
 
-          ux = (Axx * sux) + (Ayx * suy) + (Azx * suz);
-          uy = (Axy * sux) + (Ayy * suy) + (Azy * suz);
-          uz = (Axz * sux) + (Ayz * suy) + (Azz * suz);
-
-          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
-
-          u = sqrt(uSqr);
-          ux = ux / u;
-          uy = uy / u;
-          uz = uz / u;
-
-          dAtom->setSUx(ux);
-          dAtom->setSUy(uy);
-          dAtom->setSUz(uz);
         }
-        else{
-          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
-                                               info[k].getConfiguration());
-        }
-        molInfo.myAtoms[j]->setType(currentAtom->getType());
 
+        molInfo.myAtoms[j]->setType(currentAtom->getType());
 #ifdef IS_MPI
 
-        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
+        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
 
 #endif // is_mpi
       } 
@@ -284,28 +280,19 @@ void SimSetup::makeMolecules(void){
         theBonds[j].a = currentBond->getA() + atomOffset;
         theBonds[j].b = currentBond->getB() + atomOffset;
 
-        exI = theBonds[j].a;
-        exJ = theBonds[j].b;
+        tempI = theBonds[j].a;
+        tempJ = theBonds[j].b;
 
-        // exclude_I must always be the smaller of the pair
-        if (exI > exJ){
-          tempEx = exI;
-          exI = exJ;
-          exJ = tempEx;
-        }
-#ifdef IS_MPI
-        tempEx = exI;
-        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-        tempEx = exJ;
-        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+#ifdef IS_MPI
+        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
+        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
+#else
+        exI = tempI + 1;
+        exJ = tempJ + 1;
+#endif
 
-        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
-#else  // isn't MPI
-
-        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
-#endif  //is_mpi
+        info[k].excludes->addPair(exI, exJ);
       }
-      excludeOffset += molInfo.nBonds;
 
       //make the bends
       for (j = 0; j < molInfo.nBends; j++){
@@ -355,33 +342,41 @@ void SimSetup::makeMolecules(void){
           }
         }
 
-        if (!theBends[j].isGhost){
-          exI = theBends[j].a;
-          exJ = theBends[j].c;
-        }
-        else{
-          exI = theBends[j].a;
-          exJ = theBends[j].b;
-        }
-
-        // exclude_I must always be the smaller of the pair
-        if (exI > exJ){
-          tempEx = exI;
-          exI = exJ;
-          exJ = tempEx;
-        }
+        if (theBends[j].isGhost) {
+          
+          tempI = theBends[j].a;
+          tempJ = theBends[j].b;
+          
 #ifdef IS_MPI
-        tempEx = exI;
-        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-        tempEx = exJ;
-        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
+          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
+          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
+#else
+          exI = tempI + 1;
+          exJ = tempJ + 1;
+#endif          
+          info[k].excludes->addPair(exI, exJ);
 
-        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
-#else  // isn't MPI
-        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
-#endif  //is_mpi
+        } else {
+
+          tempI = theBends[j].a;
+          tempJ = theBends[j].b;
+          tempK = theBends[j].c;
+          
+#ifdef IS_MPI
+          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
+          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
+          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
+#else
+          exI = tempI + 1;
+          exJ = tempJ + 1;
+          exK = tempK + 1;
+#endif
+          
+          info[k].excludes->addPair(exI, exK);
+          info[k].excludes->addPair(exI, exJ);
+          info[k].excludes->addPair(exJ, exK);
+        }
       }
-      excludeOffset += molInfo.nBends;
 
       for (j = 0; j < molInfo.nTorsions; j++){
         currentTorsion = comp_stamps[stampID]->getTorsion(j);
@@ -390,38 +385,142 @@ void SimSetup::makeMolecules(void){
         theTorsions[j].c = currentTorsion->getC() + atomOffset;
         theTorsions[j].d = currentTorsion->getD() + atomOffset;
 
-        exI = theTorsions[j].a;
-        exJ = theTorsions[j].d;
+        tempI = theTorsions[j].a;       
+        tempJ = theTorsions[j].b;
+        tempK = theTorsions[j].c;
+        tempL = theTorsions[j].d;
 
-        // exclude_I must always be the smaller of the pair
-        if (exI > exJ){
-          tempEx = exI;
-          exI = exJ;
-          exJ = tempEx;
-        }
 #ifdef IS_MPI
-        tempEx = exI;
-        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-        tempEx = exJ;
-        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
-
-        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
-#else  // isn't MPI
-        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
-#endif  //is_mpi
-      }
-      excludeOffset += molInfo.nTorsions;
+        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
+        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
+        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
+        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
+#else
+        exI = tempI + 1;
+        exJ = tempJ + 1;
+        exK = tempK + 1;
+        exL = tempL + 1;
+#endif
 
+        info[k].excludes->addPair(exI, exJ);
+        info[k].excludes->addPair(exI, exK);
+        info[k].excludes->addPair(exI, exL);        
+        info[k].excludes->addPair(exJ, exK);
+        info[k].excludes->addPair(exJ, exL);
+        info[k].excludes->addPair(exK, exL);
+      }
 
-      // send the arrays off to the forceField for init.
+      
+      molInfo.myRigidBodies.clear();
+      
+      for (j = 0; j < molInfo.nRigidBodies; j++){
 
+        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
+        nMembers = currentRigidBody->getNMembers(); 
+
+        // Create the Rigid Body:
+
+        myRB = new RigidBody();
+
+        sprintf(rbName,"%s_RB_%d", molName, j);
+        myRB->setType(rbName);
+        
+        for (rb1 = 0; rb1 < nMembers; rb1++) {
+
+          // molI is atom numbering inside this molecule
+          molI = currentRigidBody->getMember(rb1);     
+
+          // tempI is atom numbering on local processor
+          tempI = molI + atomOffset;
+
+          // currentAtom is the AtomStamp (which we need for 
+          // rigid body reference positions)
+          currentAtom = comp_stamps[stampID]->getAtom(molI);
+
+          // When we add to the rigid body, add the atom itself and 
+          // the stamp info:
+
+          myRB->addAtom(info[k].atoms[tempI], currentAtom);
+          
+          // Add this atom to the Skip List for the integrators
+#ifdef IS_MPI
+          slI = info[k].atoms[tempI]->getGlobalIndex();
+#else
+          slI = tempI;
+#endif
+          skipList.insert(slI);
+          
+        }
+        
+        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
+          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
+            
+            tempI = currentRigidBody->getMember(rb1);
+            tempJ = currentRigidBody->getMember(rb2);
+            
+            // Some explanation is required here.
+            // Fortran indexing starts at 1, while c indexing starts at 0
+            // Also, in parallel computations, the GlobalIndex is
+            // used for the exclude list:
+            
+#ifdef IS_MPI
+            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
+            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
+#else
+            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
+            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
+#endif
+            
+            info[k].excludes->addPair(exI, exJ);
+            
+          }
+        }
+
+        molInfo.myRigidBodies.push_back(myRB);
+        info[k].rigidBodies.push_back(myRB);
+      }
+      
+
+      // After this is all set up, scan through the atoms to 
+      // see if they can be added to the integrableObjects:
+
+      molInfo.myIntegrableObjects.clear();
+      
+
+      for (j = 0; j < molInfo.nAtoms; j++){
+
+#ifdef IS_MPI
+        slJ = molInfo.myAtoms[j]->getGlobalIndex();
+#else
+        slJ = j+atomOffset;
+#endif
+
+        // if they aren't on the skip list, then they can be integrated
+
+        if (skipList.find(slJ) == skipList.end()) {
+          mySD = (StuntDouble *) molInfo.myAtoms[j];
+          info[k].integrableObjects.push_back(mySD);
+          molInfo.myIntegrableObjects.push_back(mySD);
+        }
+      }
+
+      // all rigid bodies are integrated:
+
+      for (j = 0; j < molInfo.nRigidBodies; j++) {
+        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
+        info[k].integrableObjects.push_back(mySD);      
+        molInfo.myIntegrableObjects.push_back(mySD);
+      }
+    
+      
+      // send the arrays off to the forceField for init.
+      
       the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
       the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
       the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
       the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
                                  theTorsions);
 
-
       info[k].molecules[i].initialize(molInfo);
 
 
@@ -429,7 +528,7 @@ void SimSetup::makeMolecules(void){
       delete[] theBonds;
       delete[] theBends;
       delete[] theTorsions;
-    }
+    }    
   }
 
 #ifdef IS_MPI
@@ -439,7 +538,15 @@ void SimSetup::makeMolecules(void){
 
   // clean up the forcefield
 
-  the_ff->calcRcut();
+  if (!globals->haveLJrcut()){
+
+    the_ff->calcRcut();
+
+  } else {
+    
+    the_ff->setRcut( globals->getLJrcut() );
+  }
+
   the_ff->cleanMe();
 }
 
@@ -727,7 +834,7 @@ void SimSetup::gatherInfo(void){
   }
 
   //check whether sample time, status time, thermal time and reset time are divisble by dt
-  if (!isDivisible(globals->getSampleTime(), globals->getDt())){
+  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
     sprintf(painCave.errMsg,
             "Sample time is not divisible by dt.\n"
             "\tThis will result in samples that are not uniformly\n"
@@ -737,7 +844,7 @@ void SimSetup::gatherInfo(void){
     simError();    
   }
 
-  if (globals->haveStatusTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
+  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
     sprintf(painCave.errMsg,
             "Status time is not divisible by dt.\n"
             "\tThis will result in status reports that are not uniformly\n"
@@ -773,12 +880,10 @@ void SimSetup::gatherInfo(void){
     if (globals->haveSampleTime()){
       info[i].sampleTime = globals->getSampleTime();
       info[i].statusTime = info[i].sampleTime;
-      info[i].thermalTime = info[i].sampleTime;
     }
     else{
       info[i].sampleTime = globals->getRunTime();
       info[i].statusTime = info[i].sampleTime;
-      info[i].thermalTime = info[i].sampleTime;
     }
 
     if (globals->haveStatusTime()){
@@ -787,6 +892,8 @@ void SimSetup::gatherInfo(void){
 
     if (globals->haveThermalTime()){
       info[i].thermalTime = globals->getThermalTime();
+    } else {
+      info[i].thermalTime = globals->getRunTime();
     }
 
     info[i].resetIntegrator = 0;
@@ -804,6 +911,8 @@ void SimSetup::gatherInfo(void){
 
     info[i].useInitXSstate = globals->getUseInitXSstate();
     info[i].orthoTolerance = globals->getOrthoBoxTolerance();
+    info[i].useMolecularCutoffs = globals->getUseMolecularCutoffs();
+
     
   }
   
@@ -857,6 +966,7 @@ void SimSetup::finalInfoCheck(void){
 void SimSetup::finalInfoCheck(void){
   int index;
   int usesDipoles;
+  int usesCharges;
   int i;
 
   for (i = 0; i < nInfo; i++){
@@ -868,7 +978,12 @@ void SimSetup::finalInfoCheck(void){
       usesDipoles = (info[i].atoms[index])->hasDipole();
       index++;
     }
-
+    index = 0;
+    usesCharges = 0;
+    while ((index < info[i].n_atoms) && !usesCharges){
+      usesCharges= (info[i].atoms[index])->hasCharge();
+      index++;
+    }
 #ifdef IS_MPI
     int myUse = usesDipoles;
     MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
@@ -919,7 +1034,7 @@ void SimSetup::finalInfoCheck(void){
       info[i].dielectric = globals->getDielectric();
     }
     else{
-      if (usesDipoles){
+      if (usesDipoles || usesCharges){
         if (!globals->haveECR()){
           sprintf(painCave.errMsg,
                   "SimSetup Warning: No value was set for electrostaticCutoffRadius.\n"
@@ -1221,13 +1336,15 @@ void SimSetup::calcSysValues(void){
   tot_bonds = 0;
   tot_bends = 0;
   tot_torsions = 0;
+  tot_rigid = 0;
   for (i = 0; i < n_components; i++){
     tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
     tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
     tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
     tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
+    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
   }
-
+  
   tot_SRI = tot_bonds + tot_bends + tot_torsions;
   molMembershipArray = new int[tot_atoms];
 
@@ -1249,10 +1366,14 @@ void SimSetup::mpiMolDivide(void){
   int i, j, k;
   int localMol, allMol;
   int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
+  int local_rigid;
+  vector<int> globalMolIndex;
 
   mpiSim = new mpiSimulation(info);
 
-  globalIndex = mpiSim->divideLabor();
+  mpiSim->divideLabor();
+  globalAtomIndex = mpiSim->getGlobalAtomIndex();
+  //globalMolIndex = mpiSim->getGlobalMolIndex();
 
   // set up the local variables 
 
@@ -1265,9 +1386,9 @@ void SimSetup::mpiMolDivide(void){
   local_bonds = 0;
   local_bends = 0;
   local_torsions = 0;
-  globalAtomIndex = 0;
+  local_rigid = 0;
+  globalAtomCounter = 0;
 
-
   for (i = 0; i < n_components; i++){
     for (j = 0; j < components_nmol[i]; j++){
       if (mol2proc[allMol] == worldRank){
@@ -1275,11 +1396,12 @@ void SimSetup::mpiMolDivide(void){
         local_bonds += comp_stamps[i]->getNBonds();
         local_bends += comp_stamps[i]->getNBends();
         local_torsions += comp_stamps[i]->getNTorsions();
+        local_rigid += comp_stamps[i]->getNRigidBodies();
         localMol++;
       }      
       for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
-        info[0].molMembershipArray[globalAtomIndex] = allMol;
-        globalAtomIndex++;
+        info[0].molMembershipArray[globalAtomCounter] = allMol;
+        globalAtomCounter++;
       }
 
       allMol++;
@@ -1288,6 +1410,7 @@ void SimSetup::mpiMolDivide(void){
   local_SRI = local_bonds + local_bends + local_torsions;
 
   info[0].n_atoms = mpiSim->getMyNlocal();  
+  
 
   if (local_atoms != info[0].n_atoms){
     sprintf(painCave.errMsg,
@@ -1320,9 +1443,7 @@ void SimSetup::makeSysArrays(void){
 
   Atom** the_atoms;
   Molecule* the_molecules;
-  Exclude** the_excludes;
 
-
   for (l = 0; l < nInfo; l++){
     // create the atom and short range interaction arrays
 
@@ -1348,15 +1469,15 @@ void SimSetup::makeSysArrays(void){
 #else // is_mpi
 
     molIndex = 0;
-    globalAtomIndex = 0;
+    globalAtomCounter = 0;
     for (i = 0; i < n_components; i++){
       for (j = 0; j < components_nmol[i]; j++){
         the_molecules[molIndex].setStampID(i);
         the_molecules[molIndex].setMyIndex(molIndex);
         the_molecules[molIndex].setGlobalIndex(molIndex);
         for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
-          info[l].molMembershipArray[globalAtomIndex] = molIndex;
-          globalAtomIndex++;
+          info[l].molMembershipArray[globalAtomCounter] = molIndex;
+          globalAtomCounter++;
         }
         molIndex++;
       }
@@ -1365,32 +1486,14 @@ void SimSetup::makeSysArrays(void){
 
 #endif // is_mpi
 
-
-    if (info[l].n_SRI){
-      Exclude::createArray(info[l].n_SRI);
-      the_excludes = new Exclude * [info[l].n_SRI];
-      for (int ex = 0; ex < info[l].n_SRI; ex++){
-        the_excludes[ex] = new Exclude(ex);
-      }
-      info[l].globalExcludes = new int;
-      info[l].n_exclude = info[l].n_SRI;
-    }
-    else{
-      Exclude::createArray(1);
-      the_excludes = new Exclude * ;
-      the_excludes[0] = new Exclude(0);
-      the_excludes[0]->setPair(0, 0);
-      info[l].globalExcludes = new int;
-      info[l].globalExcludes[0] = 0;
-      info[l].n_exclude = 0;
-    }
-
+    info[l].globalExcludes = new int;
+    info[l].globalExcludes[0] = 0;
+    
     // set the arrays into the SimInfo object
 
     info[l].atoms = the_atoms;
     info[l].molecules = the_molecules;
     info[l].nGlobalExcludes = 0;
-    info[l].excludes = the_excludes;
 
     the_ff->setSimInfo(info);
   }
@@ -1681,7 +1784,15 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){
     zconsFixtime->setData(globals->getZconsFixtime());
     theInfo.addProperty(zconsFixtime);  
   }
+
+  //set zconsUsingSMD
+  IntData* zconsUsingSMD = new IntData();
+  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
 
+  if (globals->haveZConsUsingSMD()){
+    zconsUsingSMD->setData(globals->getZconsUsingSMD());
+    theInfo.addProperty(zconsUsingSMD);  
+  }
 
   //Determine the name of ouput file and add it into SimInfo's property list 
   //Be careful, do not use inFileName, since it is a pointer which
@@ -1712,7 +1823,8 @@ void SimSetup::setupZConstraint(SimInfo& theInfo){
     tempParaItem.zPos = zconStamp[i]->getZpos();
     tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
     tempParaItem.kRatio = zconStamp[i]->getKratio();
-
+    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
+    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
     zconsParaData->addItem(tempParaItem);
   }