# | Line 1 | Line 1 | |
---|---|---|
1 | < | /** |
2 | < | * @file ForceDecomposition.cpp |
3 | < | * @author Charles Vardeman <cvardema.at.nd.edu> |
4 | < | * @date 08/18/2010 |
5 | < | * @time 11:56am |
6 | < | * @version 1.0 |
1 | > | /* |
2 | > | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 | * | |
8 | – | * @section LICENSE |
9 | – | * Copyright (c) 2010 The University of Notre Dame. All Rights Reserved. |
10 | – | * |
4 | * The University of Notre Dame grants you ("Licensee") a | |
5 | * non-exclusive, royalty free, license to use, modify and | |
6 | * redistribute this software in source and binary code form, provided | |
# | Line 45 | Line 38 | |
38 | * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | |
39 | * [4] Vardeman & Gezelter, in progress (2009). | |
40 | */ | |
41 | + | #include "parallel/ForceMatrixDecomposition.hpp" |
42 | + | #include "math/SquareMatrix3.hpp" |
43 | + | #include "nonbonded/NonBondedInteraction.hpp" |
44 | + | #include "brains/SnapshotManager.hpp" |
45 | + | #include "brains/PairList.hpp" |
46 | ||
47 | + | using namespace std; |
48 | + | namespace OpenMD { |
49 | ||
50 | + | /** |
51 | + | * distributeInitialData is essentially a copy of the older fortran |
52 | + | * SimulationSetup |
53 | + | */ |
54 | + | |
55 | + | void ForceMatrixDecomposition::distributeInitialData() { |
56 | + | snap_ = sman_->getCurrentSnapshot(); |
57 | + | storageLayout_ = sman_->getStorageLayout(); |
58 | + | ff_ = info_->getForceField(); |
59 | + | nLocal_ = snap_->getNumberOfAtoms(); |
60 | + | nGroups_ = snap_->getNumberOfCutoffGroups(); |
61 | ||
62 | < | /* -*- c++ -*- */ |
63 | < | #include "config.h" |
64 | < | #include <stdlib.h> |
62 | > | // gather the information for atomtype IDs (atids): |
63 | > | identsLocal = info_->getIdentArray(); |
64 | > | AtomLocalToGlobal = info_->getGlobalAtomIndices(); |
65 | > | cgLocalToGlobal = info_->getGlobalGroupIndices(); |
66 | > | vector<int> globalGroupMembership = info_->getGlobalGroupMembership(); |
67 | > | vector<RealType> massFactorsLocal = info_->getMassFactors(); |
68 | > | PairList excludes = info_->getExcludedInteractions(); |
69 | > | PairList oneTwo = info_->getOneTwoInteractions(); |
70 | > | PairList oneThree = info_->getOneThreeInteractions(); |
71 | > | PairList oneFour = info_->getOneFourInteractions(); |
72 | > | |
73 | #ifdef IS_MPI | |
74 | < | #include <mpi.h> |
75 | < | #endif |
74 | > | |
75 | > | AtomCommIntRow = new Communicator<Row,int>(nLocal_); |
76 | > | AtomCommRealRow = new Communicator<Row,RealType>(nLocal_); |
77 | > | AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_); |
78 | > | AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_); |
79 | > | AtomCommPotRow = new Communicator<Row,potVec>(nLocal_); |
80 | ||
81 | < | #include <iostream> |
82 | < | #include <vector> |
83 | < | #include <algorithm> |
84 | < | #include <cmath> |
85 | < | #include "parallel/ForceDecomposition.hpp" |
81 | > | AtomCommIntColumn = new Communicator<Column,int>(nLocal_); |
82 | > | AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_); |
83 | > | AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_); |
84 | > | AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_); |
85 | > | AtomCommPotColumn = new Communicator<Column,potVec>(nLocal_); |
86 | ||
87 | + | cgCommIntRow = new Communicator<Row,int>(nGroups_); |
88 | + | cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_); |
89 | + | cgCommIntColumn = new Communicator<Column,int>(nGroups_); |
90 | + | cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_); |
91 | ||
92 | < | using namespace std; |
93 | < | using namespace OpenMD; |
92 | > | nAtomsInRow_ = AtomCommIntRow->getSize(); |
93 | > | nAtomsInCol_ = AtomCommIntColumn->getSize(); |
94 | > | nGroupsInRow_ = cgCommIntRow->getSize(); |
95 | > | nGroupsInCol_ = cgCommIntColumn->getSize(); |
96 | ||
97 | < | //__static |
98 | < | #ifdef IS_MPI |
99 | < | static vector<MPI:Comm> communictors; |
100 | < | #endif |
97 | > | // Modify the data storage objects with the correct layouts and sizes: |
98 | > | atomRowData.resize(nAtomsInRow_); |
99 | > | atomRowData.setStorageLayout(storageLayout_); |
100 | > | atomColData.resize(nAtomsInCol_); |
101 | > | atomColData.setStorageLayout(storageLayout_); |
102 | > | cgRowData.resize(nGroupsInRow_); |
103 | > | cgRowData.setStorageLayout(DataStorage::dslPosition); |
104 | > | cgColData.resize(nGroupsInCol_); |
105 | > | cgColData.setStorageLayout(DataStorage::dslPosition); |
106 | > | |
107 | > | identsRow.reserve(nAtomsInRow_); |
108 | > | identsCol.reserve(nAtomsInCol_); |
109 | > | |
110 | > | AtomCommIntRow->gather(identsLocal, identsRow); |
111 | > | AtomCommIntColumn->gather(identsLocal, identsCol); |
112 | > | |
113 | > | AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal); |
114 | > | AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal); |
115 | > | |
116 | > | cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal); |
117 | > | cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal); |
118 | ||
119 | < | //____ MPITypeTraits |
120 | < | template<typename T> |
75 | < | struct MPITypeTraits; |
119 | > | AtomCommRealRow->gather(massFactorsLocal, massFactorsRow); |
120 | > | AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol); |
121 | ||
122 | < | #ifdef IS_MPI |
123 | < | template<> |
124 | < | struct MPITypeTraits<RealType> { |
125 | < | static const MPI::Datatype datatype; |
126 | < | }; |
127 | < | const MPI_Datatype MPITypeTraits<RealType>::datatype = MY_MPI_REAL; |
122 | > | groupListRow_.clear(); |
123 | > | groupListRow_.reserve(nGroupsInRow_); |
124 | > | for (int i = 0; i < nGroupsInRow_; i++) { |
125 | > | int gid = cgRowToGlobal[i]; |
126 | > | for (int j = 0; j < nAtomsInRow_; j++) { |
127 | > | int aid = AtomRowToGlobal[j]; |
128 | > | if (globalGroupMembership[aid] == gid) |
129 | > | groupListRow_[i].push_back(j); |
130 | > | } |
131 | > | } |
132 | ||
133 | < | template<> |
134 | < | struct MPITypeTraits<int> { |
135 | < | static const MPI::Datatype datatype; |
136 | < | }; |
137 | < | const MPI::Datatype MPITypeTraits<int>::datatype = MPI_INT; |
138 | < | #endif |
133 | > | groupListCol_.clear(); |
134 | > | groupListCol_.reserve(nGroupsInCol_); |
135 | > | for (int i = 0; i < nGroupsInCol_; i++) { |
136 | > | int gid = cgColToGlobal[i]; |
137 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
138 | > | int aid = AtomColToGlobal[j]; |
139 | > | if (globalGroupMembership[aid] == gid) |
140 | > | groupListCol_[i].push_back(j); |
141 | > | } |
142 | > | } |
143 | ||
144 | < | /** |
145 | < | * Constructor for ForceDecomposition Parallel Decomposition Method |
146 | < | * Will try to construct a symmetric grid of processors. Ideally, the |
147 | < | * number of processors will be a square ex: 4, 9, 16, 25. |
148 | < | * |
149 | < | */ |
144 | > | skipsForRowAtom.clear(); |
145 | > | skipsForRowAtom.reserve(nAtomsInRow_); |
146 | > | for (int i = 0; i < nAtomsInRow_; i++) { |
147 | > | int iglob = AtomRowToGlobal[i]; |
148 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
149 | > | int jglob = AtomColToGlobal[j]; |
150 | > | if (excludes.hasPair(iglob, jglob)) |
151 | > | skipsForRowAtom[i].push_back(j); |
152 | > | } |
153 | > | } |
154 | ||
155 | < | ForceDecomposition::ForceDecomposition() { |
155 | > | toposForRowAtom.clear(); |
156 | > | toposForRowAtom.reserve(nAtomsInRow_); |
157 | > | for (int i = 0; i < nAtomsInRow_; i++) { |
158 | > | int iglob = AtomRowToGlobal[i]; |
159 | > | int nTopos = 0; |
160 | > | for (int j = 0; j < nAtomsInCol_; j++) { |
161 | > | int jglob = AtomColToGlobal[j]; |
162 | > | if (oneTwo.hasPair(iglob, jglob)) { |
163 | > | toposForRowAtom[i].push_back(j); |
164 | > | topoDistRow[i][nTopos] = 1; |
165 | > | nTopos++; |
166 | > | } |
167 | > | if (oneThree.hasPair(iglob, jglob)) { |
168 | > | toposForRowAtom[i].push_back(j); |
169 | > | topoDistRow[i][nTopos] = 2; |
170 | > | nTopos++; |
171 | > | } |
172 | > | if (oneFour.hasPair(iglob, jglob)) { |
173 | > | toposForRowAtom[i].push_back(j); |
174 | > | topoDistRow[i][nTopos] = 3; |
175 | > | nTopos++; |
176 | > | } |
177 | > | } |
178 | > | } |
179 | ||
100 | – | #ifdef IS_MPI |
101 | – | int nProcs = MPI::COMM_WORLD.Get_size(); |
102 | – | int worldRank = MPI::COMM_WORLD.Get_rank(); |
180 | #endif | |
181 | ||
182 | < | // First time through, construct column stride. |
183 | < | if (communicators.size() == 0) |
184 | < | { |
185 | < | int nColumnsMax = (int) round(sqrt((float) nProcs)); |
186 | < | for (int i = 0; i < nProcs; ++i) |
187 | < | { |
188 | < | if (nProcs%i==0) nColumns=i; |
182 | > | groupList_.clear(); |
183 | > | groupList_.reserve(nGroups_); |
184 | > | for (int i = 0; i < nGroups_; i++) { |
185 | > | int gid = cgLocalToGlobal[i]; |
186 | > | for (int j = 0; j < nLocal_; j++) { |
187 | > | int aid = AtomLocalToGlobal[j]; |
188 | > | if (globalGroupMembership[aid] == gid) |
189 | > | groupList_[i].push_back(j); |
190 | > | } |
191 | } | |
192 | ||
193 | < | int nRows = nProcs/nColumns; |
194 | < | myRank_ = (int) worldRank%nColumns; |
116 | < | } |
117 | < | else |
118 | < | { |
119 | < | myRank_ = myRank/nColumns; |
120 | < | } |
121 | < | MPI::Comm newComm = MPI:COMM_WORLD.Split(myRank_,0); |
122 | < | |
123 | < | isColumn_ = false; |
124 | < | |
125 | < | } |
193 | > | skipsForLocalAtom.clear(); |
194 | > | skipsForLocalAtom.reserve(nLocal_); |
195 | ||
196 | < | ForceDecomposition::gather(sendbuf, receivebuf){ |
197 | < | communicators(myIndex_).Allgatherv(); |
198 | < | } |
196 | > | for (int i = 0; i < nLocal_; i++) { |
197 | > | int iglob = AtomLocalToGlobal[i]; |
198 | > | for (int j = 0; j < nLocal_; j++) { |
199 | > | int jglob = AtomLocalToGlobal[j]; |
200 | > | if (excludes.hasPair(iglob, jglob)) |
201 | > | skipsForLocalAtom[i].push_back(j); |
202 | > | } |
203 | > | } |
204 | ||
205 | < | |
206 | < | |
207 | < | ForceDecomposition::scatter(sbuffer, rbuffer){ |
208 | < | communicators(myIndex_).Reduce_scatter(sbuffer, recevbuf. recvcounts, MPI::DOUBLE, MPI::SUM); |
209 | < | } |
205 | > | toposForLocalAtom.clear(); |
206 | > | toposForLocalAtom.reserve(nLocal_); |
207 | > | for (int i = 0; i < nLocal_; i++) { |
208 | > | int iglob = AtomLocalToGlobal[i]; |
209 | > | int nTopos = 0; |
210 | > | for (int j = 0; j < nLocal_; j++) { |
211 | > | int jglob = AtomLocalToGlobal[j]; |
212 | > | if (oneTwo.hasPair(iglob, jglob)) { |
213 | > | toposForLocalAtom[i].push_back(j); |
214 | > | topoDistLocal[i][nTopos] = 1; |
215 | > | nTopos++; |
216 | > | } |
217 | > | if (oneThree.hasPair(iglob, jglob)) { |
218 | > | toposForLocalAtom[i].push_back(j); |
219 | > | topoDistLocal[i][nTopos] = 2; |
220 | > | nTopos++; |
221 | > | } |
222 | > | if (oneFour.hasPair(iglob, jglob)) { |
223 | > | toposForLocalAtom[i].push_back(j); |
224 | > | topoDistLocal[i][nTopos] = 3; |
225 | > | nTopos++; |
226 | > | } |
227 | > | } |
228 | > | } |
229 | > | } |
230 | > | |
231 | > | void ForceMatrixDecomposition::zeroWorkArrays() { |
232 | ||
233 | + | for (int j = 0; j < N_INTERACTION_FAMILIES; j++) { |
234 | + | longRangePot_[j] = 0.0; |
235 | + | } |
236 | ||
237 | + | #ifdef IS_MPI |
238 | + | if (storageLayout_ & DataStorage::dslForce) { |
239 | + | fill(atomRowData.force.begin(), atomRowData.force.end(), V3Zero); |
240 | + | fill(atomColData.force.begin(), atomColData.force.end(), V3Zero); |
241 | + | } |
242 | + | |
243 | + | if (storageLayout_ & DataStorage::dslTorque) { |
244 | + | fill(atomRowData.torque.begin(), atomRowData.torque.end(), V3Zero); |
245 | + | fill(atomColData.torque.begin(), atomColData.torque.end(), V3Zero); |
246 | + | } |
247 | + | |
248 | + | fill(pot_row.begin(), pot_row.end(), |
249 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
250 | + | |
251 | + | fill(pot_col.begin(), pot_col.end(), |
252 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
253 | + | |
254 | + | pot_local = Vector<RealType, N_INTERACTION_FAMILIES>(0.0); |
255 | + | |
256 | + | if (storageLayout_ & DataStorage::dslParticlePot) { |
257 | + | fill(atomRowData.particlePot.begin(), atomRowData.particlePot.end(), 0.0); |
258 | + | fill(atomColData.particlePot.begin(), atomColData.particlePot.end(), 0.0); |
259 | + | } |
260 | + | |
261 | + | if (storageLayout_ & DataStorage::dslDensity) { |
262 | + | fill(atomRowData.density.begin(), atomRowData.density.end(), 0.0); |
263 | + | fill(atomColData.density.begin(), atomColData.density.end(), 0.0); |
264 | + | } |
265 | + | |
266 | + | if (storageLayout_ & DataStorage::dslFunctional) { |
267 | + | fill(atomRowData.functional.begin(), atomRowData.functional.end(), 0.0); |
268 | + | fill(atomColData.functional.begin(), atomColData.functional.end(), 0.0); |
269 | + | } |
270 | + | |
271 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { |
272 | + | fill(atomRowData.functionalDerivative.begin(), |
273 | + | atomRowData.functionalDerivative.end(), 0.0); |
274 | + | fill(atomColData.functionalDerivative.begin(), |
275 | + | atomColData.functionalDerivative.end(), 0.0); |
276 | + | } |
277 | + | |
278 | + | #else |
279 | + | |
280 | + | if (storageLayout_ & DataStorage::dslParticlePot) { |
281 | + | fill(snap_->atomData.particlePot.begin(), |
282 | + | snap_->atomData.particlePot.end(), 0.0); |
283 | + | } |
284 | + | |
285 | + | if (storageLayout_ & DataStorage::dslDensity) { |
286 | + | fill(snap_->atomData.density.begin(), |
287 | + | snap_->atomData.density.end(), 0.0); |
288 | + | } |
289 | + | if (storageLayout_ & DataStorage::dslFunctional) { |
290 | + | fill(snap_->atomData.functional.begin(), |
291 | + | snap_->atomData.functional.end(), 0.0); |
292 | + | } |
293 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { |
294 | + | fill(snap_->atomData.functionalDerivative.begin(), |
295 | + | snap_->atomData.functionalDerivative.end(), 0.0); |
296 | + | } |
297 | + | #endif |
298 | + | |
299 | + | } |
300 | + | |
301 | + | |
302 | + | void ForceMatrixDecomposition::distributeData() { |
303 | + | snap_ = sman_->getCurrentSnapshot(); |
304 | + | storageLayout_ = sman_->getStorageLayout(); |
305 | + | #ifdef IS_MPI |
306 | + | |
307 | + | // gather up the atomic positions |
308 | + | AtomCommVectorRow->gather(snap_->atomData.position, |
309 | + | atomRowData.position); |
310 | + | AtomCommVectorColumn->gather(snap_->atomData.position, |
311 | + | atomColData.position); |
312 | + | |
313 | + | // gather up the cutoff group positions |
314 | + | cgCommVectorRow->gather(snap_->cgData.position, |
315 | + | cgRowData.position); |
316 | + | cgCommVectorColumn->gather(snap_->cgData.position, |
317 | + | cgColData.position); |
318 | + | |
319 | + | // if needed, gather the atomic rotation matrices |
320 | + | if (storageLayout_ & DataStorage::dslAmat) { |
321 | + | AtomCommMatrixRow->gather(snap_->atomData.aMat, |
322 | + | atomRowData.aMat); |
323 | + | AtomCommMatrixColumn->gather(snap_->atomData.aMat, |
324 | + | atomColData.aMat); |
325 | + | } |
326 | + | |
327 | + | // if needed, gather the atomic eletrostatic frames |
328 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
329 | + | AtomCommMatrixRow->gather(snap_->atomData.electroFrame, |
330 | + | atomRowData.electroFrame); |
331 | + | AtomCommMatrixColumn->gather(snap_->atomData.electroFrame, |
332 | + | atomColData.electroFrame); |
333 | + | } |
334 | + | #endif |
335 | + | } |
336 | + | |
337 | + | /* collects information obtained during the pre-pair loop onto local |
338 | + | * data structures. |
339 | + | */ |
340 | + | void ForceMatrixDecomposition::collectIntermediateData() { |
341 | + | snap_ = sman_->getCurrentSnapshot(); |
342 | + | storageLayout_ = sman_->getStorageLayout(); |
343 | + | #ifdef IS_MPI |
344 | + | |
345 | + | if (storageLayout_ & DataStorage::dslDensity) { |
346 | + | |
347 | + | AtomCommRealRow->scatter(atomRowData.density, |
348 | + | snap_->atomData.density); |
349 | + | |
350 | + | int n = snap_->atomData.density.size(); |
351 | + | vector<RealType> rho_tmp(n, 0.0); |
352 | + | AtomCommRealColumn->scatter(atomColData.density, rho_tmp); |
353 | + | for (int i = 0; i < n; i++) |
354 | + | snap_->atomData.density[i] += rho_tmp[i]; |
355 | + | } |
356 | + | #endif |
357 | + | } |
358 | + | |
359 | + | /* |
360 | + | * redistributes information obtained during the pre-pair loop out to |
361 | + | * row and column-indexed data structures |
362 | + | */ |
363 | + | void ForceMatrixDecomposition::distributeIntermediateData() { |
364 | + | snap_ = sman_->getCurrentSnapshot(); |
365 | + | storageLayout_ = sman_->getStorageLayout(); |
366 | + | #ifdef IS_MPI |
367 | + | if (storageLayout_ & DataStorage::dslFunctional) { |
368 | + | AtomCommRealRow->gather(snap_->atomData.functional, |
369 | + | atomRowData.functional); |
370 | + | AtomCommRealColumn->gather(snap_->atomData.functional, |
371 | + | atomColData.functional); |
372 | + | } |
373 | + | |
374 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { |
375 | + | AtomCommRealRow->gather(snap_->atomData.functionalDerivative, |
376 | + | atomRowData.functionalDerivative); |
377 | + | AtomCommRealColumn->gather(snap_->atomData.functionalDerivative, |
378 | + | atomColData.functionalDerivative); |
379 | + | } |
380 | + | #endif |
381 | + | } |
382 | + | |
383 | + | |
384 | + | void ForceMatrixDecomposition::collectData() { |
385 | + | snap_ = sman_->getCurrentSnapshot(); |
386 | + | storageLayout_ = sman_->getStorageLayout(); |
387 | + | #ifdef IS_MPI |
388 | + | int n = snap_->atomData.force.size(); |
389 | + | vector<Vector3d> frc_tmp(n, V3Zero); |
390 | + | |
391 | + | AtomCommVectorRow->scatter(atomRowData.force, frc_tmp); |
392 | + | for (int i = 0; i < n; i++) { |
393 | + | snap_->atomData.force[i] += frc_tmp[i]; |
394 | + | frc_tmp[i] = 0.0; |
395 | + | } |
396 | + | |
397 | + | AtomCommVectorColumn->scatter(atomColData.force, frc_tmp); |
398 | + | for (int i = 0; i < n; i++) |
399 | + | snap_->atomData.force[i] += frc_tmp[i]; |
400 | + | |
401 | + | |
402 | + | if (storageLayout_ & DataStorage::dslTorque) { |
403 | + | |
404 | + | int nt = snap_->atomData.force.size(); |
405 | + | vector<Vector3d> trq_tmp(nt, V3Zero); |
406 | + | |
407 | + | AtomCommVectorRow->scatter(atomRowData.torque, trq_tmp); |
408 | + | for (int i = 0; i < n; i++) { |
409 | + | snap_->atomData.torque[i] += trq_tmp[i]; |
410 | + | trq_tmp[i] = 0.0; |
411 | + | } |
412 | + | |
413 | + | AtomCommVectorColumn->scatter(atomColData.torque, trq_tmp); |
414 | + | for (int i = 0; i < n; i++) |
415 | + | snap_->atomData.torque[i] += trq_tmp[i]; |
416 | + | } |
417 | + | |
418 | + | nLocal_ = snap_->getNumberOfAtoms(); |
419 | + | |
420 | + | vector<potVec> pot_temp(nLocal_, |
421 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
422 | + | |
423 | + | // scatter/gather pot_row into the members of my column |
424 | + | |
425 | + | AtomCommPotRow->scatter(pot_row, pot_temp); |
426 | + | |
427 | + | for (int ii = 0; ii < pot_temp.size(); ii++ ) |
428 | + | pot_local += pot_temp[ii]; |
429 | + | |
430 | + | fill(pot_temp.begin(), pot_temp.end(), |
431 | + | Vector<RealType, N_INTERACTION_FAMILIES> (0.0)); |
432 | + | |
433 | + | AtomCommPotColumn->scatter(pot_col, pot_temp); |
434 | + | |
435 | + | for (int ii = 0; ii < pot_temp.size(); ii++ ) |
436 | + | pot_local += pot_temp[ii]; |
437 | + | |
438 | + | #endif |
439 | + | } |
440 | + | |
441 | + | int ForceMatrixDecomposition::getNAtomsInRow() { |
442 | + | #ifdef IS_MPI |
443 | + | return nAtomsInRow_; |
444 | + | #else |
445 | + | return nLocal_; |
446 | + | #endif |
447 | + | } |
448 | + | |
449 | + | /** |
450 | + | * returns the list of atoms belonging to this group. |
451 | + | */ |
452 | + | vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){ |
453 | + | #ifdef IS_MPI |
454 | + | return groupListRow_[cg1]; |
455 | + | #else |
456 | + | return groupList_[cg1]; |
457 | + | #endif |
458 | + | } |
459 | + | |
460 | + | vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){ |
461 | + | #ifdef IS_MPI |
462 | + | return groupListCol_[cg2]; |
463 | + | #else |
464 | + | return groupList_[cg2]; |
465 | + | #endif |
466 | + | } |
467 | + | |
468 | + | Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){ |
469 | + | Vector3d d; |
470 | + | |
471 | + | #ifdef IS_MPI |
472 | + | d = cgColData.position[cg2] - cgRowData.position[cg1]; |
473 | + | #else |
474 | + | d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1]; |
475 | + | #endif |
476 | + | |
477 | + | snap_->wrapVector(d); |
478 | + | return d; |
479 | + | } |
480 | + | |
481 | + | |
482 | + | Vector3d ForceMatrixDecomposition::getAtomToGroupVectorRow(int atom1, int cg1){ |
483 | + | |
484 | + | Vector3d d; |
485 | + | |
486 | + | #ifdef IS_MPI |
487 | + | d = cgRowData.position[cg1] - atomRowData.position[atom1]; |
488 | + | #else |
489 | + | d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1]; |
490 | + | #endif |
491 | + | |
492 | + | snap_->wrapVector(d); |
493 | + | return d; |
494 | + | } |
495 | + | |
496 | + | Vector3d ForceMatrixDecomposition::getAtomToGroupVectorColumn(int atom2, int cg2){ |
497 | + | Vector3d d; |
498 | + | |
499 | + | #ifdef IS_MPI |
500 | + | d = cgColData.position[cg2] - atomColData.position[atom2]; |
501 | + | #else |
502 | + | d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2]; |
503 | + | #endif |
504 | + | |
505 | + | snap_->wrapVector(d); |
506 | + | return d; |
507 | + | } |
508 | + | |
509 | + | RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) { |
510 | + | #ifdef IS_MPI |
511 | + | return massFactorsRow[atom1]; |
512 | + | #else |
513 | + | return massFactorsLocal[atom1]; |
514 | + | #endif |
515 | + | } |
516 | + | |
517 | + | RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) { |
518 | + | #ifdef IS_MPI |
519 | + | return massFactorsCol[atom2]; |
520 | + | #else |
521 | + | return massFactorsLocal[atom2]; |
522 | + | #endif |
523 | + | |
524 | + | } |
525 | + | |
526 | + | Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){ |
527 | + | Vector3d d; |
528 | + | |
529 | + | #ifdef IS_MPI |
530 | + | d = atomColData.position[atom2] - atomRowData.position[atom1]; |
531 | + | #else |
532 | + | d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1]; |
533 | + | #endif |
534 | + | |
535 | + | snap_->wrapVector(d); |
536 | + | return d; |
537 | + | } |
538 | + | |
539 | + | vector<int> ForceMatrixDecomposition::getSkipsForRowAtom(int atom1) { |
540 | + | #ifdef IS_MPI |
541 | + | return skipsForRowAtom[atom1]; |
542 | + | #else |
543 | + | return skipsForLocalAtom[atom1]; |
544 | + | #endif |
545 | + | } |
546 | + | |
547 | + | /** |
548 | + | * There are a number of reasons to skip a pair or a |
549 | + | * particle. Mostly we do this to exclude atoms who are involved in |
550 | + | * short range interactions (bonds, bends, torsions), but we also |
551 | + | * need to exclude some overcounted interactions that result from |
552 | + | * the parallel decomposition. |
553 | + | */ |
554 | + | bool ForceMatrixDecomposition::skipAtomPair(int atom1, int atom2) { |
555 | + | int unique_id_1, unique_id_2; |
556 | + | |
557 | + | #ifdef IS_MPI |
558 | + | // in MPI, we have to look up the unique IDs for each atom |
559 | + | unique_id_1 = AtomRowToGlobal[atom1]; |
560 | + | unique_id_2 = AtomColToGlobal[atom2]; |
561 | + | |
562 | + | // this situation should only arise in MPI simulations |
563 | + | if (unique_id_1 == unique_id_2) return true; |
564 | + | |
565 | + | // this prevents us from doing the pair on multiple processors |
566 | + | if (unique_id_1 < unique_id_2) { |
567 | + | if ((unique_id_1 + unique_id_2) % 2 == 0) return true; |
568 | + | } else { |
569 | + | if ((unique_id_1 + unique_id_2) % 2 == 1) return true; |
570 | + | } |
571 | + | #else |
572 | + | // in the normal loop, the atom numbers are unique |
573 | + | unique_id_1 = atom1; |
574 | + | unique_id_2 = atom2; |
575 | + | #endif |
576 | + | |
577 | + | #ifdef IS_MPI |
578 | + | for (vector<int>::iterator i = skipsForRowAtom[atom1].begin(); |
579 | + | i != skipsForRowAtom[atom1].end(); ++i) { |
580 | + | if ( (*i) == unique_id_2 ) return true; |
581 | + | } |
582 | + | #else |
583 | + | for (vector<int>::iterator i = skipsForLocalAtom[atom1].begin(); |
584 | + | i != skipsForLocalAtom[atom1].end(); ++i) { |
585 | + | if ( (*i) == unique_id_2 ) return true; |
586 | + | } |
587 | + | #endif |
588 | + | } |
589 | + | |
590 | + | int ForceMatrixDecomposition::getTopoDistance(int atom1, int atom2) { |
591 | + | |
592 | + | #ifdef IS_MPI |
593 | + | for (int i = 0; i < toposForRowAtom[atom1].size(); i++) { |
594 | + | if ( toposForRowAtom[atom1][i] == atom2 ) return topoDistRow[atom1][i]; |
595 | + | } |
596 | + | #else |
597 | + | for (int i = 0; i < toposForLocalAtom[atom1].size(); i++) { |
598 | + | if ( toposForLocalAtom[atom1][i] == atom2 ) return topoDistLocal[atom1][i]; |
599 | + | } |
600 | + | #endif |
601 | + | |
602 | + | // zero is default for unconnected (i.e. normal) pair interactions |
603 | + | return 0; |
604 | + | } |
605 | + | |
606 | + | void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){ |
607 | + | #ifdef IS_MPI |
608 | + | atomRowData.force[atom1] += fg; |
609 | + | #else |
610 | + | snap_->atomData.force[atom1] += fg; |
611 | + | #endif |
612 | + | } |
613 | + | |
614 | + | void ForceMatrixDecomposition::addForceToAtomColumn(int atom2, Vector3d fg){ |
615 | + | #ifdef IS_MPI |
616 | + | atomColData.force[atom2] += fg; |
617 | + | #else |
618 | + | snap_->atomData.force[atom2] += fg; |
619 | + | #endif |
620 | + | } |
621 | + | |
622 | + | // filling interaction blocks with pointers |
623 | + | InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) { |
624 | + | InteractionData idat; |
625 | + | |
626 | + | #ifdef IS_MPI |
627 | + | |
628 | + | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
629 | + | ff_->getAtomType(identsCol[atom2]) ); |
630 | + | |
631 | + | |
632 | + | if (storageLayout_ & DataStorage::dslAmat) { |
633 | + | idat.A1 = &(atomRowData.aMat[atom1]); |
634 | + | idat.A2 = &(atomColData.aMat[atom2]); |
635 | + | } |
636 | + | |
637 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
638 | + | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
639 | + | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
640 | + | } |
641 | + | |
642 | + | if (storageLayout_ & DataStorage::dslTorque) { |
643 | + | idat.t1 = &(atomRowData.torque[atom1]); |
644 | + | idat.t2 = &(atomColData.torque[atom2]); |
645 | + | } |
646 | + | |
647 | + | if (storageLayout_ & DataStorage::dslDensity) { |
648 | + | idat.rho1 = &(atomRowData.density[atom1]); |
649 | + | idat.rho2 = &(atomColData.density[atom2]); |
650 | + | } |
651 | + | |
652 | + | if (storageLayout_ & DataStorage::dslFunctional) { |
653 | + | idat.frho1 = &(atomRowData.functional[atom1]); |
654 | + | idat.frho2 = &(atomColData.functional[atom2]); |
655 | + | } |
656 | + | |
657 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { |
658 | + | idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]); |
659 | + | idat.dfrho2 = &(atomColData.functionalDerivative[atom2]); |
660 | + | } |
661 | + | |
662 | + | if (storageLayout_ & DataStorage::dslParticlePot) { |
663 | + | idat.particlePot1 = &(atomRowData.particlePot[atom1]); |
664 | + | idat.particlePot2 = &(atomColData.particlePot[atom2]); |
665 | + | } |
666 | + | |
667 | + | #else |
668 | + | |
669 | + | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
670 | + | ff_->getAtomType(identsLocal[atom2]) ); |
671 | + | |
672 | + | if (storageLayout_ & DataStorage::dslAmat) { |
673 | + | idat.A1 = &(snap_->atomData.aMat[atom1]); |
674 | + | idat.A2 = &(snap_->atomData.aMat[atom2]); |
675 | + | } |
676 | + | |
677 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
678 | + | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
679 | + | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
680 | + | } |
681 | + | |
682 | + | if (storageLayout_ & DataStorage::dslTorque) { |
683 | + | idat.t1 = &(snap_->atomData.torque[atom1]); |
684 | + | idat.t2 = &(snap_->atomData.torque[atom2]); |
685 | + | } |
686 | + | |
687 | + | if (storageLayout_ & DataStorage::dslDensity) { |
688 | + | idat.rho1 = &(snap_->atomData.density[atom1]); |
689 | + | idat.rho2 = &(snap_->atomData.density[atom2]); |
690 | + | } |
691 | + | |
692 | + | if (storageLayout_ & DataStorage::dslFunctional) { |
693 | + | idat.frho1 = &(snap_->atomData.functional[atom1]); |
694 | + | idat.frho2 = &(snap_->atomData.functional[atom2]); |
695 | + | } |
696 | + | |
697 | + | if (storageLayout_ & DataStorage::dslFunctionalDerivative) { |
698 | + | idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]); |
699 | + | idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]); |
700 | + | } |
701 | + | |
702 | + | if (storageLayout_ & DataStorage::dslParticlePot) { |
703 | + | idat.particlePot1 = &(snap_->atomData.particlePot[atom1]); |
704 | + | idat.particlePot2 = &(snap_->atomData.particlePot[atom2]); |
705 | + | } |
706 | + | |
707 | + | #endif |
708 | + | return idat; |
709 | + | } |
710 | + | |
711 | + | |
712 | + | void ForceMatrixDecomposition::unpackInteractionData(InteractionData idat, int atom1, int atom2) { |
713 | + | #ifdef IS_MPI |
714 | + | pot_row[atom1] += 0.5 * *(idat.pot); |
715 | + | pot_col[atom2] += 0.5 * *(idat.pot); |
716 | + | |
717 | + | atomRowData.force[atom1] += *(idat.f1); |
718 | + | atomColData.force[atom2] -= *(idat.f1); |
719 | + | #else |
720 | + | longRangePot_ += *(idat.pot); |
721 | + | |
722 | + | snap_->atomData.force[atom1] += *(idat.f1); |
723 | + | snap_->atomData.force[atom2] -= *(idat.f1); |
724 | + | #endif |
725 | + | |
726 | + | } |
727 | + | |
728 | + | |
729 | + | InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){ |
730 | + | |
731 | + | InteractionData idat; |
732 | + | #ifdef IS_MPI |
733 | + | idat.atypes = make_pair( ff_->getAtomType(identsRow[atom1]), |
734 | + | ff_->getAtomType(identsCol[atom2]) ); |
735 | + | |
736 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
737 | + | idat.eFrame1 = &(atomRowData.electroFrame[atom1]); |
738 | + | idat.eFrame2 = &(atomColData.electroFrame[atom2]); |
739 | + | } |
740 | + | if (storageLayout_ & DataStorage::dslTorque) { |
741 | + | idat.t1 = &(atomRowData.torque[atom1]); |
742 | + | idat.t2 = &(atomColData.torque[atom2]); |
743 | + | } |
744 | + | #else |
745 | + | idat.atypes = make_pair( ff_->getAtomType(identsLocal[atom1]), |
746 | + | ff_->getAtomType(identsLocal[atom2]) ); |
747 | + | |
748 | + | if (storageLayout_ & DataStorage::dslElectroFrame) { |
749 | + | idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]); |
750 | + | idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]); |
751 | + | } |
752 | + | if (storageLayout_ & DataStorage::dslTorque) { |
753 | + | idat.t1 = &(snap_->atomData.torque[atom1]); |
754 | + | idat.t2 = &(snap_->atomData.torque[atom2]); |
755 | + | } |
756 | + | #endif |
757 | + | } |
758 | + | |
759 | + | /* |
760 | + | * buildNeighborList |
761 | + | * |
762 | + | * first element of pair is row-indexed CutoffGroup |
763 | + | * second element of pair is column-indexed CutoffGroup |
764 | + | */ |
765 | + | vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() { |
766 | + | |
767 | + | vector<pair<int, int> > neighborList; |
768 | + | #ifdef IS_MPI |
769 | + | cellListRow_.clear(); |
770 | + | cellListCol_.clear(); |
771 | + | #else |
772 | + | cellList_.clear(); |
773 | + | #endif |
774 | + | |
775 | + | // dangerous to not do error checking. |
776 | + | RealType rCut_; |
777 | + | |
778 | + | RealType rList_ = (rCut_ + skinThickness_); |
779 | + | RealType rl2 = rList_ * rList_; |
780 | + | Snapshot* snap_ = sman_->getCurrentSnapshot(); |
781 | + | Mat3x3d Hmat = snap_->getHmat(); |
782 | + | Vector3d Hx = Hmat.getColumn(0); |
783 | + | Vector3d Hy = Hmat.getColumn(1); |
784 | + | Vector3d Hz = Hmat.getColumn(2); |
785 | + | |
786 | + | nCells_.x() = (int) ( Hx.length() )/ rList_; |
787 | + | nCells_.y() = (int) ( Hy.length() )/ rList_; |
788 | + | nCells_.z() = (int) ( Hz.length() )/ rList_; |
789 | + | |
790 | + | Mat3x3d invHmat = snap_->getInvHmat(); |
791 | + | Vector3d rs, scaled, dr; |
792 | + | Vector3i whichCell; |
793 | + | int cellIndex; |
794 | + | |
795 | + | #ifdef IS_MPI |
796 | + | for (int i = 0; i < nGroupsInRow_; i++) { |
797 | + | rs = cgRowData.position[i]; |
798 | + | // scaled positions relative to the box vectors |
799 | + | scaled = invHmat * rs; |
800 | + | // wrap the vector back into the unit box by subtracting integer box |
801 | + | // numbers |
802 | + | for (int j = 0; j < 3; j++) |
803 | + | scaled[j] -= roundMe(scaled[j]); |
804 | + | |
805 | + | // find xyz-indices of cell that cutoffGroup is in. |
806 | + | whichCell.x() = nCells_.x() * scaled.x(); |
807 | + | whichCell.y() = nCells_.y() * scaled.y(); |
808 | + | whichCell.z() = nCells_.z() * scaled.z(); |
809 | + | |
810 | + | // find single index of this cell: |
811 | + | cellIndex = Vlinear(whichCell, nCells_); |
812 | + | // add this cutoff group to the list of groups in this cell; |
813 | + | cellListRow_[cellIndex].push_back(i); |
814 | + | } |
815 | + | |
816 | + | for (int i = 0; i < nGroupsInCol_; i++) { |
817 | + | rs = cgColData.position[i]; |
818 | + | // scaled positions relative to the box vectors |
819 | + | scaled = invHmat * rs; |
820 | + | // wrap the vector back into the unit box by subtracting integer box |
821 | + | // numbers |
822 | + | for (int j = 0; j < 3; j++) |
823 | + | scaled[j] -= roundMe(scaled[j]); |
824 | + | |
825 | + | // find xyz-indices of cell that cutoffGroup is in. |
826 | + | whichCell.x() = nCells_.x() * scaled.x(); |
827 | + | whichCell.y() = nCells_.y() * scaled.y(); |
828 | + | whichCell.z() = nCells_.z() * scaled.z(); |
829 | + | |
830 | + | // find single index of this cell: |
831 | + | cellIndex = Vlinear(whichCell, nCells_); |
832 | + | // add this cutoff group to the list of groups in this cell; |
833 | + | cellListCol_[cellIndex].push_back(i); |
834 | + | } |
835 | + | #else |
836 | + | for (int i = 0; i < nGroups_; i++) { |
837 | + | rs = snap_->cgData.position[i]; |
838 | + | // scaled positions relative to the box vectors |
839 | + | scaled = invHmat * rs; |
840 | + | // wrap the vector back into the unit box by subtracting integer box |
841 | + | // numbers |
842 | + | for (int j = 0; j < 3; j++) |
843 | + | scaled[j] -= roundMe(scaled[j]); |
844 | + | |
845 | + | // find xyz-indices of cell that cutoffGroup is in. |
846 | + | whichCell.x() = nCells_.x() * scaled.x(); |
847 | + | whichCell.y() = nCells_.y() * scaled.y(); |
848 | + | whichCell.z() = nCells_.z() * scaled.z(); |
849 | + | |
850 | + | // find single index of this cell: |
851 | + | cellIndex = Vlinear(whichCell, nCells_); |
852 | + | // add this cutoff group to the list of groups in this cell; |
853 | + | cellList_[cellIndex].push_back(i); |
854 | + | } |
855 | + | #endif |
856 | + | |
857 | + | for (int m1z = 0; m1z < nCells_.z(); m1z++) { |
858 | + | for (int m1y = 0; m1y < nCells_.y(); m1y++) { |
859 | + | for (int m1x = 0; m1x < nCells_.x(); m1x++) { |
860 | + | Vector3i m1v(m1x, m1y, m1z); |
861 | + | int m1 = Vlinear(m1v, nCells_); |
862 | + | |
863 | + | for (vector<Vector3i>::iterator os = cellOffsets_.begin(); |
864 | + | os != cellOffsets_.end(); ++os) { |
865 | + | |
866 | + | Vector3i m2v = m1v + (*os); |
867 | + | |
868 | + | if (m2v.x() >= nCells_.x()) { |
869 | + | m2v.x() = 0; |
870 | + | } else if (m2v.x() < 0) { |
871 | + | m2v.x() = nCells_.x() - 1; |
872 | + | } |
873 | + | |
874 | + | if (m2v.y() >= nCells_.y()) { |
875 | + | m2v.y() = 0; |
876 | + | } else if (m2v.y() < 0) { |
877 | + | m2v.y() = nCells_.y() - 1; |
878 | + | } |
879 | + | |
880 | + | if (m2v.z() >= nCells_.z()) { |
881 | + | m2v.z() = 0; |
882 | + | } else if (m2v.z() < 0) { |
883 | + | m2v.z() = nCells_.z() - 1; |
884 | + | } |
885 | + | |
886 | + | int m2 = Vlinear (m2v, nCells_); |
887 | + | |
888 | + | #ifdef IS_MPI |
889 | + | for (vector<int>::iterator j1 = cellListRow_[m1].begin(); |
890 | + | j1 != cellListRow_[m1].end(); ++j1) { |
891 | + | for (vector<int>::iterator j2 = cellListCol_[m2].begin(); |
892 | + | j2 != cellListCol_[m2].end(); ++j2) { |
893 | + | |
894 | + | // Always do this if we're in different cells or if |
895 | + | // we're in the same cell and the global index of the |
896 | + | // j2 cutoff group is less than the j1 cutoff group |
897 | + | |
898 | + | if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) { |
899 | + | dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)]; |
900 | + | snap_->wrapVector(dr); |
901 | + | if (dr.lengthSquare() < rl2) { |
902 | + | neighborList.push_back(make_pair((*j1), (*j2))); |
903 | + | } |
904 | + | } |
905 | + | } |
906 | + | } |
907 | + | #else |
908 | + | for (vector<int>::iterator j1 = cellList_[m1].begin(); |
909 | + | j1 != cellList_[m1].end(); ++j1) { |
910 | + | for (vector<int>::iterator j2 = cellList_[m2].begin(); |
911 | + | j2 != cellList_[m2].end(); ++j2) { |
912 | + | |
913 | + | // Always do this if we're in different cells or if |
914 | + | // we're in the same cell and the global index of the |
915 | + | // j2 cutoff group is less than the j1 cutoff group |
916 | + | |
917 | + | if (m2 != m1 || (*j2) < (*j1)) { |
918 | + | dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)]; |
919 | + | snap_->wrapVector(dr); |
920 | + | if (dr.lengthSquare() < rl2) { |
921 | + | neighborList.push_back(make_pair((*j1), (*j2))); |
922 | + | } |
923 | + | } |
924 | + | } |
925 | + | } |
926 | + | #endif |
927 | + | } |
928 | + | } |
929 | + | } |
930 | + | } |
931 | + | |
932 | + | // save the local cutoff group positions for the check that is |
933 | + | // done on each loop: |
934 | + | saved_CG_positions_.clear(); |
935 | + | for (int i = 0; i < nGroups_; i++) |
936 | + | saved_CG_positions_.push_back(snap_->cgData.position[i]); |
937 | + | |
938 | + | return neighborList; |
939 | + | } |
940 | + | } //end namespace OpenMD |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |