| 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. | 
| 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> | 
| 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 |  |  | 
| 49 | – | namespace oopse { | 
| 53 |  |  | 
| 54 | + | namespace OpenMD { | 
| 55 |  |  | 
| 52 | – | SelectionEvaluator::SelectionEvaluator(SimInfo* si, const std::string& script) : info(si), finder(info){ | 
| 53 | – | } | 
| 56 |  |  | 
| 57 | < | bool SelectionEvaluator::loadScript(const std::string& filename, const std::string& script) { | 
| 57 | > | SelectionEvaluator::SelectionEvaluator(SimInfo* si) | 
| 58 | > | : info(si), nameFinder(info), distanceFinder(info), hullFinder(info), | 
| 59 | > | indexFinder(info), hasSurfaceArea_(false), | 
| 60 | > | isLoaded_(false){ | 
| 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 | > | } | 
| 67 | > |  | 
| 68 | > | bool SelectionEvaluator::loadScript(const std::string& filename, | 
| 69 | > | const std::string& script) { | 
| 70 | > | clearDefinitionsAndLoadPredefined(); | 
| 71 |  | this->filename = filename; | 
| 72 |  | this->script = script; | 
| 73 |  | if (! compiler.compile(filename, script)) { | 
| 74 | < | error = true; | 
| 75 | < | errorMessage = compiler.getErrorMessage(); | 
| 76 | < | return false; | 
| 74 | > | error = true; | 
| 75 | > | errorMessage = compiler.getErrorMessage(); | 
| 76 | > |  | 
| 77 | > | sprintf( painCave.errMsg, | 
| 78 | > | "SelectionCompiler Error: %s\n", errorMessage.c_str()); | 
| 79 | > | painCave.severity = OPENMD_ERROR; | 
| 80 | > | painCave.isFatal = 1; | 
| 81 | > | simError(); | 
| 82 | > | return false; | 
| 83 |  | } | 
| 84 |  |  | 
| 85 |  | pc = 0; | 
| 86 |  | aatoken = compiler.getAatokenCompiled(); | 
| 87 |  | linenumbers = compiler.getLineNumbers(); | 
| 88 |  | lineIndices = compiler.getLineIndices(); | 
| 89 | + |  | 
| 90 | + | std::vector<std::vector<Token> >::const_iterator i; | 
| 91 | + |  | 
| 92 | + | isDynamic_ = false; | 
| 93 | + | for (i = aatoken.begin(); i != aatoken.end(); ++i) { | 
| 94 | + | if (containDynamicToken(*i)) { | 
| 95 | + | isDynamic_ = true; | 
| 96 | + | break; | 
| 97 | + | } | 
| 98 | + | } | 
| 99 | + |  | 
| 100 | + | isLoaded_ = true; | 
| 101 |  | return true; | 
| 102 | < | } | 
| 102 | > | } | 
| 103 |  |  | 
| 104 | < | void SelectionEvaluator::clearState() { | 
| 72 | < | for (int i = scriptLevelMax; --i >= 0; ) | 
| 73 | < | stack[i].clear(); | 
| 74 | < | scriptLevel = 0; | 
| 104 | > | void SelectionEvaluator::clearState() { | 
| 105 |  | error = false; | 
| 106 |  | errorMessage = ""; | 
| 107 | < | } | 
| 107 | > | } | 
| 108 |  |  | 
| 109 | < | bool SelectionEvaluator::loadScriptString(const std::string& script) { | 
| 109 | > | bool SelectionEvaluator::loadScriptString(const std::string& script) { | 
| 110 |  | clearState(); | 
| 111 |  | return loadScript("", script); | 
| 112 | < | } | 
| 112 | > | } | 
| 113 |  |  | 
| 114 | < | bool SelectionEvaluator::loadScriptFile(const std::string& filename) { | 
| 114 | > | bool SelectionEvaluator::loadScriptFile(const std::string& filename) { | 
| 115 |  | clearState(); | 
| 116 |  | return loadScriptFileInternal(filename); | 
| 117 | < | } | 
| 117 | > | } | 
| 118 |  |  | 
| 119 | < | bool SelectionEvaluator::loadScriptFileInternal(const  std::string & filename) { | 
| 120 | < | return true; /**@todo */ | 
| 121 | < | } | 
| 119 | > | bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) { | 
| 120 | > | ifstrstream ifs(filename.c_str()); | 
| 121 | > | if (!ifs.is_open()) { | 
| 122 | > | return false; | 
| 123 | > | } | 
| 124 | > |  | 
| 125 | > | const int bufferSize = 65535; | 
| 126 | > | char buffer[bufferSize]; | 
| 127 | > | std::string script; | 
| 128 | > | while(ifs.getline(buffer, bufferSize)) { | 
| 129 | > | script += buffer; | 
| 130 | > | } | 
| 131 | > | return loadScript(filename, script); | 
| 132 | > | } | 
| 133 | > |  | 
| 134 | > | void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs){ | 
| 135 | > |  | 
| 136 | > | while ( pc < aatoken.size()) { | 
| 137 | > | statement = aatoken[pc++]; | 
| 138 | > | statementLength = statement.size(); | 
| 139 | > | Token token = statement[0]; | 
| 140 | > | switch (token.tok) { | 
| 141 | > | case Token::define: | 
| 142 | > | define(); | 
| 143 | > | break; | 
| 144 | > | case Token::select: | 
| 145 | > | select(bs); | 
| 146 | > | break; | 
| 147 | > | default: | 
| 148 | > | unrecognizedCommand(token); | 
| 149 | > | return; | 
| 150 | > | } | 
| 151 | > | } | 
| 152 |  |  | 
| 153 | < | void SelectionEvaluator::instructionDispatchLoop(){ | 
| 153 | > | } | 
| 154 |  |  | 
| 155 | + | void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs, int frame){ | 
| 156 | + |  | 
| 157 |  | while ( pc < aatoken.size()) { | 
| 158 | < | statement = aatoken[pc++]; | 
| 159 | < | statementLength = statement.size(); | 
| 160 | < | Token token = statement[0]; | 
| 161 | < | switch (token.tok) { | 
| 162 | < | case Token::define: | 
| 163 | < | define(); | 
| 164 | < | break; | 
| 165 | < | case Token::select: | 
| 166 | < | select(); | 
| 167 | < | break; | 
| 168 | < | default: | 
| 169 | < | unrecognizedCommand(token); | 
| 170 | < | return; | 
| 171 | < | } | 
| 158 | > | statement = aatoken[pc++]; | 
| 159 | > | statementLength = statement.size(); | 
| 160 | > | Token token = statement[0]; | 
| 161 | > | switch (token.tok) { | 
| 162 | > | case Token::define: | 
| 163 | > | define(); | 
| 164 | > | break; | 
| 165 | > | case Token::select: | 
| 166 | > | select(bs, frame); | 
| 167 | > | break; | 
| 168 | > | default: | 
| 169 | > | unrecognizedCommand(token); | 
| 170 | > | return; | 
| 171 | > | } | 
| 172 |  | } | 
| 111 | – | } | 
| 173 |  |  | 
| 174 | < | BitSet SelectionEvaluator::expression(const std::vector<Token>& code, int pcStart) { | 
| 175 | < | BitSet bs; | 
| 176 | < | std::stack<BitSet> stack; | 
| 177 | < |  | 
| 178 | < | for (int pc = pcStart; ; ++pc) { | 
| 174 | > | } | 
| 175 | > |  | 
| 176 | > | SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code, | 
| 177 | > | int pcStart) { | 
| 178 | > | SelectionSet bs = createSelectionSets(); | 
| 179 | > | std::stack<SelectionSet> stack; | 
| 180 | > | vector<int> bsSize = bs.size(); | 
| 181 | > |  | 
| 182 | > | for (unsigned int pc = pcStart; pc < code.size(); ++pc) { | 
| 183 |  | Token instruction = code[pc]; | 
| 184 |  |  | 
| 185 |  | switch (instruction.tok) { | 
| 188 |  | case Token::expressionEnd: | 
| 189 |  | break; | 
| 190 |  | case Token::all: | 
| 191 | < | bs = BitSet(nStuntDouble); | 
| 192 | < | bs.setAll(); | 
| 191 | > | bs = allInstruction(); | 
| 192 | > | stack.push(bs); | 
| 193 | > | break; | 
| 194 | > | case Token::none: | 
| 195 | > | bs = createSelectionSets(); | 
| 196 | > | stack.push(bs); | 
| 197 | > | break; | 
| 198 | > | case Token::opOr: | 
| 199 | > | bs= stack.top(); | 
| 200 | > | stack.pop(); | 
| 201 | > | stack.top() |= bs; | 
| 202 | > | break; | 
| 203 | > | case Token::opAnd: | 
| 204 | > | bs = stack.top(); | 
| 205 | > | stack.pop(); | 
| 206 | > | stack.top() &= bs; | 
| 207 | > | break; | 
| 208 | > | case Token::opNot: | 
| 209 | > | stack.top().flip(); | 
| 210 | > | break; | 
| 211 | > | case Token::within: | 
| 212 | > | withinInstruction(instruction, stack.top()); | 
| 213 | > | break; | 
| 214 | > | case Token::hull: | 
| 215 | > | stack.push(hull()); | 
| 216 | > | break; | 
| 217 | > | //case Token::selected: | 
| 218 | > | //  stack.push(getSelectionSet()); | 
| 219 | > | //  break; | 
| 220 | > | case Token::name: | 
| 221 | > | stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value))); | 
| 222 | > | break; | 
| 223 | > | case Token::index: | 
| 224 | > | stack.push(indexInstruction(instruction.value)); | 
| 225 | > | break; | 
| 226 | > | case Token::identifier: | 
| 227 | > | stack.push(lookupValue(boost::any_cast<std::string>(instruction.value))); | 
| 228 | > | break; | 
| 229 | > | case Token::opLT: | 
| 230 | > | case Token::opLE: | 
| 231 | > | case Token::opGE: | 
| 232 | > | case Token::opGT: | 
| 233 | > | case Token::opEQ: | 
| 234 | > | case Token::opNE: | 
| 235 | > | stack.push(comparatorInstruction(instruction)); | 
| 236 | > | break; | 
| 237 | > | default: | 
| 238 | > | unrecognizedExpression(); | 
| 239 | > | } | 
| 240 | > | } | 
| 241 | > | if (stack.size() != 1) | 
| 242 | > | evalError("atom expression compiler error - stack over/underflow"); | 
| 243 | > |  | 
| 244 | > | return stack.top(); | 
| 245 | > | } | 
| 246 | > |  | 
| 247 | > |  | 
| 248 | > | SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code, | 
| 249 | > | int pcStart, int frame) { | 
| 250 | > | SelectionSet bs = createSelectionSets(); | 
| 251 | > | std::stack<SelectionSet> stack; | 
| 252 | > |  | 
| 253 | > | for (unsigned int pc = pcStart; pc < code.size(); ++pc) { | 
| 254 | > | Token instruction = code[pc]; | 
| 255 | > |  | 
| 256 | > | switch (instruction.tok) { | 
| 257 | > | case Token::expressionBegin: | 
| 258 | > | break; | 
| 259 | > | case Token::expressionEnd: | 
| 260 | > | break; | 
| 261 | > | case Token::all: | 
| 262 | > | bs = allInstruction(); | 
| 263 |  | stack.push(bs); | 
| 264 |  | break; | 
| 265 |  | case Token::none: | 
| 266 | < | bs = BitSet(nStuntDouble); | 
| 266 | > | bs = SelectionSet(nObjects); | 
| 267 |  | stack.push(bs); | 
| 268 |  | break; | 
| 269 |  | case Token::opOr: | 
| 280 |  | stack.top().flip(); | 
| 281 |  | break; | 
| 282 |  | case Token::within: | 
| 283 | < |  | 
| 149 | < | withinInstruction(instruction, stack.top()); | 
| 283 | > | withinInstruction(instruction, stack.top(), frame); | 
| 284 |  | break; | 
| 285 | < | //case Token::selected: | 
| 286 | < | //  stack.push(getSelectionSet()); | 
| 287 | < | //  break; | 
| 285 | > | case Token::hull: | 
| 286 | > | stack.push(hull(frame)); | 
| 287 | > | break; | 
| 288 | > | //case Token::selected: | 
| 289 | > | //  stack.push(getSelectionSet()); | 
| 290 | > | //  break; | 
| 291 |  | case Token::name: | 
| 292 |  | stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value))); | 
| 293 |  | break; | 
| 294 | + | case Token::index: | 
| 295 | + | stack.push(indexInstruction(instruction.value)); | 
| 296 |  | break; | 
| 297 |  | case Token::identifier: | 
| 298 |  | stack.push(lookupValue(boost::any_cast<std::string>(instruction.value))); | 
| 303 |  | case Token::opGT: | 
| 304 |  | case Token::opEQ: | 
| 305 |  | case Token::opNE: | 
| 306 | < | stack.push(comparatorInstruction(instruction)); | 
| 306 | > | stack.push(comparatorInstruction(instruction, frame)); | 
| 307 |  | break; | 
| 308 |  | default: | 
| 309 |  | unrecognizedExpression(); | 
| 317 |  |  | 
| 318 |  |  | 
| 319 |  |  | 
| 320 | < | BitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { | 
| 320 | > | SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { | 
| 321 |  | int comparator = instruction.tok; | 
| 322 |  | int property = instruction.intValue; | 
| 323 |  | float comparisonValue = boost::any_cast<float>(instruction.value); | 
| 324 | < | float propertyValue; | 
| 325 | < | BitSet bs(nStuntDouble); | 
| 324 | > | SelectionSet bs = createSelectionSets(); | 
| 325 | > | bs.clearAll(); | 
| 326 |  |  | 
| 327 |  | SimInfo::MoleculeIterator mi; | 
| 328 |  | Molecule* mol; | 
| 330 |  | Atom* atom; | 
| 331 |  | Molecule::RigidBodyIterator rbIter; | 
| 332 |  | RigidBody* rb; | 
| 194 | – |  | 
| 195 | – | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 333 |  |  | 
| 334 | < | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 335 | < | compareProperty(atom, bs, property, comparator, comparisonValue); | 
| 336 | < | } | 
| 337 | < |  | 
| 338 | < | //change the positions of atoms which belong to the rigidbodies | 
| 339 | < | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 340 | < | compareProperty(rb, bs, property, comparator, comparisonValue); | 
| 341 | < | } | 
| 334 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 335 | > | mol = info->nextMolecule(mi)) { | 
| 336 | > |  | 
| 337 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 338 | > | compareProperty(atom, bs, property, comparator, comparisonValue); | 
| 339 | > | } | 
| 340 | > |  | 
| 341 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 342 | > | rb = mol->nextRigidBody(rbIter)) { | 
| 343 | > | compareProperty(rb, bs, property, comparator, comparisonValue); | 
| 344 | > | } | 
| 345 |  | } | 
| 346 |  |  | 
| 347 |  | return bs; | 
| 348 | < | } | 
| 348 | > | } | 
| 349 |  |  | 
| 350 | < | void SelectionEvaluator::compareProperty(StuntDouble* sd, BitSet& bs, int property, int comparator, float comparisonValue) { | 
| 351 | < | double propertyValue; | 
| 352 | < | switch (property) { | 
| 353 | < | case Token::mass: | 
| 354 | < | propertyValue = sd->getMass(); | 
| 355 | < | break; | 
| 216 | < | case Token::charge: | 
| 217 | < | return; | 
| 218 | < | //break; | 
| 219 | < | case Token::dipole: | 
| 220 | < | return; | 
| 221 | < | //break; | 
| 222 | < | default: | 
| 223 | < | unrecognizedAtomProperty(property); | 
| 224 | < | } | 
| 225 | < |  | 
| 226 | < | bool match = false; | 
| 227 | < | switch (comparator) { | 
| 228 | < | case Token::opLT: | 
| 229 | < | match = propertyValue < comparisonValue; | 
| 230 | < | break; | 
| 231 | < | case Token::opLE: | 
| 232 | < | match = propertyValue <= comparisonValue; | 
| 233 | < | break; | 
| 234 | < | case Token::opGE: | 
| 235 | < | match = propertyValue >= comparisonValue; | 
| 236 | < | break; | 
| 237 | < | case Token::opGT: | 
| 238 | < | match = propertyValue > comparisonValue; | 
| 239 | < | break; | 
| 240 | < | case Token::opEQ: | 
| 241 | < | match = propertyValue == comparisonValue; | 
| 242 | < | break; | 
| 243 | < | case Token::opNE: | 
| 244 | < | match = propertyValue != comparisonValue; | 
| 245 | < | break; | 
| 246 | < | } | 
| 247 | < | if (match) | 
| 248 | < | bs.setBitOn(sd->getGlobalIndex()); | 
| 249 | < |  | 
| 250 | < | } | 
| 251 | < |  | 
| 252 | < | void SelectionEvaluator::withinInstruction(const Token& instruction, BitSet& bs){ | 
| 253 | < |  | 
| 254 | < | boost::any withinSpec = instruction.value; | 
| 255 | < | if (withinSpec.type() == typeid(float)){ | 
| 256 | < | // | 
| 257 | < | return; | 
| 258 | < | } | 
| 350 | > | SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction, int frame) { | 
| 351 | > | int comparator = instruction.tok; | 
| 352 | > | int property = instruction.intValue; | 
| 353 | > | float comparisonValue = boost::any_cast<float>(instruction.value); | 
| 354 | > | SelectionSet bs = createSelectionSets(); | 
| 355 | > | bs.clearAll(); | 
| 356 |  |  | 
| 357 | < | evalError("Unrecognized within parameter"); | 
| 358 | < | } | 
| 357 | > | SimInfo::MoleculeIterator mi; | 
| 358 | > | Molecule* mol; | 
| 359 | > | Molecule::AtomIterator ai; | 
| 360 | > | Atom* atom; | 
| 361 | > | Molecule::RigidBodyIterator rbIter; | 
| 362 | > | RigidBody* rb; | 
| 363 |  |  | 
| 364 | < | void SelectionEvaluator::define() { | 
| 365 | < | assert(statement.size() >= 3); | 
| 364 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 365 | > | mol = info->nextMolecule(mi)) { | 
| 366 |  |  | 
| 367 | < | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 368 | < |  | 
| 369 | < | variables.insert(std::make_pair(variable, expression(statement, 2))); | 
| 370 | < | } | 
| 371 | < |  | 
| 372 | < |  | 
| 373 | < | /** @todo */ | 
| 374 | < | void SelectionEvaluator::predefine(const std::string& script) { | 
| 274 | < |  | 
| 275 | < | if (compiler.compile("#predefine", script)) { | 
| 276 | < | std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled(); | 
| 277 | < | if (aatoken.size() != 1) { | 
| 278 | < | evalError("predefinition does not have exactly 1 command:" | 
| 279 | < | + script); | 
| 280 | < | return; | 
| 281 | < | } | 
| 282 | < | std::vector<Token> statement = aatoken[0]; | 
| 283 | < | if (statement.size() > 2) { | 
| 284 | < | int tok = statement[1].tok; | 
| 285 | < | if (tok == Token::identifier || (tok & Token::predefinedset) == Token::predefinedset) { | 
| 286 | < | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 287 | < | variables.insert(std::make_pair(variable, statement)); | 
| 288 | < |  | 
| 289 | < | } else { | 
| 290 | < | evalError("invalid variable name:" + script); | 
| 291 | < | } | 
| 292 | < | }else { | 
| 293 | < | evalError("bad predefinition length:" + script); | 
| 294 | < | } | 
| 295 | < |  | 
| 296 | < |  | 
| 297 | < | } else { | 
| 298 | < | evalError("predefined set compile error:" + script + | 
| 299 | < | "\ncompile error:" + compiler.getErrorMessage()); | 
| 367 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 368 | > | compareProperty(atom, bs, property, comparator, comparisonValue, frame); | 
| 369 | > | } | 
| 370 | > |  | 
| 371 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 372 | > | rb = mol->nextRigidBody(rbIter)) { | 
| 373 | > | compareProperty(rb, bs, property, comparator, comparisonValue, frame); | 
| 374 | > | } | 
| 375 |  | } | 
| 376 |  |  | 
| 377 | < | } | 
| 377 | > | return bs; | 
| 378 | > | } | 
| 379 |  |  | 
| 380 | < | void SelectionEvaluator::select(){ | 
| 381 | < | //viewer.setSelectionSet(expression(statement, 1)); | 
| 382 | < | } | 
| 380 | > | void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs, | 
| 381 | > | int property, int comparator, | 
| 382 | > | float comparisonValue) { | 
| 383 | > | RealType propertyValue = 0.0; | 
| 384 | > | Vector3d pos; | 
| 385 |  |  | 
| 386 | < | BitSet SelectionEvaluator::lookupValue(const std::string& variable){ | 
| 386 | > | switch (property) { | 
| 387 | > | case Token::mass: | 
| 388 | > | propertyValue = sd->getMass(); | 
| 389 | > | break; | 
| 390 | > | case Token::charge: | 
| 391 | > | if (sd->isAtom()){ | 
| 392 | > | Atom* atom = static_cast<Atom*>(sd); | 
| 393 | > | propertyValue = getCharge(atom); | 
| 394 | > | } else if (sd->isRigidBody()) { | 
| 395 | > | RigidBody* rb = static_cast<RigidBody*>(sd); | 
| 396 | > | RigidBody::AtomIterator ai; | 
| 397 | > | Atom* atom; | 
| 398 | > | for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) { | 
| 399 | > | propertyValue+=  getCharge(atom); | 
| 400 | > | } | 
| 401 | > | } | 
| 402 | > | break; | 
| 403 | > | case Token::x: | 
| 404 | > | propertyValue = sd->getPos().x(); | 
| 405 | > | break; | 
| 406 | > | case Token::y: | 
| 407 | > | propertyValue = sd->getPos().y(); | 
| 408 | > | break; | 
| 409 | > | case Token::z: | 
| 410 | > | propertyValue = sd->getPos().z(); | 
| 411 | > | break; | 
| 412 | > | case Token::wrappedX: | 
| 413 | > | pos = sd->getPos(); | 
| 414 | > | info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); | 
| 415 | > | propertyValue = pos.x(); | 
| 416 | > | break; | 
| 417 | > | case Token::wrappedY: | 
| 418 | > | pos = sd->getPos(); | 
| 419 | > | info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); | 
| 420 | > | propertyValue = pos.y(); | 
| 421 | > | break; | 
| 422 | > | case Token::wrappedZ: | 
| 423 | > | pos = sd->getPos(); | 
| 424 | > | info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); | 
| 425 | > | propertyValue = pos.z(); | 
| 426 | > | break; | 
| 427 | > | case Token::r: | 
| 428 | > | propertyValue = sd->getPos().length(); | 
| 429 | > | break; | 
| 430 | > | default: | 
| 431 | > | unrecognizedAtomProperty(property); | 
| 432 | > | } | 
| 433 | > |  | 
| 434 | > | bool match = false; | 
| 435 | > | switch (comparator) { | 
| 436 | > | case Token::opLT: | 
| 437 | > | match = propertyValue < comparisonValue; | 
| 438 | > | break; | 
| 439 | > | case Token::opLE: | 
| 440 | > | match = propertyValue <= comparisonValue; | 
| 441 | > | break; | 
| 442 | > | case Token::opGE: | 
| 443 | > | match = propertyValue >= comparisonValue; | 
| 444 | > | break; | 
| 445 | > | case Token::opGT: | 
| 446 | > | match = propertyValue > comparisonValue; | 
| 447 | > | break; | 
| 448 | > | case Token::opEQ: | 
| 449 | > | match = propertyValue == comparisonValue; | 
| 450 | > | break; | 
| 451 | > | case Token::opNE: | 
| 452 | > | match = propertyValue != comparisonValue; | 
| 453 | > | break; | 
| 454 | > | } | 
| 455 |  |  | 
| 456 | < | std::map<std::string, boost::any>::iterator i = variables.find(variable); | 
| 456 | > | if (match) | 
| 457 | > | bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex()); | 
| 458 |  |  | 
| 459 | + |  | 
| 460 | + | } | 
| 461 | + |  | 
| 462 | + | void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs, | 
| 463 | + | int property, int comparator, | 
| 464 | + | float comparisonValue, int frame) { | 
| 465 | + | RealType propertyValue = 0.0; | 
| 466 | + | Vector3d pos; | 
| 467 | + | switch (property) { | 
| 468 | + | case Token::mass: | 
| 469 | + | propertyValue = sd->getMass(); | 
| 470 | + | break; | 
| 471 | + | case Token::charge: | 
| 472 | + | if (sd->isAtom()){ | 
| 473 | + | Atom* atom = static_cast<Atom*>(sd); | 
| 474 | + | propertyValue = getCharge(atom,frame); | 
| 475 | + | } else if (sd->isRigidBody()) { | 
| 476 | + | RigidBody* rb = static_cast<RigidBody*>(sd); | 
| 477 | + | RigidBody::AtomIterator ai; | 
| 478 | + | Atom* atom; | 
| 479 | + | for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) { | 
| 480 | + | propertyValue+=  getCharge(atom,frame); | 
| 481 | + | } | 
| 482 | + | } | 
| 483 | + | break; | 
| 484 | + | case Token::x: | 
| 485 | + | propertyValue = sd->getPos(frame).x(); | 
| 486 | + | break; | 
| 487 | + | case Token::y: | 
| 488 | + | propertyValue = sd->getPos(frame).y(); | 
| 489 | + | break; | 
| 490 | + | case Token::z: | 
| 491 | + | propertyValue = sd->getPos(frame).z(); | 
| 492 | + | break; | 
| 493 | + | case Token::wrappedX: | 
| 494 | + | pos = sd->getPos(frame); | 
| 495 | + | info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); | 
| 496 | + | propertyValue = pos.x(); | 
| 497 | + | break; | 
| 498 | + | case Token::wrappedY: | 
| 499 | + | pos = sd->getPos(frame); | 
| 500 | + | info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); | 
| 501 | + | propertyValue = pos.y(); | 
| 502 | + | break; | 
| 503 | + | case Token::wrappedZ: | 
| 504 | + | pos = sd->getPos(frame); | 
| 505 | + | info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); | 
| 506 | + | propertyValue = pos.z(); | 
| 507 | + | break; | 
| 508 | + |  | 
| 509 | + | case Token::r: | 
| 510 | + | propertyValue = sd->getPos(frame).length(); | 
| 511 | + | break; | 
| 512 | + | default: | 
| 513 | + | unrecognizedAtomProperty(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 | + | if (match) | 
| 538 | + | bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex()); | 
| 539 | + |  | 
| 540 | + |  | 
| 541 | + | } | 
| 542 | + |  | 
| 543 | + | void SelectionEvaluator::withinInstruction(const Token& instruction, | 
| 544 | + | SelectionSet& bs){ | 
| 545 | + |  | 
| 546 | + | boost::any withinSpec = instruction.value; | 
| 547 | + | float distance; | 
| 548 | + | if (withinSpec.type() == typeid(float)){ | 
| 549 | + | distance = boost::any_cast<float>(withinSpec); | 
| 550 | + | } else if (withinSpec.type() == typeid(int)) { | 
| 551 | + | distance = boost::any_cast<int>(withinSpec); | 
| 552 | + | } else { | 
| 553 | + | evalError("casting error in withinInstruction"); | 
| 554 | + | bs.clearAll(); | 
| 555 | + | } | 
| 556 | + |  | 
| 557 | + | bs = distanceFinder.find(bs, distance); | 
| 558 | + | } | 
| 559 | + |  | 
| 560 | + | void SelectionEvaluator::withinInstruction(const Token& instruction, | 
| 561 | + | SelectionSet& bs, int frame){ | 
| 562 | + |  | 
| 563 | + | boost::any withinSpec = instruction.value; | 
| 564 | + | float distance; | 
| 565 | + | if (withinSpec.type() == typeid(float)){ | 
| 566 | + | distance = boost::any_cast<float>(withinSpec); | 
| 567 | + | } else if (withinSpec.type() == typeid(int)) { | 
| 568 | + | distance = boost::any_cast<int>(withinSpec); | 
| 569 | + | } else { | 
| 570 | + | evalError("casting error in withinInstruction"); | 
| 571 | + | bs.clearAll(); | 
| 572 | + | } | 
| 573 | + |  | 
| 574 | + | bs = distanceFinder.find(bs, distance, frame); | 
| 575 | + | } | 
| 576 | + |  | 
| 577 | + | void SelectionEvaluator::define() { | 
| 578 | + | assert(statement.size() >= 3); | 
| 579 | + |  | 
| 580 | + | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 581 | + |  | 
| 582 | + | variables.insert(VariablesType::value_type(variable, | 
| 583 | + | expression(statement, 2))); | 
| 584 | + | } | 
| 585 | + |  | 
| 586 | + |  | 
| 587 | + | /** @todo */ | 
| 588 | + | void SelectionEvaluator::predefine(const std::string& script) { | 
| 589 | + |  | 
| 590 | + | if (compiler.compile("#predefine", script)) { | 
| 591 | + | std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled(); | 
| 592 | + | if (aatoken.size() != 1) { | 
| 593 | + | evalError("predefinition does not have exactly 1 command:" | 
| 594 | + | + script); | 
| 595 | + | return; | 
| 596 | + | } | 
| 597 | + | std::vector<Token> statement = aatoken[0]; | 
| 598 | + | if (statement.size() > 2) { | 
| 599 | + | int tok = statement[1].tok; | 
| 600 | + | if (tok == Token::identifier || | 
| 601 | + | (tok & Token::predefinedset) == Token::predefinedset) { | 
| 602 | + | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 603 | + | variables.insert(VariablesType::value_type(variable, statement)); | 
| 604 | + |  | 
| 605 | + | } else { | 
| 606 | + | evalError("invalid variable name:" + script); | 
| 607 | + | } | 
| 608 | + | }else { | 
| 609 | + | evalError("bad predefinition length:" + script); | 
| 610 | + | } | 
| 611 | + |  | 
| 612 | + | } else { | 
| 613 | + | evalError("predefined set compile error:" + script + | 
| 614 | + | "\ncompile error:" + compiler.getErrorMessage()); | 
| 615 | + | } | 
| 616 | + | } | 
| 617 | + |  | 
| 618 | + | void SelectionEvaluator::select(SelectionSet& bs){ | 
| 619 | + | bs = expression(statement, 1); | 
| 620 | + | } | 
| 621 | + |  | 
| 622 | + | void SelectionEvaluator::select(SelectionSet& bs, int frame){ | 
| 623 | + | bs = expression(statement, 1, frame); | 
| 624 | + | } | 
| 625 | + |  | 
| 626 | + | SelectionSet SelectionEvaluator::lookupValue(const std::string& variable){ | 
| 627 | + |  | 
| 628 | + | SelectionSet bs = createSelectionSets(); | 
| 629 | + | std::map<std::string, boost::any>::iterator i = variables.find(variable); | 
| 630 | + |  | 
| 631 |  | if (i != variables.end()) { | 
| 632 | < | if (i->second.type() == typeid(BitSet)) { | 
| 633 | < | return boost::any_cast<BitSet>(i->second); | 
| 634 | < | } else if (i->second.type() ==  typeid(std::vector<Token>)){ | 
| 635 | < | BitSet bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2); | 
| 636 | < | i->second =  bs; /**@todo fixme */ | 
| 637 | < | return bs; | 
| 638 | < | } | 
| 632 | > | if (i->second.type() == typeid(SelectionSet)) { | 
| 633 | > | return boost::any_cast<SelectionSet>(i->second); | 
| 634 | > | } else if (i->second.type() ==  typeid(std::vector<Token>)){ | 
| 635 | > | bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2); | 
| 636 | > | i->second =  bs; /**@todo fixme */ | 
| 637 | > | return bs; | 
| 638 | > | } | 
| 639 |  | } else { | 
| 640 | < | unrecognizedIdentifier(variable); | 
| 640 | > | unrecognizedIdentifier(variable); | 
| 641 |  | } | 
| 642 | < | } | 
| 642 | > |  | 
| 643 | > | return bs; | 
| 644 | > | } | 
| 645 | > |  | 
| 646 | > | SelectionSet SelectionEvaluator::nameInstruction(const std::string& name){ | 
| 647 | > | return nameFinder.match(name); | 
| 648 | > | } | 
| 649 |  |  | 
| 650 | < | BitSet SelectionEvaluator::nameInstruction(const std::string& name){ | 
| 651 | < | BitSet bs(nStuntDouble); | 
| 650 | > | bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){ | 
| 651 | > | std::vector<Token>::const_iterator i; | 
| 652 | > | for (i = tokens.begin(); i != tokens.end(); ++i) { | 
| 653 | > | if (i->tok & Token::dynamic) { | 
| 654 | > | return true; | 
| 655 | > | } | 
| 656 | > | } | 
| 657 |  |  | 
| 658 | < | bool hasError = finder.match(name, bs); | 
| 658 | > | return false; | 
| 659 | > | } | 
| 660 | > |  | 
| 661 | > | void SelectionEvaluator::clearDefinitionsAndLoadPredefined() { | 
| 662 | > | variables.clear(); | 
| 663 | > | //load predefine | 
| 664 | > | //predefine(); | 
| 665 | > | } | 
| 666 | > |  | 
| 667 | > | SelectionSet SelectionEvaluator::createSelectionSets() { | 
| 668 | > | SelectionSet ss(nObjects); | 
| 669 | > | return ss; | 
| 670 | > | } | 
| 671 | > |  | 
| 672 | > | SelectionSet SelectionEvaluator::evaluate() { | 
| 673 | > | SelectionSet bs = createSelectionSets(); | 
| 674 | > | if (isLoaded_) { | 
| 675 | > | pc = 0; | 
| 676 | > | instructionDispatchLoop(bs); | 
| 677 | > | } | 
| 678 | > | return bs; | 
| 679 | > | } | 
| 680 | > |  | 
| 681 | > | SelectionSet SelectionEvaluator::evaluate(int frame) { | 
| 682 | > | SelectionSet bs = createSelectionSets(); | 
| 683 | > | if (isLoaded_) { | 
| 684 | > | pc = 0; | 
| 685 | > | instructionDispatchLoop(bs, frame); | 
| 686 | > | } | 
| 687 | > | return bs; | 
| 688 | > | } | 
| 689 | > |  | 
| 690 | > | SelectionSet SelectionEvaluator::indexInstruction(const boost::any& value) { | 
| 691 | > | SelectionSet bs = createSelectionSets(); | 
| 692 | > |  | 
| 693 | > | if (value.type() == typeid(int)) { | 
| 694 | > | int index = boost::any_cast<int>(value); | 
| 695 | > | if (index < 0 || index >= bs.bitsets_[STUNTDOUBLE].size()) { | 
| 696 | > | invalidIndex(index); | 
| 697 | > | } else { | 
| 698 | > | bs = indexFinder.find(index); | 
| 699 | > | } | 
| 700 | > | } else if (value.type() == typeid(std::pair<int, int>)) { | 
| 701 | > | std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value); | 
| 702 | > | assert(indexRange.first <= indexRange.second); | 
| 703 | > | if (indexRange.first < 0 || | 
| 704 | > | indexRange.second >= bs.bitsets_[STUNTDOUBLE].size()) { | 
| 705 | > | invalidIndexRange(indexRange); | 
| 706 | > | }else { | 
| 707 | > | bs = indexFinder.find(indexRange.first, indexRange.second); | 
| 708 | > | } | 
| 709 | > | } | 
| 710 | > |  | 
| 711 | > | return bs; | 
| 712 | > | } | 
| 713 | > |  | 
| 714 | > | SelectionSet SelectionEvaluator::allInstruction() { | 
| 715 | > | SelectionSet ss = createSelectionSets(); | 
| 716 | > |  | 
| 717 | > | SimInfo::MoleculeIterator mi; | 
| 718 | > | Molecule::AtomIterator ai; | 
| 719 | > | Molecule::RigidBodyIterator rbIter; | 
| 720 | > | Molecule::BondIterator bondIter; | 
| 721 | > | Molecule::BendIterator bendIter; | 
| 722 | > | Molecule::TorsionIterator torsionIter; | 
| 723 | > | Molecule::InversionIterator inversionIter; | 
| 724 | > |  | 
| 725 | > | Molecule* mol; | 
| 726 | > | Atom* atom; | 
| 727 | > | RigidBody* rb; | 
| 728 | > | Bond* bond; | 
| 729 | > | Bend* bend; | 
| 730 | > | Torsion* torsion; | 
| 731 | > | Inversion* inversion; | 
| 732 | > |  | 
| 733 | > | // Doing the loop insures that we're actually on this processor. | 
| 734 | > |  | 
| 735 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 736 | > | mol = info->nextMolecule(mi)) { | 
| 737 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 738 | > | ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex()); | 
| 739 | > | } | 
| 740 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 741 | > | rb = mol->nextRigidBody(rbIter)) { | 
| 742 | > | ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex()); | 
| 743 | > | } | 
| 744 | > | for (bond = mol->beginBond(bondIter); bond != NULL; | 
| 745 | > | bond = mol->nextBond(bondIter)) { | 
| 746 | > | ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex()); | 
| 747 | > | } | 
| 748 | > | for (bend = mol->beginBend(bendIter); bend != NULL; | 
| 749 | > | bend = mol->nextBend(bendIter)) { | 
| 750 | > | ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex()); | 
| 751 | > | } | 
| 752 | > | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; | 
| 753 | > | torsion = mol->nextTorsion(torsionIter)) { | 
| 754 | > | ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex()); | 
| 755 | > | } | 
| 756 | > | for (inversion = mol->beginInversion(inversionIter); inversion != NULL; | 
| 757 | > | inversion = mol->nextInversion(inversionIter)) { | 
| 758 | > | ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex()); | 
| 759 | > | } | 
| 760 | > | } | 
| 761 | > |  | 
| 762 | > | return ss; | 
| 763 | > | } | 
| 764 | > |  | 
| 765 | > | SelectionSet SelectionEvaluator::hull() { | 
| 766 | > | SelectionSet bs = createSelectionSets(); | 
| 767 |  |  | 
| 768 | < | return bs; | 
| 769 | < | } | 
| 768 | > | bs = hullFinder.findHull(); | 
| 769 | > | surfaceArea_ = hullFinder.getSurfaceArea(); | 
| 770 | > | hasSurfaceArea_ = true; | 
| 771 | > | return bs; | 
| 772 | > | } | 
| 773 |  |  | 
| 774 |  |  | 
| 775 | + | SelectionSet SelectionEvaluator::hull(int frame) { | 
| 776 | + | SelectionSet bs = createSelectionSets(); | 
| 777 | + |  | 
| 778 | + | bs = hullFinder.findHull(frame); | 
| 779 | + |  | 
| 780 | + | return bs; | 
| 781 | + | } | 
| 782 | + |  | 
| 783 | + | RealType SelectionEvaluator::getCharge(Atom* atom) { | 
| 784 | + | RealType charge = 0.0; | 
| 785 | + | AtomType* atomType = atom->getAtomType(); | 
| 786 | + |  | 
| 787 | + | FixedChargeAdapter fca = FixedChargeAdapter(atomType); | 
| 788 | + | if ( fca.isFixedCharge() ) { | 
| 789 | + | charge = fca.getCharge(); | 
| 790 | + | } | 
| 791 | + |  | 
| 792 | + | FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); | 
| 793 | + | if ( fqa.isFluctuatingCharge() ) { | 
| 794 | + | charge += atom->getFlucQPos(); | 
| 795 | + | } | 
| 796 | + | return charge; | 
| 797 | + | } | 
| 798 | + |  | 
| 799 | + | RealType SelectionEvaluator::getCharge(Atom* atom, int frame) { | 
| 800 | + | RealType charge = 0.0; | 
| 801 | + | AtomType* atomType = atom->getAtomType(); | 
| 802 | + |  | 
| 803 | + | FixedChargeAdapter fca = FixedChargeAdapter(atomType); | 
| 804 | + | if ( fca.isFixedCharge() ) { | 
| 805 | + | charge = fca.getCharge(); | 
| 806 | + | } | 
| 807 | + |  | 
| 808 | + | FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); | 
| 809 | + | if ( fqa.isFluctuatingCharge() ) { | 
| 810 | + | charge += atom->getFlucQPos(frame); | 
| 811 | + | } | 
| 812 | + | return charge; | 
| 813 | + | } | 
| 814 | + |  | 
| 815 |  | } |