| 53 |  | // surrounding cells (not just the 14 upper triangular blocks that | 
| 54 |  | // are used when the processor can see all pairs) | 
| 55 |  | #ifdef IS_MPI | 
| 56 | < | cellOffsets_.push_back( Vector3i(-1, 0, 0) ); | 
| 57 | < | cellOffsets_.push_back( Vector3i(-1,-1, 0) ); | 
| 58 | < | cellOffsets_.push_back( Vector3i( 0,-1, 0) ); | 
| 59 | < | cellOffsets_.push_back( Vector3i( 1,-1, 0) ); | 
| 60 | < | cellOffsets_.push_back( Vector3i( 0, 0,-1) ); | 
| 61 | < | cellOffsets_.push_back( Vector3i(-1, 0, 1) ); | 
| 56 | > | cellOffsets_.clear(); | 
| 57 |  | cellOffsets_.push_back( Vector3i(-1,-1,-1) ); | 
| 58 |  | cellOffsets_.push_back( Vector3i( 0,-1,-1) ); | 
| 59 | < | cellOffsets_.push_back( Vector3i( 1,-1,-1) ); | 
| 59 | > | cellOffsets_.push_back( Vector3i( 1,-1,-1) ); | 
| 60 | > | cellOffsets_.push_back( Vector3i(-1, 0,-1) ); | 
| 61 | > | cellOffsets_.push_back( Vector3i( 0, 0,-1) ); | 
| 62 |  | cellOffsets_.push_back( Vector3i( 1, 0,-1) ); | 
| 66 | – | cellOffsets_.push_back( Vector3i( 1, 1,-1) ); | 
| 67 | – | cellOffsets_.push_back( Vector3i( 0, 1,-1) ); | 
| 63 |  | cellOffsets_.push_back( Vector3i(-1, 1,-1) ); | 
| 64 | + | cellOffsets_.push_back( Vector3i( 0, 1,-1) ); | 
| 65 | + | cellOffsets_.push_back( Vector3i( 1, 1,-1) ); | 
| 66 | + | cellOffsets_.push_back( Vector3i(-1,-1, 0) ); | 
| 67 | + | cellOffsets_.push_back( Vector3i( 0,-1, 0) ); | 
| 68 | + | cellOffsets_.push_back( Vector3i( 1,-1, 0) ); | 
| 69 | + | cellOffsets_.push_back( Vector3i(-1, 0, 0) ); | 
| 70 | + | cellOffsets_.push_back( Vector3i( 0, 0, 0) ); | 
| 71 | + | cellOffsets_.push_back( Vector3i( 1, 0, 0) ); | 
| 72 | + | cellOffsets_.push_back( Vector3i(-1, 1, 0) ); | 
| 73 | + | cellOffsets_.push_back( Vector3i( 0, 1, 0) ); | 
| 74 | + | cellOffsets_.push_back( Vector3i( 1, 1, 0) ); | 
| 75 | + | cellOffsets_.push_back( Vector3i(-1,-1, 1) ); | 
| 76 | + | cellOffsets_.push_back( Vector3i( 0,-1, 1) ); | 
| 77 | + | cellOffsets_.push_back( Vector3i( 1,-1, 1) ); | 
| 78 | + | cellOffsets_.push_back( Vector3i(-1, 0, 1) ); | 
| 79 | + | cellOffsets_.push_back( Vector3i( 0, 0, 1) ); | 
| 80 | + | cellOffsets_.push_back( Vector3i( 1, 0, 1) ); | 
| 81 | + | cellOffsets_.push_back( Vector3i(-1, 1, 1) ); | 
| 82 | + | cellOffsets_.push_back( Vector3i( 0, 1, 1) ); | 
| 83 | + | cellOffsets_.push_back( Vector3i( 1, 1, 1) ); | 
| 84 |  | #endif | 
| 85 |  | } | 
| 86 |  |  | 
| 169 |  | AtomPlanIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | 
| 170 |  | AtomPlanIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | 
| 171 |  |  | 
| 157 | – | cerr << "Atoms in Local:\n"; | 
| 158 | – | for (int i = 0; i < AtomLocalToGlobal.size(); i++) { | 
| 159 | – | cerr << "i =\t" << i << "\t localAt =\t" << AtomLocalToGlobal[i] << "\n"; | 
| 160 | – | } | 
| 161 | – | cerr << "Atoms in Row:\n"; | 
| 162 | – | for (int i = 0; i < AtomRowToGlobal.size(); i++) { | 
| 163 | – | cerr << "i =\t" << i << "\t rowAt =\t" << AtomRowToGlobal[i] << "\n"; | 
| 164 | – | } | 
| 165 | – | cerr << "Atoms in Col:\n"; | 
| 166 | – | for (int i = 0; i < AtomColToGlobal.size(); i++) { | 
| 167 | – | cerr << "i =\t" << i << "\t colAt =\t" << AtomColToGlobal[i] << "\n"; | 
| 168 | – | } | 
| 169 | – |  | 
| 172 |  | cgRowToGlobal.resize(nGroupsInRow_); | 
| 173 |  | cgColToGlobal.resize(nGroupsInCol_); | 
| 174 |  | cgPlanIntRow->gather(cgLocalToGlobal, cgRowToGlobal); | 
| 175 |  | cgPlanIntColumn->gather(cgLocalToGlobal, cgColToGlobal); | 
| 176 |  |  | 
| 175 | – | cerr << "Gruops in Local:\n"; | 
| 176 | – | for (int i = 0; i < cgLocalToGlobal.size(); i++) { | 
| 177 | – | cerr << "i =\t" << i << "\t localCG =\t" << cgLocalToGlobal[i] << "\n"; | 
| 178 | – | } | 
| 179 | – | cerr << "Groups in Row:\n"; | 
| 180 | – | for (int i = 0; i < cgRowToGlobal.size(); i++) { | 
| 181 | – | cerr << "i =\t" << i << "\t rowCG =\t" << cgRowToGlobal[i] << "\n"; | 
| 182 | – | } | 
| 183 | – | cerr << "Groups in Col:\n"; | 
| 184 | – | for (int i = 0; i < cgColToGlobal.size(); i++) { | 
| 185 | – | cerr << "i =\t" << i << "\t colCG =\t" << cgColToGlobal[i] << "\n"; | 
| 186 | – | } | 
| 187 | – |  | 
| 188 | – |  | 
| 177 |  | massFactorsRow.resize(nAtomsInRow_); | 
| 178 |  | massFactorsCol.resize(nAtomsInCol_); | 
| 179 |  | AtomPlanRealRow->gather(massFactors, massFactorsRow); | 
| 564 |  |  | 
| 565 |  | // gather up the cutoff group positions | 
| 566 |  |  | 
| 579 | – | cerr  << "before gather\n"; | 
| 580 | – | for (int i = 0; i < snap_->cgData.position.size(); i++) { | 
| 581 | – | cerr << "cgpos = " << snap_->cgData.position[i] << "\n"; | 
| 582 | – | } | 
| 583 | – |  | 
| 567 |  | cgPlanVectorRow->gather(snap_->cgData.position, | 
| 568 |  | cgRowData.position); | 
| 569 |  |  | 
| 587 | – | cerr  << "after gather\n"; | 
| 588 | – | for (int i = 0; i < cgRowData.position.size(); i++) { | 
| 589 | – | cerr << "cgRpos = " << cgRowData.position[i] << "\n"; | 
| 590 | – | } | 
| 591 | – |  | 
| 570 |  | cgPlanVectorColumn->gather(snap_->cgData.position, | 
| 571 |  | cgColData.position); | 
| 594 | – | for (int i = 0; i < cgColData.position.size(); i++) { | 
| 595 | – | cerr << "cgCpos = " << cgColData.position[i] << "\n"; | 
| 596 | – | } | 
| 572 |  |  | 
| 573 |  |  | 
| 574 |  | // if needed, gather the atomic rotation matrices | 
| 706 |  |  | 
| 707 |  | for (int ii = 0;  ii < pot_temp.size(); ii++ ) | 
| 708 |  | pairwisePot += pot_temp[ii]; | 
| 709 | + |  | 
| 710 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { | 
| 711 | + | RealType ploc1 = pairwisePot[ii]; | 
| 712 | + | RealType ploc2 = 0.0; | 
| 713 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); | 
| 714 | + | pairwisePot[ii] = ploc2; | 
| 715 | + | } | 
| 716 | + |  | 
| 717 |  | #endif | 
| 718 |  |  | 
| 736 | – | cerr << "pairwisePot = " <<  pairwisePot << "\n"; | 
| 719 |  | } | 
| 720 |  |  | 
| 721 |  | int ForceMatrixDecomposition::getNAtomsInRow() { | 
| 750 |  |  | 
| 751 |  | #ifdef IS_MPI | 
| 752 |  | d = cgColData.position[cg2] - cgRowData.position[cg1]; | 
| 771 | – | cerr << "cg1 = " << cg1 << "\tcg1p = " << cgRowData.position[cg1] << "\n"; | 
| 772 | – | cerr << "cg2 = " << cg2 << "\tcg2p = " << cgColData.position[cg2] << "\n"; | 
| 753 |  | #else | 
| 754 |  | d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; | 
| 775 | – | cerr << "cg1 = " << cg1 << "\tcg1p = " << snap_->cgData.position[cg1] << "\n"; | 
| 776 | – | cerr << "cg2 = " << cg2 << "\tcg2p = " << snap_->cgData.position[cg2] << "\n"; | 
| 755 |  | #endif | 
| 756 |  |  | 
| 757 |  | snap_->wrapVector(d); | 
| 827 |  | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { | 
| 828 |  | int unique_id_1, unique_id_2; | 
| 829 |  |  | 
| 852 | – |  | 
| 853 | – | cerr << "sap with atom1, atom2 =\t" << atom1 << "\t" << atom2 << "\n"; | 
| 830 |  | #ifdef IS_MPI | 
| 831 |  | // in MPI, we have to look up the unique IDs for each atom | 
| 832 |  | unique_id_1 = AtomRowToGlobal[atom1]; | 
| 833 |  | unique_id_2 = AtomColToGlobal[atom2]; | 
| 834 |  |  | 
| 859 | – | cerr << "sap with uid1, uid2 =\t" << unique_id_1 << "\t" << unique_id_2 << "\n"; | 
| 835 |  | // this situation should only arise in MPI simulations | 
| 836 |  | if (unique_id_1 == unique_id_2) return true; | 
| 837 |  |  | 
| 1106 |  | // add this cutoff group to the list of groups in this cell; | 
| 1107 |  | cellListCol_[cellIndex].push_back(i); | 
| 1108 |  | } | 
| 1109 | + |  | 
| 1110 |  | #else | 
| 1111 |  | for (int i = 0; i < nGroups_; i++) { | 
| 1112 |  | rs = snap_->cgData.position[i]; | 
| 1132 |  | // add this cutoff group to the list of groups in this cell; | 
| 1133 |  | cellList_[cellIndex].push_back(i); | 
| 1134 |  | } | 
| 1135 | + |  | 
| 1136 |  | #endif | 
| 1137 |  |  | 
| 1138 |  | for (int m1z = 0; m1z < nCells_.z(); m1z++) { | 
| 1145 |  | os != cellOffsets_.end(); ++os) { | 
| 1146 |  |  | 
| 1147 |  | Vector3i m2v = m1v + (*os); | 
| 1148 | < |  | 
| 1148 | > |  | 
| 1149 | > |  | 
| 1150 |  | if (m2v.x() >= nCells_.x()) { | 
| 1151 |  | m2v.x() = 0; | 
| 1152 |  | } else if (m2v.x() < 0) { | 
| 1164 |  | } else if (m2v.z() < 0) { | 
| 1165 |  | m2v.z() = nCells_.z() - 1; | 
| 1166 |  | } | 
| 1167 | < |  | 
| 1167 | > |  | 
| 1168 |  | int m2 = Vlinear (m2v, nCells_); | 
| 1169 |  |  | 
| 1170 |  | #ifdef IS_MPI | 
| 1173 |  | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); | 
| 1174 |  | j2 != cellListCol_[m2].end(); ++j2) { | 
| 1175 |  |  | 
| 1176 | < | // In parallel, we need to visit *all* pairs of row & | 
| 1177 | < | // column indicies and will truncate later on. | 
| 1176 | > | // In parallel, we need to visit *all* pairs of row | 
| 1177 | > | // & column indicies and will divide labor in the | 
| 1178 | > | // force evaluation later. | 
| 1179 |  | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; | 
| 1180 |  | snap_->wrapVector(dr); | 
| 1181 |  | cuts = getGroupCutoffs( (*j1), (*j2) ); |