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 413 by tim, Wed Mar 9 17:30:29 2005 UTC vs.
Revision 1412 by gezelter, Mon Mar 22 18:45:39 2010 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
28   * arising out of the use of or inability to use software, even if the
29   * University of Notre Dame has been advised of the possibility of
30   * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
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).                        
40   */
41  
42   #include <stack>
# Line 45 | Line 45
45   #include "primitives/DirectionalAtom.hpp"
46   #include "primitives/RigidBody.hpp"
47   #include "primitives/Molecule.hpp"
48 + #include "io/basic_ifstrstream.hpp"
49  
50 < namespace oopse {
50 > namespace OpenMD {
51  
52  
53 < SelectionEvaluator::SelectionEvaluator(SimInfo* si)
54 <    : info(si), nameFinder(info), distanceFinder(info), indexFinder(info), isLoaded_(false){
55 <    
56 <    nStuntDouble = info->getNGlobalAtoms() + info->getNRigidBodies();
57 < }            
53 >  SelectionEvaluator::SelectionEvaluator(SimInfo* si)
54 >    : info(si), nameFinder(info), distanceFinder(info), hullFinder(info),
55 >      indexFinder(info),
56 >      isLoaded_(false){    
57 >      nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies();
58 >    }            
59  
60 < bool SelectionEvaluator::loadScript(const std::string& filename, const std::string& script) {
60 >  bool SelectionEvaluator::loadScript(const std::string& filename,
61 >                                      const std::string& script) {
62      clearDefinitionsAndLoadPredefined();
63      this->filename = filename;
64      this->script = script;
65      if (! compiler.compile(filename, script)) {
66 <        error = true;
67 <        errorMessage = compiler.getErrorMessage();
68 <        std::cerr << "SelectionCompiler Error: " << errorMessage << std::endl;
69 <        return false;
66 >      error = true;
67 >      errorMessage = compiler.getErrorMessage();
68 >
69 >      sprintf( painCave.errMsg,
70 >               "SelectionCompiler Error: %s\n", errorMessage.c_str());
71 >      painCave.severity = OPENMD_ERROR;
72 >      painCave.isFatal = 1;
73 >      simError();
74 >      return false;
75      }
76  
77      pc = 0;
# Line 75 | Line 83 | bool SelectionEvaluator::loadScript(const std::string&
83  
84      isDynamic_ = false;
85      for (i = aatoken.begin(); i != aatoken.end(); ++i) {
86 <        if (containDynamicToken(*i)) {
87 <            isDynamic_ = true;
88 <            break;
89 <        }
86 >      if (containDynamicToken(*i)) {
87 >        isDynamic_ = true;
88 >        break;
89 >      }
90      }
91  
92      isLoaded_ = true;
93      return true;
94 < }
94 >  }
95  
96 < void SelectionEvaluator::clearState() {
89 <    //for (int i = scriptLevelMax; --i >= 0; )
90 <    //    stack[i].clear();
91 <    //scriptLevel = 0;
96 >  void SelectionEvaluator::clearState() {
97      error = false;
98      errorMessage = "";
99 < }
99 >  }
100  
101 < bool SelectionEvaluator::loadScriptString(const std::string& script) {
101 >  bool SelectionEvaluator::loadScriptString(const std::string& script) {
102      clearState();
103      return loadScript("", script);
104 < }
104 >  }
105  
106 < bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
106 >  bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
107      clearState();
108      return loadScriptFileInternal(filename);
109 < }
109 >  }
110  
111 < bool SelectionEvaluator::loadScriptFileInternal(const  std::string & filename) {
112 <  std::ifstream ifs(filename.c_str());
111 >  bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) {
112 >    ifstrstream ifs(filename.c_str());
113      if (!ifs.is_open()) {
114 <        return false;
114 >      return false;
115      }
116 <
116 >    
117      const int bufferSize = 65535;
118      char buffer[bufferSize];
119      std::string script;
120      while(ifs.getline(buffer, bufferSize)) {
121 <        script += buffer;
121 >      script += buffer;
122      }
123      return loadScript(filename, script);
124 < }
125 <
126 < void SelectionEvaluator::instructionDispatchLoop(BitSet& bs){
124 >  }
125 >  
126 >  void SelectionEvaluator::instructionDispatchLoop(OpenMDBitSet& bs){
127      
128      while ( pc < aatoken.size()) {
129 <        statement = aatoken[pc++];
130 <        statementLength = statement.size();
131 <        Token token = statement[0];
132 <        switch (token.tok) {
133 <            case Token::define:
134 <                define();
135 <            break;
136 <            case Token::select:
137 <                select(bs);
138 <            break;
139 <            default:
140 <                unrecognizedCommand(token);
141 <            return;
142 <        }
129 >      statement = aatoken[pc++];
130 >      statementLength = statement.size();
131 >      Token token = statement[0];
132 >      switch (token.tok) {
133 >      case Token::define:
134 >        define();
135 >        break;
136 >      case Token::select:
137 >        select(bs);
138 >        break;
139 >      default:
140 >        unrecognizedCommand(token);
141 >        return;
142 >      }
143      }
144  
145 < }
145 >  }
146  
147 <  BitSet SelectionEvaluator::expression(const std::vector<Token>& code, int pcStart) {
148 <    BitSet bs;
149 <    std::stack<BitSet> stack;
150 <    
147 >  OpenMDBitSet SelectionEvaluator::expression(const std::vector<Token>& code,
148 >                                             int pcStart) {
149 >    OpenMDBitSet bs;
150 >    std::stack<OpenMDBitSet> stack;
151 >  
152      for (int pc = pcStart; pc < code.size(); ++pc) {
153        Token instruction = code[pc];
154  
# Line 152 | Line 158 | void SelectionEvaluator::instructionDispatchLoop(BitSe
158        case Token::expressionEnd:
159          break;
160        case Token::all:
161 <        bs = BitSet(nStuntDouble);
161 >        bs = OpenMDBitSet(nStuntDouble);
162          bs.setAll();
163          stack.push(bs);            
164          break;
165        case Token::none:
166 <        bs = BitSet(nStuntDouble);
166 >        bs = OpenMDBitSet(nStuntDouble);
167          stack.push(bs);            
168          break;
169        case Token::opOr:
# Line 174 | Line 180 | void SelectionEvaluator::instructionDispatchLoop(BitSe
180          stack.top().flip();
181          break;
182        case Token::within:
177 –
183          withinInstruction(instruction, stack.top());
184          break;
185 <      //case Token::selected:
186 <      //  stack.push(getSelectionSet());
187 <      //  break;
185 >      case Token::hull:
186 >        stack.push(hull());
187 >        break;
188 >        //case Token::selected:
189 >        //  stack.push(getSelectionSet());
190 >        //  break;
191        case Token::name:
192          stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
193          break;
# Line 209 | Line 217 | void SelectionEvaluator::instructionDispatchLoop(BitSe
217  
218  
219  
220 < BitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
220 >  OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
221      int comparator = instruction.tok;
222      int property = instruction.intValue;
223      float comparisonValue = boost::any_cast<float>(instruction.value);
224      float propertyValue;
225 <    BitSet bs(nStuntDouble);
225 >    OpenMDBitSet bs(nStuntDouble);
226      bs.clearAll();
227      
228      SimInfo::MoleculeIterator mi;
# Line 223 | Line 231 | BitSet SelectionEvaluator::comparatorInstruction(const
231      Atom* atom;
232      Molecule::RigidBodyIterator rbIter;
233      RigidBody* rb;
226 –    
227 –    for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
234  
235 <        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
236 <            compareProperty(atom, bs, property, comparator, comparisonValue);
237 <        }
238 <        
239 <        //change the positions of atoms which belong to the rigidbodies
240 <        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
241 <            compareProperty(rb, bs, property, comparator, comparisonValue);
242 <        }        
235 >    for (mol = info->beginMolecule(mi); mol != NULL;
236 >         mol = info->nextMolecule(mi)) {
237 >
238 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
239 >        compareProperty(atom, bs, property, comparator, comparisonValue);
240 >      }
241 >    
242 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
243 >           rb = mol->nextRigidBody(rbIter)) {
244 >        compareProperty(rb, bs, property, comparator, comparisonValue);
245 >      }
246      }
247  
248      return bs;
249 < }
249 >  }
250  
251 < void SelectionEvaluator::compareProperty(StuntDouble* sd, BitSet& bs, int property, int comparator, float comparisonValue) {
252 <        double propertyValue;
253 <        switch (property) {
254 <        case Token::mass:
255 <            propertyValue = sd->getMass();
256 <            break;
257 <        case Token::charge:
258 <            return;
259 <            //break;
260 <        case Token::dipole:
261 <            return;
262 <            //break;
263 <        default:
264 <            unrecognizedAtomProperty(property);
265 <        }
251 >  void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs,
252 >                                           int property, int comparator,
253 >                                           float comparisonValue) {
254 >    RealType propertyValue = 0.0;
255 >    switch (property) {
256 >    case Token::mass:
257 >      propertyValue = sd->getMass();
258 >      break;
259 >    case Token::charge:
260 >      if (sd->isAtom()){
261 >        Atom* atom = static_cast<Atom*>(sd);
262 >        propertyValue = getCharge(atom);
263 >      } else if (sd->isRigidBody()) {
264 >        RigidBody* rb = static_cast<RigidBody*>(sd);
265 >        RigidBody::AtomIterator ai;
266 >        Atom* atom;
267 >        for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
268 >          propertyValue+=  getCharge(atom);
269 >        }
270 >      }
271 >      break;
272 >    case Token::x:
273 >      propertyValue = sd->getPos().x();
274 >      break;
275 >    case Token::y:
276 >      propertyValue = sd->getPos().y();
277 >      break;
278 >    case Token::z:
279 >      propertyValue = sd->getPos().z();
280 >      break;
281 >    default:
282 >      unrecognizedAtomProperty(property);
283 >    }
284          
285 <        bool match = false;
286 <        switch (comparator) {
287 <            case Token::opLT:
288 <                match = propertyValue < comparisonValue;
289 <                break;
290 <            case Token::opLE:
291 <                match = propertyValue <= comparisonValue;
292 <                break;
293 <            case Token::opGE:
294 <                match = propertyValue >= comparisonValue;
295 <                break;
296 <            case Token::opGT:
297 <                match = propertyValue > comparisonValue;
298 <                break;
299 <            case Token::opEQ:
300 <                match = propertyValue == comparisonValue;
301 <                break;
302 <            case Token::opNE:
303 <                match = propertyValue != comparisonValue;
304 <                break;
305 <        }
306 <        if (match)
307 <            bs.setBitOn(sd->getGlobalIndex());
308 <
309 < }
285 >    bool match = false;
286 >    switch (comparator) {
287 >    case Token::opLT:
288 >      match = propertyValue < comparisonValue;
289 >      break;
290 >    case Token::opLE:
291 >      match = propertyValue <= comparisonValue;
292 >      break;
293 >    case Token::opGE:
294 >      match = propertyValue >= comparisonValue;
295 >      break;
296 >    case Token::opGT:
297 >      match = propertyValue > comparisonValue;
298 >      break;
299 >    case Token::opEQ:
300 >      match = propertyValue == comparisonValue;
301 >      break;
302 >    case Token::opNE:
303 >      match = propertyValue != comparisonValue;
304 >      break;
305 >    }
306 >    if (match)
307 >      bs.setBitOn(sd->getGlobalIndex());
308 >    
309 >  }
310  
311 < void SelectionEvaluator::withinInstruction(const Token& instruction, BitSet& bs){
311 >  void SelectionEvaluator::withinInstruction(const Token& instruction,
312 >                                             OpenMDBitSet& bs){
313      
314      boost::any withinSpec = instruction.value;
315      float distance;
316      if (withinSpec.type() == typeid(float)){
317 <        distance = boost::any_cast<float>(withinSpec);
317 >      distance = boost::any_cast<float>(withinSpec);
318      } else if (withinSpec.type() == typeid(int)) {
319 <        distance = boost::any_cast<int>(withinSpec);    
319 >      distance = boost::any_cast<int>(withinSpec);    
320      } else {
321 <        evalError("casting error in withinInstruction");
322 <        bs.clearAll();
321 >      evalError("casting error in withinInstruction");
322 >      bs.clearAll();
323      }
324      
325      bs = distanceFinder.find(bs, distance);            
326 < }
327 <
328 < void SelectionEvaluator::define() {
326 >  }
327 >  
328 >  void SelectionEvaluator::define() {
329      assert(statement.size() >= 3);
330 <
330 >    
331      std::string variable = boost::any_cast<std::string>(statement[1].value);
332 +    
333 +    variables.insert(VariablesType::value_type(variable,
334 +                                               expression(statement, 2)));
335 +  }
336 +  
337  
338 <    variables.insert(VariablesType::value_type(variable, expression(statement, 2)));
339 < }
340 <
308 <
309 < /** @todo */
310 < void SelectionEvaluator::predefine(const std::string& script) {
311 <
338 >  /** @todo */
339 >  void SelectionEvaluator::predefine(const std::string& script) {
340 >    
341      if (compiler.compile("#predefine", script)) {
342 <        std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
343 <        if (aatoken.size() != 1) {
344 <            evalError("predefinition does not have exactly 1 command:"
345 <                + script);
346 <            return;
347 <        }
348 <        std::vector<Token> statement = aatoken[0];
349 <        if (statement.size() > 2) {
350 <            int tok = statement[1].tok;
351 <            if (tok == Token::identifier || (tok & Token::predefinedset) == Token::predefinedset) {
352 <                std::string variable = boost::any_cast<std::string>(statement[1].value);
353 <                variables.insert(VariablesType::value_type(variable, statement));
354 <
355 <            } else {
356 <                evalError("invalid variable name:" + script);
357 <            }
358 <        }else {
359 <            evalError("bad predefinition length:" + script);
360 <        }
361 <
342 >      std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
343 >      if (aatoken.size() != 1) {
344 >        evalError("predefinition does not have exactly 1 command:"
345 >                  + script);
346 >        return;
347 >      }
348 >      std::vector<Token> statement = aatoken[0];
349 >      if (statement.size() > 2) {
350 >        int tok = statement[1].tok;
351 >        if (tok == Token::identifier ||
352 >            (tok & Token::predefinedset) == Token::predefinedset) {
353 >          std::string variable = boost::any_cast<std::string>(statement[1].value);
354 >          variables.insert(VariablesType::value_type(variable, statement));
355 >          
356 >        } else {
357 >          evalError("invalid variable name:" + script);
358 >        }
359 >      }else {
360 >        evalError("bad predefinition length:" + script);
361 >      }      
362          
363      } else {
364 <        evalError("predefined set compile error:" + script +
365 <          "\ncompile error:" + compiler.getErrorMessage());
364 >      evalError("predefined set compile error:" + script +
365 >                "\ncompile error:" + compiler.getErrorMessage());
366      }
367 +  }
368  
369 < }
340 <
341 < void SelectionEvaluator::select(BitSet& bs){
369 >  void SelectionEvaluator::select(OpenMDBitSet& bs){
370      bs = expression(statement, 1);
371 < }
372 <
373 < BitSet SelectionEvaluator::lookupValue(const std::string& variable){
374 <
375 <    BitSet bs(nStuntDouble);
371 >  }
372 >  
373 >  OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){
374 >    
375 >    OpenMDBitSet bs(nStuntDouble);
376      std::map<std::string, boost::any>::iterator i = variables.find(variable);
377      
378      if (i != variables.end()) {
379 <        if (i->second.type() == typeid(BitSet)) {
380 <            return boost::any_cast<BitSet>(i->second);
381 <        } else if (i->second.type() ==  typeid(std::vector<Token>)){
382 <            bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
383 <            i->second =  bs; /**@todo fixme */
384 <            return bs;
385 <        }
379 >      if (i->second.type() == typeid(OpenMDBitSet)) {
380 >        return boost::any_cast<OpenMDBitSet>(i->second);
381 >      } else if (i->second.type() ==  typeid(std::vector<Token>)){
382 >        bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
383 >        i->second =  bs; /**@todo fixme */
384 >        return bs;
385 >      }
386      } else {
387 <        unrecognizedIdentifier(variable);
387 >      unrecognizedIdentifier(variable);
388      }
389 <
389 >    
390      return bs;
391 < }
391 >  }
392 >  
393 >  OpenMDBitSet SelectionEvaluator::nameInstruction(const std::string& name){    
394 >    return nameFinder.match(name);    
395 >  }    
396  
397 < BitSet SelectionEvaluator::nameInstruction(const std::string& name){
366 <    
367 <    return nameFinder.match(name);
368 <
369 < }    
370 <
371 < bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
397 >  bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
398      std::vector<Token>::const_iterator i;
399      for (i = tokens.begin(); i != tokens.end(); ++i) {
400 <        if (i->tok & Token::dynamic) {
401 <            return true;
402 <        }
400 >      if (i->tok & Token::dynamic) {
401 >        return true;
402 >      }
403      }
404 <
404 >    
405      return false;
406 < }    
406 >  }    
407  
408 < void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
408 >  void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
409      variables.clear();
410      //load predefine
411      //predefine();
412 < }
412 >  }
413  
414 < BitSet SelectionEvaluator::evaluate() {
415 <    BitSet bs(nStuntDouble);
414 >  OpenMDBitSet SelectionEvaluator::evaluate() {
415 >    OpenMDBitSet bs(nStuntDouble);
416      if (isLoaded_) {
417 <        pc = 0;
418 <        instructionDispatchLoop(bs);
417 >      pc = 0;
418 >      instructionDispatchLoop(bs);
419      }
420  
421      return bs;
422 < }
422 >  }
423  
424 < BitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
425 <    BitSet bs(nStuntDouble);
424 >  OpenMDBitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
425 >    OpenMDBitSet bs(nStuntDouble);
426  
427      if (value.type() == typeid(int)) {
428 <        int index = boost::any_cast<int>(value);
429 <        if (index < 0 || index >= bs.size()) {
430 <            invalidIndex(index);
431 <        } else {
432 <            indexFinder.find(index);
433 <        }
428 >      int index = boost::any_cast<int>(value);
429 >      if (index < 0 || index >= bs.size()) {
430 >        invalidIndex(index);
431 >      } else {
432 >        bs = indexFinder.find(index);
433 >      }
434      } else if (value.type() == typeid(std::pair<int, int>)) {
435 <        std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
436 <        assert(indexRange.first <= indexRange.second);
437 <        if (indexRange.first < 0 || indexRange.second >= bs.size()) {
438 <            invalidIndexRange(indexRange);
439 <        }else {
440 <            indexFinder.find(indexRange.first, indexRange.second);
441 <        }
435 >      std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
436 >      assert(indexRange.first <= indexRange.second);
437 >      if (indexRange.first < 0 || indexRange.second >= bs.size()) {
438 >        invalidIndexRange(indexRange);
439 >      }else {
440 >        bs = indexFinder.find(indexRange.first, indexRange.second);
441 >      }
442      }
443  
444      return bs;
445 < }
445 >  }
446  
447 +
448 +  OpenMDBitSet SelectionEvaluator::hull() {
449 +    OpenMDBitSet bs(nStuntDouble);
450 +    
451 +    bs = hullFinder.findHull();
452 +    
453 +    return bs;
454 +  }
455 +
456 +
457 +
458 +  RealType SelectionEvaluator::getCharge(Atom* atom) {
459 +    RealType charge =0.0;
460 +    AtomType* atomType = atom->getAtomType();
461 +    if (atomType->isCharge()) {
462 +      GenericData* data = atomType->getPropertyByName("Charge");
463 +      if (data != NULL) {
464 +        DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data);
465 +
466 +        if (doubleData != NULL) {
467 +          charge = doubleData->getData();
468 +
469 +        } else {
470 +          sprintf( painCave.errMsg,
471 +                   "Can not cast GenericData to DoubleGenericData\n");
472 +          painCave.severity = OPENMD_ERROR;
473 +          painCave.isFatal = 1;
474 +          simError();          
475 +        }
476 +      }
477 +    }
478 +
479 +    return charge;
480 +  }
481 +
482   }

Diff Legend

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