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 2055 by gezelter, Fri Feb 6 18:49:27 2015 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 <    }            
63 <
61 >    nObjects.push_back(info->getNGlobalAtoms() + info->getNGlobalRigidBodies());
62 >    nObjects.push_back(info->getNGlobalBonds());
63 >    nObjects.push_back(info->getNGlobalBends());
64 >    nObjects.push_back(info->getNGlobalTorsions());
65 >    nObjects.push_back(info->getNGlobalInversions());
66 >    nObjects.push_back(info->getNGlobalMolecules());    
67 >  }
68 >  
69    bool SelectionEvaluator::loadScript(const std::string& filename,
70                                        const std::string& script) {
71      clearDefinitionsAndLoadPredefined();
# Line 124 | Line 132 | namespace OpenMD {
132      return loadScript(filename, script);
133    }
134    
135 <  void SelectionEvaluator::instructionDispatchLoop(OpenMDBitSet& bs){
135 >  void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs){
136      
137      while ( pc < aatoken.size()) {
138        statement = aatoken[pc++];
# Line 145 | Line 153 | namespace OpenMD {
153  
154    }
155  
156 <  OpenMDBitSet SelectionEvaluator::expression(const std::vector<Token>& code,
156 >  void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs, int frame){
157 >    
158 >    while ( pc < aatoken.size()) {
159 >      statement = aatoken[pc++];
160 >      statementLength = statement.size();
161 >      Token token = statement[0];
162 >      switch (token.tok) {
163 >      case Token::define:
164 >        define();
165 >        break;
166 >      case Token::select:
167 >        select(bs, frame);
168 >        break;
169 >      default:
170 >        unrecognizedCommand(token);
171 >        return;
172 >      }
173 >    }
174 >
175 >  }
176 >
177 >  SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code,
178                                               int pcStart) {
179 <    OpenMDBitSet bs;
180 <    std::stack<OpenMDBitSet> stack;
179 >    SelectionSet bs = createSelectionSets();
180 >    std::stack<SelectionSet> stack;
181 >    vector<int> bsSize = bs.size();
182    
183      for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
184        Token instruction = code[pc];
# Line 160 | Line 190 | namespace OpenMD {
190          break;
191        case Token::all:
192          bs = allInstruction();
193 +        stack.push(bs);
194 +        break;
195 +      case Token::none:
196 +        bs = createSelectionSets();
197 +        stack.push(bs);
198 +        break;
199 +      case Token::opOr:
200 +          bs= stack.top();
201 +          stack.pop();
202 +          stack.top() |= bs;
203 +        break;
204 +      case Token::opAnd:
205 +          bs = stack.top();
206 +          stack.pop();
207 +          stack.top() &= bs;
208 +        break;
209 +      case Token::opNot:
210 +          stack.top().flip();
211 +        break;
212 +      case Token::within:
213 +        withinInstruction(instruction, stack.top());
214 +        break;
215 +      case Token::hull:
216 +          stack.push(hull());
217 +        break;
218 +        //case Token::selected:
219 +        //  stack.push(getSelectionSet());
220 +        //  break;
221 +      case Token::name:
222 +          stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
223 +        break;
224 +      case Token::index:
225 +          stack.push(indexInstruction(instruction.value));
226 +        break;
227 +      case Token::identifier:
228 +          stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
229 +        break;
230 +      case Token::opLT:
231 +      case Token::opLE:
232 +      case Token::opGE:
233 +      case Token::opGT:
234 +      case Token::opEQ:
235 +      case Token::opNE:
236 +        stack.push(comparatorInstruction(instruction));
237 +        break;
238 +      default:
239 +        unrecognizedExpression();
240 +      }
241 +    }
242 +    if (stack.size() != 1)
243 +      evalError("atom expression compiler error - stack over/underflow");
244 +          
245 +    return stack.top();
246 +  }
247 +
248 +
249 +  SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code,
250 +                                              int pcStart, int frame) {
251 +    SelectionSet bs = createSelectionSets();
252 +    std::stack<SelectionSet> stack;
253 +  
254 +    for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
255 +      Token instruction = code[pc];
256 +
257 +      switch (instruction.tok) {
258 +      case Token::expressionBegin:
259 +        break;
260 +      case Token::expressionEnd:
261 +        break;
262 +      case Token::all:
263 +        bs = allInstruction();
264          stack.push(bs);            
265          break;
266        case Token::none:
267 <        bs = OpenMDBitSet(nStuntDouble);
267 >        bs = SelectionSet(nObjects);
268          stack.push(bs);            
269          break;
270        case Token::opOr:
# Line 180 | Line 281 | namespace OpenMD {
281          stack.top().flip();
282          break;
283        case Token::within:
284 <        withinInstruction(instruction, stack.top());
284 >        withinInstruction(instruction, stack.top(), frame);
285          break;
286        case Token::hull:
287 <        stack.push(hull());
287 >        stack.push(hull(frame));
288          break;
289          //case Token::selected:
290          //  stack.push(getSelectionSet());
# Line 203 | Line 304 | namespace OpenMD {
304        case Token::opGT:
305        case Token::opEQ:
306        case Token::opNE:
307 <        stack.push(comparatorInstruction(instruction));
307 >        stack.push(comparatorInstruction(instruction, frame));
308          break;
309        default:
310          unrecognizedExpression();
# Line 217 | Line 318 | namespace OpenMD {
318  
319  
320  
321 <  OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
321 >  SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
322      int comparator = instruction.tok;
323      int property = instruction.intValue;
324      float comparisonValue = boost::any_cast<float>(instruction.value);
325 <    OpenMDBitSet bs(nStuntDouble);
325 >    SelectionSet bs = createSelectionSets();
326      bs.clearAll();
327      
328      SimInfo::MoleculeIterator mi;
# Line 234 | Line 335 | namespace OpenMD {
335      for (mol = info->beginMolecule(mi); mol != NULL;
336           mol = info->nextMolecule(mi)) {
337  
338 +      compareProperty(mol, bs, property, comparator, comparisonValue);
339 +
340        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
341          compareProperty(atom, bs, property, comparator, comparisonValue);
342        }
# Line 247 | Line 350 | namespace OpenMD {
350      return bs;
351    }
352  
353 <  void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs,
353 >  SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction, int frame) {
354 >    int comparator = instruction.tok;
355 >    int property = instruction.intValue;
356 >    float comparisonValue = boost::any_cast<float>(instruction.value);
357 >    SelectionSet bs = createSelectionSets();
358 >    bs.clearAll();
359 >    
360 >    SimInfo::MoleculeIterator mi;
361 >    Molecule* mol;
362 >    Molecule::AtomIterator ai;
363 >    Atom* atom;
364 >    Molecule::RigidBodyIterator rbIter;
365 >    RigidBody* rb;
366 >
367 >    for (mol = info->beginMolecule(mi); mol != NULL;
368 >         mol = info->nextMolecule(mi)) {
369 >
370 >      compareProperty(mol, bs, property, comparator, comparisonValue, frame);
371 >
372 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
373 >        compareProperty(atom, bs, property, comparator, comparisonValue, frame);
374 >      }
375 >    
376 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
377 >           rb = mol->nextRigidBody(rbIter)) {
378 >        compareProperty(rb, bs, property, comparator, comparisonValue, frame);
379 >      }
380 >    }
381 >
382 >    return bs;
383 >  }
384 >
385 >  void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs,
386                                             int property, int comparator,
387                                             float comparisonValue) {
388      RealType propertyValue = 0.0;
389 +    Vector3d pos;
390 +
391      switch (property) {
392      case Token::mass:
393        propertyValue = sd->getMass();
# Line 277 | Line 414 | namespace OpenMD {
414      case Token::z:
415        propertyValue = sd->getPos().z();
416        break;
417 +    case Token::wrappedX:    
418 +      pos = sd->getPos();
419 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
420 +      propertyValue = pos.x();
421 +      break;
422 +    case Token::wrappedY:
423 +      pos = sd->getPos();
424 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
425 +      propertyValue = pos.y();
426 +      break;
427 +    case Token::wrappedZ:
428 +      pos = sd->getPos();
429 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
430 +      propertyValue = pos.z();
431 +      break;
432      case Token::r:
433        propertyValue = sd->getPos().length();
434        break;
# Line 305 | Line 457 | namespace OpenMD {
457        match = propertyValue != comparisonValue;
458        break;
459      }
460 <    if (match)
461 <      bs.setBitOn(sd->getGlobalIndex());
460 >
461 >    if (match)
462 >      bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex());
463 >
464      
465    }
466  
467 <  void SelectionEvaluator::withinInstruction(const Token& instruction,
468 <                                             OpenMDBitSet& bs){
467 >  void SelectionEvaluator::compareProperty(Molecule* mol, SelectionSet& bs,
468 >                                           int property, int comparator,
469 >                                           float comparisonValue) {
470 >    RealType propertyValue = 0.0;
471 >    Vector3d pos;
472 >
473 >    switch (property) {
474 >    case Token::mass:
475 >      propertyValue = mol->getMass();
476 >      break;
477 >    case Token::charge:
478 >      {
479 >        Molecule::AtomIterator ai;
480 >        Atom* atom;
481 >        for (atom = mol->beginAtom(ai); atom != NULL;
482 >             atom = mol->nextAtom(ai)) {
483 >          propertyValue+=  getCharge(atom);
484 >        }
485 >      }
486 >      break;
487 >    case Token::x:
488 >      propertyValue = mol->getCom().x();
489 >      break;
490 >    case Token::y:
491 >      propertyValue = mol->getCom().y();
492 >      break;
493 >    case Token::z:
494 >      propertyValue = mol->getCom().z();
495 >      break;
496 >    case Token::wrappedX:    
497 >      pos = mol->getCom();
498 >      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
499 >      propertyValue = pos.x();
500 >      break;
501 >    case Token::wrappedY:
502 >      pos = mol->getCom();
503 >      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
504 >      propertyValue = pos.y();
505 >      break;
506 >    case Token::wrappedZ:
507 >      pos = mol->getCom();
508 >      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
509 >      propertyValue = pos.z();
510 >      break;
511 >    case Token::r:
512 >      propertyValue = mol->getCom().length();
513 >      break;
514 >    default:
515 >      unrecognizedMoleculeProperty(property);
516 >    }
517      
518 +    bool match = false;
519 +    switch (comparator) {
520 +    case Token::opLT:
521 +      match = propertyValue < comparisonValue;
522 +      break;
523 +    case Token::opLE:
524 +      match = propertyValue <= comparisonValue;
525 +      break;
526 +    case Token::opGE:
527 +      match = propertyValue >= comparisonValue;
528 +      break;
529 +    case Token::opGT:
530 +      match = propertyValue > comparisonValue;
531 +      break;
532 +    case Token::opEQ:
533 +      match = propertyValue == comparisonValue;
534 +      break;
535 +    case Token::opNE:
536 +      match = propertyValue != comparisonValue;
537 +      break;
538 +    }
539 +    
540 +    if (match)
541 +      bs.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());    
542 +  }
543 +
544 +  void SelectionEvaluator::compareProperty(Molecule* mol, SelectionSet& bs,
545 +                                           int property, int comparator,
546 +                                           float comparisonValue, int frame) {
547 +    RealType propertyValue = 0.0;
548 +    Vector3d pos;
549 +    switch (property) {
550 +    case Token::mass:
551 +      propertyValue = mol->getMass();
552 +      break;
553 +    case Token::charge:
554 +      {
555 +        Molecule::AtomIterator ai;
556 +        Atom* atom;
557 +        for (atom = mol->beginAtom(ai); atom != NULL;
558 +             atom = mol->nextAtom(ai)) {
559 +          propertyValue+=  getCharge(atom,frame);
560 +        }
561 +      }      
562 +      break;
563 +    case Token::x:
564 +      propertyValue = mol->getCom(frame).x();
565 +      break;
566 +    case Token::y:
567 +      propertyValue = mol->getCom(frame).y();
568 +      break;
569 +    case Token::z:
570 +      propertyValue = mol->getCom(frame).z();
571 +      break;
572 +    case Token::wrappedX:    
573 +      pos = mol->getCom(frame);
574 +      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
575 +      propertyValue = pos.x();
576 +      break;
577 +    case Token::wrappedY:
578 +      pos = mol->getCom(frame);
579 +      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
580 +      propertyValue = pos.y();
581 +      break;
582 +    case Token::wrappedZ:
583 +      pos = mol->getCom(frame);
584 +      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
585 +      propertyValue = pos.z();
586 +      break;
587 +
588 +    case Token::r:
589 +      propertyValue = mol->getCom(frame).length();
590 +      break;
591 +    default:
592 +      unrecognizedMoleculeProperty(property);
593 +    }
594 +        
595 +    bool match = false;
596 +    switch (comparator) {
597 +    case Token::opLT:
598 +      match = propertyValue < comparisonValue;
599 +      break;
600 +    case Token::opLE:
601 +      match = propertyValue <= comparisonValue;
602 +      break;
603 +    case Token::opGE:
604 +      match = propertyValue >= comparisonValue;
605 +      break;
606 +    case Token::opGT:
607 +      match = propertyValue > comparisonValue;
608 +      break;
609 +    case Token::opEQ:
610 +      match = propertyValue == comparisonValue;
611 +      break;
612 +    case Token::opNE:
613 +      match = propertyValue != comparisonValue;
614 +      break;
615 +    }
616 +    if (match)
617 +      bs.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
618 +    
619 +    
620 +  }
621 +  void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs,
622 +                                           int property, int comparator,
623 +                                           float comparisonValue, int frame) {
624 +    RealType propertyValue = 0.0;
625 +    Vector3d pos;
626 +    switch (property) {
627 +    case Token::mass:
628 +      propertyValue = sd->getMass();
629 +      break;
630 +    case Token::charge:
631 +      if (sd->isAtom()){
632 +        Atom* atom = static_cast<Atom*>(sd);
633 +        propertyValue = getCharge(atom,frame);
634 +      } else if (sd->isRigidBody()) {
635 +        RigidBody* rb = static_cast<RigidBody*>(sd);
636 +        RigidBody::AtomIterator ai;
637 +        Atom* atom;
638 +        for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
639 +          propertyValue+=  getCharge(atom,frame);
640 +        }
641 +      }
642 +      break;
643 +    case Token::x:
644 +      propertyValue = sd->getPos(frame).x();
645 +      break;
646 +    case Token::y:
647 +      propertyValue = sd->getPos(frame).y();
648 +      break;
649 +    case Token::z:
650 +      propertyValue = sd->getPos(frame).z();
651 +      break;
652 +    case Token::wrappedX:    
653 +      pos = sd->getPos(frame);
654 +      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
655 +      propertyValue = pos.x();
656 +      break;
657 +    case Token::wrappedY:
658 +      pos = sd->getPos(frame);
659 +      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
660 +      propertyValue = pos.y();
661 +      break;
662 +    case Token::wrappedZ:
663 +      pos = sd->getPos(frame);
664 +      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
665 +      propertyValue = pos.z();
666 +      break;
667 +
668 +    case Token::r:
669 +      propertyValue = sd->getPos(frame).length();
670 +      break;
671 +    default:
672 +      unrecognizedAtomProperty(property);
673 +    }
674 +        
675 +    bool match = false;
676 +    switch (comparator) {
677 +    case Token::opLT:
678 +      match = propertyValue < comparisonValue;
679 +      break;
680 +    case Token::opLE:
681 +      match = propertyValue <= comparisonValue;
682 +      break;
683 +    case Token::opGE:
684 +      match = propertyValue >= comparisonValue;
685 +      break;
686 +    case Token::opGT:
687 +      match = propertyValue > comparisonValue;
688 +      break;
689 +    case Token::opEQ:
690 +      match = propertyValue == comparisonValue;
691 +      break;
692 +    case Token::opNE:
693 +      match = propertyValue != comparisonValue;
694 +      break;
695 +    }
696 +    if (match)
697 +      bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex());
698 +
699 +    
700 +  }
701 +
702 +  
703 +  void SelectionEvaluator::withinInstruction(const Token& instruction,
704 +                                             SelectionSet& bs){
705 +    
706      boost::any withinSpec = instruction.value;
707      float distance;
708      if (withinSpec.type() == typeid(float)){
# Line 326 | Line 716 | namespace OpenMD {
716      
717      bs = distanceFinder.find(bs, distance);            
718    }
719 +
720 +  void SelectionEvaluator::withinInstruction(const Token& instruction,
721 +                                             SelectionSet& bs, int frame){
722 +    
723 +    boost::any withinSpec = instruction.value;
724 +    float distance;
725 +    if (withinSpec.type() == typeid(float)){
726 +      distance = boost::any_cast<float>(withinSpec);
727 +    } else if (withinSpec.type() == typeid(int)) {
728 +      distance = boost::any_cast<int>(withinSpec);    
729 +    } else {
730 +      evalError("casting error in withinInstruction");
731 +      bs.clearAll();
732 +    }
733 +    
734 +    bs = distanceFinder.find(bs, distance, frame);            
735 +  }
736    
737    void SelectionEvaluator::define() {
738      assert(statement.size() >= 3);
# Line 368 | Line 775 | namespace OpenMD {
775      }
776    }
777  
778 <  void SelectionEvaluator::select(OpenMDBitSet& bs){
778 >  void SelectionEvaluator::select(SelectionSet& bs){
779      bs = expression(statement, 1);
780    }
781 <  
782 <  OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){
781 >
782 >  void SelectionEvaluator::select(SelectionSet& bs, int frame){
783 >    bs = expression(statement, 1, frame);
784 >  }
785 >  
786 >  SelectionSet SelectionEvaluator::lookupValue(const std::string& variable){
787      
788 <    OpenMDBitSet bs(nStuntDouble);
788 >    SelectionSet bs = createSelectionSets();
789      std::map<std::string, boost::any>::iterator i = variables.find(variable);
790      
791      if (i != variables.end()) {
792 <      if (i->second.type() == typeid(OpenMDBitSet)) {
793 <        return boost::any_cast<OpenMDBitSet>(i->second);
792 >      if (i->second.type() == typeid(SelectionSet)) {
793 >        return boost::any_cast<SelectionSet>(i->second);
794        } else if (i->second.type() ==  typeid(std::vector<Token>)){
795          bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
796          i->second =  bs; /**@todo fixme */
# Line 392 | Line 803 | namespace OpenMD {
803      return bs;
804    }
805    
806 <  OpenMDBitSet SelectionEvaluator::nameInstruction(const std::string& name){    
806 >  SelectionSet SelectionEvaluator::nameInstruction(const std::string& name){    
807      return nameFinder.match(name);    
808    }    
809  
# Line 413 | Line 824 | namespace OpenMD {
824      //predefine();
825    }
826  
827 <  OpenMDBitSet SelectionEvaluator::evaluate() {
828 <    OpenMDBitSet bs(nStuntDouble);
827 >  SelectionSet SelectionEvaluator::createSelectionSets() {
828 >    SelectionSet ss(nObjects);
829 >    return ss;
830 >  }
831 >
832 >  SelectionSet SelectionEvaluator::evaluate() {
833 >    SelectionSet bs = createSelectionSets();
834      if (isLoaded_) {
835        pc = 0;
836        instructionDispatchLoop(bs);
837      }
838 +    return bs;
839 +  }
840  
841 +  SelectionSet SelectionEvaluator::evaluate(int frame) {
842 +    SelectionSet bs = createSelectionSets();
843 +    if (isLoaded_) {
844 +      pc = 0;
845 +      instructionDispatchLoop(bs, frame);
846 +    }
847      return bs;
848    }
849  
850 <  OpenMDBitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
851 <    OpenMDBitSet bs(nStuntDouble);
850 >  SelectionSet SelectionEvaluator::indexInstruction(const boost::any& value) {
851 >    SelectionSet bs = createSelectionSets();
852  
853      if (value.type() == typeid(int)) {
854        int index = boost::any_cast<int>(value);
855 <      if (index < 0 || index >= bs.size()) {
855 >      if (index < 0 || index >= bs.bitsets_[STUNTDOUBLE].size()) {
856          invalidIndex(index);
857        } else {
858          bs = indexFinder.find(index);
# Line 436 | Line 860 | namespace OpenMD {
860      } else if (value.type() == typeid(std::pair<int, int>)) {
861        std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
862        assert(indexRange.first <= indexRange.second);
863 <      if (indexRange.first < 0 || indexRange.second >= bs.size()) {
863 >      if (indexRange.first < 0 ||
864 >          indexRange.second >= bs.bitsets_[STUNTDOUBLE].size()) {
865          invalidIndexRange(indexRange);
866        }else {
867          bs = indexFinder.find(indexRange.first, indexRange.second);
# Line 446 | Line 871 | namespace OpenMD {
871      return bs;
872    }
873  
874 <  OpenMDBitSet SelectionEvaluator::allInstruction() {
875 <    OpenMDBitSet bs(nStuntDouble);
874 >  SelectionSet SelectionEvaluator::allInstruction() {
875 >    SelectionSet ss = createSelectionSets();
876  
877      SimInfo::MoleculeIterator mi;
453 –    Molecule* mol;
878      Molecule::AtomIterator ai;
455 –    Atom* atom;
879      Molecule::RigidBodyIterator rbIter;
880 +    Molecule::BondIterator bondIter;
881 +    Molecule::BendIterator bendIter;
882 +    Molecule::TorsionIterator torsionIter;
883 +    Molecule::InversionIterator inversionIter;
884 +
885 +    Molecule* mol;
886 +    Atom* atom;
887      RigidBody* rb;
888 +    Bond* bond;
889 +    Bend* bend;
890 +    Torsion* torsion;
891 +    Inversion* inversion;    
892  
893      // Doing the loop insures that we're actually on this processor.
894  
895      for (mol = info->beginMolecule(mi); mol != NULL;
896           mol = info->nextMolecule(mi)) {
463 –
897        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
898 <        bs.setBitOn(atom->getGlobalIndex());
899 <      }
467 <    
898 >        ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
899 >      }    
900        for (rb = mol->beginRigidBody(rbIter); rb != NULL;
901             rb = mol->nextRigidBody(rbIter)) {
902 <        bs.setBitOn(rb->getGlobalIndex());
902 >        ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
903        }
904 +      for (bond = mol->beginBond(bondIter); bond != NULL;
905 +           bond = mol->nextBond(bondIter)) {
906 +        ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
907 +      }  
908 +      for (bend = mol->beginBend(bendIter); bend != NULL;
909 +           bend = mol->nextBend(bendIter)) {
910 +        ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
911 +      }  
912 +      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
913 +           torsion = mol->nextTorsion(torsionIter)) {
914 +        ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
915 +      }  
916 +      for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
917 +           inversion = mol->nextInversion(inversionIter)) {
918 +        ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
919 +      }
920 +      ss.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
921      }
922  
923 <    return bs;
923 >    return ss;
924    }
925  
926 <  OpenMDBitSet SelectionEvaluator::hull() {
927 <    OpenMDBitSet bs(nStuntDouble);
926 >  SelectionSet SelectionEvaluator::hull() {
927 >    SelectionSet bs = createSelectionSets();
928      
929      bs = hullFinder.findHull();
930 +    surfaceArea_ = hullFinder.getSurfaceArea();
931 +    hasSurfaceArea_ = true;
932 +    return bs;
933 +  }
934 +
935 +
936 +  SelectionSet SelectionEvaluator::hull(int frame) {
937 +    SelectionSet bs = createSelectionSets();
938      
939 +    bs = hullFinder.findHull(frame);
940 +
941      return bs;
942    }
943  
944    RealType SelectionEvaluator::getCharge(Atom* atom) {
945 <    RealType charge =0.0;
945 >    RealType charge = 0.0;
946      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);
947  
948 <        if (doubleData != NULL) {
949 <          charge = doubleData->getData();
950 <
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 <      }
948 >    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
949 >    if ( fca.isFixedCharge() ) {
950 >      charge = fca.getCharge();
951      }
952 +    
953 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
954 +    if ( fqa.isFluctuatingCharge() ) {
955 +      charge += atom->getFlucQPos();
956 +    }
957 +    return charge;
958 +  }
959  
960 +  RealType SelectionEvaluator::getCharge(Atom* atom, int frame) {
961 +    RealType charge = 0.0;
962 +    AtomType* atomType = atom->getAtomType();    
963 +
964 +    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
965 +    if ( fca.isFixedCharge() ) {
966 +      charge = fca.getCharge();
967 +    }
968 +    
969 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
970 +    if ( fqa.isFluctuatingCharge() ) {
971 +      charge += atom->getFlucQPos(frame);
972 +    }
973      return charge;
974    }
975  

Diff Legend

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