ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/selection/DistanceFinder.cpp
(Generate patch)

Comparing branches/development/src/selection/DistanceFinder.cpp (file contents):
Revision 1465 by chuckv, Fri Jul 9 23:08:25 2010 UTC vs.
Revision 1831 by gezelter, Thu Jan 10 14:06:34 2013 UTC

# Line 36 | Line 36
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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include "selection/DistanceFinder.hpp"
44   #include "primitives/Molecule.hpp"
45 < namespace OpenMD {
45 > #ifdef IS_MPI
46 > #include <mpi.h>
47 > #endif
48  
49 + namespace OpenMD {
50 +  
51    DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) {
52 <
52 >    
53      nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies();
54      stuntdoubles_.resize(nStuntDoubles_);
55      
# Line 54 | Line 59 | namespace OpenMD {
59      Atom* atom;
60      Molecule::RigidBodyIterator rbIter;
61      RigidBody* rb;
57 –
62      
63 <    for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
63 >    
64 >    for (mol = info_->beginMolecule(mi); mol != NULL;
65 >         mol = info_->nextMolecule(mi)) {
66          
67 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
67 >      for(atom = mol->beginAtom(ai); atom != NULL;
68 >          atom = mol->nextAtom(ai)) {
69          stuntdoubles_[atom->getGlobalIndex()] = atom;
70        }
71  
72 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
72 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
73 >           rb = mol->nextRigidBody(rbIter)) {
74          stuntdoubles_[rb->getGlobalIndex()] = rb;
75        }
76 <        
69 <    }    
70 <
76 >    }
77    }
78  
79    OpenMDBitSet DistanceFinder::find(const OpenMDBitSet& bs, RealType distance) {
80      StuntDouble * center;
81      Vector3d centerPos;
82      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
83 <    OpenMDBitSet bsResult(nStuntDoubles_);
83 >    OpenMDBitSet bsResult(nStuntDoubles_);  
84      assert(bsResult.size() == bs.size());
85 <  
86 <    for (int j = 0; j < stuntdoubles_.size(); ++j) {
85 >
86 > #ifdef IS_MPI
87 >    int mol;
88 >    int proc;
89 >    RealType data[3];
90 >    int worldRank = MPI::COMM_WORLD.Get_rank();
91 > #endif
92 >
93 >    for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
94        if (stuntdoubles_[j]->isRigidBody()) {
95          RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
96          rb->updateAtoms();
97        }
98      }
99 <  
100 <    for (int i = bs.firstOnBit(); i != -1; i = bs.nextOnBit(i)) {
99 >      
100 >    OpenMDBitSet bsTemp(nStuntDoubles_);
101 >    bsTemp = bs;
102 >    bsTemp.parallelReduce();
103 >
104 >    for (int i = bsTemp.firstOnBit(); i != -1; i = bsTemp.nextOnBit(i)) {
105 >
106 >      // Now, if we own stuntdouble i, we can use the position, but in
107 >      // parallel, we'll need to let everyone else know what that
108 >      // position is!
109 >
110 > #ifdef IS_MPI
111 >      mol = info_->getGlobalMolMembership(i);
112 >      proc = info_->getMolToProc(mol);
113 >    
114 >      if (proc == worldRank) {
115 >        center = stuntdoubles_[i];
116 >        centerPos = center->getPos();
117 >        data[0] = centerPos.x();
118 >        data[1] = centerPos.y();
119 >        data[2] = centerPos.z();          
120 >        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
121 >      } else {
122 >        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
123 >        centerPos = Vector3d(data);
124 >      }
125 > #else
126        center = stuntdoubles_[i];
127        centerPos = center->getPos();
128 <      for (int j = 0; j < stuntdoubles_.size(); ++j) {
128 > #endif
129 >      
130 >      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
131          Vector3d r =centerPos - stuntdoubles_[j]->getPos();
132          currSnapshot->wrapVector(r);
133          if (r.length() <= distance) {
# Line 95 | Line 135 | namespace OpenMD {
135          }
136        }
137      }
98 –
138      return bsResult;
139    }
140  
141 +
142 + OpenMDBitSet DistanceFinder::find(const OpenMDBitSet& bs, RealType distance, int frame ) {
143 +    StuntDouble * center;
144 +    Vector3d centerPos;
145 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame);
146 +    OpenMDBitSet bsResult(nStuntDoubles_);  
147 +    assert(bsResult.size() == bs.size());
148 +
149 + #ifdef IS_MPI
150 +    int mol;
151 +    int proc;
152 +    RealType data[3];
153 +    int worldRank = MPI::COMM_WORLD.Get_rank();
154 + #endif
155 +
156 +    for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
157 +      if (stuntdoubles_[j]->isRigidBody()) {
158 +        RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
159 +        rb->updateAtoms(frame);
160 +      }
161 +    }
162 +      
163 +    OpenMDBitSet bsTemp(nStuntDoubles_);
164 +    bsTemp = bs;
165 +    bsTemp.parallelReduce();
166 +
167 +    for (int i = bsTemp.firstOnBit(); i != -1; i = bsTemp.nextOnBit(i)) {
168 +
169 +      // Now, if we own stuntdouble i, we can use the position, but in
170 +      // parallel, we'll need to let everyone else know what that
171 +      // position is!
172 +
173 + #ifdef IS_MPI
174 +      mol = info_->getGlobalMolMembership(i);
175 +      proc = info_->getMolToProc(mol);
176 +    
177 +      if (proc == worldRank) {
178 +        center = stuntdoubles_[i];
179 +        centerPos = center->getPos(frame);
180 +        data[0] = centerPos.x();
181 +        data[1] = centerPos.y();
182 +        data[2] = centerPos.z();          
183 +        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
184 +      } else {
185 +        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
186 +        centerPos = Vector3d(data);
187 +      }
188 + #else
189 +      center = stuntdoubles_[i];
190 +      centerPos = center->getPos(frame);
191 + #endif
192 +      
193 +      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
194 +        Vector3d r =centerPos - stuntdoubles_[j]->getPos(frame);
195 +        currSnapshot->wrapVector(r);
196 +        if (r.length() <= distance) {
197 +          bsResult.setBitOn(j);
198 +        }
199 +      }
200 +    }
201 +    return bsResult;
202 +  }
203   }

Diff Legend

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