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 1938 by gezelter, Thu Oct 31 15:32:17 2013 UTC vs.
Revision 2052 by gezelter, Fri Jan 9 19:06:35 2015 UTC

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

Diff Legend

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