ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/selection/DistanceFinder.cpp

Comparing trunk/src/selection/DistanceFinder.cpp (file contents):
Revision 770 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC

#	Line 1 | Line 1
1		/*
2	<	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
2	>	* Copyright (c) 2005, 2010 The University of Notre Dame. All Rights Reserved.
3		*
4		* The University of Notre Dame grants you ("Licensee") a
5		* non-exclusive, royalty free, license to use, modify and
6		* redistribute this software in source and binary code form, provided
7		* that the following conditions are met:
8		*
9	<	* 1. Acknowledgement of the program authors must be made in any
10	<	*    publication of scientific results based in part on use of the
11	<	*    program.  An acceptable form of acknowledgement is citation of
12	<	*    the article in which the program was described (Matthew
13	<	*    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	<	*    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	<	*    Parallel Simulation Engine for Molecular Dynamics,"
16	<	*    J. Comput. Chem. 26, pp. 252-271 (2005))
17	<	*
18	<	* 2. Redistributions of source code must retain the above copyright
9	>	* 1. Redistributions of source code must retain the above copyright
10		*    notice, this list of conditions and the following disclaimer.
11		*
12	<	* 3. Redistributions in binary form must reproduce the above copyright
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13		*    notice, this list of conditions and the following disclaimer in the
14		*    documentation and/or other materials provided with the
15		*    distribution.
#	Line 37 | Line 28
28		* arising out of the use of or inability to use software, even if the
29		* University of Notre Dame has been advised of the possibility of
30		* such damages.
31	+	*
32	+	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	+	* research, please cite the appropriate papers when you publish your
34	+	* work.  Good starting points are:
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, 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 oopse {
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	<	OOPSEBitSet DistanceFinder::find(const OOPSEBitSet& bs, double distance) {
112	>	SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance) {
113		StuntDouble * center;
114		Vector3d centerPos;
115		Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
116	<	OOPSEBitSet bsResult(nStuntDoubles_);
116	>	SelectionSet bsResult(nObjects_);
117		assert(bsResult.size() == bs.size());
118	<
119	<	for (int i = bs.firstOnBit(); i != -1; i = bs.nextOnBit(i)) {
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	>	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		}
91	–
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		}

Comparing trunk/src/selection/DistanceFinder.cpp (property svn:keywords):
Revision 770 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 1969 by gezelter, Wed Feb 26 14:14:50 2014 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

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