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

Diff Legend

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