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 1801 by gezelter, Mon Oct 1 18:21:15 2012 UTC vs.
Revision 1903 by gezelter, Tue Jul 16 18:58:08 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).          
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   */
# Line 47 | Line 47
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),
59 >      indexFinder(info), hasSurfaceArea_(false),
60        isLoaded_(false){    
61        nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies();
62      }            
# Line 145 | 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;
# 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;
# Line 247 | 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      switch (property) {
381      case Token::mass:
382        propertyValue = sd->getMass();
# Line 277 | Line 403 | namespace OpenMD {
403      case Token::z:
404        propertyValue = sd->getPos().z();
405        break;
406 +    case Token::wrappedX:    
407 +      pos = sd->getPos();
408 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
409 +      propertyValue = pos.x();
410 +      break;
411 +    case Token::wrappedY:
412 +      pos = sd->getPos();
413 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
414 +      propertyValue = pos.y();
415 +      break;
416 +    case Token::wrappedZ:
417 +      pos = sd->getPos();
418 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
419 +      propertyValue = pos.z();
420 +      break;
421      case Token::r:
422        propertyValue = sd->getPos().length();
423        break;
# Line 305 | Line 446 | namespace OpenMD {
446        match = propertyValue != comparisonValue;
447        break;
448      }
449 <    if (match)
449 >    if (match)
450 >      bs.setBitOn(sd->getGlobalIndex());
451 >
452 >    
453 >  }
454 >
455 >  void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs,
456 >                                           int property, int comparator,
457 >                                           float comparisonValue, int frame) {
458 >    RealType propertyValue = 0.0;
459 >    Vector3d pos;
460 >    switch (property) {
461 >    case Token::mass:
462 >      propertyValue = sd->getMass();
463 >      break;
464 >    case Token::charge:
465 >      if (sd->isAtom()){
466 >        Atom* atom = static_cast<Atom*>(sd);
467 >        propertyValue = getCharge(atom,frame);
468 >      } else if (sd->isRigidBody()) {
469 >        RigidBody* rb = static_cast<RigidBody*>(sd);
470 >        RigidBody::AtomIterator ai;
471 >        Atom* atom;
472 >        for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
473 >          propertyValue+=  getCharge(atom,frame);
474 >        }
475 >      }
476 >      break;
477 >    case Token::x:
478 >      propertyValue = sd->getPos(frame).x();
479 >      break;
480 >    case Token::y:
481 >      propertyValue = sd->getPos(frame).y();
482 >      break;
483 >    case Token::z:
484 >      propertyValue = sd->getPos(frame).z();
485 >      break;
486 >    case Token::wrappedX:    
487 >      pos = sd->getPos(frame);
488 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
489 >      propertyValue = pos.x();
490 >      break;
491 >    case Token::wrappedY:
492 >      pos = sd->getPos(frame);
493 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
494 >      propertyValue = pos.y();
495 >      break;
496 >    case Token::wrappedZ:
497 >      pos = sd->getPos(frame);
498 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
499 >      propertyValue = pos.z();
500 >      break;
501 >
502 >    case Token::r:
503 >      propertyValue = sd->getPos(frame).length();
504 >      break;
505 >    default:
506 >      unrecognizedAtomProperty(property);
507 >    }
508 >        
509 >    bool match = false;
510 >    switch (comparator) {
511 >    case Token::opLT:
512 >      match = propertyValue < comparisonValue;
513 >      break;
514 >    case Token::opLE:
515 >      match = propertyValue <= comparisonValue;
516 >      break;
517 >    case Token::opGE:
518 >      match = propertyValue >= comparisonValue;
519 >      break;
520 >    case Token::opGT:
521 >      match = propertyValue > comparisonValue;
522 >      break;
523 >    case Token::opEQ:
524 >      match = propertyValue == comparisonValue;
525 >      break;
526 >    case Token::opNE:
527 >      match = propertyValue != comparisonValue;
528 >      break;
529 >    }
530 >    if (match)
531        bs.setBitOn(sd->getGlobalIndex());
532 +
533      
534    }
535  
# Line 325 | Line 548 | namespace OpenMD {
548      }
549      
550      bs = distanceFinder.find(bs, distance);            
551 +  }
552 +
553 +  void SelectionEvaluator::withinInstruction(const Token& instruction,
554 +                                             OpenMDBitSet& bs, int frame){
555 +    
556 +    boost::any withinSpec = instruction.value;
557 +    float distance;
558 +    if (withinSpec.type() == typeid(float)){
559 +      distance = boost::any_cast<float>(withinSpec);
560 +    } else if (withinSpec.type() == typeid(int)) {
561 +      distance = boost::any_cast<int>(withinSpec);    
562 +    } else {
563 +      evalError("casting error in withinInstruction");
564 +      bs.clearAll();
565 +    }
566 +    
567 +    bs = distanceFinder.find(bs, distance, frame);            
568    }
569    
570    void SelectionEvaluator::define() {
# Line 371 | Line 611 | namespace OpenMD {
611    void SelectionEvaluator::select(OpenMDBitSet& bs){
612      bs = expression(statement, 1);
613    }
614 +
615 +  void SelectionEvaluator::select(OpenMDBitSet& bs, int frame){
616 +    bs = expression(statement, 1, frame);
617 +  }
618    
619    OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){
620      
# Line 419 | Line 663 | namespace OpenMD {
663        pc = 0;
664        instructionDispatchLoop(bs);
665      }
666 +    return bs;
667 +  }
668  
669 +  OpenMDBitSet SelectionEvaluator::evaluate(int frame) {
670 +    OpenMDBitSet bs(nStuntDouble);
671 +    if (isLoaded_) {
672 +      pc = 0;
673 +      instructionDispatchLoop(bs, frame);
674 +    }
675      return bs;
676    }
677  
# Line 478 | Line 730 | namespace OpenMD {
730      OpenMDBitSet bs(nStuntDouble);
731      
732      bs = hullFinder.findHull();
733 +    surfaceArea_ = hullFinder.getSurfaceArea();
734 +    hasSurfaceArea_ = true;
735 +    return bs;
736 +  }
737 +
738 +
739 +  OpenMDBitSet SelectionEvaluator::hull(int frame) {
740 +    OpenMDBitSet bs(nStuntDouble);
741      
742 +    bs = hullFinder.findHull(frame);
743 +
744      return bs;
745    }
746  
747    RealType SelectionEvaluator::getCharge(Atom* atom) {
748 <    RealType charge =0.0;
748 >    RealType charge = 0.0;
749      AtomType* atomType = atom->getAtomType();
488 –    if (atomType->isCharge()) {
489 –      GenericData* data = atomType->getPropertyByName("Charge");
490 –      if (data != NULL) {
491 –        DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data);
750  
751 <        if (doubleData != NULL) {
752 <          charge = doubleData->getData();
753 <
496 <        } else {
497 <          sprintf( painCave.errMsg,
498 <                   "Can not cast GenericData to DoubleGenericData\n");
499 <          painCave.severity = OPENMD_ERROR;
500 <          painCave.isFatal = 1;
501 <          simError();          
502 <        }
503 <      }
751 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
752 >    if ( fca.isFixedCharge() ) {
753 >      charge = fca.getCharge();
754      }
755 +    
756 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
757 +    if ( fqa.isFluctuatingCharge() ) {
758 +      charge += atom->getFlucQPos();
759 +    }
760 +    return charge;
761 +  }
762  
763 +  RealType SelectionEvaluator::getCharge(Atom* atom, int frame) {
764 +    RealType charge = 0.0;
765 +    AtomType* atomType = atom->getAtomType();    
766 +
767 +    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
768 +    if ( fca.isFixedCharge() ) {
769 +      charge = fca.getCharge();
770 +    }
771 +    
772 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
773 +    if ( fqa.isFluctuatingCharge() ) {
774 +      charge += atom->getFlucQPos(frame);
775 +    }
776      return charge;
777    }
778  

Diff Legend

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