| 69 |  | PairList oneTwo = info_->getOneTwoInteractions(); | 
| 70 |  | PairList oneThree = info_->getOneThreeInteractions(); | 
| 71 |  | PairList oneFour = info_->getOneFourInteractions(); | 
| 72 | – | vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0); | 
| 72 |  |  | 
| 73 |  | #ifdef IS_MPI | 
| 74 |  |  | 
| 76 |  | AtomCommRealRow = new Communicator<Row,RealType>(nLocal_); | 
| 77 |  | AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_); | 
| 78 |  | AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_); | 
| 79 | + | AtomCommPotRow = new Communicator<Row,potVec>(nLocal_); | 
| 80 |  |  | 
| 81 |  | AtomCommIntColumn = new Communicator<Column,int>(nLocal_); | 
| 82 |  | AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_); | 
| 83 |  | AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_); | 
| 84 |  | AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_); | 
| 85 | + | AtomCommPotColumn = new Communicator<Column,potVec>(nLocal_); | 
| 86 |  |  | 
| 87 |  | cgCommIntRow = new Communicator<Row,int>(nGroups_); | 
| 88 |  | cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_); | 
| 103 |  | cgRowData.setStorageLayout(DataStorage::dslPosition); | 
| 104 |  | cgColData.resize(nGroupsInCol_); | 
| 105 |  | cgColData.setStorageLayout(DataStorage::dslPosition); | 
| 106 | < |  | 
| 106 | < | vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES, | 
| 107 | < | vector<RealType> (nAtomsInRow_, 0.0)); | 
| 108 | < | vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES, | 
| 109 | < | vector<RealType> (nAtomsInCol_, 0.0)); | 
| 110 | < |  | 
| 106 | > |  | 
| 107 |  | identsRow.reserve(nAtomsInRow_); | 
| 108 |  | identsCol.reserve(nAtomsInCol_); | 
| 109 |  |  | 
| 228 |  | } | 
| 229 |  | } | 
| 230 |  |  | 
| 231 | + | void ForceMatrixDecomposition::zeroWorkArrays() { | 
| 232 | + |  | 
| 233 | + | for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { | 
| 234 | + | longRangePot_[j] = 0.0; | 
| 235 | + | } | 
| 236 | + |  | 
| 237 | + | #ifdef IS_MPI | 
| 238 | + | if (storageLayout_ & DataStorage::dslForce) { | 
| 239 | + | fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero); | 
| 240 | + | fill(atomColData.force.begin(), atomColData.force.end(), V3Zero); | 
| 241 | + | } | 
| 242 | + |  | 
| 243 | + | if (storageLayout_ & DataStorage::dslTorque) { | 
| 244 | + | fill(atomRowData.torque.begin(), atomRowData.torque.end(), V3Zero); | 
| 245 | + | fill(atomColData.torque.begin(), atomColData.torque.end(), V3Zero); | 
| 246 | + | } | 
| 247 | + |  | 
| 248 | + | fill(pot_row.begin(), pot_row.end(), | 
| 249 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 250 | + |  | 
| 251 | + | fill(pot_col.begin(), pot_col.end(), | 
| 252 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 253 | + |  | 
| 254 | + | pot_local = Vector<RealType, N_INTERACTION_FAMILIES>(0.0); | 
| 255 | + |  | 
| 256 | + | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 257 | + | fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); | 
| 258 | + | fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0); | 
| 259 | + | } | 
| 260 | + |  | 
| 261 | + | if (storageLayout_ & DataStorage::dslDensity) { | 
| 262 | + | fill(atomRowData.density.begin(), atomRowData.density.end(), 0.0); | 
| 263 | + | fill(atomColData.density.begin(), atomColData.density.end(), 0.0); | 
| 264 | + | } | 
| 265 | + |  | 
| 266 | + | if (storageLayout_ & DataStorage::dslFunctional) { | 
| 267 | + | fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0); | 
| 268 | + | fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0); | 
| 269 | + | } | 
| 270 | + |  | 
| 271 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | 
| 272 | + | fill(atomRowData.functionalDerivative.begin(), | 
| 273 | + | atomRowData.functionalDerivative.end(), 0.0); | 
| 274 | + | fill(atomColData.functionalDerivative.begin(), | 
| 275 | + | atomColData.functionalDerivative.end(), 0.0); | 
| 276 | + | } | 
| 277 | + |  | 
| 278 | + | #else | 
| 279 | + |  | 
| 280 | + | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 281 | + | fill(snap_->atomData.particlePot.begin(), | 
| 282 | + | snap_->atomData.particlePot.end(), 0.0); | 
| 283 | + | } | 
| 284 | + |  | 
| 285 | + | if (storageLayout_ & DataStorage::dslDensity) { | 
| 286 | + | fill(snap_->atomData.density.begin(), | 
| 287 | + | snap_->atomData.density.end(), 0.0); | 
| 288 | + | } | 
| 289 | + | if (storageLayout_ & DataStorage::dslFunctional) { | 
| 290 | + | fill(snap_->atomData.functional.begin(), | 
| 291 | + | snap_->atomData.functional.end(), 0.0); | 
| 292 | + | } | 
| 293 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | 
| 294 | + | fill(snap_->atomData.functionalDerivative.begin(), | 
| 295 | + | snap_->atomData.functionalDerivative.end(), 0.0); | 
| 296 | + | } | 
| 297 | + | #endif | 
| 298 | + |  | 
| 299 | + | } | 
| 300 | + |  | 
| 301 | + |  | 
| 302 |  | void ForceMatrixDecomposition::distributeData()  { | 
| 303 |  | snap_ = sman_->getCurrentSnapshot(); | 
| 304 |  | storageLayout_ = sman_->getStorageLayout(); | 
| 334 |  | #endif | 
| 335 |  | } | 
| 336 |  |  | 
| 337 | + | /* collects information obtained during the pre-pair loop onto local | 
| 338 | + | * data structures. | 
| 339 | + | */ | 
| 340 |  | void ForceMatrixDecomposition::collectIntermediateData() { | 
| 341 |  | snap_ = sman_->getCurrentSnapshot(); | 
| 342 |  | storageLayout_ = sman_->getStorageLayout(); | 
| 348 |  | snap_->atomData.density); | 
| 349 |  |  | 
| 350 |  | int n = snap_->atomData.density.size(); | 
| 351 | < | std::vector<RealType> rho_tmp(n, 0.0); | 
| 351 | > | vector<RealType> rho_tmp(n, 0.0); | 
| 352 |  | AtomCommRealColumn->scatter(atomColData.density, rho_tmp); | 
| 353 |  | for (int i = 0; i < n; i++) | 
| 354 |  | snap_->atomData.density[i] += rho_tmp[i]; | 
| 355 |  | } | 
| 356 |  | #endif | 
| 357 |  | } | 
| 358 | < |  | 
| 358 | > |  | 
| 359 | > | /* | 
| 360 | > | * redistributes information obtained during the pre-pair loop out to | 
| 361 | > | * row and column-indexed data structures | 
| 362 | > | */ | 
| 363 |  | void ForceMatrixDecomposition::distributeIntermediateData() { | 
| 364 |  | snap_ = sman_->getCurrentSnapshot(); | 
| 365 |  | storageLayout_ = sman_->getStorageLayout(); | 
| 417 |  |  | 
| 418 |  | nLocal_ = snap_->getNumberOfAtoms(); | 
| 419 |  |  | 
| 420 | < | vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES, | 
| 421 | < | vector<RealType> (nLocal_, 0.0)); | 
| 420 | > | vector<potVec> pot_temp(nLocal_, | 
| 421 | > | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 422 | > |  | 
| 423 | > | // scatter/gather pot_row into the members of my column | 
| 424 | > |  | 
| 425 | > | AtomCommPotRow->scatter(pot_row, pot_temp); | 
| 426 | > |  | 
| 427 | > | for (int ii = 0;  ii < pot_temp.size(); ii++ ) | 
| 428 | > | pot_local += pot_temp[ii]; | 
| 429 |  |  | 
| 430 | < | for (int i = 0; i < N_INTERACTION_FAMILIES; i++) { | 
| 431 | < | AtomCommRealRow->scatter(pot_row[i], pot_temp[i]); | 
| 432 | < | for (int ii = 0;  ii < pot_temp[i].size(); ii++ ) { | 
| 433 | < | pot_local[i] += pot_temp[i][ii]; | 
| 434 | < | } | 
| 435 | < | } | 
| 430 | > | fill(pot_temp.begin(), pot_temp.end(), | 
| 431 | > | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 432 | > |  | 
| 433 | > | AtomCommPotColumn->scatter(pot_col, pot_temp); | 
| 434 | > |  | 
| 435 | > | for (int ii = 0;  ii < pot_temp.size(); ii++ ) | 
| 436 | > | pot_local += pot_temp[ii]; | 
| 437 | > |  | 
| 438 |  | #endif | 
| 439 |  | } | 
| 440 |  |  | 
| 545 |  | } | 
| 546 |  |  | 
| 547 |  | /** | 
| 548 | < | * there are a number of reasons to skip a pair or a particle mostly | 
| 549 | < | * we do this to exclude atoms who are involved in short range | 
| 550 | < | * interactions (bonds, bends, torsions), but we also need to | 
| 551 | < | * exclude some overcounted interactions that result from the | 
| 552 | < | * parallel decomposition. | 
| 548 | > | * There are a number of reasons to skip a pair or a | 
| 549 | > | * particle. Mostly we do this to exclude atoms who are involved in | 
| 550 | > | * short range interactions (bonds, bends, torsions), but we also | 
| 551 | > | * need to exclude some overcounted interactions that result from | 
| 552 | > | * the parallel decomposition. | 
| 553 |  | */ | 
| 554 |  | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | 
| 555 |  | int unique_id_1, unique_id_2; | 
| 628 |  | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | 
| 629 |  | ff_->getAtomType(identsCol[atom2]) ); | 
| 630 |  |  | 
| 631 | + |  | 
| 632 |  | if (storageLayout_ & DataStorage::dslAmat) { | 
| 633 |  | idat.A1 = &(atomRowData.aMat[atom1]); | 
| 634 |  | idat.A2 = &(atomColData.aMat[atom2]); | 
| 647 |  | if (storageLayout_ & DataStorage::dslDensity) { | 
| 648 |  | idat.rho1 = &(atomRowData.density[atom1]); | 
| 649 |  | idat.rho2 = &(atomColData.density[atom2]); | 
| 650 | + | } | 
| 651 | + |  | 
| 652 | + | if (storageLayout_ & DataStorage::dslFunctional) { | 
| 653 | + | idat.frho1 = &(atomRowData.functional[atom1]); | 
| 654 | + | idat.frho2 = &(atomColData.functional[atom2]); | 
| 655 |  | } | 
| 656 |  |  | 
| 657 |  | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | 
| 659 |  | idat.dfrho2 = &(atomColData.functionalDerivative[atom2]); | 
| 660 |  | } | 
| 661 |  |  | 
| 662 | + | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 663 | + | idat.particlePot1 = &(atomRowData.particlePot[atom1]); | 
| 664 | + | idat.particlePot2 = &(atomColData.particlePot[atom2]); | 
| 665 | + | } | 
| 666 | + |  | 
| 667 |  | #else | 
| 668 |  |  | 
| 669 |  | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), | 
| 689 |  | idat.rho2 = &(snap_->atomData.density[atom2]); | 
| 690 |  | } | 
| 691 |  |  | 
| 692 | + | if (storageLayout_ & DataStorage::dslFunctional) { | 
| 693 | + | idat.frho1 = &(snap_->atomData.functional[atom1]); | 
| 694 | + | idat.frho2 = &(snap_->atomData.functional[atom2]); | 
| 695 | + | } | 
| 696 | + |  | 
| 697 |  | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { | 
| 698 |  | idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]); | 
| 699 |  | idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]); | 
| 700 |  | } | 
| 701 | + |  | 
| 702 | + | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 703 | + | idat.particlePot1 = &(snap_->atomData.particlePot[atom1]); | 
| 704 | + | idat.particlePot2 = &(snap_->atomData.particlePot[atom2]); | 
| 705 | + | } | 
| 706 | + |  | 
| 707 |  | #endif | 
| 708 |  | return idat; | 
| 709 | + | } | 
| 710 | + |  | 
| 711 | + |  | 
| 712 | + | void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { | 
| 713 | + | #ifdef IS_MPI | 
| 714 | + | pot_row[atom1] += 0.5 *  *(idat.pot); | 
| 715 | + | pot_col[atom2] += 0.5 *  *(idat.pot); | 
| 716 | + |  | 
| 717 | + | atomRowData.force[atom1] += *(idat.f1); | 
| 718 | + | atomColData.force[atom2] -= *(idat.f1); | 
| 719 | + | #else | 
| 720 | + | longRangePot_ += *(idat.pot); | 
| 721 | + |  | 
| 722 | + | snap_->atomData.force[atom1] += *(idat.f1); | 
| 723 | + | snap_->atomData.force[atom2] -= *(idat.f1); | 
| 724 | + | #endif | 
| 725 | + |  | 
| 726 |  | } | 
| 727 |  |  | 
| 728 | + |  | 
| 729 |  | InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ | 
| 730 |  |  | 
| 731 |  | InteractionData idat; | 
| 741 |  | idat.t1 = &(atomRowData.torque[atom1]); | 
| 742 |  | idat.t2 = &(atomColData.torque[atom2]); | 
| 743 |  | } | 
| 621 | – | if (storageLayout_ & DataStorage::dslForce) { | 
| 622 | – | idat.t1 = &(atomRowData.force[atom1]); | 
| 623 | – | idat.t2 = &(atomColData.force[atom2]); | 
| 624 | – | } | 
| 744 |  | #else | 
| 745 |  | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), | 
| 746 |  | ff_->getAtomType(identsLocal[atom2]) ); | 
| 753 |  | idat.t1 = &(snap_->atomData.torque[atom1]); | 
| 754 |  | idat.t2 = &(snap_->atomData.torque[atom2]); | 
| 755 |  | } | 
| 637 | – | if (storageLayout_ & DataStorage::dslForce) { | 
| 638 | – | idat.t1 = &(snap_->atomData.force[atom1]); | 
| 639 | – | idat.t2 = &(snap_->atomData.force[atom2]); | 
| 640 | – | } | 
| 756 |  | #endif | 
| 757 |  | } | 
| 758 |  |  | 
| 854 |  | } | 
| 855 |  | #endif | 
| 856 |  |  | 
| 742 | – |  | 
| 743 | – |  | 
| 857 |  | for (int m1z = 0; m1z < nCells_.z(); m1z++) { | 
| 858 |  | for (int m1y = 0; m1y < nCells_.y(); m1y++) { | 
| 859 |  | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |