--- branches/development/src/parallel/ForceDecomposition.hpp	2011/05/26 13:55:04	1569
+++ trunk/src/parallel/ForceDecomposition.hpp	2013/08/19 13:12:00	1929
@@ -35,8 +35,9 @@
  *                                                                      
  * [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
@@ -45,10 +46,15 @@
 #include "brains/SimInfo.hpp"
 #include "brains/SnapshotManager.hpp"
 #include "nonbonded/NonBondedInteraction.hpp"
+#include "nonbonded/Cutoffs.hpp"
+#include "nonbonded/InteractionManager.hpp"
+#include "utils/Tuple.hpp"
 
 using namespace std;
 namespace OpenMD {
   
+  typedef tuple3<RealType, RealType, RealType> groupCutoffs;
+
   /**
    * @class ForceDecomposition 
    *
@@ -78,6 +84,10 @@ namespace OpenMD {
    *  |  endif
    *  end
    * collectData                        (parallel communication)
+   * loop over i
+   * | localComputation
+   * end
+   * collectSelfData                    (parallel communication)
    *
    * ForceDecomposition provides the interface for ForceLoop to do the
    * communication steps and to iterate using the correct set of atoms
@@ -86,81 +96,114 @@ namespace OpenMD {
   class ForceDecomposition {
   public:
 
-    ForceDecomposition(SimInfo* info);
+    ForceDecomposition(SimInfo* info, InteractionManager* iMan);
     virtual ~ForceDecomposition() {}
     
     virtual void distributeInitialData() = 0;
     virtual void distributeData() = 0;
+    virtual void zeroWorkArrays() = 0;
     virtual void collectIntermediateData() = 0;
     virtual void distributeIntermediateData() = 0;
     virtual void collectData() = 0;
+    virtual void collectSelfData() = 0;
+    virtual potVec* getEmbeddingPotential() { return &embeddingPot; }
+    virtual potVec* getPairwisePotential() { return &pairwisePot; }
+    virtual potVec* getExcludedPotential() { return &excludedPot; }
+    virtual potVec* getExcludedSelfPotential() { return &excludedSelfPot; }
 
     // neighbor list routines
     virtual bool checkNeighborList();
-    virtual vector<pair<int, int> >  buildNeighborList() = 0;
+    virtual void buildNeighborList(vector<pair<int, int> >& neighborList) = 0;
 
+    // how to handle cutoffs:
+    void setCutoffPolicy(CutoffPolicy cp) {cutoffPolicy_ = cp;}
+    void setUserCutoff(RealType rcut) {userCutoff_ = rcut; userChoseCutoff_ = true; }
+
     // group bookkeeping
-    virtual pair<int, int> getGroupTypes(int cg1, int cg2) = 0;
+    virtual void getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq) = 0;
+    virtual Vector3d& getGroupVelocityColumn(int atom2) = 0;
 
     // Group->atom bookkeeping
-    virtual vector<int> getAtomsInGroupRow(int cg1) = 0; 
-    virtual vector<int> getAtomsInGroupColumn(int cg2) = 0;
+    virtual vector<int>& getAtomsInGroupRow(int cg1) = 0; 
+    virtual vector<int>& getAtomsInGroupColumn(int cg2) = 0;
 
     virtual Vector3d getAtomToGroupVectorRow(int atom1, int cg1) = 0;
     virtual Vector3d getAtomToGroupVectorColumn(int atom2, int cg2) = 0;
-    virtual RealType getMassFactorRow(int atom1) = 0;
-    virtual RealType getMassFactorColumn(int atom2) = 0;
+    virtual RealType& getMassFactorRow(int atom1) = 0;
+    virtual RealType& getMassFactorColumn(int atom2) = 0;
 
     // spatial data
     virtual Vector3d getIntergroupVector(int cg1, int cg2) = 0;
     virtual Vector3d getInteratomicVector(int atom1, int atom2) = 0;
        
     // atom bookkeeping
-    virtual vector<int> getAtomList() = 0;
-    virtual vector<int> getSkipsForAtom(int atom1) = 0;
-    virtual bool skipAtomPair(int atom1, int atom2) = 0;
+    virtual int& getNAtomsInRow() = 0;
+    virtual vector<int>& getExcludesForAtom(int atom1) = 0;
+    virtual bool skipAtomPair(int atom1, int atom2, int cg1, int cg2) = 0;
+    virtual bool excludeAtomPair(int atom1, int atom2) = 0;
+    virtual int getTopologicalDistance(int atom1, int atom2) = 0;
     virtual void addForceToAtomRow(int atom1, Vector3d fg) = 0;
     virtual void addForceToAtomColumn(int atom2, Vector3d fg) = 0;
+    virtual Vector3d& getAtomVelocityColumn(int atom2) = 0;
 
     // filling interaction blocks with pointers
-    virtual InteractionData fillInteractionData(int atom1, int atom2) = 0;
-    virtual InteractionData fillSkipData(int atom1, int atom2) = 0;
-    virtual SelfData fillSelfData(int atom1);
+    virtual void fillInteractionData(InteractionData &idat, int atom1, int atom2) = 0;
+    virtual void unpackInteractionData(InteractionData &idat, int atom1, int atom2) = 0;
+
+    virtual void fillSelfData(SelfData &sdat, int atom1);
+
+    virtual void addToHeatFlux(Vector3d hf);
+    virtual void setHeatFlux(Vector3d hf);
     
   protected:
     SimInfo* info_;   
     SnapshotManager* sman_;    
     Snapshot* snap_;
+    ForceField* ff_;
+    InteractionManager* interactionMan_;
+
     int storageLayout_;
+    bool needVelocities_;
+    bool usePeriodicBoundaryConditions_;
     RealType skinThickness_;   /**< Verlet neighbor list skin thickness */    
+    RealType largestRcut_;
 
-    map<pair<int, int>, int> topoDist; //< topoDist gives the
-                                       //topological distance between
-                                       //two atomic sites.  This
-                                       //declaration is agnostic
-                                       //regarding the parallel
-                                       //decomposition.  The two
-                                       //indices could be local or row
-                                       //& column.  It will be up to
-                                       //the specific decomposition
-                                       //method to fill this.
-    map<pair<int, int>, bool> exclude; //< exclude is the set of pairs
-                                       //to leave out of non-bonded
-                                       //force evaluations.  This
-                                       //declaration is agnostic
-                                       //regarding the parallel
-                                       //decomposition.  The two
-                                       //indices could be local or row
-                                       //& column.  It will be up to
-                                       //the specific decomposition
-                                       //method to fill this.
+    vector<int> idents;
+    vector<int> regions;
+    potVec pairwisePot;
+    potVec embeddingPot;
+    potVec excludedPot;
+    potVec excludedSelfPot;
 
+    /** 
+     * The topological distance between two atomic sites is handled
+     * via two vector structures for speed.  These structures agnostic
+     * regarding the parallel decomposition.  The index for
+     * toposForAtom could be local or row, while the values could be
+     * local or column.  It will be up to the specific decomposition
+     * method to fill these.
+     */
+    vector<vector<int> > toposForAtom; 
+    vector<vector<int> > topoDist;                                       
+    vector<vector<int> > excludesForAtom;
     vector<vector<int> > groupList_;
+    vector<RealType> massFactors;
+    vector<AtomType*> atypesLocal;
+
     vector<Vector3i> cellOffsets_;
     Vector3i nCells_;
     vector<vector<int> > cellList_;
     vector<Vector3d> saved_CG_positions_;
 
+    bool userChoseCutoff_;
+    RealType userCutoff_;
+    CutoffPolicy cutoffPolicy_;
+
+    //map<pair<int, int>, tuple3<RealType, RealType, RealType> > gTypeCutoffMap;
+    vector<vector<RealType> > GrCut;
+    vector<vector<RealType> > GrCutSq;
+    vector<vector<RealType> > GrlistSq;
+
   };    
 }
 #endif