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 2052 by gezelter, Fri Jan 9 19:06:35 2015 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	+	nObjects_.push_back(info_->getNGlobalMolecules());
59
60	<	nStuntDoubles_ = info_->getNGlobalAtoms() + info_->getNGlobalRigidBodies();
61	<	stuntdoubles_.resize(nStuntDoubles_);
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
67		SimInfo::MoleculeIterator mi;
52	–	Molecule* mol;
68		Molecule::AtomIterator ai;
54	–	Atom* atom;
69		Molecule::RigidBodyIterator rbIter;
70	<	RigidBody* rb;
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
83	<	for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) {
84	<
85	<	for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
83	>	for (mol = info_->beginMolecule(mi); mol != NULL;
84	>	mol = info_->nextMolecule(mi)) {
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		}
92	<
93	<	for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
92	>	for (rb = mol->beginRigidBody(rbIter); rb != NULL;
93	>	rb = mol->nextRigidBody(rbIter)) {
94		stuntdoubles_[rb->getGlobalIndex()] = rb;
95		}
96	<
97	<	}
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	<	OOPSEBitSet DistanceFinder::find(const OOPSEBitSet& bs, double distance) {
116	>	SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance) {
117		StuntDouble * center;
118		Vector3d centerPos;
119		Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
120	<	OOPSEBitSet bsResult(nStuntDoubles_);
120	>	SelectionSet bsResult(nObjects_);
121		assert(bsResult.size() == bs.size());
122	<
123	<	for (int i = bs.firstOnBit(); i != -1; i = bs.nextOnBit(i)) {
122	>
123	>	#ifdef IS_MPI
124	>	int mol;
125	>	int proc;
126	>	RealType data[3];
127	>	int worldRank;
128	>	MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
129	>	#endif
130	>
131	>	for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
132	>	if (stuntdoubles_[j]->isRigidBody()) {
133	>	RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
134	>	rb->updateAtoms();
135	>	}
136	>	}
137	>
138	>	SelectionSet bsTemp(nObjects_);
139	>	bsTemp = bs;
140	>	bsTemp.parallelReduce();
141	>
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	>
149	>	#ifdef IS_MPI
150	>	mol = info_->getGlobalMolMembership(i);
151	>	proc = info_->getMolToProc(mol);
152	>
153	>	if (proc == worldRank) {
154	>	center = stuntdoubles_[i];
155	>	centerPos = center->getPos();
156	>	data[0] = centerPos.x();
157	>	data[1] = centerPos.y();
158	>	data[2] = centerPos.z();
159	>	MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
160	>	} else {
161	>	MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
162	>	centerPos = Vector3d(data);
163	>	}
164	>	#else
165		center = stuntdoubles_[i];
166		centerPos = center->getPos();
167	<	for (int j = 0; j < stuntdoubles_.size(); ++j) {
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		}
91	–
229		return bsResult;
230		}
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	+	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;
245	+	MPI_Comm_rank( MPI_COMM_WORLD, &worldRank);
246	+	#endif
247	+
248	+	for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
249	+	if (stuntdoubles_[j]->isRigidBody()) {
250	+	RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
251	+	rb->updateAtoms(frame);
252	+	}
253	+	}
254	+
255	+	SelectionSet bsTemp(nObjects_);
256	+	bsTemp = bs;
257	+	bsTemp.parallelReduce();
258	+
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
264	+	// position is!
265	+
266	+	#ifdef IS_MPI
267	+	mol = info_->getGlobalMolMembership(i);
268	+	proc = info_->getMolToProc(mol);
269	+
270	+	if (proc == worldRank) {
271	+	center = stuntdoubles_[i];
272	+	centerPos = center->getPos(frame);
273	+	data[0] = centerPos.x();
274	+	data[1] = centerPos.y();
275	+	data[2] = centerPos.z();
276	+	MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
277	+	} else {
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.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	+	}
348		}

Comparing trunk/src/selection/DistanceFinder.cpp (property svn:keywords):
Revision 770 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 2052 by gezelter, Fri Jan 9 19:06:35 2015 UTC

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

Diff Legend

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