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 282 by tim, Thu Feb 3 14:04:59 2005 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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, 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>
44   #include "selection/SelectionEvaluator.hpp"
45 < namespace oopse {
45 > #include "primitives/Atom.hpp"
46 > #include "primitives/DirectionalAtom.hpp"
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 < bool SelectionEvaluator::loadScript(const std::string& filename, const std::string& script) {
54 > namespace OpenMD {
55 >
56 >  SelectionEvaluator::SelectionEvaluator(SimInfo* si)
57 >    : info(si), nameFinder(info), distanceFinder(info), hullFinder(info),
58 >      indexFinder(info), isLoaded_(false), hasSurfaceArea_(false) {
59 >    nObjects.push_back(info->getNGlobalAtoms() + info->getNGlobalRigidBodies());
60 >    nObjects.push_back(info->getNGlobalBonds());
61 >    nObjects.push_back(info->getNGlobalBends());
62 >    nObjects.push_back(info->getNGlobalTorsions());
63 >    nObjects.push_back(info->getNGlobalInversions());
64 >    nObjects.push_back(info->getNGlobalMolecules());    
65 >  }
66 >  
67 >  bool SelectionEvaluator::loadScript(const std::string& filename,
68 >                                      const std::string& script) {
69 >    clearDefinitionsAndLoadPredefined();
70      this->filename = filename;
71      this->script = script;
72      if (! compiler.compile(filename, script)) {
73 <        error = true;
74 <        errorMessage = compiler.getErrorMessage();
75 <        return false;
73 >      error = true;
74 >      errorMessage = compiler.getErrorMessage();
75 >
76 >      sprintf( painCave.errMsg,
77 >               "SelectionCompiler Error: %s\n", errorMessage.c_str());
78 >      painCave.severity = OPENMD_ERROR;
79 >      painCave.isFatal = 1;
80 >      simError();
81 >      return false;
82      }
83  
84      pc = 0;
85      aatoken = compiler.getAatokenCompiled();
86      linenumbers = compiler.getLineNumbers();
87      lineIndices = compiler.getLineIndices();
88 +
89 +    std::vector<std::vector<Token> >::const_iterator i;  
90 +
91 +    isDynamic_ = false;
92 +    for (i = aatoken.begin(); i != aatoken.end(); ++i) {
93 +      if (containDynamicToken(*i)) {
94 +        isDynamic_ = true;
95 +        break;
96 +      }
97 +    }
98 +
99 +    isLoaded_ = true;
100      return true;
101 < }
101 >  }
102  
103 < void SelectionEvaluator::clearState() {
63 <    for (int i = scriptLevelMax; --i >= 0; )
64 <        stack[i].clear();
65 <    scriptLevel = 0;
103 >  void SelectionEvaluator::clearState() {
104      error = false;
105      errorMessage = "";
106 < }
106 >  }
107  
108 < bool SelectionEvaluator::loadScriptString(const std::string& script) {
108 >  bool SelectionEvaluator::loadScriptString(const std::string& script) {
109      clearState();
110      return loadScript("", script);
111 < }
111 >  }
112  
113 < bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
113 >  bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
114      clearState();
115      return loadScriptFileInternal(filename);
116 < }
116 >  }
117  
118 < bool SelectionEvaluator::loadScriptFileInternal(const  string & filename) {
118 >  bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) {
119 >    ifstrstream ifs(filename.c_str());
120 >    if (!ifs.is_open()) {
121 >      return false;
122 >    }
123 >    
124 >    const int bufferSize = 65535;
125 >    char buffer[bufferSize];
126 >    std::string script;
127 >    while(ifs.getline(buffer, bufferSize)) {
128 >      script += buffer;
129 >    }
130 >    return loadScript(filename, script);
131 >  }
132 >  
133 >  void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs){
134 >    
135 >    while ( pc < aatoken.size()) {
136 >      statement = aatoken[pc++];
137 >      statementLength = statement.size();
138 >      Token token = statement[0];
139 >      switch (token.tok) {
140 >      case Token::define:
141 >        define();
142 >        break;
143 >      case Token::select:
144 >        select(bs);
145 >        break;
146 >      default:
147 >        unrecognizedCommand(token);
148 >        return;
149 >      }
150 >    }
151  
152 < }
152 >  }
153  
154 < void SelectionEvaluator::instructionDispatchLoop(){
155 <
154 >  void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs, int frame){
155 >    
156      while ( pc < aatoken.size()) {
157 <        statement = aatoken[pc++];
158 <        statementLength = statement.size();
159 <        Token token = statement[0];
160 <        switch (token.tok) {
161 <            case Token::define:
162 <                define();
163 <            break;
164 <            case Token::select:
165 <                select();
166 <            break;
167 <            default:
168 <                unrecognizedCommand(token);
169 <            return;
170 <        }
157 >      statement = aatoken[pc++];
158 >      statementLength = statement.size();
159 >      Token token = statement[0];
160 >      switch (token.tok) {
161 >      case Token::define:
162 >        define();
163 >        break;
164 >      case Token::select:
165 >        select(bs, frame);
166 >        break;
167 >      default:
168 >        unrecognizedCommand(token);
169 >        return;
170 >      }
171      }
102 }
172  
173 <  BitSet SelectionEvaluator::expression(std::vector<Token>& code, int pcStart) {
105 <    int numberOfAtoms = viewer.getAtomCount();
106 <    BitSet bs;
107 <    BitSet[] stack = new BitSet[10];
108 <    int sp = 0;
173 >  }
174  
175 <    for (int pc = pcStart; ; ++pc) {
175 >  SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code,
176 >                                             int pcStart) {
177 >    SelectionSet bs = createSelectionSets();
178 >    std::stack<SelectionSet> stack;
179 >    vector<int> bsSize = bs.size();
180 >  
181 >    for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
182        Token instruction = code[pc];
183  
184        switch (instruction.tok) {
# Line 116 | Line 187 | void SelectionEvaluator::instructionDispatchLoop(){
187        case Token::expressionEnd:
188          break;
189        case Token::all:
190 <        bs = stack[sp++] = new BitSet(numberOfAtoms);
191 <        for (int i = numberOfAtoms; --i >= 0; )
121 <          bs.set(i);
190 >        bs = allInstruction();
191 >        stack.push(bs);
192          break;
193        case Token::none:
194 <        stack[sp++] = new BitSet();
194 >        bs = createSelectionSets();
195 >        stack.push(bs);
196          break;
197        case Token::opOr:
198 <        bs = stack[--sp];
199 <        stack[sp-1].or(bs);
198 >          bs= stack.top();
199 >          stack.pop();
200 >          stack.top() |= bs;
201          break;
202        case Token::opAnd:
203 <        bs = stack[--sp];
204 <        stack[sp-1].and(bs);
203 >          bs = stack.top();
204 >          stack.pop();
205 >          stack.top() &= bs;
206          break;
207        case Token::opNot:
208 <        bs = stack[sp - 1];
136 <        notSet(bs);
208 >          stack.top().flip();
209          break;
210        case Token::within:
211 <        bs = stack[sp - 1];
140 <        stack[sp - 1] = new BitSet();
141 <        withinInstruction(instruction, bs, stack[sp - 1]);
211 >        withinInstruction(instruction, stack.top());
212          break;
213 <      case Token::selected:
214 <        stack[sp++] = copyBitSet(viewer.getSelectionSet());
213 >      case Token::hull:
214 >          stack.push(hull());
215          break;
216 +        //case Token::selected:
217 +        //  stack.push(getSelectionSet());
218 +        //  break;
219        case Token::name:
220 <
220 >          stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
221          break;
222 <      case  Token::index:
223 <        
222 >      case Token::index:
223 >          stack.push(indexInstruction(instruction.value));
224          break;
152      case Token::molname:  
153
154        break;
155      case Token::molindex:
156        break;
225        case Token::identifier:
226 <        stack[sp++] = lookupIdentifierValue((std::string)instruction.value);
226 >          stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
227          break;
228        case Token::opLT:
229        case Token::opLE:
# Line 163 | Line 231 | void SelectionEvaluator::instructionDispatchLoop(){
231        case Token::opGT:
232        case Token::opEQ:
233        case Token::opNE:
234 <        bs = stack[sp++] = new BitSet();
167 <        comparatorInstruction(instruction, bs);
234 >        stack.push(comparatorInstruction(instruction));
235          break;
236        default:
237          unrecognizedExpression();
238        }
239      }
240 <    if (sp != 1)
240 >    if (stack.size() != 1)
241        evalError("atom expression compiler error - stack over/underflow");
242 <    return stack[0];
242 >          
243 >    return stack.top();
244    }
245  
246  
247 +  SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code,
248 +                                              int pcStart, int frame) {
249 +    SelectionSet bs = createSelectionSets();
250 +    std::stack<SelectionSet> stack;
251 +  
252 +    for (unsigned int pc = pcStart; pc < code.size(); ++pc) {
253 +      Token instruction = code[pc];
254  
255 <  void SelectionEvaluator::comparatorInstruction(Token instruction, BitSet bs) {
256 <    int comparator = instruction.tok;
182 <    int property = instruction.intValue;
183 <    float propertyValue = 0; // just for temperature
184 <    int comparisonValue = ((Integer)instruction.value).intValue();
185 <    int numberOfAtoms = viewer.getAtomCount();
186 <    Frame frame = viewer.getFrame();
187 <    for (int i = 0; i < numberOfAtoms; ++i) {
188 <      Atom atom = frame.getAtomAt(i);
189 <      switch (property) {
190 <      case Token::mass:
191 <        //propertyValue = atom.getAtomNumber();
255 >      switch (instruction.tok) {
256 >      case Token::expressionBegin:
257          break;
258 <      case Token::charge:
194 <
258 >      case Token::expressionEnd:
259          break;
260 <      case Token::dipole:
261 <
262 <        break;
199 <      default:
200 <        unrecognizedAtomProperty(property);
201 <      }
202 <      bool match = false;
203 <      switch (comparator) {
204 <      case Token::opLT:
205 <        match = propertyValue < comparisonValue;
260 >      case Token::all:
261 >        bs = allInstruction();
262 >        stack.push(bs);            
263          break;
264 <      case Token::opLE:
265 <        match = propertyValue <= comparisonValue;
264 >      case Token::none:
265 >        bs = SelectionSet(nObjects);
266 >        stack.push(bs);            
267          break;
268 <      case Token::opGE:
269 <        match = propertyValue >= comparisonValue;
268 >      case Token::opOr:
269 >        bs = stack.top();
270 >        stack.pop();
271 >        stack.top() |= bs;
272          break;
273 <      case Token::opGT:
274 <        match = propertyValue > comparisonValue;
273 >      case Token::opAnd:
274 >        bs = stack.top();
275 >        stack.pop();
276 >        stack.top() &= bs;
277          break;
278 <      case Token::opEQ:
279 <        match = propertyValue == comparisonValue;
278 >      case Token::opNot:
279 >        stack.top().flip();
280          break;
281 +      case Token::within:
282 +        withinInstruction(instruction, stack.top(), frame);
283 +        break;
284 +      case Token::hull:
285 +        stack.push(hull(frame));
286 +        break;
287 +        //case Token::selected:
288 +        //  stack.push(getSelectionSet());
289 +        //  break;
290 +      case Token::name:
291 +        stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
292 +        break;
293 +      case Token::index:
294 +        stack.push(indexInstruction(instruction.value));
295 +        break;
296 +      case Token::identifier:
297 +        stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
298 +        break;
299 +      case Token::opLT:
300 +      case Token::opLE:
301 +      case Token::opGE:
302 +      case Token::opGT:
303 +      case Token::opEQ:
304        case Token::opNE:
305 <        match = propertyValue != comparisonValue;
305 >        stack.push(comparatorInstruction(instruction, frame));
306          break;
307 +      default:
308 +        unrecognizedExpression();
309        }
223      if (match)
224        bs.set(i);
310      }
311 +    if (stack.size() != 1)
312 +      evalError("atom expression compiler error - stack over/underflow");
313 +          
314 +    return stack.top();
315    }
316  
228 void SelectionEvaluator::withinInstruction(const Token& instruction, BitSet& bs, BitSet& bsResult)
317  
318 <    boost::any withinSpec = instruction.value;
319 <    if (withinSpec.type() == typeid(float)){
320 <        withinDistance(boost::any_cast<float>(withinSpec), bs, bsResult);
321 <        return;
322 <    }
318 >
319 >  SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
320 >    int comparator = instruction.tok;
321 >    int property = instruction.intValue;
322 >    float comparisonValue = boost::any_cast<float>(instruction.value);
323 >    SelectionSet bs = createSelectionSets();
324 >    bs.clearAll();
325      
326 <    evalError("Unrecognized within parameter:" + withinSpec);
327 < }
326 >    SimInfo::MoleculeIterator mi;
327 >    Molecule* mol;
328 >    Molecule::AtomIterator ai;
329 >    Atom* atom;
330 >    Molecule::RigidBodyIterator rbIter;
331 >    RigidBody* rb;
332  
333 <  void SelectionEvaluator::withinDistance(float distance, const BitSet& bs, const BitSet& bsResult) {
334 <    Frame frame = viewer.getFrame();
335 <    for (int i = frame.getAtomCount(); --i >= 0; ) {
336 <      if (bs.get(i)) {
337 <        Atom atom = frame.getAtomAt(i);
338 <        AtomIterator iterWithin =
339 <          frame.getWithinIterator(atom, distance);
246 <        while (iterWithin.hasNext())
247 <          bsResult.set(iterWithin.next().getAtomIndex());
333 >    for (mol = info->beginMolecule(mi); mol != NULL;
334 >         mol = info->nextMolecule(mi)) {
335 >
336 >      compareProperty(mol, bs, property, comparator, comparisonValue);
337 >
338 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
339 >        compareProperty(atom, bs, property, comparator, comparisonValue);
340        }
341 +    
342 +      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
343 +           rb = mol->nextRigidBody(rbIter)) {
344 +        compareProperty(rb, bs, property, comparator, comparisonValue);
345 +      }
346      }
347 +
348 +    return bs;
349    }
350  
351 <  void SelectionEvaluator::define() {
352 <    assert(statement.size() >= 3);
353 <
354 <    std::string variable = boost::any_cast<std::string>(statement[1].value);
351 >  SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction, int frame) {
352 >    int comparator = instruction.tok;
353 >    int property = instruction.intValue;
354 >    float comparisonValue = boost::any_cast<float>(instruction.value);
355 >    SelectionSet bs = createSelectionSets();
356 >    bs.clearAll();
357      
358 <    variables.insert(std::make_pair(variable, expression(statement, 2)));
359 <  }
360 < }
358 >    SimInfo::MoleculeIterator mi;
359 >    Molecule* mol;
360 >    Molecule::AtomIterator ai;
361 >    Atom* atom;
362 >    Molecule::RigidBodyIterator rbIter;
363 >    RigidBody* rb;
364  
365 < /** @todo */
366 < void SelectionEvaluator::predefine(const std::string& script) {
365 >    for (mol = info->beginMolecule(mi); mol != NULL;
366 >         mol = info->nextMolecule(mi)) {
367  
368 <    if (compiler.compile("#predefine", script)) {
265 <        std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
266 <        if (aatoken.size() != 1) {
267 <            evalError("predefinition does not have exactly 1 command:"
268 <                + script);
269 <            return;
270 <        }
271 <        std::vector<Token> statement = aatoken[0];
272 <        if (statement.size() > 2) {
273 <            int tok = statement[1].tok;
274 <            if (tok == Token::identifier || (tok & Token::predefinedset) == Token::predefinedset) {
275 <                std::string variable = (std::string)statement[1].value;
276 <                variables.insert(std::make_pair(variable, statement));
368 >      compareProperty(mol, bs, property, comparator, comparisonValue, frame);
369  
370 <            } else {
371 <                evalError("invalid variable name:" + script);
372 <            }
373 <        }else {
374 <            evalError("bad predefinition length:" + script);
375 <        }
370 >      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
371 >        compareProperty(atom, bs, property, comparator, comparisonValue, frame);
372 >      }
373 >    
374 >      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
375 >           rb = mol->nextRigidBody(rbIter)) {
376 >        compareProperty(rb, bs, property, comparator, comparisonValue, frame);
377 >      }
378 >    }
379  
380 +    return bs;
381 +  }
382 +
383 +  void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs,
384 +                                           int property, int comparator,
385 +                                           float comparisonValue) {
386 +    RealType propertyValue = 0.0;
387 +    Vector3d pos;
388 +
389 +    switch (property) {
390 +    case Token::mass:
391 +      propertyValue = sd->getMass();
392 +      break;
393 +    case Token::charge:
394 +      if (sd->isAtom()){
395 +        Atom* atom = static_cast<Atom*>(sd);
396 +        propertyValue = getCharge(atom);
397 +      } else if (sd->isRigidBody()) {
398 +        RigidBody* rb = static_cast<RigidBody*>(sd);
399 +        RigidBody::AtomIterator ai;
400 +        Atom* atom;
401 +        for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
402 +          propertyValue+=  getCharge(atom);
403 +        }
404 +      }
405 +      break;
406 +    case Token::x:
407 +      propertyValue = sd->getPos().x();
408 +      break;
409 +    case Token::y:
410 +      propertyValue = sd->getPos().y();
411 +      break;
412 +    case Token::z:
413 +      propertyValue = sd->getPos().z();
414 +      break;
415 +    case Token::wrappedX:    
416 +      pos = sd->getPos();
417 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
418 +      propertyValue = pos.x();
419 +      break;
420 +    case Token::wrappedY:
421 +      pos = sd->getPos();
422 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
423 +      propertyValue = pos.y();
424 +      break;
425 +    case Token::wrappedZ:
426 +      pos = sd->getPos();
427 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
428 +      propertyValue = pos.z();
429 +      break;
430 +    case Token::r:
431 +      propertyValue = sd->getPos().length();
432 +      break;
433 +    default:
434 +      unrecognizedAtomProperty(property);
435 +    }
436          
437 <    } else {
438 <        evalError("predefined set compile error:" + script +
439 <          "\ncompile error:" + compiler.getErrorMessage());
437 >    bool match = false;
438 >    switch (comparator) {
439 >    case Token::opLT:
440 >      match = propertyValue < comparisonValue;
441 >      break;
442 >    case Token::opLE:
443 >      match = propertyValue <= comparisonValue;
444 >      break;
445 >    case Token::opGE:
446 >      match = propertyValue >= comparisonValue;
447 >      break;
448 >    case Token::opGT:
449 >      match = propertyValue > comparisonValue;
450 >      break;
451 >    case Token::opEQ:
452 >      match = propertyValue == comparisonValue;
453 >      break;
454 >    case Token::opNE:
455 >      match = propertyValue != comparisonValue;
456 >      break;
457      }
458  
459 < }
459 >    if (match)
460 >      bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex());
461  
462 < void SelectionEvaluator::select(){
463 <    viewer.setSelectionSet(expression(statement, 1));
295 < }
462 >    
463 >  }
464  
465 < BitSet SelectionEvaluator::lookupValue(const std::string& variable){
466 <
467 <    std::map<std::string, boost::any>::iterator i = variables.find(variable);
465 >  void SelectionEvaluator::compareProperty(Molecule* mol, SelectionSet& bs,
466 >                                           int property, int comparator,
467 >                                           float comparisonValue) {
468 >    RealType propertyValue = 0.0;
469 >    Vector3d pos;
470  
471 <    if (i != variables.end()) {
472 <        if (i->second.type() == typeid(BitSet)) {
473 <            return boost::any_cast<BitSet>(i->second);
474 <        } else if (i->second.type() ==  typeid(std::vector<Token>)){
475 <            BitSet bs = expression(boost::any_cast(i->second), 2);
476 <            i->second =  bs; /**@todo fixme */
477 <            return bs;
471 >    switch (property) {
472 >    case Token::mass:
473 >      propertyValue = mol->getMass();
474 >      break;
475 >    case Token::charge:
476 >      {
477 >        Molecule::AtomIterator ai;
478 >        Atom* atom;
479 >        for (atom = mol->beginAtom(ai); atom != NULL;
480 >             atom = mol->nextAtom(ai)) {
481 >          propertyValue+=  getCharge(atom);
482          }
483 +      }
484 +      break;
485 +    case Token::x:
486 +      propertyValue = mol->getCom().x();
487 +      break;
488 +    case Token::y:
489 +      propertyValue = mol->getCom().y();
490 +      break;
491 +    case Token::z:
492 +      propertyValue = mol->getCom().z();
493 +      break;
494 +    case Token::wrappedX:    
495 +      pos = mol->getCom();
496 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
497 +      propertyValue = pos.x();
498 +      break;
499 +    case Token::wrappedY:
500 +      pos = mol->getCom();
501 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
502 +      propertyValue = pos.y();
503 +      break;
504 +    case Token::wrappedZ:
505 +      pos = mol->getCom();
506 +      info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos);
507 +      propertyValue = pos.z();
508 +      break;
509 +    case Token::r:
510 +      propertyValue = mol->getCom().length();
511 +      break;
512 +    default:
513 +      unrecognizedMoleculeProperty(property);
514      }
515 +    
516 +    bool match = false;
517 +    switch (comparator) {
518 +    case Token::opLT:
519 +      match = propertyValue < comparisonValue;
520 +      break;
521 +    case Token::opLE:
522 +      match = propertyValue <= comparisonValue;
523 +      break;
524 +    case Token::opGE:
525 +      match = propertyValue >= comparisonValue;
526 +      break;
527 +    case Token::opGT:
528 +      match = propertyValue > comparisonValue;
529 +      break;
530 +    case Token::opEQ:
531 +      match = propertyValue == comparisonValue;
532 +      break;
533 +    case Token::opNE:
534 +      match = propertyValue != comparisonValue;
535 +      break;
536 +    }
537 +    
538 +    if (match)
539 +      bs.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());    
540 +  }
541  
542 < }
542 >  void SelectionEvaluator::compareProperty(Molecule* mol, SelectionSet& bs,
543 >                                           int property, int comparator,
544 >                                           float comparisonValue, int frame) {
545 >    RealType propertyValue = 0.0;
546 >    Vector3d pos;
547 >    switch (property) {
548 >    case Token::mass:
549 >      propertyValue = mol->getMass();
550 >      break;
551 >    case Token::charge:
552 >      {
553 >        Molecule::AtomIterator ai;
554 >        Atom* atom;
555 >        for (atom = mol->beginAtom(ai); atom != NULL;
556 >             atom = mol->nextAtom(ai)) {
557 >          propertyValue+=  getCharge(atom,frame);
558 >        }
559 >      }      
560 >      break;
561 >    case Token::x:
562 >      propertyValue = mol->getCom(frame).x();
563 >      break;
564 >    case Token::y:
565 >      propertyValue = mol->getCom(frame).y();
566 >      break;
567 >    case Token::z:
568 >      propertyValue = mol->getCom(frame).z();
569 >      break;
570 >    case Token::wrappedX:    
571 >      pos = mol->getCom(frame);
572 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
573 >      propertyValue = pos.x();
574 >      break;
575 >    case Token::wrappedY:
576 >      pos = mol->getCom(frame);
577 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
578 >      propertyValue = pos.y();
579 >      break;
580 >    case Token::wrappedZ:
581 >      pos = mol->getCom(frame);
582 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
583 >      propertyValue = pos.z();
584 >      break;
585 >
586 >    case Token::r:
587 >      propertyValue = mol->getCom(frame).length();
588 >      break;
589 >    default:
590 >      unrecognizedMoleculeProperty(property);
591 >    }
592 >        
593 >    bool match = false;
594 >    switch (comparator) {
595 >    case Token::opLT:
596 >      match = propertyValue < comparisonValue;
597 >      break;
598 >    case Token::opLE:
599 >      match = propertyValue <= comparisonValue;
600 >      break;
601 >    case Token::opGE:
602 >      match = propertyValue >= comparisonValue;
603 >      break;
604 >    case Token::opGT:
605 >      match = propertyValue > comparisonValue;
606 >      break;
607 >    case Token::opEQ:
608 >      match = propertyValue == comparisonValue;
609 >      break;
610 >    case Token::opNE:
611 >      match = propertyValue != comparisonValue;
612 >      break;
613 >    }
614 >    if (match)
615 >      bs.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
616 >    
617 >    
618 >  }
619 >  void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs,
620 >                                           int property, int comparator,
621 >                                           float comparisonValue, int frame) {
622 >    RealType propertyValue = 0.0;
623 >    Vector3d pos;
624 >    switch (property) {
625 >    case Token::mass:
626 >      propertyValue = sd->getMass();
627 >      break;
628 >    case Token::charge:
629 >      if (sd->isAtom()){
630 >        Atom* atom = static_cast<Atom*>(sd);
631 >        propertyValue = getCharge(atom,frame);
632 >      } else if (sd->isRigidBody()) {
633 >        RigidBody* rb = static_cast<RigidBody*>(sd);
634 >        RigidBody::AtomIterator ai;
635 >        Atom* atom;
636 >        for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
637 >          propertyValue+=  getCharge(atom,frame);
638 >        }
639 >      }
640 >      break;
641 >    case Token::x:
642 >      propertyValue = sd->getPos(frame).x();
643 >      break;
644 >    case Token::y:
645 >      propertyValue = sd->getPos(frame).y();
646 >      break;
647 >    case Token::z:
648 >      propertyValue = sd->getPos(frame).z();
649 >      break;
650 >    case Token::wrappedX:    
651 >      pos = sd->getPos(frame);
652 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
653 >      propertyValue = pos.x();
654 >      break;
655 >    case Token::wrappedY:
656 >      pos = sd->getPos(frame);
657 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
658 >      propertyValue = pos.y();
659 >      break;
660 >    case Token::wrappedZ:
661 >      pos = sd->getPos(frame);
662 >      info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos);
663 >      propertyValue = pos.z();
664 >      break;
665 >
666 >    case Token::r:
667 >      propertyValue = sd->getPos(frame).length();
668 >      break;
669 >    default:
670 >      unrecognizedAtomProperty(property);
671 >    }
672 >        
673 >    bool match = false;
674 >    switch (comparator) {
675 >    case Token::opLT:
676 >      match = propertyValue < comparisonValue;
677 >      break;
678 >    case Token::opLE:
679 >      match = propertyValue <= comparisonValue;
680 >      break;
681 >    case Token::opGE:
682 >      match = propertyValue >= comparisonValue;
683 >      break;
684 >    case Token::opGT:
685 >      match = propertyValue > comparisonValue;
686 >      break;
687 >    case Token::opEQ:
688 >      match = propertyValue == comparisonValue;
689 >      break;
690 >    case Token::opNE:
691 >      match = propertyValue != comparisonValue;
692 >      break;
693 >    }
694 >    if (match)
695 >      bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex());    
696 >  }
697 >  
698 >  
699 >  void SelectionEvaluator::withinInstruction(const Token& instruction,
700 >                                             SelectionSet& bs){
701 >    
702 >    boost::any withinSpec = instruction.value;
703 >    float distance;
704 >    if (withinSpec.type() == typeid(float)){
705 >      distance = boost::any_cast<float>(withinSpec);
706 >    } else if (withinSpec.type() == typeid(int)) {
707 >      distance = boost::any_cast<int>(withinSpec);    
708 >    } else {
709 >      evalError("casting error in withinInstruction");
710 >      bs.clearAll();
711 >    }
712 >    
713 >    bs = distanceFinder.find(bs, distance);            
714 >  }
715 >
716 >  void SelectionEvaluator::withinInstruction(const Token& instruction,
717 >                                             SelectionSet& bs, int frame){
718 >    
719 >    boost::any withinSpec = instruction.value;
720 >    float distance;
721 >    if (withinSpec.type() == typeid(float)){
722 >      distance = boost::any_cast<float>(withinSpec);
723 >    } else if (withinSpec.type() == typeid(int)) {
724 >      distance = boost::any_cast<int>(withinSpec);    
725 >    } else {
726 >      evalError("casting error in withinInstruction");
727 >      bs.clearAll();
728 >    }
729 >    
730 >    bs = distanceFinder.find(bs, distance, frame);            
731 >  }
732 >  
733 >  void SelectionEvaluator::define() {
734 >    assert(statement.size() >= 3);
735 >    
736 >    std::string variable = boost::any_cast<std::string>(statement[1].value);
737 >    
738 >    variables.insert(VariablesType::value_type(variable,
739 >                                               expression(statement, 2)));
740 >  }
741 >  
742  
743 +  /** @todo */
744 +  void SelectionEvaluator::predefine(const std::string& script) {
745 +    
746 +    if (compiler.compile("#predefine", script)) {
747 +      std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
748 +      if (aatoken.size() != 1) {
749 +        evalError("predefinition does not have exactly 1 command:"
750 +                  + script);
751 +        return;
752 +      }
753 +      std::vector<Token> statement = aatoken[0];
754 +      if (statement.size() > 2) {
755 +        int tok = statement[1].tok;
756 +        if (tok == Token::identifier ||
757 +            (tok & Token::predefinedset) == Token::predefinedset) {
758 +          std::string variable = boost::any_cast<std::string>(statement[1].value);
759 +          variables.insert(VariablesType::value_type(variable, statement));
760 +          
761 +        } else {
762 +          evalError("invalid variable name:" + script);
763 +        }
764 +      }else {
765 +        evalError("bad predefinition length:" + script);
766 +      }      
767 +        
768 +    } else {
769 +      evalError("predefined set compile error:" + script +
770 +                "\ncompile error:" + compiler.getErrorMessage());
771 +    }
772 +  }
773 +
774 +  void SelectionEvaluator::select(SelectionSet& bs){
775 +    bs = expression(statement, 1);
776 +  }
777 +
778 +  void SelectionEvaluator::select(SelectionSet& bs, int frame){
779 +    bs = expression(statement, 1, frame);
780 +  }
781 +  
782 +  SelectionSet SelectionEvaluator::lookupValue(const std::string& variable){
783 +    
784 +    SelectionSet bs = createSelectionSets();
785 +    std::map<std::string, boost::any>::iterator i = variables.find(variable);
786 +    
787 +    if (i != variables.end()) {
788 +      if (i->second.type() == typeid(SelectionSet)) {
789 +        return boost::any_cast<SelectionSet>(i->second);
790 +      } else if (i->second.type() ==  typeid(std::vector<Token>)){
791 +        bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
792 +        i->second =  bs; /**@todo fixme */
793 +        return bs;
794 +      }
795 +    } else {
796 +      unrecognizedIdentifier(variable);
797 +    }
798 +    
799 +    return bs;
800 +  }
801 +  
802 +  SelectionSet SelectionEvaluator::nameInstruction(const std::string& name){    
803 +    return nameFinder.match(name);    
804 +  }    
805 +
806 +  bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
807 +    std::vector<Token>::const_iterator i;
808 +    for (i = tokens.begin(); i != tokens.end(); ++i) {
809 +      if (i->tok & Token::dynamic) {
810 +        return true;
811 +      }
812 +    }
813 +    
814 +    return false;
815 +  }    
816 +
817 +  void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
818 +    variables.clear();
819 +    //load predefine
820 +    //predefine();
821 +  }
822 +
823 +  SelectionSet SelectionEvaluator::createSelectionSets() {
824 +    SelectionSet ss(nObjects);
825 +    return ss;
826 +  }
827 +
828 +  SelectionSet SelectionEvaluator::evaluate() {
829 +    SelectionSet bs = createSelectionSets();
830 +    if (isLoaded_) {
831 +      pc = 0;
832 +      instructionDispatchLoop(bs);
833 +    }
834 +    return bs;
835 +  }
836 +
837 +  SelectionSet SelectionEvaluator::evaluate(int frame) {
838 +    SelectionSet bs = createSelectionSets();
839 +    if (isLoaded_) {
840 +      pc = 0;
841 +      instructionDispatchLoop(bs, frame);
842 +    }
843 +    return bs;
844 +  }
845 +
846 +  SelectionSet SelectionEvaluator::indexInstruction(const boost::any& value) {
847 +    SelectionSet bs = createSelectionSets();
848 +
849 +    if (value.type() == typeid(int)) {
850 +      int index = boost::any_cast<int>(value);
851 +      if (index < 0 || index >= bs.bitsets_[STUNTDOUBLE].size()) {
852 +        invalidIndex(index);
853 +      } else {
854 +        bs = indexFinder.find(index);
855 +      }
856 +    } else if (value.type() == typeid(std::pair<int, int>)) {
857 +      std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
858 +      assert(indexRange.first <= indexRange.second);
859 +      if (indexRange.first < 0 ||
860 +          indexRange.second >= bs.bitsets_[STUNTDOUBLE].size()) {
861 +        invalidIndexRange(indexRange);
862 +      }else {
863 +        bs = indexFinder.find(indexRange.first, indexRange.second);
864 +      }
865 +    }
866 +
867 +    return bs;
868 +  }
869 +
870 +  SelectionSet SelectionEvaluator::allInstruction() {
871 +    SelectionSet ss = createSelectionSets();
872 +
873 +    SimInfo::MoleculeIterator mi;
874 +    Molecule::AtomIterator ai;
875 +    Molecule::RigidBodyIterator rbIter;
876 +    Molecule::BondIterator bondIter;
877 +    Molecule::BendIterator bendIter;
878 +    Molecule::TorsionIterator torsionIter;
879 +    Molecule::InversionIterator inversionIter;
880 +
881 +    Molecule* mol;
882 +    Atom* atom;
883 +    RigidBody* rb;
884 +    Bond* bond;
885 +    Bend* bend;
886 +    Torsion* torsion;
887 +    Inversion* inversion;    
888 +
889 +    // Doing the loop insures that we're actually on this processor.
890 +
891 +    for (mol = info->beginMolecule(mi); mol != NULL;
892 +         mol = info->nextMolecule(mi)) {
893 +      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
894 +        ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
895 +      }    
896 +      for (rb = mol->beginRigidBody(rbIter); rb != NULL;
897 +           rb = mol->nextRigidBody(rbIter)) {
898 +        ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
899 +      }
900 +      for (bond = mol->beginBond(bondIter); bond != NULL;
901 +           bond = mol->nextBond(bondIter)) {
902 +        ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
903 +      }  
904 +      for (bend = mol->beginBend(bendIter); bend != NULL;
905 +           bend = mol->nextBend(bendIter)) {
906 +        ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
907 +      }  
908 +      for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
909 +           torsion = mol->nextTorsion(torsionIter)) {
910 +        ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
911 +      }  
912 +      for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
913 +           inversion = mol->nextInversion(inversionIter)) {
914 +        ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
915 +      }
916 +      ss.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex());
917 +    }
918 +
919 +    return ss;
920 +  }
921 +
922 +  SelectionSet SelectionEvaluator::hull() {
923 +    SelectionSet bs = createSelectionSets();
924 +    
925 +    bs = hullFinder.findHull();
926 +    surfaceArea_ = hullFinder.getSurfaceArea();
927 +    hasSurfaceArea_ = true;
928 +    return bs;
929 +  }
930 +
931 +
932 +  SelectionSet SelectionEvaluator::hull(int frame) {
933 +    SelectionSet bs = createSelectionSets();
934 +    
935 +    bs = hullFinder.findHull(frame);
936 +
937 +    return bs;
938 +  }
939 +
940 +  RealType SelectionEvaluator::getCharge(Atom* atom) {
941 +    RealType charge = 0.0;
942 +    AtomType* atomType = atom->getAtomType();
943 +
944 +    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
945 +    if ( fca.isFixedCharge() ) {
946 +      charge = fca.getCharge();
947 +    }
948 +    
949 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
950 +    if ( fqa.isFluctuatingCharge() ) {
951 +      charge += atom->getFlucQPos();
952 +    }
953 +    return charge;
954 +  }
955 +
956 +  RealType SelectionEvaluator::getCharge(Atom* atom, int frame) {
957 +    RealType charge = 0.0;
958 +    AtomType* atomType = atom->getAtomType();    
959 +
960 +    FixedChargeAdapter fca = FixedChargeAdapter(atomType);
961 +    if ( fca.isFixedCharge() ) {
962 +      charge = fca.getCharge();
963 +    }
964 +    
965 +    FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
966 +    if ( fqa.isFluctuatingCharge() ) {
967 +      charge += atom->getFlucQPos(frame);
968 +    }
969 +    return charge;
970 +  }
971 +
972   }

Comparing trunk/src/selection/SelectionEvaluator.cpp (property svn:keywords):
Revision 282 by tim, Thu Feb 3 14:04:59 2005 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines