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 295 by tim, Mon Feb 7 19:13:18 2005 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 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	<	DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) {
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)) {
85	<	stuntdoubles_[atom->getGlobalIndex()] = atom;
86	<	}
84	>	for(atom = mol->beginAtom(ai); atom != NULL;
85	>	atom = mol->nextAtom(ai)) {
86	>	stuntdoubles_[atom->getGlobalIndex()] = atom;
87	>	}
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	<	for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
110	<	stuntdoubles_[rb->getGlobalIndex()] = rb;
67	<	}
68	<
69	<	}
109	>	}
110	>	}
111
112	<	}
72	<
73	<	BitSet DistanceFinder::find(const BitSet& 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	<	BitSet bsResult(nStuntDoubles_);
116	>	SelectionSet bsResult(nObjects_);
117		assert(bsResult.size() == bs.size());
118	<
119	<	for (int i = bs.nextOnBit(-1); i != -1; i = bs.nextOnBit(i)) {
118	>
119	>	#ifdef IS_MPI
120	>	int mol;
121	>	int proc;
122	>	RealType data[3];
123	>	int worldRank = MPI::COMM_WORLD.Get_rank();
124	>	#endif
125	>
126	>	for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
127	>	if (stuntdoubles_[j]->isRigidBody()) {
128	>	RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
129	>	rb->updateAtoms();
130	>	}
131	>	}
132	>
133	>	SelectionSet bsTemp(nObjects_);
134	>	bsTemp = bs;
135	>	bsTemp.parallelReduce();
136	>
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	>
144	>	#ifdef IS_MPI
145	>	mol = info_->getGlobalMolMembership(i);
146	>	proc = info_->getMolToProc(mol);
147	>
148	>	if (proc == worldRank) {
149		center = stuntdoubles_[i];
150		centerPos = center->getPos();
151	<	for (int j = 0; j < stuntdoubles_.size(); ++j) {
152	<	Vector3d r =centerPos - stuntdoubles_[j]->getPos();
153	<	currSnapshot->wrapVector(r);
154	<	if (r.length() <= distance) {
155	<	bsResult.setBitOn(j);
156	<	}
157	<	}
151	>	data[0] = centerPos.x();
152	>	data[1] = centerPos.y();
153	>	data[2] = centerPos.z();
154	>	MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
155	>	} else {
156	>	MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
157	>	centerPos = Vector3d(data);
158	>	}
159	>	#else
160	>	center = stuntdoubles_[i];
161	>	centerPos = center->getPos();
162	>	#endif
163	>
164	>	for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
165	>	Vector3d r =centerPos - stuntdoubles_[j]->getPos();
166	>	currSnapshot->wrapVector(r);
167	>	if (r.length() <= distance) {
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	+	SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance, int frame ) {
222	+	StuntDouble * center;
223	+	Vector3d centerPos;
224	+	Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame);
225	+	SelectionSet bsResult(nObjects_);
226	+	assert(bsResult.size() == bs.size());
227
228	+	#ifdef IS_MPI
229	+	int mol;
230	+	int proc;
231	+	RealType data[3];
232	+	int worldRank = MPI::COMM_WORLD.Get_rank();
233	+	#endif
234	+
235	+	for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
236	+	if (stuntdoubles_[j]->isRigidBody()) {
237	+	RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
238	+	rb->updateAtoms(frame);
239	+	}
240	+	}
241	+
242	+	SelectionSet bsTemp(nObjects_);
243	+	bsTemp = bs;
244	+	bsTemp.parallelReduce();
245	+
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
251	+	// position is!
252	+
253	+	#ifdef IS_MPI
254	+	mol = info_->getGlobalMolMembership(i);
255	+	proc = info_->getMolToProc(mol);
256	+
257	+	if (proc == worldRank) {
258	+	center = stuntdoubles_[i];
259	+	centerPos = center->getPos(frame);
260	+	data[0] = centerPos.x();
261	+	data[1] = centerPos.y();
262	+	data[2] = centerPos.z();
263	+	MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
264	+	} else {
265	+	MPI::COMM_WORLD.Bcast(data, 3, MPI::REALTYPE, proc);
266	+	centerPos = Vector3d(data);
267	+	}
268	+	#else
269	+	center = stuntdoubles_[i];
270	+	centerPos = center->getPos(frame);
271	+	#endif
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.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	+	}
327		}
94	–
95	–	}

Comparing trunk/src/selection/DistanceFinder.cpp (property svn:keywords):
Revision 295 by tim, Mon Feb 7 19:13:18 2005 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

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

Diff Legend

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