ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/selection/SelectionEvaluator.cpp
Revision: 1360
Committed: Mon Sep 7 16:31:51 2009 UTC (15 years, 7 months ago) by cli2
Original Path: trunk/src/selection/SelectionEvaluator.cpp
File size: 13337 byte(s)
Log Message:
Added new restraint infrastructure
Added MolecularRestraints
Added ObjectRestraints
Added RestraintStamp
Updated thermodynamic integration to use ObjectRestraints
Added Quaternion mathematics for twist swing decompositions
Significantly updated RestWriter and RestReader to use dump-like files
Added selections for x, y, and z coordinates of atoms
Removed monolithic Restraints class
Fixed a few bugs in gradients of Euler angles in DirectionalAtom and RigidBody
Added some rotational capabilities to prinicpalAxisCalculator

File Contents

# User Rev Content
1 tim 277 /*
2     * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3     *
4     * The University of Notre Dame grants you ("Licensee") a
5     * non-exclusive, royalty free, license to use, modify and
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
19     * notice, this list of conditions and the following disclaimer.
20     *
21     * 3. Redistributions in binary form must reproduce the above copyright
22     * notice, this list of conditions and the following disclaimer in the
23     * documentation and/or other materials provided with the
24     * distribution.
25     *
26     * This software is provided "AS IS," without a warranty of any
27     * kind. All express or implied conditions, representations and
28     * warranties, including any implied warranty of merchantability,
29     * fitness for a particular purpose or non-infringement, are hereby
30     * excluded. The University of Notre Dame and its licensors shall not
31     * be liable for any damages suffered by licensee as a result of
32     * using, modifying or distributing the software or its
33     * derivatives. In no event will the University of Notre Dame or its
34     * licensors be liable for any lost revenue, profit or data, or for
35     * direct, indirect, special, consequential, incidental or punitive
36     * damages, however caused and regardless of the theory of liability,
37     * arising out of the use of or inability to use software, even if the
38     * University of Notre Dame has been advised of the possibility of
39     * such damages.
40     */
41    
42 tim 283 #include <stack>
43 tim 277 #include "selection/SelectionEvaluator.hpp"
44 tim 283 #include "primitives/Atom.hpp"
45     #include "primitives/DirectionalAtom.hpp"
46     #include "primitives/RigidBody.hpp"
47     #include "primitives/Molecule.hpp"
48    
49 tim 277 namespace oopse {
50    
51    
52 gezelter 507 SelectionEvaluator::SelectionEvaluator(SimInfo* si)
53 tim 413 : info(si), nameFinder(info), distanceFinder(info), indexFinder(info), isLoaded_(false){
54    
55 gezelter 507 nStuntDouble = info->getNGlobalAtoms() + info->getNRigidBodies();
56     }
57 tim 283
58 gezelter 507 bool SelectionEvaluator::loadScript(const std::string& filename, const std::string& script) {
59 tim 288 clearDefinitionsAndLoadPredefined();
60 tim 278 this->filename = filename;
61     this->script = script;
62     if (! compiler.compile(filename, script)) {
63 gezelter 507 error = true;
64     errorMessage = compiler.getErrorMessage();
65     std::cerr << "SelectionCompiler Error: " << errorMessage << std::endl;
66     return false;
67 tim 278 }
68    
69     pc = 0;
70     aatoken = compiler.getAatokenCompiled();
71     linenumbers = compiler.getLineNumbers();
72     lineIndices = compiler.getLineIndices();
73 tim 288
74     std::vector<std::vector<Token> >::const_iterator i;
75    
76     isDynamic_ = false;
77     for (i = aatoken.begin(); i != aatoken.end(); ++i) {
78 gezelter 507 if (containDynamicToken(*i)) {
79     isDynamic_ = true;
80     break;
81     }
82 tim 288 }
83    
84     isLoaded_ = true;
85 tim 278 return true;
86 gezelter 507 }
87 tim 278
88 gezelter 507 void SelectionEvaluator::clearState() {
89 tim 278 error = false;
90 tim 281 errorMessage = "";
91 gezelter 507 }
92 tim 278
93 gezelter 507 bool SelectionEvaluator::loadScriptString(const std::string& script) {
94 tim 278 clearState();
95 tim 281 return loadScript("", script);
96 gezelter 507 }
97 tim 278
98 gezelter 507 bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
99 tim 278 clearState();
100     return loadScriptFileInternal(filename);
101 gezelter 507 }
102 tim 278
103 gezelter 507 bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) {
104     std::ifstream ifs(filename.c_str());
105 tim 288 if (!ifs.is_open()) {
106 gezelter 507 return false;
107 tim 288 }
108    
109     const int bufferSize = 65535;
110     char buffer[bufferSize];
111     std::string script;
112     while(ifs.getline(buffer, bufferSize)) {
113 gezelter 507 script += buffer;
114 tim 288 }
115     return loadScript(filename, script);
116 gezelter 507 }
117 tim 282
118 tim 770 void SelectionEvaluator::instructionDispatchLoop(OOPSEBitSet& bs){
119 tim 288
120 tim 281 while ( pc < aatoken.size()) {
121 gezelter 507 statement = aatoken[pc++];
122     statementLength = statement.size();
123     Token token = statement[0];
124     switch (token.tok) {
125     case Token::define:
126     define();
127     break;
128     case Token::select:
129     select(bs);
130     break;
131     default:
132     unrecognizedCommand(token);
133     return;
134     }
135 tim 278 }
136 tim 288
137 gezelter 507 }
138 tim 278
139 tim 770 OOPSEBitSet SelectionEvaluator::expression(const std::vector<Token>& code, int pcStart) {
140     OOPSEBitSet bs;
141     std::stack<OOPSEBitSet> stack;
142 tim 283
143 tim 288 for (int pc = pcStart; pc < code.size(); ++pc) {
144 tim 278 Token instruction = code[pc];
145 tim 282
146 tim 278 switch (instruction.tok) {
147 tim 281 case Token::expressionBegin:
148 tim 278 break;
149 tim 281 case Token::expressionEnd:
150 tim 282 break;
151 tim 281 case Token::all:
152 tim 770 bs = OOPSEBitSet(nStuntDouble);
153 tim 283 bs.setAll();
154     stack.push(bs);
155 tim 278 break;
156 tim 281 case Token::none:
157 tim 770 bs = OOPSEBitSet(nStuntDouble);
158 tim 283 stack.push(bs);
159 tim 278 break;
160 tim 281 case Token::opOr:
161 tim 283 bs = stack.top();
162     stack.pop();
163     stack.top() |= bs;
164 tim 278 break;
165 tim 281 case Token::opAnd:
166 tim 283 bs = stack.top();
167     stack.pop();
168     stack.top() &= bs;
169 tim 278 break;
170 tim 281 case Token::opNot:
171 tim 283 stack.top().flip();
172 tim 278 break;
173 tim 281 case Token::within:
174 tim 283
175     withinInstruction(instruction, stack.top());
176 tim 278 break;
177 gezelter 507 //case Token::selected:
178     // stack.push(getSelectionSet());
179     // break;
180 tim 281 case Token::name:
181 tim 283 stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
182 tim 278 break;
183 tim 295 case Token::index:
184     stack.push(indexInstruction(instruction.value));
185 tim 281 break;
186     case Token::identifier:
187 tim 283 stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
188 tim 281 break;
189     case Token::opLT:
190     case Token::opLE:
191     case Token::opGE:
192     case Token::opGT:
193     case Token::opEQ:
194     case Token::opNE:
195 tim 283 stack.push(comparatorInstruction(instruction));
196 tim 278 break;
197     default:
198     unrecognizedExpression();
199     }
200     }
201 tim 283 if (stack.size() != 1)
202 tim 278 evalError("atom expression compiler error - stack over/underflow");
203 tim 283
204     return stack.top();
205 tim 278 }
206    
207    
208    
209 tim 770 OOPSEBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
210 tim 278 int comparator = instruction.tok;
211     int property = instruction.intValue;
212 tim 283 float comparisonValue = boost::any_cast<float>(instruction.value);
213     float propertyValue;
214 tim 770 OOPSEBitSet bs(nStuntDouble);
215 tim 295 bs.clearAll();
216 tim 283
217     SimInfo::MoleculeIterator mi;
218     Molecule* mol;
219     Molecule::AtomIterator ai;
220     Atom* atom;
221     Molecule::RigidBodyIterator rbIter;
222     RigidBody* rb;
223    
224     for (mol = info->beginMolecule(mi); mol != NULL; mol = info->nextMolecule(mi)) {
225 tim 281
226 gezelter 507 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
227     compareProperty(atom, bs, property, comparator, comparisonValue);
228     }
229 tim 283
230 gezelter 507 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
231     compareProperty(rb, bs, property, comparator, comparisonValue);
232     }
233 tim 278 }
234    
235 tim 283 return bs;
236 gezelter 507 }
237 tim 278
238 tim 770 void SelectionEvaluator::compareProperty(StuntDouble* sd, OOPSEBitSet& bs, int property, int comparator, float comparisonValue) {
239 tim 963 RealType propertyValue = 0.0;
240 gezelter 507 switch (property) {
241     case Token::mass:
242     propertyValue = sd->getMass();
243     break;
244     case Token::charge:
245     if (sd->isAtom()){
246     Atom* atom = static_cast<Atom*>(sd);
247     propertyValue = getCharge(atom);
248     } else if (sd->isRigidBody()) {
249     RigidBody* rb = static_cast<RigidBody*>(sd);
250     RigidBody::AtomIterator ai;
251     Atom* atom;
252     for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
253     propertyValue+= getCharge(atom);
254     }
255     }
256     break;
257 cli2 1360 case Token::x:
258     propertyValue = sd->getPos().x();
259     break;
260     case Token::y:
261     propertyValue = sd->getPos().y();
262     break;
263     case Token::z:
264     propertyValue = sd->getPos().z();
265     break;
266 gezelter 507 default:
267     unrecognizedAtomProperty(property);
268     }
269 tim 283
270 gezelter 507 bool match = false;
271     switch (comparator) {
272     case Token::opLT:
273     match = propertyValue < comparisonValue;
274     break;
275     case Token::opLE:
276     match = propertyValue <= comparisonValue;
277     break;
278     case Token::opGE:
279     match = propertyValue >= comparisonValue;
280     break;
281     case Token::opGT:
282     match = propertyValue > comparisonValue;
283     break;
284     case Token::opEQ:
285     match = propertyValue == comparisonValue;
286     break;
287     case Token::opNE:
288     match = propertyValue != comparisonValue;
289     break;
290     }
291     if (match)
292     bs.setBitOn(sd->getGlobalIndex());
293 tim 283
294 gezelter 507 }
295 tim 283
296 tim 770 void SelectionEvaluator::withinInstruction(const Token& instruction, OOPSEBitSet& bs){
297 tim 295
298 tim 281 boost::any withinSpec = instruction.value;
299 tim 295 float distance;
300 tim 282 if (withinSpec.type() == typeid(float)){
301 gezelter 507 distance = boost::any_cast<float>(withinSpec);
302 tim 295 } else if (withinSpec.type() == typeid(int)) {
303 gezelter 507 distance = boost::any_cast<int>(withinSpec);
304 tim 295 } else {
305 gezelter 507 evalError("casting error in withinInstruction");
306     bs.clearAll();
307 tim 278 }
308 tim 282
309 tim 295 bs = distanceFinder.find(bs, distance);
310 gezelter 507 }
311 tim 278
312 gezelter 507 void SelectionEvaluator::define() {
313 tim 282 assert(statement.size() >= 3);
314    
315     std::string variable = boost::any_cast<std::string>(statement[1].value);
316 tim 283
317 tim 351 variables.insert(VariablesType::value_type(variable, expression(statement, 2)));
318 gezelter 507 }
319 tim 278
320 tim 283
321 gezelter 507 /** @todo */
322     void SelectionEvaluator::predefine(const std::string& script) {
323 tim 278
324 tim 282 if (compiler.compile("#predefine", script)) {
325 gezelter 507 std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
326     if (aatoken.size() != 1) {
327     evalError("predefinition does not have exactly 1 command:"
328     + script);
329     return;
330     }
331     std::vector<Token> statement = aatoken[0];
332     if (statement.size() > 2) {
333     int tok = statement[1].tok;
334     if (tok == Token::identifier || (tok & Token::predefinedset) == Token::predefinedset) {
335     std::string variable = boost::any_cast<std::string>(statement[1].value);
336     variables.insert(VariablesType::value_type(variable, statement));
337 tim 282
338 gezelter 507 } else {
339     evalError("invalid variable name:" + script);
340     }
341     }else {
342     evalError("bad predefinition length:" + script);
343     }
344 tim 282
345    
346     } else {
347 gezelter 507 evalError("predefined set compile error:" + script +
348     "\ncompile error:" + compiler.getErrorMessage());
349 tim 282 }
350    
351 gezelter 507 }
352 tim 282
353 tim 770 void SelectionEvaluator::select(OOPSEBitSet& bs){
354 tim 288 bs = expression(statement, 1);
355 gezelter 507 }
356 tim 282
357 tim 770 OOPSEBitSet SelectionEvaluator::lookupValue(const std::string& variable){
358 tim 282
359 tim 770 OOPSEBitSet bs(nStuntDouble);
360 tim 282 std::map<std::string, boost::any>::iterator i = variables.find(variable);
361 tim 295
362 tim 282 if (i != variables.end()) {
363 tim 770 if (i->second.type() == typeid(OOPSEBitSet)) {
364     return boost::any_cast<OOPSEBitSet>(i->second);
365 gezelter 507 } else if (i->second.type() == typeid(std::vector<Token>)){
366     bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
367     i->second = bs; /**@todo fixme */
368     return bs;
369     }
370 tim 283 } else {
371 gezelter 507 unrecognizedIdentifier(variable);
372 tim 282 }
373 tim 295
374     return bs;
375 gezelter 507 }
376 tim 282
377 tim 770 OOPSEBitSet SelectionEvaluator::nameInstruction(const std::string& name){
378 tim 283
379 tim 295 return nameFinder.match(name);
380 tim 288
381 gezelter 507 }
382 tim 283
383 gezelter 507 bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
384 tim 288 std::vector<Token>::const_iterator i;
385     for (i = tokens.begin(); i != tokens.end(); ++i) {
386 gezelter 507 if (i->tok & Token::dynamic) {
387     return true;
388     }
389 tim 288 }
390 tim 283
391 tim 288 return false;
392 gezelter 507 }
393 tim 288
394 gezelter 507 void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
395 tim 288 variables.clear();
396     //load predefine
397     //predefine();
398 gezelter 507 }
399 tim 288
400 tim 770 OOPSEBitSet SelectionEvaluator::evaluate() {
401     OOPSEBitSet bs(nStuntDouble);
402 tim 288 if (isLoaded_) {
403 gezelter 507 pc = 0;
404     instructionDispatchLoop(bs);
405 tim 288 }
406    
407     return bs;
408 gezelter 507 }
409 tim 295
410 tim 770 OOPSEBitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
411     OOPSEBitSet bs(nStuntDouble);
412 tim 295
413     if (value.type() == typeid(int)) {
414 gezelter 507 int index = boost::any_cast<int>(value);
415     if (index < 0 || index >= bs.size()) {
416     invalidIndex(index);
417     } else {
418     bs = indexFinder.find(index);
419     }
420 tim 295 } else if (value.type() == typeid(std::pair<int, int>)) {
421 gezelter 507 std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
422     assert(indexRange.first <= indexRange.second);
423     if (indexRange.first < 0 || indexRange.second >= bs.size()) {
424     invalidIndexRange(indexRange);
425     }else {
426     bs = indexFinder.find(indexRange.first, indexRange.second);
427     }
428 tim 295 }
429    
430     return bs;
431 gezelter 507 }
432 tim 295
433 tim 432
434 tim 963 RealType SelectionEvaluator::getCharge(Atom* atom) {
435     RealType charge =0.0;
436 tim 432 AtomType* atomType = atom->getAtomType();
437     if (atomType->isCharge()) {
438 gezelter 507 GenericData* data = atomType->getPropertyByName("Charge");
439     if (data != NULL) {
440     DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data);
441 tim 432
442 gezelter 507 if (doubleData != NULL) {
443     charge = doubleData->getData();
444 tim 432
445 gezelter 507 } else {
446     sprintf( painCave.errMsg,
447     "Can not cast GenericData to DoubleGenericData\n");
448     painCave.severity = OOPSE_ERROR;
449     painCave.isFatal = 1;
450     simError();
451     }
452     }
453 tim 432 }
454    
455     return charge;
456 gezelter 507 }
457 tim 432
458     }