OpenMD 3.1
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
DistanceFinder.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "selection/DistanceFinder.hpp"
46
47#ifdef IS_MPI
48#include <mpi.h>
49#endif
50
51#include "primitives/Molecule.hpp"
52
53namespace OpenMD {
54
55 DistanceFinder::DistanceFinder(SimInfo* info) : info_(info) {
56 nObjects_.push_back(info_->getNGlobalAtoms() +
57 info_->getNGlobalRigidBodies());
58 nObjects_.push_back(info_->getNGlobalBonds());
59 nObjects_.push_back(info_->getNGlobalBends());
60 nObjects_.push_back(info_->getNGlobalTorsions());
61 nObjects_.push_back(info_->getNGlobalInversions());
62 nObjects_.push_back(info_->getNGlobalMolecules());
63
64 stuntdoubles_.resize(nObjects_[STUNTDOUBLE]);
65 bonds_.resize(nObjects_[BOND]);
66 bends_.resize(nObjects_[BEND]);
67 torsions_.resize(nObjects_[TORSION]);
68 inversions_.resize(nObjects_[INVERSION]);
69 molecules_.resize(nObjects_[MOLECULE]);
70
71 SimInfo::MoleculeIterator mi;
72 Molecule::AtomIterator ai;
73 Molecule::RigidBodyIterator rbIter;
74 Molecule::BondIterator bondIter;
75 Molecule::BendIterator bendIter;
76 Molecule::TorsionIterator torsionIter;
77 Molecule::InversionIterator inversionIter;
78
79 Molecule* mol;
80 Atom* atom;
81 RigidBody* rb;
82 Bond* bond;
83 Bend* bend;
84 Torsion* torsion;
85 Inversion* inversion;
86
87 for (mol = info_->beginMolecule(mi); mol != NULL;
88 mol = info_->nextMolecule(mi)) {
89 molecules_[mol->getGlobalIndex()] = mol;
90
91 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
92 stuntdoubles_[atom->getGlobalIndex()] = atom;
93 }
94 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
95 rb = mol->nextRigidBody(rbIter)) {
96 stuntdoubles_[rb->getGlobalIndex()] = rb;
97 }
98 for (bond = mol->beginBond(bondIter); bond != NULL;
99 bond = mol->nextBond(bondIter)) {
100 bonds_[bond->getGlobalIndex()] = bond;
101 }
102 for (bend = mol->beginBend(bendIter); bend != NULL;
103 bend = mol->nextBend(bendIter)) {
104 bends_[bend->getGlobalIndex()] = bend;
105 }
106 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
107 torsion = mol->nextTorsion(torsionIter)) {
108 torsions_[torsion->getGlobalIndex()] = torsion;
109 }
110 for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
111 inversion = mol->nextInversion(inversionIter)) {
112 inversions_[inversion->getGlobalIndex()] = inversion;
113 }
114 }
115 }
116
117 SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance) {
118 StuntDouble* center;
119 Vector3d centerPos;
120 Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
121 SelectionSet bsResult(nObjects_);
122 assert(bsResult.size() == bs.size());
123
124#ifdef IS_MPI
125 int mol;
126 int proc;
127 RealType data[3];
128 int worldRank;
129 MPI_Comm_rank(MPI_COMM_WORLD, &worldRank);
130#endif
131
132 for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
133 if (stuntdoubles_[j] != NULL) {
134 if (stuntdoubles_[j]->isRigidBody()) {
135 RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
136 rb->updateAtoms();
137 }
138 }
139 }
140
141 SelectionSet bsTemp = bs;
142 bsTemp = bsTemp.parallelReduce();
143
144 for (size_t i = 0; i < bsTemp.bitsets_[STUNTDOUBLE].size(); ++i) {
145 if (bsTemp.bitsets_[STUNTDOUBLE][i]) {
146#ifdef IS_MPI
147
148 // Now, if we own stuntdouble i, we can use the position, but in
149 // parallel, we'll need to let everyone else know what that
150 // position is!
151
152 mol = info_->getGlobalMolMembership(i);
153 proc = info_->getMolToProc(mol);
154
155 if (proc == worldRank) {
156 center = stuntdoubles_[i];
157 centerPos = center->getPos();
158 data[0] = centerPos.x();
159 data[1] = centerPos.y();
160 data[2] = centerPos.z();
161 MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
162 } else {
163 MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
164 centerPos = Vector3d(data);
165 }
166#else
167 center = stuntdoubles_[i];
168 centerPos = center->getPos();
169#endif
170 for (size_t j = 0; j < molecules_.size(); ++j) {
171 if (molecules_[j] != NULL) {
172 Vector3d r = centerPos - molecules_[j]->getCom();
173 currSnapshot->wrapVector(r);
174 if (r.length() <= distance) {
175 bsResult.bitsets_[MOLECULE].setBitOn(j);
176 }
177 }
178 }
179 for (size_t j = 0; j < stuntdoubles_.size(); ++j) {
180 if (stuntdoubles_[j] != NULL) {
181 Vector3d r = centerPos - stuntdoubles_[j]->getPos();
182 currSnapshot->wrapVector(r);
183 if (r.length() <= distance) {
184 bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
185 }
186 }
187 }
188 for (size_t j = 0; j < bonds_.size(); ++j) {
189 if (bonds_[j] != NULL) {
190 Vector3d loc = bonds_[j]->getAtomA()->getPos();
191 loc += bonds_[j]->getAtomB()->getPos();
192 loc = loc / 2.0;
193 Vector3d r = centerPos - loc;
194 currSnapshot->wrapVector(r);
195 if (r.length() <= distance) { bsResult.bitsets_[BOND].setBitOn(j); }
196 }
197 }
198 for (size_t j = 0; j < bends_.size(); ++j) {
199 if (bends_[j] != NULL) {
200 Vector3d loc = bends_[j]->getAtomA()->getPos();
201 loc += bends_[j]->getAtomB()->getPos();
202 loc += bends_[j]->getAtomC()->getPos();
203 loc = loc / 3.0;
204 Vector3d r = centerPos - loc;
205 currSnapshot->wrapVector(r);
206 if (r.length() <= distance) { bsResult.bitsets_[BEND].setBitOn(j); }
207 }
208 }
209 for (size_t j = 0; j < torsions_.size(); ++j) {
210 if (torsions_[j] != NULL) {
211 Vector3d loc = torsions_[j]->getAtomA()->getPos();
212 loc += torsions_[j]->getAtomB()->getPos();
213 loc += torsions_[j]->getAtomC()->getPos();
214 loc += torsions_[j]->getAtomD()->getPos();
215 loc = loc / 4.0;
216 Vector3d r = centerPos - loc;
217 currSnapshot->wrapVector(r);
218 if (r.length() <= distance) {
219 bsResult.bitsets_[TORSION].setBitOn(j);
220 }
221 }
222 }
223 for (size_t j = 0; j < inversions_.size(); ++j) {
224 if (inversions_[j] != NULL) {
225 Vector3d loc = inversions_[j]->getAtomA()->getPos();
226 loc += inversions_[j]->getAtomB()->getPos();
227 loc += inversions_[j]->getAtomC()->getPos();
228 loc += inversions_[j]->getAtomD()->getPos();
229 loc = loc / 4.0;
230 Vector3d r = centerPos - loc;
231 currSnapshot->wrapVector(r);
232 if (r.length() <= distance) {
233 bsResult.bitsets_[INVERSION].setBitOn(j);
234 }
235 }
236 }
237 }
238 }
239 return bsResult;
240 }
241
242 SelectionSet DistanceFinder::find(const SelectionSet& bs, RealType distance,
243 int frame) {
244 StuntDouble* center;
245 Vector3d centerPos;
246 Snapshot* currSnapshot = info_->getSnapshotManager()->getSnapshot(frame);
247 SelectionSet bsResult(nObjects_);
248 assert(bsResult.size() == bs.size());
249
250#ifdef IS_MPI
251 int mol;
252 int proc;
253 RealType data[3];
254 int worldRank;
255 MPI_Comm_rank(MPI_COMM_WORLD, &worldRank);
256#endif
257
258 for (unsigned int j = 0; j < stuntdoubles_.size(); ++j) {
259 if (stuntdoubles_[j] != NULL) {
260 if (stuntdoubles_[j]->isRigidBody()) {
261 RigidBody* rb = static_cast<RigidBody*>(stuntdoubles_[j]);
262 rb->updateAtoms(frame);
263 }
264 }
265 }
266
267 SelectionSet bsTemp = bs;
268 bsTemp = bsTemp.parallelReduce();
269
270 for (size_t i = 0; i < bsTemp.bitsets_[STUNTDOUBLE].size(); ++i) {
271 if (bsTemp.bitsets_[STUNTDOUBLE][i]) {
272 // Now, if we own stuntdouble i, we can use the position, but in
273 // parallel, we'll need to let everyone else know what that
274 // position is!
275
276#ifdef IS_MPI
277 mol = info_->getGlobalMolMembership(i);
278 proc = info_->getMolToProc(mol);
279
280 if (proc == worldRank) {
281 center = stuntdoubles_[i];
282 centerPos = center->getPos(frame);
283 data[0] = centerPos.x();
284 data[1] = centerPos.y();
285 data[2] = centerPos.z();
286 MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
287 } else {
288 MPI_Bcast(data, 3, MPI_REALTYPE, proc, MPI_COMM_WORLD);
289 centerPos = Vector3d(data);
290 }
291#else
292 center = stuntdoubles_[i];
293 centerPos = center->getPos(frame);
294#endif
295 for (size_t j = 0; j < molecules_.size(); ++j) {
296 if (molecules_[j] != NULL) {
297 Vector3d r = centerPos - molecules_[j]->getCom(frame);
298 currSnapshot->wrapVector(r);
299 if (r.length() <= distance) {
300 bsResult.bitsets_[MOLECULE].setBitOn(j);
301 }
302 }
303 }
304
305 for (size_t j = 0; j < stuntdoubles_.size(); ++j) {
306 if (stuntdoubles_[j] != NULL) {
307 Vector3d r = centerPos - stuntdoubles_[j]->getPos(frame);
308 currSnapshot->wrapVector(r);
309 if (r.length() <= distance) {
310 bsResult.bitsets_[STUNTDOUBLE].setBitOn(j);
311 }
312 }
313 }
314 for (size_t j = 0; j < bonds_.size(); ++j) {
315 if (bonds_[j] != NULL) {
316 Vector3d loc = bonds_[j]->getAtomA()->getPos(frame);
317 loc += bonds_[j]->getAtomB()->getPos(frame);
318 loc = loc / 2.0;
319 Vector3d r = centerPos - loc;
320 currSnapshot->wrapVector(r);
321 if (r.length() <= distance) { bsResult.bitsets_[BOND].setBitOn(j); }
322 }
323 }
324 for (size_t j = 0; j < bends_.size(); ++j) {
325 if (bends_[j] != NULL) {
326 Vector3d loc = bends_[j]->getAtomA()->getPos(frame);
327 loc += bends_[j]->getAtomB()->getPos(frame);
328 loc += bends_[j]->getAtomC()->getPos(frame);
329 loc = loc / 3.0;
330 Vector3d r = centerPos - loc;
331 currSnapshot->wrapVector(r);
332 if (r.length() <= distance) { bsResult.bitsets_[BEND].setBitOn(j); }
333 }
334 }
335 for (size_t j = 0; j < torsions_.size(); ++j) {
336 if (torsions_[j] != NULL) {
337 Vector3d loc = torsions_[j]->getAtomA()->getPos(frame);
338 loc += torsions_[j]->getAtomB()->getPos(frame);
339 loc += torsions_[j]->getAtomC()->getPos(frame);
340 loc += torsions_[j]->getAtomD()->getPos(frame);
341 loc = loc / 4.0;
342 Vector3d r = centerPos - loc;
343 currSnapshot->wrapVector(r);
344 if (r.length() <= distance) {
345 bsResult.bitsets_[TORSION].setBitOn(j);
346 }
347 }
348 }
349 for (size_t j = 0; j < inversions_.size(); ++j) {
350 if (inversions_[j] != NULL) {
351 Vector3d loc = inversions_[j]->getAtomA()->getPos(frame);
352 loc += inversions_[j]->getAtomB()->getPos(frame);
353 loc += inversions_[j]->getAtomC()->getPos(frame);
354 loc += inversions_[j]->getAtomD()->getPos(frame);
355 loc = loc / 4.0;
356 Vector3d r = centerPos - loc;
357 currSnapshot->wrapVector(r);
358 if (r.length() <= distance) {
359 bsResult.bitsets_[INVERSION].setBitOn(j);
360 }
361 }
362 }
363 }
364 }
365 return bsResult;
366 }
367} // namespace OpenMD
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:60
OpenMD::STUNTDOUBLE
@ STUNTDOUBLE
StuntDoubles (Atoms & RigidBodies)
Definition SelectionSet.hpp:61