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 1893 by gezelter, Wed Jun 19 17:19:07 2013 UTC vs.
Revision 2061 by gezelter, Tue Mar 3 16:24:44 2015 UTC

# Line 50 | Line 50 | namespace OpenMD {
50  
51    ForceMatrixDecomposition::ForceMatrixDecomposition(SimInfo* info, InteractionManager* iMan) : ForceDecomposition(info, iMan) {
52  
53 <    // In a parallel computation, row and colum scans must visit all
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
53 >    // Row and colum scans must visit all surrounding cells
54      cellOffsets_.clear();
55      cellOffsets_.push_back( Vector3i(-1,-1,-1) );
56      cellOffsets_.push_back( Vector3i( 0,-1,-1) );
# Line 82 | Line 79 | namespace OpenMD {
79      cellOffsets_.push_back( Vector3i(-1, 1, 1) );
80      cellOffsets_.push_back( Vector3i( 0, 1, 1) );
81      cellOffsets_.push_back( Vector3i( 1, 1, 1) );
85 – #endif    
82    }
83  
84  
# Line 99 | Line 95 | namespace OpenMD {
95      nGroups_ = info_->getNLocalCutoffGroups();
96      // gather the information for atomtype IDs (atids):
97      idents = info_->getIdentArray();
98 +    regions = info_->getRegions();
99      AtomLocalToGlobal = info_->getGlobalAtomIndices();
100      cgLocalToGlobal = info_->getGlobalGroupIndices();
101      vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
# Line 118 | Line 115 | namespace OpenMD {
115      
116   #ifdef IS_MPI
117  
118 <    MPI::Intracomm row = rowComm.getComm();
119 <    MPI::Intracomm col = colComm.getComm();
118 >    MPI_Comm row = rowComm.getComm();
119 >    MPI_Comm col = colComm.getComm();
120  
121      AtomPlanIntRow = new Plan<int>(row, nLocal_);
122      AtomPlanRealRow = new Plan<RealType>(row, nLocal_);
# Line 163 | Line 160 | namespace OpenMD {
160      
161      AtomPlanIntRow->gather(idents, identsRow);
162      AtomPlanIntColumn->gather(idents, identsCol);
163 +
164 +    regionsRow.resize(nAtomsInRow_);
165 +    regionsCol.resize(nAtomsInCol_);
166      
167 +    AtomPlanIntRow->gather(regions, regionsRow);
168 +    AtomPlanIntColumn->gather(regions, regionsCol);
169 +    
170      // allocate memory for the parallel objects
171      atypesRow.resize(nAtomsInRow_);
172      atypesCol.resize(nAtomsInCol_);
# Line 299 | Line 302 | namespace OpenMD {
302            groupList_[i].push_back(j);
303          }
304        }      
305 <    }
303 <
304 <
305 <    createGtypeCutoffMap();
306 <
305 >    }    
306    }
308 –  
309 –  void ForceMatrixDecomposition::createGtypeCutoffMap() {
310 –    
311 –    RealType tol = 1e-6;
312 –    largestRcut_ = 0.0;
313 –    int atid;
314 –    set<AtomType*> atypes = info_->getSimulatedAtomTypes();
315 –    
316 –    map<int, RealType> atypeCutoff;
317 –      
318 –    for (set<AtomType*>::iterator at = atypes.begin();
319 –         at != atypes.end(); ++at){
320 –      atid = (*at)->getIdent();
321 –      if (userChoseCutoff_)
322 –        atypeCutoff[atid] = userCutoff_;
323 –      else
324 –        atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at);
325 –    }
326 –    
327 –    vector<RealType> gTypeCutoffs;
328 –    // first we do a single loop over the cutoff groups to find the
329 –    // largest cutoff for any atypes present in this group.
330 – #ifdef IS_MPI
331 –    vector<RealType> groupCutoffRow(nGroupsInRow_, 0.0);
332 –    groupRowToGtype.resize(nGroupsInRow_);
333 –    for (int cg1 = 0; cg1 < nGroupsInRow_; cg1++) {
334 –      vector<int> atomListRow = getAtomsInGroupRow(cg1);
335 –      for (vector<int>::iterator ia = atomListRow.begin();
336 –           ia != atomListRow.end(); ++ia) {            
337 –        int atom1 = (*ia);
338 –        atid = identsRow[atom1];
339 –        if (atypeCutoff[atid] > groupCutoffRow[cg1]) {
340 –          groupCutoffRow[cg1] = atypeCutoff[atid];
341 –        }
342 –      }
343 –
344 –      bool gTypeFound = false;
345 –      for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
346 –        if (abs(groupCutoffRow[cg1] - gTypeCutoffs[gt]) < tol) {
347 –          groupRowToGtype[cg1] = gt;
348 –          gTypeFound = true;
349 –        }
350 –      }
351 –      if (!gTypeFound) {
352 –        gTypeCutoffs.push_back( groupCutoffRow[cg1] );
353 –        groupRowToGtype[cg1] = gTypeCutoffs.size() - 1;
354 –      }
355 –      
356 –    }
357 –    vector<RealType> groupCutoffCol(nGroupsInCol_, 0.0);
358 –    groupColToGtype.resize(nGroupsInCol_);
359 –    for (int cg2 = 0; cg2 < nGroupsInCol_; cg2++) {
360 –      vector<int> atomListCol = getAtomsInGroupColumn(cg2);
361 –      for (vector<int>::iterator jb = atomListCol.begin();
362 –           jb != atomListCol.end(); ++jb) {            
363 –        int atom2 = (*jb);
364 –        atid = identsCol[atom2];
365 –        if (atypeCutoff[atid] > groupCutoffCol[cg2]) {
366 –          groupCutoffCol[cg2] = atypeCutoff[atid];
367 –        }
368 –      }
369 –      bool gTypeFound = false;
370 –      for (int gt = 0; gt < gTypeCutoffs.size(); gt++) {
371 –        if (abs(groupCutoffCol[cg2] - gTypeCutoffs[gt]) < tol) {
372 –          groupColToGtype[cg2] = gt;
373 –          gTypeFound = true;
374 –        }
375 –      }
376 –      if (!gTypeFound) {
377 –        gTypeCutoffs.push_back( groupCutoffCol[cg2] );
378 –        groupColToGtype[cg2] = gTypeCutoffs.size() - 1;
379 –      }
380 –    }
381 – #else
382 –
383 –    vector<RealType> groupCutoff(nGroups_, 0.0);
384 –    groupToGtype.resize(nGroups_);
385 –    for (int cg1 = 0; cg1 < nGroups_; cg1++) {
386 –      groupCutoff[cg1] = 0.0;
387 –      vector<int> atomList = getAtomsInGroupRow(cg1);
388 –      for (vector<int>::iterator ia = atomList.begin();
389 –           ia != atomList.end(); ++ia) {            
390 –        int atom1 = (*ia);
391 –        atid = idents[atom1];
392 –        if (atypeCutoff[atid] > groupCutoff[cg1])
393 –          groupCutoff[cg1] = atypeCutoff[atid];
394 –      }
395 –      
396 –      bool gTypeFound = false;
397 –      for (unsigned int gt = 0; gt < gTypeCutoffs.size(); gt++) {
398 –        if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) {
399 –          groupToGtype[cg1] = gt;
400 –          gTypeFound = true;
401 –        }
402 –      }
403 –      if (!gTypeFound) {      
404 –        gTypeCutoffs.push_back( groupCutoff[cg1] );
405 –        groupToGtype[cg1] = gTypeCutoffs.size() - 1;
406 –      }      
407 –    }
408 – #endif
409 –
410 –    // Now we find the maximum group cutoff value present in the simulation
411 –
412 –    RealType groupMax = *max_element(gTypeCutoffs.begin(),
413 –                                     gTypeCutoffs.end());
414 –
415 – #ifdef IS_MPI
416 –    MPI::COMM_WORLD.Allreduce(&groupMax, &groupMax, 1, MPI::REALTYPE,
417 –                              MPI::MAX);
418 – #endif
307      
420 –    RealType tradRcut = groupMax;
421 –
422 –    for (unsigned int i = 0; i < gTypeCutoffs.size();  i++) {
423 –      for (unsigned int j = 0; j < gTypeCutoffs.size();  j++) {      
424 –        RealType thisRcut;
425 –        switch(cutoffPolicy_) {
426 –        case TRADITIONAL:
427 –          thisRcut = tradRcut;
428 –          break;
429 –        case MIX:
430 –          thisRcut = 0.5 * (gTypeCutoffs[i] + gTypeCutoffs[j]);
431 –          break;
432 –        case MAX:
433 –          thisRcut = max(gTypeCutoffs[i], gTypeCutoffs[j]);
434 –          break;
435 –        default:
436 –          sprintf(painCave.errMsg,
437 –                  "ForceMatrixDecomposition::createGtypeCutoffMap "
438 –                  "hit an unknown cutoff policy!\n");
439 –          painCave.severity = OPENMD_ERROR;
440 –          painCave.isFatal = 1;
441 –          simError();
442 –          break;
443 –        }
444 –
445 –        pair<int,int> key = make_pair(i,j);
446 –        gTypeCutoffMap[key].first = thisRcut;
447 –        if (thisRcut > largestRcut_) largestRcut_ = thisRcut;
448 –        gTypeCutoffMap[key].second = thisRcut*thisRcut;
449 –        gTypeCutoffMap[key].third = pow(thisRcut + skinThickness_, 2);
450 –        // sanity check
451 –        
452 –        if (userChoseCutoff_) {
453 –          if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) {
454 –            sprintf(painCave.errMsg,
455 –                    "ForceMatrixDecomposition::createGtypeCutoffMap "
456 –                    "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_);
457 –            painCave.severity = OPENMD_ERROR;
458 –            painCave.isFatal = 1;
459 –            simError();            
460 –          }
461 –        }
462 –      }
463 –    }
464 –  }
465 –
466 –  groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) {
467 –    int i, j;  
468 – #ifdef IS_MPI
469 –    i = groupRowToGtype[cg1];
470 –    j = groupColToGtype[cg2];
471 – #else
472 –    i = groupToGtype[cg1];
473 –    j = groupToGtype[cg2];
474 – #endif    
475 –    return gTypeCutoffMap[make_pair(i,j)];
476 –  }
477 –
308    int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) {
309      for (unsigned int j = 0; j < toposForAtom[atom1].size(); j++) {
310        if (toposForAtom[atom1][j] == atom2)
# Line 559 | Line 389 | namespace OpenMD {
389             atomColData.electricField.end(), V3Zero);
390      }
391  
392 +    if (storageLayout_ & DataStorage::dslSitePotential) {    
393 +      fill(atomRowData.sitePotential.begin(),
394 +           atomRowData.sitePotential.end(), 0.0);
395 +      fill(atomColData.sitePotential.begin(),
396 +           atomColData.sitePotential.end(), 0.0);
397 +    }
398 +
399   #endif
400      // even in parallel, we need to zero out the local arrays:
401  
# Line 591 | Line 428 | namespace OpenMD {
428        fill(snap_->atomData.electricField.begin(),
429             snap_->atomData.electricField.end(), V3Zero);
430      }
431 +    if (storageLayout_ & DataStorage::dslSitePotential) {      
432 +      fill(snap_->atomData.sitePotential.begin(),
433 +           snap_->atomData.sitePotential.end(), 0.0);
434 +    }
435    }
436  
437  
438    void ForceMatrixDecomposition::distributeData()  {
439 +  
440 + #ifdef IS_MPI
441 +
442      snap_ = sman_->getCurrentSnapshot();
443      storageLayout_ = sman_->getStorageLayout();
600 – #ifdef IS_MPI
444      
445 +    bool needsCG = true;
446 +    if(info_->getNCutoffGroups() != info_->getNAtoms())
447 +      needsCG = false;
448 +
449      // gather up the atomic positions
450      AtomPlanVectorRow->gather(snap_->atomData.position,
451                                atomRowData.position);
# Line 607 | Line 454 | namespace OpenMD {
454      
455      // gather up the cutoff group positions
456  
457 <    cgPlanVectorRow->gather(snap_->cgData.position,
458 <                            cgRowData.position);
457 >    if (needsCG) {
458 >      cgPlanVectorRow->gather(snap_->cgData.position,
459 >                              cgRowData.position);
460 >      
461 >      cgPlanVectorColumn->gather(snap_->cgData.position,
462 >                                 cgColData.position);
463 >    }
464  
613 –    cgPlanVectorColumn->gather(snap_->cgData.position,
614 –                               cgColData.position);
465  
616 –
617 –
466      if (needVelocities_) {
467        // gather up the atomic velocities
468        AtomPlanVectorColumn->gather(snap_->atomData.velocity,
469                                     atomColData.velocity);
470 <      
471 <      cgPlanVectorColumn->gather(snap_->cgData.velocity,
472 <                                 cgColData.velocity);
470 >
471 >      if (needsCG) {        
472 >        cgPlanVectorColumn->gather(snap_->cgData.velocity,
473 >                                   cgColData.velocity);
474 >      }
475      }
476  
477      
# Line 805 | Line 655 | namespace OpenMD {
655          snap_->atomData.electricField[i] += efield_tmp[i];
656      }
657  
658 +    if (storageLayout_ & DataStorage::dslSitePotential) {
659  
660 +      int nsp = snap_->atomData.sitePotential.size();
661 +      vector<RealType> sp_tmp(nsp, 0.0);
662 +
663 +      AtomPlanRealRow->scatter(atomRowData.sitePotential, sp_tmp);
664 +      for (int i = 0; i < nsp; i++) {
665 +        snap_->atomData.sitePotential[i] += sp_tmp[i];
666 +        sp_tmp[i] = 0.0;
667 +      }
668 +      
669 +      AtomPlanRealColumn->scatter(atomColData.sitePotential, sp_tmp);
670 +      for (int i = 0; i < nsp; i++)
671 +        snap_->atomData.sitePotential[i] += sp_tmp[i];
672 +    }
673 +
674      nLocal_ = snap_->getNumberOfAtoms();
675  
676      vector<potVec> pot_temp(nLocal_,
# Line 889 | Line 754 | namespace OpenMD {
754      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
755        RealType ploc1 = pairwisePot[ii];
756        RealType ploc2 = 0.0;
757 <      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
757 >      MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
758        pairwisePot[ii] = ploc2;
759      }
760  
761      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
762        RealType ploc1 = excludedPot[ii];
763        RealType ploc2 = 0.0;
764 <      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
764 >      MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
765        excludedPot[ii] = ploc2;
766      }
767  
768      // Here be dragons.
769 <    MPI::Intracomm col = colComm.getComm();
769 >    MPI_Comm col = colComm.getComm();
770  
771 <    col.Allreduce(MPI::IN_PLACE,
771 >    MPI_Allreduce(MPI_IN_PLACE,
772                    &snap_->frameData.conductiveHeatFlux[0], 3,
773 <                  MPI::REALTYPE, MPI::SUM);
773 >                  MPI_REALTYPE, MPI_SUM, col);
774  
775  
776   #endif
# Line 924 | Line 789 | namespace OpenMD {
789      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
790        RealType ploc1 = embeddingPot[ii];
791        RealType ploc2 = 0.0;
792 <      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
792 >      MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
793        embeddingPot[ii] = ploc2;
794      }    
795      for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) {
796        RealType ploc1 = excludedSelfPot[ii];
797        RealType ploc2 = 0.0;
798 <      MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM);
798 >      MPI_Allreduce(&ploc1, &ploc2, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD);
799        excludedSelfPot[ii] = ploc2;
800      }    
801   #endif
# Line 1143 | Line 1008 | namespace OpenMD {
1008  
1009      // filling interaction blocks with pointers
1010    void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat,
1011 <                                                     int atom1, int atom2) {
1011 >                                                     int atom1, int atom2,
1012 >                                                     bool newAtom1) {
1013  
1014      idat.excluded = excludeAtomPair(atom1, atom2);
1015 <  
1015 >
1016 >    if (newAtom1) {
1017 >      
1018   #ifdef IS_MPI
1019 <    idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]);
1020 <    //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]),
1021 <    //                         ff_->getAtomType(identsCol[atom2]) );
1022 <    
1019 >      idat.atid1 = identsRow[atom1];
1020 >      idat.atid2 = identsCol[atom2];
1021 >      
1022 >      if (regionsRow[atom1] >= 0 && regionsCol[atom2] >= 0) {
1023 >        idat.sameRegion = (regionsRow[atom1] == regionsCol[atom2]);
1024 >      } else {
1025 >        idat.sameRegion = false;
1026 >      }
1027 >      
1028 >      if (storageLayout_ & DataStorage::dslAmat) {
1029 >        idat.A1 = &(atomRowData.aMat[atom1]);
1030 >        idat.A2 = &(atomColData.aMat[atom2]);
1031 >      }
1032 >      
1033 >      if (storageLayout_ & DataStorage::dslTorque) {
1034 >        idat.t1 = &(atomRowData.torque[atom1]);
1035 >        idat.t2 = &(atomColData.torque[atom2]);
1036 >      }
1037 >      
1038 >      if (storageLayout_ & DataStorage::dslDipole) {
1039 >        idat.dipole1 = &(atomRowData.dipole[atom1]);
1040 >        idat.dipole2 = &(atomColData.dipole[atom2]);
1041 >      }
1042 >      
1043 >      if (storageLayout_ & DataStorage::dslQuadrupole) {
1044 >        idat.quadrupole1 = &(atomRowData.quadrupole[atom1]);
1045 >        idat.quadrupole2 = &(atomColData.quadrupole[atom2]);
1046 >      }
1047 >      
1048 >      if (storageLayout_ & DataStorage::dslDensity) {
1049 >        idat.rho1 = &(atomRowData.density[atom1]);
1050 >        idat.rho2 = &(atomColData.density[atom2]);
1051 >      }
1052 >      
1053 >      if (storageLayout_ & DataStorage::dslFunctional) {
1054 >        idat.frho1 = &(atomRowData.functional[atom1]);
1055 >        idat.frho2 = &(atomColData.functional[atom2]);
1056 >      }
1057 >      
1058 >      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
1059 >        idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
1060 >        idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
1061 >      }
1062 >      
1063 >      if (storageLayout_ & DataStorage::dslParticlePot) {
1064 >        idat.particlePot1 = &(atomRowData.particlePot[atom1]);
1065 >        idat.particlePot2 = &(atomColData.particlePot[atom2]);
1066 >      }
1067 >      
1068 >      if (storageLayout_ & DataStorage::dslSkippedCharge) {              
1069 >        idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]);
1070 >        idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
1071 >      }
1072 >      
1073 >      if (storageLayout_ & DataStorage::dslFlucQPosition) {              
1074 >        idat.flucQ1 = &(atomRowData.flucQPos[atom1]);
1075 >        idat.flucQ2 = &(atomColData.flucQPos[atom2]);
1076 >      }
1077 >      
1078 > #else
1079 >      
1080 >      idat.atid1 = idents[atom1];
1081 >      idat.atid2 = idents[atom2];
1082 >      
1083 >      if (regions[atom1] >= 0 && regions[atom2] >= 0) {
1084 >        idat.sameRegion = (regions[atom1] == regions[atom2]);
1085 >      } else {
1086 >        idat.sameRegion = false;
1087 >      }
1088 >      
1089 >      if (storageLayout_ & DataStorage::dslAmat) {
1090 >        idat.A1 = &(snap_->atomData.aMat[atom1]);
1091 >        idat.A2 = &(snap_->atomData.aMat[atom2]);
1092 >      }
1093 >      
1094 >      if (storageLayout_ & DataStorage::dslTorque) {
1095 >        idat.t1 = &(snap_->atomData.torque[atom1]);
1096 >        idat.t2 = &(snap_->atomData.torque[atom2]);
1097 >      }
1098 >      
1099 >      if (storageLayout_ & DataStorage::dslDipole) {
1100 >        idat.dipole1 = &(snap_->atomData.dipole[atom1]);
1101 >        idat.dipole2 = &(snap_->atomData.dipole[atom2]);
1102 >      }
1103 >      
1104 >      if (storageLayout_ & DataStorage::dslQuadrupole) {
1105 >        idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]);
1106 >        idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]);
1107 >      }
1108 >      
1109 >      if (storageLayout_ & DataStorage::dslDensity) {    
1110 >        idat.rho1 = &(snap_->atomData.density[atom1]);
1111 >        idat.rho2 = &(snap_->atomData.density[atom2]);
1112 >      }
1113 >      
1114 >      if (storageLayout_ & DataStorage::dslFunctional) {
1115 >        idat.frho1 = &(snap_->atomData.functional[atom1]);
1116 >        idat.frho2 = &(snap_->atomData.functional[atom2]);
1117 >      }
1118 >      
1119 >      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
1120 >        idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
1121 >        idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
1122 >      }
1123 >      
1124 >      if (storageLayout_ & DataStorage::dslParticlePot) {
1125 >        idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
1126 >        idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
1127 >      }
1128 >      
1129 >      if (storageLayout_ & DataStorage::dslSkippedCharge) {
1130 >        idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
1131 >        idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
1132 >      }
1133 >      
1134 >      if (storageLayout_ & DataStorage::dslFlucQPosition) {              
1135 >        idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]);
1136 >        idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]);
1137 >      }
1138 > #endif
1139 >      
1140 >    } else {
1141 >      // atom1 is not new, so don't bother updating properties of that atom:
1142 > #ifdef IS_MPI
1143 >    idat.atid2 = identsCol[atom2];
1144 >
1145 >    if (regionsRow[atom1] >= 0 && regionsCol[atom2] >= 0) {
1146 >      idat.sameRegion = (regionsRow[atom1] == regionsCol[atom2]);
1147 >    } else {
1148 >      idat.sameRegion = false;
1149 >    }
1150 >
1151      if (storageLayout_ & DataStorage::dslAmat) {
1156 –      idat.A1 = &(atomRowData.aMat[atom1]);
1152        idat.A2 = &(atomColData.aMat[atom2]);
1153      }
1154      
1155      if (storageLayout_ & DataStorage::dslTorque) {
1161 –      idat.t1 = &(atomRowData.torque[atom1]);
1156        idat.t2 = &(atomColData.torque[atom2]);
1157      }
1158  
1159      if (storageLayout_ & DataStorage::dslDipole) {
1166 –      idat.dipole1 = &(atomRowData.dipole[atom1]);
1160        idat.dipole2 = &(atomColData.dipole[atom2]);
1161      }
1162  
1163      if (storageLayout_ & DataStorage::dslQuadrupole) {
1171 –      idat.quadrupole1 = &(atomRowData.quadrupole[atom1]);
1164        idat.quadrupole2 = &(atomColData.quadrupole[atom2]);
1165      }
1166  
1167      if (storageLayout_ & DataStorage::dslDensity) {
1176 –      idat.rho1 = &(atomRowData.density[atom1]);
1168        idat.rho2 = &(atomColData.density[atom2]);
1169      }
1170  
1171      if (storageLayout_ & DataStorage::dslFunctional) {
1181 –      idat.frho1 = &(atomRowData.functional[atom1]);
1172        idat.frho2 = &(atomColData.functional[atom2]);
1173      }
1174  
1175      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
1186 –      idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
1176        idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
1177      }
1178  
1179      if (storageLayout_ & DataStorage::dslParticlePot) {
1191 –      idat.particlePot1 = &(atomRowData.particlePot[atom1]);
1180        idat.particlePot2 = &(atomColData.particlePot[atom2]);
1181      }
1182  
1183      if (storageLayout_ & DataStorage::dslSkippedCharge) {              
1196 –      idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]);
1184        idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]);
1185      }
1186  
1187 <    if (storageLayout_ & DataStorage::dslFlucQPosition) {              
1201 <      idat.flucQ1 = &(atomRowData.flucQPos[atom1]);
1187 >    if (storageLayout_ & DataStorage::dslFlucQPosition) {
1188        idat.flucQ2 = &(atomColData.flucQPos[atom2]);
1189      }
1190  
1191 < #else
1192 <    
1207 <    idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]);
1191 > #else  
1192 >    idat.atid2 = idents[atom2];
1193  
1194 +    if (regions[atom1] >= 0 && regions[atom2] >= 0) {
1195 +      idat.sameRegion = (regions[atom1] == regions[atom2]);
1196 +    } else {
1197 +      idat.sameRegion = false;
1198 +    }
1199 +
1200      if (storageLayout_ & DataStorage::dslAmat) {
1210 –      idat.A1 = &(snap_->atomData.aMat[atom1]);
1201        idat.A2 = &(snap_->atomData.aMat[atom2]);
1202      }
1203  
1204      if (storageLayout_ & DataStorage::dslTorque) {
1215 –      idat.t1 = &(snap_->atomData.torque[atom1]);
1205        idat.t2 = &(snap_->atomData.torque[atom2]);
1206      }
1207  
1208      if (storageLayout_ & DataStorage::dslDipole) {
1220 –      idat.dipole1 = &(snap_->atomData.dipole[atom1]);
1209        idat.dipole2 = &(snap_->atomData.dipole[atom2]);
1210      }
1211  
1212      if (storageLayout_ & DataStorage::dslQuadrupole) {
1225 –      idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]);
1213        idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]);
1214      }
1215  
1216      if (storageLayout_ & DataStorage::dslDensity) {    
1230 –      idat.rho1 = &(snap_->atomData.density[atom1]);
1217        idat.rho2 = &(snap_->atomData.density[atom2]);
1218      }
1219  
1220      if (storageLayout_ & DataStorage::dslFunctional) {
1235 –      idat.frho1 = &(snap_->atomData.functional[atom1]);
1221        idat.frho2 = &(snap_->atomData.functional[atom2]);
1222      }
1223  
1224      if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
1240 –      idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
1225        idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
1226      }
1227  
1228      if (storageLayout_ & DataStorage::dslParticlePot) {
1245 –      idat.particlePot1 = &(snap_->atomData.particlePot[atom1]);
1229        idat.particlePot2 = &(snap_->atomData.particlePot[atom2]);
1230      }
1231  
1232      if (storageLayout_ & DataStorage::dslSkippedCharge) {
1250 –      idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]);
1233        idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]);
1234      }
1235  
1236      if (storageLayout_ & DataStorage::dslFlucQPosition) {              
1255 –      idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]);
1237        idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]);
1238      }
1239  
1240   #endif
1241 +    }
1242    }
1261 –
1243    
1244 <  void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) {    
1244 >  void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat,
1245 >                                                       int atom1, int atom2) {  
1246   #ifdef IS_MPI
1247      pot_row[atom1] += RealType(0.5) *  *(idat.pot);
1248      pot_col[atom2] += RealType(0.5) *  *(idat.pot);
# Line 1280 | Line 1262 | namespace OpenMD {
1262        atomColData.electricField[atom2] += *(idat.eField2);
1263      }
1264  
1265 +    if (storageLayout_ & DataStorage::dslSitePotential) {              
1266 +      atomRowData.sitePotential[atom1] += *(idat.sPot1);
1267 +      atomColData.sitePotential[atom2] += *(idat.sPot2);
1268 +    }
1269 +
1270   #else
1271      pairwisePot += *(idat.pot);
1272      excludedPot += *(idat.excludedPot);
# Line 1306 | Line 1293 | namespace OpenMD {
1293        snap_->atomData.electricField[atom2] += *(idat.eField2);
1294      }
1295  
1296 +    if (storageLayout_ & DataStorage::dslSitePotential) {              
1297 +      snap_->atomData.sitePotential[atom1] += *(idat.sPot1);
1298 +      snap_->atomData.sitePotential[atom2] += *(idat.sPot2);
1299 +    }
1300 +
1301   #endif
1302      
1303    }
# Line 1313 | Line 1305 | namespace OpenMD {
1305    /*
1306     * buildNeighborList
1307     *
1308 <   * first element of pair is row-indexed CutoffGroup
1309 <   * second element of pair is column-indexed CutoffGroup
1308 >   * Constructs the Verlet neighbor list for a force-matrix
1309 >   * decomposition.  In this case, each processor is responsible for
1310 >   * row-site interactions with column-sites.
1311 >   *
1312 >   * neighborList is returned as a packed array of neighboring
1313 >   * column-ordered CutoffGroups.  The starting position in
1314 >   * neighborList for each row-ordered CutoffGroup is given by the
1315 >   * returned vector point.
1316     */
1317 <  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
1318 <      
1319 <    vector<pair<int, int> > neighborList;
1320 <    groupCutoffs cuts;
1317 >  void ForceMatrixDecomposition::buildNeighborList(vector<int>& neighborList,
1318 >                                                   vector<int>& point) {
1319 >    neighborList.clear();
1320 >    point.clear();
1321 >    int len = 0;
1322 >    
1323      bool doAllPairs = false;
1324  
1325 –    RealType rList_ = (largestRcut_ + skinThickness_);
1325      Snapshot* snap_ = sman_->getCurrentSnapshot();
1326      Mat3x3d box;
1327      Mat3x3d invBox;
# Line 1334 | Line 1333 | namespace OpenMD {
1333   #ifdef IS_MPI
1334      cellListRow_.clear();
1335      cellListCol_.clear();
1336 +    point.resize(nGroupsInRow_+1);
1337   #else
1338      cellList_.clear();
1339 +    point.resize(nGroups_+1);
1340   #endif
1341      
1342      if (!usePeriodicBoundaryConditions_) {
# Line 1346 | Line 1347 | namespace OpenMD {
1347        invBox = snap_->getInvHmat();
1348      }
1349      
1350 <    Vector3d boxX = box.getColumn(0);
1351 <    Vector3d boxY = box.getColumn(1);
1352 <    Vector3d boxZ = box.getColumn(2);
1350 >    Vector3d A = box.getColumn(0);
1351 >    Vector3d B = box.getColumn(1);
1352 >    Vector3d C = box.getColumn(2);
1353 >
1354 >    // Required for triclinic cells
1355 >    Vector3d AxB = cross(A, B);
1356 >    Vector3d BxC = cross(B, C);
1357 >    Vector3d CxA = cross(C, A);
1358 >
1359 >    // unit vectors perpendicular to the faces of the triclinic cell:
1360 >    AxB.normalize();
1361 >    BxC.normalize();
1362 >    CxA.normalize();
1363 >
1364 >    // A set of perpendicular lengths in triclinic cells:
1365 >    RealType Wa = abs(dot(A, BxC));
1366 >    RealType Wb = abs(dot(B, CxA));
1367 >    RealType Wc = abs(dot(C, AxB));
1368      
1369 <    nCells_.x() = (int) ( boxX.length() )/ rList_;
1370 <    nCells_.y() = (int) ( boxY.length() )/ rList_;
1371 <    nCells_.z() = (int) ( boxZ.length() )/ rList_;
1369 >    nCells_.x() = int( Wa / rList_ );
1370 >    nCells_.y() = int( Wb / rList_ );
1371 >    nCells_.z() = int( Wc / rList_ );
1372      
1373      // handle small boxes where the cell offsets can end up repeating cells
1358 –    
1374      if (nCells_.x() < 3) doAllPairs = true;
1375      if (nCells_.y() < 3) doAllPairs = true;
1376      if (nCells_.z() < 3) doAllPairs = true;
# Line 1370 | Line 1385 | namespace OpenMD {
1385   #endif
1386      
1387      if (!doAllPairs) {
1388 +      
1389   #ifdef IS_MPI
1390        
1391        for (int i = 0; i < nGroupsInRow_; i++) {
# Line 1428 | Line 1444 | namespace OpenMD {
1444          // add this cutoff group to the list of groups in this cell;
1445          cellListCol_[cellIndex].push_back(i);
1446        }
1447 <      
1447 >            
1448   #else
1449        for (int i = 0; i < nGroups_; i++) {
1450          rs = snap_->cgData.position[i];
# Line 1448 | Line 1464 | namespace OpenMD {
1464          }
1465          
1466          // find xyz-indices of cell that cutoffGroup is in.
1467 <        whichCell.x() = nCells_.x() * scaled.x();
1468 <        whichCell.y() = nCells_.y() * scaled.y();
1469 <        whichCell.z() = nCells_.z() * scaled.z();
1467 >        whichCell.x() = int(nCells_.x() * scaled.x());
1468 >        whichCell.y() = int(nCells_.y() * scaled.y());
1469 >        whichCell.z() = int(nCells_.z() * scaled.z());
1470          
1471          // find single index of this cell:
1472          cellIndex = Vlinear(whichCell, nCells_);
# Line 1461 | Line 1477 | namespace OpenMD {
1477  
1478   #endif
1479  
1480 <      for (int m1z = 0; m1z < nCells_.z(); m1z++) {
1481 <        for (int m1y = 0; m1y < nCells_.y(); m1y++) {
1482 <          for (int m1x = 0; m1x < nCells_.x(); m1x++) {
1483 <            Vector3i m1v(m1x, m1y, m1z);
1468 <            int m1 = Vlinear(m1v, nCells_);
1469 <            
1470 <            for (vector<Vector3i>::iterator os = cellOffsets_.begin();
1471 <                 os != cellOffsets_.end(); ++os) {
1472 <              
1473 <              Vector3i m2v = m1v + (*os);
1474 <            
1480 > #ifdef IS_MPI
1481 >      for (int j1 = 0; j1 < nGroupsInRow_; j1++) {
1482 >        rs = cgRowData.position[j1];
1483 > #else
1484  
1485 <              if (m2v.x() >= nCells_.x()) {
1486 <                m2v.x() = 0;          
1487 <              } else if (m2v.x() < 0) {
1488 <                m2v.x() = nCells_.x() - 1;
1489 <              }
1490 <              
1491 <              if (m2v.y() >= nCells_.y()) {
1492 <                m2v.y() = 0;          
1493 <              } else if (m2v.y() < 0) {
1494 <                m2v.y() = nCells_.y() - 1;
1495 <              }
1496 <              
1497 <              if (m2v.z() >= nCells_.z()) {
1498 <                m2v.z() = 0;          
1499 <              } else if (m2v.z() < 0) {
1500 <                m2v.z() = nCells_.z() - 1;
1501 <              }
1485 >      for (int j1 = 0; j1 < nGroups_; j1++) {
1486 >        rs = snap_->cgData.position[j1];
1487 > #endif
1488 >        point[j1] = len;
1489 >        
1490 >        // scaled positions relative to the box vectors
1491 >        scaled = invBox * rs;
1492 >        
1493 >        // wrap the vector back into the unit box by subtracting integer box
1494 >        // numbers
1495 >        for (int j = 0; j < 3; j++) {
1496 >          scaled[j] -= roundMe(scaled[j]);
1497 >          scaled[j] += 0.5;
1498 >          // Handle the special case when an object is exactly on the
1499 >          // boundary (a scaled coordinate of 1.0 is the same as
1500 >          // scaled coordinate of 0.0)
1501 >          if (scaled[j] >= 1.0) scaled[j] -= 1.0;
1502 >        }
1503 >        
1504 >        // find xyz-indices of cell that cutoffGroup is in.
1505 >        whichCell.x() = nCells_.x() * scaled.x();
1506 >        whichCell.y() = nCells_.y() * scaled.y();
1507 >        whichCell.z() = nCells_.z() * scaled.z();
1508 >        
1509 >        // find single index of this cell:
1510 >        int m1 = Vlinear(whichCell, nCells_);
1511  
1512 <              int m2 = Vlinear (m2v, nCells_);
1512 >        for (vector<Vector3i>::iterator os = cellOffsets_.begin();
1513 >             os != cellOffsets_.end(); ++os) {
1514                
1515 +          Vector3i m2v = whichCell + (*os);
1516 +
1517 +          if (m2v.x() >= nCells_.x()) {
1518 +            m2v.x() = 0;          
1519 +          } else if (m2v.x() < 0) {
1520 +            m2v.x() = nCells_.x() - 1;
1521 +          }
1522 +          
1523 +          if (m2v.y() >= nCells_.y()) {
1524 +            m2v.y() = 0;          
1525 +          } else if (m2v.y() < 0) {
1526 +            m2v.y() = nCells_.y() - 1;
1527 +          }
1528 +          
1529 +          if (m2v.z() >= nCells_.z()) {
1530 +            m2v.z() = 0;          
1531 +          } else if (m2v.z() < 0) {
1532 +            m2v.z() = nCells_.z() - 1;
1533 +          }
1534 +          int m2 = Vlinear (m2v, nCells_);                                      
1535   #ifdef IS_MPI
1536 <              for (vector<int>::iterator j1 = cellListRow_[m1].begin();
1537 <                   j1 != cellListRow_[m1].end(); ++j1) {
1538 <                for (vector<int>::iterator j2 = cellListCol_[m2].begin();
1539 <                     j2 != cellListCol_[m2].end(); ++j2) {
1540 <                  
1541 <                  // In parallel, we need to visit *all* pairs of row
1542 <                  // & column indicies and will divide labor in the
1543 <                  // force evaluation later.
1544 <                  dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)];
1545 <                  if (usePeriodicBoundaryConditions_) {
1546 <                    snap_->wrapVector(dr);
1547 <                  }
1548 <                  cuts = getGroupCutoffs( (*j1), (*j2) );
1549 <                  if (dr.lengthSquare() < cuts.third) {
1550 <                    neighborList.push_back(make_pair((*j1), (*j2)));
1512 <                  }                  
1513 <                }
1514 <              }
1536 >          for (vector<int>::iterator j2 = cellListCol_[m2].begin();
1537 >               j2 != cellListCol_[m2].end(); ++j2) {
1538 >            
1539 >            // In parallel, we need to visit *all* pairs of row
1540 >            // & column indicies and will divide labor in the
1541 >            // force evaluation later.
1542 >            dr = cgColData.position[(*j2)] - rs;
1543 >            if (usePeriodicBoundaryConditions_) {
1544 >              snap_->wrapVector(dr);
1545 >            }
1546 >            if (dr.lengthSquare() < rListSq_) {
1547 >              neighborList.push_back( (*j2) );
1548 >              ++len;
1549 >            }                
1550 >          }        
1551   #else
1552 <              for (vector<int>::iterator j1 = cellList_[m1].begin();
1553 <                   j1 != cellList_[m1].end(); ++j1) {
1554 <                for (vector<int>::iterator j2 = cellList_[m2].begin();
1555 <                     j2 != cellList_[m2].end(); ++j2) {
1556 <    
1557 <                  // Always do this if we're in different cells or if
1558 <                  // we're in the same cell and the global index of
1559 <                  // the j2 cutoff group is greater than or equal to
1560 <                  // the j1 cutoff group.  Note that Rappaport's code
1561 <                  // has a "less than" conditional here, but that
1562 <                  // deals with atom-by-atom computation.  OpenMD
1563 <                  // allows atoms within a single cutoff group to
1564 <                  // interact with each other.
1565 <
1566 <                  if (m2 != m1 || (*j2) >= (*j1) ) {
1567 <
1568 <                    dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
1533 <                    if (usePeriodicBoundaryConditions_) {
1534 <                      snap_->wrapVector(dr);
1535 <                    }
1536 <                    cuts = getGroupCutoffs( (*j1), (*j2) );
1537 <                    if (dr.lengthSquare() < cuts.third) {
1538 <                      neighborList.push_back(make_pair((*j1), (*j2)));
1539 <                    }
1540 <                  }
1541 <                }
1552 >          for (vector<int>::iterator j2 = cellList_[m2].begin();
1553 >               j2 != cellList_[m2].end(); ++j2) {
1554 >          
1555 >            // Always do this if we're in different cells or if
1556 >            // we're in the same cell and the global index of
1557 >            // the j2 cutoff group is greater than or equal to
1558 >            // the j1 cutoff group.  Note that Rappaport's code
1559 >            // has a "less than" conditional here, but that
1560 >            // deals with atom-by-atom computation.  OpenMD
1561 >            // allows atoms within a single cutoff group to
1562 >            // interact with each other.
1563 >            
1564 >            if ( (*j2) >= j1 ) {
1565 >              
1566 >              dr = snap_->cgData.position[(*j2)] - rs;
1567 >              if (usePeriodicBoundaryConditions_) {
1568 >                snap_->wrapVector(dr);
1569                }
1570 < #endif
1570 >              if ( dr.lengthSquare() < rListSq_) {
1571 >                neighborList.push_back( (*j2) );
1572 >                ++len;
1573 >              }
1574              }
1575 <          }
1575 >          }                
1576 > #endif
1577          }
1578 <      }
1578 >      }      
1579      } else {
1580        // branch to do all cutoff group pairs
1581   #ifdef IS_MPI
1582        for (int j1 = 0; j1 < nGroupsInRow_; j1++) {
1583 +        point[j1] = len;
1584 +        rs = cgRowData.position[j1];
1585          for (int j2 = 0; j2 < nGroupsInCol_; j2++) {    
1586 <          dr = cgColData.position[j2] - cgRowData.position[j1];
1586 >          dr = cgColData.position[j2] - rs;
1587            if (usePeriodicBoundaryConditions_) {
1588              snap_->wrapVector(dr);
1589            }
1590 <          cuts = getGroupCutoffs( j1, j2 );
1591 <          if (dr.lengthSquare() < cuts.third) {
1592 <            neighborList.push_back(make_pair(j1, j2));
1590 >          if (dr.lengthSquare() < rListSq_) {
1591 >            neighborList.push_back( j2 );
1592 >            ++len;
1593            }
1594          }
1595        }      
1596   #else
1597        // include all groups here.
1598        for (int j1 = 0; j1 < nGroups_; j1++) {
1599 +        point[j1] = len;
1600 +        rs = snap_->cgData.position[j1];
1601          // include self group interactions j2 == j1
1602          for (int j2 = j1; j2 < nGroups_; j2++) {
1603 <          dr = snap_->cgData.position[j2] - snap_->cgData.position[j1];
1603 >          dr = snap_->cgData.position[j2] - rs;
1604            if (usePeriodicBoundaryConditions_) {
1605              snap_->wrapVector(dr);
1606            }
1607 <          cuts = getGroupCutoffs( j1, j2 );
1608 <          if (dr.lengthSquare() < cuts.third) {
1609 <            neighborList.push_back(make_pair(j1, j2));
1607 >          if (dr.lengthSquare() < rListSq_) {
1608 >            neighborList.push_back( j2 );
1609 >            ++len;
1610            }
1611          }    
1612        }
1613   #endif
1614      }
1615 <      
1615 >
1616 > #ifdef IS_MPI
1617 >    point[nGroupsInRow_] = len;
1618 > #else
1619 >    point[nGroups_] = len;
1620 > #endif
1621 >  
1622      // save the local cutoff group positions for the check that is
1623      // done on each loop:
1624      saved_CG_positions_.clear();
1625 +    saved_CG_positions_.reserve(nGroups_);
1626      for (int i = 0; i < nGroups_; i++)
1627        saved_CG_positions_.push_back(snap_->cgData.position[i]);
1586 –    
1587 –    return neighborList;
1628    }
1629   } //end namespace OpenMD

Diff Legend

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