95 |
|
storageLayout_ = sman_->getStorageLayout(); |
96 |
|
ff_ = info_->getForceField(); |
97 |
|
nLocal_ = snap_->getNumberOfAtoms(); |
98 |
< |
|
98 |
> |
|
99 |
|
nGroups_ = info_->getNLocalCutoffGroups(); |
100 |
|
// gather the information for atomtype IDs (atids): |
101 |
|
idents = info_->getIdentArray(); |
109 |
|
PairList* oneTwo = info_->getOneTwoInteractions(); |
110 |
|
PairList* oneThree = info_->getOneThreeInteractions(); |
111 |
|
PairList* oneFour = info_->getOneFourInteractions(); |
112 |
< |
|
112 |
> |
|
113 |
> |
if (needVelocities_) |
114 |
> |
snap_->cgData.setStorageLayout(DataStorage::dslPosition | |
115 |
> |
DataStorage::dslVelocity); |
116 |
> |
else |
117 |
> |
snap_->cgData.setStorageLayout(DataStorage::dslPosition); |
118 |
> |
|
119 |
|
#ifdef IS_MPI |
120 |
|
|
121 |
|
MPI::Intracomm row = rowComm.getComm(); |
151 |
|
cgRowData.resize(nGroupsInRow_); |
152 |
|
cgRowData.setStorageLayout(DataStorage::dslPosition); |
153 |
|
cgColData.resize(nGroupsInCol_); |
154 |
< |
cgColData.setStorageLayout(DataStorage::dslPosition); |
155 |
< |
|
154 |
> |
if (needVelocities_) |
155 |
> |
// we only need column velocities if we need them. |
156 |
> |
cgColData.setStorageLayout(DataStorage::dslPosition | |
157 |
> |
DataStorage::dslVelocity); |
158 |
> |
else |
159 |
> |
cgColData.setStorageLayout(DataStorage::dslPosition); |
160 |
> |
|
161 |
|
identsRow.resize(nAtomsInRow_); |
162 |
|
identsCol.resize(nAtomsInCol_); |
163 |
|
|
536 |
|
atomColData.skippedCharge.end(), 0.0); |
537 |
|
} |
538 |
|
|
539 |
+ |
if (storageLayout_ & DataStorage::dslFlucQForce) { |
540 |
+ |
fill(atomRowData.flucQFrc.begin(), |
541 |
+ |
atomRowData.flucQFrc.end(), 0.0); |
542 |
+ |
fill(atomColData.flucQFrc.begin(), |
543 |
+ |
atomColData.flucQFrc.end(), 0.0); |
544 |
+ |
} |
545 |
+ |
|
546 |
|
if (storageLayout_ & DataStorage::dslElectricField) { |
547 |
|
fill(atomRowData.electricField.begin(), |
548 |
|
atomRowData.electricField.end(), V3Zero); |
549 |
|
fill(atomColData.electricField.begin(), |
550 |
|
atomColData.electricField.end(), V3Zero); |
551 |
|
} |
552 |
+ |
|
553 |
|
if (storageLayout_ & DataStorage::dslFlucQForce) { |
554 |
|
fill(atomRowData.flucQFrc.begin(), atomRowData.flucQFrc.end(), |
555 |
|
0.0); |
611 |
|
cgPlanVectorColumn->gather(snap_->cgData.position, |
612 |
|
cgColData.position); |
613 |
|
|
614 |
+ |
|
615 |
+ |
|
616 |
+ |
if (needVelocities_) { |
617 |
+ |
// gather up the atomic velocities |
618 |
+ |
AtomPlanVectorColumn->gather(snap_->atomData.velocity, |
619 |
+ |
atomColData.velocity); |
620 |
+ |
|
621 |
+ |
cgPlanVectorColumn->gather(snap_->cgData.velocity, |
622 |
+ |
cgColData.velocity); |
623 |
+ |
} |
624 |
+ |
|
625 |
|
|
626 |
|
// if needed, gather the atomic rotation matrices |
627 |
|
if (storageLayout_ & DataStorage::dslAmat) { |
788 |
|
|
789 |
|
for (int ii = 0; ii < pot_temp.size(); ii++ ) |
790 |
|
pairwisePot += pot_temp[ii]; |
791 |
< |
|
791 |
> |
|
792 |
> |
if (storageLayout_ & DataStorage::dslParticlePot) { |
793 |
> |
// This is the pairwise contribution to the particle pot. The |
794 |
> |
// embedding contribution is added in each of the low level |
795 |
> |
// non-bonded routines. In single processor, this is done in |
796 |
> |
// unpackInteractionData, not in collectData. |
797 |
> |
for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
798 |
> |
for (int i = 0; i < nLocal_; i++) { |
799 |
> |
// factor of two is because the total potential terms are divided |
800 |
> |
// by 2 in parallel due to row/ column scatter |
801 |
> |
snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii); |
802 |
> |
} |
803 |
> |
} |
804 |
> |
} |
805 |
> |
|
806 |
|
fill(pot_temp.begin(), pot_temp.end(), |
807 |
|
Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
808 |
|
|
810 |
|
|
811 |
|
for (int ii = 0; ii < pot_temp.size(); ii++ ) |
812 |
|
pairwisePot += pot_temp[ii]; |
813 |
+ |
|
814 |
+ |
if (storageLayout_ & DataStorage::dslParticlePot) { |
815 |
+ |
// This is the pairwise contribution to the particle pot. The |
816 |
+ |
// embedding contribution is added in each of the low level |
817 |
+ |
// non-bonded routines. In single processor, this is done in |
818 |
+ |
// unpackInteractionData, not in collectData. |
819 |
+ |
for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
820 |
+ |
for (int i = 0; i < nLocal_; i++) { |
821 |
+ |
// factor of two is because the total potential terms are divided |
822 |
+ |
// by 2 in parallel due to row/ column scatter |
823 |
+ |
snap_->atomData.particlePot[i] += 2.0 * pot_temp[i](ii); |
824 |
+ |
} |
825 |
+ |
} |
826 |
+ |
} |
827 |
|
|
828 |
+ |
if (storageLayout_ & DataStorage::dslParticlePot) { |
829 |
+ |
int npp = snap_->atomData.particlePot.size(); |
830 |
+ |
vector<RealType> ppot_temp(npp, 0.0); |
831 |
+ |
|
832 |
+ |
// This is the direct or embedding contribution to the particle |
833 |
+ |
// pot. |
834 |
+ |
|
835 |
+ |
AtomPlanRealRow->scatter(atomRowData.particlePot, ppot_temp); |
836 |
+ |
for (int i = 0; i < npp; i++) { |
837 |
+ |
snap_->atomData.particlePot[i] += ppot_temp[i]; |
838 |
+ |
} |
839 |
+ |
|
840 |
+ |
fill(ppot_temp.begin(), ppot_temp.end(), 0.0); |
841 |
+ |
|
842 |
+ |
AtomPlanRealColumn->scatter(atomColData.particlePot, ppot_temp); |
843 |
+ |
for (int i = 0; i < npp; i++) { |
844 |
+ |
snap_->atomData.particlePot[i] += ppot_temp[i]; |
845 |
+ |
} |
846 |
+ |
} |
847 |
+ |
|
848 |
|
for (int ii = 0; ii < N_INTERACTION_FAMILIES; ii++) { |
849 |
|
RealType ploc1 = pairwisePot[ii]; |
850 |
|
RealType ploc2 = 0.0; |
858 |
|
MPI::COMM_WORLD.Allreduce(&ploc1, &ploc2, 1, MPI::REALTYPE, MPI::SUM); |
859 |
|
embeddingPot[ii] = ploc2; |
860 |
|
} |
861 |
+ |
|
862 |
+ |
// Here be dragons. |
863 |
+ |
MPI::Intracomm col = colComm.getComm(); |
864 |
|
|
865 |
+ |
col.Allreduce(MPI::IN_PLACE, |
866 |
+ |
&snap_->frameData.conductiveHeatFlux[0], 3, |
867 |
+ |
MPI::REALTYPE, MPI::SUM); |
868 |
+ |
|
869 |
+ |
|
870 |
|
#endif |
871 |
|
|
872 |
|
} |
909 |
|
|
910 |
|
snap_->wrapVector(d); |
911 |
|
return d; |
912 |
+ |
} |
913 |
+ |
|
914 |
+ |
Vector3d ForceMatrixDecomposition::getGroupVelocityColumn(int cg2){ |
915 |
+ |
#ifdef IS_MPI |
916 |
+ |
return cgColData.velocity[cg2]; |
917 |
+ |
#else |
918 |
+ |
return snap_->cgData.velocity[cg2]; |
919 |
+ |
#endif |
920 |
|
} |
921 |
|
|
922 |
+ |
Vector3d ForceMatrixDecomposition::getAtomVelocityColumn(int atom2){ |
923 |
+ |
#ifdef IS_MPI |
924 |
+ |
return atomColData.velocity[atom2]; |
925 |
+ |
#else |
926 |
+ |
return snap_->atomData.velocity[atom2]; |
927 |
+ |
#endif |
928 |
+ |
} |
929 |
|
|
930 |
+ |
|
931 |
|
Vector3d ForceMatrixDecomposition::getAtomToGroupVectorRow(int atom1, int cg1){ |
932 |
|
|
933 |
|
Vector3d d; |
1107 |
|
if (storageLayout_ & DataStorage::dslSkippedCharge) { |
1108 |
|
idat.skippedCharge1 = &(atomRowData.skippedCharge[atom1]); |
1109 |
|
idat.skippedCharge2 = &(atomColData.skippedCharge[atom2]); |
1110 |
+ |
} |
1111 |
+ |
|
1112 |
+ |
if (storageLayout_ & DataStorage::dslFlucQPosition) { |
1113 |
+ |
idat.flucQ1 = &(atomRowData.flucQPos[atom1]); |
1114 |
+ |
idat.flucQ2 = &(atomColData.flucQPos[atom2]); |
1115 |
|
} |
1116 |
|
|
1117 |
|
#else |
1164 |
|
idat.skippedCharge1 = &(snap_->atomData.skippedCharge[atom1]); |
1165 |
|
idat.skippedCharge2 = &(snap_->atomData.skippedCharge[atom2]); |
1166 |
|
} |
1167 |
+ |
|
1168 |
+ |
if (storageLayout_ & DataStorage::dslFlucQPosition) { |
1169 |
+ |
idat.flucQ1 = &(snap_->atomData.flucQPos[atom1]); |
1170 |
+ |
idat.flucQ2 = &(snap_->atomData.flucQPos[atom2]); |
1171 |
+ |
} |
1172 |
+ |
|
1173 |
|
#endif |
1174 |
|
} |
1175 |
|
|
1182 |
|
atomRowData.force[atom1] += *(idat.f1); |
1183 |
|
atomColData.force[atom2] -= *(idat.f1); |
1184 |
|
|
1185 |
< |
// should particle pot be done here also? |
1185 |
> |
if (storageLayout_ & DataStorage::dslFlucQForce) { |
1186 |
> |
atomRowData.flucQFrc[atom1] += *(idat.dVdFQ1); |
1187 |
> |
atomColData.flucQFrc[atom2] += *(idat.dVdFQ2); |
1188 |
> |
} |
1189 |
> |
|
1190 |
> |
if (storageLayout_ & DataStorage::dslElectricField) { |
1191 |
> |
atomRowData.electricField[atom1] += *(idat.eField1); |
1192 |
> |
atomColData.electricField[atom2] += *(idat.eField2); |
1193 |
> |
} |
1194 |
> |
|
1195 |
|
#else |
1196 |
|
pairwisePot += *(idat.pot); |
1197 |
|
|
1199 |
|
snap_->atomData.force[atom2] -= *(idat.f1); |
1200 |
|
|
1201 |
|
if (idat.doParticlePot) { |
1202 |
< |
snap_->atomData.particlePot[atom1] += *(idat.vpair) * *(idat.sw); |
1203 |
< |
snap_->atomData.particlePot[atom2] -= *(idat.vpair) * *(idat.sw); |
1202 |
> |
// This is the pairwise contribution to the particle pot. The |
1203 |
> |
// embedding contribution is added in each of the low level |
1204 |
> |
// non-bonded routines. In parallel, this calculation is done |
1205 |
> |
// in collectData, not in unpackInteractionData. |
1206 |
> |
snap_->atomData.particlePot[atom1] += *(idat.vpair) * *(idat.sw); |
1207 |
> |
snap_->atomData.particlePot[atom2] += *(idat.vpair) * *(idat.sw); |
1208 |
|
} |
1209 |
< |
|
1209 |
> |
|
1210 |
> |
if (storageLayout_ & DataStorage::dslFlucQForce) { |
1211 |
> |
snap_->atomData.flucQFrc[atom1] += *(idat.dVdFQ1); |
1212 |
> |
snap_->atomData.flucQFrc[atom2] -= *(idat.dVdFQ2); |
1213 |
> |
} |
1214 |
> |
|
1215 |
> |
if (storageLayout_ & DataStorage::dslElectricField) { |
1216 |
> |
snap_->atomData.electricField[atom1] += *(idat.eField1); |
1217 |
> |
snap_->atomData.electricField[atom2] += *(idat.eField2); |
1218 |
> |
} |
1219 |
> |
|
1220 |
|
#endif |
1221 |
|
|
1222 |
|
} |