--- branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/07/08 20:25:32	1587
+++ branches/development/src/parallel/ForceMatrixDecomposition.cpp	2011/07/12 20:33:14	1592
@@ -112,12 +112,30 @@ namespace OpenMD {
     AtomCommIntRow->gather(idents, identsRow);
     AtomCommIntColumn->gather(idents, identsCol);
     
+    // allocate memory for the parallel objects
+    atypesRow.resize(nAtomsInRow_);
+    atypesCol.resize(nAtomsInCol_);
+
+    for (int i = 0; i < nAtomsInRow_; i++) 
+      atypesRow[i] = ff_->getAtomType(identsRow[i]);
+    for (int i = 0; i < nAtomsInCol_; i++) 
+      atypesCol[i] = ff_->getAtomType(identsCol[i]);         
+
+    pot_row.resize(nAtomsInRow_);
+    pot_col.resize(nAtomsInCol_);
+
+    AtomRowToGlobal.resize(nAtomsInRow_);
+    AtomColToGlobal.resize(nAtomsInCol_);
     AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
     AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
     
+    cgRowToGlobal.resize(nGroupsInRow_);
+    cgColToGlobal.resize(nGroupsInCol_);
     cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
     cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
 
+    massFactorsRow.resize(nAtomsInRow_);
+    massFactorsCol.resize(nAtomsInCol_);
     AtomCommRealRow->gather(massFactors, massFactorsRow);
     AtomCommRealColumn->gather(massFactors, massFactorsCol);
 
@@ -176,6 +194,12 @@ namespace OpenMD {
     }
 
 #endif
+
+    // allocate memory for the parallel objects
+    atypesLocal.resize(nLocal_);
+
+    for (int i = 0; i < nLocal_; i++) 
+      atypesLocal[i] = ff_->getAtomType(idents[i]);
 
     groupList_.clear();
     groupList_.resize(nGroups_);
@@ -229,9 +253,11 @@ namespace OpenMD {
   void ForceMatrixDecomposition::createGtypeCutoffMap() {
     
     RealType tol = 1e-6;
+    largestRcut_ = 0.0;
     RealType rc;
     int atid;
     set<AtomType*> atypes = info_->getSimulatedAtomTypes();
+    
     map<int, RealType> atypeCutoff;
       
     for (set<AtomType*>::iterator at = atypes.begin(); 
@@ -239,10 +265,10 @@ namespace OpenMD {
       atid = (*at)->getIdent();
       if (userChoseCutoff_) 
         atypeCutoff[atid] = userCutoff_;
-      else 
+      else
         atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at);
     }
-
+    
     vector<RealType> gTypeCutoffs;
     // first we do a single loop over the cutoff groups to find the
     // largest cutoff for any atypes present in this group.
@@ -302,19 +328,16 @@ namespace OpenMD {
     vector<RealType> groupCutoff(nGroups_, 0.0);
     groupToGtype.resize(nGroups_);
     for (int cg1 = 0; cg1 < nGroups_; cg1++) {
-
       groupCutoff[cg1] = 0.0;
       vector<int> atomList = getAtomsInGroupRow(cg1);
-
       for (vector<int>::iterator ia = atomList.begin(); 
            ia != atomList.end(); ++ia) {            
         int atom1 = (*ia);
         atid = idents[atom1];
-        if (atypeCutoff[atid] > groupCutoff[cg1]) {
+        if (atypeCutoff[atid] > groupCutoff[cg1]) 
           groupCutoff[cg1] = atypeCutoff[atid];
-        }
       }
-
+      
       bool gTypeFound = false;
       for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
         if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
@@ -322,7 +345,7 @@ namespace OpenMD {
           gTypeFound = true;
         } 
       }
-      if (!gTypeFound) {
+      if (!gTypeFound) {      
         gTypeCutoffs.push_back( groupCutoff[cg1] );
         groupToGtype[cg1] = gTypeCutoffs.size() - 1;
       }      
@@ -331,10 +354,12 @@ namespace OpenMD {
 
     // Now we find the maximum group cutoff value present in the simulation
 
-    RealType groupMax = *max_element(gTypeCutoffs.begin(), gTypeCutoffs.end());
+    RealType groupMax = *max_element(gTypeCutoffs.begin(), 
+                                     gTypeCutoffs.end());
 
 #ifdef IS_MPI
-    MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, MPI::MAX);
+    MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE, 
+                              MPI::MAX);
 #endif
     
     RealType tradRcut = groupMax;
@@ -364,13 +389,9 @@ namespace OpenMD {
 
         pair<int,int> key = make_pair(i,j);
         gTypeCutoffMap[key].first = thisRcut;
-
         if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
-
         gTypeCutoffMap[key].second = thisRcut*thisRcut;
-        
         gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
-
         // sanity check
         
         if (userChoseCutoff_) {
@@ -430,8 +451,10 @@ namespace OpenMD {
          Vector<RealType, N_INTERACTION_FAMILIES> (0.0));   
 
     if (storageLayout_ & DataStorage::dslParticlePot) {    
-      fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0);
-      fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0);
+      fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(),
+           0.0);
+      fill(atomColData.particlePot.begin(), atomColData.particlePot.end(),
+           0.0);
     }
 
     if (storageLayout_ & DataStorage::dslDensity) {      
@@ -440,8 +463,10 @@ namespace OpenMD {
     }
 
     if (storageLayout_ & DataStorage::dslFunctional) {   
-      fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0);
-      fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0);
+      fill(atomRowData.functional.begin(), atomRowData.functional.end(),
+           0.0);
+      fill(atomColData.functional.begin(), atomColData.functional.end(),
+           0.0);
     }
 
     if (storageLayout_ & DataStorage::dslFunctionalDerivative) {      
@@ -458,8 +483,9 @@ namespace OpenMD {
            atomColData.skippedCharge.end(), 0.0);
     }
 
-#else
-    
+#endif
+    // even in parallel, we need to zero out the local arrays:
+
     if (storageLayout_ & DataStorage::dslParticlePot) {      
       fill(snap_->atomData.particlePot.begin(), 
            snap_->atomData.particlePot.end(), 0.0);
@@ -481,7 +507,6 @@ namespace OpenMD {
       fill(snap_->atomData.skippedCharge.begin(), 
            snap_->atomData.skippedCharge.end(), 0.0);
     }
-#endif
     
   }
 
@@ -518,6 +543,7 @@ namespace OpenMD {
       AtomCommMatrixColumn->gather(snap_->atomData.electroFrame, 
                                    atomColData.electroFrame);
     }
+
 #endif      
   }
   
@@ -584,8 +610,7 @@ namespace OpenMD {
     AtomCommVectorColumn->scatter(atomColData.force, frc_tmp);
     for (int i = 0; i < n; i++)
       snap_->atomData.force[i] += frc_tmp[i];
-    
-    
+        
     if (storageLayout_ & DataStorage::dslTorque) {
 
       int nt = snap_->atomData.torque.size();
@@ -609,7 +634,7 @@ namespace OpenMD {
 
       AtomCommRealRow->scatter(atomRowData.skippedCharge, skch_tmp);
       for (int i = 0; i < ns; i++) {
-        snap_->atomData.skippedCharge[i] = skch_tmp[i];
+        snap_->atomData.skippedCharge[i] += skch_tmp[i];
         skch_tmp[i] = 0.0;
       }
       
@@ -820,9 +845,9 @@ namespace OpenMD {
     idat.excluded = excludeAtomPair(atom1, atom2);
    
 #ifdef IS_MPI
-    
-    idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), 
-                             ff_->getAtomType(identsCol[atom2]) );
+    idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
+    //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), 
+    //                         ff_->getAtomType(identsCol[atom2]) );
     
     if (storageLayout_ & DataStorage::dslAmat) {
       idat.A1 = &(atomRowData.aMat[atom1]);
@@ -866,8 +891,9 @@ namespace OpenMD {
 
 #else
 
-    idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), 
-                             ff_->getAtomType(idents[atom2]) );
+    idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
+    //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), 
+    //                         ff_->getAtomType(idents[atom2]) );
 
     if (storageLayout_ & DataStorage::dslAmat) {
       idat.A1 = &(snap_->atomData.aMat[atom1]);