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 1554 by gezelter, Sat Apr 30 02:54:02 2011 UTC vs.
Revision 1568 by gezelter, Wed May 25 16:20:37 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 > #ifdef IS_MPI
61 >
62 >    AtomCommIntRow = new Communicator<Row,int>(nLocal_);
63 >    AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
64 >    AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
65 >    AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
66  
67 <    cgCommIntRow = new Communicator<Row,int>(nGroups);
68 <    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups);
69 <    cgCommIntColumn = new Communicator<Column,int>(nGroups);
70 <    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups);
67 >    AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
68 >    AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
69 >    AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
70 >    AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
71  
72 <    int nAtomsInRow = AtomCommIntRow->getSize();
73 <    int nAtomsInCol = AtomCommIntColumn->getSize();
74 <    int nGroupsInRow = cgCommIntRow->getSize();
75 <    int nGroupsInCol = cgCommIntColumn->getSize();
72 >    cgCommIntRow = new Communicator<Row,int>(nGroups_);
73 >    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
74 >    cgCommIntColumn = new Communicator<Column,int>(nGroups_);
75 >    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_);
76  
77 +    nAtomsInRow_ = AtomCommIntRow->getSize();
78 +    nAtomsInCol_ = AtomCommIntColumn->getSize();
79 +    nGroupsInRow_ = cgCommIntRow->getSize();
80 +    nGroupsInCol_ = cgCommIntColumn->getSize();
81 +
82      // Modify the data storage objects with the correct layouts and sizes:
83 <    atomRowData.resize(nAtomsInRow);
83 >    atomRowData.resize(nAtomsInRow_);
84      atomRowData.setStorageLayout(storageLayout_);
85 <    atomColData.resize(nAtomsInCol);
85 >    atomColData.resize(nAtomsInCol_);
86      atomColData.setStorageLayout(storageLayout_);
87 <    cgRowData.resize(nGroupsInRow);
87 >    cgRowData.resize(nGroupsInRow_);
88      cgRowData.setStorageLayout(DataStorage::dslPosition);
89 <    cgColData.resize(nGroupsInCol);
89 >    cgColData.resize(nGroupsInCol_);
90      cgColData.setStorageLayout(DataStorage::dslPosition);
91      
92      vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
93 <                                      vector<RealType> (nAtomsInRow, 0.0));
93 >                                      vector<RealType> (nAtomsInRow_, 0.0));
94      vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
95 <                                      vector<RealType> (nAtomsInCol, 0.0));
95 >                                      vector<RealType> (nAtomsInCol_, 0.0));
96  
97  
98      vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
99      
100      // gather the information for atomtype IDs (atids):
101      vector<int> identsLocal = info_->getIdentArray();
102 <    identsRow.reserve(nAtomsInRow);
103 <    identsCol.reserve(nAtomsInCol);
102 >    identsRow.reserve(nAtomsInRow_);
103 >    identsCol.reserve(nAtomsInCol_);
104      
105      AtomCommIntRow->gather(identsLocal, identsRow);
106      AtomCommIntColumn->gather(identsLocal, identsCol);
# Line 229 | Line 230 | namespace OpenMD {
230          snap_->atomData.torque[i] += trq_tmp[i];
231      }
232      
233 <    int nLocal = snap_->getNumberOfAtoms();
233 >    nLocal_ = snap_->getNumberOfAtoms();
234  
235      vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES,
236 <                                       vector<RealType> (nLocal, 0.0));
236 >                                       vector<RealType> (nLocal_, 0.0));
237      
238      for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
239        AtomCommRealRow->scatter(pot_row[i], pot_temp[i]);
# Line 312 | Line 313 | namespace OpenMD {
313   #else
314      snap_->atomData.force[atom2] += fg;
315   #endif
315 –
316    }
317  
318      // filling interaction blocks with pointers
319    InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {    
320 –
320      InteractionData idat;
321 +
322   #ifdef IS_MPI
323      if (storageLayout_ & DataStorage::dslAmat) {
324        idat.A1 = &(atomRowData.aMat[atom1]);
325        idat.A2 = &(atomColData.aMat[atom2]);
326      }
327 <
327 >    
328      if (storageLayout_ & DataStorage::dslElectroFrame) {
329        idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
330        idat.eFrame2 = &(atomColData.electroFrame[atom2]);
# Line 343 | Line 343 | namespace OpenMD {
343      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
344        idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
345        idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
346 +    }
347 + #else
348 +    if (storageLayout_ & DataStorage::dslAmat) {
349 +      idat.A1 = &(snap_->atomData.aMat[atom1]);
350 +      idat.A2 = &(snap_->atomData.aMat[atom2]);
351 +    }
352 +
353 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
354 +      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
355 +      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
356 +    }
357 +
358 +    if (storageLayout_ & DataStorage::dslTorque) {
359 +      idat.t1 = &(snap_->atomData.torque[atom1]);
360 +      idat.t2 = &(snap_->atomData.torque[atom2]);
361      }
362 +
363 +    if (storageLayout_ & DataStorage::dslDensity) {
364 +      idat.rho1 = &(snap_->atomData.density[atom1]);
365 +      idat.rho2 = &(snap_->atomData.density[atom2]);
366 +    }
367 +
368 +    if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
369 +      idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
370 +      idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
371 +    }
372   #endif
373 <    
373 >    return idat;
374    }
375 +
376    InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
377 +
378 +    InteractionData idat;
379 + #ifdef IS_MPI
380 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
381 +      idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
382 +      idat.eFrame2 = &(atomColData.electroFrame[atom2]);
383 +    }
384 +    if (storageLayout_ & DataStorage::dslTorque) {
385 +      idat.t1 = &(atomRowData.torque[atom1]);
386 +      idat.t2 = &(atomColData.torque[atom2]);
387 +    }
388 +    if (storageLayout_ & DataStorage::dslForce) {
389 +      idat.t1 = &(atomRowData.force[atom1]);
390 +      idat.t2 = &(atomColData.force[atom2]);
391 +    }
392 + #else
393 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
394 +      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
395 +      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
396 +    }
397 +    if (storageLayout_ & DataStorage::dslTorque) {
398 +      idat.t1 = &(snap_->atomData.torque[atom1]);
399 +      idat.t2 = &(snap_->atomData.torque[atom2]);
400 +    }
401 +    if (storageLayout_ & DataStorage::dslForce) {
402 +      idat.t1 = &(snap_->atomData.force[atom1]);
403 +      idat.t2 = &(snap_->atomData.force[atom2]);
404 +    }
405 + #endif
406 +    
407    }
352 –  SelfData ForceMatrixDecomposition::fillSelfData(int atom1) {
353 –  }
408  
409 <  
409 >
410 >
411 >
412 >  /*
413 >   * buildNeighborList
414 >   *
415 >   * first element of pair is row-indexed CutoffGroup
416 >   * second element of pair is column-indexed CutoffGroup
417 >   */
418 >  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
419 >      
420 >    vector<pair<int, int> > neighborList;
421 > #ifdef IS_MPI
422 >    cellListRow_.clear();
423 >    cellListCol_.clear();
424 > #else
425 >    cellList_.clear();
426 > #endif
427 >
428 >    // dangerous to not do error checking.
429 >    RealType rCut_;
430 >
431 >    RealType rList_ = (rCut_ + skinThickness_);
432 >    RealType rl2 = rList_ * rList_;
433 >    Snapshot* snap_ = sman_->getCurrentSnapshot();
434 >    Mat3x3d Hmat = snap_->getHmat();
435 >    Vector3d Hx = Hmat.getColumn(0);
436 >    Vector3d Hy = Hmat.getColumn(1);
437 >    Vector3d Hz = Hmat.getColumn(2);
438 >
439 >    nCells_.x() = (int) ( Hx.length() )/ rList_;
440 >    nCells_.y() = (int) ( Hy.length() )/ rList_;
441 >    nCells_.z() = (int) ( Hz.length() )/ rList_;
442 >
443 >    Mat3x3d invHmat = snap_->getInvHmat();
444 >    Vector3d rs, scaled, dr;
445 >    Vector3i whichCell;
446 >    int cellIndex;
447 >
448 > #ifdef IS_MPI
449 >    for (int i = 0; i < nGroupsInRow_; i++) {
450 >      rs = cgRowData.position[i];
451 >      // scaled positions relative to the box vectors
452 >      scaled = invHmat * rs;
453 >      // wrap the vector back into the unit box by subtracting integer box
454 >      // numbers
455 >      for (int j = 0; j < 3; j++)
456 >        scaled[j] -= roundMe(scaled[j]);
457 >    
458 >      // find xyz-indices of cell that cutoffGroup is in.
459 >      whichCell.x() = nCells_.x() * scaled.x();
460 >      whichCell.y() = nCells_.y() * scaled.y();
461 >      whichCell.z() = nCells_.z() * scaled.z();
462 >
463 >      // find single index of this cell:
464 >      cellIndex = Vlinear(whichCell, nCells_);
465 >      // add this cutoff group to the list of groups in this cell;
466 >      cellListRow_[cellIndex].push_back(i);
467 >    }
468 >
469 >    for (int i = 0; i < nGroupsInCol_; i++) {
470 >      rs = cgColData.position[i];
471 >      // scaled positions relative to the box vectors
472 >      scaled = invHmat * rs;
473 >      // wrap the vector back into the unit box by subtracting integer box
474 >      // numbers
475 >      for (int j = 0; j < 3; j++)
476 >        scaled[j] -= roundMe(scaled[j]);
477 >
478 >      // find xyz-indices of cell that cutoffGroup is in.
479 >      whichCell.x() = nCells_.x() * scaled.x();
480 >      whichCell.y() = nCells_.y() * scaled.y();
481 >      whichCell.z() = nCells_.z() * scaled.z();
482 >
483 >      // find single index of this cell:
484 >      cellIndex = Vlinear(whichCell, nCells_);
485 >      // add this cutoff group to the list of groups in this cell;
486 >      cellListCol_[cellIndex].push_back(i);
487 >    }
488 > #else
489 >    for (int i = 0; i < nGroups_; i++) {
490 >      rs = snap_->cgData.position[i];
491 >      // scaled positions relative to the box vectors
492 >      scaled = invHmat * rs;
493 >      // wrap the vector back into the unit box by subtracting integer box
494 >      // numbers
495 >      for (int j = 0; j < 3; j++)
496 >        scaled[j] -= roundMe(scaled[j]);
497 >
498 >      // find xyz-indices of cell that cutoffGroup is in.
499 >      whichCell.x() = nCells_.x() * scaled.x();
500 >      whichCell.y() = nCells_.y() * scaled.y();
501 >      whichCell.z() = nCells_.z() * scaled.z();
502 >
503 >      // find single index of this cell:
504 >      cellIndex = Vlinear(whichCell, nCells_);
505 >      // add this cutoff group to the list of groups in this cell;
506 >      cellList_[cellIndex].push_back(i);
507 >    }
508 > #endif
509 >
510 >
511 >
512 >    for (int m1z = 0; m1z < nCells_.z(); m1z++) {
513 >      for (int m1y = 0; m1y < nCells_.y(); m1y++) {
514 >        for (int m1x = 0; m1x < nCells_.x(); m1x++) {
515 >          Vector3i m1v(m1x, m1y, m1z);
516 >          int m1 = Vlinear(m1v, nCells_);
517 >
518 >          for (vector<Vector3i>::iterator os = cellOffsets_.begin();
519 >               os != cellOffsets_.end(); ++os) {
520 >            
521 >            Vector3i m2v = m1v + (*os);
522 >            
523 >            if (m2v.x() >= nCells_.x()) {
524 >              m2v.x() = 0;          
525 >            } else if (m2v.x() < 0) {
526 >              m2v.x() = nCells_.x() - 1;
527 >            }
528 >            
529 >            if (m2v.y() >= nCells_.y()) {
530 >              m2v.y() = 0;          
531 >            } else if (m2v.y() < 0) {
532 >              m2v.y() = nCells_.y() - 1;
533 >            }
534 >            
535 >            if (m2v.z() >= nCells_.z()) {
536 >              m2v.z() = 0;          
537 >            } else if (m2v.z() < 0) {
538 >              m2v.z() = nCells_.z() - 1;
539 >            }
540 >            
541 >            int m2 = Vlinear (m2v, nCells_);
542 >
543 > #ifdef IS_MPI
544 >            for (vector<int>::iterator j1 = cellListRow_[m1].begin();
545 >                 j1 != cellListRow_[m1].end(); ++j1) {
546 >              for (vector<int>::iterator j2 = cellListCol_[m2].begin();
547 >                   j2 != cellListCol_[m2].end(); ++j2) {
548 >                              
549 >                // Always do this if we're in different cells or if
550 >                // we're in the same cell and the global index of the
551 >                // j2 cutoff group is less than the j1 cutoff group
552 >
553 >                if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) {
554 >                  dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)];
555 >                  snap_->wrapVector(dr);
556 >                  if (dr.lengthSquare() < rl2) {
557 >                    neighborList.push_back(make_pair((*j1), (*j2)));
558 >                  }
559 >                }
560 >              }
561 >            }
562 > #else
563 >            for (vector<int>::iterator j1 = cellList_[m1].begin();
564 >                 j1 != cellList_[m1].end(); ++j1) {
565 >              for (vector<int>::iterator j2 = cellList_[m2].begin();
566 >                   j2 != cellList_[m2].end(); ++j2) {
567 >                              
568 >                // Always do this if we're in different cells or if
569 >                // we're in the same cell and the global index of the
570 >                // j2 cutoff group is less than the j1 cutoff group
571 >
572 >                if (m2 != m1 || (*j2) < (*j1)) {
573 >                  dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
574 >                  snap_->wrapVector(dr);
575 >                  if (dr.lengthSquare() < rl2) {
576 >                    neighborList.push_back(make_pair((*j1), (*j2)));
577 >                  }
578 >                }
579 >              }
580 >            }
581 > #endif
582 >          }
583 >        }
584 >      }
585 >    }
586 >
587 >    // save the local cutoff group positions for the check that is
588 >    // done on each loop:
589 >    saved_CG_positions_.clear();
590 >    for (int i = 0; i < nGroups_; i++)
591 >      saved_CG_positions_.push_back(snap_->cgData.position[i]);
592 >
593 >    return neighborList;
594 >  }
595   } //end namespace OpenMD

Diff Legend

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