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

Comparing trunk/src/selection/DistanceFinder.cpp (file contents):
Revision 1801 by gezelter, Mon Oct 1 18:21:15 2012 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

# 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   */
42  
43 + #ifdef IS_MPI
44 + #include <mpi.h>
45 + #endif
46 +
47   #include "selection/DistanceFinder.hpp"
48   #include "primitives/Molecule.hpp"
45 – namespace OpenMD {
49  
50 + namespace OpenMD {
51 +  
52    DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) {
53 +    nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies());
54 +    nObjects_.push_back(info_->getNGlobalBonds());
55 +    nObjects_.push_back(info_->getNGlobalBends());
56 +    nObjects_.push_back(info_->getNGlobalTorsions());
57 +    nObjects_.push_back(info_->getNGlobalInversions());
58  
59 <    nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies();
60 <    stuntdoubles_.resize(nStuntDoubles_);
59 >    stuntdoubles_.resize(nObjects_[STUNTDOUBLE]);
60 >    bonds_.resize(nObjects_[BOND]);
61 >    bends_.resize(nObjects_[BEND]);
62 >    torsions_.resize(nObjects_[TORSION]);
63 >    inversions_.resize(nObjects_[INVERSION]);
64      
65      SimInfo::MoleculeIterator mi;
53 –    Molecule* mol;
66      Molecule::AtomIterator ai;
55 –    Atom* atom;
67      Molecule::RigidBodyIterator rbIter;
68 <    RigidBody* rb;
68 >    Molecule::BondIterator bondIter;
69 >    Molecule::BendIterator bendIter;
70 >    Molecule::TorsionIterator torsionIter;
71 >    Molecule::InversionIterator inversionIter;
72  
73 +    Molecule* mol;
74 +    Atom* atom;
75 +    RigidBody* rb;
76 +    Bond* bond;
77 +    Bend* bend;
78 +    Torsion* torsion;
79 +    Inversion* inversion;    
80      
81      for (mol = info_->beginMolecule(mi); mol != NULL;
82           mol = info_->nextMolecule(mi)) {
# Line 64 | Line 85 | namespace OpenMD {
85            atom = mol->nextAtom(ai)) {
86          stuntdoubles_[atom->getGlobalIndex()] = atom;
87        }
67 –
88        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
89             rb = mol->nextRigidBody(rbIter)) {
90          stuntdoubles_[rb->getGlobalIndex()] = rb;
91        }
92 +      for (bond = mol->beginBond(bondIter); bond != NULL;
93 +           bond = mol->nextBond(bondIter)) {
94 +        bonds_[bond->getGlobalIndex()] = bond;
95 +      }  
96 +      for (bend = mol->beginBend(bendIter); bend != NULL;
97 +           bend = mol->nextBend(bendIter)) {
98 +        bends_[bend->getGlobalIndex()] = bend;
99 +      }  
100 +      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
101 +           torsion = mol->nextTorsion(torsionIter)) {
102 +        torsions_[torsion->getGlobalIndex()] = torsion;
103 +      }  
104 +      for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
105 +           inversion = mol->nextInversion(inversionIter)) {
106 +        inversions_[inversion->getGlobalIndex()] = inversion;
107 +      }  
108 +
109      }
110    }
111  
112 <  OpenMDBitSet DistanceFinder::find(const OpenMDBitSet& bs, RealType distance) {
112 >  SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance) {
113      StuntDouble * center;
114      Vector3d centerPos;
115      Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
116 <    OpenMDBitSet bsResult(nStuntDoubles_);
116 >    SelectionSet bsResult(nObjects_);  
117      assert(bsResult.size() == bs.size());
118 <  
118 >
119 > #ifdef IS_MPI
120 >    int mol;
121 >    int proc;
122 >    RealType data[3];
123 >    int worldRank = MPI::COMM_WORLD.Get_rank();
124 > #endif
125 >
126      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
127        if (stuntdoubles_[j]->isRigidBody()) {
128          RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
129          rb->updateAtoms();
130        }
131      }
132 <  
133 <    // This will fail in parallel because i might not be on this processor.
132 >            
133 >    SelectionSet bsTemp(nObjects_);
134 >    bsTemp = bs;
135 >    bsTemp.parallelReduce();
136  
137 <    for (int i = bs.firstOnBit(); i != -1; i = bs.nextOnBit(i)) {
137 >    for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1;
138 >         i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) {
139 >      
140 >      // Now, if we own stuntdouble i, we can use the position, but in
141 >      // parallel, we'll need to let everyone else know what that
142 >      // position is!
143 >      
144 > #ifdef IS_MPI
145 >      mol = info_->getGlobalMolMembership(i);
146 >      proc = info_->getMolToProc(mol);
147 >    
148 >      if (proc == worldRank) {
149 >        center = stuntdoubles_[i];
150 >        centerPos = center->getPos();
151 >        data[0] = centerPos.x();
152 >        data[1] = centerPos.y();
153 >        data[2] = centerPos.z();          
154 >        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
155 >      } else {
156 >        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
157 >        centerPos = Vector3d(data);
158 >      }
159 > #else
160        center = stuntdoubles_[i];
161        centerPos = center->getPos();
162 + #endif
163 +      
164        for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
165          Vector3d r =centerPos - stuntdoubles_[j]->getPos();
166          currSnapshot->wrapVector(r);
167          if (r.length() <= distance) {
168 <          bsResult.setBitOn(j);
168 >          bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
169          }
170        }
171 +      for (unsigned int j = 0; j < bonds_.size(); ++j) {
172 +        Vector3d loc = bonds_[j]->getAtomA()->getPos();
173 +        loc += bonds_[j]->getAtomB()->getPos();
174 +        loc = loc / 2.0;
175 +        Vector3d r = centerPos - loc;
176 +        currSnapshot->wrapVector(r);
177 +        if (r.length() <= distance) {
178 +          bsResult.bitsets_[BOND].setBitOn(j);
179 +        }
180 +      }
181 +      for (unsigned int j = 0; j < bends_.size(); ++j) {
182 +        Vector3d loc = bends_[j]->getAtomA()->getPos();
183 +        loc += bends_[j]->getAtomB()->getPos();
184 +        loc += bends_[j]->getAtomC()->getPos();
185 +        loc = loc / 3.0;
186 +        Vector3d r = centerPos - loc;
187 +        currSnapshot->wrapVector(r);
188 +        if (r.length() <= distance) {
189 +          bsResult.bitsets_[BEND].setBitOn(j);
190 +        }
191 +      }
192 +      for (unsigned int j = 0; j < torsions_.size(); ++j) {
193 +        Vector3d loc = torsions_[j]->getAtomA()->getPos();
194 +        loc += torsions_[j]->getAtomB()->getPos();
195 +        loc += torsions_[j]->getAtomC()->getPos();
196 +        loc += torsions_[j]->getAtomD()->getPos();
197 +        loc = loc / 4.0;
198 +        Vector3d r = centerPos - loc;
199 +        currSnapshot->wrapVector(r);
200 +        if (r.length() <= distance) {
201 +          bsResult.bitsets_[TORSION].setBitOn(j);
202 +        }
203 +      }
204 +      for (unsigned int j = 0; j < inversions_.size(); ++j) {
205 +        Vector3d loc = inversions_[j]->getAtomA()->getPos();
206 +        loc += inversions_[j]->getAtomB()->getPos();
207 +        loc += inversions_[j]->getAtomC()->getPos();
208 +        loc += inversions_[j]->getAtomD()->getPos();
209 +        loc = loc / 4.0;
210 +        Vector3d r = centerPos - loc;
211 +        currSnapshot->wrapVector(r);
212 +        if (r.length() <= distance) {
213 +          bsResult.bitsets_[INVERSION].setBitOn(j);
214 +        }
215 +      }
216      }
102 –
217      return bsResult;
218    }
219 +  
220 +  
221 +  SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance, int frame ) {
222 +    StuntDouble * center;
223 +    Vector3d centerPos;
224 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame);
225 +    SelectionSet bsResult(nObjects_);  
226 +    assert(bsResult.size() == bs.size());
227  
228 + #ifdef IS_MPI
229 +    int mol;
230 +    int proc;
231 +    RealType data[3];
232 +    int worldRank = MPI::COMM_WORLD.Get_rank();
233 + #endif
234 +
235 +    for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
236 +      if (stuntdoubles_[j]->isRigidBody()) {
237 +        RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
238 +        rb->updateAtoms(frame);
239 +      }
240 +    }
241 +      
242 +    SelectionSet bsTemp(nObjects_);
243 +    bsTemp = bs;
244 +    bsTemp.parallelReduce();
245 +
246 +    for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1;
247 +         i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) {
248 +
249 +      // Now, if we own stuntdouble i, we can use the position, but in
250 +      // parallel, we'll need to let everyone else know what that
251 +      // position is!
252 +
253 + #ifdef IS_MPI
254 +      mol = info_->getGlobalMolMembership(i);
255 +      proc = info_->getMolToProc(mol);
256 +    
257 +      if (proc == worldRank) {
258 +        center = stuntdoubles_[i];
259 +        centerPos = center->getPos(frame);
260 +        data[0] = centerPos.x();
261 +        data[1] = centerPos.y();
262 +        data[2] = centerPos.z();          
263 +        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
264 +      } else {
265 +        MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
266 +        centerPos = Vector3d(data);
267 +      }
268 + #else
269 +      center = stuntdoubles_[i];
270 +      centerPos = center->getPos(frame);
271 + #endif
272 +      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
273 +        Vector3d r =centerPos - stuntdoubles_[j]->getPos(frame);
274 +        currSnapshot->wrapVector(r);
275 +        if (r.length() <= distance) {
276 +          bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
277 +        }
278 +      }
279 +      for (unsigned int j = 0; j < bonds_.size(); ++j) {      
280 +        Vector3d loc = bonds_[j]->getAtomA()->getPos(frame);
281 +        loc += bonds_[j]->getAtomB()->getPos(frame);
282 +        loc = loc / 2.0;
283 +        Vector3d r = centerPos - loc;
284 +        currSnapshot->wrapVector(r);
285 +        if (r.length() <= distance) {
286 +          bsResult.bitsets_[BOND].setBitOn(j);
287 +        }
288 +      }
289 +      for (unsigned int j = 0; j < bends_.size(); ++j) {
290 +        Vector3d loc = bends_[j]->getAtomA()->getPos(frame);
291 +        loc += bends_[j]->getAtomB()->getPos(frame);
292 +        loc += bends_[j]->getAtomC()->getPos(frame);
293 +        loc = loc / 3.0;
294 +        Vector3d r = centerPos - loc;
295 +        currSnapshot->wrapVector(r);
296 +        if (r.length() <= distance) {
297 +          bsResult.bitsets_[BEND].setBitOn(j);
298 +        }
299 +      }
300 +      for (unsigned int j = 0; j < torsions_.size(); ++j) {
301 +        Vector3d loc = torsions_[j]->getAtomA()->getPos(frame);
302 +        loc += torsions_[j]->getAtomB()->getPos(frame);
303 +        loc += torsions_[j]->getAtomC()->getPos(frame);
304 +        loc += torsions_[j]->getAtomD()->getPos(frame);
305 +        loc = loc / 4.0;
306 +        Vector3d r = centerPos - loc;
307 +        currSnapshot->wrapVector(r);
308 +        if (r.length() <= distance) {
309 +          bsResult.bitsets_[TORSION].setBitOn(j);
310 +        }
311 +      }
312 +      for (unsigned int j = 0; j < inversions_.size(); ++j) {
313 +        Vector3d loc = inversions_[j]->getAtomA()->getPos(frame);
314 +        loc += inversions_[j]->getAtomB()->getPos(frame);
315 +        loc += inversions_[j]->getAtomC()->getPos(frame);
316 +        loc += inversions_[j]->getAtomD()->getPos(frame);
317 +        loc = loc / 4.0;
318 +        Vector3d r = centerPos - loc;
319 +        currSnapshot->wrapVector(r);
320 +        if (r.length() <= distance) {
321 +          bsResult.bitsets_[INVERSION].setBitOn(j);
322 +        }
323 +      }
324 +    }
325 +    return bsResult;
326 +  }
327   }

Diff Legend

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