# | Line 35 | Line 35 | |
---|---|---|
35 | * | |
36 | * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | |
37 | * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | |
38 | < | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
38 | > | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 | * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). | |
40 | * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | |
41 | */ | |
# | Line 175 | Line 175 | namespace OpenMD { | |
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_); | |
# | Line 307 | Line 310 | namespace OpenMD { | |
310 | ||
311 | RealType tol = 1e-6; | |
312 | largestRcut_ = 0.0; | |
310 | – | RealType rc; |
313 | int atid; | |
314 | set<AtomType*> atypes = info_->getSimulatedAtomTypes(); | |
315 | ||
# | Line 392 | Line 394 | namespace OpenMD { | |
394 | } | |
395 | ||
396 | bool gTypeFound = false; | |
397 | < | for (int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
397 | > | for (unsigned int gt = 0; gt < gTypeCutoffs.size(); gt++) { |
398 | if (abs(groupCutoff[cg1] - gTypeCutoffs[gt]) < tol) { | |
399 | groupToGtype[cg1] = gt; | |
400 | gTypeFound = true; | |
# | Line 417 | Line 419 | namespace OpenMD { | |
419 | ||
420 | RealType tradRcut = groupMax; | |
421 | ||
422 | < | for (int i = 0; i < gTypeCutoffs.size(); i++) { |
423 | < | for (int j = 0; j < gTypeCutoffs.size(); j++) { |
422 | > | for (unsigned int i = 0; i < gTypeCutoffs.size(); i++) { |
423 | > | for (unsigned int j = 0; j < gTypeCutoffs.size(); j++) { |
424 | RealType thisRcut; | |
425 | switch(cutoffPolicy_) { | |
426 | case TRADITIONAL: | |
# | Line 461 | Line 463 | namespace OpenMD { | |
463 | } | |
464 | } | |
465 | ||
464 | – | |
466 | groupCutoffs ForceMatrixDecomposition::getGroupCutoffs(int cg1, int cg2) { | |
467 | int i, j; | |
468 | #ifdef IS_MPI | |
# | Line 475 | Line 476 | namespace OpenMD { | |
476 | } | |
477 | ||
478 | int ForceMatrixDecomposition::getTopologicalDistance(int atom1, int atom2) { | |
479 | < | for (int j = 0; j < toposForAtom[atom1].size(); j++) { |
479 | > | for (unsigned int j = 0; j < toposForAtom[atom1].size(); j++) { |
480 | if (toposForAtom[atom1][j] == atom2) | |
481 | return topoDist[atom1][j]; | |
482 | < | } |
482 | > | } |
483 | return 0; | |
484 | } | |
485 | ||
486 | void ForceMatrixDecomposition::zeroWorkArrays() { | |
487 | pairwisePot = 0.0; | |
488 | embeddingPot = 0.0; | |
489 | + | excludedPot = 0.0; |
490 | + | excludedSelfPot = 0.0; |
491 | ||
492 | #ifdef IS_MPI | |
493 | if (storageLayout_ & DataStorage::dslForce) { | |
# | Line 501 | Line 504 | namespace OpenMD { | |
504 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
505 | ||
506 | fill(pot_col.begin(), pot_col.end(), | |
507 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
508 | + | |
509 | + | fill(expot_row.begin(), expot_row.end(), |
510 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
511 | + | |
512 | + | fill(expot_col.begin(), expot_col.end(), |
513 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
514 | ||
515 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
# | Line 550 | Line 559 | namespace OpenMD { | |
559 | atomColData.electricField.end(), V3Zero); | |
560 | } | |
561 | ||
553 | – | if (storageLayout_ & DataStorage::dslFlucQForce) { |
554 | – | fill(atomRowData.flucQFrc.begin(), atomRowData.flucQFrc.end(), |
555 | – | 0.0); |
556 | – | fill(atomColData.flucQFrc.begin(), atomColData.flucQFrc.end(), |
557 | – | 0.0); |
558 | – | } |
559 | – | |
562 | #endif | |
563 | // even in parallel, we need to zero out the local arrays: | |
564 | ||
# | Line 630 | Line 632 | namespace OpenMD { | |
632 | AtomPlanMatrixColumn->gather(snap_->atomData.aMat, | |
633 | atomColData.aMat); | |
634 | } | |
635 | < | |
636 | < | // if needed, gather the atomic eletrostatic frames |
637 | < | if (storageLayout_ & DataStorage::dslElectroFrame) { |
638 | < | AtomPlanMatrixRow->gather(snap_->atomData.electroFrame, |
639 | < | atomRowData.electroFrame); |
640 | < | AtomPlanMatrixColumn->gather(snap_->atomData.electroFrame, |
641 | < | atomColData.electroFrame); |
635 | > | |
636 | > | // if needed, gather the atomic eletrostatic information |
637 | > | if (storageLayout_ & DataStorage::dslDipole) { |
638 | > | AtomPlanVectorRow->gather(snap_->atomData.dipole, |
639 | > | atomRowData.dipole); |
640 | > | AtomPlanVectorColumn->gather(snap_->atomData.dipole, |
641 | > | atomColData.dipole); |
642 | } | |
643 | ||
644 | + | if (storageLayout_ & DataStorage::dslQuadrupole) { |
645 | + | AtomPlanMatrixRow->gather(snap_->atomData.quadrupole, |
646 | + | atomRowData.quadrupole); |
647 | + | AtomPlanMatrixColumn->gather(snap_->atomData.quadrupole, |
648 | + | atomColData.quadrupole); |
649 | + | } |
650 | + | |
651 | // if needed, gather the atomic fluctuating charge values | |
652 | if (storageLayout_ & DataStorage::dslFlucQPosition) { | |
653 | AtomPlanRealRow->gather(snap_->atomData.flucQPos, | |
# | Line 670 | Line 679 | namespace OpenMD { | |
679 | snap_->atomData.density[i] += rho_tmp[i]; | |
680 | } | |
681 | ||
682 | + | // this isn't necessary if we don't have polarizable atoms, but |
683 | + | // we'll leave it here for now. |
684 | if (storageLayout_ & DataStorage::dslElectricField) { | |
685 | ||
686 | AtomPlanVectorRow->scatter(atomRowData.electricField, | |
# | Line 677 | Line 688 | namespace OpenMD { | |
688 | ||
689 | int n = snap_->atomData.electricField.size(); | |
690 | vector<Vector3d> field_tmp(n, V3Zero); | |
691 | < | AtomPlanVectorColumn->scatter(atomColData.electricField, field_tmp); |
691 | > | AtomPlanVectorColumn->scatter(atomColData.electricField, |
692 | > | field_tmp); |
693 | for (int i = 0; i < n; i++) | |
694 | snap_->atomData.electricField[i] += field_tmp[i]; | |
695 | } | |
# | Line 777 | Line 789 | namespace OpenMD { | |
789 | ||
790 | } | |
791 | ||
792 | + | if (storageLayout_ & DataStorage::dslElectricField) { |
793 | + | |
794 | + | int nef = snap_->atomData.electricField.size(); |
795 | + | vector<Vector3d> efield_tmp(nef, V3Zero); |
796 | + | |
797 | + | AtomPlanVectorRow->scatter(atomRowData.electricField, efield_tmp); |
798 | + | for (int i = 0; i < nef; i++) { |
799 | + | snap_->atomData.electricField[i] += efield_tmp[i]; |
800 | + | efield_tmp[i] = 0.0; |
801 | + | } |
802 | + | |
803 | + | AtomPlanVectorColumn->scatter(atomColData.electricField, efield_tmp); |
804 | + | for (int i = 0; i < nef; i++) |
805 | + | snap_->atomData.electricField[i] += efield_tmp[i]; |
806 | + | } |
807 | + | |
808 | + | |
809 | nLocal_ = snap_->getNumberOfAtoms(); | |
810 | ||
811 | vector<potVec> pot_temp(nLocal_, | |
812 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
813 | + | vector<potVec> expot_temp(nLocal_, |
814 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
815 | ||
816 | // scatter/gather pot_row into the members of my column | |
817 | ||
818 | AtomPlanPotRow->scatter(pot_row, pot_temp); | |
819 | + | AtomPlanPotRow->scatter(expot_row, expot_temp); |
820 | ||
821 | < | for (int ii = 0; ii < pot_temp.size(); ii++ ) |
821 | > | for (int ii = 0; ii < pot_temp.size(); ii++ ) |
822 | pairwisePot += pot_temp[ii]; | |
823 | + | |
824 | + | for (int ii = 0; ii < expot_temp.size(); ii++ ) |
825 | + | excludedPot += expot_temp[ii]; |
826 | ||
827 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
828 | // This is the pairwise contribution to the particle pot. The | |
# | Line 805 | Line 840 | namespace OpenMD { | |
840 | ||
841 | fill(pot_temp.begin(), pot_temp.end(), | |
842 | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); | |
843 | + | fill(expot_temp.begin(), expot_temp.end(), |
844 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
845 | ||
846 | AtomPlanPotColumn->scatter(pot_col, pot_temp); | |
847 | + | AtomPlanPotColumn->scatter(expot_col, expot_temp); |
848 | ||
849 | for (int ii = 0; ii < pot_temp.size(); ii++ ) | |
850 | pairwisePot += pot_temp[ii]; | |
851 | ||
852 | + | for (int ii = 0; ii < expot_temp.size(); ii++ ) |
853 | + | excludedPot += expot_temp[ii]; |
854 | + | |
855 | if (storageLayout_ & DataStorage::dslParticlePot) { | |
856 | // This is the pairwise contribution to the particle pot. The | |
857 | // embedding contribution is added in each of the low level | |
# | Line 853 | Line 894 | namespace OpenMD { | |
894 | } | |
895 | ||
896 | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { | |
897 | < | RealType ploc1 = embeddingPot[ii]; |
897 | > | RealType ploc1 = excludedPot[ii]; |
898 | RealType ploc2 = 0.0; | |
899 | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); | |
900 | < | embeddingPot[ii] = ploc2; |
900 | > | excludedPot[ii] = ploc2; |
901 | } | |
902 | < | |
902 | > | |
903 | // Here be dragons. | |
904 | MPI::Intracomm col = colComm.getComm(); | |
905 | ||
# | Line 871 | Line 912 | namespace OpenMD { | |
912 | ||
913 | } | |
914 | ||
915 | < | int ForceMatrixDecomposition::getNAtomsInRow() { |
915 | > | /** |
916 | > | * Collects information obtained during the post-pair (and embedding |
917 | > | * functional) loops onto local data structures. |
918 | > | */ |
919 | > | void ForceMatrixDecomposition::collectSelfData() { |
920 | > | snap_ = sman_->getCurrentSnapshot(); |
921 | > | storageLayout_ = sman_->getStorageLayout(); |
922 | > | |
923 | #ifdef IS_MPI | |
924 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
925 | + | RealType ploc1 = embeddingPot[ii]; |
926 | + | RealType ploc2 = 0.0; |
927 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
928 | + | embeddingPot[ii] = ploc2; |
929 | + | } |
930 | + | for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
931 | + | RealType ploc1 = excludedSelfPot[ii]; |
932 | + | RealType ploc2 = 0.0; |
933 | + | MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
934 | + | excludedSelfPot[ii] = ploc2; |
935 | + | } |
936 | + | #endif |
937 | + | |
938 | + | } |
939 | + | |
940 | + | |
941 | + | |
942 | + | int& ForceMatrixDecomposition::getNAtomsInRow() { |
943 | + | #ifdef IS_MPI |
944 | return nAtomsInRow_; | |
945 | #else | |
946 | return nLocal_; | |
# | Line 882 | Line 950 | namespace OpenMD { | |
950 | /** | |
951 | * returns the list of atoms belonging to this group. | |
952 | */ | |
953 | < | vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
953 | > | vector<int>& ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
954 | #ifdef IS_MPI | |
955 | return groupListRow_[cg1]; | |
956 | #else | |
# | Line 890 | Line 958 | namespace OpenMD { | |
958 | #endif | |
959 | } | |
960 | ||
961 | < | vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
961 | > | vector<int>& ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
962 | #ifdef IS_MPI | |
963 | return groupListCol_[cg2]; | |
964 | #else | |
# | Line 907 | Line 975 | namespace OpenMD { | |
975 | d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; | |
976 | #endif | |
977 | ||
978 | < | snap_->wrapVector(d); |
978 | > | if (usePeriodicBoundaryConditions_) { |
979 | > | snap_->wrapVector(d); |
980 | > | } |
981 | return d; | |
982 | } | |
983 | ||
984 | < | Vector3d ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ |
984 | > | Vector3d& ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ |
985 | #ifdef IS_MPI | |
986 | return cgColData.velocity[cg2]; | |
987 | #else | |
# | Line 919 | Line 989 | namespace OpenMD { | |
989 | #endif | |
990 | } | |
991 | ||
992 | < | Vector3d ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ |
992 | > | Vector3d& ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ |
993 | #ifdef IS_MPI | |
994 | return atomColData.velocity[atom2]; | |
995 | #else | |
# | Line 937 | Line 1007 | namespace OpenMD { | |
1007 | #else | |
1008 | d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1]; | |
1009 | #endif | |
1010 | < | |
1011 | < | snap_->wrapVector(d); |
1010 | > | if (usePeriodicBoundaryConditions_) { |
1011 | > | snap_->wrapVector(d); |
1012 | > | } |
1013 | return d; | |
1014 | } | |
1015 | ||
# | Line 950 | Line 1021 | namespace OpenMD { | |
1021 | #else | |
1022 | d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2]; | |
1023 | #endif | |
1024 | < | |
1025 | < | snap_->wrapVector(d); |
1024 | > | if (usePeriodicBoundaryConditions_) { |
1025 | > | snap_->wrapVector(d); |
1026 | > | } |
1027 | return d; | |
1028 | } | |
1029 | ||
1030 | < | RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
1030 | > | RealType& ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
1031 | #ifdef IS_MPI | |
1032 | return massFactorsRow[atom1]; | |
1033 | #else | |
# | Line 963 | Line 1035 | namespace OpenMD { | |
1035 | #endif | |
1036 | } | |
1037 | ||
1038 | < | RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
1038 | > | RealType& ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
1039 | #ifdef IS_MPI | |
1040 | return massFactorsCol[atom2]; | |
1041 | #else | |
# | Line 980 | Line 1052 | namespace OpenMD { | |
1052 | #else | |
1053 | d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1]; | |
1054 | #endif | |
1055 | < | |
1056 | < | snap_->wrapVector(d); |
1055 | > | if (usePeriodicBoundaryConditions_) { |
1056 | > | snap_->wrapVector(d); |
1057 | > | } |
1058 | return d; | |
1059 | } | |
1060 | ||
1061 | < | vector<int> ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
1061 | > | vector<int>& ForceMatrixDecomposition::getExcludesForAtom(int atom1) { |
1062 | return excludesForAtom[atom1]; | |
1063 | } | |
1064 | ||
# | Line 993 | Line 1066 | namespace OpenMD { | |
1066 | * We need to exclude some overcounted interactions that result from | |
1067 | * the parallel decomposition. | |
1068 | */ | |
1069 | < | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { |
1069 | > | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2, int cg1, int cg2) { |
1070 | int unique_id_1, unique_id_2; | |
1071 | ||
1072 | #ifdef IS_MPI | |
1073 | // in MPI, we have to look up the unique IDs for each atom | |
1074 | unique_id_1 = AtomRowToGlobal[atom1]; | |
1075 | unique_id_2 = AtomColToGlobal[atom2]; | |
1076 | + | // group1 = cgRowToGlobal[cg1]; |
1077 | + | // group2 = cgColToGlobal[cg2]; |
1078 | #else | |
1079 | unique_id_1 = AtomLocalToGlobal[atom1]; | |
1080 | unique_id_2 = AtomLocalToGlobal[atom2]; | |
1081 | + | int group1 = cgLocalToGlobal[cg1]; |
1082 | + | int group2 = cgLocalToGlobal[cg2]; |
1083 | #endif | |
1084 | ||
1085 | if (unique_id_1 == unique_id_2) return true; | |
# | Line 1014 | Line 1091 | namespace OpenMD { | |
1091 | } else { | |
1092 | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; | |
1093 | } | |
1094 | + | #endif |
1095 | + | |
1096 | + | #ifndef IS_MPI |
1097 | + | if (group1 == group2) { |
1098 | + | if (unique_id_1 < unique_id_2) return true; |
1099 | + | } |
1100 | #endif | |
1101 | ||
1102 | return false; | |
# | Line 1066 | Line 1149 | namespace OpenMD { | |
1149 | ||
1150 | #ifdef IS_MPI | |
1151 | idat.atypes = make_pair( atypesRow[atom1], atypesCol[atom2]); | |
1152 | + | idat.atid1 = identsRow[atom1]; |
1153 | + | idat.atid2 = identsCol[atom2]; |
1154 | //idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), | |
1155 | // ff_->getAtomType(identsCol[atom2]) ); | |
1156 | ||
# | Line 1074 | Line 1159 | namespace OpenMD { | |
1159 | idat.A2 = &(atomColData.aMat[atom2]); | |
1160 | } | |
1161 | ||
1077 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
1078 | – | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
1079 | – | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
1080 | – | } |
1081 | – | |
1162 | if (storageLayout_ & DataStorage::dslTorque) { | |
1163 | idat.t1 = &(atomRowData.torque[atom1]); | |
1164 | idat.t2 = &(atomColData.torque[atom2]); | |
1165 | } | |
1166 | ||
1167 | + | if (storageLayout_ & DataStorage::dslDipole) { |
1168 | + | idat.dipole1 = &(atomRowData.dipole[atom1]); |
1169 | + | idat.dipole2 = &(atomColData.dipole[atom2]); |
1170 | + | } |
1171 | + | |
1172 | + | if (storageLayout_ & DataStorage::dslQuadrupole) { |
1173 | + | idat.quadrupole1 = &(atomRowData.quadrupole[atom1]); |
1174 | + | idat.quadrupole2 = &(atomColData.quadrupole[atom2]); |
1175 | + | } |
1176 | + | |
1177 | if (storageLayout_ & DataStorage::dslDensity) { | |
1178 | idat.rho1 = &(atomRowData.density[atom1]); | |
1179 | idat.rho2 = &(atomColData.density[atom2]); | |
# | Line 1116 | Line 1206 | namespace OpenMD { | |
1206 | ||
1207 | #else | |
1208 | ||
1119 | – | |
1120 | – | // cerr << "atoms = " << atom1 << " " << atom2 << "\n"; |
1121 | – | // cerr << "pos1 = " << snap_->atomData.position[atom1] << "\n"; |
1122 | – | // cerr << "pos2 = " << snap_->atomData.position[atom2] << "\n"; |
1123 | – | |
1209 | idat.atypes = make_pair( atypesLocal[atom1], atypesLocal[atom2]); | |
1210 | < | //idat.atypes = make_pair( ff_->getAtomType(idents[atom1]), |
1211 | < | // ff_->getAtomType(idents[atom2]) ); |
1210 | > | idat.atid1 = idents[atom1]; |
1211 | > | idat.atid2 = idents[atom2]; |
1212 | ||
1213 | if (storageLayout_ & DataStorage::dslAmat) { | |
1214 | idat.A1 = &(snap_->atomData.aMat[atom1]); | |
1215 | idat.A2 = &(snap_->atomData.aMat[atom2]); | |
1131 | – | } |
1132 | – | |
1133 | – | if (storageLayout_ & DataStorage::dslElectroFrame) { |
1134 | – | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
1135 | – | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
1216 | } | |
1217 | ||
1218 | if (storageLayout_ & DataStorage::dslTorque) { | |
# | Line 1140 | Line 1220 | namespace OpenMD { | |
1220 | idat.t2 = &(snap_->atomData.torque[atom2]); | |
1221 | } | |
1222 | ||
1223 | + | if (storageLayout_ & DataStorage::dslDipole) { |
1224 | + | idat.dipole1 = &(snap_->atomData.dipole[atom1]); |
1225 | + | idat.dipole2 = &(snap_->atomData.dipole[atom2]); |
1226 | + | } |
1227 | + | |
1228 | + | if (storageLayout_ & DataStorage::dslQuadrupole) { |
1229 | + | idat.quadrupole1 = &(snap_->atomData.quadrupole[atom1]); |
1230 | + | idat.quadrupole2 = &(snap_->atomData.quadrupole[atom2]); |
1231 | + | } |
1232 | + | |
1233 | if (storageLayout_ & DataStorage::dslDensity) { | |
1234 | idat.rho1 = &(snap_->atomData.density[atom1]); | |
1235 | idat.rho2 = &(snap_->atomData.density[atom2]); | |
# | Line 1178 | Line 1268 | namespace OpenMD { | |
1268 | #ifdef IS_MPI | |
1269 | pot_row[atom1] += RealType(0.5) * *(idat.pot); | |
1270 | pot_col[atom2] += RealType(0.5) * *(idat.pot); | |
1271 | + | expot_row[atom1] += RealType(0.5) * *(idat.excludedPot); |
1272 | + | expot_col[atom2] += RealType(0.5) * *(idat.excludedPot); |
1273 | ||
1274 | atomRowData.force[atom1] += *(idat.f1); | |
1275 | atomColData.force[atom2] -= *(idat.f1); | |
# | Line 1194 | Line 1286 | namespace OpenMD { | |
1286 | ||
1287 | #else | |
1288 | pairwisePot += *(idat.pot); | |
1289 | + | excludedPot += *(idat.excludedPot); |
1290 | ||
1291 | snap_->atomData.force[atom1] += *(idat.f1); | |
1292 | snap_->atomData.force[atom2] -= *(idat.f1); | |
# | Line 1227 | Line 1320 | namespace OpenMD { | |
1320 | * first element of pair is row-indexed CutoffGroup | |
1321 | * second element of pair is column-indexed CutoffGroup | |
1322 | */ | |
1323 | < | vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() { |
1324 | < | |
1325 | < | vector<pair<int, int> > neighborList; |
1323 | > | void ForceMatrixDecomposition::buildNeighborList(vector<pair<int,int> >& neighborList) { |
1324 | > | |
1325 | > | neighborList.clear(); |
1326 | groupCutoffs cuts; | |
1327 | bool doAllPairs = false; | |
1328 | ||
1329 | + | RealType rList_ = (largestRcut_ + skinThickness_); |
1330 | + | Snapshot* snap_ = sman_->getCurrentSnapshot(); |
1331 | + | Mat3x3d box; |
1332 | + | Mat3x3d invBox; |
1333 | + | |
1334 | + | Vector3d rs, scaled, dr; |
1335 | + | Vector3i whichCell; |
1336 | + | int cellIndex; |
1337 | + | |
1338 | #ifdef IS_MPI | |
1339 | cellListRow_.clear(); | |
1340 | cellListCol_.clear(); | |
1341 | #else | |
1342 | cellList_.clear(); | |
1343 | #endif | |
1344 | < | |
1345 | < | RealType rList_ = (largestRcut_ + skinThickness_); |
1346 | < | RealType rl2 = rList_ * rList_; |
1347 | < | Snapshot* snap_ = sman_->getCurrentSnapshot(); |
1348 | < | Mat3x3d Hmat = snap_->getHmat(); |
1349 | < | Vector3d Hx = Hmat.getColumn(0); |
1350 | < | Vector3d Hy = Hmat.getColumn(1); |
1351 | < | Vector3d Hz = Hmat.getColumn(2); |
1352 | < | |
1353 | < | nCells_.x() = (int) ( Hx.length() )/ rList_; |
1354 | < | nCells_.y() = (int) ( Hy.length() )/ rList_; |
1355 | < | nCells_.z() = (int) ( Hz.length() )/ rList_; |
1356 | < | |
1344 | > | |
1345 | > | if (!usePeriodicBoundaryConditions_) { |
1346 | > | box = snap_->getBoundingBox(); |
1347 | > | invBox = snap_->getInvBoundingBox(); |
1348 | > | } else { |
1349 | > | box = snap_->getHmat(); |
1350 | > | invBox = snap_->getInvHmat(); |
1351 | > | } |
1352 | > | |
1353 | > | Vector3d boxX = box.getColumn(0); |
1354 | > | Vector3d boxY = box.getColumn(1); |
1355 | > | Vector3d boxZ = box.getColumn(2); |
1356 | > | |
1357 | > | nCells_.x() = (int) ( boxX.length() )/ rList_; |
1358 | > | nCells_.y() = (int) ( boxY.length() )/ rList_; |
1359 | > | nCells_.z() = (int) ( boxZ.length() )/ rList_; |
1360 | > | |
1361 | // handle small boxes where the cell offsets can end up repeating cells | |
1362 | ||
1363 | if (nCells_.x() < 3) doAllPairs = true; | |
1364 | if (nCells_.y() < 3) doAllPairs = true; | |
1365 | if (nCells_.z() < 3) doAllPairs = true; | |
1366 | < | |
1261 | < | Mat3x3d invHmat = snap_->getInvHmat(); |
1262 | < | Vector3d rs, scaled, dr; |
1263 | < | Vector3i whichCell; |
1264 | < | int cellIndex; |
1366 | > | |
1367 | int nCtot = nCells_.x() * nCells_.y() * nCells_.z(); | |
1368 | < | |
1368 | > | |
1369 | #ifdef IS_MPI | |
1370 | cellListRow_.resize(nCtot); | |
1371 | cellListCol_.resize(nCtot); | |
1372 | #else | |
1373 | cellList_.resize(nCtot); | |
1374 | #endif | |
1375 | < | |
1375 | > | |
1376 | if (!doAllPairs) { | |
1377 | #ifdef IS_MPI | |
1378 | < | |
1378 | > | |
1379 | for (int i = 0; i < nGroupsInRow_; i++) { | |
1380 | rs = cgRowData.position[i]; | |
1381 | ||
1382 | // scaled positions relative to the box vectors | |
1383 | < | scaled = invHmat * rs; |
1383 | > | scaled = invBox * rs; |
1384 | ||
1385 | // wrap the vector back into the unit box by subtracting integer box | |
1386 | // numbers | |
1387 | for (int j = 0; j < 3; j++) { | |
1388 | scaled[j] -= roundMe(scaled[j]); | |
1389 | scaled[j] += 0.5; | |
1390 | + | // Handle the special case when an object is exactly on the |
1391 | + | // boundary (a scaled coordinate of 1.0 is the same as |
1392 | + | // scaled coordinate of 0.0) |
1393 | + | if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
1394 | } | |
1395 | ||
1396 | // find xyz-indices of cell that cutoffGroup is in. | |
# | Line 1302 | Line 1408 | namespace OpenMD { | |
1408 | rs = cgColData.position[i]; | |
1409 | ||
1410 | // scaled positions relative to the box vectors | |
1411 | < | scaled = invHmat * rs; |
1411 | > | scaled = invBox * rs; |
1412 | ||
1413 | // wrap the vector back into the unit box by subtracting integer box | |
1414 | // numbers | |
1415 | for (int j = 0; j < 3; j++) { | |
1416 | scaled[j] -= roundMe(scaled[j]); | |
1417 | scaled[j] += 0.5; | |
1418 | + | // Handle the special case when an object is exactly on the |
1419 | + | // boundary (a scaled coordinate of 1.0 is the same as |
1420 | + | // scaled coordinate of 0.0) |
1421 | + | if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
1422 | } | |
1423 | ||
1424 | // find xyz-indices of cell that cutoffGroup is in. | |
# | Line 1322 | Line 1432 | namespace OpenMD { | |
1432 | // add this cutoff group to the list of groups in this cell; | |
1433 | cellListCol_[cellIndex].push_back(i); | |
1434 | } | |
1435 | < | |
1435 | > | |
1436 | #else | |
1437 | for (int i = 0; i < nGroups_; i++) { | |
1438 | rs = snap_->cgData.position[i]; | |
1439 | ||
1440 | // scaled positions relative to the box vectors | |
1441 | < | scaled = invHmat * rs; |
1441 | > | scaled = invBox * rs; |
1442 | ||
1443 | // wrap the vector back into the unit box by subtracting integer box | |
1444 | // numbers | |
1445 | for (int j = 0; j < 3; j++) { | |
1446 | scaled[j] -= roundMe(scaled[j]); | |
1447 | scaled[j] += 0.5; | |
1448 | + | // Handle the special case when an object is exactly on the |
1449 | + | // boundary (a scaled coordinate of 1.0 is the same as |
1450 | + | // scaled coordinate of 0.0) |
1451 | + | if (scaled[j] >= 1.0) scaled[j] -= 1.0; |
1452 | } | |
1453 | ||
1454 | // find xyz-indices of cell that cutoffGroup is in. | |
# | Line 1393 | Line 1507 | namespace OpenMD { | |
1507 | // & column indicies and will divide labor in the | |
1508 | // force evaluation later. | |
1509 | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; | |
1510 | < | snap_->wrapVector(dr); |
1510 | > | if (usePeriodicBoundaryConditions_) { |
1511 | > | snap_->wrapVector(dr); |
1512 | > | } |
1513 | cuts = getGroupCutoffs( (*j1), (*j2) ); | |
1514 | if (dr.lengthSquare() < cuts.third) { | |
1515 | neighborList.push_back(make_pair((*j1), (*j2))); | |
# | Line 1415 | Line 1531 | namespace OpenMD { | |
1531 | // allows atoms within a single cutoff group to | |
1532 | // interact with each other. | |
1533 | ||
1418 | – | |
1419 | – | |
1534 | if (m2 != m1 || (*j2) >= (*j1) ) { | |
1535 | ||
1536 | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; | |
1537 | < | snap_->wrapVector(dr); |
1537 | > | if (usePeriodicBoundaryConditions_) { |
1538 | > | snap_->wrapVector(dr); |
1539 | > | } |
1540 | cuts = getGroupCutoffs( (*j1), (*j2) ); | |
1541 | if (dr.lengthSquare() < cuts.third) { | |
1542 | neighborList.push_back(make_pair((*j1), (*j2))); | |
# | Line 1439 | Line 1555 | namespace OpenMD { | |
1555 | for (int j1 = 0; j1 < nGroupsInRow_; j1++) { | |
1556 | for (int j2 = 0; j2 < nGroupsInCol_; j2++) { | |
1557 | dr = cgColData.position[j2] - cgRowData.position[j1]; | |
1558 | < | snap_->wrapVector(dr); |
1558 | > | if (usePeriodicBoundaryConditions_) { |
1559 | > | snap_->wrapVector(dr); |
1560 | > | } |
1561 | cuts = getGroupCutoffs( j1, j2 ); | |
1562 | if (dr.lengthSquare() < cuts.third) { | |
1563 | neighborList.push_back(make_pair(j1, j2)); | |
# | Line 1452 | Line 1570 | namespace OpenMD { | |
1570 | // include self group interactions j2 == j1 | |
1571 | for (int j2 = j1; j2 < nGroups_; j2++) { | |
1572 | dr = snap_->cgData.position[j2] - snap_->cgData.position[j1]; | |
1573 | < | snap_->wrapVector(dr); |
1573 | > | if (usePeriodicBoundaryConditions_) { |
1574 | > | snap_->wrapVector(dr); |
1575 | > | } |
1576 | cuts = getGroupCutoffs( j1, j2 ); | |
1577 | if (dr.lengthSquare() < cuts.third) { | |
1578 | neighborList.push_back(make_pair(j1, j2)); | |
# | Line 1467 | Line 1587 | namespace OpenMD { | |
1587 | saved_CG_positions_.clear(); | |
1588 | for (int i = 0; i < nGroups_; i++) | |
1589 | saved_CG_positions_.push_back(snap_->cgData.position[i]); | |
1470 | – | |
1471 | – | return neighborList; |
1590 | } | |
1591 | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |