| 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 |  |  | 
| 56 |  | SelectionEvaluator::SelectionEvaluator(SimInfo* si) | 
| 57 | < | : info(si), nameFinder(info), distanceFinder(info), indexFinder(info), isLoaded_(false){ | 
| 58 | < |  | 
| 59 | < | nStuntDouble = info->getNGlobalAtoms() + info->getNRigidBodies(); | 
| 60 | < | } | 
| 61 | < |  | 
| 62 | < | bool SelectionEvaluator::loadScript(const std::string& filename, const std::string& script) { | 
| 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 | < | std::cerr << "SelectionCompiler Error: " << errorMessage << std::endl; | 
| 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 |  |  | 
| 115 |  | return loadScriptFileInternal(filename); | 
| 116 |  | } | 
| 117 |  |  | 
| 118 | < | bool SelectionEvaluator::loadScriptFileInternal(const  std::string & filename) { | 
| 119 | < | std::ifstream ifs(filename.c_str()); | 
| 118 | > | bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) { | 
| 119 | > | ifstrstream ifs(filename.c_str()); | 
| 120 |  | if (!ifs.is_open()) { | 
| 121 |  | return false; | 
| 122 |  | } | 
| 123 | < |  | 
| 123 | > |  | 
| 124 |  | const int bufferSize = 65535; | 
| 125 |  | char buffer[bufferSize]; | 
| 126 |  | std::string script; | 
| 129 |  | } | 
| 130 |  | return loadScript(filename, script); | 
| 131 |  | } | 
| 132 | < |  | 
| 133 | < | void SelectionEvaluator::instructionDispatchLoop(BitSet& bs){ | 
| 132 | > |  | 
| 133 | > | void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs){ | 
| 134 |  |  | 
| 135 |  | while ( pc < aatoken.size()) { | 
| 136 |  | statement = aatoken[pc++]; | 
| 151 |  |  | 
| 152 |  | } | 
| 153 |  |  | 
| 154 | < | BitSet SelectionEvaluator::expression(const std::vector<Token>& code, int pcStart) { | 
| 140 | < | BitSet bs; | 
| 141 | < | std::stack<BitSet> stack; | 
| 154 | > | void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs, int frame){ | 
| 155 |  |  | 
| 156 | < | for (int pc = pcStart; pc < code.size(); ++pc) { | 
| 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(bs, frame); | 
| 166 | > | break; | 
| 167 | > | default: | 
| 168 | > | unrecognizedCommand(token); | 
| 169 | > | return; | 
| 170 | > | } | 
| 171 | > | } | 
| 172 | > |  | 
| 173 | > | } | 
| 174 | > |  | 
| 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) { | 
| 187 |  | case Token::expressionEnd: | 
| 188 |  | break; | 
| 189 |  | case Token::all: | 
| 190 | < | bs = BitSet(nStuntDouble); | 
| 191 | < | bs.setAll(); | 
| 190 | > | bs = allInstruction(); | 
| 191 | > | stack.push(bs); | 
| 192 | > | break; | 
| 193 | > | case Token::none: | 
| 194 | > | bs = createSelectionSets(); | 
| 195 | > | stack.push(bs); | 
| 196 | > | break; | 
| 197 | > | case Token::opOr: | 
| 198 | > | bs= stack.top(); | 
| 199 | > | stack.pop(); | 
| 200 | > | stack.top() |= bs; | 
| 201 | > | break; | 
| 202 | > | case Token::opAnd: | 
| 203 | > | bs = stack.top(); | 
| 204 | > | stack.pop(); | 
| 205 | > | stack.top() &= bs; | 
| 206 | > | break; | 
| 207 | > | case Token::opNot: | 
| 208 | > | stack.top().flip(); | 
| 209 | > | break; | 
| 210 | > | case Token::within: | 
| 211 | > | withinInstruction(instruction, stack.top()); | 
| 212 | > | break; | 
| 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 | > | stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value))); | 
| 221 | > | break; | 
| 222 | > | case Token::index: | 
| 223 | > | stack.push(indexInstruction(instruction.value)); | 
| 224 | > | break; | 
| 225 | > | case Token::identifier: | 
| 226 | > | stack.push(lookupValue(boost::any_cast<std::string>(instruction.value))); | 
| 227 | > | break; | 
| 228 | > | case Token::opLT: | 
| 229 | > | case Token::opLE: | 
| 230 | > | case Token::opGE: | 
| 231 | > | case Token::opGT: | 
| 232 | > | case Token::opEQ: | 
| 233 | > | case Token::opNE: | 
| 234 | > | stack.push(comparatorInstruction(instruction)); | 
| 235 | > | break; | 
| 236 | > | default: | 
| 237 | > | unrecognizedExpression(); | 
| 238 | > | } | 
| 239 | > | } | 
| 240 | > | if (stack.size() != 1) | 
| 241 | > | evalError("atom expression compiler error - stack over/underflow"); | 
| 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 | > | switch (instruction.tok) { | 
| 256 | > | case Token::expressionBegin: | 
| 257 | > | break; | 
| 258 | > | case Token::expressionEnd: | 
| 259 | > | break; | 
| 260 | > | case Token::all: | 
| 261 | > | bs = allInstruction(); | 
| 262 |  | stack.push(bs); | 
| 263 |  | break; | 
| 264 |  | case Token::none: | 
| 265 | < | bs = BitSet(nStuntDouble); | 
| 265 | > | bs = SelectionSet(nObjects); | 
| 266 |  | stack.push(bs); | 
| 267 |  | break; | 
| 268 |  | case Token::opOr: | 
| 279 |  | stack.top().flip(); | 
| 280 |  | break; | 
| 281 |  | case Token::within: | 
| 282 | < |  | 
| 175 | < | withinInstruction(instruction, stack.top()); | 
| 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; | 
| 302 |  | case Token::opGT: | 
| 303 |  | case Token::opEQ: | 
| 304 |  | case Token::opNE: | 
| 305 | < | stack.push(comparatorInstruction(instruction)); | 
| 305 | > | stack.push(comparatorInstruction(instruction, frame)); | 
| 306 |  | break; | 
| 307 |  | default: | 
| 308 |  | unrecognizedExpression(); | 
| 316 |  |  | 
| 317 |  |  | 
| 318 |  |  | 
| 319 | < | BitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { | 
| 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 | < | float propertyValue; | 
| 214 | < | BitSet bs(nStuntDouble); | 
| 323 | > | SelectionSet bs = createSelectionSets(); | 
| 324 |  | bs.clearAll(); | 
| 325 |  |  | 
| 326 |  | SimInfo::MoleculeIterator mi; | 
| 329 |  | Atom* atom; | 
| 330 |  | Molecule::RigidBodyIterator rbIter; | 
| 331 |  | RigidBody* rb; | 
| 223 | – |  | 
| 224 | – | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 332 |  |  | 
| 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; rb = mol->nextRigidBody(rbIter)) { | 
| 343 | < | compareProperty(rb, bs, property, comparator, comparisonValue); | 
| 344 | < | } | 
| 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::compareProperty(StuntDouble* sd, BitSet& bs, int property, int comparator, float comparisonValue) { | 
| 352 | < | double propertyValue = 0.0; | 
| 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 | > | SimInfo::MoleculeIterator mi; | 
| 359 | > | Molecule* mol; | 
| 360 | > | Molecule::AtomIterator ai; | 
| 361 | > | Atom* atom; | 
| 362 | > | Molecule::RigidBodyIterator rbIter; | 
| 363 | > | RigidBody* rb; | 
| 364 | > |  | 
| 365 | > | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 366 | > | mol = info->nextMolecule(mi)) { | 
| 367 | > |  | 
| 368 | > | compareProperty(mol, bs, property, comparator, comparisonValue, frame); | 
| 369 | > |  | 
| 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(); | 
| 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 |  | } | 
| 455 |  | match = propertyValue != comparisonValue; | 
| 456 |  | break; | 
| 457 |  | } | 
| 282 | – | if (match) | 
| 283 | – | bs.setBitOn(sd->getGlobalIndex()); | 
| 458 |  |  | 
| 459 | + | if (match) | 
| 460 | + | bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex()); | 
| 461 | + |  | 
| 462 | + |  | 
| 463 |  | } | 
| 464 |  |  | 
| 465 | < | void SelectionEvaluator::withinInstruction(const Token& instruction, BitSet& bs){ | 
| 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 | > | 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 | + | 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)){ | 
| 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 | < |  | 
| 735 | > |  | 
| 736 |  | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 737 | < |  | 
| 738 | < | variables.insert(VariablesType::value_type(variable, expression(statement, 2))); | 
| 737 | > |  | 
| 738 | > | variables.insert(VariablesType::value_type(variable, | 
| 739 | > | expression(statement, 2))); | 
| 740 |  | } | 
| 741 | + |  | 
| 742 |  |  | 
| 311 | – |  | 
| 743 |  | /** @todo */ | 
| 744 |  | void SelectionEvaluator::predefine(const std::string& script) { | 
| 745 | < |  | 
| 745 | > |  | 
| 746 |  | if (compiler.compile("#predefine", script)) { | 
| 747 |  | std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled(); | 
| 748 |  | if (aatoken.size() != 1) { | 
| 753 |  | std::vector<Token> statement = aatoken[0]; | 
| 754 |  | if (statement.size() > 2) { | 
| 755 |  | int tok = statement[1].tok; | 
| 756 | < | if (tok == Token::identifier || (tok & Token::predefinedset) == Token::predefinedset) { | 
| 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 | < |  | 
| 760 | > |  | 
| 761 |  | } else { | 
| 762 |  | evalError("invalid variable name:" + script); | 
| 763 |  | } | 
| 764 |  | }else { | 
| 765 |  | evalError("bad predefinition length:" + script); | 
| 766 | < | } | 
| 335 | < |  | 
| 766 | > | } | 
| 767 |  |  | 
| 768 |  | } else { | 
| 769 |  | evalError("predefined set compile error:" + script + | 
| 770 |  | "\ncompile error:" + compiler.getErrorMessage()); | 
| 771 |  | } | 
| 341 | – |  | 
| 772 |  | } | 
| 773 |  |  | 
| 774 | < | void SelectionEvaluator::select(BitSet& bs){ | 
| 774 | > | void SelectionEvaluator::select(SelectionSet& bs){ | 
| 775 |  | bs = expression(statement, 1); | 
| 776 |  | } | 
| 777 |  |  | 
| 778 | < | BitSet SelectionEvaluator::lookupValue(const std::string& variable){ | 
| 779 | < |  | 
| 780 | < | BitSet bs(nStuntDouble); | 
| 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(BitSet)) { | 
| 789 | < | return boost::any_cast<BitSet>(i->second); | 
| 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 */ | 
| 795 |  | } else { | 
| 796 |  | unrecognizedIdentifier(variable); | 
| 797 |  | } | 
| 798 | < |  | 
| 798 | > |  | 
| 799 |  | return bs; | 
| 800 |  | } | 
| 801 | < |  | 
| 802 | < | BitSet SelectionEvaluator::nameInstruction(const std::string& name){ | 
| 803 | < |  | 
| 370 | < | return nameFinder.match(name); | 
| 371 | < |  | 
| 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){ | 
| 810 |  | return true; | 
| 811 |  | } | 
| 812 |  | } | 
| 813 | < |  | 
| 813 | > |  | 
| 814 |  | return false; | 
| 815 |  | } | 
| 816 |  |  | 
| 820 |  | //predefine(); | 
| 821 |  | } | 
| 822 |  |  | 
| 823 | < | BitSet SelectionEvaluator::evaluate() { | 
| 824 | < | BitSet bs(nStuntDouble); | 
| 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 | < | BitSet SelectionEvaluator::indexInstruction(const boost::any& value) { | 
| 847 | < | BitSet bs(nStuntDouble); | 
| 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.size()) { | 
| 851 | > | if (index < 0 || index >= bs.bitsets_[STUNTDOUBLE].size()) { | 
| 852 |  | invalidIndex(index); | 
| 853 |  | } else { | 
| 854 |  | bs = indexFinder.find(index); | 
| 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 || indexRange.second >= bs.size()) { | 
| 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); | 
| 867 |  | return bs; | 
| 868 |  | } | 
| 869 |  |  | 
| 870 | + | SelectionSet SelectionEvaluator::allInstruction() { | 
| 871 | + | SelectionSet ss = createSelectionSets(); | 
| 872 |  |  | 
| 873 | < | double SelectionEvaluator::getCharge(Atom* atom) { | 
| 874 | < | double charge =0.0; | 
| 875 | < | AtomType* atomType = atom->getAtomType(); | 
| 876 | < | if (atomType->isCharge()) { | 
| 877 | < | GenericData* data = atomType->getPropertyByName("Charge"); | 
| 878 | < | if (data != NULL) { | 
| 879 | < | DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data); | 
| 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 | < | if (doubleData != NULL) { | 
| 882 | < | charge = doubleData->getData(); | 
| 881 | > | Molecule* mol; | 
| 882 | > | Atom* atom; | 
| 883 | > | RigidBody* rb; | 
| 884 | > | Bond* bond; | 
| 885 | > | Bend* bend; | 
| 886 | > | Torsion* torsion; | 
| 887 | > | Inversion* inversion; | 
| 888 |  |  | 
| 889 | < | } else { | 
| 890 | < | sprintf( painCave.errMsg, | 
| 891 | < | "Can not cast GenericData to DoubleGenericData\n"); | 
| 892 | < | painCave.severity = OOPSE_ERROR; | 
| 893 | < | painCave.isFatal = 1; | 
| 894 | < | simError(); | 
| 895 | < | } | 
| 896 | < | } | 
| 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 |  | } |