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

Comparing trunk/src/parallel/ForceMatrixDecomposition.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 1896 by gezelter, Tue Jul 2 20:02:31 2013 UTC

# Line 308 | Line 308 | namespace OpenMD {
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;
317      int atid;
# Line 419 | Line 423 | namespace OpenMD {
423      
424      RealType tradRcut = groupMax;
425  
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_) {
# Line 442 | Line 455 | namespace OpenMD {
455            break;
456          }
457  
458 <        pair<int,int> key = make_pair(i,j);
446 <        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 463 | Line 479 | namespace OpenMD {
479      }
480    }
481  
482 <  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 472 | 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 (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  
# Line 939 | Line 959 | namespace OpenMD {
959  
960  
961  
962 <  int ForceMatrixDecomposition::getNAtomsInRow() {  
962 >  int& ForceMatrixDecomposition::getNAtomsInRow() {  
963   #ifdef IS_MPI
964      return nAtomsInRow_;
965   #else
# Line 950 | 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 958 | 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 981 | Line 1001 | namespace OpenMD {
1001      return d;    
1002    }
1003  
1004 <  Vector3d ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){
1004 >  Vector3d& ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){
1005   #ifdef IS_MPI
1006      return cgColData.velocity[cg2];
1007   #else
# Line 989 | Line 1009 | namespace OpenMD {
1009   #endif
1010    }
1011  
1012 <  Vector3d ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){
1012 >  Vector3d& ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){
1013   #ifdef IS_MPI
1014      return atomColData.velocity[atom2];
1015   #else
# Line 1027 | Line 1047 | namespace OpenMD {
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 1035 | 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 1058 | Line 1078 | namespace OpenMD {
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 1149 | 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 1205 | Line 1227 | namespace OpenMD {
1227   #else
1228      
1229      idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
1230 +    idat.atid1 = idents[atom1];
1231 +    idat.atid2 = idents[atom2];
1232  
1233      if (storageLayout_ & DataStorage::dslAmat) {
1234        idat.A1 = &(snap_->atomData.aMat[atom1]);
# Line 1316 | 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;
# Line 1506 | Line 1531 | namespace OpenMD {
1531                    if (usePeriodicBoundaryConditions_) {
1532                      snap_->wrapVector(dr);
1533                    }
1534 <                  cuts = getGroupCutoffs( (*j1), (*j2) );
1535 <                  if (dr.lengthSquare() < cuts.third) {
1534 >                  getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq );
1535 >                  if (dr.lengthSquare() < rlistsq) {
1536                      neighborList.push_back(make_pair((*j1), (*j2)));
1537                    }                  
1538                  }
# Line 1533 | Line 1558 | namespace OpenMD {
1558                      if (usePeriodicBoundaryConditions_) {
1559                        snap_->wrapVector(dr);
1560                      }
1561 <                    cuts = getGroupCutoffs( (*j1), (*j2) );
1562 <                    if (dr.lengthSquare() < cuts.third) {
1561 >                    getGroupCutoffs( (*j1), (*j2), rcut, rcutsq, rlistsq );
1562 >                    if (dr.lengthSquare() < rlistsq) {
1563                        neighborList.push_back(make_pair((*j1), (*j2)));
1564                      }
1565                    }
# Line 1554 | Line 1579 | namespace OpenMD {
1579            if (usePeriodicBoundaryConditions_) {
1580              snap_->wrapVector(dr);
1581            }
1582 <          cuts = getGroupCutoffs( j1, j2 );
1583 <          if (dr.lengthSquare() < cuts.third) {
1582 >          getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq);
1583 >          if (dr.lengthSquare() < rlistsq) {
1584              neighborList.push_back(make_pair(j1, j2));
1585            }
1586          }
# Line 1569 | Line 1594 | namespace OpenMD {
1594            if (usePeriodicBoundaryConditions_) {
1595              snap_->wrapVector(dr);
1596            }
1597 <          cuts = getGroupCutoffs( j1, j2 );
1598 <          if (dr.lengthSquare() < cuts.third) {
1597 >          getGroupCutoffs( j1, j2, rcut, rcutsq, rlistsq );
1598 >          if (dr.lengthSquare() < rlistsq) {
1599              neighborList.push_back(make_pair(j1, j2));
1600            }
1601          }    
# Line 1583 | 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]);
1586 –    
1587 –    return neighborList;
1611    }
1612   } //end namespace OpenMD

Diff Legend

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