--- branches/development/src/parallel/ForceDecomposition.hpp	2011/01/11 18:58:12	1538
+++ branches/development/src/parallel/ForceMatrixDecomposition.hpp	2011/08/12 19:59:56	1612
@@ -1,13 +1,6 @@
-/**
- * @file ForceDecomposition.cpp
- * @author Charles Vardeman <cvardema.at.nd.edu>
- * @date 08/18/2010
- * @time 11:56am
- * @version 1.0
+/*
+ * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
- * @section LICENSE
- * Copyright (c) 2010 The University of Notre Dame. All Rights Reserved.
- *
  * The University of Notre Dame grants you ("Licensee") a
  * non-exclusive, royalty free, license to use, modify and
  * redistribute this software in source and binary code form, provided
@@ -45,41 +38,139 @@
  * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
  * [4]  Vardeman & Gezelter, in progress (2009).                        
  */
+ 
+#ifndef PARALLEL_FORCEMATRIXDECOMPOSITION_HPP
+#define PARALLEL_FORCEMATRIXDECOMPOSITION_HPP
 
+#include "parallel/ForceDecomposition.hpp"
+#include "math/SquareMatrix3.hpp"
+#include "brains/Snapshot.hpp"
 
-/*  -*- c++ -*-  */
-#ifndef FORCEDECOMPOSITION_PARALLEL_HPP
-#define FORCEDECOMPOSITION_PARALLEL_HPP
+#ifdef IS_MPI
+#include "parallel/Communicator.hpp"
+#endif
 
+using namespace std;
+namespace OpenMD {
+  
+  class ForceMatrixDecomposition : public ForceDecomposition {
+  public:
+    ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan);
 
-namespace OpenMD{
+    void distributeInitialData();
+    void zeroWorkArrays();
+    void distributeData();
+    void collectIntermediateData();
+    void distributeIntermediateData();
+    void collectData();
 
-class ForceDecomposition
-{
- public:
-  ForceDecomposition();
-  virtual ~ForceDecomposition();
-  void gather();
-  void scatter();
-  
- protected:
- private:
-  int myRank_;
-  static bool isColumn_;
-  static int myIndex_:
-};
+    // neighbor list routines
+    vector<pair<int, int> >  buildNeighborList();
 
+    // group bookkeeping
+    groupCutoffs getGroupCutoffs(int cg1, int cg2);
 
+    // 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);
 
-}
+    // 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);
+    bool excludeAtomPair(int atom1, int atom2);
+    void addForceToAtomRow(int atom1, Vector3d fg);
+    void addForceToAtomColumn(int atom2, Vector3d fg);
 
+    // filling interaction blocks with pointers
+    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> 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<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_;
+
+    vector<vector<int> > groupListRow_;
+    vector<vector<int> > groupListCol_;
+
+    vector<RealType> massFactorsRow;
+    vector<RealType> massFactorsCol;
 #endif
 
+  };
 
+}
+#endif
+