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 1513 by gezelter, Tue Oct 19 18:40:54 2010 UTC vs.
Revision 1931 by gezelter, Mon Aug 19 19:20:32 2013 UTC

# Line 35 | Line 35
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]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (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>
# 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), hullFinder(info),
59 <      indexFinder(info),
59 >      indexFinder(info), hasSurfaceArea_(false),
60        isLoaded_(false){    
61        nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies();
62      }            
# Line 144 | 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 158 | Line 183 | namespace OpenMD {
183        case Token::expressionEnd:
184          break;
185        case Token::all:
186 <        bs = OpenMDBitSet(nStuntDouble);
162 <        bs.setAll();
186 >        bs = allInstruction();
187          stack.push(bs);            
188          break;
189        case Token::none:
# Line 216 | 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);
224 –    float propertyValue;
319      OpenMDBitSet bs(nStuntDouble);
320      bs.clearAll();
321      
# Line 248 | 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) {
378      RealType propertyValue = 0.0;
379 +    Vector3d pos;
380 +
381      switch (property) {
382      case Token::mass:
383        propertyValue = sd->getMass();
# Line 278 | Line 404 | namespace OpenMD {
404      case Token::z:
405        propertyValue = sd->getPos().z();
406        break;
407 +    case Token::wrappedX:    
408 +      pos = sd->getPos();
409 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
410 +      propertyValue = pos.x();
411 +      break;
412 +    case Token::wrappedY:
413 +      pos = sd->getPos();
414 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
415 +      propertyValue = pos.y();
416 +      break;
417 +    case Token::wrappedZ:
418 +      pos = sd->getPos();
419 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
420 +      propertyValue = pos.z();
421 +      break;
422      case Token::r:
423        propertyValue = sd->getPos().length();
424        break;
# Line 306 | Line 447 | namespace OpenMD {
447        match = propertyValue != comparisonValue;
448        break;
449      }
450 <    if (match)
450 >
451 >    if (match)
452 >      bs.setBitOn(sd->getGlobalIndex());
453 >
454 >    
455 >  }
456 >
457 >  void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs,
458 >                                           int property, int comparator,
459 >                                           float comparisonValue, int frame) {
460 >    RealType propertyValue = 0.0;
461 >    Vector3d pos;
462 >    switch (property) {
463 >    case Token::mass:
464 >      propertyValue = sd->getMass();
465 >      break;
466 >    case Token::charge:
467 >      if (sd->isAtom()){
468 >        Atom* atom = static_cast<Atom*>(sd);
469 >        propertyValue = getCharge(atom,frame);
470 >      } else if (sd->isRigidBody()) {
471 >        RigidBody* rb = static_cast<RigidBody*>(sd);
472 >        RigidBody::AtomIterator ai;
473 >        Atom* atom;
474 >        for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
475 >          propertyValue+=  getCharge(atom,frame);
476 >        }
477 >      }
478 >      break;
479 >    case Token::x:
480 >      propertyValue = sd->getPos(frame).x();
481 >      break;
482 >    case Token::y:
483 >      propertyValue = sd->getPos(frame).y();
484 >      break;
485 >    case Token::z:
486 >      propertyValue = sd->getPos(frame).z();
487 >      break;
488 >    case Token::wrappedX:    
489 >      pos = sd->getPos(frame);
490 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
491 >      propertyValue = pos.x();
492 >      break;
493 >    case Token::wrappedY:
494 >      pos = sd->getPos(frame);
495 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
496 >      propertyValue = pos.y();
497 >      break;
498 >    case Token::wrappedZ:
499 >      pos = sd->getPos(frame);
500 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
501 >      propertyValue = pos.z();
502 >      break;
503 >
504 >    case Token::r:
505 >      propertyValue = sd->getPos(frame).length();
506 >      break;
507 >    default:
508 >      unrecognizedAtomProperty(property);
509 >    }
510 >        
511 >    bool match = false;
512 >    switch (comparator) {
513 >    case Token::opLT:
514 >      match = propertyValue < comparisonValue;
515 >      break;
516 >    case Token::opLE:
517 >      match = propertyValue <= comparisonValue;
518 >      break;
519 >    case Token::opGE:
520 >      match = propertyValue >= comparisonValue;
521 >      break;
522 >    case Token::opGT:
523 >      match = propertyValue > comparisonValue;
524 >      break;
525 >    case Token::opEQ:
526 >      match = propertyValue == comparisonValue;
527 >      break;
528 >    case Token::opNE:
529 >      match = propertyValue != comparisonValue;
530 >      break;
531 >    }
532 >    if (match)
533        bs.setBitOn(sd->getGlobalIndex());
534 +
535      
536    }
537  
# Line 327 | Line 551 | namespace OpenMD {
551      
552      bs = distanceFinder.find(bs, distance);            
553    }
554 +
555 +  void SelectionEvaluator::withinInstruction(const Token& instruction,
556 +                                             OpenMDBitSet& bs, int frame){
557 +    
558 +    boost::any withinSpec = instruction.value;
559 +    float distance;
560 +    if (withinSpec.type() == typeid(float)){
561 +      distance = boost::any_cast<float>(withinSpec);
562 +    } else if (withinSpec.type() == typeid(int)) {
563 +      distance = boost::any_cast<int>(withinSpec);    
564 +    } else {
565 +      evalError("casting error in withinInstruction");
566 +      bs.clearAll();
567 +    }
568 +    
569 +    bs = distanceFinder.find(bs, distance, frame);            
570 +  }
571    
572    void SelectionEvaluator::define() {
573      assert(statement.size() >= 3);
# Line 372 | Line 613 | namespace OpenMD {
613    void SelectionEvaluator::select(OpenMDBitSet& bs){
614      bs = expression(statement, 1);
615    }
616 +
617 +  void SelectionEvaluator::select(OpenMDBitSet& bs, int frame){
618 +    bs = expression(statement, 1, frame);
619 +  }
620    
621    OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){
622      
# Line 420 | Line 665 | namespace OpenMD {
665        pc = 0;
666        instructionDispatchLoop(bs);
667      }
668 +    return bs;
669 +  }
670  
671 +  OpenMDBitSet SelectionEvaluator::evaluate(int frame) {
672 +    OpenMDBitSet bs(nStuntDouble);
673 +    if (isLoaded_) {
674 +      pc = 0;
675 +      instructionDispatchLoop(bs, frame);
676 +    }
677      return bs;
678    }
679  
# Line 447 | Line 700 | namespace OpenMD {
700      return bs;
701    }
702  
703 +  OpenMDBitSet SelectionEvaluator::allInstruction() {
704 +    OpenMDBitSet bs(nStuntDouble);
705  
706 +    SimInfo::MoleculeIterator mi;
707 +    Molecule* mol;
708 +    Molecule::AtomIterator ai;
709 +    Atom* atom;
710 +    Molecule::RigidBodyIterator rbIter;
711 +    RigidBody* rb;
712 +
713 +    // Doing the loop insures that we're actually on this processor.
714 +
715 +    for (mol = info->beginMolecule(mi); mol != NULL;
716 +         mol = info->nextMolecule(mi)) {
717 +
718 +      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
719 +        bs.setBitOn(atom->getGlobalIndex());
720 +      }
721 +    
722 +      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
723 +           rb = mol->nextRigidBody(rbIter)) {
724 +        bs.setBitOn(rb->getGlobalIndex());
725 +      }
726 +    }
727 +
728 +    return bs;
729 +  }
730 +
731    OpenMDBitSet SelectionEvaluator::hull() {
732      OpenMDBitSet bs(nStuntDouble);
733      
734      bs = hullFinder.findHull();
735 <    
735 >    surfaceArea_ = hullFinder.getSurfaceArea();
736 >    hasSurfaceArea_ = true;
737      return bs;
738    }
739  
740  
741 +  OpenMDBitSet SelectionEvaluator::hull(int frame) {
742 +    OpenMDBitSet bs(nStuntDouble);
743 +    
744 +    bs = hullFinder.findHull(frame);
745  
746 +    return bs;
747 +  }
748 +
749    RealType SelectionEvaluator::getCharge(Atom* atom) {
750 <    RealType charge =0.0;
750 >    RealType charge = 0.0;
751      AtomType* atomType = atom->getAtomType();
464 –    if (atomType->isCharge()) {
465 –      GenericData* data = atomType->getPropertyByName("Charge");
466 –      if (data != NULL) {
467 –        DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data);
752  
753 <        if (doubleData != NULL) {
754 <          charge = doubleData->getData();
755 <
472 <        } else {
473 <          sprintf( painCave.errMsg,
474 <                   "Can not cast GenericData to DoubleGenericData\n");
475 <          painCave.severity = OPENMD_ERROR;
476 <          painCave.isFatal = 1;
477 <          simError();          
478 <        }
479 <      }
753 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
754 >    if ( fca.isFixedCharge() ) {
755 >      charge = fca.getCharge();
756      }
757 +    
758 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
759 +    if ( fqa.isFluctuatingCharge() ) {
760 +      charge += atom->getFlucQPos();
761 +    }
762 +    return charge;
763 +  }
764  
765 +  RealType SelectionEvaluator::getCharge(Atom* atom, int frame) {
766 +    RealType charge = 0.0;
767 +    AtomType* atomType = atom->getAtomType();    
768 +
769 +    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
770 +    if ( fca.isFixedCharge() ) {
771 +      charge = fca.getCharge();
772 +    }
773 +    
774 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
775 +    if ( fqa.isFluctuatingCharge() ) {
776 +      charge += atom->getFlucQPos(frame);
777 +    }
778      return charge;
779    }
780  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines