--- trunk/OOPSE/libmdtools/SimSetup.cpp	2004/06/03 20:02:25	1229
+++ trunk/OOPSE/libmdtools/SimSetup.cpp	2004/07/16 16:29:44	1330
@@ -11,8 +11,9 @@
 #include "simError.h"
 #include "RigidBody.hpp"
 #include "OOPSEMinimizer.hpp"
-//#include "ConstraintElement.hpp"
-//#include "ConstraintPair.hpp"
+#include "ConstraintElement.hpp"
+#include "ConstraintPair.hpp"
+#include "ConstraintManager.hpp"
 
 #ifdef IS_MPI
 #include "mpiBASS.h"
@@ -157,6 +158,10 @@ void SimSetup::createSim(void){
 
   initFortran();
 
+  //creat constraint manager
+  for(int i = 0; i < nInfo; i++)
+    info[i].consMan = new ConstraintManager(&info[i]);
+
   if (globals->haveMinimizer())
     // make minimizer
     makeMinimizer();
@@ -199,12 +204,13 @@ void SimSetup::makeMolecules(void){
   char* molName;
   char rbName[100];
 
-  //ConstraintPair* consPair; //constraint pair
-  //ConstraintElement* consElement1;  //first element of constraint pair
-  //ConstraintElement* consElement2;  //second element of constraint pair
-  //int whichRigidBody; 
-  //int consAtomIndex;  //index of constraint atom in rigid body's atom array
-  //vector<pair<int, int> > jointAtoms;
+  ConstraintPair* consPair; //constraint pair
+  ConstraintElement* consElement1;  //first element of constraint pair
+  ConstraintElement* consElement2;  //second element of constraint pair
+  int whichRigidBody; 
+  int consAtomIndex;  //index of constraint atom in rigid body's atom array
+  vector<pair<int, int> > jointAtoms;
+  double bondLength2;
   //init the forceField paramters
 
   the_ff->readParams();
@@ -599,37 +605,50 @@ void SimSetup::makeMolecules(void){
         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);
 
-    /*
 
       //creat ConstraintPair.
-      molInfo.myConstraintPair.clear();
+      molInfo.myConstraintPairs.clear();
       
       for (j = 0; j < molInfo.nBonds; j++){
 
-        //if both atoms are in the same rigid body, just skip it
-        currentBond = comp_stamps[stampID]->getBond(j);
-        if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
+        //if bond is constrained bond, add it into constraint pair
+        if(molInfo.myBonds[j]->is_constrained()){
 
-          tempI = currentBond->getA() + atomOffset;
-          if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
-            consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
-          else
-             consElement1 = new ConstraintAtom(info[k].atoms[tempI]);       
+          //if both atoms are in the same rigid body, just skip it
+          currentBond = comp_stamps[stampID]->getBond(j);
+          
+          if(!comp_stamps[stampID]->isBondInSameRigidBody(currentBond)){
 
-          tempJ =  currentBond->getB() + atomOffset;
-          if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
-            consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
-          else
-             consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);     
+            tempI = currentBond->getA() + atomOffset;
+            if( comp_stamps[stampID]->isAtomInRigidBody(currentBond->getA(), whichRigidBody, consAtomIndex))
+              consElement1 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
+            else
+               consElement1 = new ConstraintAtom(info[k].atoms[tempI]);       
 
-          consPair = new DistanceConstraintPair(consElement1, consElement2);
-          molInfo.myConstraintPairs.push_back(consPair);
-        }
+            tempJ =  currentBond->getB() + atomOffset;
+            if(comp_stamps[stampID]->isAtomInRigidBody(currentBond->getB(), whichRigidBody, consAtomIndex))
+              consElement2 = new ConstraintRigidBody(molInfo.myRigidBodies[whichRigidBody], consAtomIndex);
+            else
+               consElement2 = new ConstraintAtom(info[k].atoms[tempJ]);     
+
+            bondLength2 = molInfo.myBonds[j]->get_constraint()->get_dsqr();            
+            consPair = new DistanceConstraintPair(consElement1, consElement2, bondLength2);
+
+            molInfo.myConstraintPairs.push_back(consPair);
+          }
+        }//end if(molInfo.myBonds[j]->is_constrained())
       }  
       
-      //loop over rigid bodies, if two rigid bodies share same joint, creat a HingeConstraintPair
+      //loop over rigid bodies, if two rigid bodies share same joint, creat a JointConstraintPair
       for (int rb1 = 0; rb1 < molInfo.nRigidBodies -1 ; rb1++){
         for (int rb2 = rb1 + 1; rb2 < molInfo.nRigidBodies ; rb2++){
           
@@ -646,14 +665,6 @@ void SimSetup::makeMolecules(void){
         }
       }
       
-*/      
-      // 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);
       
@@ -905,8 +916,12 @@ void SimSetup::gatherInfo(void){
          painCave.isFatal = 1;
          simError();
   }
-
-    // get the ensemble
+  if (globals->haveForceFieldVariant()) {
+    strcpy(forcefield_variant, globals->getForceFieldVariant());
+    has_forcefield_variant = 1;
+  }
+  
+  // get the ensemble
 
   strcpy(ensemble, globals->getEnsemble());
 
@@ -1503,7 +1518,7 @@ void SimSetup::createFF(void){
 void SimSetup::createFF(void){
   switch (ffCase){
     case FF_DUFF:
-      the_ff = new DUFF();
+        the_ff = new DUFF();
       break;
 
     case FF_LJ:
@@ -1511,7 +1526,10 @@ void SimSetup::createFF(void){
       break;
 
     case FF_EAM:
-      the_ff = new EAM_FF();
+      if (has_forcefield_variant) 
+        the_ff = new EAM_FF(forcefield_variant);
+      else
+        the_ff = new EAM_FF();
       break;
 
     case FF_H2O:
@@ -1525,6 +1543,7 @@ void SimSetup::createFF(void){
       simError();
   }
 
+
 #ifdef IS_MPI
   strcpy(checkPointMsg, "ForceField creation successful");
   MPIcheckPoint();
@@ -1801,11 +1820,11 @@ void SimSetup::makeIntegrator(void){
 void SimSetup::makeIntegrator(void){
   int k;
 
-  NVE<RealIntegrator>* myNVE = NULL;
-  NVT<RealIntegrator>* myNVT = NULL;
-  NPTi<NPT<RealIntegrator> >* myNPTi = NULL;
-  NPTf<NPT<RealIntegrator> >* myNPTf = NULL;
-  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL;
+  NVE<Integrator<BaseIntegrator> >* myNVE = NULL;
+  NVT<Integrator<BaseIntegrator> >* myNVT = NULL;
+  NPTi<NPT<Integrator<BaseIntegrator> > >* myNPTi = NULL;
+  NPTf<NPT<Integrator<BaseIntegrator> > >* myNPTf = NULL;
+  NPTxyz<NPT<Integrator<BaseIntegrator> > >* myNPTxyz = NULL;
   
   for (k = 0; k < nInfo; k++){
     switch (ensembleCase){