ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/selection/SelectionEvaluator.cpp
(Generate patch)

Comparing trunk/src/selection/SelectionEvaluator.cpp (file contents):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1816 by gezelter, Fri Dec 7 18:33:33 2012 UTC

# Line 36 | Line 36
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>
# Line 45 | Line 46
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      }            
# Line 143 | Line 148 | namespace OpenMD {
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) {
# Line 157 | Line 183 | namespace OpenMD {
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:
# Line 181 | Line 206 | namespace OpenMD {
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;
# Line 212 | Line 240 | namespace OpenMD {
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      
# Line 244 | Line 342 | namespace OpenMD {
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) {
# Line 274 | Line 402 | namespace OpenMD {
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      }
# Line 299 | Line 430 | namespace OpenMD {
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      
# Line 320 | Line 516 | namespace OpenMD {
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);
# Line 365 | Line 578 | namespace OpenMD {
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      
# Line 413 | Line 630 | namespace OpenMD {
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  
# Line 440 | Line 665 | namespace OpenMD {
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  

Comparing trunk/src/selection/SelectionEvaluator.cpp (property svn:keywords):
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC vs.
Revision 1816 by gezelter, Fri Dec 7 18:33:33 2012 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines