OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
FragmentStamp.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#include "types/FragmentStamp.hpp"
45
46#include <algorithm>
47#include <functional>
48#include <iostream>
49#include <sstream>
50#include <tuple>
51
52#include "utils/MemoryUtils.hpp"
53
54using namespace std;
55
56namespace OpenMD {
57
58 template<class ContainerType>
59 bool hasDuplicateElement(const ContainerType& cont) {
60 ContainerType tmp = cont;
61 std::sort(tmp.begin(), tmp.end());
62 tmp.erase(std::unique(tmp.begin(), tmp.end()), tmp.end());
63 return tmp.size() != cont.size();
64 }
65
66 FragmentStamp::FragmentStamp() { DefineParameter(Name, "name"); }
67
68 FragmentStamp::~FragmentStamp() {
69 Utils::deletePointers(atomStamps_);
70 Utils::deletePointers(bondStamps_);
71 Utils::deletePointers(bendStamps_);
72 Utils::deletePointers(torsionStamps_);
73 Utils::deletePointers(inversionStamps_);
74 Utils::deletePointers(rigidBodyStamps_);
75 Utils::deletePointers(cutoffGroupStamps_);
76 Utils::deletePointers(nodesStamps_);
77 Utils::deletePointers(constraintStamps_);
78 }
79
80 bool FragmentStamp::addAtomStamp(AtomStamp* atom) {
81 bool ret = addIndexSensitiveStamp(atomStamps_, atom);
82 if (!ret) {
83 std::ostringstream oss;
84 oss << "Error in Fragment " << getName()
85 << ": multiple atoms have the same indices" << atom->getIndex()
86 << "\n";
87 throw OpenMDException(oss.str());
88 }
89 return ret;
90 }
91
92 bool FragmentStamp::addBondStamp(BondStamp* bond) {
93 bondStamps_.push_back(bond);
94 return true;
95 }
96
97 bool FragmentStamp::addBendStamp(BendStamp* bend) {
98 bendStamps_.push_back(bend);
99 return true;
100 }
101
102 bool FragmentStamp::addTorsionStamp(TorsionStamp* torsion) {
103 torsionStamps_.push_back(torsion);
104 return true;
105 }
106 bool FragmentStamp::addInversionStamp(InversionStamp* inversion) {
107 inversionStamps_.push_back(inversion);
108 return true;
109 }
110
111 bool FragmentStamp::addRigidBodyStamp(RigidBodyStamp* rigidbody) {
112 bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
113 if (!ret) {
114 std::ostringstream oss;
115 oss << "Error in Fragment " << getName()
116 << ": multiple rigidbodies have the same indices: "
117 << rigidbody->getIndex() << "\n";
118 throw OpenMDException(oss.str());
119 }
120 return ret;
121 }
122
123 bool FragmentStamp::addCutoffGroupStamp(CutoffGroupStamp* cutoffgroup) {
124 cutoffGroupStamps_.push_back(cutoffgroup);
125 return true;
126 }
127
128 bool FragmentStamp::addNodesStamp(NodesStamp* nodes) {
129 nodesStamps_.push_back(nodes);
130 return true;
131 }
132
133 bool FragmentStamp::addConstraintStamp(ConstraintStamp* constraint) {
134 constraintStamps_.push_back(constraint);
135 return true;
136 }
137
138 void FragmentStamp::validate() {
139 DataHolder::validate();
140 CheckParameter(Name, isNotEmpty());
141
142 atom2Rigidbody.resize(getNAtoms());
143
144 // A negative number means the atom is a free atom, and does not
145 // belong to rigidbody. Every element in atom2Rigidbody has unique
146 // negative number at the very beginning
147
148 for (unsigned int i = 0; i < atom2Rigidbody.size(); ++i) {
149 atom2Rigidbody[i] = -1 - int(i);
150 }
151 for (std::size_t i = 0; i < getNRigidBodies(); ++i) {
152 RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
153 std::vector<int> members = rbStamp->getMembers();
154 for (std::vector<int>::iterator j = members.begin(); j != members.end();
155 ++j) {
156 atom2Rigidbody[*j] = i;
157 }
158 }
159 checkAtoms();
160 checkBonds();
161 fillBondInfo();
162 checkBends();
163 checkTorsions();
164 checkInversions();
165 checkRigidBodies();
166 checkCutoffGroups();
167 checkConstraints();
168 checkNodes();
169 }
170
171 void FragmentStamp::checkAtoms() {
172 std::vector<AtomStamp*>::iterator ai = std::find(
173 atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
174 if (ai != atomStamps_.end()) {
175 std::ostringstream oss;
176 oss << "Error in Fragment " << getName() << ": atom["
177 << ai - atomStamps_.begin() << "] is missing\n";
178 throw OpenMDException(oss.str());
179 }
180 }
181
182 void FragmentStamp::checkBonds() {
183 std::ostringstream oss;
184 // make sure index is not out of range
185 int natoms = getNAtoms();
186 for (std::size_t i = 0; i < getNBonds(); ++i) {
187 BondStamp* bondStamp = getBondStamp(i);
188 if (bondStamp->getA() > natoms - 1 || bondStamp->getA() < 0 ||
189 bondStamp->getB() > natoms - 1 || bondStamp->getB() < 0 ||
190 bondStamp->getA() == bondStamp->getB()) {
191 oss << "Error in Fragment " << getName() << ": bond("
192 << bondStamp->getA() << ", " << bondStamp->getB()
193 << ") is invalid\n";
194 throw OpenMDException(oss.str());
195 }
196 }
197
198 // make sure bonds are unique
199 std::set<std::pair<int, int>> allBonds;
200 for (std::size_t i = 0; i < getNBonds(); ++i) {
201 BondStamp* bondStamp = getBondStamp(i);
202 std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
203 // make sure bondPair.first is always less than or equal to
204 // bondPair.third
205 if (bondPair.first > bondPair.second) {
206 std::swap(bondPair.first, bondPair.second);
207 }
208
209 std::set<std::pair<int, int>>::iterator iter = allBonds.find(bondPair);
210 if (iter != allBonds.end()) {
211 oss << "Error in Fragment " << getName() << ": "
212 << "bond(" << iter->first << ", " << iter->second
213 << ") appears multiple times\n";
214 throw OpenMDException(oss.str());
215 } else {
216 allBonds.insert(bondPair);
217 }
218 }
219
220 // make sure atoms belong to same rigidbody do not bond to each other
221 for (std::size_t i = 0; i < getNBonds(); ++i) {
222 BondStamp* bondStamp = getBondStamp(i);
223 if (atom2Rigidbody[bondStamp->getA()] ==
224 atom2Rigidbody[bondStamp->getB()]) {
225 oss << "Error in Fragment " << getName() << ": "
226 << "bond(" << bondStamp->getA() << ", " << bondStamp->getB()
227 << ") belong to same rigidbody "
228 << atom2Rigidbody[bondStamp->getA()] << "\n";
229 throw OpenMDException(oss.str());
230 }
231 }
232 }
233
234 void FragmentStamp::checkBends() {
235 std::ostringstream oss;
236 for (std::size_t i = 0; i < getNBends(); ++i) {
237 BendStamp* bendStamp = getBendStamp(i);
238 std::vector<int> bendAtoms = bendStamp->getMembers();
239 std::vector<int>::iterator j = std::find_if(
240 bendAtoms.begin(), bendAtoms.end(),
241 std::bind(std::greater<int>(), placeholders::_1, getNAtoms() - 1));
242 std::vector<int>::iterator k =
243 std::find_if(bendAtoms.begin(), bendAtoms.end(),
244 std::bind(std::less<int>(), placeholders::_1, 0));
245
246 if (j != bendAtoms.end() || k != bendAtoms.end()) {
247 oss << "Error in Fragment " << getName() << " : atoms of bend"
248 << containerToString(bendAtoms) << " have invalid indices\n";
249 throw OpenMDException(oss.str());
250 }
251
252 if (hasDuplicateElement(bendAtoms)) {
253 oss << "Error in Fragment " << getName() << " : atoms of bend"
254 << containerToString(bendAtoms) << " have duplicated indices\n";
255 throw OpenMDException(oss.str());
256 }
257
258 if (bendAtoms.size() == 2) {
259 if (!bendStamp->haveGhostVectorSource()) {
260 oss << "Error in Fragment " << getName()
261 << ": ghostVectorSouce is missing\n";
262 throw OpenMDException(oss.str());
263 } else {
264 std::size_t ghostIndex = bendStamp->getGhostVectorSource();
265 if (ghostIndex < getNAtoms()) {
266 if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) ==
267 bendAtoms.end()) {
268 oss << "Error in Fragment " << getName() << ": ghostVectorSouce "
269 << ghostIndex << "is invalid\n";
270 throw OpenMDException(oss.str());
271 }
272 if (!getAtomStamp(ghostIndex)->haveOrientation()) {
273 oss << "Error in Fragment " << getName()
274 << ": ghost atom must be a directional atom\n";
275 throw OpenMDException(oss.str());
276 }
277 } else {
278 oss << "Error in Fragment " << getName() << ": ghostVectorSource "
279 << ghostIndex << " is invalid\n";
280 throw OpenMDException(oss.str());
281 }
282 }
283 } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
284 oss << "Error in Fragment " << getName()
285 << ": normal bend should not have ghostVectorSouce\n";
286 throw OpenMDException(oss.str());
287 }
288 }
289
290 for (std::size_t i = 0; i < getNBends(); ++i) {
291 BendStamp* bendStamp = getBendStamp(i);
292 std::vector<int> bendAtoms = bendStamp->getMembers();
293 std::vector<int> rigidSet(getNRigidBodies(), 0);
294 std::vector<int>::iterator j;
295 for (j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
296 int rigidbodyIndex = atom2Rigidbody[*j];
297 if (rigidbodyIndex >= 0) {
298 ++rigidSet[rigidbodyIndex];
299 if (rigidSet[rigidbodyIndex] > 2) {
300 oss << "Error in Fragment " << getName() << ": bend"
301 << containerToString(bendAtoms)
302 << "has three atoms on the same rigid body\n";
303 throw OpenMDException(oss.str());
304 }
305 }
306 }
307 }
308
309 std::set<std::tuple<int, int, int>> allBends;
310 std::set<std::tuple<int, int, int>>::iterator iter;
311 for (std::size_t i = 0; i < getNBends(); ++i) {
312 BendStamp* bendStamp = getBendStamp(i);
313 std::vector<int> bend = bendStamp->getMembers();
314 if (bend.size() == 2) {
315 // in case we have two ghost bend. For example,
316 // bend {
317 // members (0, 1);
318 // ghostVectorSource = 0;
319 // }
320 // and
321 // bend {
322 // members (0, 1);
323 // ghostVectorSource = 0;
324 // }
325 // In order to distinguish them. we expand them to Tuple3.
326 // the first one is expanded to (0, 0, 1) while the second one
327 // is expaned to (0, 1, 1)
328 int ghostIndex = bendStamp->getGhostVectorSource();
329 std::vector<int>::iterator j =
330 std::find(bend.begin(), bend.end(), ghostIndex);
331 if (j != bend.end()) { bend.insert(j, ghostIndex); }
332 }
333
334 std::tuple<int, int, int> bendTuple {bend[0], bend[1], bend[2]};
335 auto& [first, second, third] = bendTuple;
336
337 // make sure bendTuple.first is always less than or equal to
338 // bendTuple.third
339 if (first > third) { std::swap(first, third); }
340
341 iter = allBends.find(bendTuple);
342 if (iter != allBends.end()) {
343 oss << "Error in Fragment " << getName() << ": "
344 << "Bend" << containerToString(bend) << " appears multiple times\n";
345 throw OpenMDException(oss.str());
346 } else {
347 allBends.insert(bendTuple);
348 }
349 }
350 }
351
352 void FragmentStamp::checkTorsions() {
353 std::ostringstream oss;
354 for (std::size_t i = 0; i < getNTorsions(); ++i) {
355 TorsionStamp* torsionStamp = getTorsionStamp(i);
356 std::vector<int> torsionAtoms = torsionStamp->getMembers();
357 std::vector<int>::iterator j = std::find_if(
358 torsionAtoms.begin(), torsionAtoms.end(),
359 std::bind(std::greater<int>(), placeholders::_1, getNAtoms() - 1));
360 std::vector<int>::iterator k =
361 std::find_if(torsionAtoms.begin(), torsionAtoms.end(),
362 std::bind(std::less<int>(), placeholders::_1, 0));
363
364 if (j != torsionAtoms.end() || k != torsionAtoms.end()) {
365 oss << "Error in Fragment " << getName() << ": atoms of torsion"
366 << containerToString(torsionAtoms) << " have invalid indices\n";
367 throw OpenMDException(oss.str());
368 }
369 if (hasDuplicateElement(torsionAtoms)) {
370 oss << "Error in Fragment " << getName() << " : atoms of torsion"
371 << containerToString(torsionAtoms) << " have duplicated indices\n";
372 throw OpenMDException(oss.str());
373 }
374 }
375
376 for (std::size_t i = 0; i < getNTorsions(); ++i) {
377 TorsionStamp* torsionStamp = getTorsionStamp(i);
378 std::vector<int> torsionAtoms = torsionStamp->getMembers();
379 std::vector<int> rigidSet(getNRigidBodies(), 0);
380 std::vector<int>::iterator j;
381 for (j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
382 int rigidbodyIndex = atom2Rigidbody[*j];
383 if (rigidbodyIndex >= 0) {
384 ++rigidSet[rigidbodyIndex];
385 if (rigidSet[rigidbodyIndex] > 3) {
386 oss << "Error in Fragment " << getName() << ": torsion"
387 << containerToString(torsionAtoms)
388 << "has four atoms on the same rigid body\n";
389 throw OpenMDException(oss.str());
390 }
391 }
392 }
393 }
394
395 std::set<std::tuple<int, int, int, int>> allTorsions;
396 std::set<std::tuple<int, int, int, int>>::iterator iter;
397 for (std::size_t i = 0; i < getNTorsions(); ++i) {
398 TorsionStamp* torsionStamp = getTorsionStamp(i);
399 std::vector<int> torsion = torsionStamp->getMembers();
400 if (torsion.size() == 3) {
401 int ghostIndex = torsionStamp->getGhostVectorSource();
402 std::vector<int>::iterator j =
403 std::find(torsion.begin(), torsion.end(), ghostIndex);
404 if (j != torsion.end()) { torsion.insert(j, ghostIndex); }
405 }
406
407 std::tuple<int, int, int, int> torsionTuple(torsion[0], torsion[1],
408 torsion[2], torsion[3]);
409 auto& [first, second, third, fourth] = torsionTuple;
410
411 if (first > fourth) {
412 std::swap(first, fourth);
413 std::swap(second, third);
414 }
415
416 iter = allTorsions.find(torsionTuple);
417 if (iter == allTorsions.end()) {
418 allTorsions.insert(torsionTuple);
419 } else {
420 oss << "Error in Fragment " << getName() << ": "
421 << "Torsion" << containerToString(torsion)
422 << " appears multiple times\n";
423 throw OpenMDException(oss.str());
424 }
425 }
426 }
427
428 void FragmentStamp::checkInversions() {
429 std::ostringstream oss;
430
431 // first we automatically find the other three atoms that
432 // are satellites of an inversion center:
433
434 for (std::size_t i = 0; i < getNInversions(); ++i) {
435 InversionStamp* inversionStamp = getInversionStamp(i);
436 int center = inversionStamp->getCenter();
437 std::vector<int> satellites;
438
439 // Some inversions come pre-programmed with the satellites. If
440 // so, don't add the satellites as they are already there.
441
442 if (inversionStamp->getNSatellites() != 3) {
443 for (std::size_t j = 0; j < getNBonds(); ++j) {
444 BondStamp* bondStamp = getBondStamp(j);
445 int a = bondStamp->getA();
446 int b = bondStamp->getB();
447
448 if (a == center) { satellites.push_back(b); }
449 if (b == center) { satellites.push_back(a); }
450 }
451
452 if (satellites.size() == 3) {
453 std::sort(satellites.begin(), satellites.end());
454 inversionStamp->setSatellites(satellites);
455 } else {
456 oss << "Error in Fragment " << getName() << ": found wrong number"
457 << " of bonds for inversion center " << center;
458 throw OpenMDException(oss.str());
459 }
460 }
461 }
462
463 // then we do some sanity checking on the inversions:
464
465 for (std::size_t i = 0; i < getNInversions(); ++i) {
466 InversionStamp* inversionStamp = getInversionStamp(i);
467
468 std::vector<int> inversionAtoms = inversionStamp->getSatellites();
469 // add the central atom to the beginning of the list:
470 inversionAtoms.insert(inversionAtoms.begin(),
471 inversionStamp->getCenter());
472
473 std::vector<int>::iterator j = std::find_if(
474 inversionAtoms.begin(), inversionAtoms.end(),
475 std::bind(std::greater<int>(), placeholders::_1, getNAtoms() - 1));
476 std::vector<int>::iterator k =
477 std::find_if(inversionAtoms.begin(), inversionAtoms.end(),
478 std::bind(std::less<int>(), placeholders::_1, 0));
479
480 if (j != inversionAtoms.end() || k != inversionAtoms.end()) {
481 oss << "Error in Fragment " << getName() << ": atoms of inversion"
482 << containerToString(inversionAtoms) << " have invalid indices\n";
483 throw OpenMDException(oss.str());
484 }
485
486 if (hasDuplicateElement(inversionAtoms)) {
487 oss << "Error in Fragment " << getName() << " : atoms of inversion"
488 << containerToString(inversionAtoms)
489 << " have duplicated indices\n";
490 throw OpenMDException(oss.str());
491 }
492 }
493
494 for (std::size_t i = 0; i < getNInversions(); ++i) {
495 InversionStamp* inversionStamp = getInversionStamp(i);
496 std::vector<int> inversionAtoms = inversionStamp->getSatellites();
497 inversionAtoms.push_back(inversionStamp->getCenter());
498 std::vector<int> rigidSet(getNRigidBodies(), 0);
499 std::vector<int>::iterator j;
500 for (j = inversionAtoms.begin(); j != inversionAtoms.end(); ++j) {
501 int rigidbodyIndex = atom2Rigidbody[*j];
502 if (rigidbodyIndex >= 0) {
503 ++rigidSet[rigidbodyIndex];
504 if (rigidSet[rigidbodyIndex] > 3) {
505 oss << "Error in Fragment " << getName()
506 << ": inversion centered on atom "
507 << inversionStamp->getCenter()
508 << " has four atoms that belong to same rigidbody "
509 << rigidbodyIndex << "\n";
510 throw OpenMDException(oss.str());
511 }
512 }
513 }
514 }
515
516 std::set<std::tuple<int, int, int, int>> allInversions;
517 std::set<std::tuple<int, int, int, int>>::iterator iter;
518 for (std::size_t i = 0; i < getNInversions(); ++i) {
519 InversionStamp* inversionStamp = getInversionStamp(i);
520 int cent = inversionStamp->getCenter();
521 std::vector<int> inversion = inversionStamp->getSatellites();
522
523 std::tuple<int, int, int, int> inversionTuple(cent, inversion[0],
524 inversion[1], inversion[2]);
525 auto& [first, second, third, fourth] = inversionTuple;
526
527 // In OpenMD, the Central atom in an inversion comes first, and
528 // has a special position. The other three atoms can come in
529 // random order, and should be sorted in increasing numerical
530 // order to check for duplicates. This requires three pairwise
531 // swaps:
532 if (third > fourth) std::swap(third, fourth);
533 if (second > third) std::swap(second, third);
534 if (third > fourth) std::swap(third, fourth);
535
536 iter = allInversions.find(inversionTuple);
537 if (iter == allInversions.end()) {
538 allInversions.insert(inversionTuple);
539 } else {
540 oss << "Error in Fragment " << getName() << ": "
541 << "Inversion" << containerToString(inversion)
542 << " appears multiple times\n";
543 throw OpenMDException(oss.str());
544 }
545 }
546 }
547
548 void FragmentStamp::checkRigidBodies() {
549 std::ostringstream oss;
550 std::vector<RigidBodyStamp*>::iterator ri =
551 std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(),
552 static_cast<RigidBodyStamp*>(NULL));
553 if (ri != rigidBodyStamps_.end()) {
554 oss << "Error in Fragment " << getName() << ":rigidBody["
555 << ri - rigidBodyStamps_.begin() << "] is missing\n";
556 throw OpenMDException(oss.str());
557 }
558
559 for (std::size_t i = 0; i < getNRigidBodies(); ++i) {
560 RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
561 std::vector<int> rigidAtoms = rbStamp->getMembers();
562 std::vector<int>::iterator j = std::find_if(
563 rigidAtoms.begin(), rigidAtoms.end(),
564 std::bind(std::greater<int>(), placeholders::_1, getNAtoms() - 1));
565 if (j != rigidAtoms.end()) {
566 oss << "Error in Fragment " << getName();
567 throw OpenMDException(oss.str());
568 }
569 }
570 }
571
572 void FragmentStamp::checkCutoffGroups() {
573 std::vector<AtomStamp*>::iterator ai;
574 std::vector<int>::iterator fai;
575
576 // add all atoms into freeAtoms_ set
577 for (ai = atomStamps_.begin(); ai != atomStamps_.end(); ++ai) {
578 freeAtoms_.push_back((*ai)->getIndex());
579 }
580
581 for (std::size_t i = 0; i < getNCutoffGroups(); ++i) {
582 CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
583 std::vector<int> cutoffGroupAtoms = cutoffGroupStamp->getMembers();
584 std::vector<int>::iterator j = std::find_if(
585 cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(),
586 std::bind(std::greater<int>(), placeholders::_1, getNAtoms() - 1));
587 if (j != cutoffGroupAtoms.end()) {
588 std::ostringstream oss;
589 oss << "Error in Fragment " << getName() << ": cutoffGroup"
590 << " is out of range\n";
591 throw OpenMDException(oss.str());
592 }
593
594 for (fai = cutoffGroupAtoms.begin(); fai != cutoffGroupAtoms.end();
595 ++fai) {
596 // erase the atoms belonging to cutoff groups from freeAtoms_ vector
597 freeAtoms_.erase(
598 std::remove(freeAtoms_.begin(), freeAtoms_.end(), (*fai)),
599 freeAtoms_.end());
600 }
601 }
602 }
603
604 void FragmentStamp::checkConstraints() {
605 std::ostringstream oss;
606 // make sure index is not out of range
607 int natoms = getNAtoms();
608 for (std::size_t i = 0; i < getNConstraints(); ++i) {
609 ConstraintStamp* constraintStamp = getConstraintStamp(i);
610 if (constraintStamp->getA() > natoms - 1 || constraintStamp->getA() < 0 ||
611 constraintStamp->getB() > natoms - 1 || constraintStamp->getB() < 0 ||
612 constraintStamp->getA() == constraintStamp->getB()) {
613 oss << "Error in Fragment " << getName() << ": constraint("
614 << constraintStamp->getA() << ", " << constraintStamp->getB()
615 << ") is invalid\n";
616 throw OpenMDException(oss.str());
617 }
618 }
619
620 // make sure constraints are unique
621 std::set<std::pair<int, int>> allConstraints;
622 for (std::size_t i = 0; i < getNConstraints(); ++i) {
623 ConstraintStamp* constraintStamp = getConstraintStamp(i);
624 std::pair<int, int> constraintPair(constraintStamp->getA(),
625 constraintStamp->getB());
626 // make sure constraintPair.first is always less than or equal to
627 // constraintPair.third
628 if (constraintPair.first > constraintPair.second) {
629 std::swap(constraintPair.first, constraintPair.second);
630 }
631
632 std::set<std::pair<int, int>>::iterator iter =
633 allConstraints.find(constraintPair);
634 if (iter != allConstraints.end()) {
635 oss << "Error in Fragment " << getName() << ": "
636 << "constraint(" << iter->first << ", " << iter->second
637 << ") appears multiple times\n";
638 throw OpenMDException(oss.str());
639 } else {
640 allConstraints.insert(constraintPair);
641 }
642 }
643
644 // make sure atoms belong to same rigidbody are not constrained to
645 // each other
646 for (std::size_t i = 0; i < getNConstraints(); ++i) {
647 ConstraintStamp* constraintStamp = getConstraintStamp(i);
648 if (atom2Rigidbody[constraintStamp->getA()] ==
649 atom2Rigidbody[constraintStamp->getB()]) {
650 oss << "Error in Fragment " << getName() << ": "
651 << "constraint(" << constraintStamp->getA() << ", "
652 << constraintStamp->getB() << ") belong to same rigidbody "
653 << atom2Rigidbody[constraintStamp->getA()] << "\n";
654 throw OpenMDException(oss.str());
655 }
656 }
657 }
658
659 void FragmentStamp::checkNodes() {
660 std::ostringstream oss;
661 std::vector<NodesStamp*>::iterator ni =
662 std::find(nodesStamps_.begin(), nodesStamps_.end(),
663 static_cast<NodesStamp*>(NULL));
664 if (ni != nodesStamps_.end()) {
665 oss << "Error in Molecule " << getName() << ":nodes["
666 << ni - nodesStamps_.begin() << "] is missing\n";
667 throw OpenMDException(oss.str());
668 }
669
670 for (std::size_t i = 0; i < getNNodes(); ++i) {
671 NodesStamp* nStamp = getNodesStamp(i);
672 std::vector<int> nodeAtoms = nStamp->getMembers();
673 std::vector<int>::iterator j = std::find_if(
674 nodeAtoms.begin(), nodeAtoms.end(),
675 std::bind(std::greater<int>(), placeholders::_1, getNAtoms() - 1));
676 if (j != nodeAtoms.end()) {
677 oss << "Error in Fragment " << getName();
678 throw OpenMDException(oss.str());
679 }
680 }
681 }
682
683 void FragmentStamp::fillBondInfo() {
684 for (std::size_t i = 0; i < getNBonds(); ++i) {
685 BondStamp* bondStamp = getBondStamp(i);
686 int a = bondStamp->getA();
687 int b = bondStamp->getB();
688 AtomStamp* atomA = getAtomStamp(a);
689 AtomStamp* atomB = getAtomStamp(b);
690 atomA->addBond(i);
691 atomA->addBondedAtom(b);
692 atomB->addBond(i);
693 atomB->addBondedAtom(a);
694 }
695 }
696
697 // Function Name: isBondInSameRigidBody
698 // Returns true is both atoms of the bond belong to the same rigid
699 // body, otherwise return false
700 bool FragmentStamp::isBondInSameRigidBody(BondStamp* bond) {
701 int rbA;
702 int rbB;
703 int consAtomA;
704 int consAtomB;
705
706 if (!isAtomInRigidBody(bond->getA(), rbA, consAtomA)) return false;
707
708 if (!isAtomInRigidBody(bond->getB(), rbB, consAtomB)) return false;
709
710 if (rbB == rbA)
711 return true;
712 else
713 return false;
714 }
715
716 // Function Name: isAtomInRigidBody
717 // Returns false if atom does not belong to a rigid body, otherwise
718 // returns true
719 bool FragmentStamp::isAtomInRigidBody(int atomIndex) {
720 return atom2Rigidbody[atomIndex] >= 0;
721 }
722
723 // Function Name: isAtomInRigidBody
724 // Returns false if atom does not belong to a rigid body otherwise
725 // returns true and sets whichRigidBody and consAtomIndex
726 // atomIndex : the index of atom in component
727 // whichRigidBody: the index of the rigidbody in the component
728 // consAtomIndex: the position the joint atom appears in the rigidbody's
729 // definition
730 bool FragmentStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody,
731 int& consAtomIndex) {
732 whichRigidBody = -1;
733 consAtomIndex = -1;
734
735 if (atom2Rigidbody[atomIndex] >= 0) {
736 whichRigidBody = atom2Rigidbody[atomIndex];
737 RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
738 int numAtom = rbStamp->getNMembers();
739 for (int j = 0; j < numAtom; j++) {
740 if (rbStamp->getMemberAt(j) == atomIndex) {
741 consAtomIndex = j;
742 return true;
743 }
744 }
745 }
746
747 return false;
748 }
749
750 // Returns the position of joint atoms apearing in a rigidbody's definition
751 // For the time being, we will use the most inefficient algorithm,
752 // the complexity is O(N^2). We could improve the
753 // complexity to O(NlogN) by sorting the atom index in rigid body
754 // first
755 std::vector<std::pair<int, int>> FragmentStamp::getJointAtoms(int rb1,
756 int rb2) {
757 RigidBodyStamp* rbStamp1;
758 RigidBodyStamp* rbStamp2;
759 int natomInRb1;
760 int natomInRb2;
761 int atomIndex1;
762 int atomIndex2;
763 std::vector<std::pair<int, int>> jointAtomIndexPair;
764
765 rbStamp1 = this->getRigidBodyStamp(rb1);
766 natomInRb1 = rbStamp1->getNMembers();
767
768 rbStamp2 = this->getRigidBodyStamp(rb2);
769 natomInRb2 = rbStamp2->getNMembers();
770
771 for (int i = 0; i < natomInRb1; i++) {
772 atomIndex1 = rbStamp1->getMemberAt(i);
773
774 for (int j = 0; j < natomInRb2; j++) {
775 atomIndex2 = rbStamp2->getMemberAt(j);
776
777 if (atomIndex1 == atomIndex2) {
778 jointAtomIndexPair.push_back(std::make_pair(i, j));
779 break;
780 }
781 }
782 }
783
784 return jointAtomIndexPair;
785 }
786
787} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.