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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 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).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
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"
44 – 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;
52 –    Molecule* mol;
66      Molecule::AtomIterator ai;
54 –    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; mol = info_->nextMolecule(mi)) {
81 >    for (mol = info_->beginMolecule(mi); mol != NULL;
82 >         mol = info_->nextMolecule(mi)) {
83          
84 <      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
84 >      for(atom = mol->beginAtom(ai); atom != NULL;
85 >          atom = mol->nextAtom(ai)) {
86          stuntdoubles_[atom->getGlobalIndex()] = atom;
87        }
88 <
89 <      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
88 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
89 >           rb = mol->nextRigidBody(rbIter)) {
90          stuntdoubles_[rb->getGlobalIndex()] = rb;
91        }
92 <        
93 <    }    
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 <  
119 <    for (int j = 0; j < stuntdoubles_.size(); ++j) {
118 >
119 > #ifdef IS_MPI
120 >    int mol;
121 >    int proc;
122 >    RealType data[3];
123 >    int worldRank;
124 >    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
125 > #endif
126 >
127 >    for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
128        if (stuntdoubles_[j]->isRigidBody()) {
129          RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
130          rb->updateAtoms();
131        }
132      }
133 <  
134 <    for (int i = bs.firstOnBit(); i != -1; i = bs.nextOnBit(i)) {
133 >            
134 >    SelectionSet bsTemp(nObjects_);
135 >    bsTemp = bs;
136 >    bsTemp.parallelReduce();
137 >
138 >    for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1;
139 >         i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) {
140 >      
141 >      // Now, if we own stuntdouble i, we can use the position, but in
142 >      // parallel, we'll need to let everyone else know what that
143 >      // position is!
144 >      
145 > #ifdef IS_MPI
146 >      mol = info_->getGlobalMolMembership(i);
147 >      proc = info_->getMolToProc(mol);
148 >    
149 >      if (proc == worldRank) {
150 >        center = stuntdoubles_[i];
151 >        centerPos = center->getPos();
152 >        data[0] = centerPos.x();
153 >        data[1] = centerPos.y();
154 >        data[2] = centerPos.z();          
155 >        MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
156 >      } else {
157 >        MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
158 >        centerPos = Vector3d(data);
159 >      }
160 > #else
161        center = stuntdoubles_[i];
162        centerPos = center->getPos();
163 <      for (int j = 0; j < stuntdoubles_.size(); ++j) {
163 > #endif
164 >      
165 >      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
166          Vector3d r =centerPos - stuntdoubles_[j]->getPos();
167          currSnapshot->wrapVector(r);
168          if (r.length() <= distance) {
169 <          bsResult.setBitOn(j);
169 >          bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
170          }
171        }
172 +      for (unsigned int j = 0; j < bonds_.size(); ++j) {
173 +        Vector3d loc = bonds_[j]->getAtomA()->getPos();
174 +        loc += bonds_[j]->getAtomB()->getPos();
175 +        loc = loc / 2.0;
176 +        Vector3d r = centerPos - loc;
177 +        currSnapshot->wrapVector(r);
178 +        if (r.length() <= distance) {
179 +          bsResult.bitsets_[BOND].setBitOn(j);
180 +        }
181 +      }
182 +      for (unsigned int j = 0; j < bends_.size(); ++j) {
183 +        Vector3d loc = bends_[j]->getAtomA()->getPos();
184 +        loc += bends_[j]->getAtomB()->getPos();
185 +        loc += bends_[j]->getAtomC()->getPos();
186 +        loc = loc / 3.0;
187 +        Vector3d r = centerPos - loc;
188 +        currSnapshot->wrapVector(r);
189 +        if (r.length() <= distance) {
190 +          bsResult.bitsets_[BEND].setBitOn(j);
191 +        }
192 +      }
193 +      for (unsigned int j = 0; j < torsions_.size(); ++j) {
194 +        Vector3d loc = torsions_[j]->getAtomA()->getPos();
195 +        loc += torsions_[j]->getAtomB()->getPos();
196 +        loc += torsions_[j]->getAtomC()->getPos();
197 +        loc += torsions_[j]->getAtomD()->getPos();
198 +        loc = loc / 4.0;
199 +        Vector3d r = centerPos - loc;
200 +        currSnapshot->wrapVector(r);
201 +        if (r.length() <= distance) {
202 +          bsResult.bitsets_[TORSION].setBitOn(j);
203 +        }
204 +      }
205 +      for (unsigned int j = 0; j < inversions_.size(); ++j) {
206 +        Vector3d loc = inversions_[j]->getAtomA()->getPos();
207 +        loc += inversions_[j]->getAtomB()->getPos();
208 +        loc += inversions_[j]->getAtomC()->getPos();
209 +        loc += inversions_[j]->getAtomD()->getPos();
210 +        loc = loc / 4.0;
211 +        Vector3d r = centerPos - loc;
212 +        currSnapshot->wrapVector(r);
213 +        if (r.length() <= distance) {
214 +          bsResult.bitsets_[INVERSION].setBitOn(j);
215 +        }
216 +      }
217      }
98 –
218      return bsResult;
219    }
220 +  
221 +  
222 +  SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance, int frame ) {
223 +    StuntDouble * center;
224 +    Vector3d centerPos;
225 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame);
226 +    SelectionSet bsResult(nObjects_);  
227 +    assert(bsResult.size() == bs.size());
228  
229 + #ifdef IS_MPI
230 +    int mol;
231 +    int proc;
232 +    RealType data[3];
233 +    int worldRank;
234 +    MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
235 + #endif
236 +
237 +    for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
238 +      if (stuntdoubles_[j]->isRigidBody()) {
239 +        RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
240 +        rb->updateAtoms(frame);
241 +      }
242 +    }
243 +      
244 +    SelectionSet bsTemp(nObjects_);
245 +    bsTemp = bs;
246 +    bsTemp.parallelReduce();
247 +
248 +    for (int i = bsTemp.bitsets_[STUNTDOUBLE].firstOnBit(); i != -1;
249 +         i = bsTemp.bitsets_[STUNTDOUBLE].nextOnBit(i)) {
250 +
251 +      // Now, if we own stuntdouble i, we can use the position, but in
252 +      // parallel, we'll need to let everyone else know what that
253 +      // position is!
254 +
255 + #ifdef IS_MPI
256 +      mol = info_->getGlobalMolMembership(i);
257 +      proc = info_->getMolToProc(mol);
258 +    
259 +      if (proc == worldRank) {
260 +        center = stuntdoubles_[i];
261 +        centerPos = center->getPos(frame);
262 +        data[0] = centerPos.x();
263 +        data[1] = centerPos.y();
264 +        data[2] = centerPos.z();          
265 +        MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
266 +      } else {
267 +        MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
268 +        centerPos = Vector3d(data);
269 +      }
270 + #else
271 +      center = stuntdoubles_[i];
272 +      centerPos = center->getPos(frame);
273 + #endif
274 +      for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
275 +        Vector3d r =centerPos - stuntdoubles_[j]->getPos(frame);
276 +        currSnapshot->wrapVector(r);
277 +        if (r.length() <= distance) {
278 +          bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
279 +        }
280 +      }
281 +      for (unsigned int j = 0; j < bonds_.size(); ++j) {      
282 +        Vector3d loc = bonds_[j]->getAtomA()->getPos(frame);
283 +        loc += bonds_[j]->getAtomB()->getPos(frame);
284 +        loc = loc / 2.0;
285 +        Vector3d r = centerPos - loc;
286 +        currSnapshot->wrapVector(r);
287 +        if (r.length() <= distance) {
288 +          bsResult.bitsets_[BOND].setBitOn(j);
289 +        }
290 +      }
291 +      for (unsigned int j = 0; j < bends_.size(); ++j) {
292 +        Vector3d loc = bends_[j]->getAtomA()->getPos(frame);
293 +        loc += bends_[j]->getAtomB()->getPos(frame);
294 +        loc += bends_[j]->getAtomC()->getPos(frame);
295 +        loc = loc / 3.0;
296 +        Vector3d r = centerPos - loc;
297 +        currSnapshot->wrapVector(r);
298 +        if (r.length() <= distance) {
299 +          bsResult.bitsets_[BEND].setBitOn(j);
300 +        }
301 +      }
302 +      for (unsigned int j = 0; j < torsions_.size(); ++j) {
303 +        Vector3d loc = torsions_[j]->getAtomA()->getPos(frame);
304 +        loc += torsions_[j]->getAtomB()->getPos(frame);
305 +        loc += torsions_[j]->getAtomC()->getPos(frame);
306 +        loc += torsions_[j]->getAtomD()->getPos(frame);
307 +        loc = loc / 4.0;
308 +        Vector3d r = centerPos - loc;
309 +        currSnapshot->wrapVector(r);
310 +        if (r.length() <= distance) {
311 +          bsResult.bitsets_[TORSION].setBitOn(j);
312 +        }
313 +      }
314 +      for (unsigned int j = 0; j < inversions_.size(); ++j) {
315 +        Vector3d loc = inversions_[j]->getAtomA()->getPos(frame);
316 +        loc += inversions_[j]->getAtomB()->getPos(frame);
317 +        loc += inversions_[j]->getAtomC()->getPos(frame);
318 +        loc += inversions_[j]->getAtomD()->getPos(frame);
319 +        loc = loc / 4.0;
320 +        Vector3d r = centerPos - loc;
321 +        currSnapshot->wrapVector(r);
322 +        if (r.length() <= distance) {
323 +          bsResult.bitsets_[INVERSION].setBitOn(j);
324 +        }
325 +      }
326 +    }
327 +    return bsResult;
328 +  }
329   }

Diff Legend

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