--- branches/development/src/parallel/ForceMatrixDecomposition.hpp	2011/05/26 13:55:04	1569
+++ branches/development/src/parallel/ForceMatrixDecomposition.hpp	2012/06/18 18:23:20	1756
@@ -36,7 +36,8 @@
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
  * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
 #ifndef PARALLEL_FORCEMATRIXDECOMPOSITION_HPP
@@ -55,19 +56,22 @@ namespace OpenMD {
   
   class ForceMatrixDecomposition : public ForceDecomposition {
   public:
-    ForceMatrixDecomposition(SimInfo* info) : ForceDecomposition(info) {};
+    ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan);
 
     void distributeInitialData();
+    void zeroWorkArrays();
     void distributeData();
     void collectIntermediateData();
     void distributeIntermediateData();
+    void collectSelfData();
     void collectData();
 
     // neighbor list routines
     vector<pair<int, int> >  buildNeighborList();
 
     // group bookkeeping
-    pair<int, int> getGroupTypes(int cg1, int cg2);
+    groupCutoffs getGroupCutoffs(int cg1, int cg2);
+    Vector3d getGroupVelocityColumn(int cg2);
 
     // Group->atom bookkeeping
     vector<int> getAtomsInGroupRow(int cg1);
@@ -82,27 +86,31 @@ namespace OpenMD {
     Vector3d getInteratomicVector(int atom1, int atom2);
        
     // atom bookkeeping
-    vector<int> getAtomList();
-    vector<int> getSkipsForAtom(int atom1);
-    bool skipAtomPair(int atom1, int atom2);
+    int getNAtomsInRow();
+    int getTopologicalDistance(int atom1, int atom2);
+    vector<int> getExcludesForAtom(int atom1); 
+    bool skipAtomPair(int atom1, int atom2, int cg1, int cg2);
+    bool excludeAtomPair(int atom1, int atom2);
     void addForceToAtomRow(int atom1, Vector3d fg);
     void addForceToAtomColumn(int atom2, Vector3d fg);
+    Vector3d getAtomVelocityColumn(int atom2);
 
     // filling interaction blocks with pointers
-    InteractionData fillInteractionData(int atom1, int atom2);
-    InteractionData fillSkipData(int atom1, int atom2);
+    void fillInteractionData(InteractionData &idat, int atom1, int atom2);
+    void unpackInteractionData(InteractionData &idat, int atom1, int atom2);
 
   private: 
+    void createGtypeCutoffMap();
+
     int nLocal_;
     int nGroups_;
     vector<int> AtomLocalToGlobal;
     vector<int> cgLocalToGlobal;
-    vector<RealType> pot_local;
-    vector<RealType> massFactorsLocal;
+    vector<RealType> groupCutoff;
+    vector<int> groupToGtype;
 
 #ifdef IS_MPI    
-    DataStorage atomRowData;
-    DataStorage atomColData;
+    DataStorage atomRowData;    DataStorage atomColData;
     DataStorage cgRowData;
     DataStorage cgColData;
 
@@ -111,33 +119,49 @@ namespace OpenMD {
     int nGroupsInRow_;
     int nGroupsInCol_;
 
-    Communicator<Row, int>* AtomCommIntRow;
-    Communicator<Row, RealType>* AtomCommRealRow; 
-    Communicator<Row, Vector3d>* AtomCommVectorRow; 
-    Communicator<Row, Mat3x3d>*  AtomCommMatrixRow; 
+    Communicator<Row> rowComm;
+    Communicator<Column> colComm;
 
-    Communicator<Column, int>* AtomCommIntColumn;
-    Communicator<Column, RealType>* AtomCommRealColumn; 
-    Communicator<Column, Vector3d>* AtomCommVectorColumn; 
-    Communicator<Column, Mat3x3d>*  AtomCommMatrixColumn; 
+    Plan<int>* AtomPlanIntRow;
+    Plan<RealType>* AtomPlanRealRow; 
+    Plan<Vector3d>* AtomPlanVectorRow; 
+    Plan<Mat3x3d>*  AtomPlanMatrixRow; 
+    Plan<potVec>* AtomPlanPotRow;
 
-    Communicator<Row, int>* cgCommIntRow;
-    Communicator<Row, Vector3d>* cgCommVectorRow; 
-    Communicator<Column, int>* cgCommIntColumn;
-    Communicator<Column, Vector3d>* cgCommVectorColumn; 
+    Plan<int>* AtomPlanIntColumn;
+    Plan<RealType>* AtomPlanRealColumn; 
+    Plan<Vector3d>* AtomPlanVectorColumn; 
+    Plan<Mat3x3d>*  AtomPlanMatrixColumn;
+    Plan<potVec>* AtomPlanPotColumn; 
 
-    vector<vector<RealType> > pot_row;
-    vector<vector<RealType> > pot_col;
+    Plan<int>* cgPlanIntRow;
+    Plan<Vector3d>* cgPlanVectorRow; 
+    Plan<int>* cgPlanIntColumn;
+    Plan<Vector3d>* cgPlanVectorColumn; 
 
+    // work arrays for assembling potential energy
+    vector<potVec> pot_row;
+    vector<potVec> pot_col;
+
     vector<int> identsRow;
     vector<int> identsCol;
 
+    vector<AtomType*> atypesRow;
+    vector<AtomType*> atypesCol;
+
     vector<int> AtomRowToGlobal;
     vector<int> AtomColToGlobal;
 
+public:
     vector<int> cgRowToGlobal;
     vector<int> cgColToGlobal;
 
+private:
+    vector<RealType> groupCutoffRow;
+    vector<RealType> groupCutoffCol;
+    vector<int> groupColToGtype;
+    vector<int> groupRowToGtype;
+
     vector<vector<int> > cellListRow_;
     vector<vector<int> > cellListCol_;
 
@@ -146,7 +170,6 @@ namespace OpenMD {
 
     vector<RealType> massFactorsRow;
     vector<RealType> massFactorsCol;
-
 #endif
 
   };