--- branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/04/30 02:54:02	1554
+++ branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/05/12 17:00:14	1562
@@ -343,13 +343,149 @@ namespace OpenMD {
     if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
       idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
       idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
+    }
+#else
+    if (storageLayout_ & DataStorage::dslAmat) {
+      idat.A1 = &(snap_->atomData.aMat[atom1]);
+      idat.A2 = &(snap_->atomData.aMat[atom2]);
+    }
+
+    if (storageLayout_ & DataStorage::dslElectroFrame) {
+      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
+      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
+    }
+
+    if (storageLayout_ & DataStorage::dslTorque) {
+      idat.t1 = &(snap_->atomData.torque[atom1]);
+      idat.t2 = &(snap_->atomData.torque[atom2]);
+    }
+
+    if (storageLayout_ & DataStorage::dslDensity) {
+      idat.rho1 = &(snap_->atomData.density[atom1]);
+      idat.rho2 = &(snap_->atomData.density[atom2]);
     }
+
+    if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
+      idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
+      idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
+    }
 #endif
     
   }
   InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
+    InteractionData idat;
+    skippedCharge1
+      skippedCharge2
+      rij
+      d
+    electroMult
+    sw
+    f
+#ifdef IS_MPI
+
+    if (storageLayout_ & DataStorage::dslElectroFrame) {
+      idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
+      idat.eFrame2 = &(atomColData.electroFrame[atom2]);
+    }
+    if (storageLayout_ & DataStorage::dslTorque) {
+      idat.t1 = &(atomRowData.torque[atom1]);
+      idat.t2 = &(atomColData.torque[atom2]);
+    }
+
+    
   }
   SelfData ForceMatrixDecomposition::fillSelfData(int atom1) {
+  }
+
+
+  /*
+   * buildNeighborList
+   *
+   * first element of pair is row-indexed CutoffGroup
+   * second element of pair is column-indexed CutoffGroup
+   */
+  vector<pair<int, int> >  buildNeighborList() {
+    Vector3d dr, invWid, rs, shift;
+    Vector3i cc, m1v, m2s;
+    RealType rrNebr;
+    int c, j1, j2, m1, m1x, m1y, m1z, m2, n, offset;
+
+
+    vector<pair<int, int> > neighborList;   
+    Vector3i nCells;
+    Vector3d invWid, r;
+
+    rList_ = (rCut_ + skinThickness_);
+    rl2 = rList_ * rList_;
+
+    snap_ = sman_->getCurrentSnapshot();
+    Mat3x3d Hmat = snap_->getHmat();
+    Vector3d Hx = Hmat.getColumn(0);
+    Vector3d Hy = Hmat.getColumn(1);
+    Vector3d Hz = Hmat.getColumn(2);
+
+    nCells.x() = (int) ( Hx.length() )/ rList_;
+    nCells.y() = (int) ( Hy.length() )/ rList_;
+    nCells.z() = (int) ( Hz.length() )/ rList_;
+
+    for (i = 0; i < nGroupsInRow; i++) {
+      rs = cgRowData.position[i];
+      snap_->scaleVector(rs);     
+    }
+    
+
+    VDiv (invWid, cells, region);
+    for (n = nMol; n < nMol + cells.componentProduct(); n ++) cellList[n] = -1;
+    for (n = 0; n < nMol; n ++) {
+      VSAdd (rs, mol[n].r, 0.5, region);
+      VMul (cc, rs, invWid);
+      c = VLinear (cc, cells) + nMol;
+      cellList[n] = cellList[c];
+      cellList[c] = n;
+    }
+    nebrTabLen = 0;
+    for (m1z = 0; m1z < cells.z(); m1z++) {
+      for (m1y = 0; m1y < cells.y(); m1y++) {
+        for (m1x = 0; m1x < cells.x(); m1x++) {
+          Vector3i m1v(m1x, m1y, m1z);
+          m1 = VLinear(m1v, cells) + nMol;
+          for (offset = 0; offset < nOffset_; offset++) {
+            m2v = m1v + cellOffsets_[offset];
+            shift = V3Zero();
+
+            if (m2v.x() >= cells.x) {
+              m2v.x() = 0;           
+              shift.x() = region.x();  
+            } else if (m2v.x() < 0) {
+              m2v.x() = cells.x() - 1; 
+              shift.x() = - region.x();
+            }
+
+            if (m2v.y() >= cells.y()) {
+              m2v.y() = 0;           
+              shift.y() = region.y();  
+            } else if (m2v.y() < 0) {
+              m2v.y() = cells.y() - 1; 
+              shift.y() = - region.y();
+            }
+
+            m2 = VLinear (m2v, cells) + nMol;
+            for (j1 = cellList[m1]; j1 >= 0; j1 = cellList[j1]) {
+              for (j2 = cellList[m2]; j2 >= 0; j2 = cellList[j2]) {
+                if (m1 != m2 || j2 < j1) {
+                  dr = mol[j1].r - mol[j2].r;
+                  VSub (dr, mol[j1].r, mol[j2].r);
+                  VVSub (dr, shift);
+                  if (VLenSq (dr) < rrNebr) {
+                    neighborList.push_back(make_pair(j1, j2));
+                  }
+                }
+              }
+            }
+          }
+        }
+      }
+    }
   }