OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
SelectionEvaluator.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 following paper when you publish your work:
33 *
34 * [1] Drisko et al., J. Open Source Softw. 9, 7004 (2024).
35 *
36 * Good starting points for code and simulation methodology are:
37 *
38 * [2] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
39 * [3] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
40 * [4] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
41 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
42 * [6] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
43 * [7] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
44 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
45 * [9] Drisko & Gezelter, J. Chem. Theory Comput. 20, 4986-4997 (2024).
46 */
47
48#include "selection/SelectionEvaluator.hpp"
49
50#include <any>
51#include <map>
52#include <stack>
53#include <string>
54#include <utility>
55
56#include "io/ifstrstream.hpp"
57#include "primitives/Atom.hpp"
58#include "primitives/DirectionalAtom.hpp"
59#include "primitives/Molecule.hpp"
60#include "primitives/RigidBody.hpp"
61#include "types/FixedChargeAdapter.hpp"
62#include "types/FluctuatingChargeAdapter.hpp"
63
64namespace OpenMD {
65
66 SelectionEvaluator::SelectionEvaluator(SimInfo* si) :
67 info(si), nameFinder(info), distanceFinder(info), hullFinder(info),
68 alphaHullFinder(info), indexFinder(info), isLoaded_(false),
69 hasSurfaceArea_(false), hasVolume_(false) {
70 nObjects.push_back(info->getNGlobalAtoms() + info->getNGlobalRigidBodies());
71 nObjects.push_back(info->getNGlobalBonds());
72 nObjects.push_back(info->getNGlobalBends());
73 nObjects.push_back(info->getNGlobalTorsions());
74 nObjects.push_back(info->getNGlobalInversions());
75 nObjects.push_back(info->getNGlobalMolecules());
76 }
77
78 bool SelectionEvaluator::loadScript(const std::string& filename,
79 const std::string& script) {
80 clearDefinitionsAndLoadPredefined();
81 this->filename = filename;
82 this->script = script;
83 if (!compiler.compile(filename, script)) {
84 error = true;
85 errorMessage = compiler.getErrorMessage();
86
87 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
88 "SelectionCompiler Error: %s\n", errorMessage.c_str());
89 painCave.severity = OPENMD_ERROR;
90 painCave.isFatal = 1;
91 simError();
92 return false;
93 }
94
95 pc = 0;
96 aatoken = compiler.getAatokenCompiled();
97 linenumbers = compiler.getLineNumbers();
98 lineIndices = compiler.getLineIndices();
99
100 std::vector<std::vector<Token>>::const_iterator i;
101
102 isDynamic_ = false;
103 for (i = aatoken.begin(); i != aatoken.end(); ++i) {
104 if (containDynamicToken(*i)) {
105 isDynamic_ = true;
106 break;
107 }
108 }
109
110 isLoaded_ = true;
111 return true;
112 }
113
114 void SelectionEvaluator::clearState() {
115 error = false;
116 errorMessage = "";
117 }
118
119 bool SelectionEvaluator::loadScriptString(const std::string& script) {
120 clearState();
121 return loadScript("", script);
122 }
123
124 bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
125 clearState();
126 return loadScriptFileInternal(filename);
127 }
128
129 bool SelectionEvaluator::loadScriptFileInternal(const std::string& filename) {
130 ifstrstream ifs(filename.c_str());
131 if (!ifs.is_open()) { return false; }
132
133 const int bufferSize = 65535;
134 char buffer[bufferSize];
135 std::string script;
136 while (ifs.getline(buffer, bufferSize)) {
137 script += buffer;
138 }
139 return loadScript(filename, script);
140 }
141
142 void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs) {
143 while (pc < aatoken.size()) {
144 statement = aatoken[pc++];
145 statementLength = statement.size();
146 Token token = statement[0];
147 switch (token.tok) {
148 case Token::define:
149 define();
150 break;
151 case Token::select:
152 select(bs);
153 break;
154 default:
155 unrecognizedCommand(token);
156 return;
157 }
158 }
159 }
160
161 void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs,
162 int frame) {
163 while (pc < aatoken.size()) {
164 statement = aatoken[pc++];
165 statementLength = statement.size();
166 Token token = statement[0];
167 switch (token.tok) {
168 case Token::define:
169 define();
170 break;
171 case Token::select:
172 select(bs, frame);
173 break;
174 default:
175 unrecognizedCommand(token);
176 return;
177 }
178 }
179 }
180
181 SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code,
182 int pcStart) {
183 SelectionSet bs = createSelectionSets();
184 std::stack<SelectionSet> stack;
185 vector<int> bsSize = bs.size();
186
187 for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
188 Token instruction = code[pc];
189
190 switch (instruction.tok) {
191 case Token::expressionBegin:
192 break;
193 case Token::expressionEnd:
194 break;
195 case Token::all:
196 bs = allInstruction();
197 stack.push(bs);
198 break;
199 case Token::none:
200 bs = createSelectionSets();
201 stack.push(bs);
202 break;
203 case Token::opOr:
204 bs = stack.top();
205 stack.pop();
206 stack.top() |= bs;
207 break;
208 case Token::opAnd:
209 bs = stack.top();
210 stack.pop();
211 stack.top() &= bs;
212 break;
213 case Token::opNot:
214 stack.top().flip();
215 break;
216 case Token::within:
217 withinInstruction(instruction, stack.top());
218 break;
219 case Token::alphahull:
220 stack.push(alphaHullInstruction(instruction));
221 break;
222 case Token::hull:
223 stack.push(hull());
224 break;
225 // case Token::selected:
226 // stack.push(getSelectionSet());
227 // break;
228 case Token::name:
229 stack.push(
230 nameInstruction(std::any_cast<std::string>(instruction.value)));
231 break;
232 case Token::index:
233 stack.push(indexInstruction(instruction.value));
234 break;
235 case Token::identifier:
236 stack.push(lookupValue(std::any_cast<std::string>(instruction.value)));
237 break;
238 case Token::opLT:
239 case Token::opLE:
240 case Token::opGE:
241 case Token::opGT:
242 case Token::opEQ:
243 case Token::opNE:
244 stack.push(comparatorInstruction(instruction));
245 break;
246 default:
247 unrecognizedExpression();
248 }
249 }
250 if (stack.size() != 1)
251 evalError("atom expression compiler error - stack over/underflow");
252
253 return stack.top();
254 }
255
256 SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code,
257 int pcStart, int frame) {
258 SelectionSet bs = createSelectionSets();
259 std::stack<SelectionSet> stack;
260
261 for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
262 Token instruction = code[pc];
263
264 switch (instruction.tok) {
265 case Token::expressionBegin:
266 break;
267 case Token::expressionEnd:
268 break;
269 case Token::all:
270 bs = allInstruction();
271 stack.push(bs);
272 break;
273 case Token::none:
274 bs = SelectionSet(nObjects);
275 stack.push(bs);
276 break;
277 case Token::opOr:
278 bs = stack.top();
279 stack.pop();
280 stack.top() |= bs;
281 break;
282 case Token::opAnd:
283 bs = stack.top();
284 stack.pop();
285 stack.top() &= bs;
286 break;
287 case Token::opNot:
288 stack.top().flip();
289 break;
290 case Token::within:
291 withinInstruction(instruction, stack.top(), frame);
292 break;
293 case Token::alphahull:
294 stack.push(alphaHullInstruction(instruction, frame));
295 break;
296 case Token::hull:
297 stack.push(hull(frame));
298 break;
299 // case Token::selected:
300 // stack.push(getSelectionSet());
301 // break;
302 case Token::name:
303 stack.push(
304 nameInstruction(std::any_cast<std::string>(instruction.value)));
305 break;
306 case Token::index:
307 stack.push(indexInstruction(instruction.value));
308 break;
309 case Token::identifier:
310 stack.push(lookupValue(std::any_cast<std::string>(instruction.value)));
311 break;
312 case Token::opLT:
313 case Token::opLE:
314 case Token::opGE:
315 case Token::opGT:
316 case Token::opEQ:
317 case Token::opNE:
318 stack.push(comparatorInstruction(instruction, frame));
319 break;
320 default:
321 unrecognizedExpression();
322 }
323 }
324 if (stack.size() != 1)
325 evalError("atom expression compiler error - stack over/underflow");
326
327 return stack.top();
328 }
329
330 SelectionSet SelectionEvaluator::comparatorInstruction(
331 const Token& instruction) {
332 int comparator = instruction.tok;
333 int property = instruction.intValue;
334 float comparisonValue = std::any_cast<float>(instruction.value);
335 SelectionSet bs = createSelectionSets();
336 bs.clearAll();
337
338 SimInfo::MoleculeIterator mi;
339 Molecule* mol;
340 Molecule::AtomIterator ai;
341 Atom* atom;
342 Molecule::RigidBodyIterator rbIter;
343 RigidBody* rb;
344
345 for (mol = info->beginMolecule(mi); mol != NULL;
346 mol = info->nextMolecule(mi)) {
347 compareProperty(mol, bs, property, comparator, comparisonValue);
348
349 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
350 compareProperty(atom, bs, property, comparator, comparisonValue);
351 }
352
353 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
354 rb = mol->nextRigidBody(rbIter)) {
355 compareProperty(rb, bs, property, comparator, comparisonValue);
356 }
357 }
358
359 return bs.parallelReduce();
360 }
361
362 SelectionSet SelectionEvaluator::comparatorInstruction(
363 const Token& instruction, int frame) {
364 int comparator = instruction.tok;
365 int property = instruction.intValue;
366 float comparisonValue = std::any_cast<float>(instruction.value);
367 SelectionSet bs = createSelectionSets();
368 bs.clearAll();
369
370 SimInfo::MoleculeIterator mi;
371 Molecule* mol;
372 Molecule::AtomIterator ai;
373 Atom* atom;
374 Molecule::RigidBodyIterator rbIter;
375 RigidBody* rb;
376
377 for (mol = info->beginMolecule(mi); mol != NULL;
378 mol = info->nextMolecule(mi)) {
379 compareProperty(mol, bs, property, comparator, comparisonValue, frame);
380
381 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
382 compareProperty(atom, bs, property, comparator, comparisonValue, frame);
383 }
384
385 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
386 rb = mol->nextRigidBody(rbIter)) {
387 compareProperty(rb, bs, property, comparator, comparisonValue, frame);
388 }
389 }
390
391 return bs.parallelReduce();
392 }
393
394 void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs,
395 int property, int comparator,
396 float comparisonValue) {
397 RealType propertyValue = 0.0;
398 Vector3d pos;
399
400 switch (property) {
401 case Token::atomno:
402 if (sd->isAtom()) {
403 Atom* atom = static_cast<Atom*>(sd);
404 propertyValue = atom->getGlobalIndex();
405 }
406 break;
407 case Token::mass:
408 propertyValue = sd->getMass();
409 break;
410 case Token::charge:
411 if (sd->isAtom()) {
412 Atom* atom = static_cast<Atom*>(sd);
413 propertyValue = getCharge(atom);
414 } else if (sd->isRigidBody()) {
415 RigidBody* rb = static_cast<RigidBody*>(sd);
416 RigidBody::AtomIterator ai;
417 Atom* atom;
418 for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
419 propertyValue += getCharge(atom);
420 }
421 }
422 break;
423 case Token::x:
424 propertyValue = sd->getPos().x();
425 break;
426 case Token::y:
427 propertyValue = sd->getPos().y();
428 break;
429 case Token::z:
430 propertyValue = sd->getPos().z();
431 break;
432 case Token::wrappedX:
433 pos = sd->getPos();
434 info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
435 propertyValue = pos.x();
436 break;
437 case Token::wrappedY:
438 pos = sd->getPos();
439 info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
440 propertyValue = pos.y();
441 break;
442 case Token::wrappedZ:
443 pos = sd->getPos();
444 info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
445 propertyValue = pos.z();
446 break;
447 case Token::r:
448 propertyValue = sd->getPos().length();
449 break;
450 default:
451 unrecognizedAtomProperty(property);
452 }
453
454 bool match = false;
455 switch (comparator) {
456 case Token::opLT:
457 match = propertyValue < comparisonValue;
458 break;
459 case Token::opLE:
460 match = propertyValue <= comparisonValue;
461 break;
462 case Token::opGE:
463 match = propertyValue >= comparisonValue;
464 break;
465 case Token::opGT:
466 match = propertyValue > comparisonValue;
467 break;
468 case Token::opEQ:
469 match = fabs(propertyValue - comparisonValue) <= OpenMD::epsilon;
470 break;
471 case Token::opNE:
472 match = propertyValue != comparisonValue;
473 break;
474 }
475
476 if (match) bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex());
477 }
478
479 void SelectionEvaluator::compareProperty(Molecule* mol, SelectionSet& bs,
480 int property, int comparator,
481 float comparisonValue) {
482 RealType propertyValue = 0.0;
483 Vector3d pos;
484
485 switch (property) {
486 case Token::mass:
487 propertyValue = mol->getMass();
488 break;
489 case Token::charge: {
490 Molecule::AtomIterator ai;
491 Atom* atom;
492 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
493 propertyValue += getCharge(atom);
494 }
495 } break;
496 case Token::x:
497 propertyValue = mol->getCom().x();
498 break;
499 case Token::y:
500 propertyValue = mol->getCom().y();
501 break;
502 case Token::z:
503 propertyValue = mol->getCom().z();
504 break;
505 case Token::wrappedX:
506 pos = mol->getCom();
507 info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
508 propertyValue = pos.x();
509 break;
510 case Token::wrappedY:
511 pos = mol->getCom();
512 info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
513 propertyValue = pos.y();
514 break;
515 case Token::wrappedZ:
516 pos = mol->getCom();
517 info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
518 propertyValue = pos.z();
519 break;
520 case Token::r:
521 propertyValue = mol->getCom().length();
522 break;
523 default:
524 unrecognizedMoleculeProperty(property);
525 }
526
527 bool match = false;
528 switch (comparator) {
529 case Token::opLT:
530 match = propertyValue < comparisonValue;
531 break;
532 case Token::opLE:
533 match = propertyValue <= comparisonValue;
534 break;
535 case Token::opGE:
536 match = propertyValue >= comparisonValue;
537 break;
538 case Token::opGT:
539 match = propertyValue > comparisonValue;
540 break;
541 case Token::opEQ:
542 match = fabs(propertyValue - comparisonValue) <= OpenMD::epsilon;
543 break;
544 case Token::opNE:
545 match = propertyValue != comparisonValue;
546 break;
547 }
548
549 if (match) bs.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
550 }
551
552 void SelectionEvaluator::compareProperty(Molecule* mol, SelectionSet& bs,
553 int property, int comparator,
554 float comparisonValue, int frame) {
555 RealType propertyValue = 0.0;
556 Vector3d pos;
557 switch (property) {
558 case Token::mass:
559 propertyValue = mol->getMass();
560 break;
561 case Token::charge: {
562 Molecule::AtomIterator ai;
563 Atom* atom;
564 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
565 propertyValue += getCharge(atom, frame);
566 }
567 } break;
568 case Token::x:
569 propertyValue = mol->getCom(frame).x();
570 break;
571 case Token::y:
572 propertyValue = mol->getCom(frame).y();
573 break;
574 case Token::z:
575 propertyValue = mol->getCom(frame).z();
576 break;
577 case Token::wrappedX:
578 pos = mol->getCom(frame);
579 info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
580 propertyValue = pos.x();
581 break;
582 case Token::wrappedY:
583 pos = mol->getCom(frame);
584 info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
585 propertyValue = pos.y();
586 break;
587 case Token::wrappedZ:
588 pos = mol->getCom(frame);
589 info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
590 propertyValue = pos.z();
591 break;
592
593 case Token::r:
594 propertyValue = mol->getCom(frame).length();
595 break;
596 default:
597 unrecognizedMoleculeProperty(property);
598 }
599
600 bool match = false;
601 switch (comparator) {
602 case Token::opLT:
603 match = propertyValue < comparisonValue;
604 break;
605 case Token::opLE:
606 match = propertyValue <= comparisonValue;
607 break;
608 case Token::opGE:
609 match = propertyValue >= comparisonValue;
610 break;
611 case Token::opGT:
612 match = propertyValue > comparisonValue;
613 break;
614 case Token::opEQ:
615 match = fabs(propertyValue - comparisonValue) <= OpenMD::epsilon;
616 break;
617 case Token::opNE:
618 match = propertyValue != comparisonValue;
619 break;
620 }
621 if (match) bs.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
622 }
623 void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs,
624 int property, int comparator,
625 float comparisonValue, int frame) {
626 RealType propertyValue = 0.0;
627 Vector3d pos;
628 switch (property) {
629 case Token::atomno:
630 if (sd->isAtom()) {
631 Atom* atom = static_cast<Atom*>(sd);
632 propertyValue = atom->getGlobalIndex();
633 }
634 break;
635 case Token::mass:
636 propertyValue = sd->getMass();
637 break;
638 case Token::charge:
639 if (sd->isAtom()) {
640 Atom* atom = static_cast<Atom*>(sd);
641 propertyValue = getCharge(atom, frame);
642 } else if (sd->isRigidBody()) {
643 RigidBody* rb = static_cast<RigidBody*>(sd);
644 RigidBody::AtomIterator ai;
645 Atom* atom;
646 for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
647 propertyValue += getCharge(atom, frame);
648 }
649 }
650 break;
651 case Token::x:
652 propertyValue = sd->getPos(frame).x();
653 break;
654 case Token::y:
655 propertyValue = sd->getPos(frame).y();
656 break;
657 case Token::z:
658 propertyValue = sd->getPos(frame).z();
659 break;
660 case Token::wrappedX:
661 pos = sd->getPos(frame);
662 info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
663 propertyValue = pos.x();
664 break;
665 case Token::wrappedY:
666 pos = sd->getPos(frame);
667 info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
668 propertyValue = pos.y();
669 break;
670 case Token::wrappedZ:
671 pos = sd->getPos(frame);
672 info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
673 propertyValue = pos.z();
674 break;
675
676 case Token::r:
677 propertyValue = sd->getPos(frame).length();
678 break;
679 default:
680 unrecognizedAtomProperty(property);
681 }
682
683 bool match = false;
684 switch (comparator) {
685 case Token::opLT:
686 match = propertyValue < comparisonValue;
687 break;
688 case Token::opLE:
689 match = propertyValue <= comparisonValue;
690 break;
691 case Token::opGE:
692 match = propertyValue >= comparisonValue;
693 break;
694 case Token::opGT:
695 match = propertyValue > comparisonValue;
696 break;
697 case Token::opEQ:
698 match = fabs(propertyValue - comparisonValue) <= OpenMD::epsilon;
699 break;
700 case Token::opNE:
701 match = propertyValue != comparisonValue;
702 break;
703 }
704 if (match) bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex());
705 }
706
707 void SelectionEvaluator::withinInstruction(const Token& instruction,
708 SelectionSet& bs) {
709 std::any withinSpec = instruction.value;
710 float distance(0.0);
711 if (withinSpec.type() == typeid(float)) {
712 distance = std::any_cast<float>(withinSpec);
713 } else if (withinSpec.type() == typeid(int)) {
714 distance = std::any_cast<int>(withinSpec);
715 } else {
716 evalError("casting error in withinInstruction");
717 bs.clearAll();
718 }
719
720 bs = distanceFinder.find(bs, distance);
721 }
722
723 void SelectionEvaluator::withinInstruction(const Token& instruction,
724 SelectionSet& bs, int frame) {
725 std::any withinSpec = instruction.value;
726 float distance(0.0);
727 if (withinSpec.type() == typeid(float)) {
728 distance = std::any_cast<float>(withinSpec);
729 } else if (withinSpec.type() == typeid(int)) {
730 distance = std::any_cast<int>(withinSpec);
731 } else {
732 evalError("casting error in withinInstruction");
733 bs.clearAll();
734 }
735
736 bs = distanceFinder.find(bs, distance, frame);
737 }
738
739 SelectionSet SelectionEvaluator::alphaHullInstruction(
740 const Token& instruction) {
741 SelectionSet bs = createSelectionSets();
742
743 std::any alphaSpec = instruction.value;
744 float alpha(0.0);
745 if (alphaSpec.type() == typeid(float)) {
746 alpha = std::any_cast<float>(alphaSpec);
747 } else if (alphaSpec.type() == typeid(int)) {
748 alpha = std::any_cast<int>(alphaSpec);
749 } else {
750 evalError("casting error in alphaHullInstruction");
751 bs.clearAll();
752 }
753
754 alphaHullFinder.setAlpha(alpha);
755 bs = alphaHullFinder.findHull();
756 surfaceArea_ = alphaHullFinder.getSurfaceArea();
757 hasSurfaceArea_ = true;
758 volume_ = alphaHullFinder.getVolume();
759 hasVolume_ = true;
760
761 return bs.parallelReduce();
762 }
763
764 SelectionSet SelectionEvaluator::alphaHullInstruction(
765 const Token& instruction, int frame) {
766 SelectionSet bs = createSelectionSets();
767
768 std::any alphaSpec = instruction.value;
769 float alpha(0.0);
770 if (alphaSpec.type() == typeid(float)) {
771 alpha = std::any_cast<float>(alphaSpec);
772 } else if (alphaSpec.type() == typeid(int)) {
773 alpha = std::any_cast<int>(alphaSpec);
774 } else {
775 evalError("casting error in alphaHullInstruction");
776 bs.clearAll();
777 }
778
779 alphaHullFinder.setAlpha(alpha);
780 bs = alphaHullFinder.findHull(frame);
781 surfaceArea_ = alphaHullFinder.getSurfaceArea();
782 hasSurfaceArea_ = true;
783 volume_ = alphaHullFinder.getVolume();
784 hasVolume_ = true;
785
786 return bs.parallelReduce();
787 }
788
789 void SelectionEvaluator::define() {
790 assert(statement.size() >= 3);
791
792 std::string variable = std::any_cast<std::string>(statement[1].value);
793
794 variables.insert(
795 VariablesType::value_type(variable, expression(statement, 2)));
796 }
797
798 /** @todo */
799 void SelectionEvaluator::predefine(const std::string& script) {
800 if (compiler.compile("#predefine", script)) {
801 std::vector<std::vector<Token>> aatoken = compiler.getAatokenCompiled();
802 if (aatoken.size() != 1) {
803 evalError("predefinition does not have exactly 1 command:" + script);
804 return;
805 }
806 std::vector<Token> statement = aatoken[0];
807 if (statement.size() > 2) {
808 int tok = statement[1].tok;
809 if (tok == Token::identifier ||
810 (tok & Token::predefinedset) == Token::predefinedset) {
811 std::string variable = std::any_cast<std::string>(statement[1].value);
812 variables.insert(VariablesType::value_type(variable, statement));
813
814 } else {
815 evalError("invalid variable name:" + script);
816 }
817 } else {
818 evalError("bad predefinition length:" + script);
819 }
820
821 } else {
822 evalError("predefined set compile error:" + script +
823 "\ncompile error:" + compiler.getErrorMessage());
824 }
825 }
826
827 void SelectionEvaluator::select(SelectionSet& bs) {
828 bs = expression(statement, 1);
829 }
830
831 void SelectionEvaluator::select(SelectionSet& bs, int frame) {
832 bs = expression(statement, 1, frame);
833 }
834
835 SelectionSet SelectionEvaluator::lookupValue(const std::string& variable) {
836 SelectionSet bs = createSelectionSets();
837 std::map<std::string, std::any>::iterator i = variables.find(variable);
838
839 if (i != variables.end()) {
840 if (i->second.type() == typeid(SelectionSet)) {
841 return std::any_cast<SelectionSet>(i->second);
842 } else if (i->second.type() == typeid(std::vector<Token>)) {
843 bs = expression(std::any_cast<std::vector<Token>>(i->second), 2);
844 i->second = bs; /**@todo fixme */
845 return bs.parallelReduce();
846 }
847 } else {
848 unrecognizedIdentifier(variable);
849 }
850
851 return bs.parallelReduce();
852 }
853
854 SelectionSet SelectionEvaluator::nameInstruction(const std::string& name) {
855 return nameFinder.match(name);
856 }
857
858 bool SelectionEvaluator::containDynamicToken(
859 const std::vector<Token>& tokens) {
860 std::vector<Token>::const_iterator i;
861 for (i = tokens.begin(); i != tokens.end(); ++i) {
862 if (i->tok & Token::dynamic) { return true; }
863 }
864
865 return false;
866 }
867
868 void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
869 variables.clear();
870 // load predefine
871 // predefine();
872 }
873
874 SelectionSet SelectionEvaluator::createSelectionSets() {
875 SelectionSet ss(nObjects);
876 return ss;
877 }
878
879 SelectionSet SelectionEvaluator::evaluate() {
880 SelectionSet bs = createSelectionSets();
881 if (isLoaded_) {
882 pc = 0;
883 instructionDispatchLoop(bs);
884 }
885 return bs.parallelReduce();
886 }
887
888 SelectionSet SelectionEvaluator::evaluate(int frame) {
889 SelectionSet bs = createSelectionSets();
890 if (isLoaded_) {
891 pc = 0;
892 instructionDispatchLoop(bs, frame);
893 }
894 return bs.parallelReduce();
895 }
896
897 SelectionSet SelectionEvaluator::indexInstruction(const std::any& value) {
898 SelectionSet bs = createSelectionSets();
899
900 if (value.type() == typeid(int)) {
901 int index = std::any_cast<int>(value);
902 if (index < 0 ||
903 index >= static_cast<int>(bs.bitsets_[STUNTDOUBLE].size())) {
904 invalidIndex(index);
905 } else {
906 bs = indexFinder.find(index);
907 }
908 } else if (value.type() == typeid(std::pair<int, int>)) {
909 std::pair<int, int> indexRange =
910 std::any_cast<std::pair<int, int>>(value);
911 assert(indexRange.first <= indexRange.second);
912 if (indexRange.first < 0 ||
913 indexRange.second >=
914 static_cast<int>(bs.bitsets_[STUNTDOUBLE].size())) {
915 invalidIndexRange(indexRange);
916 } else {
917 bs = indexFinder.find(indexRange.first, indexRange.second);
918 }
919 }
920
921 return bs.parallelReduce();
922 }
923
924 SelectionSet SelectionEvaluator::allInstruction() {
925 SelectionSet ss = createSelectionSets();
926
927 SimInfo::MoleculeIterator mi;
928 Molecule::AtomIterator ai;
929 Molecule::RigidBodyIterator rbIter;
930 Molecule::BondIterator bondIter;
931 Molecule::BendIterator bendIter;
932 Molecule::TorsionIterator torsionIter;
933 Molecule::InversionIterator inversionIter;
934
935 Molecule* mol;
936 Atom* atom;
937 RigidBody* rb;
938 Bond* bond;
939 Bend* bend;
940 Torsion* torsion;
941 Inversion* inversion;
942
943 // Doing the loop insures that we're actually on this processor.
944
945 for (mol = info->beginMolecule(mi); mol != NULL;
946 mol = info->nextMolecule(mi)) {
947 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
948 ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
949 }
950 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
951 rb = mol->nextRigidBody(rbIter)) {
952 ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
953 }
954 for (bond = mol->beginBond(bondIter); bond != NULL;
955 bond = mol->nextBond(bondIter)) {
956 ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
957 }
958 for (bend = mol->beginBend(bendIter); bend != NULL;
959 bend = mol->nextBend(bendIter)) {
960 ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
961 }
962 for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
963 torsion = mol->nextTorsion(torsionIter)) {
964 ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
965 }
966 for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
967 inversion = mol->nextInversion(inversionIter)) {
968 ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
969 }
970 ss.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
971 }
972
973 return ss;
974 }
975
976 SelectionSet SelectionEvaluator::hull() {
977 SelectionSet bs = createSelectionSets();
978
979 bs = hullFinder.findHull();
980 surfaceArea_ = hullFinder.getSurfaceArea();
981 hasSurfaceArea_ = true;
982 volume_ = hullFinder.getVolume();
983 hasVolume_ = true;
984
985 return bs.parallelReduce();
986 }
987
988 SelectionSet SelectionEvaluator::hull(int frame) {
989 SelectionSet bs = createSelectionSets();
990
991 bs = hullFinder.findHull(frame);
992
993 return bs.parallelReduce();
994 }
995
996 RealType SelectionEvaluator::getCharge(Atom* atom) {
997 RealType charge = 0.0;
998 AtomType* atomType = atom->getAtomType();
999
1000 FixedChargeAdapter fca = FixedChargeAdapter(atomType);
1001 if (fca.isFixedCharge()) { charge = fca.getCharge(); }
1002
1003 FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
1004 if (fqa.isFluctuatingCharge()) { charge += atom->getFlucQPos(); }
1005 return charge;
1006 }
1007
1008 RealType SelectionEvaluator::getCharge(Atom* atom, int frame) {
1009 RealType charge = 0.0;
1010 AtomType* atomType = atom->getAtomType();
1011
1012 FixedChargeAdapter fca = FixedChargeAdapter(atomType);
1013 if (fca.isFixedCharge()) { charge = fca.getCharge(); }
1014
1015 FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
1016 if (fqa.isFluctuatingCharge()) { charge += atom->getFlucQPos(frame); }
1017 return charge;
1018 }
1019} // namespace OpenMD
RigidBody.hpp
OpenMD::SimInfo
One of the heavy-weight classes of OpenMD, SimInfo maintains objects and variables relating to the cu...
Definition SimInfo.hpp:96
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:63
OpenMD::INVERSION
@ INVERSION
Inversions.
Definition SelectionSet.hpp:68
OpenMD::STUNTDOUBLE
@ STUNTDOUBLE
StuntDoubles (Atoms & RigidBodies).
Definition SelectionSet.hpp:64
OpenMD::TORSION
@ TORSION
Torsions.
Definition SelectionSet.hpp:67
OpenMD::BEND
@ BEND
Bends.
Definition SelectionSet.hpp:66
OpenMD::BOND
@ BOND
Bonds.
Definition SelectionSet.hpp:65
OpenMD::MOLECULE
@ MOLECULE
Molecules.
Definition SelectionSet.hpp:69
OpenMD::distance
Real distance(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the distance between two DynamicVectors.
Definition DynamicVector.hpp:470