ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/parallel/ForceMatrixDecomposition.cpp
(Generate patch)

Comparing:
branches/development/src/parallel/ForceDecomposition.cpp (file contents), Revision 1541 by gezelter, Fri Feb 4 20:04:56 2011 UTC vs.
branches/development/src/parallel/ForceMatrixDecomposition.cpp (file contents), Revision 1569 by gezelter, Thu May 26 13:55:04 2011 UTC

# Line 38 | Line 38
38   * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39   * [4]  Vardeman & Gezelter, in progress (2009).                        
40   */
41 < #include "parallel/ForceDecomposition.hpp"
42 < #include "parallel/Communicator.hpp"
41 > #include "parallel/ForceMatrixDecomposition.hpp"
42   #include "math/SquareMatrix3.hpp"
43 + #include "nonbonded/NonBondedInteraction.hpp"
44 + #include "brains/SnapshotManager.hpp"
45  
46   using namespace std;
47   namespace OpenMD {
48  
49 <  void ForceDecomposition::distributeInitialData() {
49 >  /**
50 >   * distributeInitialData is essentially a copy of the older fortran
51 >   * SimulationSetup
52 >   */
53 >  
54 >  void ForceMatrixDecomposition::distributeInitialData() {
55 >    snap_ = sman_->getCurrentSnapshot();
56 >    storageLayout_ = sman_->getStorageLayout();
57 >    nLocal_ = snap_->getNumberOfAtoms();
58 >    nGroups_ = snap_->getNumberOfCutoffGroups();
59 >
60 >    // gather the information for atomtype IDs (atids):
61 >    vector<int> identsLocal = info_->getIdentArray();
62 >    AtomLocalToGlobal = info_->getGlobalAtomIndices();
63 >    cgLocalToGlobal = info_->getGlobalGroupIndices();
64 >    vector<int> globalGroupMembership = info_->getGlobalGroupMembership();
65 >    vector<RealType> massFactorsLocal = info_->getMassFactors();
66 >    vector<RealType> pot_local(N_INTERACTION_FAMILIES, 0.0);
67 >
68   #ifdef IS_MPI
69 <    Snapshot* snap = sman_->getCurrentSnapshot();
70 <    int nAtoms = snap->getNumberOfAtoms();
71 <    int nGroups = snap->getNumberOfCutoffGroups();
69 >
70 >    AtomCommIntRow = new Communicator<Row,int>(nLocal_);
71 >    AtomCommRealRow = new Communicator<Row,RealType>(nLocal_);
72 >    AtomCommVectorRow = new Communicator<Row,Vector3d>(nLocal_);
73 >    AtomCommMatrixRow = new Communicator<Row,Mat3x3d>(nLocal_);
74  
75 <    AtomCommRealI = new Communicator<Row,RealType>(nAtoms);
76 <    AtomCommVectorI = new Communicator<Row,Vector3d>(nAtoms);
77 <    AtomCommMatrixI = new Communicator<Row,Mat3x3d>(nAtoms);
75 >    AtomCommIntColumn = new Communicator<Column,int>(nLocal_);
76 >    AtomCommRealColumn = new Communicator<Column,RealType>(nLocal_);
77 >    AtomCommVectorColumn = new Communicator<Column,Vector3d>(nLocal_);
78 >    AtomCommMatrixColumn = new Communicator<Column,Mat3x3d>(nLocal_);
79  
80 <    AtomCommRealJ = new Communicator<Column,RealType>(nAtoms);
81 <    AtomCommVectorJ = new Communicator<Column,Vector3d>(nAtoms);
82 <    AtomCommMatrixJ = new Communicator<Column,Mat3x3d>(nAtoms);
80 >    cgCommIntRow = new Communicator<Row,int>(nGroups_);
81 >    cgCommVectorRow = new Communicator<Row,Vector3d>(nGroups_);
82 >    cgCommIntColumn = new Communicator<Column,int>(nGroups_);
83 >    cgCommVectorColumn = new Communicator<Column,Vector3d>(nGroups_);
84  
85 <    cgCommVectorI = new Communicator<Row,Vector3d>(nGroups);
86 <    cgCommVectorJ = new Communicator<Column,Vector3d>(nGroups);
85 >    nAtomsInRow_ = AtomCommIntRow->getSize();
86 >    nAtomsInCol_ = AtomCommIntColumn->getSize();
87 >    nGroupsInRow_ = cgCommIntRow->getSize();
88 >    nGroupsInCol_ = cgCommIntColumn->getSize();
89  
90 <    int nInRow = AtomCommRealI.getSize();
91 <    int nInCol = AtomCommRealJ.getSize();
90 >    // Modify the data storage objects with the correct layouts and sizes:
91 >    atomRowData.resize(nAtomsInRow_);
92 >    atomRowData.setStorageLayout(storageLayout_);
93 >    atomColData.resize(nAtomsInCol_);
94 >    atomColData.setStorageLayout(storageLayout_);
95 >    cgRowData.resize(nGroupsInRow_);
96 >    cgRowData.setStorageLayout(DataStorage::dslPosition);
97 >    cgColData.resize(nGroupsInCol_);
98 >    cgColData.setStorageLayout(DataStorage::dslPosition);
99 >    
100 >    vector<vector<RealType> > pot_row(N_INTERACTION_FAMILIES,
101 >                                      vector<RealType> (nAtomsInRow_, 0.0));
102 >    vector<vector<RealType> > pot_col(N_INTERACTION_FAMILIES,
103 >                                      vector<RealType> (nAtomsInCol_, 0.0));
104 >    
105 >    identsRow.reserve(nAtomsInRow_);
106 >    identsCol.reserve(nAtomsInCol_);
107 >    
108 >    AtomCommIntRow->gather(identsLocal, identsRow);
109 >    AtomCommIntColumn->gather(identsLocal, identsCol);
110 >    
111 >    AtomCommIntRow->gather(AtomLocalToGlobal, AtomRowToGlobal);
112 >    AtomCommIntColumn->gather(AtomLocalToGlobal, AtomColToGlobal);
113 >    
114 >    cgCommIntRow->gather(cgLocalToGlobal, cgRowToGlobal);
115 >    cgCommIntColumn->gather(cgLocalToGlobal, cgColToGlobal);
116  
117 <    vector<vector<RealType> > pot_row(LR_POT_TYPES,
118 <                                      vector<RealType> (nInRow, 0.0));
70 <    vector<vector<RealType> > pot_col(LR_POT_TYPES,
71 <                                      vector<RealType> (nInCol, 0.0));
117 >    AtomCommRealRow->gather(massFactorsLocal, massFactorsRow);
118 >    AtomCommRealColumn->gather(massFactorsLocal, massFactorsCol);
119  
120 <    vector<vector<RealType> > pot_local(LR_POT_TYPES,
121 <                                        vector<RealType> (nAtoms, 0.0));
120 >    groupListRow_.clear();
121 >    groupListRow_.reserve(nGroupsInRow_);
122 >    for (int i = 0; i < nGroupsInRow_; i++) {
123 >      int gid = cgRowToGlobal[i];
124 >      for (int j = 0; j < nAtomsInRow_; j++) {
125 >        int aid = AtomRowToGlobal[j];
126 >        if (globalGroupMembership[aid] == gid)
127 >          groupListRow_[i].push_back(j);
128 >      }      
129 >    }
130  
131 +    groupListCol_.clear();
132 +    groupListCol_.reserve(nGroupsInCol_);
133 +    for (int i = 0; i < nGroupsInCol_; i++) {
134 +      int gid = cgColToGlobal[i];
135 +      for (int j = 0; j < nAtomsInCol_; j++) {
136 +        int aid = AtomColToGlobal[j];
137 +        if (globalGroupMembership[aid] == gid)
138 +          groupListCol_[i].push_back(j);
139 +      }      
140 +    }
141 +
142   #endif
143 +
144 +    groupList_.clear();
145 +    groupList_.reserve(nGroups_);
146 +    for (int i = 0; i < nGroups_; i++) {
147 +      int gid = cgLocalToGlobal[i];
148 +      for (int j = 0; j < nLocal_; j++) {
149 +        int aid = AtomLocalToGlobal[j];
150 +        if (globalGroupMembership[aid] == gid)
151 +          groupList_[i].push_back(j);
152 +      }      
153 +    }
154 +
155 +  
156 +    // still need:
157 +    // topoDist
158 +    // exclude
159 +
160    }
161      
162  
163  
164 <  void ForceDecomposition::distributeData()  {
164 >  void ForceMatrixDecomposition::distributeData()  {
165 >    snap_ = sman_->getCurrentSnapshot();
166 >    storageLayout_ = sman_->getStorageLayout();
167   #ifdef IS_MPI
83 –    Snapshot* snap = sman_->getCurrentSnapshot();
168      
169      // gather up the atomic positions
170 <    AtomCommVectorI->gather(snap->atomData.position,
171 <                            snap->atomIData.position);
172 <    AtomCommVectorJ->gather(snap->atomData.position,
173 <                            snap->atomJData.position);
170 >    AtomCommVectorRow->gather(snap_->atomData.position,
171 >                              atomRowData.position);
172 >    AtomCommVectorColumn->gather(snap_->atomData.position,
173 >                                 atomColData.position);
174      
175      // gather up the cutoff group positions
176 <    cgCommVectorI->gather(snap->cgData.position,
177 <                          snap->cgIData.position);
178 <    cgCommVectorJ->gather(snap->cgData.position,
179 <                          snap->cgJData.position);
176 >    cgCommVectorRow->gather(snap_->cgData.position,
177 >                            cgRowData.position);
178 >    cgCommVectorColumn->gather(snap_->cgData.position,
179 >                               cgColData.position);
180      
181      // if needed, gather the atomic rotation matrices
182 <    if (snap->atomData.getStorageLayout() & DataStorage::dslAmat) {
183 <      AtomCommMatrixI->gather(snap->atomData.aMat,
184 <                              snap->atomIData.aMat);
185 <      AtomCommMatrixJ->gather(snap->atomData.aMat,
186 <                              snap->atomJData.aMat);
182 >    if (storageLayout_ & DataStorage::dslAmat) {
183 >      AtomCommMatrixRow->gather(snap_->atomData.aMat,
184 >                                atomRowData.aMat);
185 >      AtomCommMatrixColumn->gather(snap_->atomData.aMat,
186 >                                   atomColData.aMat);
187      }
188      
189      // if needed, gather the atomic eletrostatic frames
190 <    if (snap->atomData.getStorageLayout() & DataStorage::dslElectroFrame) {
191 <      AtomCommMatrixI->gather(snap->atomData.electroFrame,
192 <                              snap->atomIData.electroFrame);
193 <      AtomCommMatrixJ->gather(snap->atomData.electroFrame,
194 <                              snap->atomJData.electroFrame);
190 >    if (storageLayout_ & DataStorage::dslElectroFrame) {
191 >      AtomCommMatrixRow->gather(snap_->atomData.electroFrame,
192 >                                atomRowData.electroFrame);
193 >      AtomCommMatrixColumn->gather(snap_->atomData.electroFrame,
194 >                                   atomColData.electroFrame);
195      }
196   #endif      
197    }
198    
199 <  void ForceDecomposition::collectIntermediateData() {
200 < #ifdef IS_MPI
201 <    Snapshot* snap = sman_->getCurrentSnapshot();
199 >  void ForceMatrixDecomposition::collectIntermediateData() {
200 >    snap_ = sman_->getCurrentSnapshot();
201 >    storageLayout_ = sman_->getStorageLayout();
202 > #ifdef IS_MPI
203      
204 <    if (snap->atomData.getStorageLayout() & DataStorage::dslDensity) {
205 <
206 <      AtomCommRealI->scatter(snap->atomIData.density,
207 <                             snap->atomData.density);
208 <
209 <      int n = snap->atomData.density.size();
204 >    if (storageLayout_ & DataStorage::dslDensity) {
205 >      
206 >      AtomCommRealRow->scatter(atomRowData.density,
207 >                               snap_->atomData.density);
208 >      
209 >      int n = snap_->atomData.density.size();
210        std::vector<RealType> rho_tmp(n, 0.0);
211 <      AtomCommRealJ->scatter(snap->atomJData.density, rho_tmp);
211 >      AtomCommRealColumn->scatter(atomColData.density, rho_tmp);
212        for (int i = 0; i < n; i++)
213 <        snap->atomData.density[i] += rho_tmp[i];
213 >        snap_->atomData.density[i] += rho_tmp[i];
214      }
215   #endif
216    }
217    
218 <  void ForceDecomposition::distributeIntermediateData() {
218 >  void ForceMatrixDecomposition::distributeIntermediateData() {
219 >    snap_ = sman_->getCurrentSnapshot();
220 >    storageLayout_ = sman_->getStorageLayout();
221   #ifdef IS_MPI
222 <    Snapshot* snap = sman_->getCurrentSnapshot();
223 <    if (snap->atomData.getStorageLayout() & DataStorage::dslFunctional) {
224 <      AtomCommRealI->gather(snap->atomData.functional,
225 <                            snap->atomIData.functional);
226 <      AtomCommRealJ->gather(snap->atomData.functional,
140 <                            snap->atomJData.functional);
222 >    if (storageLayout_ & DataStorage::dslFunctional) {
223 >      AtomCommRealRow->gather(snap_->atomData.functional,
224 >                              atomRowData.functional);
225 >      AtomCommRealColumn->gather(snap_->atomData.functional,
226 >                                 atomColData.functional);
227      }
228      
229 <    if (snap->atomData.getStorageLayout() & DataStorage::dslFunctionalDerivative) {
230 <      AtomCommRealI->gather(snap->atomData.functionalDerivative,
231 <                            snap->atomIData.functionalDerivative);
232 <      AtomCommRealJ->gather(snap->atomData.functionalDerivative,
233 <                            snap->atomJData.functionalDerivative);
229 >    if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
230 >      AtomCommRealRow->gather(snap_->atomData.functionalDerivative,
231 >                              atomRowData.functionalDerivative);
232 >      AtomCommRealColumn->gather(snap_->atomData.functionalDerivative,
233 >                                 atomColData.functionalDerivative);
234      }
235   #endif
236    }
237    
238    
239 <  void ForceDecomposition::collectData() {
240 < #ifdef IS_MPI
241 <    Snapshot* snap = sman_->getCurrentSnapshot();
239 >  void ForceMatrixDecomposition::collectData() {
240 >    snap_ = sman_->getCurrentSnapshot();
241 >    storageLayout_ = sman_->getStorageLayout();
242 > #ifdef IS_MPI    
243 >    int n = snap_->atomData.force.size();
244 >    vector<Vector3d> frc_tmp(n, V3Zero);
245      
246 <    int n = snap->atomData.force.size();
158 <    std::vector<Vector3d> frc_tmp(n, 0.0);
159 <    
160 <    AtomCommVectorI->scatter(snap->atomIData.force, frc_tmp);
246 >    AtomCommVectorRow->scatter(atomRowData.force, frc_tmp);
247      for (int i = 0; i < n; i++) {
248 <      snap->atomData.force[i] += frc_tmp[i];
248 >      snap_->atomData.force[i] += frc_tmp[i];
249        frc_tmp[i] = 0.0;
250      }
251      
252 <    AtomCommVectorJ->scatter(snap->atomJData.force, frc_tmp);
252 >    AtomCommVectorColumn->scatter(atomColData.force, frc_tmp);
253      for (int i = 0; i < n; i++)
254 <      snap->atomData.force[i] += frc_tmp[i];
254 >      snap_->atomData.force[i] += frc_tmp[i];
255      
256      
257 <    if (snap->atomData.getStorageLayout() & DataStorage::dslTorque) {
257 >    if (storageLayout_ & DataStorage::dslTorque) {
258  
259 <      int nt = snap->atomData.force.size();
260 <      std::vector<Vector3d> trq_tmp(nt, 0.0);
259 >      int nt = snap_->atomData.force.size();
260 >      vector<Vector3d> trq_tmp(nt, V3Zero);
261  
262 <      AtomCommVectorI->scatter(snap->atomIData.torque, trq_tmp);
262 >      AtomCommVectorRow->scatter(atomRowData.torque, trq_tmp);
263        for (int i = 0; i < n; i++) {
264 <        snap->atomData.torque[i] += trq_tmp[i];
264 >        snap_->atomData.torque[i] += trq_tmp[i];
265          trq_tmp[i] = 0.0;
266        }
267        
268 <      AtomCommVectorJ->scatter(snap->atomJData.torque, trq_tmp);
268 >      AtomCommVectorColumn->scatter(atomColData.torque, trq_tmp);
269        for (int i = 0; i < n; i++)
270 <        snap->atomData.torque[i] += trq_tmp[i];
270 >        snap_->atomData.torque[i] += trq_tmp[i];
271      }
272      
273 +    nLocal_ = snap_->getNumberOfAtoms();
274 +
275 +    vector<vector<RealType> > pot_temp(N_INTERACTION_FAMILIES,
276 +                                       vector<RealType> (nLocal_, 0.0));
277      
278 <    vector<vector<RealType> > pot_temp(LR_POT_TYPES,
279 <                                       vector<RealType> (nAtoms, 0.0));
190 <    
191 <    for (int i = 0; i < LR_POT_TYPES; i++) {
192 <      AtomCommRealI->scatter(pot_row[i], pot_temp[i]);
278 >    for (int i = 0; i < N_INTERACTION_FAMILIES; i++) {
279 >      AtomCommRealRow->scatter(pot_row[i], pot_temp[i]);
280        for (int ii = 0;  ii < pot_temp[i].size(); ii++ ) {
281          pot_local[i] += pot_temp[i][ii];
282        }
283      }
284 + #endif
285 +  }
286 +
287 +  /**
288 +   * returns the list of atoms belonging to this group.  
289 +   */
290 +  vector<int> ForceMatrixDecomposition::getAtomsInGroupRow(int cg1){
291 + #ifdef IS_MPI
292 +    return groupListRow_[cg1];
293 + #else
294 +    return groupList_[cg1];
295 + #endif
296 +  }
297 +
298 +  vector<int> ForceMatrixDecomposition::getAtomsInGroupColumn(int cg2){
299 + #ifdef IS_MPI
300 +    return groupListCol_[cg2];
301 + #else
302 +    return groupList_[cg2];
303 + #endif
304 +  }
305 +  
306 +  Vector3d ForceMatrixDecomposition::getIntergroupVector(int cg1, int cg2){
307 +    Vector3d d;
308      
309 + #ifdef IS_MPI
310 +    d = cgColData.position[cg2] - cgRowData.position[cg1];
311 + #else
312 +    d = snap_->cgData.position[cg2] - snap_->cgData.position[cg1];
313 + #endif
314 +    
315 +    snap_->wrapVector(d);
316 +    return d;    
317 +  }
318  
319  
320 +  Vector3d ForceMatrixDecomposition::getAtomToGroupVectorRow(int atom1, int cg1){
321 +
322 +    Vector3d d;
323 +    
324 + #ifdef IS_MPI
325 +    d = cgRowData.position[cg1] - atomRowData.position[atom1];
326 + #else
327 +    d = snap_->cgData.position[cg1] - snap_->atomData.position[atom1];
328   #endif
329 +
330 +    snap_->wrapVector(d);
331 +    return d;    
332    }
333    
334 +  Vector3d ForceMatrixDecomposition::getAtomToGroupVectorColumn(int atom2, int cg2){
335 +    Vector3d d;
336 +    
337 + #ifdef IS_MPI
338 +    d = cgColData.position[cg2] - atomColData.position[atom2];
339 + #else
340 +    d = snap_->cgData.position[cg2] - snap_->atomData.position[atom2];
341 + #endif
342 +    
343 +    snap_->wrapVector(d);
344 +    return d;    
345 +  }
346 +
347 +  RealType ForceMatrixDecomposition::getMassFactorRow(int atom1) {
348 + #ifdef IS_MPI
349 +    return massFactorsRow[atom1];
350 + #else
351 +    return massFactorsLocal[atom1];
352 + #endif
353 +  }
354 +
355 +  RealType ForceMatrixDecomposition::getMassFactorColumn(int atom2) {
356 + #ifdef IS_MPI
357 +    return massFactorsCol[atom2];
358 + #else
359 +    return massFactorsLocal[atom2];
360 + #endif
361 +
362 +  }
363 +    
364 +  Vector3d ForceMatrixDecomposition::getInteratomicVector(int atom1, int atom2){
365 +    Vector3d d;
366 +    
367 + #ifdef IS_MPI
368 +    d = atomColData.position[atom2] - atomRowData.position[atom1];
369 + #else
370 +    d = snap_->atomData.position[atom2] - snap_->atomData.position[atom1];
371 + #endif
372 +
373 +    snap_->wrapVector(d);
374 +    return d;    
375 +  }
376 +
377 +  void ForceMatrixDecomposition::addForceToAtomRow(int atom1, Vector3d fg){
378 + #ifdef IS_MPI
379 +    atomRowData.force[atom1] += fg;
380 + #else
381 +    snap_->atomData.force[atom1] += fg;
382 + #endif
383 +  }
384 +
385 +  void ForceMatrixDecomposition::addForceToAtomColumn(int atom2, Vector3d fg){
386 + #ifdef IS_MPI
387 +    atomColData.force[atom2] += fg;
388 + #else
389 +    snap_->atomData.force[atom2] += fg;
390 + #endif
391 +  }
392 +
393 +    // filling interaction blocks with pointers
394 +  InteractionData ForceMatrixDecomposition::fillInteractionData(int atom1, int atom2) {    
395 +    InteractionData idat;
396 +
397 + #ifdef IS_MPI
398 +    if (storageLayout_ & DataStorage::dslAmat) {
399 +      idat.A1 = &(atomRowData.aMat[atom1]);
400 +      idat.A2 = &(atomColData.aMat[atom2]);
401 +    }
402 +    
403 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
404 +      idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
405 +      idat.eFrame2 = &(atomColData.electroFrame[atom2]);
406 +    }
407 +
408 +    if (storageLayout_ & DataStorage::dslTorque) {
409 +      idat.t1 = &(atomRowData.torque[atom1]);
410 +      idat.t2 = &(atomColData.torque[atom2]);
411 +    }
412 +
413 +    if (storageLayout_ & DataStorage::dslDensity) {
414 +      idat.rho1 = &(atomRowData.density[atom1]);
415 +      idat.rho2 = &(atomColData.density[atom2]);
416 +    }
417 +
418 +    if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
419 +      idat.dfrho1 = &(atomRowData.functionalDerivative[atom1]);
420 +      idat.dfrho2 = &(atomColData.functionalDerivative[atom2]);
421 +    }
422 + #else
423 +    if (storageLayout_ & DataStorage::dslAmat) {
424 +      idat.A1 = &(snap_->atomData.aMat[atom1]);
425 +      idat.A2 = &(snap_->atomData.aMat[atom2]);
426 +    }
427 +
428 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
429 +      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
430 +      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
431 +    }
432 +
433 +    if (storageLayout_ & DataStorage::dslTorque) {
434 +      idat.t1 = &(snap_->atomData.torque[atom1]);
435 +      idat.t2 = &(snap_->atomData.torque[atom2]);
436 +    }
437 +
438 +    if (storageLayout_ & DataStorage::dslDensity) {
439 +      idat.rho1 = &(snap_->atomData.density[atom1]);
440 +      idat.rho2 = &(snap_->atomData.density[atom2]);
441 +    }
442 +
443 +    if (storageLayout_ & DataStorage::dslFunctionalDerivative) {
444 +      idat.dfrho1 = &(snap_->atomData.functionalDerivative[atom1]);
445 +      idat.dfrho2 = &(snap_->atomData.functionalDerivative[atom2]);
446 +    }
447 + #endif
448 +    return idat;
449 +  }
450 +
451 +  InteractionData ForceMatrixDecomposition::fillSkipData(int atom1, int atom2){
452 +
453 +    InteractionData idat;
454 + #ifdef IS_MPI
455 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
456 +      idat.eFrame1 = &(atomRowData.electroFrame[atom1]);
457 +      idat.eFrame2 = &(atomColData.electroFrame[atom2]);
458 +    }
459 +    if (storageLayout_ & DataStorage::dslTorque) {
460 +      idat.t1 = &(atomRowData.torque[atom1]);
461 +      idat.t2 = &(atomColData.torque[atom2]);
462 +    }
463 +    if (storageLayout_ & DataStorage::dslForce) {
464 +      idat.t1 = &(atomRowData.force[atom1]);
465 +      idat.t2 = &(atomColData.force[atom2]);
466 +    }
467 + #else
468 +    if (storageLayout_ & DataStorage::dslElectroFrame) {
469 +      idat.eFrame1 = &(snap_->atomData.electroFrame[atom1]);
470 +      idat.eFrame2 = &(snap_->atomData.electroFrame[atom2]);
471 +    }
472 +    if (storageLayout_ & DataStorage::dslTorque) {
473 +      idat.t1 = &(snap_->atomData.torque[atom1]);
474 +      idat.t2 = &(snap_->atomData.torque[atom2]);
475 +    }
476 +    if (storageLayout_ & DataStorage::dslForce) {
477 +      idat.t1 = &(snap_->atomData.force[atom1]);
478 +      idat.t2 = &(snap_->atomData.force[atom2]);
479 +    }
480 + #endif
481 +    
482 +  }
483 +
484 +
485 +
486 +
487 +  /*
488 +   * buildNeighborList
489 +   *
490 +   * first element of pair is row-indexed CutoffGroup
491 +   * second element of pair is column-indexed CutoffGroup
492 +   */
493 +  vector<pair<int, int> > ForceMatrixDecomposition::buildNeighborList() {
494 +      
495 +    vector<pair<int, int> > neighborList;
496 + #ifdef IS_MPI
497 +    cellListRow_.clear();
498 +    cellListCol_.clear();
499 + #else
500 +    cellList_.clear();
501 + #endif
502 +
503 +    // dangerous to not do error checking.
504 +    RealType rCut_;
505 +
506 +    RealType rList_ = (rCut_ + skinThickness_);
507 +    RealType rl2 = rList_ * rList_;
508 +    Snapshot* snap_ = sman_->getCurrentSnapshot();
509 +    Mat3x3d Hmat = snap_->getHmat();
510 +    Vector3d Hx = Hmat.getColumn(0);
511 +    Vector3d Hy = Hmat.getColumn(1);
512 +    Vector3d Hz = Hmat.getColumn(2);
513 +
514 +    nCells_.x() = (int) ( Hx.length() )/ rList_;
515 +    nCells_.y() = (int) ( Hy.length() )/ rList_;
516 +    nCells_.z() = (int) ( Hz.length() )/ rList_;
517 +
518 +    Mat3x3d invHmat = snap_->getInvHmat();
519 +    Vector3d rs, scaled, dr;
520 +    Vector3i whichCell;
521 +    int cellIndex;
522 +
523 + #ifdef IS_MPI
524 +    for (int i = 0; i < nGroupsInRow_; i++) {
525 +      rs = cgRowData.position[i];
526 +      // scaled positions relative to the box vectors
527 +      scaled = invHmat * rs;
528 +      // wrap the vector back into the unit box by subtracting integer box
529 +      // numbers
530 +      for (int j = 0; j < 3; j++)
531 +        scaled[j] -= roundMe(scaled[j]);
532 +    
533 +      // find xyz-indices of cell that cutoffGroup is in.
534 +      whichCell.x() = nCells_.x() * scaled.x();
535 +      whichCell.y() = nCells_.y() * scaled.y();
536 +      whichCell.z() = nCells_.z() * scaled.z();
537 +
538 +      // find single index of this cell:
539 +      cellIndex = Vlinear(whichCell, nCells_);
540 +      // add this cutoff group to the list of groups in this cell;
541 +      cellListRow_[cellIndex].push_back(i);
542 +    }
543 +
544 +    for (int i = 0; i < nGroupsInCol_; i++) {
545 +      rs = cgColData.position[i];
546 +      // scaled positions relative to the box vectors
547 +      scaled = invHmat * rs;
548 +      // wrap the vector back into the unit box by subtracting integer box
549 +      // numbers
550 +      for (int j = 0; j < 3; j++)
551 +        scaled[j] -= roundMe(scaled[j]);
552 +
553 +      // find xyz-indices of cell that cutoffGroup is in.
554 +      whichCell.x() = nCells_.x() * scaled.x();
555 +      whichCell.y() = nCells_.y() * scaled.y();
556 +      whichCell.z() = nCells_.z() * scaled.z();
557 +
558 +      // find single index of this cell:
559 +      cellIndex = Vlinear(whichCell, nCells_);
560 +      // add this cutoff group to the list of groups in this cell;
561 +      cellListCol_[cellIndex].push_back(i);
562 +    }
563 + #else
564 +    for (int i = 0; i < nGroups_; i++) {
565 +      rs = snap_->cgData.position[i];
566 +      // scaled positions relative to the box vectors
567 +      scaled = invHmat * rs;
568 +      // wrap the vector back into the unit box by subtracting integer box
569 +      // numbers
570 +      for (int j = 0; j < 3; j++)
571 +        scaled[j] -= roundMe(scaled[j]);
572 +
573 +      // find xyz-indices of cell that cutoffGroup is in.
574 +      whichCell.x() = nCells_.x() * scaled.x();
575 +      whichCell.y() = nCells_.y() * scaled.y();
576 +      whichCell.z() = nCells_.z() * scaled.z();
577 +
578 +      // find single index of this cell:
579 +      cellIndex = Vlinear(whichCell, nCells_);
580 +      // add this cutoff group to the list of groups in this cell;
581 +      cellList_[cellIndex].push_back(i);
582 +    }
583 + #endif
584 +
585 +
586 +
587 +    for (int m1z = 0; m1z < nCells_.z(); m1z++) {
588 +      for (int m1y = 0; m1y < nCells_.y(); m1y++) {
589 +        for (int m1x = 0; m1x < nCells_.x(); m1x++) {
590 +          Vector3i m1v(m1x, m1y, m1z);
591 +          int m1 = Vlinear(m1v, nCells_);
592 +
593 +          for (vector<Vector3i>::iterator os = cellOffsets_.begin();
594 +               os != cellOffsets_.end(); ++os) {
595 +            
596 +            Vector3i m2v = m1v + (*os);
597 +            
598 +            if (m2v.x() >= nCells_.x()) {
599 +              m2v.x() = 0;          
600 +            } else if (m2v.x() < 0) {
601 +              m2v.x() = nCells_.x() - 1;
602 +            }
603 +            
604 +            if (m2v.y() >= nCells_.y()) {
605 +              m2v.y() = 0;          
606 +            } else if (m2v.y() < 0) {
607 +              m2v.y() = nCells_.y() - 1;
608 +            }
609 +            
610 +            if (m2v.z() >= nCells_.z()) {
611 +              m2v.z() = 0;          
612 +            } else if (m2v.z() < 0) {
613 +              m2v.z() = nCells_.z() - 1;
614 +            }
615 +            
616 +            int m2 = Vlinear (m2v, nCells_);
617 +
618 + #ifdef IS_MPI
619 +            for (vector<int>::iterator j1 = cellListRow_[m1].begin();
620 +                 j1 != cellListRow_[m1].end(); ++j1) {
621 +              for (vector<int>::iterator j2 = cellListCol_[m2].begin();
622 +                   j2 != cellListCol_[m2].end(); ++j2) {
623 +                              
624 +                // Always do this if we're in different cells or if
625 +                // we're in the same cell and the global index of the
626 +                // j2 cutoff group is less than the j1 cutoff group
627 +
628 +                if (m2 != m1 || cgColToGlobal[(*j2)] < cgRowToGlobal[(*j1)]) {
629 +                  dr = cgColData.position[(*j2)] - cgRowData.position[(*j1)];
630 +                  snap_->wrapVector(dr);
631 +                  if (dr.lengthSquare() < rl2) {
632 +                    neighborList.push_back(make_pair((*j1), (*j2)));
633 +                  }
634 +                }
635 +              }
636 +            }
637 + #else
638 +            for (vector<int>::iterator j1 = cellList_[m1].begin();
639 +                 j1 != cellList_[m1].end(); ++j1) {
640 +              for (vector<int>::iterator j2 = cellList_[m2].begin();
641 +                   j2 != cellList_[m2].end(); ++j2) {
642 +                              
643 +                // Always do this if we're in different cells or if
644 +                // we're in the same cell and the global index of the
645 +                // j2 cutoff group is less than the j1 cutoff group
646 +
647 +                if (m2 != m1 || (*j2) < (*j1)) {
648 +                  dr = snap_->cgData.position[(*j2)] - snap_->cgData.position[(*j1)];
649 +                  snap_->wrapVector(dr);
650 +                  if (dr.lengthSquare() < rl2) {
651 +                    neighborList.push_back(make_pair((*j1), (*j2)));
652 +                  }
653 +                }
654 +              }
655 +            }
656 + #endif
657 +          }
658 +        }
659 +      }
660 +    }
661 +
662 +    // save the local cutoff group positions for the check that is
663 +    // done on each loop:
664 +    saved_CG_positions_.clear();
665 +    for (int i = 0; i < nGroups_; i++)
666 +      saved_CG_positions_.push_back(snap_->cgData.position[i]);
667 +
668 +    return neighborList;
669 +  }
670   } //end namespace OpenMD

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines