ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/selection/SelectionEvaluator.cpp
Revision: 1831
Committed: Thu Jan 10 14:06:34 2013 UTC (12 years, 3 months ago) by gezelter
File size: 21507 byte(s)
Log Message:
Merging trunk changes 1812:1830 into development branch.

File Contents

# Content
1 /*
2 * Copyright (c) 2005 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. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
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.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
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, 24107 (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 #include <stack>
44 #include "selection/SelectionEvaluator.hpp"
45 #include "primitives/Atom.hpp"
46 #include "primitives/DirectionalAtom.hpp"
47 #include "primitives/RigidBody.hpp"
48 #include "primitives/Molecule.hpp"
49 #include "io/ifstrstream.hpp"
50 #include "types/FixedChargeAdapter.hpp"
51 #include "types/FluctuatingChargeAdapter.hpp"
52
53
54 namespace OpenMD {
55
56
57 SelectionEvaluator::SelectionEvaluator(SimInfo* si)
58 : info(si), nameFinder(info), distanceFinder(info), hullFinder(info),
59 indexFinder(info),
60 isLoaded_(false){
61 nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies();
62 }
63
64 bool SelectionEvaluator::loadScript(const std::string& filename,
65 const std::string& script) {
66 clearDefinitionsAndLoadPredefined();
67 this->filename = filename;
68 this->script = script;
69 if (! compiler.compile(filename, script)) {
70 error = true;
71 errorMessage = compiler.getErrorMessage();
72
73 sprintf( painCave.errMsg,
74 "SelectionCompiler Error: %s\n", errorMessage.c_str());
75 painCave.severity = OPENMD_ERROR;
76 painCave.isFatal = 1;
77 simError();
78 return false;
79 }
80
81 pc = 0;
82 aatoken = compiler.getAatokenCompiled();
83 linenumbers = compiler.getLineNumbers();
84 lineIndices = compiler.getLineIndices();
85
86 std::vector<std::vector<Token> >::const_iterator i;
87
88 isDynamic_ = false;
89 for (i = aatoken.begin(); i != aatoken.end(); ++i) {
90 if (containDynamicToken(*i)) {
91 isDynamic_ = true;
92 break;
93 }
94 }
95
96 isLoaded_ = true;
97 return true;
98 }
99
100 void SelectionEvaluator::clearState() {
101 error = false;
102 errorMessage = "";
103 }
104
105 bool SelectionEvaluator::loadScriptString(const std::string& script) {
106 clearState();
107 return loadScript("", script);
108 }
109
110 bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
111 clearState();
112 return loadScriptFileInternal(filename);
113 }
114
115 bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) {
116 ifstrstream ifs(filename.c_str());
117 if (!ifs.is_open()) {
118 return false;
119 }
120
121 const int bufferSize = 65535;
122 char buffer[bufferSize];
123 std::string script;
124 while(ifs.getline(buffer, bufferSize)) {
125 script += buffer;
126 }
127 return loadScript(filename, script);
128 }
129
130 void SelectionEvaluator::instructionDispatchLoop(OpenMDBitSet& bs){
131
132 while ( pc < aatoken.size()) {
133 statement = aatoken[pc++];
134 statementLength = statement.size();
135 Token token = statement[0];
136 switch (token.tok) {
137 case Token::define:
138 define();
139 break;
140 case Token::select:
141 select(bs);
142 break;
143 default:
144 unrecognizedCommand(token);
145 return;
146 }
147 }
148
149 }
150
151 void SelectionEvaluator::instructionDispatchLoop(OpenMDBitSet& bs, int frame){
152
153 while ( pc < aatoken.size()) {
154 statement = aatoken[pc++];
155 statementLength = statement.size();
156 Token token = statement[0];
157 switch (token.tok) {
158 case Token::define:
159 define();
160 break;
161 case Token::select:
162 select(bs, frame);
163 break;
164 default:
165 unrecognizedCommand(token);
166 return;
167 }
168 }
169
170 }
171
172 OpenMDBitSet SelectionEvaluator::expression(const std::vector<Token>& code,
173 int pcStart) {
174 OpenMDBitSet bs;
175 std::stack<OpenMDBitSet> stack;
176
177 for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
178 Token instruction = code[pc];
179
180 switch (instruction.tok) {
181 case Token::expressionBegin:
182 break;
183 case Token::expressionEnd:
184 break;
185 case Token::all:
186 bs = allInstruction();
187 stack.push(bs);
188 break;
189 case Token::none:
190 bs = OpenMDBitSet(nStuntDouble);
191 stack.push(bs);
192 break;
193 case Token::opOr:
194 bs = stack.top();
195 stack.pop();
196 stack.top() |= bs;
197 break;
198 case Token::opAnd:
199 bs = stack.top();
200 stack.pop();
201 stack.top() &= bs;
202 break;
203 case Token::opNot:
204 stack.top().flip();
205 break;
206 case Token::within:
207 withinInstruction(instruction, stack.top());
208 break;
209 case Token::hull:
210 stack.push(hull());
211 break;
212 //case Token::selected:
213 // stack.push(getSelectionSet());
214 // break;
215 case Token::name:
216 stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
217 break;
218 case Token::index:
219 stack.push(indexInstruction(instruction.value));
220 break;
221 case Token::identifier:
222 stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
223 break;
224 case Token::opLT:
225 case Token::opLE:
226 case Token::opGE:
227 case Token::opGT:
228 case Token::opEQ:
229 case Token::opNE:
230 stack.push(comparatorInstruction(instruction));
231 break;
232 default:
233 unrecognizedExpression();
234 }
235 }
236 if (stack.size() != 1)
237 evalError("atom expression compiler error - stack over/underflow");
238
239 return stack.top();
240 }
241
242
243 OpenMDBitSet SelectionEvaluator::expression(const std::vector<Token>& code,
244 int pcStart, int frame) {
245 OpenMDBitSet bs;
246 std::stack<OpenMDBitSet> stack;
247
248 for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
249 Token instruction = code[pc];
250
251 switch (instruction.tok) {
252 case Token::expressionBegin:
253 break;
254 case Token::expressionEnd:
255 break;
256 case Token::all:
257 bs = allInstruction();
258 stack.push(bs);
259 break;
260 case Token::none:
261 bs = OpenMDBitSet(nStuntDouble);
262 stack.push(bs);
263 break;
264 case Token::opOr:
265 bs = stack.top();
266 stack.pop();
267 stack.top() |= bs;
268 break;
269 case Token::opAnd:
270 bs = stack.top();
271 stack.pop();
272 stack.top() &= bs;
273 break;
274 case Token::opNot:
275 stack.top().flip();
276 break;
277 case Token::within:
278 withinInstruction(instruction, stack.top(), frame);
279 break;
280 case Token::hull:
281 stack.push(hull(frame));
282 break;
283 //case Token::selected:
284 // stack.push(getSelectionSet());
285 // break;
286 case Token::name:
287 stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
288 break;
289 case Token::index:
290 stack.push(indexInstruction(instruction.value));
291 break;
292 case Token::identifier:
293 stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
294 break;
295 case Token::opLT:
296 case Token::opLE:
297 case Token::opGE:
298 case Token::opGT:
299 case Token::opEQ:
300 case Token::opNE:
301 stack.push(comparatorInstruction(instruction, frame));
302 break;
303 default:
304 unrecognizedExpression();
305 }
306 }
307 if (stack.size() != 1)
308 evalError("atom expression compiler error - stack over/underflow");
309
310 return stack.top();
311 }
312
313
314
315 OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
316 int comparator = instruction.tok;
317 int property = instruction.intValue;
318 float comparisonValue = boost::any_cast<float>(instruction.value);
319 OpenMDBitSet bs(nStuntDouble);
320 bs.clearAll();
321
322 SimInfo::MoleculeIterator mi;
323 Molecule* mol;
324 Molecule::AtomIterator ai;
325 Atom* atom;
326 Molecule::RigidBodyIterator rbIter;
327 RigidBody* rb;
328
329 for (mol = info->beginMolecule(mi); mol != NULL;
330 mol = info->nextMolecule(mi)) {
331
332 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
333 compareProperty(atom, bs, property, comparator, comparisonValue);
334 }
335
336 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
337 rb = mol->nextRigidBody(rbIter)) {
338 compareProperty(rb, bs, property, comparator, comparisonValue);
339 }
340 }
341
342 return bs;
343 }
344
345 OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction, int frame) {
346 int comparator = instruction.tok;
347 int property = instruction.intValue;
348 float comparisonValue = boost::any_cast<float>(instruction.value);
349 OpenMDBitSet bs(nStuntDouble);
350 bs.clearAll();
351
352 SimInfo::MoleculeIterator mi;
353 Molecule* mol;
354 Molecule::AtomIterator ai;
355 Atom* atom;
356 Molecule::RigidBodyIterator rbIter;
357 RigidBody* rb;
358
359 for (mol = info->beginMolecule(mi); mol != NULL;
360 mol = info->nextMolecule(mi)) {
361
362 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
363 compareProperty(atom, bs, property, comparator, comparisonValue, frame);
364 }
365
366 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
367 rb = mol->nextRigidBody(rbIter)) {
368 compareProperty(rb, bs, property, comparator, comparisonValue, frame);
369 }
370 }
371
372 return bs;
373 }
374
375 void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs,
376 int property, int comparator,
377 float comparisonValue) {
378 RealType propertyValue = 0.0;
379 switch (property) {
380 case Token::mass:
381 propertyValue = sd->getMass();
382 break;
383 case Token::charge:
384 if (sd->isAtom()){
385 Atom* atom = static_cast<Atom*>(sd);
386 propertyValue = getCharge(atom);
387 } else if (sd->isRigidBody()) {
388 RigidBody* rb = static_cast<RigidBody*>(sd);
389 RigidBody::AtomIterator ai;
390 Atom* atom;
391 for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
392 propertyValue+= getCharge(atom);
393 }
394 }
395 break;
396 case Token::x:
397 propertyValue = sd->getPos().x();
398 break;
399 case Token::y:
400 propertyValue = sd->getPos().y();
401 break;
402 case Token::z:
403 propertyValue = sd->getPos().z();
404 break;
405 case Token::r:
406 propertyValue = sd->getPos().length();
407 break;
408 default:
409 unrecognizedAtomProperty(property);
410 }
411
412 bool match = false;
413 switch (comparator) {
414 case Token::opLT:
415 match = propertyValue < comparisonValue;
416 break;
417 case Token::opLE:
418 match = propertyValue <= comparisonValue;
419 break;
420 case Token::opGE:
421 match = propertyValue >= comparisonValue;
422 break;
423 case Token::opGT:
424 match = propertyValue > comparisonValue;
425 break;
426 case Token::opEQ:
427 match = propertyValue == comparisonValue;
428 break;
429 case Token::opNE:
430 match = propertyValue != comparisonValue;
431 break;
432 }
433 if (match)
434 bs.setBitOn(sd->getGlobalIndex());
435
436
437 }
438
439 void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs,
440 int property, int comparator,
441 float comparisonValue, int frame) {
442 RealType propertyValue = 0.0;
443 switch (property) {
444 case Token::mass:
445 propertyValue = sd->getMass();
446 break;
447 case Token::charge:
448 if (sd->isAtom()){
449 Atom* atom = static_cast<Atom*>(sd);
450 propertyValue = getCharge(atom,frame);
451 } else if (sd->isRigidBody()) {
452 RigidBody* rb = static_cast<RigidBody*>(sd);
453 RigidBody::AtomIterator ai;
454 Atom* atom;
455 for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
456 propertyValue+= getCharge(atom,frame);
457 }
458 }
459 break;
460 case Token::x:
461 propertyValue = sd->getPos(frame).x();
462 break;
463 case Token::y:
464 propertyValue = sd->getPos(frame).y();
465 break;
466 case Token::z:
467 propertyValue = sd->getPos(frame).z();
468 break;
469 case Token::r:
470 propertyValue = sd->getPos(frame).length();
471 break;
472 default:
473 unrecognizedAtomProperty(property);
474 }
475
476 bool match = false;
477 switch (comparator) {
478 case Token::opLT:
479 match = propertyValue < comparisonValue;
480 break;
481 case Token::opLE:
482 match = propertyValue <= comparisonValue;
483 break;
484 case Token::opGE:
485 match = propertyValue >= comparisonValue;
486 break;
487 case Token::opGT:
488 match = propertyValue > comparisonValue;
489 break;
490 case Token::opEQ:
491 match = propertyValue == comparisonValue;
492 break;
493 case Token::opNE:
494 match = propertyValue != comparisonValue;
495 break;
496 }
497 if (match)
498 bs.setBitOn(sd->getGlobalIndex());
499
500
501 }
502
503 void SelectionEvaluator::withinInstruction(const Token& instruction,
504 OpenMDBitSet& bs){
505
506 boost::any withinSpec = instruction.value;
507 float distance;
508 if (withinSpec.type() == typeid(float)){
509 distance = boost::any_cast<float>(withinSpec);
510 } else if (withinSpec.type() == typeid(int)) {
511 distance = boost::any_cast<int>(withinSpec);
512 } else {
513 evalError("casting error in withinInstruction");
514 bs.clearAll();
515 }
516
517 bs = distanceFinder.find(bs, distance);
518 }
519
520 void SelectionEvaluator::withinInstruction(const Token& instruction,
521 OpenMDBitSet& bs, int frame){
522
523 boost::any withinSpec = instruction.value;
524 float distance;
525 if (withinSpec.type() == typeid(float)){
526 distance = boost::any_cast<float>(withinSpec);
527 } else if (withinSpec.type() == typeid(int)) {
528 distance = boost::any_cast<int>(withinSpec);
529 } else {
530 evalError("casting error in withinInstruction");
531 bs.clearAll();
532 }
533
534 bs = distanceFinder.find(bs, distance, frame);
535 }
536
537 void SelectionEvaluator::define() {
538 assert(statement.size() >= 3);
539
540 std::string variable = boost::any_cast<std::string>(statement[1].value);
541
542 variables.insert(VariablesType::value_type(variable,
543 expression(statement, 2)));
544 }
545
546
547 /** @todo */
548 void SelectionEvaluator::predefine(const std::string& script) {
549
550 if (compiler.compile("#predefine", script)) {
551 std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
552 if (aatoken.size() != 1) {
553 evalError("predefinition does not have exactly 1 command:"
554 + script);
555 return;
556 }
557 std::vector<Token> statement = aatoken[0];
558 if (statement.size() > 2) {
559 int tok = statement[1].tok;
560 if (tok == Token::identifier ||
561 (tok & Token::predefinedset) == Token::predefinedset) {
562 std::string variable = boost::any_cast<std::string>(statement[1].value);
563 variables.insert(VariablesType::value_type(variable, statement));
564
565 } else {
566 evalError("invalid variable name:" + script);
567 }
568 }else {
569 evalError("bad predefinition length:" + script);
570 }
571
572 } else {
573 evalError("predefined set compile error:" + script +
574 "\ncompile error:" + compiler.getErrorMessage());
575 }
576 }
577
578 void SelectionEvaluator::select(OpenMDBitSet& bs){
579 bs = expression(statement, 1);
580 }
581
582 void SelectionEvaluator::select(OpenMDBitSet& bs, int frame){
583 bs = expression(statement, 1, frame);
584 }
585
586 OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){
587
588 OpenMDBitSet bs(nStuntDouble);
589 std::map<std::string, boost::any>::iterator i = variables.find(variable);
590
591 if (i != variables.end()) {
592 if (i->second.type() == typeid(OpenMDBitSet)) {
593 return boost::any_cast<OpenMDBitSet>(i->second);
594 } else if (i->second.type() == typeid(std::vector<Token>)){
595 bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
596 i->second = bs; /**@todo fixme */
597 return bs;
598 }
599 } else {
600 unrecognizedIdentifier(variable);
601 }
602
603 return bs;
604 }
605
606 OpenMDBitSet SelectionEvaluator::nameInstruction(const std::string& name){
607 return nameFinder.match(name);
608 }
609
610 bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
611 std::vector<Token>::const_iterator i;
612 for (i = tokens.begin(); i != tokens.end(); ++i) {
613 if (i->tok & Token::dynamic) {
614 return true;
615 }
616 }
617
618 return false;
619 }
620
621 void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
622 variables.clear();
623 //load predefine
624 //predefine();
625 }
626
627 OpenMDBitSet SelectionEvaluator::evaluate() {
628 OpenMDBitSet bs(nStuntDouble);
629 if (isLoaded_) {
630 pc = 0;
631 instructionDispatchLoop(bs);
632 }
633 return bs;
634 }
635
636 OpenMDBitSet SelectionEvaluator::evaluate(int frame) {
637 OpenMDBitSet bs(nStuntDouble);
638 if (isLoaded_) {
639 pc = 0;
640 instructionDispatchLoop(bs, frame);
641 }
642 return bs;
643 }
644
645 OpenMDBitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
646 OpenMDBitSet bs(nStuntDouble);
647
648 if (value.type() == typeid(int)) {
649 int index = boost::any_cast<int>(value);
650 if (index < 0 || index >= bs.size()) {
651 invalidIndex(index);
652 } else {
653 bs = indexFinder.find(index);
654 }
655 } else if (value.type() == typeid(std::pair<int, int>)) {
656 std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
657 assert(indexRange.first <= indexRange.second);
658 if (indexRange.first < 0 || indexRange.second >= bs.size()) {
659 invalidIndexRange(indexRange);
660 }else {
661 bs = indexFinder.find(indexRange.first, indexRange.second);
662 }
663 }
664
665 return bs;
666 }
667
668 OpenMDBitSet SelectionEvaluator::allInstruction() {
669 OpenMDBitSet bs(nStuntDouble);
670
671 SimInfo::MoleculeIterator mi;
672 Molecule* mol;
673 Molecule::AtomIterator ai;
674 Atom* atom;
675 Molecule::RigidBodyIterator rbIter;
676 RigidBody* rb;
677
678 // Doing the loop insures that we're actually on this processor.
679
680 for (mol = info->beginMolecule(mi); mol != NULL;
681 mol = info->nextMolecule(mi)) {
682
683 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
684 bs.setBitOn(atom->getGlobalIndex());
685 }
686
687 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
688 rb = mol->nextRigidBody(rbIter)) {
689 bs.setBitOn(rb->getGlobalIndex());
690 }
691 }
692
693 return bs;
694 }
695
696 OpenMDBitSet SelectionEvaluator::hull() {
697 OpenMDBitSet bs(nStuntDouble);
698
699 bs = hullFinder.findHull();
700
701 return bs;
702 }
703
704
705 OpenMDBitSet SelectionEvaluator::hull(int frame) {
706 OpenMDBitSet bs(nStuntDouble);
707
708 bs = hullFinder.findHull(frame);
709
710 return bs;
711 }
712
713 RealType SelectionEvaluator::getCharge(Atom* atom) {
714 RealType charge = 0.0;
715 AtomType* atomType = atom->getAtomType();
716
717 FixedChargeAdapter fca = FixedChargeAdapter(atomType);
718 if ( fca.isFixedCharge() ) {
719 charge = fca.getCharge();
720 }
721
722 FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
723 if ( fqa.isFluctuatingCharge() ) {
724 charge += atom->getFlucQPos();
725 }
726 return charge;
727 }
728
729 RealType SelectionEvaluator::getCharge(Atom* atom, int frame) {
730 RealType charge = 0.0;
731 AtomType* atomType = atom->getAtomType();
732
733 FixedChargeAdapter fca = FixedChargeAdapter(atomType);
734 if ( fca.isFixedCharge() ) {
735 charge = fca.getCharge();
736 }
737
738 FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
739 if ( fqa.isFluctuatingCharge() ) {
740 charge += atom->getFlucQPos(frame);
741 }
742 return charge;
743 }
744
745 }

Properties

Name Value
svn:keywords Author Id Revision Date