| 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) | 
| 53 | – | : info(si), nameFinder(info), distanceFinder(info), indexFinder(info), isLoaded_(false){ | 
| 54 | – |  | 
| 55 | – | nStuntDouble = info->getNGlobalAtoms() + info->getNRigidBodies(); | 
| 56 | – | } | 
| 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 | > | nObjects.push_back(info->getNGlobalMolecules()); | 
| 67 | > | } | 
| 68 | > |  | 
| 69 | > | bool SelectionEvaluator::loadScript(const std::string& filename, | 
| 70 | > | const std::string& script) { | 
| 71 |  | clearDefinitionsAndLoadPredefined(); | 
| 72 |  | this->filename = filename; | 
| 73 |  | this->script = script; | 
| 74 |  | if (! compiler.compile(filename, script)) { | 
| 75 |  | error = true; | 
| 76 |  | errorMessage = compiler.getErrorMessage(); | 
| 77 | < | std::cerr << "SelectionCompiler Error: " << errorMessage << std::endl; | 
| 77 | > |  | 
| 78 | > | sprintf( painCave.errMsg, | 
| 79 | > | "SelectionCompiler Error: %s\n", errorMessage.c_str()); | 
| 80 | > | painCave.severity = OPENMD_ERROR; | 
| 81 | > | painCave.isFatal = 1; | 
| 82 | > | simError(); | 
| 83 |  | return false; | 
| 84 |  | } | 
| 85 |  |  | 
| 117 |  | return loadScriptFileInternal(filename); | 
| 118 |  | } | 
| 119 |  |  | 
| 120 | < | bool SelectionEvaluator::loadScriptFileInternal(const  std::string & filename) { | 
| 121 | < | std::ifstream ifs(filename.c_str()); | 
| 120 | > | bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) { | 
| 121 | > | ifstrstream ifs(filename.c_str()); | 
| 122 |  | if (!ifs.is_open()) { | 
| 123 |  | return false; | 
| 124 |  | } | 
| 125 | < |  | 
| 125 | > |  | 
| 126 |  | const int bufferSize = 65535; | 
| 127 |  | char buffer[bufferSize]; | 
| 128 |  | std::string script; | 
| 131 |  | } | 
| 132 |  | return loadScript(filename, script); | 
| 133 |  | } | 
| 134 | < |  | 
| 135 | < | void SelectionEvaluator::instructionDispatchLoop(OOPSEBitSet& bs){ | 
| 134 | > |  | 
| 135 | > | void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs){ | 
| 136 |  |  | 
| 137 |  | while ( pc < aatoken.size()) { | 
| 138 |  | statement = aatoken[pc++]; | 
| 153 |  |  | 
| 154 |  | } | 
| 155 |  |  | 
| 156 | < | OOPSEBitSet SelectionEvaluator::expression(const std::vector<Token>& code, int pcStart) { | 
| 140 | < | OOPSEBitSet bs; | 
| 141 | < | std::stack<OOPSEBitSet> stack; | 
| 156 | > | void SelectionEvaluator::instructionDispatchLoop(SelectionSet& bs, int frame){ | 
| 157 |  |  | 
| 158 | < | for (int pc = pcStart; pc < code.size(); ++pc) { | 
| 158 | > | while ( pc < aatoken.size()) { | 
| 159 | > | statement = aatoken[pc++]; | 
| 160 | > | statementLength = statement.size(); | 
| 161 | > | Token token = statement[0]; | 
| 162 | > | switch (token.tok) { | 
| 163 | > | case Token::define: | 
| 164 | > | define(); | 
| 165 | > | break; | 
| 166 | > | case Token::select: | 
| 167 | > | select(bs, frame); | 
| 168 | > | break; | 
| 169 | > | default: | 
| 170 | > | unrecognizedCommand(token); | 
| 171 | > | return; | 
| 172 | > | } | 
| 173 | > | } | 
| 174 | > |  | 
| 175 | > | } | 
| 176 | > |  | 
| 177 | > | SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code, | 
| 178 | > | int pcStart) { | 
| 179 | > | SelectionSet bs = createSelectionSets(); | 
| 180 | > | std::stack<SelectionSet> stack; | 
| 181 | > | vector<int> bsSize = bs.size(); | 
| 182 | > |  | 
| 183 | > | for (unsigned int pc = pcStart; pc < code.size(); ++pc) { | 
| 184 |  | Token instruction = code[pc]; | 
| 185 |  |  | 
| 186 |  | switch (instruction.tok) { | 
| 189 |  | case Token::expressionEnd: | 
| 190 |  | break; | 
| 191 |  | case Token::all: | 
| 192 | < | bs = OOPSEBitSet(nStuntDouble); | 
| 193 | < | bs.setAll(); | 
| 192 | > | bs = allInstruction(); | 
| 193 | > | stack.push(bs); | 
| 194 | > | break; | 
| 195 | > | case Token::none: | 
| 196 | > | bs = createSelectionSets(); | 
| 197 | > | stack.push(bs); | 
| 198 | > | break; | 
| 199 | > | case Token::opOr: | 
| 200 | > | bs= stack.top(); | 
| 201 | > | stack.pop(); | 
| 202 | > | stack.top() |= bs; | 
| 203 | > | break; | 
| 204 | > | case Token::opAnd: | 
| 205 | > | bs = stack.top(); | 
| 206 | > | stack.pop(); | 
| 207 | > | stack.top() &= bs; | 
| 208 | > | break; | 
| 209 | > | case Token::opNot: | 
| 210 | > | stack.top().flip(); | 
| 211 | > | break; | 
| 212 | > | case Token::within: | 
| 213 | > | withinInstruction(instruction, stack.top()); | 
| 214 | > | break; | 
| 215 | > | case Token::hull: | 
| 216 | > | stack.push(hull()); | 
| 217 | > | break; | 
| 218 | > | //case Token::selected: | 
| 219 | > | //  stack.push(getSelectionSet()); | 
| 220 | > | //  break; | 
| 221 | > | case Token::name: | 
| 222 | > | stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value))); | 
| 223 | > | break; | 
| 224 | > | case Token::index: | 
| 225 | > | stack.push(indexInstruction(instruction.value)); | 
| 226 | > | break; | 
| 227 | > | case Token::identifier: | 
| 228 | > | stack.push(lookupValue(boost::any_cast<std::string>(instruction.value))); | 
| 229 | > | break; | 
| 230 | > | case Token::opLT: | 
| 231 | > | case Token::opLE: | 
| 232 | > | case Token::opGE: | 
| 233 | > | case Token::opGT: | 
| 234 | > | case Token::opEQ: | 
| 235 | > | case Token::opNE: | 
| 236 | > | stack.push(comparatorInstruction(instruction)); | 
| 237 | > | break; | 
| 238 | > | default: | 
| 239 | > | unrecognizedExpression(); | 
| 240 | > | } | 
| 241 | > | } | 
| 242 | > | if (stack.size() != 1) | 
| 243 | > | evalError("atom expression compiler error - stack over/underflow"); | 
| 244 | > |  | 
| 245 | > | return stack.top(); | 
| 246 | > | } | 
| 247 | > |  | 
| 248 | > |  | 
| 249 | > | SelectionSet SelectionEvaluator::expression(const std::vector<Token>& code, | 
| 250 | > | int pcStart, int frame) { | 
| 251 | > | SelectionSet bs = createSelectionSets(); | 
| 252 | > | std::stack<SelectionSet> stack; | 
| 253 | > |  | 
| 254 | > | for (unsigned int pc = pcStart; pc < code.size(); ++pc) { | 
| 255 | > | Token instruction = code[pc]; | 
| 256 | > |  | 
| 257 | > | switch (instruction.tok) { | 
| 258 | > | case Token::expressionBegin: | 
| 259 | > | break; | 
| 260 | > | case Token::expressionEnd: | 
| 261 | > | break; | 
| 262 | > | case Token::all: | 
| 263 | > | bs = allInstruction(); | 
| 264 |  | stack.push(bs); | 
| 265 |  | break; | 
| 266 |  | case Token::none: | 
| 267 | < | bs = OOPSEBitSet(nStuntDouble); | 
| 267 | > | bs = SelectionSet(nObjects); | 
| 268 |  | stack.push(bs); | 
| 269 |  | break; | 
| 270 |  | case Token::opOr: | 
| 281 |  | stack.top().flip(); | 
| 282 |  | break; | 
| 283 |  | case Token::within: | 
| 284 | < |  | 
| 175 | < | withinInstruction(instruction, stack.top()); | 
| 284 | > | withinInstruction(instruction, stack.top(), frame); | 
| 285 |  | break; | 
| 286 | + | case Token::hull: | 
| 287 | + | stack.push(hull(frame)); | 
| 288 | + | break; | 
| 289 |  | //case Token::selected: | 
| 290 |  | //  stack.push(getSelectionSet()); | 
| 291 |  | //  break; | 
| 304 |  | case Token::opGT: | 
| 305 |  | case Token::opEQ: | 
| 306 |  | case Token::opNE: | 
| 307 | < | stack.push(comparatorInstruction(instruction)); | 
| 307 | > | stack.push(comparatorInstruction(instruction, frame)); | 
| 308 |  | break; | 
| 309 |  | default: | 
| 310 |  | unrecognizedExpression(); | 
| 318 |  |  | 
| 319 |  |  | 
| 320 |  |  | 
| 321 | < | OOPSEBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { | 
| 321 | > | SelectionSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { | 
| 322 |  | int comparator = instruction.tok; | 
| 323 |  | int property = instruction.intValue; | 
| 324 |  | float comparisonValue = boost::any_cast<float>(instruction.value); | 
| 325 | < | float propertyValue; | 
| 214 | < | OOPSEBitSet bs(nStuntDouble); | 
| 325 | > | SelectionSet bs = createSelectionSets(); | 
| 326 |  | bs.clearAll(); | 
| 327 |  |  | 
| 328 |  | SimInfo::MoleculeIterator mi; | 
| 331 |  | Atom* atom; | 
| 332 |  | Molecule::RigidBodyIterator rbIter; | 
| 333 |  | RigidBody* rb; | 
| 223 | – |  | 
| 224 | – | for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) { | 
| 334 |  |  | 
| 335 | + | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 336 | + | mol = info->nextMolecule(mi)) { | 
| 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, OOPSEBitSet& bs, int property, int comparator, float comparisonValue) { | 
| 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 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 369 | > | compareProperty(atom, bs, property, comparator, comparisonValue, frame); | 
| 370 | > | } | 
| 371 | > |  | 
| 372 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 373 | > | rb = mol->nextRigidBody(rbIter)) { | 
| 374 | > | compareProperty(rb, bs, property, comparator, comparisonValue, frame); | 
| 375 | > | } | 
| 376 | > | } | 
| 377 | > |  | 
| 378 | > | return bs; | 
| 379 | > | } | 
| 380 | > |  | 
| 381 | > | void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs, | 
| 382 | > | int property, int comparator, | 
| 383 | > | float comparisonValue) { | 
| 384 |  | RealType propertyValue = 0.0; | 
| 385 | + | Vector3d pos; | 
| 386 | + |  | 
| 387 |  | switch (property) { | 
| 388 |  | case Token::mass: | 
| 389 |  | propertyValue = sd->getMass(); | 
| 401 |  | } | 
| 402 |  | } | 
| 403 |  | break; | 
| 404 | + | case Token::x: | 
| 405 | + | propertyValue = sd->getPos().x(); | 
| 406 | + | break; | 
| 407 | + | case Token::y: | 
| 408 | + | propertyValue = sd->getPos().y(); | 
| 409 | + | break; | 
| 410 | + | case Token::z: | 
| 411 | + | propertyValue = sd->getPos().z(); | 
| 412 | + | break; | 
| 413 | + | case Token::wrappedX: | 
| 414 | + | pos = sd->getPos(); | 
| 415 | + | info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); | 
| 416 | + | propertyValue = pos.x(); | 
| 417 | + | break; | 
| 418 | + | case Token::wrappedY: | 
| 419 | + | pos = sd->getPos(); | 
| 420 | + | info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); | 
| 421 | + | propertyValue = pos.y(); | 
| 422 | + | break; | 
| 423 | + | case Token::wrappedZ: | 
| 424 | + | pos = sd->getPos(); | 
| 425 | + | info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); | 
| 426 | + | propertyValue = pos.z(); | 
| 427 | + | break; | 
| 428 | + | case Token::r: | 
| 429 | + | propertyValue = sd->getPos().length(); | 
| 430 | + | break; | 
| 431 |  | default: | 
| 432 |  | unrecognizedAtomProperty(property); | 
| 433 |  | } | 
| 453 |  | match = propertyValue != comparisonValue; | 
| 454 |  | break; | 
| 455 |  | } | 
| 282 | – | if (match) | 
| 283 | – | bs.setBitOn(sd->getGlobalIndex()); | 
| 456 |  |  | 
| 457 | + | if (match) | 
| 458 | + | bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex()); | 
| 459 | + |  | 
| 460 | + |  | 
| 461 |  | } | 
| 462 |  |  | 
| 463 | < | void SelectionEvaluator::withinInstruction(const Token& instruction, OOPSEBitSet& bs){ | 
| 463 | > | void SelectionEvaluator::compareProperty(StuntDouble* sd, SelectionSet& bs, | 
| 464 | > | int property, int comparator, | 
| 465 | > | float comparisonValue, int frame) { | 
| 466 | > | RealType propertyValue = 0.0; | 
| 467 | > | Vector3d pos; | 
| 468 | > | switch (property) { | 
| 469 | > | case Token::mass: | 
| 470 | > | propertyValue = sd->getMass(); | 
| 471 | > | break; | 
| 472 | > | case Token::charge: | 
| 473 | > | if (sd->isAtom()){ | 
| 474 | > | Atom* atom = static_cast<Atom*>(sd); | 
| 475 | > | propertyValue = getCharge(atom,frame); | 
| 476 | > | } else if (sd->isRigidBody()) { | 
| 477 | > | RigidBody* rb = static_cast<RigidBody*>(sd); | 
| 478 | > | RigidBody::AtomIterator ai; | 
| 479 | > | Atom* atom; | 
| 480 | > | for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) { | 
| 481 | > | propertyValue+=  getCharge(atom,frame); | 
| 482 | > | } | 
| 483 | > | } | 
| 484 | > | break; | 
| 485 | > | case Token::x: | 
| 486 | > | propertyValue = sd->getPos(frame).x(); | 
| 487 | > | break; | 
| 488 | > | case Token::y: | 
| 489 | > | propertyValue = sd->getPos(frame).y(); | 
| 490 | > | break; | 
| 491 | > | case Token::z: | 
| 492 | > | propertyValue = sd->getPos(frame).z(); | 
| 493 | > | break; | 
| 494 | > | case Token::wrappedX: | 
| 495 | > | pos = sd->getPos(frame); | 
| 496 | > | info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); | 
| 497 | > | propertyValue = pos.x(); | 
| 498 | > | break; | 
| 499 | > | case Token::wrappedY: | 
| 500 | > | pos = sd->getPos(frame); | 
| 501 | > | info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); | 
| 502 | > | propertyValue = pos.y(); | 
| 503 | > | break; | 
| 504 | > | case Token::wrappedZ: | 
| 505 | > | pos = sd->getPos(frame); | 
| 506 | > | info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); | 
| 507 | > | propertyValue = pos.z(); | 
| 508 | > | break; | 
| 509 | > |  | 
| 510 | > | case Token::r: | 
| 511 | > | propertyValue = sd->getPos(frame).length(); | 
| 512 | > | break; | 
| 513 | > | default: | 
| 514 | > | unrecognizedAtomProperty(property); | 
| 515 | > | } | 
| 516 | > |  | 
| 517 | > | bool match = false; | 
| 518 | > | switch (comparator) { | 
| 519 | > | case Token::opLT: | 
| 520 | > | match = propertyValue < comparisonValue; | 
| 521 | > | break; | 
| 522 | > | case Token::opLE: | 
| 523 | > | match = propertyValue <= comparisonValue; | 
| 524 | > | break; | 
| 525 | > | case Token::opGE: | 
| 526 | > | match = propertyValue >= comparisonValue; | 
| 527 | > | break; | 
| 528 | > | case Token::opGT: | 
| 529 | > | match = propertyValue > comparisonValue; | 
| 530 | > | break; | 
| 531 | > | case Token::opEQ: | 
| 532 | > | match = propertyValue == comparisonValue; | 
| 533 | > | break; | 
| 534 | > | case Token::opNE: | 
| 535 | > | match = propertyValue != comparisonValue; | 
| 536 | > | break; | 
| 537 | > | } | 
| 538 | > | if (match) | 
| 539 | > | bs.bitsets_[STUNTDOUBLE].setBitOn(sd->getGlobalIndex()); | 
| 540 | > |  | 
| 541 |  |  | 
| 542 | + | } | 
| 543 | + |  | 
| 544 | + | void SelectionEvaluator::withinInstruction(const Token& instruction, | 
| 545 | + | SelectionSet& bs){ | 
| 546 | + |  | 
| 547 |  | boost::any withinSpec = instruction.value; | 
| 548 |  | float distance; | 
| 549 |  | if (withinSpec.type() == typeid(float)){ | 
| 558 |  | bs = distanceFinder.find(bs, distance); | 
| 559 |  | } | 
| 560 |  |  | 
| 561 | + | void SelectionEvaluator::withinInstruction(const Token& instruction, | 
| 562 | + | SelectionSet& bs, int frame){ | 
| 563 | + |  | 
| 564 | + | boost::any withinSpec = instruction.value; | 
| 565 | + | float distance; | 
| 566 | + | if (withinSpec.type() == typeid(float)){ | 
| 567 | + | distance = boost::any_cast<float>(withinSpec); | 
| 568 | + | } else if (withinSpec.type() == typeid(int)) { | 
| 569 | + | distance = boost::any_cast<int>(withinSpec); | 
| 570 | + | } else { | 
| 571 | + | evalError("casting error in withinInstruction"); | 
| 572 | + | bs.clearAll(); | 
| 573 | + | } | 
| 574 | + |  | 
| 575 | + | bs = distanceFinder.find(bs, distance, frame); | 
| 576 | + | } | 
| 577 | + |  | 
| 578 |  | void SelectionEvaluator::define() { | 
| 579 |  | assert(statement.size() >= 3); | 
| 580 | < |  | 
| 580 | > |  | 
| 581 |  | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 582 | < |  | 
| 583 | < | variables.insert(VariablesType::value_type(variable, expression(statement, 2))); | 
| 582 | > |  | 
| 583 | > | variables.insert(VariablesType::value_type(variable, | 
| 584 | > | expression(statement, 2))); | 
| 585 |  | } | 
| 586 | + |  | 
| 587 |  |  | 
| 311 | – |  | 
| 588 |  | /** @todo */ | 
| 589 |  | void SelectionEvaluator::predefine(const std::string& script) { | 
| 590 | < |  | 
| 590 | > |  | 
| 591 |  | if (compiler.compile("#predefine", script)) { | 
| 592 |  | std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled(); | 
| 593 |  | if (aatoken.size() != 1) { | 
| 598 |  | std::vector<Token> statement = aatoken[0]; | 
| 599 |  | if (statement.size() > 2) { | 
| 600 |  | int tok = statement[1].tok; | 
| 601 | < | if (tok == Token::identifier || (tok & Token::predefinedset) == Token::predefinedset) { | 
| 601 | > | if (tok == Token::identifier || | 
| 602 | > | (tok & Token::predefinedset) == Token::predefinedset) { | 
| 603 |  | std::string variable = boost::any_cast<std::string>(statement[1].value); | 
| 604 |  | variables.insert(VariablesType::value_type(variable, statement)); | 
| 605 | < |  | 
| 605 | > |  | 
| 606 |  | } else { | 
| 607 |  | evalError("invalid variable name:" + script); | 
| 608 |  | } | 
| 609 |  | }else { | 
| 610 |  | evalError("bad predefinition length:" + script); | 
| 611 | < | } | 
| 335 | < |  | 
| 611 | > | } | 
| 612 |  |  | 
| 613 |  | } else { | 
| 614 |  | evalError("predefined set compile error:" + script + | 
| 615 |  | "\ncompile error:" + compiler.getErrorMessage()); | 
| 616 |  | } | 
| 341 | – |  | 
| 617 |  | } | 
| 618 |  |  | 
| 619 | < | void SelectionEvaluator::select(OOPSEBitSet& bs){ | 
| 619 | > | void SelectionEvaluator::select(SelectionSet& bs){ | 
| 620 |  | bs = expression(statement, 1); | 
| 621 |  | } | 
| 622 |  |  | 
| 623 | < | OOPSEBitSet SelectionEvaluator::lookupValue(const std::string& variable){ | 
| 624 | < |  | 
| 625 | < | OOPSEBitSet bs(nStuntDouble); | 
| 623 | > | void SelectionEvaluator::select(SelectionSet& bs, int frame){ | 
| 624 | > | bs = expression(statement, 1, frame); | 
| 625 | > | } | 
| 626 | > |  | 
| 627 | > | SelectionSet SelectionEvaluator::lookupValue(const std::string& variable){ | 
| 628 | > |  | 
| 629 | > | SelectionSet bs = createSelectionSets(); | 
| 630 |  | std::map<std::string, boost::any>::iterator i = variables.find(variable); | 
| 631 |  |  | 
| 632 |  | if (i != variables.end()) { | 
| 633 | < | if (i->second.type() == typeid(OOPSEBitSet)) { | 
| 634 | < | return boost::any_cast<OOPSEBitSet>(i->second); | 
| 633 | > | if (i->second.type() == typeid(SelectionSet)) { | 
| 634 | > | return boost::any_cast<SelectionSet>(i->second); | 
| 635 |  | } else if (i->second.type() ==  typeid(std::vector<Token>)){ | 
| 636 |  | bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2); | 
| 637 |  | i->second =  bs; /**@todo fixme */ | 
| 640 |  | } else { | 
| 641 |  | unrecognizedIdentifier(variable); | 
| 642 |  | } | 
| 643 | < |  | 
| 643 | > |  | 
| 644 |  | return bs; | 
| 645 |  | } | 
| 646 | < |  | 
| 647 | < | OOPSEBitSet SelectionEvaluator::nameInstruction(const std::string& name){ | 
| 648 | < |  | 
| 370 | < | return nameFinder.match(name); | 
| 371 | < |  | 
| 646 | > |  | 
| 647 | > | SelectionSet SelectionEvaluator::nameInstruction(const std::string& name){ | 
| 648 | > | return nameFinder.match(name); | 
| 649 |  | } | 
| 650 |  |  | 
| 651 |  | bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){ | 
| 655 |  | return true; | 
| 656 |  | } | 
| 657 |  | } | 
| 658 | < |  | 
| 658 | > |  | 
| 659 |  | return false; | 
| 660 |  | } | 
| 661 |  |  | 
| 665 |  | //predefine(); | 
| 666 |  | } | 
| 667 |  |  | 
| 668 | < | OOPSEBitSet SelectionEvaluator::evaluate() { | 
| 669 | < | OOPSEBitSet bs(nStuntDouble); | 
| 668 | > | SelectionSet SelectionEvaluator::createSelectionSets() { | 
| 669 | > | SelectionSet ss(nObjects); | 
| 670 | > | return ss; | 
| 671 | > | } | 
| 672 | > |  | 
| 673 | > | SelectionSet SelectionEvaluator::evaluate() { | 
| 674 | > | SelectionSet bs = createSelectionSets(); | 
| 675 |  | if (isLoaded_) { | 
| 676 |  | pc = 0; | 
| 677 |  | instructionDispatchLoop(bs); | 
| 678 |  | } | 
| 679 | + | return bs; | 
| 680 | + | } | 
| 681 |  |  | 
| 682 | + | SelectionSet SelectionEvaluator::evaluate(int frame) { | 
| 683 | + | SelectionSet bs = createSelectionSets(); | 
| 684 | + | if (isLoaded_) { | 
| 685 | + | pc = 0; | 
| 686 | + | instructionDispatchLoop(bs, frame); | 
| 687 | + | } | 
| 688 |  | return bs; | 
| 689 |  | } | 
| 690 |  |  | 
| 691 | < | OOPSEBitSet SelectionEvaluator::indexInstruction(const boost::any& value) { | 
| 692 | < | OOPSEBitSet bs(nStuntDouble); | 
| 691 | > | SelectionSet SelectionEvaluator::indexInstruction(const boost::any& value) { | 
| 692 | > | SelectionSet bs = createSelectionSets(); | 
| 693 |  |  | 
| 694 |  | if (value.type() == typeid(int)) { | 
| 695 |  | int index = boost::any_cast<int>(value); | 
| 696 | < | if (index < 0 || index >= bs.size()) { | 
| 696 | > | if (index < 0 || index >= bs.bitsets_[STUNTDOUBLE].size()) { | 
| 697 |  | invalidIndex(index); | 
| 698 |  | } else { | 
| 699 |  | bs = indexFinder.find(index); | 
| 701 |  | } else if (value.type() == typeid(std::pair<int, int>)) { | 
| 702 |  | std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value); | 
| 703 |  | assert(indexRange.first <= indexRange.second); | 
| 704 | < | if (indexRange.first < 0 || indexRange.second >= bs.size()) { | 
| 704 | > | if (indexRange.first < 0 || | 
| 705 | > | indexRange.second >= bs.bitsets_[STUNTDOUBLE].size()) { | 
| 706 |  | invalidIndexRange(indexRange); | 
| 707 |  | }else { | 
| 708 |  | bs = indexFinder.find(indexRange.first, indexRange.second); | 
| 712 |  | return bs; | 
| 713 |  | } | 
| 714 |  |  | 
| 715 | + | SelectionSet SelectionEvaluator::allInstruction() { | 
| 716 | + | SelectionSet ss = createSelectionSets(); | 
| 717 |  |  | 
| 718 | + | SimInfo::MoleculeIterator mi; | 
| 719 | + | Molecule::AtomIterator ai; | 
| 720 | + | Molecule::RigidBodyIterator rbIter; | 
| 721 | + | Molecule::BondIterator bondIter; | 
| 722 | + | Molecule::BendIterator bendIter; | 
| 723 | + | Molecule::TorsionIterator torsionIter; | 
| 724 | + | Molecule::InversionIterator inversionIter; | 
| 725 | + |  | 
| 726 | + | Molecule* mol; | 
| 727 | + | Atom* atom; | 
| 728 | + | RigidBody* rb; | 
| 729 | + | Bond* bond; | 
| 730 | + | Bend* bend; | 
| 731 | + | Torsion* torsion; | 
| 732 | + | Inversion* inversion; | 
| 733 | + |  | 
| 734 | + | // Doing the loop insures that we're actually on this processor. | 
| 735 | + |  | 
| 736 | + | for (mol = info->beginMolecule(mi); mol != NULL; | 
| 737 | + | mol = info->nextMolecule(mi)) { | 
| 738 | + | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 739 | + | ss.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex()); | 
| 740 | + | } | 
| 741 | + | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 742 | + | rb = mol->nextRigidBody(rbIter)) { | 
| 743 | + | ss.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex()); | 
| 744 | + | } | 
| 745 | + | for (bond = mol->beginBond(bondIter); bond != NULL; | 
| 746 | + | bond = mol->nextBond(bondIter)) { | 
| 747 | + | ss.bitsets_[BOND].setBitOn(bond->getGlobalIndex()); | 
| 748 | + | } | 
| 749 | + | for (bend = mol->beginBend(bendIter); bend != NULL; | 
| 750 | + | bend = mol->nextBend(bendIter)) { | 
| 751 | + | ss.bitsets_[BEND].setBitOn(bend->getGlobalIndex()); | 
| 752 | + | } | 
| 753 | + | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; | 
| 754 | + | torsion = mol->nextTorsion(torsionIter)) { | 
| 755 | + | ss.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex()); | 
| 756 | + | } | 
| 757 | + | for (inversion = mol->beginInversion(inversionIter); inversion != NULL; | 
| 758 | + | inversion = mol->nextInversion(inversionIter)) { | 
| 759 | + | ss.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex()); | 
| 760 | + | } | 
| 761 | + | ss.bitsets_[MOLECULE].setBitOn(mol->getGlobalIndex()); | 
| 762 | + | } | 
| 763 | + |  | 
| 764 | + | return ss; | 
| 765 | + | } | 
| 766 | + |  | 
| 767 | + | SelectionSet SelectionEvaluator::hull() { | 
| 768 | + | SelectionSet bs = createSelectionSets(); | 
| 769 | + |  | 
| 770 | + | bs = hullFinder.findHull(); | 
| 771 | + | surfaceArea_ = hullFinder.getSurfaceArea(); | 
| 772 | + | hasSurfaceArea_ = true; | 
| 773 | + | return bs; | 
| 774 | + | } | 
| 775 | + |  | 
| 776 | + |  | 
| 777 | + | SelectionSet SelectionEvaluator::hull(int frame) { | 
| 778 | + | SelectionSet bs = createSelectionSets(); | 
| 779 | + |  | 
| 780 | + | bs = hullFinder.findHull(frame); | 
| 781 | + |  | 
| 782 | + | return bs; | 
| 783 | + | } | 
| 784 | + |  | 
| 785 |  | RealType SelectionEvaluator::getCharge(Atom* atom) { | 
| 786 | < | RealType charge =0.0; | 
| 786 | > | RealType charge = 0.0; | 
| 787 |  | AtomType* atomType = atom->getAtomType(); | 
| 428 | – | if (atomType->isCharge()) { | 
| 429 | – | GenericData* data = atomType->getPropertyByName("Charge"); | 
| 430 | – | if (data != NULL) { | 
| 431 | – | DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data); | 
| 788 |  |  | 
| 789 | < | if (doubleData != NULL) { | 
| 790 | < | charge = doubleData->getData(); | 
| 791 | < |  | 
| 436 | < | } else { | 
| 437 | < | sprintf( painCave.errMsg, | 
| 438 | < | "Can not cast GenericData to DoubleGenericData\n"); | 
| 439 | < | painCave.severity = OOPSE_ERROR; | 
| 440 | < | painCave.isFatal = 1; | 
| 441 | < | simError(); | 
| 442 | < | } | 
| 443 | < | } | 
| 789 | > | FixedChargeAdapter fca = FixedChargeAdapter(atomType); | 
| 790 | > | if ( fca.isFixedCharge() ) { | 
| 791 | > | charge = fca.getCharge(); | 
| 792 |  | } | 
| 793 | + |  | 
| 794 | + | FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); | 
| 795 | + | if ( fqa.isFluctuatingCharge() ) { | 
| 796 | + | charge += atom->getFlucQPos(); | 
| 797 | + | } | 
| 798 | + | return charge; | 
| 799 | + | } | 
| 800 |  |  | 
| 801 | + | RealType SelectionEvaluator::getCharge(Atom* atom, int frame) { | 
| 802 | + | RealType charge = 0.0; | 
| 803 | + | AtomType* atomType = atom->getAtomType(); | 
| 804 | + |  | 
| 805 | + | FixedChargeAdapter fca = FixedChargeAdapter(atomType); | 
| 806 | + | if ( fca.isFixedCharge() ) { | 
| 807 | + | charge = fca.getCharge(); | 
| 808 | + | } | 
| 809 | + |  | 
| 810 | + | FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); | 
| 811 | + | if ( fqa.isFluctuatingCharge() ) { | 
| 812 | + | charge += atom->getFlucQPos(frame); | 
| 813 | + | } | 
| 814 |  | return charge; | 
| 815 |  | } | 
| 816 |  |  |