--- trunk/src/parallel/ForceDecomposition.cpp	2013/08/19 13:12:00	1929
+++ trunk/src/parallel/ForceDecomposition.cpp	2015/03/07 21:41:51	2071
@@ -39,12 +39,14 @@
  * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
  * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
-#include "parallel/ForceDecomposition.hpp"
-#include "math/Vector3.hpp"
+
 #ifdef IS_MPI
 #include <mpi.h>
 #endif
 
+#include "parallel/ForceDecomposition.hpp"
+#include "math/Vector3.hpp"
+
 using namespace std;
 namespace OpenMD {
 
@@ -53,7 +55,6 @@ namespace OpenMD {
     sman_ = info_->getSnapshotManager();
     storageLayout_ = sman_->getStorageLayout();
     ff_ = info_->getForceField();
-    userChoseCutoff_ = false;
 
     usePeriodicBoundaryConditions_ = info->getSimParams()->getUsePeriodicBoundaryConditions();
 
@@ -96,9 +97,14 @@ namespace OpenMD {
     cellOffsets_.push_back( Vector3i(1, -1,1) );
   }
 
+  void ForceDecomposition::setCutoffRadius(RealType rcut) {
+    rCut_ = rcut;
+    rList_ = rCut_ + skinThickness_;
+    rListSq_ = rList_ * rList_;
+  }
+
   void ForceDecomposition::fillSelfData(SelfData &sdat, int atom1) {
 
-    sdat.atype = atypesLocal[atom1];
     sdat.atid = idents[atom1];
 
     sdat.pot = &embeddingPot;
@@ -146,37 +152,31 @@ namespace OpenMD {
   }
 
   bool ForceDecomposition::checkNeighborList() {
-
-    int nGroups = snap_->cgData.position.size();
+    RealType st2 = pow( skinThickness_ / 2.0, 2);
+    std::size_t nGroups = snap_->cgData.position.size();
     if (needVelocities_) 
-      snap_->cgData.setStorageLayout(DataStorage::dslPosition | DataStorage::dslVelocity);
+      snap_->cgData.setStorageLayout(DataStorage::dslPosition |
+                                     DataStorage::dslVelocity);
     
     // if we have changed the group identities or haven't set up the
     // saved positions we automatically will need a neighbor list update:
-
+    
     if ( saved_CG_positions_.size() != nGroups ) return true;
 
     RealType dispmax = 0.0;
     Vector3d disp;    
 
-    for (int i = 0; i < nGroups; i++) {
+    for (std::size_t i = 0; i < nGroups; i++) {
       disp = snap_->cgData.position[i]  - saved_CG_positions_[i];
-      for (int j = 0; j < 3; j++)
-        dispmax = max( abs(disp[j]), dispmax);
+      dispmax = max(dispmax, disp.lengthSquare());
     }
 
 #ifdef IS_MPI
-    MPI::COMM_WORLD.Allreduce(&dispmax, &dispmax, 1, MPI::REALTYPE, MPI::MAX);
+    MPI_Allreduce(MPI_IN_PLACE, &dispmax, 1, MPI_REALTYPE, MPI_MAX, 
+                  MPI_COMM_WORLD);
 #endif
 
-    // a conservative test of list skin crossings
-    dispmax = 2.0 * sqrt (3.0 * dispmax * dispmax);
-
-
-    if (dispmax > skinThickness_) 
-      return (dispmax > skinThickness_);   
-
-    return false;
+    return (dispmax > st2) ? true : false;
   }
 
   void ForceDecomposition::addToHeatFlux(Vector3d hf) {