--- trunk/OOPSE/libmdtools/ForceFields.cpp	2004/05/11 20:33:41	1157
+++ trunk/OOPSE/libmdtools/ForceFields.cpp	2004/06/01 18:07:34	1213
@@ -81,6 +81,7 @@ void ForceFields::doForces( int calcPot, int calcStres
   double* u_l;
   double* rc;
   double* massRatio;
+  double factor;
   SimState* config;
 
   Molecule* myMols;
@@ -92,11 +93,8 @@ void ForceFields::doForces( int calcPot, int calcStres
   int numCutoffGroups;
   CutoffGroup* myCutoffGroup;
   vector<CutoffGroup*>::iterator iterCutoff;
-  Atom* cutoffAtom;
-  vector<Atom*>::iterator iterAtom;  
   double com[3];
-  double tempPos[3];
-  int atomIndex;
+  vector<double> rcGroup;
   
   short int passedCalcPot = (short int)calcPot;
   short int passedCalcStress = (short int)calcStress;
@@ -131,46 +129,25 @@ void ForceFields::doForces( int calcPot, int calcStres
   u_l = config->getUlArray();
 
   if(entry_plug->haveCutoffGroups){
-    //if 
     myMols = entry_plug->molecules;
     numMol = entry_plug->n_mol;
     for(int i  = 0; i < numMol; i++){
-      numAtom = myMols[i].getNAtoms();
-      myAtoms = myMols[i].getMyAtoms();
-
-      
-      for(int j = 0; j < numAtom; j++){
-#ifdef IS_MPI
-        atomIndex = myAtoms[j]->getGlobalIndex();
-#else
-        atomIndex = myAtoms[j]->getIndex();
-#endif
-
-        if(myMols[i].belongToCutoffGroup(atomIndex))
-          continue;
-        else{
-          myAtoms[j]->getPos(tempPos);
-          myAtoms[j]->setRc(tempPos);
-        }
-          
-      }
         
       numCutoffGroups = myMols[i].getNCutoffGroups();
       for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); myCutoffGroup != NULL; 
                                                     myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
         //get center of mass of the cutoff group
-        myCutoffGroup->getCOM(com);
+        myCutoffGroup->getCOM(com); 
 
-        for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL; 
-                                             cutoffAtom = myCutoffGroup->beginAtom(iterAtom)){
-          cutoffAtom->setRc(com);
-        }  
-                                
+        rcGroup.push_back(com[0]);
+        rcGroup.push_back(com[1]);
+        rcGroup.push_back(com[2]);
+        
       }// end for(myCutoffGroup)
       
     }//end for(int i = 0)
 
-    rc = config->getRcArray();
+    rc = &rcGroup[0];
   }
   else{
     // center of mass of the group is the same as position of the atom  if cutoff group does not exist
@@ -218,9 +195,41 @@ void ForceFields::doForces( int calcPot, int calcStres
 
   for(i=0; i<entry_plug->n_mol; i++ ){
     entry_plug->molecules[i].atoms2rigidBodies();
+  }
+
+
+  if (entry_plug->useSolidThermInt && !entry_plug->useLiquidThermInt) {
+
+    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
+    for (i=0; i < entry_plug->integrableObjects.size(); i++) {
+      for (j=0; j< 3; j++) 
+	frc[3*i + j] *= factor; 
+      if (entry_plug->integrableObjects[i]->isDirectional()) {
+	for (j=0; j< 3; j++) 
+	  trq[3*i + j] *= factor;
+      }
+    }
+    entry_plug->vRaw = entry_plug->lrPot;
+    entry_plug->lrPot *= factor;
+    entry_plug->lrPot += entry_plug->restraint->Calc_Restraint_Forces(entry_plug->integrableObjects);
+    entry_plug->vHarm = entry_plug->restraint->getVharm();
   }
 
+  if (entry_plug->useLiquidThermInt) {
 
+    factor = pow(entry_plug->thermIntLambda, entry_plug->thermIntK);
+    for (i=0; i < entry_plug->integrableObjects.size(); i++) {
+      for (j=0; j< 3; j++) 
+	frc[3*i + j] *= factor; 
+      if (entry_plug->integrableObjects[i]->isDirectional()) {
+	for (j=0; j< 3; j++) 
+	  trq[3*i + j] *= factor;
+      }
+    }
+    entry_plug->vRaw = entry_plug->lrPot;
+    entry_plug->lrPot *= factor;
+  }
+
 #ifdef IS_MPI
   sprintf( checkPointMsg,
 	   "returned from the force calculation.\n" );
@@ -252,3 +261,22 @@ void ForceFields::initFortran(int ljMixPolicy, int use
 #endif // is_mpi
   
 }
+
+
+void ForceFields::initRestraints(){
+  int i;
+  // store the initial info.
+  // set the omega values to zero
+  for (i=0; i<entry_plug->integrableObjects.size(); i++)
+    entry_plug->integrableObjects[i]->setZangle( 0.0 );
+
+  entry_plug->restraint->Store_Init_Info(entry_plug->integrableObjects);
+
+}
+
+void ForceFields::dumpzAngle(){
+
+  // store the initial info.
+  entry_plug->restraint->Write_zAngle_File(entry_plug->integrableObjects);
+
+}