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 1816 by gezelter, Fri Dec 7 18:33:33 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 – #include "selection/DistanceFinder.hpp"
44 – #include "primitives/Molecule.hpp"
43   #ifdef IS_MPI
44   #include <mpi.h>
45   #endif
46  
47 + #include "selection/DistanceFinder.hpp"
48 + #include "primitives/Molecule.hpp"
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 +    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      
53 –    nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies();
54 –    stuntdoubles_.resize(nStuntDoubles_);
55 –    
65      SimInfo::MoleculeIterator mi;
57 –    Molecule* mol;
66      Molecule::AtomIterator ai;
59 –    Atom* atom;
67      Molecule::RigidBodyIterator rbIter;
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      
63 –    
81      for (mol = info_->beginMolecule(mi); mol != NULL;
82           mol = info_->nextMolecule(mi)) {
83          
# Line 68 | Line 85 | namespace OpenMD {
85            atom = mol->nextAtom(ai)) {
86          stuntdoubles_[atom->getGlobalIndex()] = atom;
87        }
71 –
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  
119   #ifdef IS_MPI
# Line 96 | Line 129 | namespace OpenMD {
129          rb->updateAtoms();
130        }
131      }
132 <      
133 <    OpenMDBitSet bsTemp(nStuntDoubles_);
132 >            
133 >    SelectionSet bsTemp(nObjects_);
134      bsTemp = bs;
135      bsTemp.parallelReduce();
136  
137 <    for (int i = bsTemp.firstOnBit(); i != -1; i = bsTemp.nextOnBit(i)) {
138 <
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 <
143 >      
144   #ifdef IS_MPI
145        mol = info_->getGlobalMolMembership(i);
146        proc = info_->getMolToProc(mol);
# Line 131 | Line 165 | namespace OpenMD {
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      }
217      return bsResult;
218    }
219 <
220 <
221 < OpenMDBitSet DistanceFinder::find(const OpenMDBitSet& bs, RealType distance, int frame ) {
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 <    OpenMDBitSet bsResult(nStuntDoubles_);  
225 >    SelectionSet bsResult(nObjects_);  
226      assert(bsResult.size() == bs.size());
227  
228   #ifdef IS_MPI
# Line 160 | Line 239 | OpenMDBitSet DistanceFinder::find(const OpenMDBitSet&
239        }
240      }
241        
242 <    OpenMDBitSet bsTemp(nStuntDoubles_);
242 >    SelectionSet bsTemp(nObjects_);
243      bsTemp = bs;
244      bsTemp.parallelReduce();
245  
246 <    for (int i = bsTemp.firstOnBit(); i != -1; i = bsTemp.nextOnBit(i)) {
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
# Line 189 | Line 269 | OpenMDBitSet DistanceFinder::find(const OpenMDBitSet&
269        center = stuntdoubles_[i];
270        centerPos = center->getPos(frame);
271   #endif
192 –      
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.setBitOn(j);
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    }

Diff Legend

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