--- trunk/src/brains/ForceManager.cpp	2005/10/04 19:34:03	645
+++ trunk/src/brains/ForceManager.cpp	2006/05/17 21:51:42	963
@@ -50,9 +50,31 @@
 #include "brains/ForceManager.hpp"
 #include "primitives/Molecule.hpp"
 #include "UseTheForce/doForces_interface.h"
+#define __C
+#include "UseTheForce/DarkSide/fInteractionMap.h"
 #include "utils/simError.h"
+#include "primitives/Bend.hpp"
+#include "primitives/Bend.hpp"
 namespace oopse {
 
+/*
+  struct BendOrderStruct {
+    Bend* bend;
+    BendDataSet dataSet;
+  };
+  struct TorsionOrderStruct {
+    Torsion* torsion;
+    TorsionDataSet dataSet;
+  };
+
+  bool  BendSortFunctor(const BendOrderStruct& b1, const BendOrderStruct& b2) {
+    return b1.dataSet.deltaV < b2.dataSet.deltaV;
+  }
+
+  bool  TorsionSortFunctor(const TorsionOrderStruct& t1, const TorsionOrderStruct& t2) {
+    return t1.dataSet.deltaV < t2.dataSet.deltaV;
+  }
+  */
   void ForceManager::calcForces(bool needPotential, bool needStress) {
 
     if (!info_->isFortranInitialized()) {
@@ -66,7 +88,37 @@ namespace oopse {
     calcLongRangeInteraction(needPotential, needStress);
 
     postCalculation();
-        
+
+/*
+    std::vector<BendOrderStruct> bendOrderStruct;
+    for(std::map<Bend*, BendDataSet>::iterator i = bendDataSets.begin(); i != bendDataSets.end(); ++i) {
+        BendOrderStruct tmp;
+        tmp.bend= const_cast<Bend*>(i->first);
+        tmp.dataSet = i->second;
+        bendOrderStruct.push_back(tmp);
+    }
+
+    std::vector<TorsionOrderStruct> torsionOrderStruct;
+    for(std::map<Torsion*, TorsionDataSet>::iterator j = torsionDataSets.begin(); j != torsionDataSets.end(); ++j) {
+        TorsionOrderStruct tmp;
+        tmp.torsion = const_cast<Torsion*>(j->first);
+        tmp.dataSet = j->second;
+        torsionOrderStruct.push_back(tmp);
+    }
+    
+    std::sort(bendOrderStruct.begin(), bendOrderStruct.end(), std::ptr_fun(BendSortFunctor));
+    std::sort(torsionOrderStruct.begin(), torsionOrderStruct.end(), std::ptr_fun(TorsionSortFunctor));
+    for (std::vector<BendOrderStruct>::iterator k = bendOrderStruct.begin(); k != bendOrderStruct.end(); ++k) {
+        Bend* bend = k->bend;
+        std::cout << "Bend: atom1=" <<bend->getAtomA()->getGlobalIndex() << ",atom2 = "<< bend->getAtomB()->getGlobalIndex() << ",atom3="<<bend->getAtomC()->getGlobalIndex() << " ";
+        std::cout << "deltaV=" << k->dataSet.deltaV << ",p_theta=" << k->dataSet.prev.angle <<",p_pot=" << k->dataSet.prev.potential<< ",c_theta=" << k->dataSet.curr.angle << ", c_pot = " << k->dataSet.curr.potential <<std::endl;
+    }
+    for (std::vector<TorsionOrderStruct>::iterator l = torsionOrderStruct.begin(); l != torsionOrderStruct.end(); ++l) {
+        Torsion* torsion = l->torsion;
+        std::cout << "Torsion: atom1=" <<torsion->getAtomA()->getGlobalIndex() << ",atom2 = "<< torsion->getAtomB()->getGlobalIndex() << ",atom3="<<torsion->getAtomC()->getGlobalIndex() << ",atom4="<<torsion->getAtomD()->getGlobalIndex()<< " ";
+        std::cout << "deltaV=" << l->dataSet.deltaV << ",p_theta=" << l->dataSet.prev.angle <<",p_pot=" << l->dataSet.prev.potential<< ",c_theta=" << l->dataSet.curr.angle << ", c_pot = " << l->dataSet.curr.potential <<std::endl;
+    }
+   */ 
   }
 
   void ForceManager::preCalculation() {
@@ -103,47 +155,87 @@ namespace oopse {
     Molecule::BondIterator bondIter;;
     Molecule::BendIterator  bendIter;
     Molecule::TorsionIterator  torsionIter;
+    RealType bondPotential = 0.0;
+    RealType bendPotential = 0.0;
+    RealType torsionPotential = 0.0;
 
     //calculate short range interactions    
     for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
 
       //change the positions of atoms which belong to the rigidbodies
       for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
-	rb->updateAtoms();
+  	  rb->updateAtoms();
       }
 
       for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
-	bond->calcForce();
+        bond->calcForce();
+        bondPotential += bond->getPotential();
       }
 
+
       for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
-	bend->calcForce();
+
+          RealType angle;
+	    bend->calcForce(angle);
+          RealType currBendPot = bend->getPotential();          
+  	    bendPotential += bend->getPotential();
+          std::map<Bend*, BendDataSet>::iterator i = bendDataSets.find(bend);
+          if (i == bendDataSets.end()) {
+            BendDataSet dataSet;
+            dataSet.prev.angle = dataSet.curr.angle = angle;
+            dataSet.prev.potential = dataSet.curr.potential = currBendPot;
+            dataSet.deltaV = 0.0;
+            bendDataSets.insert(std::map<Bend*, BendDataSet>::value_type(bend, dataSet));
+          }else {
+            i->second.prev.angle = i->second.curr.angle;
+            i->second.prev.potential = i->second.curr.potential;
+            i->second.curr.angle = angle;
+            i->second.curr.potential = currBendPot;
+            i->second.deltaV =  fabs(i->second.curr.potential -  i->second.prev.potential);
+          }
       }
 
       for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
-	torsion->calcForce();
+        RealType angle;
+  	  torsion->calcForce(angle);
+        RealType currTorsionPot = torsion->getPotential();
+	  torsionPotential += torsion->getPotential();
+          std::map<Torsion*, TorsionDataSet>::iterator i = torsionDataSets.find(torsion);
+          if (i == torsionDataSets.end()) {
+            TorsionDataSet dataSet;
+            dataSet.prev.angle = dataSet.curr.angle = angle;
+            dataSet.prev.potential = dataSet.curr.potential = currTorsionPot;
+            dataSet.deltaV = 0.0;
+            torsionDataSets.insert(std::map<Torsion*, TorsionDataSet>::value_type(torsion, dataSet));
+          }else {
+            i->second.prev.angle = i->second.curr.angle;
+            i->second.prev.potential = i->second.curr.potential;
+            i->second.curr.angle = angle;
+            i->second.curr.potential = currTorsionPot;
+            i->second.deltaV =  fabs(i->second.curr.potential -  i->second.prev.potential);
+          }      
       }
 
     }
     
-    double  shortRangePotential = 0.0;
-    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
-      shortRangePotential += mol->getPotential();
-    }
-
+    RealType  shortRangePotential = bondPotential + bendPotential + torsionPotential;    
     Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
     curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential;
+    curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential;
+    curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential;
+    curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential;
+    
   }
 
   void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) {
     Snapshot* curSnapshot;
     DataStorage* config;
-    double* frc;
-    double* pos;
-    double* trq;
-    double* A;
-    double* electroFrame;
-    double* rc;
+    RealType* frc;
+    RealType* pos;
+    RealType* trq;
+    RealType* A;
+    RealType* electroFrame;
+    RealType* rc;
     
     //get current snapshot from SimInfo
     curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
@@ -180,12 +272,18 @@ namespace oopse {
     }
   
     //initialize data before passing to fortran
-    double longRangePotential = 0.0;
+    RealType longRangePotential[LR_POT_TYPES];
+    RealType lrPot = 0.0;
+    
     Mat3x3d tau;
     short int passedCalcPot = needPotential;
     short int passedCalcStress = needStress;
     int isError = 0;
 
+    for (int i=0; i<LR_POT_TYPES;i++){
+      longRangePotential[i]=0.0; //Initialize array
+    }
+
     doForceLoop( pos,
 		 rc,
 		 A,
@@ -193,7 +291,7 @@ namespace oopse {
 		 frc,
 		 trq,
 		 tau.getArrayPointer(),
-		 &longRangePotential, 
+		 longRangePotential, 
 		 &passedCalcPot,
 		 &passedCalcStress,
 		 &isError );
@@ -204,9 +302,15 @@ namespace oopse {
       painCave.isFatal = 1;
       simError();
     }
+    for (int i=0; i<LR_POT_TYPES;i++){
+      lrPot += longRangePotential[i]; //Quick hack
+    }
 
     //store the tau and long range potential    
-    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential;
+    curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot;
+    curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT];
+    curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT];
+
     curSnapshot->statData.setTau(tau);
   }