--- branches/development/src/parallel/ForceDecomposition.hpp	2011/01/14 22:31:31	1539
+++ trunk/src/parallel/ForceMatrixDecomposition.hpp	2013/08/19 13:12:00	1929
@@ -35,51 +35,149 @@
  *                                                                      
  * [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).                        
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
-#ifndef PARALLEL_FORCEDECOMPOSITION_HPP
-#define PARALLEL_FORCEDECOMPOSITION_HPP
+#ifndef PARALLEL_FORCEMATRIXDECOMPOSITION_HPP
+#define PARALLEL_FORCEMATRIXDECOMPOSITION_HPP
 
-#include "Parallel/Decomposition.hpp"
-#include "Parallel/Communicator.hpp"
+#include "parallel/ForceDecomposition.hpp"
 #include "math/SquareMatrix3.hpp"
+#include "brains/Snapshot.hpp"
 
+#ifdef IS_MPI
+#include "parallel/Communicator.hpp"
+#endif
+
+using namespace std;
 namespace OpenMD {
   
-  class ForceDecomposition : public Decomposition {
+  class ForceMatrixDecomposition : public ForceDecomposition {
   public:
-    ForceDecomposition(Snapshot* sman);
+    ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan);
+
     void distributeInitialData();
+    void zeroWorkArrays();
     void distributeData();
     void collectIntermediateData();
     void distributeIntermediateData();
+    void collectSelfData();
     void collectData();
 
-    unsigned int getNcutoffGroupsI();
-    unsigned int getNcutoffGroupsJ();
+    // neighbor list routines
+    void buildNeighborList(vector<pair<int, int> >& neighborList);
 
-    vector<int> getAtomsInGroupI(int whichCGI);
-    vector<int> getAtomsInGroupJ(int whichCGJ);
+    // group bookkeeping
+    void getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq);
+    Vector3d& getGroupVelocityColumn(int cg2);
 
-    AtomType* getAtomTypeI(int whichAtomI);
-    AtomType* getAtomTypeJ(int whichAtomJ);   
+    // Group->atom bookkeeping
+    vector<int>& getAtomsInGroupRow(int cg1);
+    vector<int>& getAtomsInGroupColumn(int cg2);
+    Vector3d getAtomToGroupVectorRow(int atom1, int cg1);
+    Vector3d getAtomToGroupVectorColumn(int atom2, int cg2);
+    RealType& getMassFactorRow(int atom1);
+    RealType& getMassFactorColumn(int atom2);
 
-#ifdef IS_MPI
-    
-    Comm<I, RealType>* AtomCommRealI; 
-    Comm<I, Vector3d>* AtomCommVectorI; 
-    Comm<I, Mat3x3d>*  AtomCommMatrixI; 
+    // spatial data
+    Vector3d getIntergroupVector(int cg1, int cg2);
+    Vector3d getInteratomicVector(int atom1, int atom2);
+       
+    // atom bookkeeping
+    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);
 
-    Comm<J, RealType>* AtomCommRealJ; 
-    Comm<J, Vector3d>* AtomCommVectorJ; 
-    Comm<J, Mat3x3d>*  AtomCommMatrixJ; 
+    // filling interaction blocks with pointers
+    void fillInteractionData(InteractionData &idat, int atom1, int atom2);
+    void unpackInteractionData(InteractionData &idat, int atom1, int atom2);
 
-    Comm<I, Vector3d>* cgCommVectorI; 
-    Comm<J, Vector3d>* cgCommVectorJ; 
-    
+  private: 
+    void createGtypeCutoffMap();
+
+    int nLocal_;
+    int nGroups_;
+    vector<int> AtomLocalToGlobal;
+    vector<int> cgLocalToGlobal;
+    vector<RealType> groupCutoff;
+    vector<int> groupToGtype;
+
+#ifdef IS_MPI    
+    DataStorage atomRowData;    DataStorage atomColData;
+    DataStorage cgRowData;
+    DataStorage cgColData;
+
+    int nAtomsInRow_; 
+    int nAtomsInCol_;
+    int nGroupsInRow_;
+    int nGroupsInCol_;
+
+    Communicator<Row> rowComm;
+    Communicator<Column> colComm;
+
+    Plan<int>* AtomPlanIntRow;
+    Plan<RealType>* AtomPlanRealRow; 
+    Plan<Vector3d>* AtomPlanVectorRow; 
+    Plan<Mat3x3d>*  AtomPlanMatrixRow; 
+    Plan<potVec>* AtomPlanPotRow;
+
+    Plan<int>* AtomPlanIntColumn;
+    Plan<RealType>* AtomPlanRealColumn; 
+    Plan<Vector3d>* AtomPlanVectorColumn; 
+    Plan<Mat3x3d>*  AtomPlanMatrixColumn;
+    Plan<potVec>* AtomPlanPotColumn; 
+
+    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<potVec> expot_row;
+    vector<potVec> expot_col;
+
+    vector<int> identsRow;
+    vector<int> identsCol;
+
+    vector<int> regionsRow;
+    vector<int> regionsCol;
+
+    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_;
+
+    vector<vector<int> > groupListRow_;
+    vector<vector<int> > groupListCol_;
+
+    vector<RealType> massFactorsRow;
+    vector<RealType> massFactorsCol;
 #endif
+
   };
 
 }