| 175 |  |  | 
| 176 |  | pot_row.resize(nAtomsInRow_); | 
| 177 |  | pot_col.resize(nAtomsInCol_); | 
| 178 | + |  | 
| 179 | + | expot_row.resize(nAtomsInRow_); | 
| 180 | + | expot_col.resize(nAtomsInCol_); | 
| 181 |  |  | 
| 182 |  | AtomRowToGlobal.resize(nAtomsInRow_); | 
| 183 |  | AtomColToGlobal.resize(nAtomsInCol_); | 
| 464 |  | } | 
| 465 |  | } | 
| 466 |  |  | 
| 464 | – |  | 
| 467 |  | groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { | 
| 468 |  | int i, j; | 
| 469 |  | #ifdef IS_MPI | 
| 487 |  | void ForceMatrixDecomposition::zeroWorkArrays() { | 
| 488 |  | pairwisePot = 0.0; | 
| 489 |  | embeddingPot = 0.0; | 
| 490 | + | excludedPot = 0.0; | 
| 491 | + | excludedSelfPot = 0.0; | 
| 492 |  |  | 
| 493 |  | #ifdef IS_MPI | 
| 494 |  | if (storageLayout_ & DataStorage::dslForce) { | 
| 505 |  | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 506 |  |  | 
| 507 |  | fill(pot_col.begin(), pot_col.end(), | 
| 508 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 509 | + |  | 
| 510 | + | fill(expot_row.begin(), expot_row.end(), | 
| 511 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 512 | + |  | 
| 513 | + | fill(expot_col.begin(), expot_col.end(), | 
| 514 |  | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 515 |  |  | 
| 516 |  | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 791 |  |  | 
| 792 |  | vector<potVec> pot_temp(nLocal_, | 
| 793 |  | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 794 | + | vector<potVec> expot_temp(nLocal_, | 
| 795 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 796 |  |  | 
| 797 |  | // scatter/gather pot_row into the members of my column | 
| 798 |  |  | 
| 799 |  | AtomPlanPotRow->scatter(pot_row, pot_temp); | 
| 800 | + | AtomPlanPotRow->scatter(expot_row, expot_temp); | 
| 801 |  |  | 
| 802 | < | for (int ii = 0;  ii < pot_temp.size(); ii++ ) | 
| 802 | > | for (int ii = 0;  ii < pot_temp.size(); ii++ ) | 
| 803 |  | pairwisePot += pot_temp[ii]; | 
| 804 | + |  | 
| 805 | + | for (int ii = 0;  ii < expot_temp.size(); ii++ ) | 
| 806 | + | excludedPot += expot_temp[ii]; | 
| 807 |  |  | 
| 808 |  | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 809 |  | // This is the pairwise contribution to the particle pot.  The | 
| 821 |  |  | 
| 822 |  | fill(pot_temp.begin(), pot_temp.end(), | 
| 823 |  | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 824 | + | fill(expot_temp.begin(), expot_temp.end(), | 
| 825 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | 
| 826 |  |  | 
| 827 |  | AtomPlanPotColumn->scatter(pot_col, pot_temp); | 
| 828 | + | AtomPlanPotColumn->scatter(expot_col, expot_temp); | 
| 829 |  |  | 
| 830 |  | for (int ii = 0;  ii < pot_temp.size(); ii++ ) | 
| 831 |  | pairwisePot += pot_temp[ii]; | 
| 832 |  |  | 
| 833 | + | for (int ii = 0;  ii < expot_temp.size(); ii++ ) | 
| 834 | + | excludedPot += expot_temp[ii]; | 
| 835 | + |  | 
| 836 |  | if (storageLayout_ & DataStorage::dslParticlePot) { | 
| 837 |  | // This is the pairwise contribution to the particle pot.  The | 
| 838 |  | // embedding contribution is added in each of the low level | 
| 875 |  | } | 
| 876 |  |  | 
| 877 |  | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { | 
| 878 | < | RealType ploc1 = embeddingPot[ii]; | 
| 878 | > | RealType ploc1 = excludedPot[ii]; | 
| 879 |  | RealType ploc2 = 0.0; | 
| 880 |  | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); | 
| 881 | < | embeddingPot[ii] = ploc2; | 
| 881 | > | excludedPot[ii] = ploc2; | 
| 882 |  | } | 
| 883 | < |  | 
| 883 | > |  | 
| 884 |  | // Here be dragons. | 
| 885 |  | MPI::Intracomm col = colComm.getComm(); | 
| 886 |  |  | 
| 890 |  |  | 
| 891 |  |  | 
| 892 |  | #endif | 
| 893 | + |  | 
| 894 | + | } | 
| 895 | + |  | 
| 896 | + | /** | 
| 897 | + | * Collects information obtained during the post-pair (and embedding | 
| 898 | + | * functional) loops onto local data structures. | 
| 899 | + | */ | 
| 900 | + | void ForceMatrixDecomposition::collectSelfData() { | 
| 901 | + | snap_ = sman_->getCurrentSnapshot(); | 
| 902 | + | storageLayout_ = sman_->getStorageLayout(); | 
| 903 |  |  | 
| 904 | + | #ifdef IS_MPI | 
| 905 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { | 
| 906 | + | RealType ploc1 = embeddingPot[ii]; | 
| 907 | + | RealType ploc2 = 0.0; | 
| 908 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); | 
| 909 | + | embeddingPot[ii] = ploc2; | 
| 910 | + | } | 
| 911 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { | 
| 912 | + | RealType ploc1 = excludedSelfPot[ii]; | 
| 913 | + | RealType ploc2 = 0.0; | 
| 914 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); | 
| 915 | + | excludedSelfPot[ii] = ploc2; | 
| 916 | + | } | 
| 917 | + | #endif | 
| 918 | + |  | 
| 919 |  | } | 
| 920 |  |  | 
| 921 | + |  | 
| 922 | + |  | 
| 923 |  | int ForceMatrixDecomposition::getNAtomsInRow() { | 
| 924 |  | #ifdef IS_MPI | 
| 925 |  | return nAtomsInRow_; | 
| 1042 |  | * We need to exclude some overcounted interactions that result from | 
| 1043 |  | * the parallel decomposition. | 
| 1044 |  | */ | 
| 1045 | < | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | 
| 1046 | < | int unique_id_1, unique_id_2; | 
| 1045 | > | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) { | 
| 1046 | > | int unique_id_1, unique_id_2, group1, group2; | 
| 1047 |  |  | 
| 1048 |  | #ifdef IS_MPI | 
| 1049 |  | // in MPI, we have to look up the unique IDs for each atom | 
| 1050 |  | unique_id_1 = AtomRowToGlobal[atom1]; | 
| 1051 |  | unique_id_2 = AtomColToGlobal[atom2]; | 
| 1052 | + | group1 = cgRowToGlobal[cg1]; | 
| 1053 | + | group2 = cgColToGlobal[cg2]; | 
| 1054 |  | #else | 
| 1055 |  | unique_id_1 = AtomLocalToGlobal[atom1]; | 
| 1056 |  | unique_id_2 = AtomLocalToGlobal[atom2]; | 
| 1057 | + | group1 = cgLocalToGlobal[cg1]; | 
| 1058 | + | group2 = cgLocalToGlobal[cg2]; | 
| 1059 |  | #endif | 
| 1060 |  |  | 
| 1061 |  | if (unique_id_1 == unique_id_2) return true; | 
| 1067 |  | } else { | 
| 1068 |  | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; | 
| 1069 |  | } | 
| 1070 | + | #endif | 
| 1071 | + |  | 
| 1072 | + | #ifndef IS_MPI | 
| 1073 | + | if (group1 == group2) { | 
| 1074 | + | if (unique_id_1 < unique_id_2) return true; | 
| 1075 | + | } | 
| 1076 |  | #endif | 
| 1077 |  |  | 
| 1078 |  | return false; | 
| 1230 |  | #ifdef IS_MPI | 
| 1231 |  | pot_row[atom1] += RealType(0.5) *  *(idat.pot); | 
| 1232 |  | pot_col[atom2] += RealType(0.5) *  *(idat.pot); | 
| 1233 | + | expot_row[atom1] += RealType(0.5) *  *(idat.excludedPot); | 
| 1234 | + | expot_col[atom2] += RealType(0.5) *  *(idat.excludedPot); | 
| 1235 |  |  | 
| 1236 |  | atomRowData.force[atom1] += *(idat.f1); | 
| 1237 |  | atomColData.force[atom2] -= *(idat.f1); | 
| 1248 |  |  | 
| 1249 |  | #else | 
| 1250 |  | pairwisePot += *(idat.pot); | 
| 1251 | + | excludedPot += *(idat.excludedPot); | 
| 1252 |  |  | 
| 1253 |  | snap_->atomData.force[atom1] += *(idat.f1); | 
| 1254 |  | snap_->atomData.force[atom2] -= *(idat.f1); |