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 1569 by gezelter, Thu May 26 13:55:04 2011 UTC

# Line 54 | Line 54 | namespace OpenMD {
54    void ForceMatrixDecomposition::distributeInitialData() {
55      snap_ = sman_->getCurrentSnapshot();
56      storageLayout_ = sman_->getStorageLayout();
57 < #ifdef IS_MPI    
58 <    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);
57 >    nLocal_ = snap_->getNumberOfAtoms();
58 >    nGroups_ = snap_->getNumberOfCutoffGroups();
59  
60 <    AtomCommIntColumn = new Communicator<Column,int>(nLocal);
61 <    AtomCommRealColumn = new Communicator<Column,RealType>(nLocal);
62 <    AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal);
63 <    AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal);
60 >    // gather the information for atomtype IDs (atids):
61 >    vector<int> identsLocal = info_->getIdentArray();
62 >    AtomLocalToGlobal = info_->getGlobalAtomIndices();
63 >    cgLocalToGlobal = info_->getGlobalGroupIndices();
64 >    vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
65 >    vector<RealType> massFactorsLocal = info_->getMassFactors();
66 >    vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
67  
68 <    cgCommIntRow = new Communicator<Row,int>(nGroups);
69 <    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups);
70 <    cgCommIntColumn = new Communicator<Column,int>(nGroups);
71 <    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups);
68 > #ifdef IS_MPI
69 >
70 >    AtomCommIntRow = new Communicator<Row,int>(nLocal_);
71 >    AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
72 >    AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
73 >    AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
74  
75 <    int nAtomsInRow = AtomCommIntRow->getSize();
76 <    int nAtomsInCol = AtomCommIntColumn->getSize();
77 <    int nGroupsInRow = cgCommIntRow->getSize();
78 <    int nGroupsInCol = cgCommIntColumn->getSize();
75 >    AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
76 >    AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
77 >    AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
78 >    AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
79  
80 +    cgCommIntRow = new Communicator<Row,int>(nGroups_);
81 +    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
82 +    cgCommIntColumn = new Communicator<Column,int>(nGroups_);
83 +    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_);
84 +
85 +    nAtomsInRow_ = AtomCommIntRow->getSize();
86 +    nAtomsInCol_ = AtomCommIntColumn->getSize();
87 +    nGroupsInRow_ = cgCommIntRow->getSize();
88 +    nGroupsInCol_ = cgCommIntColumn->getSize();
89 +
90      // Modify the data storage objects with the correct layouts and sizes:
91 <    atomRowData.resize(nAtomsInRow);
91 >    atomRowData.resize(nAtomsInRow_);
92      atomRowData.setStorageLayout(storageLayout_);
93 <    atomColData.resize(nAtomsInCol);
93 >    atomColData.resize(nAtomsInCol_);
94      atomColData.setStorageLayout(storageLayout_);
95 <    cgRowData.resize(nGroupsInRow);
95 >    cgRowData.resize(nGroupsInRow_);
96      cgRowData.setStorageLayout(DataStorage::dslPosition);
97 <    cgColData.resize(nGroupsInCol);
97 >    cgColData.resize(nGroupsInCol_);
98      cgColData.setStorageLayout(DataStorage::dslPosition);
99      
100      vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
101 <                                      vector<RealType> (nAtomsInRow, 0.0));
101 >                                      vector<RealType> (nAtomsInRow_, 0.0));
102      vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
103 <                                      vector<RealType> (nAtomsInCol, 0.0));
95 <
96 <
97 <    vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
103 >                                      vector<RealType> (nAtomsInCol_, 0.0));
104      
105 <    // gather the information for atomtype IDs (atids):
106 <    vector<int> identsLocal = info_->getIdentArray();
101 <    identsRow.reserve(nAtomsInRow);
102 <    identsCol.reserve(nAtomsInCol);
105 >    identsRow.reserve(nAtomsInRow_);
106 >    identsCol.reserve(nAtomsInCol_);
107      
108      AtomCommIntRow->gather(identsLocal, identsRow);
109      AtomCommIntColumn->gather(identsLocal, identsCol);
110      
107 –    AtomLocalToGlobal = info_->getGlobalAtomIndices();
111      AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
112      AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
113      
111 –    cgLocalToGlobal = info_->getGlobalGroupIndices();
114      cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
115      cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
116 +
117 +    AtomCommRealRow->gather(massFactorsLocal, massFactorsRow);
118 +    AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol);
119 +
120 +    groupListRow_.clear();
121 +    groupListRow_.reserve(nGroupsInRow_);
122 +    for (int i = 0; i < nGroupsInRow_; i++) {
123 +      int gid = cgRowToGlobal[i];
124 +      for (int j = 0; j < nAtomsInRow_; j++) {
125 +        int aid = AtomRowToGlobal[j];
126 +        if (globalGroupMembership[aid] == gid)
127 +          groupListRow_[i].push_back(j);
128 +      }      
129 +    }
130 +
131 +    groupListCol_.clear();
132 +    groupListCol_.reserve(nGroupsInCol_);
133 +    for (int i = 0; i < nGroupsInCol_; i++) {
134 +      int gid = cgColToGlobal[i];
135 +      for (int j = 0; j < nAtomsInCol_; j++) {
136 +        int aid = AtomColToGlobal[j];
137 +        if (globalGroupMembership[aid] == gid)
138 +          groupListCol_[i].push_back(j);
139 +      }      
140 +    }
141 +
142 + #endif
143 +
144 +    groupList_.clear();
145 +    groupList_.reserve(nGroups_);
146 +    for (int i = 0; i < nGroups_; i++) {
147 +      int gid = cgLocalToGlobal[i];
148 +      for (int j = 0; j < nLocal_; j++) {
149 +        int aid = AtomLocalToGlobal[j];
150 +        if (globalGroupMembership[aid] == gid)
151 +          groupList_[i].push_back(j);
152 +      }      
153 +    }
154  
155 +  
156      // still need:
157      // topoDist
158      // exclude
159 < #endif
159 >
160    }
161      
162  
# Line 229 | Line 270 | namespace OpenMD {
270          snap_->atomData.torque[i] += trq_tmp[i];
271      }
272      
273 <    int nLocal = snap_->getNumberOfAtoms();
273 >    nLocal_ = snap_->getNumberOfAtoms();
274  
275      vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES,
276 <                                       vector<RealType> (nLocal, 0.0));
276 >                                       vector<RealType> (nLocal_, 0.0));
277      
278      for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
279        AtomCommRealRow->scatter(pot_row[i], pot_temp[i]);
# Line 243 | Line 284 | namespace OpenMD {
284   #endif
285    }
286  
287 +  /**
288 +   * returns the list of atoms belonging to this group.  
289 +   */
290 +  vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
291 + #ifdef IS_MPI
292 +    return groupListRow_[cg1];
293 + #else
294 +    return groupList_[cg1];
295 + #endif
296 +  }
297 +
298 +  vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
299 + #ifdef IS_MPI
300 +    return groupListCol_[cg2];
301 + #else
302 +    return groupList_[cg2];
303 + #endif
304 +  }
305    
306    Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){
307      Vector3d d;
# Line 284 | Line 343 | namespace OpenMD {
343      snap_->wrapVector(d);
344      return d;    
345    }
346 +
347 +  RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) {
348 + #ifdef IS_MPI
349 +    return massFactorsRow[atom1];
350 + #else
351 +    return massFactorsLocal[atom1];
352 + #endif
353 +  }
354 +
355 +  RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
356 + #ifdef IS_MPI
357 +    return massFactorsCol[atom2];
358 + #else
359 +    return massFactorsLocal[atom2];
360 + #endif
361 +
362 +  }
363      
364    Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){
365      Vector3d d;
# Line 312 | Line 388 | namespace OpenMD {
388   #else
389      snap_->atomData.force[atom2] += fg;
390   #endif
315 –
391    }
392  
393      // filling interaction blocks with pointers
394    InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {    
320 –
395      InteractionData idat;
396 +
397   #ifdef IS_MPI
398      if (storageLayout_ & DataStorage::dslAmat) {
399        idat.A1 = &(atomRowData.aMat[atom1]);
400        idat.A2 = &(atomColData.aMat[atom2]);
401      }
402 <
402 >    
403      if (storageLayout_ & DataStorage::dslElectroFrame) {
404        idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
405        idat.eFrame2 = &(atomColData.electroFrame[atom2]);
# Line 370 | Line 445 | namespace OpenMD {
445        idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
446      }
447   #endif
448 <    
448 >    return idat;
449    }
450 +
451    InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
452 +
453      InteractionData idat;
377 –    skippedCharge1
378 –      skippedCharge2
379 –      rij
380 –      d
381 –    electroMult
382 –    sw
383 –    f
454   #ifdef IS_MPI
385 –
455      if (storageLayout_ & DataStorage::dslElectroFrame) {
456        idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
457        idat.eFrame2 = &(atomColData.electroFrame[atom2]);
# Line 391 | Line 460 | namespace OpenMD {
460        idat.t1 = &(atomRowData.torque[atom1]);
461        idat.t2 = &(atomColData.torque[atom2]);
462      }
463 <
463 >    if (storageLayout_ & DataStorage::dslForce) {
464 >      idat.t1 = &(atomRowData.force[atom1]);
465 >      idat.t2 = &(atomColData.force[atom2]);
466 >    }
467 > #else
468 >    if (storageLayout_ & DataStorage::dslElectroFrame) {
469 >      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
470 >      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
471 >    }
472 >    if (storageLayout_ & DataStorage::dslTorque) {
473 >      idat.t1 = &(snap_->atomData.torque[atom1]);
474 >      idat.t2 = &(snap_->atomData.torque[atom2]);
475 >    }
476 >    if (storageLayout_ & DataStorage::dslForce) {
477 >      idat.t1 = &(snap_->atomData.force[atom1]);
478 >      idat.t2 = &(snap_->atomData.force[atom2]);
479 >    }
480 > #endif
481      
482    }
397 –  SelfData ForceMatrixDecomposition::fillSelfData(int atom1) {
398 –  }
483  
484  
485 +
486 +
487    /*
488     * buildNeighborList
489     *
490     * first element of pair is row-indexed CutoffGroup
491     * second element of pair is column-indexed CutoffGroup
492     */
493 <  vector<pair<int, int> >  buildNeighborList() {
494 <    Vector3d dr, invWid, rs, shift;
495 <    Vector3i cc, m1v, m2s;
496 <    RealType rrNebr;
497 <    int c, j1, j2, m1, m1x, m1y, m1z, m2, n, offset;
498 <
499 <
500 <    vector<pair<int, int> > neighborList;  
501 <    Vector3i nCells;
416 <    Vector3d invWid, r;
417 <
418 <    rList_ = (rCut_ + skinThickness_);
419 <    rl2 = rList_ * rList_;
493 >  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
494 >      
495 >    vector<pair<int, int> > neighborList;
496 > #ifdef IS_MPI
497 >    cellListRow_.clear();
498 >    cellListCol_.clear();
499 > #else
500 >    cellList_.clear();
501 > #endif
502  
503 <    snap_ = sman_->getCurrentSnapshot();
503 >    // dangerous to not do error checking.
504 >    RealType rCut_;
505 >
506 >    RealType rList_ = (rCut_ + skinThickness_);
507 >    RealType rl2 = rList_ * rList_;
508 >    Snapshot* snap_ = sman_->getCurrentSnapshot();
509      Mat3x3d Hmat = snap_->getHmat();
510      Vector3d Hx = Hmat.getColumn(0);
511      Vector3d Hy = Hmat.getColumn(1);
512      Vector3d Hz = Hmat.getColumn(2);
513  
514 <    nCells.x() = (int) ( Hx.length() )/ rList_;
515 <    nCells.y() = (int) ( Hy.length() )/ rList_;
516 <    nCells.z() = (int) ( Hz.length() )/ rList_;
514 >    nCells_.x() = (int) ( Hx.length() )/ rList_;
515 >    nCells_.y() = (int) ( Hy.length() )/ rList_;
516 >    nCells_.z() = (int) ( Hz.length() )/ rList_;
517  
518 <    for (i = 0; i < nGroupsInRow; i++) {
518 >    Mat3x3d invHmat = snap_->getInvHmat();
519 >    Vector3d rs, scaled, dr;
520 >    Vector3i whichCell;
521 >    int cellIndex;
522 >
523 > #ifdef IS_MPI
524 >    for (int i = 0; i < nGroupsInRow_; i++) {
525        rs = cgRowData.position[i];
526 <      snap_->scaleVector(rs);    
526 >      // scaled positions relative to the box vectors
527 >      scaled = invHmat * rs;
528 >      // wrap the vector back into the unit box by subtracting integer box
529 >      // numbers
530 >      for (int j = 0; j < 3; j++)
531 >        scaled[j] -= roundMe(scaled[j]);
532 >    
533 >      // find xyz-indices of cell that cutoffGroup is in.
534 >      whichCell.x() = nCells_.x() * scaled.x();
535 >      whichCell.y() = nCells_.y() * scaled.y();
536 >      whichCell.z() = nCells_.z() * scaled.z();
537 >
538 >      // find single index of this cell:
539 >      cellIndex = Vlinear(whichCell, nCells_);
540 >      // add this cutoff group to the list of groups in this cell;
541 >      cellListRow_[cellIndex].push_back(i);
542      }
435 –    
543  
544 <    VDiv (invWid, cells, region);
545 <    for (n = nMol; n < nMol + cells.componentProduct(); n ++) cellList[n] = -1;
546 <    for (n = 0; n < nMol; n ++) {
547 <      VSAdd (rs, mol[n].r, 0.5, region);
548 <      VMul (cc, rs, invWid);
549 <      c = VLinear (cc, cells) + nMol;
550 <      cellList[n] = cellList[c];
551 <      cellList[c] = n;
544 >    for (int i = 0; i < nGroupsInCol_; i++) {
545 >      rs = cgColData.position[i];
546 >      // scaled positions relative to the box vectors
547 >      scaled = invHmat * rs;
548 >      // wrap the vector back into the unit box by subtracting integer box
549 >      // numbers
550 >      for (int j = 0; j < 3; j++)
551 >        scaled[j] -= roundMe(scaled[j]);
552 >
553 >      // find xyz-indices of cell that cutoffGroup is in.
554 >      whichCell.x() = nCells_.x() * scaled.x();
555 >      whichCell.y() = nCells_.y() * scaled.y();
556 >      whichCell.z() = nCells_.z() * scaled.z();
557 >
558 >      // find single index of this cell:
559 >      cellIndex = Vlinear(whichCell, nCells_);
560 >      // add this cutoff group to the list of groups in this cell;
561 >      cellListCol_[cellIndex].push_back(i);
562      }
563 <    nebrTabLen = 0;
564 <    for (m1z = 0; m1z < cells.z(); m1z++) {
565 <      for (m1y = 0; m1y < cells.y(); m1y++) {
566 <        for (m1x = 0; m1x < cells.x(); m1x++) {
563 > #else
564 >    for (int i = 0; i < nGroups_; i++) {
565 >      rs = snap_->cgData.position[i];
566 >      // scaled positions relative to the box vectors
567 >      scaled = invHmat * rs;
568 >      // wrap the vector back into the unit box by subtracting integer box
569 >      // numbers
570 >      for (int j = 0; j < 3; j++)
571 >        scaled[j] -= roundMe(scaled[j]);
572 >
573 >      // find xyz-indices of cell that cutoffGroup is in.
574 >      whichCell.x() = nCells_.x() * scaled.x();
575 >      whichCell.y() = nCells_.y() * scaled.y();
576 >      whichCell.z() = nCells_.z() * scaled.z();
577 >
578 >      // find single index of this cell:
579 >      cellIndex = Vlinear(whichCell, nCells_);
580 >      // add this cutoff group to the list of groups in this cell;
581 >      cellList_[cellIndex].push_back(i);
582 >    }
583 > #endif
584 >
585 >
586 >
587 >    for (int m1z = 0; m1z < nCells_.z(); m1z++) {
588 >      for (int m1y = 0; m1y < nCells_.y(); m1y++) {
589 >        for (int m1x = 0; m1x < nCells_.x(); m1x++) {
590            Vector3i m1v(m1x, m1y, m1z);
591 <          m1 = VLinear(m1v, cells) + nMol;
452 <          for (offset = 0; offset < nOffset_; offset++) {
453 <            m2v = m1v + cellOffsets_[offset];
454 <            shift = V3Zero();
591 >          int m1 = Vlinear(m1v, nCells_);
592  
593 <            if (m2v.x() >= cells.x) {
593 >          for (vector<Vector3i>::iterator os = cellOffsets_.begin();
594 >               os != cellOffsets_.end(); ++os) {
595 >            
596 >            Vector3i m2v = m1v + (*os);
597 >            
598 >            if (m2v.x() >= nCells_.x()) {
599                m2v.x() = 0;          
458 –              shift.x() = region.x();  
600              } else if (m2v.x() < 0) {
601 <              m2v.x() = cells.x() - 1;
461 <              shift.x() = - region.x();
601 >              m2v.x() = nCells_.x() - 1;
602              }
603 <
604 <            if (m2v.y() >= cells.y()) {
603 >            
604 >            if (m2v.y() >= nCells_.y()) {
605                m2v.y() = 0;          
466 –              shift.y() = region.y();  
606              } else if (m2v.y() < 0) {
607 <              m2v.y() = cells.y() - 1;
469 <              shift.y() = - region.y();
607 >              m2v.y() = nCells_.y() - 1;
608              }
609 +            
610 +            if (m2v.z() >= nCells_.z()) {
611 +              m2v.z() = 0;          
612 +            } else if (m2v.z() < 0) {
613 +              m2v.z() = nCells_.z() - 1;
614 +            }
615 +            
616 +            int m2 = Vlinear (m2v, nCells_);
617  
618 <            m2 = VLinear (m2v, cells) + nMol;
619 <            for (j1 = cellList[m1]; j1 >= 0; j1 = cellList[j1]) {
620 <              for (j2 = cellList[m2]; j2 >= 0; j2 = cellList[j2]) {
621 <                if (m1 != m2 || j2 < j1) {
622 <                  dr = mol[j1].r - mol[j2].r;
623 <                  VSub (dr, mol[j1].r, mol[j2].r);
624 <                  VVSub (dr, shift);
625 <                  if (VLenSq (dr) < rrNebr) {
626 <                    neighborList.push_back(make_pair(j1, j2));
618 > #ifdef IS_MPI
619 >            for (vector<int>::iterator j1 = cellListRow_[m1].begin();
620 >                 j1 != cellListRow_[m1].end(); ++j1) {
621 >              for (vector<int>::iterator j2 = cellListCol_[m2].begin();
622 >                   j2 != cellListCol_[m2].end(); ++j2) {
623 >                              
624 >                // Always do this if we're in different cells or if
625 >                // we're in the same cell and the global index of the
626 >                // j2 cutoff group is less than the j1 cutoff group
627 >
628 >                if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) {
629 >                  dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)];
630 >                  snap_->wrapVector(dr);
631 >                  if (dr.lengthSquare() < rl2) {
632 >                    neighborList.push_back(make_pair((*j1), (*j2)));
633                    }
634                  }
635                }
636              }
637 + #else
638 +            for (vector<int>::iterator j1 = cellList_[m1].begin();
639 +                 j1 != cellList_[m1].end(); ++j1) {
640 +              for (vector<int>::iterator j2 = cellList_[m2].begin();
641 +                   j2 != cellList_[m2].end(); ++j2) {
642 +                              
643 +                // Always do this if we're in different cells or if
644 +                // we're in the same cell and the global index of the
645 +                // j2 cutoff group is less than the j1 cutoff group
646 +
647 +                if (m2 != m1 || (*j2) < (*j1)) {
648 +                  dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
649 +                  snap_->wrapVector(dr);
650 +                  if (dr.lengthSquare() < rl2) {
651 +                    neighborList.push_back(make_pair((*j1), (*j2)));
652 +                  }
653 +                }
654 +              }
655 +            }
656 + #endif
657            }
658          }
659        }
660      }
489 –  }
661  
662 <  
662 >    // save the local cutoff group positions for the check that is
663 >    // done on each loop:
664 >    saved_CG_positions_.clear();
665 >    for (int i = 0; i < nGroups_; i++)
666 >      saved_CG_positions_.push_back(snap_->cgData.position[i]);
667 >
668 >    return neighborList;
669 >  }
670   } //end namespace OpenMD

Diff Legend

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