--- trunk/OOPSE/libmdtools/ForceFields.cpp	2004/05/01 18:52:38	1144
+++ trunk/OOPSE/libmdtools/ForceFields.cpp	2004/05/11 20:33:41	1157
@@ -73,7 +73,7 @@ void ForceFields::doForces( int calcPot, int calcStres
 
 void ForceFields::doForces( int calcPot, int calcStress ){
 
-  int i, isError;
+  int i, j, isError;
   double* frc;
   double* pos;
   double* trq;
@@ -83,6 +83,21 @@ void ForceFields::doForces( int calcPot, int calcStres
   double* massRatio;
   SimState* config;
 
+  Molecule* myMols;
+  Atom** myAtoms;
+  int numAtom;
+  int curIndex;
+  double mtot;
+  int numMol;
+  int numCutoffGroups;
+  CutoffGroup* myCutoffGroup;
+  vector<CutoffGroup*>::iterator iterCutoff;
+  Atom* cutoffAtom;
+  vector<Atom*>::iterator iterAtom;  
+  double com[3];
+  double tempPos[3];
+  int atomIndex;
+  
   short int passedCalcPot = (short int)calcPot;
   short int passedCalcStress = (short int)calcStress;
 
@@ -100,7 +115,7 @@ void ForceFields::doForces( int calcPot, int calcStres
   for(i=0; i<entry_plug->n_mol; i++ ){
     // CalcForces in molecules takes care of mapping rigid body coordinates
     // into atomic coordinates
-    entry_plug->molecules[i].calcForces();
+    entry_plug->molecules[i].calcForces();    
   }
 
 #ifdef PROFILE
@@ -114,9 +129,56 @@ void ForceFields::doForces( int calcPot, int calcStres
   trq = config->getTrqArray();
   A   = config->getAmatArray();
   u_l = config->getUlArray();
-  rc = config->getRcArray();
-  massRatio = config->getMassRatioArray();
 
+  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);
+
+        for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); cutoffAtom != NULL; 
+                                             cutoffAtom = myCutoffGroup->beginAtom(iterAtom)){
+          cutoffAtom->setRc(com);
+        }  
+                                
+      }// end for(myCutoffGroup)
+      
+    }//end for(int i = 0)
+
+    rc = config->getRcArray();
+  }
+  else{
+    // center of mass of the group is the same as position of the atom  if cutoff group does not exist
+    rc = pos;
+  }
+  
+
+
   isError = 0;
   entry_plug->lrPot = 0.0;