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] Vardeman & Gezelter, in progress (2009). |
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> |
46 |
|
#include "primitives/DirectionalAtom.hpp" |
47 |
|
#include "primitives/RigidBody.hpp" |
48 |
|
#include "primitives/Molecule.hpp" |
49 |
< |
#include "io/basic_ifstrstream.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), indexFinder(info), |
58 |
> |
: info(si), nameFinder(info), distanceFinder(info), hullFinder(info), |
59 |
> |
indexFinder(info), |
60 |
|
isLoaded_(false){ |
61 |
|
nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies(); |
62 |
|
} |
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 (int pc = pcStart; pc < code.size(); ++pc) { |
177 |
> |
for (unsigned int pc = pcStart; pc < code.size(); ++pc) { |
178 |
|
Token instruction = code[pc]; |
179 |
|
|
180 |
|
switch (instruction.tok) { |
183 |
|
case Token::expressionEnd: |
184 |
|
break; |
185 |
|
case Token::all: |
186 |
< |
bs = OpenMDBitSet(nStuntDouble); |
161 |
< |
bs.setAll(); |
186 |
> |
bs = allInstruction(); |
187 |
|
stack.push(bs); |
188 |
|
break; |
189 |
|
case Token::none: |
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; |
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); |
220 |
– |
float propertyValue; |
319 |
|
OpenMDBitSet bs(nStuntDouble); |
320 |
|
bs.clearAll(); |
321 |
|
|
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) { |
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 |
|
} |
430 |
|
match = propertyValue != comparisonValue; |
431 |
|
break; |
432 |
|
} |
433 |
< |
if (match) |
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 |
|
|
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); |
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 |
|
|
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 |
|
|
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; |
714 |
> |
RealType charge = 0.0; |
715 |
|
AtomType* atomType = atom->getAtomType(); |
447 |
– |
if (atomType->isCharge()) { |
448 |
– |
GenericData* data = atomType->getPropertyByName("Charge"); |
449 |
– |
if (data != NULL) { |
450 |
– |
DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data); |
716 |
|
|
717 |
< |
if (doubleData != NULL) { |
718 |
< |
charge = doubleData->getData(); |
719 |
< |
|
455 |
< |
} else { |
456 |
< |
sprintf( painCave.errMsg, |
457 |
< |
"Can not cast GenericData to DoubleGenericData\n"); |
458 |
< |
painCave.severity = OPENMD_ERROR; |
459 |
< |
painCave.isFatal = 1; |
460 |
< |
simError(); |
461 |
< |
} |
462 |
< |
} |
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 |
|
|