# | Line 61 | Line 61 | namespace OpenMD { | |
---|---|---|
61 | nGroups_ = info_->getNLocalCutoffGroups(); | |
62 | cerr << "in dId, nGroups = " << nGroups_ << "\n"; | |
63 | // gather the information for atomtype IDs (atids): | |
64 | < | identsLocal = info_->getIdentArray(); |
64 | > | idents = info_->getIdentArray(); |
65 | AtomLocalToGlobal = info_->getGlobalAtomIndices(); | |
66 | cgLocalToGlobal = info_->getGlobalGroupIndices(); | |
67 | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); | |
68 | < | vector<RealType> massFactorsLocal = info_->getMassFactors(); |
68 | > | massFactors = info_->getMassFactors(); |
69 | > | |
70 | PairList excludes = info_->getExcludedInteractions(); | |
71 | PairList oneTwo = info_->getOneTwoInteractions(); | |
72 | PairList oneThree = info_->getOneThreeInteractions(); | |
# | Line 108 | Line 109 | namespace OpenMD { | |
109 | identsRow.resize(nAtomsInRow_); | |
110 | identsCol.resize(nAtomsInCol_); | |
111 | ||
112 | < | AtomCommIntRow->gather(identsLocal, identsRow); |
113 | < | AtomCommIntColumn->gather(identsLocal, identsCol); |
112 | > | AtomCommIntRow->gather(idents, identsRow); |
113 | > | AtomCommIntColumn->gather(idents, identsCol); |
114 | ||
115 | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); | |
116 | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); | |
# | Line 117 | Line 118 | namespace OpenMD { | |
118 | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); | |
119 | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); | |
120 | ||
121 | < | AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); |
122 | < | AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); |
121 | > | AtomCommRealRow->gather(massFactors, massFactorsRow); |
122 | > | AtomCommRealColumn->gather(massFactors, massFactorsCol); |
123 | ||
124 | groupListRow_.clear(); | |
125 | groupListRow_.resize(nGroupsInRow_); | |
# | Line 232 | Line 233 | namespace OpenMD { | |
233 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
234 | vector<RealType> atypeCutoff; | |
235 | atypeCutoff.resize( atypes.size() ); | |
236 | < | |
236 | > | |
237 | for (set<AtomType*>::iterator at = atypes.begin(); | |
238 | at != atypes.end(); ++at){ | |
238 | – | rc = interactionMan_->getSuggestedCutoffRadius(*at); |
239 | atid = (*at)->getIdent(); | |
240 | < | atypeCutoff[atid] = rc; |
240 | > | |
241 | > | if (userChoseCutoff_) |
242 | > | atypeCutoff[atid] = userCutoff_; |
243 | > | else |
244 | > | atypeCutoff[atid] = interactionMan_->getSuggestedCutoffRadius(*at); |
245 | } | |
246 | ||
247 | vector<RealType> gTypeCutoffs; | |
# | Line 309 | Line 313 | namespace OpenMD { | |
313 | for (vector<int>::iterator ia = atomList.begin(); | |
314 | ia != atomList.end(); ++ia) { | |
315 | int atom1 = (*ia); | |
316 | < | atid = identsLocal[atom1]; |
316 | > | atid = idents[atom1]; |
317 | if (atypeCutoff[atid] > groupCutoff[cg1]) { | |
318 | groupCutoff[cg1] = atypeCutoff[atid]; | |
319 | } | |
# | Line 378 | Line 382 | namespace OpenMD { | |
382 | if (abs(gTypeCutoffMap[key].first - userCutoff_) > 0.0001) { | |
383 | sprintf(painCave.errMsg, | |
384 | "ForceMatrixDecomposition::createGtypeCutoffMap " | |
385 | < | "user-specified rCut does not match computed group Cutoff\n"); |
385 | > | "user-specified rCut (%lf) does not match computed group Cutoff\n", userCutoff_); |
386 | painCave.severity = OPENMD_ERROR; | |
387 | painCave.isFatal = 1; | |
388 | simError(); | |
# | Line 410 | Line 414 | namespace OpenMD { | |
414 | } | |
415 | ||
416 | void ForceMatrixDecomposition::zeroWorkArrays() { | |
417 | < | |
418 | < | for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { |
415 | < | longRangePot_[j] = 0.0; |
416 | < | } |
417 | > | pairwisePot = 0.0; |
418 | > | embeddingPot = 0.0; |
419 | ||
420 | #ifdef IS_MPI | |
421 | if (storageLayout_ & DataStorage::dslForce) { | |
# | Line 430 | Line 432 | namespace OpenMD { | |
432 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
433 | ||
434 | fill(pot_col.begin(), pot_col.end(), | |
435 | < | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
434 | < | |
435 | < | pot_local = Vector<RealType, N_INTERACTION_FAMILIES>(0.0); |
435 | > | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
436 | ||
437 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
438 | fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); | |
# | Line 606 | Line 606 | namespace OpenMD { | |
606 | AtomCommPotRow->scatter(pot_row, pot_temp); | |
607 | ||
608 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
609 | < | pot_local += pot_temp[ii]; |
609 | > | pairwisePot += pot_temp[ii]; |
610 | ||
611 | fill(pot_temp.begin(), pot_temp.end(), | |
612 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
# | Line 614 | Line 614 | namespace OpenMD { | |
614 | AtomCommPotColumn->scatter(pot_col, pot_temp); | |
615 | ||
616 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
617 | < | pot_local += pot_temp[ii]; |
618 | < | |
617 | > | pairwisePot += pot_temp[ii]; |
618 | #endif | |
619 | + | |
620 | } | |
621 | ||
622 | int ForceMatrixDecomposition::getNAtomsInRow() { | |
# | Line 691 | Line 691 | namespace OpenMD { | |
691 | #ifdef IS_MPI | |
692 | return massFactorsRow[atom1]; | |
693 | #else | |
694 | < | return massFactorsLocal[atom1]; |
694 | > | cerr << "mfs = " << massFactors.size() << " atom1 = " << atom1 << "\n"; |
695 | > | return massFactors[atom1]; |
696 | #endif | |
697 | } | |
698 | ||
# | Line 699 | Line 700 | namespace OpenMD { | |
700 | #ifdef IS_MPI | |
701 | return massFactorsCol[atom2]; | |
702 | #else | |
703 | < | return massFactorsLocal[atom2]; |
703 | > | return massFactors[atom2]; |
704 | #endif | |
705 | ||
706 | } | |
# | Line 754 | Line 755 | namespace OpenMD { | |
755 | for (vector<int>::iterator i = skipsForAtom[atom1].begin(); | |
756 | i != skipsForAtom[atom1].end(); ++i) { | |
757 | if ( (*i) == unique_id_2 ) return true; | |
758 | < | } |
758 | > | } |
759 | ||
760 | + | return false; |
761 | } | |
762 | ||
763 | ||
# | Line 776 | Line 778 | namespace OpenMD { | |
778 | } | |
779 | ||
780 | // filling interaction blocks with pointers | |
781 | < | InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) { |
782 | < | InteractionData idat; |
781 | < | |
781 | > | void ForceMatrixDecomposition::fillInteractionData(InteractionData &idat, |
782 | > | int atom1, int atom2) { |
783 | #ifdef IS_MPI | |
784 | ||
785 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
786 | ff_->getAtomType(identsCol[atom2]) ); | |
786 | – | |
787 | ||
788 | if (storageLayout_ & DataStorage::dslAmat) { | |
789 | idat.A1 = &(atomRowData.aMat[atom1]); | |
# | Line 822 | Line 822 | namespace OpenMD { | |
822 | ||
823 | #else | |
824 | ||
825 | < | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
826 | < | ff_->getAtomType(identsLocal[atom2]) ); |
825 | > | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
826 | > | ff_->getAtomType(idents[atom2]) ); |
827 | ||
828 | if (storageLayout_ & DataStorage::dslAmat) { | |
829 | idat.A1 = &(snap_->atomData.aMat[atom1]); | |
# | Line 840 | Line 840 | namespace OpenMD { | |
840 | idat.t2 = &(snap_->atomData.torque[atom2]); | |
841 | } | |
842 | ||
843 | < | if (storageLayout_ & DataStorage::dslDensity) { |
843 | > | if (storageLayout_ & DataStorage::dslDensity) { |
844 | idat.rho1 = &(snap_->atomData.density[atom1]); | |
845 | idat.rho2 = &(snap_->atomData.density[atom2]); | |
846 | } | |
# | Line 861 | Line 861 | namespace OpenMD { | |
861 | } | |
862 | ||
863 | #endif | |
864 | – | return idat; |
864 | } | |
865 | ||
866 | ||
867 | < | void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { |
867 | > | void ForceMatrixDecomposition::unpackInteractionData(InteractionData &idat, int atom1, int atom2) { |
868 | #ifdef IS_MPI | |
869 | pot_row[atom1] += 0.5 * *(idat.pot); | |
870 | pot_col[atom2] += 0.5 * *(idat.pot); | |
# | Line 873 | Line 872 | namespace OpenMD { | |
872 | atomRowData.force[atom1] += *(idat.f1); | |
873 | atomColData.force[atom2] -= *(idat.f1); | |
874 | #else | |
875 | < | longRangePot_ += *(idat.pot); |
876 | < | |
875 | > | pairwisePot += *(idat.pot); |
876 | > | |
877 | snap_->atomData.force[atom1] += *(idat.f1); | |
878 | snap_->atomData.force[atom2] -= *(idat.f1); | |
879 | #endif | |
# | Line 882 | Line 881 | namespace OpenMD { | |
881 | } | |
882 | ||
883 | ||
884 | < | InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ |
885 | < | |
887 | < | InteractionData idat; |
884 | > | void ForceMatrixDecomposition::fillSkipData(InteractionData &idat, |
885 | > | int atom1, int atom2) { |
886 | #ifdef IS_MPI | |
887 | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
888 | ff_->getAtomType(identsCol[atom2]) ); | |
# | Line 893 | Line 891 | namespace OpenMD { | |
891 | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); | |
892 | idat.eFrame2 = &(atomColData.electroFrame[atom2]); | |
893 | } | |
894 | + | |
895 | if (storageLayout_ & DataStorage::dslTorque) { | |
896 | idat.t1 = &(atomRowData.torque[atom1]); | |
897 | idat.t2 = &(atomColData.torque[atom2]); | |
898 | } | |
899 | + | |
900 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
901 | + | idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
902 | + | idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
903 | + | } |
904 | #else | |
905 | < | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
906 | < | ff_->getAtomType(identsLocal[atom2]) ); |
905 | > | idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
906 | > | ff_->getAtomType(idents[atom2]) ); |
907 | ||
908 | if (storageLayout_ & DataStorage::dslElectroFrame) { | |
909 | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); | |
910 | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); | |
911 | } | |
912 | + | |
913 | if (storageLayout_ & DataStorage::dslTorque) { | |
914 | idat.t1 = &(snap_->atomData.torque[atom1]); | |
915 | idat.t2 = &(snap_->atomData.torque[atom2]); | |
916 | } | |
917 | + | |
918 | + | if (storageLayout_ & DataStorage::dslSkippedCharge) { |
919 | + | idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
920 | + | idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
921 | + | } |
922 | #endif | |
923 | } | |
924 | ||
925 | + | |
926 | + | void ForceMatrixDecomposition::unpackSkipData(InteractionData &idat, int atom1, int atom2) { |
927 | + | #ifdef IS_MPI |
928 | + | pot_row[atom1] += 0.5 * *(idat.pot); |
929 | + | pot_col[atom2] += 0.5 * *(idat.pot); |
930 | + | #else |
931 | + | pairwisePot += *(idat.pot); |
932 | + | #endif |
933 | + | |
934 | + | } |
935 | + | |
936 | + | |
937 | /* | |
938 | * buildNeighborList | |
939 | * | |
# | Line 957 | Line 979 | namespace OpenMD { | |
979 | #ifdef IS_MPI | |
980 | for (int i = 0; i < nGroupsInRow_; i++) { | |
981 | rs = cgRowData.position[i]; | |
982 | + | |
983 | // scaled positions relative to the box vectors | |
984 | scaled = invHmat * rs; | |
985 | + | |
986 | // wrap the vector back into the unit box by subtracting integer box | |
987 | // numbers | |
988 | < | for (int j = 0; j < 3; j++) |
988 | > | for (int j = 0; j < 3; j++) { |
989 | scaled[j] -= roundMe(scaled[j]); | |
990 | + | scaled[j] += 0.5; |
991 | + | } |
992 | ||
993 | // find xyz-indices of cell that cutoffGroup is in. | |
994 | whichCell.x() = nCells_.x() * scaled.x(); | |
# | Line 971 | Line 997 | namespace OpenMD { | |
997 | ||
998 | // find single index of this cell: | |
999 | cellIndex = Vlinear(whichCell, nCells_); | |
1000 | + | |
1001 | // add this cutoff group to the list of groups in this cell; | |
1002 | cellListRow_[cellIndex].push_back(i); | |
1003 | } | |
1004 | ||
1005 | for (int i = 0; i < nGroupsInCol_; i++) { | |
1006 | rs = cgColData.position[i]; | |
1007 | + | |
1008 | // scaled positions relative to the box vectors | |
1009 | scaled = invHmat * rs; | |
1010 | + | |
1011 | // wrap the vector back into the unit box by subtracting integer box | |
1012 | // numbers | |
1013 | < | for (int j = 0; j < 3; j++) |
1013 | > | for (int j = 0; j < 3; j++) { |
1014 | scaled[j] -= roundMe(scaled[j]); | |
1015 | + | scaled[j] += 0.5; |
1016 | + | } |
1017 | ||
1018 | // find xyz-indices of cell that cutoffGroup is in. | |
1019 | whichCell.x() = nCells_.x() * scaled.x(); | |
# | Line 991 | Line 1022 | namespace OpenMD { | |
1022 | ||
1023 | // find single index of this cell: | |
1024 | cellIndex = Vlinear(whichCell, nCells_); | |
1025 | + | |
1026 | // add this cutoff group to the list of groups in this cell; | |
1027 | cellListCol_[cellIndex].push_back(i); | |
1028 | } | |
1029 | #else | |
1030 | for (int i = 0; i < nGroups_; i++) { | |
1031 | rs = snap_->cgData.position[i]; | |
1032 | + | |
1033 | // scaled positions relative to the box vectors | |
1034 | scaled = invHmat * rs; | |
1035 | + | |
1036 | // wrap the vector back into the unit box by subtracting integer box | |
1037 | // numbers | |
1038 | < | for (int j = 0; j < 3; j++) |
1038 | > | for (int j = 0; j < 3; j++) { |
1039 | scaled[j] -= roundMe(scaled[j]); | |
1040 | + | scaled[j] += 0.5; |
1041 | + | } |
1042 | ||
1043 | // find xyz-indices of cell that cutoffGroup is in. | |
1044 | whichCell.x() = nCells_.x() * scaled.x(); | |
# | Line 1010 | Line 1046 | namespace OpenMD { | |
1046 | whichCell.z() = nCells_.z() * scaled.z(); | |
1047 | ||
1048 | // find single index of this cell: | |
1049 | < | cellIndex = Vlinear(whichCell, nCells_); |
1049 | > | cellIndex = Vlinear(whichCell, nCells_); |
1050 | > | |
1051 | // add this cutoff group to the list of groups in this cell; | |
1052 | cellList_[cellIndex].push_back(i); | |
1053 | } | |
# | Line 1068 | Line 1105 | namespace OpenMD { | |
1105 | } | |
1106 | } | |
1107 | #else | |
1108 | + | |
1109 | for (vector<int>::iterator j1 = cellList_[m1].begin(); | |
1110 | j1 != cellList_[m1].end(); ++j1) { | |
1111 | for (vector<int>::iterator j2 = cellList_[m2].begin(); | |
1112 | j2 != cellList_[m2].end(); ++j2) { | |
1113 | < | |
1113 | > | |
1114 | // Always do this if we're in different cells or if | |
1115 | // we're in the same cell and the global index of the | |
1116 | // j2 cutoff group is less than the j1 cutoff group | |
# | Line 1092 | Line 1130 | namespace OpenMD { | |
1130 | } | |
1131 | } | |
1132 | } | |
1133 | < | |
1133 | > | |
1134 | // save the local cutoff group positions for the check that is | |
1135 | // done on each loop: | |
1136 | saved_CG_positions_.clear(); | |
1137 | for (int i = 0; i < nGroups_; i++) | |
1138 | saved_CG_positions_.push_back(snap_->cgData.position[i]); | |
1139 | < | |
1139 | > | |
1140 | return neighborList; | |
1141 | } | |
1142 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |