--- trunk/src/selection/SelectionEvaluator.cpp 2010/03/22 18:45:39 1412 +++ trunk/src/selection/SelectionEvaluator.cpp 2013/07/16 18:58:08 1903 @@ -35,8 +35,9 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). - * [4] Vardeman & Gezelter, in progress (2009). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ #include @@ -45,14 +46,17 @@ #include "primitives/DirectionalAtom.hpp" #include "primitives/RigidBody.hpp" #include "primitives/Molecule.hpp" -#include "io/basic_ifstrstream.hpp" +#include "io/ifstrstream.hpp" +#include "types/FixedChargeAdapter.hpp" +#include "types/FluctuatingChargeAdapter.hpp" + namespace OpenMD { SelectionEvaluator::SelectionEvaluator(SimInfo* si) : info(si), nameFinder(info), distanceFinder(info), hullFinder(info), - indexFinder(info), + indexFinder(info), hasSurfaceArea_(false), isLoaded_(false){ nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies(); } @@ -144,12 +148,33 @@ namespace OpenMD { } + void SelectionEvaluator::instructionDispatchLoop(OpenMDBitSet& bs, int frame){ + + while ( pc < aatoken.size()) { + statement = aatoken[pc++]; + statementLength = statement.size(); + Token token = statement[0]; + switch (token.tok) { + case Token::define: + define(); + break; + case Token::select: + select(bs, frame); + break; + default: + unrecognizedCommand(token); + return; + } + } + + } + OpenMDBitSet SelectionEvaluator::expression(const std::vector& code, int pcStart) { OpenMDBitSet bs; std::stack stack; - for (int pc = pcStart; pc < code.size(); ++pc) { + for (unsigned int pc = pcStart; pc < code.size(); ++pc) { Token instruction = code[pc]; switch (instruction.tok) { @@ -158,8 +183,7 @@ namespace OpenMD { case Token::expressionEnd: break; case Token::all: - bs = OpenMDBitSet(nStuntDouble); - bs.setAll(); + bs = allInstruction(); stack.push(bs); break; case Token::none: @@ -216,12 +240,82 @@ namespace OpenMD { } + OpenMDBitSet SelectionEvaluator::expression(const std::vector& code, + int pcStart, int frame) { + OpenMDBitSet bs; + std::stack stack; + + for (unsigned int pc = pcStart; pc < code.size(); ++pc) { + Token instruction = code[pc]; + switch (instruction.tok) { + case Token::expressionBegin: + break; + case Token::expressionEnd: + break; + case Token::all: + bs = allInstruction(); + stack.push(bs); + break; + case Token::none: + bs = OpenMDBitSet(nStuntDouble); + stack.push(bs); + break; + case Token::opOr: + bs = stack.top(); + stack.pop(); + stack.top() |= bs; + break; + case Token::opAnd: + bs = stack.top(); + stack.pop(); + stack.top() &= bs; + break; + case Token::opNot: + stack.top().flip(); + break; + case Token::within: + withinInstruction(instruction, stack.top(), frame); + break; + case Token::hull: + stack.push(hull(frame)); + break; + //case Token::selected: + // stack.push(getSelectionSet()); + // break; + case Token::name: + stack.push(nameInstruction(boost::any_cast(instruction.value))); + break; + case Token::index: + stack.push(indexInstruction(instruction.value)); + break; + case Token::identifier: + stack.push(lookupValue(boost::any_cast(instruction.value))); + break; + case Token::opLT: + case Token::opLE: + case Token::opGE: + case Token::opGT: + case Token::opEQ: + case Token::opNE: + stack.push(comparatorInstruction(instruction, frame)); + break; + default: + unrecognizedExpression(); + } + } + if (stack.size() != 1) + evalError("atom expression compiler error - stack over/underflow"); + + return stack.top(); + } + + + OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) { int comparator = instruction.tok; int property = instruction.intValue; float comparisonValue = boost::any_cast(instruction.value); - float propertyValue; OpenMDBitSet bs(nStuntDouble); bs.clearAll(); @@ -248,10 +342,41 @@ namespace OpenMD { return bs; } + OpenMDBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction, int frame) { + int comparator = instruction.tok; + int property = instruction.intValue; + float comparisonValue = boost::any_cast(instruction.value); + OpenMDBitSet bs(nStuntDouble); + bs.clearAll(); + + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + + for (mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { + + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + compareProperty(atom, bs, property, comparator, comparisonValue, frame); + } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + compareProperty(rb, bs, property, comparator, comparisonValue, frame); + } + } + + return bs; + } + void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs, int property, int comparator, float comparisonValue) { RealType propertyValue = 0.0; + Vector3d pos; switch (property) { case Token::mass: propertyValue = sd->getMass(); @@ -278,6 +403,24 @@ namespace OpenMD { case Token::z: propertyValue = sd->getPos().z(); break; + case Token::wrappedX: + pos = sd->getPos(); + info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); + propertyValue = pos.x(); + break; + case Token::wrappedY: + pos = sd->getPos(); + info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); + propertyValue = pos.y(); + break; + case Token::wrappedZ: + pos = sd->getPos(); + info->getSnapshotManager()->getCurrentSnapshot()->wrapVector(pos); + propertyValue = pos.z(); + break; + case Token::r: + propertyValue = sd->getPos().length(); + break; default: unrecognizedAtomProperty(property); } @@ -303,11 +446,93 @@ namespace OpenMD { match = propertyValue != comparisonValue; break; } - if (match) + if (match) bs.setBitOn(sd->getGlobalIndex()); + } + void SelectionEvaluator::compareProperty(StuntDouble* sd, OpenMDBitSet& bs, + int property, int comparator, + float comparisonValue, int frame) { + RealType propertyValue = 0.0; + Vector3d pos; + switch (property) { + case Token::mass: + propertyValue = sd->getMass(); + break; + case Token::charge: + if (sd->isAtom()){ + Atom* atom = static_cast(sd); + propertyValue = getCharge(atom,frame); + } else if (sd->isRigidBody()) { + RigidBody* rb = static_cast(sd); + RigidBody::AtomIterator ai; + Atom* atom; + for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) { + propertyValue+= getCharge(atom,frame); + } + } + break; + case Token::x: + propertyValue = sd->getPos(frame).x(); + break; + case Token::y: + propertyValue = sd->getPos(frame).y(); + break; + case Token::z: + propertyValue = sd->getPos(frame).z(); + break; + case Token::wrappedX: + pos = sd->getPos(frame); + info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); + propertyValue = pos.x(); + break; + case Token::wrappedY: + pos = sd->getPos(frame); + info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); + propertyValue = pos.y(); + break; + case Token::wrappedZ: + pos = sd->getPos(frame); + info->getSnapshotManager()->getSnapshot(frame)->wrapVector(pos); + propertyValue = pos.z(); + break; + + case Token::r: + propertyValue = sd->getPos(frame).length(); + break; + default: + unrecognizedAtomProperty(property); + } + + bool match = false; + switch (comparator) { + case Token::opLT: + match = propertyValue < comparisonValue; + break; + case Token::opLE: + match = propertyValue <= comparisonValue; + break; + case Token::opGE: + match = propertyValue >= comparisonValue; + break; + case Token::opGT: + match = propertyValue > comparisonValue; + break; + case Token::opEQ: + match = propertyValue == comparisonValue; + break; + case Token::opNE: + match = propertyValue != comparisonValue; + break; + } + if (match) + bs.setBitOn(sd->getGlobalIndex()); + + + } + void SelectionEvaluator::withinInstruction(const Token& instruction, OpenMDBitSet& bs){ @@ -324,6 +549,23 @@ namespace OpenMD { bs = distanceFinder.find(bs, distance); } + + void SelectionEvaluator::withinInstruction(const Token& instruction, + OpenMDBitSet& bs, int frame){ + + boost::any withinSpec = instruction.value; + float distance; + if (withinSpec.type() == typeid(float)){ + distance = boost::any_cast(withinSpec); + } else if (withinSpec.type() == typeid(int)) { + distance = boost::any_cast(withinSpec); + } else { + evalError("casting error in withinInstruction"); + bs.clearAll(); + } + + bs = distanceFinder.find(bs, distance, frame); + } void SelectionEvaluator::define() { assert(statement.size() >= 3); @@ -369,6 +611,10 @@ namespace OpenMD { void SelectionEvaluator::select(OpenMDBitSet& bs){ bs = expression(statement, 1); } + + void SelectionEvaluator::select(OpenMDBitSet& bs, int frame){ + bs = expression(statement, 1, frame); + } OpenMDBitSet SelectionEvaluator::lookupValue(const std::string& variable){ @@ -417,7 +663,15 @@ namespace OpenMD { pc = 0; instructionDispatchLoop(bs); } + return bs; + } + OpenMDBitSet SelectionEvaluator::evaluate(int frame) { + OpenMDBitSet bs(nStuntDouble); + if (isLoaded_) { + pc = 0; + instructionDispatchLoop(bs, frame); + } return bs; } @@ -444,38 +698,81 @@ namespace OpenMD { return bs; } + OpenMDBitSet SelectionEvaluator::allInstruction() { + OpenMDBitSet bs(nStuntDouble); + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + + // Doing the loop insures that we're actually on this processor. + + for (mol = info->beginMolecule(mi); mol != NULL; + mol = info->nextMolecule(mi)) { + + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + bs.setBitOn(atom->getGlobalIndex()); + } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + bs.setBitOn(rb->getGlobalIndex()); + } + } + + return bs; + } + OpenMDBitSet SelectionEvaluator::hull() { OpenMDBitSet bs(nStuntDouble); bs = hullFinder.findHull(); - + surfaceArea_ = hullFinder.getSurfaceArea(); + hasSurfaceArea_ = true; return bs; } + OpenMDBitSet SelectionEvaluator::hull(int frame) { + OpenMDBitSet bs(nStuntDouble); + + bs = hullFinder.findHull(frame); + return bs; + } + RealType SelectionEvaluator::getCharge(Atom* atom) { - RealType charge =0.0; + RealType charge = 0.0; AtomType* atomType = atom->getAtomType(); - if (atomType->isCharge()) { - GenericData* data = atomType->getPropertyByName("Charge"); - if (data != NULL) { - DoubleGenericData* doubleData= dynamic_cast(data); - if (doubleData != NULL) { - charge = doubleData->getData(); - - } else { - sprintf( painCave.errMsg, - "Can not cast GenericData to DoubleGenericData\n"); - painCave.severity = OPENMD_ERROR; - painCave.isFatal = 1; - simError(); - } - } + FixedChargeAdapter fca = FixedChargeAdapter(atomType); + if ( fca.isFixedCharge() ) { + charge = fca.getCharge(); } + + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); + if ( fqa.isFluctuatingCharge() ) { + charge += atom->getFlucQPos(); + } + return charge; + } + RealType SelectionEvaluator::getCharge(Atom* atom, int frame) { + RealType charge = 0.0; + AtomType* atomType = atom->getAtomType(); + + FixedChargeAdapter fca = FixedChargeAdapter(atomType); + if ( fca.isFixedCharge() ) { + charge = fca.getCharge(); + } + + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType); + if ( fqa.isFluctuatingCharge() ) { + charge += atom->getFlucQPos(frame); + } return charge; }