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 1562 by gezelter, Thu May 12 17:00:14 2011 UTC vs.
Revision 1570 by gezelter, Thu May 26 21:56:04 2011 UTC

# Line 42 | Line 42
42   #include "math/SquareMatrix3.hpp"
43   #include "nonbonded/NonBondedInteraction.hpp"
44   #include "brains/SnapshotManager.hpp"
45 + #include "brains/PairList.hpp"
46  
47   using namespace std;
48   namespace OpenMD {
# Line 54 | Line 55 | namespace OpenMD {
55    void ForceMatrixDecomposition::distributeInitialData() {
56      snap_ = sman_->getCurrentSnapshot();
57      storageLayout_ = sman_->getStorageLayout();
58 < #ifdef IS_MPI    
59 <    int nLocal = snap_->getNumberOfAtoms();
59 <    int nGroups = snap_->getNumberOfCutoffGroups();
60 <    
61 <    AtomCommIntRow = new Communicator<Row,int>(nLocal);
62 <    AtomCommRealRow = new Communicator<Row,RealType>(nLocal);
63 <    AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal);
64 <    AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal);
58 >    nLocal_ = snap_->getNumberOfAtoms();
59 >    nGroups_ = snap_->getNumberOfCutoffGroups();
60  
61 <    AtomCommIntColumn = new Communicator<Column,int>(nLocal);
62 <    AtomCommRealColumn = new Communicator<Column,RealType>(nLocal);
63 <    AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal);
64 <    AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal);
61 >    // gather the information for atomtype IDs (atids):
62 >    vector<int> identsLocal = info_->getIdentArray();
63 >    AtomLocalToGlobal = info_->getGlobalAtomIndices();
64 >    cgLocalToGlobal = info_->getGlobalGroupIndices();
65 >    vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
66 >    vector<RealType> massFactorsLocal = info_->getMassFactors();
67 >    PairList excludes = info_->getExcludedInteractions();
68 >    PairList oneTwo = info_->getOneTwoInteractions();
69 >    PairList oneThree = info_->getOneThreeInteractions();
70 >    PairList oneFour = info_->getOneFourInteractions();
71 >    vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
72  
73 <    cgCommIntRow = new Communicator<Row,int>(nGroups);
74 <    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups);
75 <    cgCommIntColumn = new Communicator<Column,int>(nGroups);
76 <    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups);
73 > #ifdef IS_MPI
74 >
75 >    AtomCommIntRow = new Communicator<Row,int>(nLocal_);
76 >    AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
77 >    AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
78 >    AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
79  
80 <    int nAtomsInRow = AtomCommIntRow->getSize();
81 <    int nAtomsInCol = AtomCommIntColumn->getSize();
82 <    int nGroupsInRow = cgCommIntRow->getSize();
83 <    int nGroupsInCol = cgCommIntColumn->getSize();
80 >    AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
81 >    AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
82 >    AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
83 >    AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
84  
85 +    cgCommIntRow = new Communicator<Row,int>(nGroups_);
86 +    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
87 +    cgCommIntColumn = new Communicator<Column,int>(nGroups_);
88 +    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_);
89 +
90 +    nAtomsInRow_ = AtomCommIntRow->getSize();
91 +    nAtomsInCol_ = AtomCommIntColumn->getSize();
92 +    nGroupsInRow_ = cgCommIntRow->getSize();
93 +    nGroupsInCol_ = cgCommIntColumn->getSize();
94 +
95      // Modify the data storage objects with the correct layouts and sizes:
96 <    atomRowData.resize(nAtomsInRow);
96 >    atomRowData.resize(nAtomsInRow_);
97      atomRowData.setStorageLayout(storageLayout_);
98 <    atomColData.resize(nAtomsInCol);
98 >    atomColData.resize(nAtomsInCol_);
99      atomColData.setStorageLayout(storageLayout_);
100 <    cgRowData.resize(nGroupsInRow);
100 >    cgRowData.resize(nGroupsInRow_);
101      cgRowData.setStorageLayout(DataStorage::dslPosition);
102 <    cgColData.resize(nGroupsInCol);
102 >    cgColData.resize(nGroupsInCol_);
103      cgColData.setStorageLayout(DataStorage::dslPosition);
104      
105      vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
106 <                                      vector<RealType> (nAtomsInRow, 0.0));
106 >                                      vector<RealType> (nAtomsInRow_, 0.0));
107      vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
108 <                                      vector<RealType> (nAtomsInCol, 0.0));
95 <
96 <
97 <    vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
108 >                                      vector<RealType> (nAtomsInCol_, 0.0));
109      
110 <    // gather the information for atomtype IDs (atids):
111 <    vector<int> identsLocal = info_->getIdentArray();
101 <    identsRow.reserve(nAtomsInRow);
102 <    identsCol.reserve(nAtomsInCol);
110 >    identsRow.reserve(nAtomsInRow_);
111 >    identsCol.reserve(nAtomsInCol_);
112      
113      AtomCommIntRow->gather(identsLocal, identsRow);
114      AtomCommIntColumn->gather(identsLocal, identsCol);
115      
107 –    AtomLocalToGlobal = info_->getGlobalAtomIndices();
116      AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
117      AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
118      
111 –    cgLocalToGlobal = info_->getGlobalGroupIndices();
119      cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
120      cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
121  
122 <    // still need:
123 <    // topoDist
124 <    // exclude
122 >    AtomCommRealRow->gather(massFactorsLocal, massFactorsRow);
123 >    AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol);
124 >
125 >    groupListRow_.clear();
126 >    groupListRow_.reserve(nGroupsInRow_);
127 >    for (int i = 0; i < nGroupsInRow_; i++) {
128 >      int gid = cgRowToGlobal[i];
129 >      for (int j = 0; j < nAtomsInRow_; j++) {
130 >        int aid = AtomRowToGlobal[j];
131 >        if (globalGroupMembership[aid] == gid)
132 >          groupListRow_[i].push_back(j);
133 >      }      
134 >    }
135 >
136 >    groupListCol_.clear();
137 >    groupListCol_.reserve(nGroupsInCol_);
138 >    for (int i = 0; i < nGroupsInCol_; i++) {
139 >      int gid = cgColToGlobal[i];
140 >      for (int j = 0; j < nAtomsInCol_; j++) {
141 >        int aid = AtomColToGlobal[j];
142 >        if (globalGroupMembership[aid] == gid)
143 >          groupListCol_[i].push_back(j);
144 >      }      
145 >    }
146 >
147 >    skipsForRowAtom.clear();
148 >    skipsForRowAtom.reserve(nAtomsInRow_);
149 >    for (int i = 0; i < nAtomsInRow_; i++) {
150 >      int iglob = AtomColToGlobal[i];
151 >      for (int j = 0; j < nAtomsInCol_; j++) {
152 >        int jglob = AtomRowToGlobal[j];        
153 >        if (excludes.hasPair(iglob, jglob))
154 >          skipsForRowAtom[i].push_back(j);      
155 >      }      
156 >    }
157 >
158 >    toposForRowAtom.clear();
159 >    toposForRowAtom.reserve(nAtomsInRow_);
160 >    for (int i = 0; i < nAtomsInRow_; i++) {
161 >      int iglob = AtomColToGlobal[i];
162 >      int nTopos = 0;
163 >      for (int j = 0; j < nAtomsInCol_; j++) {
164 >        int jglob = AtomRowToGlobal[j];        
165 >        if (oneTwo.hasPair(iglob, jglob)) {
166 >          toposForRowAtom[i].push_back(j);
167 >          topoDistRow[i][nTopos] = 1;
168 >          nTopos++;
169 >        }
170 >        if (oneThree.hasPair(iglob, jglob)) {
171 >          toposForRowAtom[i].push_back(j);
172 >          topoDistRow[i][nTopos] = 2;
173 >          nTopos++;
174 >        }
175 >        if (oneFour.hasPair(iglob, jglob)) {
176 >          toposForRowAtom[i].push_back(j);
177 >          topoDistRow[i][nTopos] = 3;
178 >          nTopos++;
179 >        }
180 >      }      
181 >    }
182 >
183   #endif
184 <  }
185 <    
184 >
185 >    groupList_.clear();
186 >    groupList_.reserve(nGroups_);
187 >    for (int i = 0; i < nGroups_; i++) {
188 >      int gid = cgLocalToGlobal[i];
189 >      for (int j = 0; j < nLocal_; j++) {
190 >        int aid = AtomLocalToGlobal[j];
191 >        if (globalGroupMembership[aid] == gid)
192 >          groupList_[i].push_back(j);
193 >      }      
194 >    }
195  
196 +    skipsForLocalAtom.clear();
197 +    skipsForLocalAtom.reserve(nLocal_);
198  
199 +    for (int i = 0; i < nLocal_; i++) {
200 +      int iglob = AtomLocalToGlobal[i];
201 +      for (int j = 0; j < nLocal_; j++) {
202 +        int jglob = AtomLocalToGlobal[j];        
203 +        if (excludes.hasPair(iglob, jglob))
204 +          skipsForLocalAtom[i].push_back(j);      
205 +      }      
206 +    }
207 +
208 +    toposForLocalAtom.clear();
209 +    toposForLocalAtom.reserve(nLocal_);
210 +    for (int i = 0; i < nLocal_; i++) {
211 +      int iglob = AtomLocalToGlobal[i];
212 +      int nTopos = 0;
213 +      for (int j = 0; j < nLocal_; j++) {
214 +        int jglob = AtomLocalToGlobal[j];        
215 +        if (oneTwo.hasPair(iglob, jglob)) {
216 +          toposForLocalAtom[i].push_back(j);
217 +          topoDistLocal[i][nTopos] = 1;
218 +          nTopos++;
219 +        }
220 +        if (oneThree.hasPair(iglob, jglob)) {
221 +          toposForLocalAtom[i].push_back(j);
222 +          topoDistLocal[i][nTopos] = 2;
223 +          nTopos++;
224 +        }
225 +        if (oneFour.hasPair(iglob, jglob)) {
226 +          toposForLocalAtom[i].push_back(j);
227 +          topoDistLocal[i][nTopos] = 3;
228 +          nTopos++;
229 +        }
230 +      }      
231 +    }
232 +  }
233 +  
234    void ForceMatrixDecomposition::distributeData()  {
235      snap_ = sman_->getCurrentSnapshot();
236      storageLayout_ = sman_->getStorageLayout();
# Line 229 | Line 340 | namespace OpenMD {
340          snap_->atomData.torque[i] += trq_tmp[i];
341      }
342      
343 <    int nLocal = snap_->getNumberOfAtoms();
343 >    nLocal_ = snap_->getNumberOfAtoms();
344  
345      vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES,
346 <                                       vector<RealType> (nLocal, 0.0));
346 >                                       vector<RealType> (nLocal_, 0.0));
347      
348      for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
349        AtomCommRealRow->scatter(pot_row[i], pot_temp[i]);
# Line 243 | Line 354 | namespace OpenMD {
354   #endif
355    }
356  
357 +  int ForceMatrixDecomposition::getNAtomsInRow() {  
358 + #ifdef IS_MPI
359 +    return nAtomsInRow_;
360 + #else
361 +    return nLocal_;
362 + #endif
363 +  }
364 +
365 +  /**
366 +   * returns the list of atoms belonging to this group.  
367 +   */
368 +  vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
369 + #ifdef IS_MPI
370 +    return groupListRow_[cg1];
371 + #else
372 +    return groupList_[cg1];
373 + #endif
374 +  }
375 +
376 +  vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
377 + #ifdef IS_MPI
378 +    return groupListCol_[cg2];
379 + #else
380 +    return groupList_[cg2];
381 + #endif
382 +  }
383    
384    Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){
385      Vector3d d;
# Line 284 | Line 421 | namespace OpenMD {
421      snap_->wrapVector(d);
422      return d;    
423    }
424 +
425 +  RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) {
426 + #ifdef IS_MPI
427 +    return massFactorsRow[atom1];
428 + #else
429 +    return massFactorsLocal[atom1];
430 + #endif
431 +  }
432 +
433 +  RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
434 + #ifdef IS_MPI
435 +    return massFactorsCol[atom2];
436 + #else
437 +    return massFactorsLocal[atom2];
438 + #endif
439 +
440 +  }
441      
442    Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){
443      Vector3d d;
# Line 296 | Line 450 | namespace OpenMD {
450  
451      snap_->wrapVector(d);
452      return d;    
453 +  }
454 +
455 +  vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) {
456 + #ifdef IS_MPI
457 +    return skipsForRowAtom[atom1];
458 + #else
459 +    return skipsForLocalAtom[atom1];
460 + #endif
461 +  }
462 +
463 +  /**
464 +   * there are a number of reasons to skip a pair or a particle mostly
465 +   * we do this to exclude atoms who are involved in short range
466 +   * interactions (bonds, bends, torsions), but we also need to
467 +   * exclude some overcounted interactions that result from the
468 +   * parallel decomposition.
469 +   */
470 +  bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) {
471 +    int unique_id_1, unique_id_2;
472 +
473 + #ifdef IS_MPI
474 +    // in MPI, we have to look up the unique IDs for each atom
475 +    unique_id_1 = AtomRowToGlobal[atom1];
476 +    unique_id_2 = AtomColToGlobal[atom2];
477 +
478 +    // this situation should only arise in MPI simulations
479 +    if (unique_id_1 == unique_id_2) return true;
480 +    
481 +    // this prevents us from doing the pair on multiple processors
482 +    if (unique_id_1 < unique_id_2) {
483 +      if ((unique_id_1 + unique_id_2) % 2 == 0) return true;
484 +    } else {
485 +      if ((unique_id_1 + unique_id_2) % 2 == 1) return true;
486 +    }
487 + #else
488 +    // in the normal loop, the atom numbers are unique
489 +    unique_id_1 = atom1;
490 +    unique_id_2 = atom2;
491 + #endif
492 +    
493 + #ifdef IS_MPI
494 +    for (vector<int>::iterator i = skipsForRowAtom[atom1].begin();
495 +         i != skipsForRowAtom[atom1].end(); ++i) {
496 +      if ( (*i) == unique_id_2 ) return true;
497 +    }    
498 + #else
499 +    for (vector<int>::iterator i = skipsForLocalAtom[atom1].begin();
500 +         i != skipsForLocalAtom[atom1].end(); ++i) {
501 +      if ( (*i) == unique_id_2 ) return true;
502 +    }    
503 + #endif
504 +  }
505 +
506 +  int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) {
507 +    
508 + #ifdef IS_MPI
509 +    for (int i = 0; i < toposForRowAtom[atom1].size(); i++) {
510 +      if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i];
511 +    }
512 + #else
513 +    for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) {
514 +      if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i];
515 +    }
516 + #endif
517 +
518 +    // zero is default for unconnected (i.e. normal) pair interactions
519 +    return 0;
520    }
521  
522    void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){
# Line 312 | Line 533 | namespace OpenMD {
533   #else
534      snap_->atomData.force[atom2] += fg;
535   #endif
315 –
536    }
537  
538      // filling interaction blocks with pointers
539    InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {    
320 –
540      InteractionData idat;
541 +
542   #ifdef IS_MPI
543      if (storageLayout_ & DataStorage::dslAmat) {
544        idat.A1 = &(atomRowData.aMat[atom1]);
545        idat.A2 = &(atomColData.aMat[atom2]);
546      }
547 <
547 >    
548      if (storageLayout_ & DataStorage::dslElectroFrame) {
549        idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
550        idat.eFrame2 = &(atomColData.electroFrame[atom2]);
# Line 344 | Line 564 | namespace OpenMD {
564        idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
565        idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
566      }
567 +
568   #else
569      if (storageLayout_ & DataStorage::dslAmat) {
570        idat.A1 = &(snap_->atomData.aMat[atom1]);
# Line 370 | Line 591 | namespace OpenMD {
591        idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
592      }
593   #endif
594 <    
594 >    return idat;
595    }
596 +
597    InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
598 +
599      InteractionData idat;
377 –    skippedCharge1
378 –      skippedCharge2
379 –      rij
380 –      d
381 –    electroMult
382 –    sw
383 –    f
600   #ifdef IS_MPI
385 –
601      if (storageLayout_ & DataStorage::dslElectroFrame) {
602        idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
603        idat.eFrame2 = &(atomColData.electroFrame[atom2]);
# Line 391 | Line 606 | namespace OpenMD {
606        idat.t1 = &(atomRowData.torque[atom1]);
607        idat.t2 = &(atomColData.torque[atom2]);
608      }
609 <
609 >    if (storageLayout_ & DataStorage::dslForce) {
610 >      idat.t1 = &(atomRowData.force[atom1]);
611 >      idat.t2 = &(atomColData.force[atom2]);
612 >    }
613 > #else
614 >    if (storageLayout_ & DataStorage::dslElectroFrame) {
615 >      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
616 >      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
617 >    }
618 >    if (storageLayout_ & DataStorage::dslTorque) {
619 >      idat.t1 = &(snap_->atomData.torque[atom1]);
620 >      idat.t2 = &(snap_->atomData.torque[atom2]);
621 >    }
622 >    if (storageLayout_ & DataStorage::dslForce) {
623 >      idat.t1 = &(snap_->atomData.force[atom1]);
624 >      idat.t2 = &(snap_->atomData.force[atom2]);
625 >    }
626 > #endif
627      
628    }
397 –  SelfData ForceMatrixDecomposition::fillSelfData(int atom1) {
398 –  }
629  
630  
631 +
632 +
633    /*
634     * buildNeighborList
635     *
636     * first element of pair is row-indexed CutoffGroup
637     * second element of pair is column-indexed CutoffGroup
638     */
639 <  vector<pair<int, int> >  buildNeighborList() {
640 <    Vector3d dr, invWid, rs, shift;
641 <    Vector3i cc, m1v, m2s;
642 <    RealType rrNebr;
643 <    int c, j1, j2, m1, m1x, m1y, m1z, m2, n, offset;
639 >  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
640 >      
641 >    vector<pair<int, int> > neighborList;
642 > #ifdef IS_MPI
643 >    cellListRow_.clear();
644 >    cellListCol_.clear();
645 > #else
646 >    cellList_.clear();
647 > #endif
648  
649 <
650 <    vector<pair<int, int> > neighborList;  
651 <    Vector3i nCells;
652 <    Vector3d invWid, r;
653 <
654 <    rList_ = (rCut_ + skinThickness_);
419 <    rl2 = rList_ * rList_;
420 <
421 <    snap_ = sman_->getCurrentSnapshot();
649 >    // dangerous to not do error checking.
650 >    RealType rCut_;
651 >
652 >    RealType rList_ = (rCut_ + skinThickness_);
653 >    RealType rl2 = rList_ * rList_;
654 >    Snapshot* snap_ = sman_->getCurrentSnapshot();
655      Mat3x3d Hmat = snap_->getHmat();
656      Vector3d Hx = Hmat.getColumn(0);
657      Vector3d Hy = Hmat.getColumn(1);
658      Vector3d Hz = Hmat.getColumn(2);
659  
660 <    nCells.x() = (int) ( Hx.length() )/ rList_;
661 <    nCells.y() = (int) ( Hy.length() )/ rList_;
662 <    nCells.z() = (int) ( Hz.length() )/ rList_;
660 >    nCells_.x() = (int) ( Hx.length() )/ rList_;
661 >    nCells_.y() = (int) ( Hy.length() )/ rList_;
662 >    nCells_.z() = (int) ( Hz.length() )/ rList_;
663  
664 <    for (i = 0; i < nGroupsInRow; i++) {
664 >    Mat3x3d invHmat = snap_->getInvHmat();
665 >    Vector3d rs, scaled, dr;
666 >    Vector3i whichCell;
667 >    int cellIndex;
668 >
669 > #ifdef IS_MPI
670 >    for (int i = 0; i < nGroupsInRow_; i++) {
671        rs = cgRowData.position[i];
672 <      snap_->scaleVector(rs);    
672 >      // scaled positions relative to the box vectors
673 >      scaled = invHmat * rs;
674 >      // wrap the vector back into the unit box by subtracting integer box
675 >      // numbers
676 >      for (int j = 0; j < 3; j++)
677 >        scaled[j] -= roundMe(scaled[j]);
678 >    
679 >      // find xyz-indices of cell that cutoffGroup is in.
680 >      whichCell.x() = nCells_.x() * scaled.x();
681 >      whichCell.y() = nCells_.y() * scaled.y();
682 >      whichCell.z() = nCells_.z() * scaled.z();
683 >
684 >      // find single index of this cell:
685 >      cellIndex = Vlinear(whichCell, nCells_);
686 >      // add this cutoff group to the list of groups in this cell;
687 >      cellListRow_[cellIndex].push_back(i);
688      }
435 –    
689  
690 <    VDiv (invWid, cells, region);
691 <    for (n = nMol; n < nMol + cells.componentProduct(); n ++) cellList[n] = -1;
692 <    for (n = 0; n < nMol; n ++) {
693 <      VSAdd (rs, mol[n].r, 0.5, region);
694 <      VMul (cc, rs, invWid);
695 <      c = VLinear (cc, cells) + nMol;
696 <      cellList[n] = cellList[c];
697 <      cellList[c] = n;
690 >    for (int i = 0; i < nGroupsInCol_; i++) {
691 >      rs = cgColData.position[i];
692 >      // scaled positions relative to the box vectors
693 >      scaled = invHmat * rs;
694 >      // wrap the vector back into the unit box by subtracting integer box
695 >      // numbers
696 >      for (int j = 0; j < 3; j++)
697 >        scaled[j] -= roundMe(scaled[j]);
698 >
699 >      // find xyz-indices of cell that cutoffGroup is in.
700 >      whichCell.x() = nCells_.x() * scaled.x();
701 >      whichCell.y() = nCells_.y() * scaled.y();
702 >      whichCell.z() = nCells_.z() * scaled.z();
703 >
704 >      // find single index of this cell:
705 >      cellIndex = Vlinear(whichCell, nCells_);
706 >      // add this cutoff group to the list of groups in this cell;
707 >      cellListCol_[cellIndex].push_back(i);
708      }
709 <    nebrTabLen = 0;
710 <    for (m1z = 0; m1z < cells.z(); m1z++) {
711 <      for (m1y = 0; m1y < cells.y(); m1y++) {
712 <        for (m1x = 0; m1x < cells.x(); m1x++) {
709 > #else
710 >    for (int i = 0; i < nGroups_; i++) {
711 >      rs = snap_->cgData.position[i];
712 >      // scaled positions relative to the box vectors
713 >      scaled = invHmat * rs;
714 >      // wrap the vector back into the unit box by subtracting integer box
715 >      // numbers
716 >      for (int j = 0; j < 3; j++)
717 >        scaled[j] -= roundMe(scaled[j]);
718 >
719 >      // find xyz-indices of cell that cutoffGroup is in.
720 >      whichCell.x() = nCells_.x() * scaled.x();
721 >      whichCell.y() = nCells_.y() * scaled.y();
722 >      whichCell.z() = nCells_.z() * scaled.z();
723 >
724 >      // find single index of this cell:
725 >      cellIndex = Vlinear(whichCell, nCells_);
726 >      // add this cutoff group to the list of groups in this cell;
727 >      cellList_[cellIndex].push_back(i);
728 >    }
729 > #endif
730 >
731 >
732 >
733 >    for (int m1z = 0; m1z < nCells_.z(); m1z++) {
734 >      for (int m1y = 0; m1y < nCells_.y(); m1y++) {
735 >        for (int m1x = 0; m1x < nCells_.x(); m1x++) {
736            Vector3i m1v(m1x, m1y, m1z);
737 <          m1 = VLinear(m1v, cells) + nMol;
452 <          for (offset = 0; offset < nOffset_; offset++) {
453 <            m2v = m1v + cellOffsets_[offset];
454 <            shift = V3Zero();
737 >          int m1 = Vlinear(m1v, nCells_);
738  
739 <            if (m2v.x() >= cells.x) {
739 >          for (vector<Vector3i>::iterator os = cellOffsets_.begin();
740 >               os != cellOffsets_.end(); ++os) {
741 >            
742 >            Vector3i m2v = m1v + (*os);
743 >            
744 >            if (m2v.x() >= nCells_.x()) {
745                m2v.x() = 0;          
458 –              shift.x() = region.x();  
746              } else if (m2v.x() < 0) {
747 <              m2v.x() = cells.x() - 1;
461 <              shift.x() = - region.x();
747 >              m2v.x() = nCells_.x() - 1;
748              }
749 <
750 <            if (m2v.y() >= cells.y()) {
749 >            
750 >            if (m2v.y() >= nCells_.y()) {
751                m2v.y() = 0;          
466 –              shift.y() = region.y();  
752              } else if (m2v.y() < 0) {
753 <              m2v.y() = cells.y() - 1;
469 <              shift.y() = - region.y();
753 >              m2v.y() = nCells_.y() - 1;
754              }
755 +            
756 +            if (m2v.z() >= nCells_.z()) {
757 +              m2v.z() = 0;          
758 +            } else if (m2v.z() < 0) {
759 +              m2v.z() = nCells_.z() - 1;
760 +            }
761 +            
762 +            int m2 = Vlinear (m2v, nCells_);
763  
764 <            m2 = VLinear (m2v, cells) + nMol;
765 <            for (j1 = cellList[m1]; j1 >= 0; j1 = cellList[j1]) {
766 <              for (j2 = cellList[m2]; j2 >= 0; j2 = cellList[j2]) {
767 <                if (m1 != m2 || j2 < j1) {
768 <                  dr = mol[j1].r - mol[j2].r;
769 <                  VSub (dr, mol[j1].r, mol[j2].r);
770 <                  VVSub (dr, shift);
771 <                  if (VLenSq (dr) < rrNebr) {
772 <                    neighborList.push_back(make_pair(j1, j2));
764 > #ifdef IS_MPI
765 >            for (vector<int>::iterator j1 = cellListRow_[m1].begin();
766 >                 j1 != cellListRow_[m1].end(); ++j1) {
767 >              for (vector<int>::iterator j2 = cellListCol_[m2].begin();
768 >                   j2 != cellListCol_[m2].end(); ++j2) {
769 >                              
770 >                // Always do this if we're in different cells or if
771 >                // we're in the same cell and the global index of the
772 >                // j2 cutoff group is less than the j1 cutoff group
773 >
774 >                if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) {
775 >                  dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)];
776 >                  snap_->wrapVector(dr);
777 >                  if (dr.lengthSquare() < rl2) {
778 >                    neighborList.push_back(make_pair((*j1), (*j2)));
779                    }
780                  }
781                }
782              }
783 + #else
784 +            for (vector<int>::iterator j1 = cellList_[m1].begin();
785 +                 j1 != cellList_[m1].end(); ++j1) {
786 +              for (vector<int>::iterator j2 = cellList_[m2].begin();
787 +                   j2 != cellList_[m2].end(); ++j2) {
788 +                              
789 +                // Always do this if we're in different cells or if
790 +                // we're in the same cell and the global index of the
791 +                // j2 cutoff group is less than the j1 cutoff group
792 +
793 +                if (m2 != m1 || (*j2) < (*j1)) {
794 +                  dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
795 +                  snap_->wrapVector(dr);
796 +                  if (dr.lengthSquare() < rl2) {
797 +                    neighborList.push_back(make_pair((*j1), (*j2)));
798 +                  }
799 +                }
800 +              }
801 +            }
802 + #endif
803            }
804          }
805        }
806      }
489 –  }
807  
808 <  
808 >    // save the local cutoff group positions for the check that is
809 >    // done on each loop:
810 >    saved_CG_positions_.clear();
811 >    for (int i = 0; i < nGroups_; i++)
812 >      saved_CG_positions_.push_back(snap_->cgData.position[i]);
813 >
814 >    return neighborList;
815 >  }
816   } //end namespace OpenMD

Diff Legend

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