| 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, 24107 (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 |  |  | 
| 51 | + | namespace OpenMD { | 
| 52 |  |  | 
| 53 | + |  | 
| 54 | + | SelectionEvaluator::SelectionEvaluator(SimInfo* si) | 
| 55 | + | : info(si), nameFinder(info), distanceFinder(info), hullFinder(info), | 
| 56 | + | indexFinder(info), | 
| 57 | + | isLoaded_(false){ | 
| 58 | + | nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies(); | 
| 59 | + | } | 
| 60 | + |  | 
| 61 | + | bool SelectionEvaluator::loadScript(const std::string& filename, | 
| 62 | + | const std::string& script) { | 
| 63 | + | clearDefinitionsAndLoadPredefined(); | 
| 64 | + | this->filename = filename; | 
| 65 | + | this->script = script; | 
| 66 | + | if (! compiler.compile(filename, script)) { | 
| 67 | + | error = true; | 
| 68 | + | errorMessage = compiler.getErrorMessage(); | 
| 69 | + |  | 
| 70 | + | sprintf( painCave.errMsg, | 
| 71 | + | "SelectionCompiler Error: %s\n", errorMessage.c_str()); | 
| 72 | + | painCave.severity = OPENMD_ERROR; | 
| 73 | + | painCave.isFatal = 1; | 
| 74 | + | simError(); | 
| 75 | + | return false; | 
| 76 | + | } | 
| 77 | + |  | 
| 78 | + | pc = 0; | 
| 79 | + | aatoken = compiler.getAatokenCompiled(); | 
| 80 | + | linenumbers = compiler.getLineNumbers(); | 
| 81 | + | lineIndices = compiler.getLineIndices(); | 
| 82 | + |  | 
| 83 | + | std::vector<std::vector<Token> >::const_iterator i; | 
| 84 | + |  | 
| 85 | + | isDynamic_ = false; | 
| 86 | + | for (i = aatoken.begin(); i != aatoken.end(); ++i) { | 
| 87 | + | if (containDynamicToken(*i)) { | 
| 88 | + | isDynamic_ = true; | 
| 89 | + | break; | 
| 90 | + | } | 
| 91 | + | } | 
| 92 | + |  | 
| 93 | + | isLoaded_ = true; | 
| 94 | + | return true; | 
| 95 | + | } | 
| 96 | + |  | 
| 97 | + | void SelectionEvaluator::clearState() { | 
| 98 | + | error = false; | 
| 99 | + | errorMessage = ""; | 
| 100 | + | } | 
| 101 | + |  | 
| 102 | + | bool SelectionEvaluator::loadScriptString(const std::string& script) { | 
| 103 | + | clearState(); | 
| 104 | + | return loadScript("", script); | 
| 105 | + | } | 
| 106 | + |  | 
| 107 | + | bool SelectionEvaluator::loadScriptFile(const std::string& filename) { | 
| 108 | + | clearState(); | 
| 109 | + | return loadScriptFileInternal(filename); | 
| 110 | + | } | 
| 111 | + |  | 
| 112 | + | bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) { | 
| 113 | + | ifstrstream ifs(filename.c_str()); | 
| 114 | + | if (!ifs.is_open()) { | 
| 115 | + | return false; | 
| 116 | + | } | 
| 117 | + |  | 
| 118 | + | const int bufferSize = 65535; | 
| 119 | + | char buffer[bufferSize]; | 
| 120 | + | std::string script; | 
| 121 | + | while(ifs.getline(buffer, bufferSize)) { | 
| 122 | + | script += buffer; | 
| 123 | + | } | 
| 124 | + | return loadScript(filename, script); | 
| 125 | + | } | 
| 126 | + |  | 
| 127 | + | void SelectionEvaluator::instructionDispatchLoop(OpenMDBitSet& bs){ | 
| 128 | + |  | 
| 129 | + | while ( pc < aatoken.size()) { | 
| 130 | + | statement = aatoken[pc++]; | 
| 131 | + | statementLength = statement.size(); | 
| 132 | + | Token token = statement[0]; | 
| 133 | + | switch (token.tok) { | 
| 134 | + | case Token::define: | 
| 135 | + | define(); | 
| 136 | + | break; | 
| 137 | + | case Token::select: | 
| 138 | + | select(bs); | 
| 139 | + | break; | 
| 140 | + | default: | 
| 141 | + | unrecognizedCommand(token); | 
| 142 | + | return; | 
| 143 | + | } | 
| 144 | + | } | 
| 145 | + |  | 
| 146 | + | } | 
| 147 | + |  | 
| 148 | + | OpenMDBitSet SelectionEvaluator::expression(const std::vector<Token>& code, | 
| 149 | + | int pcStart) { | 
| 150 | + | OpenMDBitSet bs; | 
| 151 | + | std::stack<OpenMDBitSet> stack; | 
| 152 | + |  | 
| 153 | + | for (unsigned int pc = pcStart; pc < code.size(); ++pc) { | 
| 154 | + | Token instruction = code[pc]; | 
| 155 | + |  | 
| 156 | + | switch (instruction.tok) { | 
| 157 | + | case Token::expressionBegin: | 
| 158 | + | break; | 
| 159 | + | case Token::expressionEnd: | 
| 160 | + | break; | 
| 161 | + | case Token::all: | 
| 162 | + | bs = OpenMDBitSet(nStuntDouble); | 
| 163 | + | bs.setAll(); | 
| 164 | + | stack.push(bs); | 
| 165 | + | break; | 
| 166 | + | case Token::none: | 
| 167 | + | bs = OpenMDBitSet(nStuntDouble); | 
| 168 | + | stack.push(bs); | 
| 169 | + | break; | 
| 170 | + | case Token::opOr: | 
| 171 | + | bs = stack.top(); | 
| 172 | + | stack.pop(); | 
| 173 | + | stack.top() |= bs; | 
| 174 | + | break; | 
| 175 | + | case Token::opAnd: | 
| 176 | + | bs = stack.top(); | 
| 177 | + | stack.pop(); | 
| 178 | + | stack.top() &= bs; | 
| 179 | + | break; | 
| 180 | + | case Token::opNot: | 
| 181 | + | stack.top().flip(); | 
| 182 | + | break; | 
| 183 | + | case Token::within: | 
| 184 | + | withinInstruction(instruction, stack.top()); | 
| 185 | + | break; | 
| 186 | + | case Token::hull: | 
| 187 | + | stack.push(hull()); | 
| 188 | + | break; | 
| 189 | + | //case Token::selected: | 
| 190 | + | //  stack.push(getSelectionSet()); | 
| 191 | + | //  break; | 
| 192 | + | case Token::name: | 
| 193 | + | stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value))); | 
| 194 | + | break; | 
| 195 | + | case Token::index: | 
| 196 | + | stack.push(indexInstruction(instruction.value)); | 
| 197 | + | break; | 
| 198 | + | case Token::identifier: | 
| 199 | + | stack.push(lookupValue(boost::any_cast<std::string>(instruction.value))); | 
| 200 | + | break; | 
| 201 | + | case Token::opLT: | 
| 202 | + | case Token::opLE: | 
| 203 | + | case Token::opGE: | 
| 204 | + | case Token::opGT: | 
| 205 | + | case Token::opEQ: | 
| 206 | + | case Token::opNE: | 
| 207 | + | stack.push(comparatorInstruction(instruction)); | 
| 208 | + | break; | 
| 209 | + | default: | 
| 210 | + | unrecognizedExpression(); | 
| 211 | + | } | 
| 212 | + | } | 
| 213 | + | if (stack.size() != 1) | 
| 214 | + | evalError("atom expression compiler error - stack over/underflow"); | 
| 215 | + |  | 
| 216 | + | return stack.top(); | 
| 217 | + | } | 
| 218 | + |  | 
| 219 | + |  | 
| 220 | + |  | 
| 221 | + | OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { | 
| 222 | + | int comparator = instruction.tok; | 
| 223 | + | int property = instruction.intValue; | 
| 224 | + | float comparisonValue = boost::any_cast<float>(instruction.value); | 
| 225 | + | OpenMDBitSet bs(nStuntDouble); | 
| 226 | + | bs.clearAll(); | 
| 227 | + |  | 
| 228 | + | SimInfo::MoleculeIterator mi; | 
| 229 | + | Molecule* mol; | 
| 230 | + | Molecule::AtomIterator ai; | 
| 231 | + | Atom* atom; | 
| 232 | + | Molecule::RigidBodyIterator rbIter; | 
| 233 | + | RigidBody* rb; | 
| 234 | + |  | 
| 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 | + | } | 
| 250 | + |  | 
| 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 | + | case Token::r: | 
| 282 | + | propertyValue = sd->getPos().length(); | 
| 283 | + | break; | 
| 284 | + | default: | 
| 285 | + | unrecognizedAtomProperty(property); | 
| 286 | + | } | 
| 287 | + |  | 
| 288 | + | bool match = false; | 
| 289 | + | switch (comparator) { | 
| 290 | + | case Token::opLT: | 
| 291 | + | match = propertyValue < comparisonValue; | 
| 292 | + | break; | 
| 293 | + | case Token::opLE: | 
| 294 | + | match = propertyValue <= comparisonValue; | 
| 295 | + | break; | 
| 296 | + | case Token::opGE: | 
| 297 | + | match = propertyValue >= comparisonValue; | 
| 298 | + | break; | 
| 299 | + | case Token::opGT: | 
| 300 | + | match = propertyValue > comparisonValue; | 
| 301 | + | break; | 
| 302 | + | case Token::opEQ: | 
| 303 | + | match = propertyValue == comparisonValue; | 
| 304 | + | break; | 
| 305 | + | case Token::opNE: | 
| 306 | + | match = propertyValue != comparisonValue; | 
| 307 | + | break; | 
| 308 | + | } | 
| 309 | + | if (match) | 
| 310 | + | bs.setBitOn(sd->getGlobalIndex()); | 
| 311 | + |  | 
| 312 | + | } | 
| 313 | + |  | 
| 314 | + | void SelectionEvaluator::withinInstruction(const Token& instruction, | 
| 315 | + | OpenMDBitSet& bs){ | 
| 316 | + |  | 
| 317 | + | boost::any withinSpec = instruction.value; | 
| 318 | + | float distance; | 
| 319 | + | if (withinSpec.type() == typeid(float)){ | 
| 320 | + | distance = boost::any_cast<float>(withinSpec); | 
| 321 | + | } else if (withinSpec.type() == typeid(int)) { | 
| 322 | + | distance = boost::any_cast<int>(withinSpec); | 
| 323 | + | } else { | 
| 324 | + | evalError("casting error in withinInstruction"); | 
| 325 | + | bs.clearAll(); | 
| 326 | + | } | 
| 327 | + |  | 
| 328 | + | bs = distanceFinder.find(bs, distance); | 
| 329 | + | } | 
| 330 | + |  | 
| 331 | + | void SelectionEvaluator::define() { | 
| 332 | + | assert(statement.size() >= 3); | 
| 333 | + |  | 
| 334 | + | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 335 | + |  | 
| 336 | + | variables.insert(VariablesType::value_type(variable, | 
| 337 | + | expression(statement, 2))); | 
| 338 | + | } | 
| 339 | + |  | 
| 340 | + |  | 
| 341 | + | /** @todo */ | 
| 342 | + | void SelectionEvaluator::predefine(const std::string& script) { | 
| 343 | + |  | 
| 344 | + | if (compiler.compile("#predefine", script)) { | 
| 345 | + | std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled(); | 
| 346 | + | if (aatoken.size() != 1) { | 
| 347 | + | evalError("predefinition does not have exactly 1 command:" | 
| 348 | + | + script); | 
| 349 | + | return; | 
| 350 | + | } | 
| 351 | + | std::vector<Token> statement = aatoken[0]; | 
| 352 | + | if (statement.size() > 2) { | 
| 353 | + | int tok = statement[1].tok; | 
| 354 | + | if (tok == Token::identifier || | 
| 355 | + | (tok & Token::predefinedset) == Token::predefinedset) { | 
| 356 | + | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 357 | + | variables.insert(VariablesType::value_type(variable, statement)); | 
| 358 | + |  | 
| 359 | + | } else { | 
| 360 | + | evalError("invalid variable name:" + script); | 
| 361 | + | } | 
| 362 | + | }else { | 
| 363 | + | evalError("bad predefinition length:" + script); | 
| 364 | + | } | 
| 365 | + |  | 
| 366 | + | } else { | 
| 367 | + | evalError("predefined set compile error:" + script + | 
| 368 | + | "\ncompile error:" + compiler.getErrorMessage()); | 
| 369 | + | } | 
| 370 | + | } | 
| 371 | + |  | 
| 372 | + | void SelectionEvaluator::select(OpenMDBitSet& bs){ | 
| 373 | + | bs = expression(statement, 1); | 
| 374 | + | } | 
| 375 | + |  | 
| 376 | + | OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){ | 
| 377 | + |  | 
| 378 | + | OpenMDBitSet bs(nStuntDouble); | 
| 379 | + | std::map<std::string, boost::any>::iterator i = variables.find(variable); | 
| 380 | + |  | 
| 381 | + | if (i != variables.end()) { | 
| 382 | + | if (i->second.type() == typeid(OpenMDBitSet)) { | 
| 383 | + | return boost::any_cast<OpenMDBitSet>(i->second); | 
| 384 | + | } else if (i->second.type() ==  typeid(std::vector<Token>)){ | 
| 385 | + | bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2); | 
| 386 | + | i->second =  bs; /**@todo fixme */ | 
| 387 | + | return bs; | 
| 388 | + | } | 
| 389 | + | } else { | 
| 390 | + | unrecognizedIdentifier(variable); | 
| 391 | + | } | 
| 392 | + |  | 
| 393 | + | return bs; | 
| 394 | + | } | 
| 395 | + |  | 
| 396 | + | OpenMDBitSet SelectionEvaluator::nameInstruction(const std::string& name){ | 
| 397 | + | return nameFinder.match(name); | 
| 398 | + | } | 
| 399 | + |  | 
| 400 | + | bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){ | 
| 401 | + | std::vector<Token>::const_iterator i; | 
| 402 | + | for (i = tokens.begin(); i != tokens.end(); ++i) { | 
| 403 | + | if (i->tok & Token::dynamic) { | 
| 404 | + | return true; | 
| 405 | + | } | 
| 406 | + | } | 
| 407 | + |  | 
| 408 | + | return false; | 
| 409 | + | } | 
| 410 | + |  | 
| 411 | + | void SelectionEvaluator::clearDefinitionsAndLoadPredefined() { | 
| 412 | + | variables.clear(); | 
| 413 | + | //load predefine | 
| 414 | + | //predefine(); | 
| 415 | + | } | 
| 416 | + |  | 
| 417 | + | OpenMDBitSet SelectionEvaluator::evaluate() { | 
| 418 | + | OpenMDBitSet bs(nStuntDouble); | 
| 419 | + | if (isLoaded_) { | 
| 420 | + | pc = 0; | 
| 421 | + | instructionDispatchLoop(bs); | 
| 422 | + | } | 
| 423 | + |  | 
| 424 | + | return bs; | 
| 425 | + | } | 
| 426 | + |  | 
| 427 | + | OpenMDBitSet SelectionEvaluator::indexInstruction(const boost::any& value) { | 
| 428 | + | OpenMDBitSet bs(nStuntDouble); | 
| 429 | + |  | 
| 430 | + | if (value.type() == typeid(int)) { | 
| 431 | + | int index = boost::any_cast<int>(value); | 
| 432 | + | if (index < 0 || index >= bs.size()) { | 
| 433 | + | invalidIndex(index); | 
| 434 | + | } else { | 
| 435 | + | bs = indexFinder.find(index); | 
| 436 | + | } | 
| 437 | + | } else if (value.type() == typeid(std::pair<int, int>)) { | 
| 438 | + | std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value); | 
| 439 | + | assert(indexRange.first <= indexRange.second); | 
| 440 | + | if (indexRange.first < 0 || indexRange.second >= bs.size()) { | 
| 441 | + | invalidIndexRange(indexRange); | 
| 442 | + | }else { | 
| 443 | + | bs = indexFinder.find(indexRange.first, indexRange.second); | 
| 444 | + | } | 
| 445 | + | } | 
| 446 | + |  | 
| 447 | + | return bs; | 
| 448 | + | } | 
| 449 | + |  | 
| 450 | + |  | 
| 451 | + | OpenMDBitSet SelectionEvaluator::hull() { | 
| 452 | + | OpenMDBitSet bs(nStuntDouble); | 
| 453 | + |  | 
| 454 | + | bs = hullFinder.findHull(); | 
| 455 | + |  | 
| 456 | + | return bs; | 
| 457 | + | } | 
| 458 | + |  | 
| 459 | + |  | 
| 460 | + |  | 
| 461 | + | RealType SelectionEvaluator::getCharge(Atom* atom) { | 
| 462 | + | RealType charge =0.0; | 
| 463 | + | AtomType* atomType = atom->getAtomType(); | 
| 464 | + | if (atomType->isCharge()) { | 
| 465 | + | GenericData* data = atomType->getPropertyByName("Charge"); | 
| 466 | + | if (data != NULL) { | 
| 467 | + | DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data); | 
| 468 | + |  | 
| 469 | + | if (doubleData != NULL) { | 
| 470 | + | charge = doubleData->getData(); | 
| 471 | + |  | 
| 472 | + | } else { | 
| 473 | + | sprintf( painCave.errMsg, | 
| 474 | + | "Can not cast GenericData to DoubleGenericData\n"); | 
| 475 | + | painCave.severity = OPENMD_ERROR; | 
| 476 | + | painCave.isFatal = 1; | 
| 477 | + | simError(); | 
| 478 | + | } | 
| 479 | + | } | 
| 480 | + | } | 
| 481 | + |  | 
| 482 | + | return charge; | 
| 483 | + | } | 
| 484 | + |  | 
| 485 |  | } |