ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/parallel/ForceMatrixDecomposition.cpp
(Generate patch)

Comparing:
branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents), Revision 1601 by gezelter, Thu Aug 4 20:04:35 2011 UTC vs.
trunk/src/parallel/ForceMatrixDecomposition.cpp (file contents), Revision 1896 by gezelter, Tue Jul 2 20:02:31 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42   #include "parallel/ForceMatrixDecomposition.hpp"
43   #include "math/SquareMatrix3.hpp"
# Line 53 | Line 54 | namespace OpenMD {
54      // surrounding cells (not just the 14 upper triangular blocks that
55      // are used when the processor can see all pairs)
56   #ifdef IS_MPI
57 <    cellOffsets_.push_back( Vector3i(-1, 0, 0) );
57 <    cellOffsets_.push_back( Vector3i(-1,-1, 0) );
58 <    cellOffsets_.push_back( Vector3i( 0,-1, 0) );
59 <    cellOffsets_.push_back( Vector3i( 1,-1, 0) );
60 <    cellOffsets_.push_back( Vector3i( 0, 0,-1) );
61 <    cellOffsets_.push_back( Vector3i(-1, 0, 1) );
57 >    cellOffsets_.clear();
58      cellOffsets_.push_back( Vector3i(-1,-1,-1) );
59      cellOffsets_.push_back( Vector3i( 0,-1,-1) );
60 <    cellOffsets_.push_back( Vector3i( 1,-1,-1) );
60 >    cellOffsets_.push_back( Vector3i( 1,-1,-1) );                          
61 >    cellOffsets_.push_back( Vector3i(-1, 0,-1) );
62 >    cellOffsets_.push_back( Vector3i( 0, 0,-1) );
63      cellOffsets_.push_back( Vector3i( 1, 0,-1) );
66 –    cellOffsets_.push_back( Vector3i( 1, 1,-1) );
67 –    cellOffsets_.push_back( Vector3i( 0, 1,-1) );
64      cellOffsets_.push_back( Vector3i(-1, 1,-1) );
65 +    cellOffsets_.push_back( Vector3i( 0, 1,-1) );      
66 +    cellOffsets_.push_back( Vector3i( 1, 1,-1) );
67 +    cellOffsets_.push_back( Vector3i(-1,-1, 0) );
68 +    cellOffsets_.push_back( Vector3i( 0,-1, 0) );
69 +    cellOffsets_.push_back( Vector3i( 1,-1, 0) );
70 +    cellOffsets_.push_back( Vector3i(-1, 0, 0) );      
71 +    cellOffsets_.push_back( Vector3i( 0, 0, 0) );
72 +    cellOffsets_.push_back( Vector3i( 1, 0, 0) );
73 +    cellOffsets_.push_back( Vector3i(-1, 1, 0) );
74 +    cellOffsets_.push_back( Vector3i( 0, 1, 0) );
75 +    cellOffsets_.push_back( Vector3i( 1, 1, 0) );
76 +    cellOffsets_.push_back( Vector3i(-1,-1, 1) );
77 +    cellOffsets_.push_back( Vector3i( 0,-1, 1) );
78 +    cellOffsets_.push_back( Vector3i( 1,-1, 1) );
79 +    cellOffsets_.push_back( Vector3i(-1, 0, 1) );
80 +    cellOffsets_.push_back( Vector3i( 0, 0, 1) );
81 +    cellOffsets_.push_back( Vector3i( 1, 0, 1) );
82 +    cellOffsets_.push_back( Vector3i(-1, 1, 1) );
83 +    cellOffsets_.push_back( Vector3i( 0, 1, 1) );
84 +    cellOffsets_.push_back( Vector3i( 1, 1, 1) );
85   #endif    
86    }
87  
# Line 79 | Line 95 | namespace OpenMD {
95      storageLayout_ = sman_->getStorageLayout();
96      ff_ = info_->getForceField();
97      nLocal_ = snap_->getNumberOfAtoms();
98 <    
98 >  
99      nGroups_ = info_->getNLocalCutoffGroups();
100      // gather the information for atomtype IDs (atids):
101      idents = info_->getIdentArray();
# Line 93 | Line 109 | namespace OpenMD {
109      PairList* oneTwo = info_->getOneTwoInteractions();
110      PairList* oneThree = info_->getOneThreeInteractions();
111      PairList* oneFour = info_->getOneFourInteractions();
112 <
112 >    
113 >    if (needVelocities_)
114 >      snap_->cgData.setStorageLayout(DataStorage::dslPosition |
115 >                                     DataStorage::dslVelocity);
116 >    else
117 >      snap_->cgData.setStorageLayout(DataStorage::dslPosition);
118 >    
119   #ifdef IS_MPI
120  
121      MPI::Intracomm row = rowComm.getComm();
# Line 129 | Line 151 | namespace OpenMD {
151      cgRowData.resize(nGroupsInRow_);
152      cgRowData.setStorageLayout(DataStorage::dslPosition);
153      cgColData.resize(nGroupsInCol_);
154 <    cgColData.setStorageLayout(DataStorage::dslPosition);
155 <        
154 >    if (needVelocities_)
155 >      // we only need column velocities if we need them.
156 >      cgColData.setStorageLayout(DataStorage::dslPosition |
157 >                                 DataStorage::dslVelocity);
158 >    else    
159 >      cgColData.setStorageLayout(DataStorage::dslPosition);
160 >      
161      identsRow.resize(nAtomsInRow_);
162      identsCol.resize(nAtomsInCol_);
163      
# Line 149 | Line 176 | namespace OpenMD {
176      pot_row.resize(nAtomsInRow_);
177      pot_col.resize(nAtomsInCol_);
178  
179 +    expot_row.resize(nAtomsInRow_);
180 +    expot_col.resize(nAtomsInCol_);
181 +
182      AtomRowToGlobal.resize(nAtomsInRow_);
183      AtomColToGlobal.resize(nAtomsInCol_);
184      AtomPlanIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
# Line 218 | Line 248 | namespace OpenMD {
248        }      
249      }
250  
251 < #endif
222 <
223 <    // allocate memory for the parallel objects
224 <    atypesLocal.resize(nLocal_);
225 <
226 <    for (int i = 0; i < nLocal_; i++)
227 <      atypesLocal[i] = ff_->getAtomType(idents[i]);
228 <
229 <    groupList_.clear();
230 <    groupList_.resize(nGroups_);
231 <    for (int i = 0; i < nGroups_; i++) {
232 <      int gid = cgLocalToGlobal[i];
233 <      for (int j = 0; j < nLocal_; j++) {
234 <        int aid = AtomLocalToGlobal[j];
235 <        if (globalGroupMembership[aid] == gid) {
236 <          groupList_[i].push_back(j);
237 <        }
238 <      }      
239 <    }
240 <
251 > #else
252      excludesForAtom.clear();
253      excludesForAtom.resize(nLocal_);
254      toposForAtom.clear();
# Line 270 | Line 281 | namespace OpenMD {
281          }
282        }      
283      }
284 <    
284 > #endif
285 >
286 >    // allocate memory for the parallel objects
287 >    atypesLocal.resize(nLocal_);
288 >
289 >    for (int i = 0; i < nLocal_; i++)
290 >      atypesLocal[i] = ff_->getAtomType(idents[i]);
291 >
292 >    groupList_.clear();
293 >    groupList_.resize(nGroups_);
294 >    for (int i = 0; i < nGroups_; i++) {
295 >      int gid = cgLocalToGlobal[i];
296 >      for (int j = 0; j < nLocal_; j++) {
297 >        int aid = AtomLocalToGlobal[j];
298 >        if (globalGroupMembership[aid] == gid) {
299 >          groupList_[i].push_back(j);
300 >        }
301 >      }      
302 >    }
303 >
304 >
305      createGtypeCutoffMap();
306  
307    }
308    
309    void ForceMatrixDecomposition::createGtypeCutoffMap() {
310      
311 +    GrCut.clear();
312 +    GrCutSq.clear();
313 +    GrlistSq.clear();
314 +
315      RealType tol = 1e-6;
316      largestRcut_ = 0.0;
282 –    RealType rc;
317      int atid;
318      set<AtomType*> atypes = info_->getSimulatedAtomTypes();
319      
# Line 364 | Line 398 | namespace OpenMD {
398        }
399        
400        bool gTypeFound = false;
401 <      for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
401 >      for (unsigned int gt = 0; gt < gTypeCutoffs.size(); gt++) {
402          if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
403            groupToGtype[cg1] = gt;
404            gTypeFound = true;
# Line 389 | Line 423 | namespace OpenMD {
423      
424      RealType tradRcut = groupMax;
425  
426 <    for (int i = 0; i < gTypeCutoffs.size();  i++) {
427 <      for (int j = 0; j < gTypeCutoffs.size();  j++) {      
426 >    GrCut.resize( gTypeCutoffs.size() );
427 >    GrCutSq.resize( gTypeCutoffs.size() );
428 >    GrlistSq.resize( gTypeCutoffs.size() );
429 >
430 >
431 >    for (unsigned int i = 0; i < gTypeCutoffs.size();  i++) {
432 >      GrCut[i].resize( gTypeCutoffs.size() , 0.0);
433 >      GrCutSq[i].resize( gTypeCutoffs.size(), 0.0 );
434 >      GrlistSq[i].resize( gTypeCutoffs.size(), 0.0 );
435 >
436 >      for (unsigned int j = 0; j < gTypeCutoffs.size();  j++) {      
437          RealType thisRcut;
438          switch(cutoffPolicy_) {
439          case TRADITIONAL:
# Line 412 | Line 455 | namespace OpenMD {
455            break;
456          }
457  
458 <        pair<int,int> key = make_pair(i,j);
416 <        gTypeCutoffMap[key].first = thisRcut;
458 >        GrCut[i][j] = thisRcut;
459          if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
460 <        gTypeCutoffMap[key].second = thisRcut*thisRcut;
461 <        gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
460 >        GrCutSq[i][j] = thisRcut * thisRcut;
461 >        GrlistSq[i][j] = pow(thisRcut + skinThickness_, 2);
462 >
463 >        // pair<int,int> key = make_pair(i,j);
464 >        // gTypeCutoffMap[key].first = thisRcut;
465 >        // gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
466          // sanity check
467          
468          if (userChoseCutoff_) {
469 <          if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) {
469 >          if (abs(GrCut[i][j] - userCutoff_) > 0.0001) {
470              sprintf(painCave.errMsg,
471                      "ForceMatrixDecomposition::createGtypeCutoffMap "
472                      "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_);
# Line 433 | Line 479 | namespace OpenMD {
479      }
480    }
481  
482 <
437 <  groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
482 >  void ForceMatrixDecomposition::getGroupCutoffs(int &cg1, int &cg2, RealType &rcut, RealType &rcutsq, RealType &rlistsq) {
483      int i, j;  
484   #ifdef IS_MPI
485      i = groupRowToGtype[cg1];
# Line 443 | Line 488 | namespace OpenMD {
488      i = groupToGtype[cg1];
489      j = groupToGtype[cg2];
490   #endif    
491 <    return gTypeCutoffMap[make_pair(i,j)];
491 >    rcut = GrCut[i][j];
492 >    rcutsq = GrCutSq[i][j];
493 >    rlistsq = GrlistSq[i][j];
494 >    return;
495 >    //return gTypeCutoffMap[make_pair(i,j)];
496    }
497  
498    int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
499 <    for (int j = 0; j < toposForAtom[atom1].size(); j++) {
499 >    for (unsigned int j = 0; j < toposForAtom[atom1].size(); j++) {
500        if (toposForAtom[atom1][j] == atom2)
501          return topoDist[atom1][j];
502 <    }
502 >    }                                          
503      return 0;
504    }
505  
506    void ForceMatrixDecomposition::zeroWorkArrays() {
507      pairwisePot = 0.0;
508      embeddingPot = 0.0;
509 +    excludedPot = 0.0;
510 +    excludedSelfPot = 0.0;
511  
512   #ifdef IS_MPI
513      if (storageLayout_ & DataStorage::dslForce) {
# Line 475 | Line 526 | namespace OpenMD {
526      fill(pot_col.begin(), pot_col.end(),
527           Vector<RealType, N_INTERACTION_FAMILIES> (0.0));  
528  
529 +    fill(expot_row.begin(), expot_row.end(),
530 +         Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
531 +
532 +    fill(expot_col.begin(), expot_col.end(),
533 +         Vector<RealType, N_INTERACTION_FAMILIES> (0.0));  
534 +
535      if (storageLayout_ & DataStorage::dslParticlePot) {    
536        fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(),
537             0.0);
# Line 508 | Line 565 | namespace OpenMD {
565             atomColData.skippedCharge.end(), 0.0);
566      }
567  
568 +    if (storageLayout_ & DataStorage::dslFlucQForce) {      
569 +      fill(atomRowData.flucQFrc.begin(),
570 +           atomRowData.flucQFrc.end(), 0.0);
571 +      fill(atomColData.flucQFrc.begin(),
572 +           atomColData.flucQFrc.end(), 0.0);
573 +    }
574 +
575 +    if (storageLayout_ & DataStorage::dslElectricField) {    
576 +      fill(atomRowData.electricField.begin(),
577 +           atomRowData.electricField.end(), V3Zero);
578 +      fill(atomColData.electricField.begin(),
579 +           atomColData.electricField.end(), V3Zero);
580 +    }
581 +
582   #endif
583      // even in parallel, we need to zero out the local arrays:
584  
# Line 520 | Line 591 | namespace OpenMD {
591        fill(snap_->atomData.density.begin(),
592             snap_->atomData.density.end(), 0.0);
593      }
594 +
595      if (storageLayout_ & DataStorage::dslFunctional) {
596        fill(snap_->atomData.functional.begin(),
597             snap_->atomData.functional.end(), 0.0);
598      }
599 +
600      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {      
601        fill(snap_->atomData.functionalDerivative.begin(),
602             snap_->atomData.functionalDerivative.end(), 0.0);
603      }
604 +
605      if (storageLayout_ & DataStorage::dslSkippedCharge) {      
606        fill(snap_->atomData.skippedCharge.begin(),
607             snap_->atomData.skippedCharge.end(), 0.0);
608      }
609 <    
609 >
610 >    if (storageLayout_ & DataStorage::dslElectricField) {      
611 >      fill(snap_->atomData.electricField.begin(),
612 >           snap_->atomData.electricField.end(), V3Zero);
613 >    }
614    }
615  
616  
# Line 554 | Line 632 | namespace OpenMD {
632  
633      cgPlanVectorColumn->gather(snap_->cgData.position,
634                                 cgColData.position);
635 +
636 +
637 +
638 +    if (needVelocities_) {
639 +      // gather up the atomic velocities
640 +      AtomPlanVectorColumn->gather(snap_->atomData.velocity,
641 +                                   atomColData.velocity);
642 +      
643 +      cgPlanVectorColumn->gather(snap_->cgData.velocity,
644 +                                 cgColData.velocity);
645 +    }
646  
647      
648      // if needed, gather the atomic rotation matrices
# Line 563 | Line 652 | namespace OpenMD {
652        AtomPlanMatrixColumn->gather(snap_->atomData.aMat,
653                                     atomColData.aMat);
654      }
655 <    
656 <    // if needed, gather the atomic eletrostatic frames
657 <    if (storageLayout_ & DataStorage::dslElectroFrame) {
658 <      AtomPlanMatrixRow->gather(snap_->atomData.electroFrame,
659 <                                atomRowData.electroFrame);
660 <      AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame,
661 <                                   atomColData.electroFrame);
655 >
656 >    // if needed, gather the atomic eletrostatic information
657 >    if (storageLayout_ & DataStorage::dslDipole) {
658 >      AtomPlanVectorRow->gather(snap_->atomData.dipole,
659 >                                atomRowData.dipole);
660 >      AtomPlanVectorColumn->gather(snap_->atomData.dipole,
661 >                                   atomColData.dipole);
662      }
663  
664 +    if (storageLayout_ & DataStorage::dslQuadrupole) {
665 +      AtomPlanMatrixRow->gather(snap_->atomData.quadrupole,
666 +                                atomRowData.quadrupole);
667 +      AtomPlanMatrixColumn->gather(snap_->atomData.quadrupole,
668 +                                   atomColData.quadrupole);
669 +    }
670 +        
671 +    // if needed, gather the atomic fluctuating charge values
672 +    if (storageLayout_ & DataStorage::dslFlucQPosition) {
673 +      AtomPlanRealRow->gather(snap_->atomData.flucQPos,
674 +                              atomRowData.flucQPos);
675 +      AtomPlanRealColumn->gather(snap_->atomData.flucQPos,
676 +                                 atomColData.flucQPos);
677 +    }
678 +
679   #endif      
680    }
681    
# Line 594 | Line 698 | namespace OpenMD {
698        for (int i = 0; i < n; i++)
699          snap_->atomData.density[i] += rho_tmp[i];
700      }
701 +
702 +    // this isn't necessary if we don't have polarizable atoms, but
703 +    // we'll leave it here for now.
704 +    if (storageLayout_ & DataStorage::dslElectricField) {
705 +      
706 +      AtomPlanVectorRow->scatter(atomRowData.electricField,
707 +                                 snap_->atomData.electricField);
708 +      
709 +      int n = snap_->atomData.electricField.size();
710 +      vector<Vector3d> field_tmp(n, V3Zero);
711 +      AtomPlanVectorColumn->scatter(atomColData.electricField,
712 +                                    field_tmp);
713 +      for (int i = 0; i < n; i++)
714 +        snap_->atomData.electricField[i] += field_tmp[i];
715 +    }
716   #endif
717    }
718  
# Line 668 | Line 787 | namespace OpenMD {
787        }
788        
789        AtomPlanRealColumn->scatter(atomColData.skippedCharge, skch_tmp);
790 <      for (int i = 0; i < ns; i++)
790 >      for (int i = 0; i < ns; i++)
791          snap_->atomData.skippedCharge[i] += skch_tmp[i];
792 +            
793      }
794      
795 +    if (storageLayout_ & DataStorage::dslFlucQForce) {
796 +
797 +      int nq = snap_->atomData.flucQFrc.size();
798 +      vector<RealType> fqfrc_tmp(nq, 0.0);
799 +
800 +      AtomPlanRealRow->scatter(atomRowData.flucQFrc, fqfrc_tmp);
801 +      for (int i = 0; i < nq; i++) {
802 +        snap_->atomData.flucQFrc[i] += fqfrc_tmp[i];
803 +        fqfrc_tmp[i] = 0.0;
804 +      }
805 +      
806 +      AtomPlanRealColumn->scatter(atomColData.flucQFrc, fqfrc_tmp);
807 +      for (int i = 0; i < nq; i++)
808 +        snap_->atomData.flucQFrc[i] += fqfrc_tmp[i];
809 +            
810 +    }
811 +
812 +    if (storageLayout_ & DataStorage::dslElectricField) {
813 +
814 +      int nef = snap_->atomData.electricField.size();
815 +      vector<Vector3d> efield_tmp(nef, V3Zero);
816 +
817 +      AtomPlanVectorRow->scatter(atomRowData.electricField, efield_tmp);
818 +      for (int i = 0; i < nef; i++) {
819 +        snap_->atomData.electricField[i] += efield_tmp[i];
820 +        efield_tmp[i] = 0.0;
821 +      }
822 +      
823 +      AtomPlanVectorColumn->scatter(atomColData.electricField, efield_tmp);
824 +      for (int i = 0; i < nef; i++)
825 +        snap_->atomData.electricField[i] += efield_tmp[i];
826 +    }
827 +
828 +
829      nLocal_ = snap_->getNumberOfAtoms();
830  
831      vector<potVec> pot_temp(nLocal_,
832                              Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
833 +    vector<potVec> expot_temp(nLocal_,
834 +                              Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
835  
836      // scatter/gather pot_row into the members of my column
837            
838      AtomPlanPotRow->scatter(pot_row, pot_temp);
839 +    AtomPlanPotRow->scatter(expot_row, expot_temp);
840  
841 <    for (int ii = 0;  ii < pot_temp.size(); ii++ )
841 >    for (int ii = 0;  ii < pot_temp.size(); ii++ )
842        pairwisePot += pot_temp[ii];
843 <    
843 >
844 >    for (int ii = 0;  ii < expot_temp.size(); ii++ )
845 >      excludedPot += expot_temp[ii];
846 >        
847 >    if (storageLayout_ & DataStorage::dslParticlePot) {
848 >      // This is the pairwise contribution to the particle pot.  The
849 >      // embedding contribution is added in each of the low level
850 >      // non-bonded routines.  In single processor, this is done in
851 >      // unpackInteractionData, not in collectData.
852 >      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
853 >        for (int i = 0; i < nLocal_; i++) {
854 >          // factor of two is because the total potential terms are divided
855 >          // by 2 in parallel due to row/ column scatter      
856 >          snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii);
857 >        }
858 >      }
859 >    }
860 >
861      fill(pot_temp.begin(), pot_temp.end(),
862           Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
863 +    fill(expot_temp.begin(), expot_temp.end(),
864 +         Vector<RealType, N_INTERACTION_FAMILIES> (0.0));
865        
866      AtomPlanPotColumn->scatter(pot_col, pot_temp);    
867 +    AtomPlanPotColumn->scatter(expot_col, expot_temp);    
868      
869      for (int ii = 0;  ii < pot_temp.size(); ii++ )
870        pairwisePot += pot_temp[ii];    
871 +
872 +    for (int ii = 0;  ii < expot_temp.size(); ii++ )
873 +      excludedPot += expot_temp[ii];    
874 +
875 +    if (storageLayout_ & DataStorage::dslParticlePot) {
876 +      // This is the pairwise contribution to the particle pot.  The
877 +      // embedding contribution is added in each of the low level
878 +      // non-bonded routines.  In single processor, this is done in
879 +      // unpackInteractionData, not in collectData.
880 +      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
881 +        for (int i = 0; i < nLocal_; i++) {
882 +          // factor of two is because the total potential terms are divided
883 +          // by 2 in parallel due to row/ column scatter      
884 +          snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii);
885 +        }
886 +      }
887 +    }
888      
889 +    if (storageLayout_ & DataStorage::dslParticlePot) {
890 +      int npp = snap_->atomData.particlePot.size();
891 +      vector<RealType> ppot_temp(npp, 0.0);
892 +
893 +      // This is the direct or embedding contribution to the particle
894 +      // pot.
895 +      
896 +      AtomPlanRealRow->scatter(atomRowData.particlePot, ppot_temp);
897 +      for (int i = 0; i < npp; i++) {
898 +        snap_->atomData.particlePot[i] += ppot_temp[i];
899 +      }
900 +
901 +      fill(ppot_temp.begin(), ppot_temp.end(), 0.0);
902 +      
903 +      AtomPlanRealColumn->scatter(atomColData.particlePot, ppot_temp);
904 +      for (int i = 0; i < npp; i++) {
905 +        snap_->atomData.particlePot[i] += ppot_temp[i];
906 +      }
907 +    }
908 +
909      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
910        RealType ploc1 = pairwisePot[ii];
911        RealType ploc2 = 0.0;
# Line 699 | Line 913 | namespace OpenMD {
913        pairwisePot[ii] = ploc2;
914      }
915  
916 +    for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
917 +      RealType ploc1 = excludedPot[ii];
918 +      RealType ploc2 = 0.0;
919 +      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
920 +      excludedPot[ii] = ploc2;
921 +    }
922 +
923 +    // Here be dragons.
924 +    MPI::Intracomm col = colComm.getComm();
925 +
926 +    col.Allreduce(MPI::IN_PLACE,
927 +                  &snap_->frameData.conductiveHeatFlux[0], 3,
928 +                  MPI::REALTYPE, MPI::SUM);
929 +
930 +
931   #endif
932  
933    }
934  
935 <  int ForceMatrixDecomposition::getNAtomsInRow() {  
935 >  /**
936 >   * Collects information obtained during the post-pair (and embedding
937 >   * functional) loops onto local data structures.
938 >   */
939 >  void ForceMatrixDecomposition::collectSelfData() {
940 >    snap_ = sman_->getCurrentSnapshot();
941 >    storageLayout_ = sman_->getStorageLayout();
942 >
943   #ifdef IS_MPI
944 +    for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
945 +      RealType ploc1 = embeddingPot[ii];
946 +      RealType ploc2 = 0.0;
947 +      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
948 +      embeddingPot[ii] = ploc2;
949 +    }    
950 +    for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
951 +      RealType ploc1 = excludedSelfPot[ii];
952 +      RealType ploc2 = 0.0;
953 +      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
954 +      excludedSelfPot[ii] = ploc2;
955 +    }    
956 + #endif
957 +    
958 +  }
959 +
960 +
961 +
962 +  int& ForceMatrixDecomposition::getNAtomsInRow() {  
963 + #ifdef IS_MPI
964      return nAtomsInRow_;
965   #else
966      return nLocal_;
# Line 714 | Line 970 | namespace OpenMD {
970    /**
971     * returns the list of atoms belonging to this group.  
972     */
973 <  vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
973 >  vector<int>& ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
974   #ifdef IS_MPI
975      return groupListRow_[cg1];
976   #else
# Line 722 | Line 978 | namespace OpenMD {
978   #endif
979    }
980  
981 <  vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
981 >  vector<int>& ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
982   #ifdef IS_MPI
983      return groupListCol_[cg2];
984   #else
# Line 739 | Line 995 | namespace OpenMD {
995      d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1];
996   #endif
997      
998 <    snap_->wrapVector(d);
998 >    if (usePeriodicBoundaryConditions_) {
999 >      snap_->wrapVector(d);
1000 >    }
1001      return d;    
1002    }
1003  
1004 +  Vector3d& ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){
1005 + #ifdef IS_MPI
1006 +    return cgColData.velocity[cg2];
1007 + #else
1008 +    return snap_->cgData.velocity[cg2];
1009 + #endif
1010 +  }
1011  
1012 +  Vector3d& ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){
1013 + #ifdef IS_MPI
1014 +    return atomColData.velocity[atom2];
1015 + #else
1016 +    return snap_->atomData.velocity[atom2];
1017 + #endif
1018 +  }
1019 +
1020 +
1021    Vector3d ForceMatrixDecomposition::getAtomToGroupVectorRow(int atom1, int cg1){
1022  
1023      Vector3d d;
# Line 753 | Line 1027 | namespace OpenMD {
1027   #else
1028      d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1];
1029   #endif
1030 <
1031 <    snap_->wrapVector(d);
1030 >    if (usePeriodicBoundaryConditions_) {
1031 >      snap_->wrapVector(d);
1032 >    }
1033      return d;    
1034    }
1035    
# Line 766 | Line 1041 | namespace OpenMD {
1041   #else
1042      d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2];
1043   #endif
1044 <    
1045 <    snap_->wrapVector(d);
1044 >    if (usePeriodicBoundaryConditions_) {
1045 >      snap_->wrapVector(d);
1046 >    }
1047      return d;    
1048    }
1049  
1050 <  RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) {
1050 >  RealType& ForceMatrixDecomposition::getMassFactorRow(int atom1) {
1051   #ifdef IS_MPI
1052      return massFactorsRow[atom1];
1053   #else
# Line 779 | Line 1055 | namespace OpenMD {
1055   #endif
1056    }
1057  
1058 <  RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
1058 >  RealType& ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
1059   #ifdef IS_MPI
1060      return massFactorsCol[atom2];
1061   #else
# Line 796 | Line 1072 | namespace OpenMD {
1072   #else
1073      d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1];
1074   #endif
1075 <
1076 <    snap_->wrapVector(d);
1075 >    if (usePeriodicBoundaryConditions_) {
1076 >      snap_->wrapVector(d);
1077 >    }
1078      return d;    
1079    }
1080  
1081 <  vector<int> ForceMatrixDecomposition::getExcludesForAtom(int atom1) {
1081 >  vector<int>& ForceMatrixDecomposition::getExcludesForAtom(int atom1) {
1082      return excludesForAtom[atom1];
1083    }
1084  
# Line 809 | Line 1086 | namespace OpenMD {
1086     * We need to exclude some overcounted interactions that result from
1087     * the parallel decomposition.
1088     */
1089 <  bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
1089 >  bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) {
1090      int unique_id_1, unique_id_2;
1091 <    
1091 >        
1092   #ifdef IS_MPI
1093      // in MPI, we have to look up the unique IDs for each atom
1094      unique_id_1 = AtomRowToGlobal[atom1];
1095      unique_id_2 = AtomColToGlobal[atom2];
1096 +    // group1 = cgRowToGlobal[cg1];
1097 +    // group2 = cgColToGlobal[cg2];
1098 + #else
1099 +    unique_id_1 = AtomLocalToGlobal[atom1];
1100 +    unique_id_2 = AtomLocalToGlobal[atom2];
1101 +    int group1 = cgLocalToGlobal[cg1];
1102 +    int group2 = cgLocalToGlobal[cg2];
1103 + #endif  
1104  
820 –    // this situation should only arise in MPI simulations
1105      if (unique_id_1 == unique_id_2) return true;
1106 <    
1106 >
1107 > #ifdef IS_MPI
1108      // this prevents us from doing the pair on multiple processors
1109      if (unique_id_1 < unique_id_2) {
1110        if ((unique_id_1 + unique_id_2) % 2 == 0) return true;
1111      } else {
1112 <      if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
1112 >      if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
1113      }
1114 + #endif    
1115 +
1116 + #ifndef IS_MPI
1117 +    if (group1 == group2) {
1118 +      if (unique_id_1 < unique_id_2) return true;
1119 +    }
1120   #endif
1121 +    
1122      return false;
1123    }
1124  
# Line 840 | Line 1132 | namespace OpenMD {
1132     * field) must still be handled for these pairs.
1133     */
1134    bool ForceMatrixDecomposition::excludeAtomPair(int atom1, int atom2) {
1135 <    int unique_id_2;
1136 < #ifdef IS_MPI
1137 <    // in MPI, we have to look up the unique IDs for the row atom.
846 <    unique_id_2 = AtomColToGlobal[atom2];
847 < #else
848 <    // in the normal loop, the atom numbers are unique
849 <    unique_id_2 = atom2;
850 < #endif
1135 >
1136 >    // excludesForAtom was constructed to use row/column indices in the MPI
1137 >    // version, and to use local IDs in the non-MPI version:
1138      
1139      for (vector<int>::iterator i = excludesForAtom[atom1].begin();
1140           i != excludesForAtom[atom1].end(); ++i) {
1141 <      if ( (*i) == unique_id_2 ) return true;
1141 >      if ( (*i) == atom2 ) return true;
1142      }
1143  
1144      return false;
# Line 882 | Line 1169 | namespace OpenMD {
1169    
1170   #ifdef IS_MPI
1171      idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
1172 +    idat.atid1 = identsRow[atom1];
1173 +    idat.atid2 = identsCol[atom2];
1174      //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
1175      //                         ff_->getAtomType(identsCol[atom2]) );
1176      
# Line 890 | Line 1179 | namespace OpenMD {
1179        idat.A2 = &(atomColData.aMat[atom2]);
1180      }
1181      
893 –    if (storageLayout_ & DataStorage::dslElectroFrame) {
894 –      idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
895 –      idat.eFrame2 = &(atomColData.electroFrame[atom2]);
896 –    }
897 –
1182      if (storageLayout_ & DataStorage::dslTorque) {
1183        idat.t1 = &(atomRowData.torque[atom1]);
1184        idat.t2 = &(atomColData.torque[atom2]);
1185      }
1186  
1187 +    if (storageLayout_ & DataStorage::dslDipole) {
1188 +      idat.dipole1 = &(atomRowData.dipole[atom1]);
1189 +      idat.dipole2 = &(atomColData.dipole[atom2]);
1190 +    }
1191 +
1192 +    if (storageLayout_ & DataStorage::dslQuadrupole) {
1193 +      idat.quadrupole1 = &(atomRowData.quadrupole[atom1]);
1194 +      idat.quadrupole2 = &(atomColData.quadrupole[atom2]);
1195 +    }
1196 +
1197      if (storageLayout_ & DataStorage::dslDensity) {
1198        idat.rho1 = &(atomRowData.density[atom1]);
1199        idat.rho2 = &(atomColData.density[atom2]);
# Line 925 | Line 1219 | namespace OpenMD {
1219        idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
1220      }
1221  
1222 < #else
1222 >    if (storageLayout_ & DataStorage::dslFlucQPosition) {              
1223 >      idat.flucQ1 = &(atomRowData.flucQPos[atom1]);
1224 >      idat.flucQ2 = &(atomColData.flucQPos[atom2]);
1225 >    }
1226  
1227 + #else
1228 +    
1229      idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
1230 <    //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]),
1231 <    //                         ff_->getAtomType(idents[atom2]) );
1230 >    idat.atid1 = idents[atom1];
1231 >    idat.atid2 = idents[atom2];
1232  
1233      if (storageLayout_ & DataStorage::dslAmat) {
1234        idat.A1 = &(snap_->atomData.aMat[atom1]);
1235        idat.A2 = &(snap_->atomData.aMat[atom2]);
1236      }
1237  
939 –    if (storageLayout_ & DataStorage::dslElectroFrame) {
940 –      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
941 –      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
942 –    }
943 –
1238      if (storageLayout_ & DataStorage::dslTorque) {
1239        idat.t1 = &(snap_->atomData.torque[atom1]);
1240        idat.t2 = &(snap_->atomData.torque[atom2]);
1241      }
1242  
1243 +    if (storageLayout_ & DataStorage::dslDipole) {
1244 +      idat.dipole1 = &(snap_->atomData.dipole[atom1]);
1245 +      idat.dipole2 = &(snap_->atomData.dipole[atom2]);
1246 +    }
1247 +
1248 +    if (storageLayout_ & DataStorage::dslQuadrupole) {
1249 +      idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]);
1250 +      idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]);
1251 +    }
1252 +
1253      if (storageLayout_ & DataStorage::dslDensity) {    
1254        idat.rho1 = &(snap_->atomData.density[atom1]);
1255        idat.rho2 = &(snap_->atomData.density[atom2]);
# Line 970 | Line 1274 | namespace OpenMD {
1274        idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
1275        idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
1276      }
1277 +
1278 +    if (storageLayout_ & DataStorage::dslFlucQPosition) {              
1279 +      idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]);
1280 +      idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]);
1281 +    }
1282 +
1283   #endif
1284    }
1285  
1286    
1287    void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) {    
1288   #ifdef IS_MPI
1289 <    pot_row[atom1] += 0.5 *  *(idat.pot);
1290 <    pot_col[atom2] += 0.5 *  *(idat.pot);
1289 >    pot_row[atom1] += RealType(0.5) *  *(idat.pot);
1290 >    pot_col[atom2] += RealType(0.5) *  *(idat.pot);
1291 >    expot_row[atom1] += RealType(0.5) *  *(idat.excludedPot);
1292 >    expot_col[atom2] += RealType(0.5) *  *(idat.excludedPot);
1293  
1294      atomRowData.force[atom1] += *(idat.f1);
1295      atomColData.force[atom2] -= *(idat.f1);
1296 +
1297 +    if (storageLayout_ & DataStorage::dslFlucQForce) {              
1298 +      atomRowData.flucQFrc[atom1] -= *(idat.dVdFQ1);
1299 +      atomColData.flucQFrc[atom2] -= *(idat.dVdFQ2);
1300 +    }
1301 +
1302 +    if (storageLayout_ & DataStorage::dslElectricField) {              
1303 +      atomRowData.electricField[atom1] += *(idat.eField1);
1304 +      atomColData.electricField[atom2] += *(idat.eField2);
1305 +    }
1306 +
1307   #else
1308      pairwisePot += *(idat.pot);
1309 +    excludedPot += *(idat.excludedPot);
1310  
1311      snap_->atomData.force[atom1] += *(idat.f1);
1312      snap_->atomData.force[atom2] -= *(idat.f1);
1313 +
1314 +    if (idat.doParticlePot) {
1315 +      // This is the pairwise contribution to the particle pot.  The
1316 +      // embedding contribution is added in each of the low level
1317 +      // non-bonded routines.  In parallel, this calculation is done
1318 +      // in collectData, not in unpackInteractionData.
1319 +      snap_->atomData.particlePot[atom1] += *(idat.vpair) * *(idat.sw);
1320 +      snap_->atomData.particlePot[atom2] += *(idat.vpair) * *(idat.sw);
1321 +    }
1322 +    
1323 +    if (storageLayout_ & DataStorage::dslFlucQForce) {              
1324 +      snap_->atomData.flucQFrc[atom1] -= *(idat.dVdFQ1);
1325 +      snap_->atomData.flucQFrc[atom2] -= *(idat.dVdFQ2);
1326 +    }
1327 +
1328 +    if (storageLayout_ & DataStorage::dslElectricField) {              
1329 +      snap_->atomData.electricField[atom1] += *(idat.eField1);
1330 +      snap_->atomData.electricField[atom2] += *(idat.eField2);
1331 +    }
1332 +
1333   #endif
1334      
1335    }
# Line 996 | Line 1340 | namespace OpenMD {
1340     * first element of pair is row-indexed CutoffGroup
1341     * second element of pair is column-indexed CutoffGroup
1342     */
1343 <  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
1344 <      
1345 <    vector<pair<int, int> > neighborList;
1343 >  void ForceMatrixDecomposition::buildNeighborList(vector<pair<int,int> >& neighborList) {
1344 >    
1345 >    neighborList.clear();
1346      groupCutoffs cuts;
1347      bool doAllPairs = false;
1348  
1349 +    RealType rList_ = (largestRcut_ + skinThickness_);
1350 +    RealType rcut, rcutsq, rlistsq;
1351 +    Snapshot* snap_ = sman_->getCurrentSnapshot();
1352 +    Mat3x3d box;
1353 +    Mat3x3d invBox;
1354 +
1355 +    Vector3d rs, scaled, dr;
1356 +    Vector3i whichCell;
1357 +    int cellIndex;
1358 +
1359   #ifdef IS_MPI
1360      cellListRow_.clear();
1361      cellListCol_.clear();
1362   #else
1363      cellList_.clear();
1364   #endif
1365 <
1366 <    RealType rList_ = (largestRcut_ + skinThickness_);
1367 <    RealType rl2 = rList_ * rList_;
1368 <    Snapshot* snap_ = sman_->getCurrentSnapshot();
1369 <    Mat3x3d Hmat = snap_->getHmat();
1370 <    Vector3d Hx = Hmat.getColumn(0);
1371 <    Vector3d Hy = Hmat.getColumn(1);
1372 <    Vector3d Hz = Hmat.getColumn(2);
1373 <
1374 <    nCells_.x() = (int) ( Hx.length() )/ rList_;
1375 <    nCells_.y() = (int) ( Hy.length() )/ rList_;
1376 <    nCells_.z() = (int) ( Hz.length() )/ rList_;
1377 <
1365 >    
1366 >    if (!usePeriodicBoundaryConditions_) {
1367 >      box = snap_->getBoundingBox();
1368 >      invBox = snap_->getInvBoundingBox();
1369 >    } else {
1370 >      box = snap_->getHmat();
1371 >      invBox = snap_->getInvHmat();
1372 >    }
1373 >    
1374 >    Vector3d boxX = box.getColumn(0);
1375 >    Vector3d boxY = box.getColumn(1);
1376 >    Vector3d boxZ = box.getColumn(2);
1377 >    
1378 >    nCells_.x() = (int) ( boxX.length() )/ rList_;
1379 >    nCells_.y() = (int) ( boxY.length() )/ rList_;
1380 >    nCells_.z() = (int) ( boxZ.length() )/ rList_;
1381 >    
1382      // handle small boxes where the cell offsets can end up repeating cells
1383      
1384      if (nCells_.x() < 3) doAllPairs = true;
1385      if (nCells_.y() < 3) doAllPairs = true;
1386      if (nCells_.z() < 3) doAllPairs = true;
1387 <
1030 <    Mat3x3d invHmat = snap_->getInvHmat();
1031 <    Vector3d rs, scaled, dr;
1032 <    Vector3i whichCell;
1033 <    int cellIndex;
1387 >    
1388      int nCtot = nCells_.x() * nCells_.y() * nCells_.z();
1389 <
1389 >    
1390   #ifdef IS_MPI
1391      cellListRow_.resize(nCtot);
1392      cellListCol_.resize(nCtot);
1393   #else
1394      cellList_.resize(nCtot);
1395   #endif
1396 <
1396 >    
1397      if (!doAllPairs) {
1398   #ifdef IS_MPI
1399 <
1399 >      
1400        for (int i = 0; i < nGroupsInRow_; i++) {
1401          rs = cgRowData.position[i];
1402          
1403          // scaled positions relative to the box vectors
1404 <        scaled = invHmat * rs;
1404 >        scaled = invBox * rs;
1405          
1406          // wrap the vector back into the unit box by subtracting integer box
1407          // numbers
1408          for (int j = 0; j < 3; j++) {
1409            scaled[j] -= roundMe(scaled[j]);
1410            scaled[j] += 0.5;
1411 +          // Handle the special case when an object is exactly on the
1412 +          // boundary (a scaled coordinate of 1.0 is the same as
1413 +          // scaled coordinate of 0.0)
1414 +          if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1415          }
1416          
1417          // find xyz-indices of cell that cutoffGroup is in.
# Line 1071 | Line 1429 | namespace OpenMD {
1429          rs = cgColData.position[i];
1430          
1431          // scaled positions relative to the box vectors
1432 <        scaled = invHmat * rs;
1432 >        scaled = invBox * rs;
1433          
1434          // wrap the vector back into the unit box by subtracting integer box
1435          // numbers
1436          for (int j = 0; j < 3; j++) {
1437            scaled[j] -= roundMe(scaled[j]);
1438            scaled[j] += 0.5;
1439 +          // Handle the special case when an object is exactly on the
1440 +          // boundary (a scaled coordinate of 1.0 is the same as
1441 +          // scaled coordinate of 0.0)
1442 +          if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1443          }
1444          
1445          // find xyz-indices of cell that cutoffGroup is in.
# Line 1091 | Line 1453 | namespace OpenMD {
1453          // add this cutoff group to the list of groups in this cell;
1454          cellListCol_[cellIndex].push_back(i);
1455        }
1456 +      
1457   #else
1458        for (int i = 0; i < nGroups_; i++) {
1459          rs = snap_->cgData.position[i];
1460          
1461          // scaled positions relative to the box vectors
1462 <        scaled = invHmat * rs;
1462 >        scaled = invBox * rs;
1463          
1464          // wrap the vector back into the unit box by subtracting integer box
1465          // numbers
1466          for (int j = 0; j < 3; j++) {
1467            scaled[j] -= roundMe(scaled[j]);
1468            scaled[j] += 0.5;
1469 +          // Handle the special case when an object is exactly on the
1470 +          // boundary (a scaled coordinate of 1.0 is the same as
1471 +          // scaled coordinate of 0.0)
1472 +          if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1473          }
1474          
1475          // find xyz-indices of cell that cutoffGroup is in.
# Line 1116 | Line 1483 | namespace OpenMD {
1483          // add this cutoff group to the list of groups in this cell;
1484          cellList_[cellIndex].push_back(i);
1485        }
1486 +
1487   #endif
1488  
1489        for (int m1z = 0; m1z < nCells_.z(); m1z++) {
# Line 1128 | Line 1496 | namespace OpenMD {
1496                   os != cellOffsets_.end(); ++os) {
1497                
1498                Vector3i m2v = m1v + (*os);
1499 <              
1499 >            
1500 >
1501                if (m2v.x() >= nCells_.x()) {
1502                  m2v.x() = 0;          
1503                } else if (m2v.x() < 0) {
# Line 1146 | Line 1515 | namespace OpenMD {
1515                } else if (m2v.z() < 0) {
1516                  m2v.z() = nCells_.z() - 1;
1517                }
1518 <              
1518 >
1519                int m2 = Vlinear (m2v, nCells_);
1520                
1521   #ifdef IS_MPI
# Line 1155 | Line 1524 | namespace OpenMD {
1524                  for (vector<int>::iterator j2 = cellListCol_[m2].begin();
1525                       j2 != cellListCol_[m2].end(); ++j2) {
1526                    
1527 <                  // In parallel, we need to visit *all* pairs of row &
1528 <                  // column indicies and will truncate later on.
1527 >                  // In parallel, we need to visit *all* pairs of row
1528 >                  // & column indicies and will divide labor in the
1529 >                  // force evaluation later.
1530                    dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)];
1531 <                  snap_->wrapVector(dr);
1532 <                  cuts = getGroupCutoffs( (*j1), (*j2) );
1533 <                  if (dr.lengthSquare() < cuts.third) {
1531 >                  if (usePeriodicBoundaryConditions_) {
1532 >                    snap_->wrapVector(dr);
1533 >                  }
1534 >                  getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq );
1535 >                  if (dr.lengthSquare() < rlistsq) {
1536                      neighborList.push_back(make_pair((*j1), (*j2)));
1537                    }                  
1538                  }
1539                }
1540   #else
1169 –              
1541                for (vector<int>::iterator j1 = cellList_[m1].begin();
1542                     j1 != cellList_[m1].end(); ++j1) {
1543                  for (vector<int>::iterator j2 = cellList_[m2].begin();
1544                       j2 != cellList_[m2].end(); ++j2) {
1545 <                  
1545 >    
1546                    // Always do this if we're in different cells or if
1547 <                  // we're in the same cell and the global index of the
1548 <                  // j2 cutoff group is less than the j1 cutoff group
1549 <                  
1550 <                  if (m2 != m1 || (*j2) < (*j1)) {
1547 >                  // we're in the same cell and the global index of
1548 >                  // the j2 cutoff group is greater than or equal to
1549 >                  // the j1 cutoff group.  Note that Rappaport's code
1550 >                  // has a "less than" conditional here, but that
1551 >                  // deals with atom-by-atom computation.  OpenMD
1552 >                  // allows atoms within a single cutoff group to
1553 >                  // interact with each other.
1554 >
1555 >                  if (m2 != m1 || (*j2) >= (*j1) ) {
1556 >
1557                      dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
1558 <                    snap_->wrapVector(dr);
1559 <                    cuts = getGroupCutoffs( (*j1), (*j2) );
1560 <                    if (dr.lengthSquare() < cuts.third) {
1558 >                    if (usePeriodicBoundaryConditions_) {
1559 >                      snap_->wrapVector(dr);
1560 >                    }
1561 >                    getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq );
1562 >                    if (dr.lengthSquare() < rlistsq) {
1563                        neighborList.push_back(make_pair((*j1), (*j2)));
1564                      }
1565                    }
# Line 1195 | Line 1574 | namespace OpenMD {
1574        // branch to do all cutoff group pairs
1575   #ifdef IS_MPI
1576        for (int j1 = 0; j1 < nGroupsInRow_; j1++) {
1577 <        for (int j2 = 0; j2 < nGroupsInCol_; j2++) {      
1577 >        for (int j2 = 0; j2 < nGroupsInCol_; j2++) {    
1578            dr = cgColData.position[j2] - cgRowData.position[j1];
1579 <          snap_->wrapVector(dr);
1580 <          cuts = getGroupCutoffs( j1, j2 );
1581 <          if (dr.lengthSquare() < cuts.third) {
1579 >          if (usePeriodicBoundaryConditions_) {
1580 >            snap_->wrapVector(dr);
1581 >          }
1582 >          getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq);
1583 >          if (dr.lengthSquare() < rlistsq) {
1584              neighborList.push_back(make_pair(j1, j2));
1585            }
1586          }
1587 <      }
1587 >      }      
1588   #else
1589 <      for (int j1 = 0; j1 < nGroups_ - 1; j1++) {
1590 <        for (int j2 = j1 + 1; j2 < nGroups_; j2++) {
1589 >      // include all groups here.
1590 >      for (int j1 = 0; j1 < nGroups_; j1++) {
1591 >        // include self group interactions j2 == j1
1592 >        for (int j2 = j1; j2 < nGroups_; j2++) {
1593            dr = snap_->cgData.position[j2] - snap_->cgData.position[j1];
1594 <          snap_->wrapVector(dr);
1595 <          cuts = getGroupCutoffs( j1, j2 );
1596 <          if (dr.lengthSquare() < cuts.third) {
1594 >          if (usePeriodicBoundaryConditions_) {
1595 >            snap_->wrapVector(dr);
1596 >          }
1597 >          getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq );
1598 >          if (dr.lengthSquare() < rlistsq) {
1599              neighborList.push_back(make_pair(j1, j2));
1600            }
1601 <        }
1602 <      }        
1601 >        }    
1602 >      }
1603   #endif
1604      }
1605        
# Line 1223 | Line 1608 | namespace OpenMD {
1608      saved_CG_positions_.clear();
1609      for (int i = 0; i < nGroups_; i++)
1610        saved_CG_positions_.push_back(snap_->cgData.position[i]);
1226 –    
1227 –    return neighborList;
1611    }
1612   } //end namespace OpenMD

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines