1 |
/* |
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. Redistributions of source code must retain the above copyright |
10 |
* notice, this list of conditions and the following disclaimer. |
11 |
* |
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. |
16 |
* |
17 |
* This software is provided "AS IS," without a warranty of any |
18 |
* kind. All express or implied conditions, representations and |
19 |
* warranties, including any implied warranty of merchantability, |
20 |
* fitness for a particular purpose or non-infringement, are hereby |
21 |
* excluded. The University of Notre Dame and its licensors shall not |
22 |
* be liable for any damages suffered by licensee as a result of |
23 |
* using, modifying or distributing the software or its |
24 |
* derivatives. In no event will the University of Notre Dame or its |
25 |
* licensors be liable for any lost revenue, profit or data, or for |
26 |
* direct, indirect, special, consequential, incidental or punitive |
27 |
* damages, however caused and regardless of the theory of liability, |
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 |
#include "selection/NameFinder.hpp" |
43 |
#include "utils/wildcards.hpp" |
44 |
#include "utils/StringTokenizer.hpp" |
45 |
#include "primitives/Molecule.hpp" |
46 |
#include "utils/StringUtils.hpp" |
47 |
namespace OpenMD { |
48 |
|
49 |
TreeNode::~TreeNode(){ |
50 |
std::map<std::string, TreeNode*>::iterator i; |
51 |
for ( i = children.begin(); i != children.end(); ++i) { |
52 |
i->second->~TreeNode(); |
53 |
} |
54 |
children.clear(); |
55 |
} |
56 |
|
57 |
NameFinder::NameFinder(SimInfo* info) : info_(info), root_(NULL){ |
58 |
nObjects_.push_back(info_->getNGlobalAtoms()+info_->getNGlobalRigidBodies()); |
59 |
nObjects_.push_back(info_->getNGlobalBonds()); |
60 |
nObjects_.push_back(info_->getNGlobalBends()); |
61 |
nObjects_.push_back(info_->getNGlobalTorsions()); |
62 |
nObjects_.push_back(info_->getNGlobalInversions()); |
63 |
loadNames(); |
64 |
} |
65 |
|
66 |
NameFinder::~NameFinder(){ |
67 |
delete root_; |
68 |
} |
69 |
|
70 |
void NameFinder::loadNames() { |
71 |
SimInfo::MoleculeIterator mi; |
72 |
Molecule::AtomIterator ai; |
73 |
Molecule::RigidBodyIterator rbIter; |
74 |
Molecule::BondIterator bondIter; |
75 |
Molecule::BendIterator bendIter; |
76 |
Molecule::TorsionIterator torsionIter; |
77 |
Molecule::InversionIterator inversionIter; |
78 |
|
79 |
Molecule* mol; |
80 |
Atom* atom; |
81 |
RigidBody* rb; |
82 |
Bond* bond; |
83 |
Bend* bend; |
84 |
Torsion* torsion; |
85 |
Inversion* inversion; |
86 |
|
87 |
root_ = new TreeNode; |
88 |
root_->bs.resize(nObjects_); |
89 |
root_->bs.setAll(); // |
90 |
|
91 |
for (mol = info_->beginMolecule(mi); mol != NULL; |
92 |
mol = info_->nextMolecule(mi)) { |
93 |
|
94 |
std::string molName = mol->getMoleculeName(); |
95 |
TreeNode* molNode = createNode(root_, molName); |
96 |
|
97 |
for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
98 |
std::string atomName = atom->getType(); |
99 |
TreeNode* atomNode = createNode(molNode, atomName); |
100 |
|
101 |
molNode->bs.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex()); |
102 |
atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex()); |
103 |
} |
104 |
|
105 |
for (rb = mol->beginRigidBody(rbIter); rb != NULL; |
106 |
rb = mol->nextRigidBody(rbIter)) { |
107 |
std::string rbName = rb->getType(); |
108 |
TreeNode* rbNode = createNode(molNode, rbName); |
109 |
|
110 |
molNode->bs.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex()); |
111 |
rbNode->bs.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex()); |
112 |
|
113 |
// COMMENTED OUT because rigid bodies are IntegrableObjects |
114 |
// (e.g. they are independently mobile, so selecting their |
115 |
// member atoms will give some odd results if we are computing |
116 |
// degrees of freedom elsewhere. |
117 |
|
118 |
// //create nodes for atoms belong to this rigidbody |
119 |
// for(atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) { |
120 |
// std::string rbAtomName = atom->getType(); |
121 |
// TreeNode* rbAtomNode = createNode(rbNode, rbAtomName); |
122 |
|
123 |
// rbAtomNode->bs.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex()); |
124 |
// } |
125 |
} |
126 |
|
127 |
for (bond = mol->beginBond(bondIter); bond != NULL; |
128 |
bond = mol->nextBond(bondIter)) { |
129 |
|
130 |
std::string bondName = bond->getName(); |
131 |
TreeNode* bondNode = createNode(molNode, bondName); |
132 |
|
133 |
molNode->bs.bitsets_[BOND].setBitOn(bond->getGlobalIndex()); |
134 |
bondNode->bs.bitsets_[BOND].setBitOn(bond->getGlobalIndex()); |
135 |
|
136 |
std::vector<Atom*> atoms = bond->getAtoms(); |
137 |
std::vector<Atom*>::iterator ai; |
138 |
|
139 |
for (ai = atoms.begin(); ai != atoms.end(); ++ai) { |
140 |
std::string atomName = (*ai)->getType(); |
141 |
TreeNode* atomNode = createNode(bondNode, atomName); |
142 |
atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex()); |
143 |
} |
144 |
} |
145 |
for (bend = mol->beginBend(bendIter); bend != NULL; |
146 |
bend = mol->nextBend(bendIter)) { |
147 |
|
148 |
std::string bendName = bend->getName(); |
149 |
TreeNode* bendNode = createNode(molNode, bendName); |
150 |
|
151 |
molNode->bs.bitsets_[BEND].setBitOn(bend->getGlobalIndex()); |
152 |
bendNode->bs.bitsets_[BEND].setBitOn(bend->getGlobalIndex()); |
153 |
|
154 |
std::vector<Atom*> atoms = bend->getAtoms(); |
155 |
std::vector<Atom*>::iterator ai; |
156 |
|
157 |
for (ai = atoms.begin(); ai != atoms.end(); ++ai) { |
158 |
std::string atomName = (*ai)->getType(); |
159 |
TreeNode* atomNode = createNode(bendNode, atomName); |
160 |
atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex()); |
161 |
} |
162 |
|
163 |
} |
164 |
for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; |
165 |
torsion = mol->nextTorsion(torsionIter)) { |
166 |
|
167 |
std::string torsionName = torsion->getName(); |
168 |
TreeNode* torsionNode = createNode(molNode, torsionName); |
169 |
|
170 |
molNode->bs.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex()); |
171 |
torsionNode->bs.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex()); |
172 |
|
173 |
std::vector<Atom*> atoms = torsion->getAtoms(); |
174 |
std::vector<Atom*>::iterator ai; |
175 |
|
176 |
for (ai = atoms.begin(); ai != atoms.end(); ++ai) { |
177 |
std::string atomName = (*ai)->getType(); |
178 |
TreeNode* atomNode = createNode(torsionNode, atomName); |
179 |
atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex()); |
180 |
} |
181 |
|
182 |
} |
183 |
for (inversion = mol->beginInversion(inversionIter); inversion != NULL; |
184 |
inversion = mol->nextInversion(inversionIter)) { |
185 |
|
186 |
std::string inversionName = inversion->getName(); |
187 |
TreeNode* inversionNode = createNode(molNode, inversionName); |
188 |
|
189 |
molNode->bs.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex()); |
190 |
inversionNode->bs.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex()); |
191 |
std::vector<Atom*> atoms = inversion->getAtoms(); |
192 |
std::vector<Atom*>::iterator ai; |
193 |
|
194 |
for (ai = atoms.begin(); ai != atoms.end(); ++ai) { |
195 |
std::string atomName = (*ai)->getType(); |
196 |
TreeNode* atomNode = createNode(inversionNode, atomName); |
197 |
atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex()); |
198 |
} |
199 |
} |
200 |
} |
201 |
} |
202 |
|
203 |
TreeNode* NameFinder::createNode(TreeNode* parent, const std::string& name) { |
204 |
TreeNode* node; |
205 |
std::map<std::string, TreeNode*>::iterator foundIter; |
206 |
foundIter = parent->children.find(name); |
207 |
if ( foundIter == parent->children.end()) { |
208 |
node = new TreeNode; |
209 |
node->name = name; |
210 |
node->bs.resize(nObjects_); |
211 |
parent->children.insert(std::make_pair(name, node)); |
212 |
}else { |
213 |
node = foundIter->second; |
214 |
} |
215 |
return node; |
216 |
} |
217 |
|
218 |
SelectionSet NameFinder::match(const std::string& name){ |
219 |
SelectionSet bs(nObjects_); |
220 |
|
221 |
StringTokenizer tokenizer(name, "."); |
222 |
|
223 |
std::vector<std::string> names; |
224 |
while(tokenizer.hasMoreTokens()) { |
225 |
names.push_back(tokenizer.nextToken()); |
226 |
} |
227 |
|
228 |
int size = names.size(); |
229 |
|
230 |
switch(size) { |
231 |
case 1 : |
232 |
//could be molecule name, atom name and rigidbody name |
233 |
matchMolecule(names[0], bs); |
234 |
matchStuntDouble("*", names[0], bs); |
235 |
matchBond("*", names[0], bs); |
236 |
matchBend("*", names[0], bs); |
237 |
matchTorsion("*", names[0], bs); |
238 |
matchInversion("*", names[0], bs); |
239 |
|
240 |
break; |
241 |
case 2: |
242 |
//could be molecule.*(include atoms and rigidbodies) or rigidbody.*(atoms belong to rigidbody) |
243 |
|
244 |
if (!isInteger(names[1])){ |
245 |
matchRigidAtoms("*", names[0], names[1], bs); |
246 |
matchStuntDouble(names[0], names[1], bs); |
247 |
} else { |
248 |
int internalIndex = lexi_cast<int>(names[1]); |
249 |
if (internalIndex < 0) { |
250 |
std::cerr << names[0] << ". " << names[1] << " is an invalid name" << std::endl; |
251 |
} else { |
252 |
matchInternalIndex(names[0], internalIndex, bs); |
253 |
} |
254 |
} |
255 |
|
256 |
break; |
257 |
case 3: |
258 |
//must be molecule.rigidbody.* |
259 |
matchRigidAtoms(names[0], names[1], names[2], bs); |
260 |
break; |
261 |
default: |
262 |
std::cerr << "invalid name: " << name << std::endl; |
263 |
break; |
264 |
} |
265 |
|
266 |
return bs; |
267 |
} |
268 |
|
269 |
void NameFinder::matchMolecule(const std::string& molName, SelectionSet& bs) { |
270 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
271 |
std::vector<TreeNode*>::iterator i; |
272 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
273 |
bs |= (*i)->bs; |
274 |
} |
275 |
} |
276 |
|
277 |
void NameFinder::matchStuntDouble(const std::string& molName, const std::string& sdName, SelectionSet& bs){ |
278 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
279 |
std::vector<TreeNode*>::iterator i; |
280 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
281 |
std::vector<TreeNode*> sdNodes = getMatchedChildren(*i, sdName); |
282 |
std::vector<TreeNode*>::iterator j; |
283 |
for (j = sdNodes.begin(); j != sdNodes.end(); ++j) { |
284 |
bs |= (*j)->bs; |
285 |
} |
286 |
} |
287 |
|
288 |
} |
289 |
|
290 |
void NameFinder::matchBond(const std::string& molName, |
291 |
const std::string& bondName, SelectionSet& bs){ |
292 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
293 |
std::vector<TreeNode*>::iterator i; |
294 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
295 |
std::vector<TreeNode*> bondNodes = getMatchedChildren(*i, bondName); |
296 |
std::vector<TreeNode*>::iterator j; |
297 |
for (j = bondNodes.begin(); j != bondNodes.end(); ++j) { |
298 |
bs |= (*j)->bs; |
299 |
std::vector<TreeNode*> bondAtomNodes = getAllChildren(*j); |
300 |
std::vector<TreeNode*>::iterator k; |
301 |
for(k = bondAtomNodes.begin(); k != bondAtomNodes.end(); ++k){ |
302 |
bs |= (*k)->bs; |
303 |
} |
304 |
} |
305 |
} |
306 |
} |
307 |
|
308 |
void NameFinder::matchBend(const std::string& molName, const std::string& bendName, SelectionSet& bs){ |
309 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
310 |
std::vector<TreeNode*>::iterator i; |
311 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
312 |
std::vector<TreeNode*> bendNodes = getMatchedChildren(*i, bendName); |
313 |
std::vector<TreeNode*>::iterator j; |
314 |
for (j = bendNodes.begin(); j != bendNodes.end(); ++j) { |
315 |
std::vector<TreeNode*> bendAtomNodes = getAllChildren(*j); |
316 |
std::vector<TreeNode*>::iterator k; |
317 |
for(k = bendAtomNodes.begin(); k != bendAtomNodes.end(); ++k){ |
318 |
bs |= (*k)->bs; |
319 |
} |
320 |
} |
321 |
} |
322 |
} |
323 |
void NameFinder::matchTorsion(const std::string& molName, const std::string& torsionName, SelectionSet& bs){ |
324 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
325 |
std::vector<TreeNode*>::iterator i; |
326 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
327 |
std::vector<TreeNode*> torsionNodes = getMatchedChildren(*i, torsionName); |
328 |
std::vector<TreeNode*>::iterator j; |
329 |
for (j = torsionNodes.begin(); j != torsionNodes.end(); ++j) { |
330 |
std::vector<TreeNode*> torsionAtomNodes = getAllChildren(*j); |
331 |
std::vector<TreeNode*>::iterator k; |
332 |
for(k = torsionAtomNodes.begin(); k != torsionAtomNodes.end(); ++k){ |
333 |
bs |= (*k)->bs; |
334 |
} |
335 |
} |
336 |
} |
337 |
} |
338 |
void NameFinder::matchInversion(const std::string& molName, const std::string& inversionName, SelectionSet& bs){ |
339 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
340 |
std::vector<TreeNode*>::iterator i; |
341 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
342 |
std::vector<TreeNode*> inversionNodes = getMatchedChildren(*i, inversionName); |
343 |
std::vector<TreeNode*>::iterator j; |
344 |
for (j = inversionNodes.begin(); j != inversionNodes.end(); ++j) { |
345 |
std::vector<TreeNode*> inversionAtomNodes = getAllChildren(*j); |
346 |
std::vector<TreeNode*>::iterator k; |
347 |
for(k = inversionAtomNodes.begin(); k != inversionAtomNodes.end(); ++k){ |
348 |
bs |= (*k)->bs; |
349 |
} |
350 |
} |
351 |
} |
352 |
} |
353 |
|
354 |
void NameFinder::matchRigidAtoms(const std::string& molName, const std::string& rbName, const std::string& rbAtomName, SelectionSet& bs){ |
355 |
std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName); |
356 |
std::vector<TreeNode*>::iterator i; |
357 |
for( i = molNodes.begin(); i != molNodes.end(); ++i ) { |
358 |
std::vector<TreeNode*> rbNodes = getMatchedChildren(*i, rbName); |
359 |
std::vector<TreeNode*>::iterator j; |
360 |
for (j = rbNodes.begin(); j != rbNodes.end(); ++j) { |
361 |
std::vector<TreeNode*> rbAtomNodes = getMatchedChildren(*j, rbAtomName); |
362 |
std::vector<TreeNode*>::iterator k; |
363 |
for(k = rbAtomNodes.begin(); k != rbAtomNodes.end(); ++k){ |
364 |
bs |= (*k)->bs; |
365 |
} |
366 |
} |
367 |
} |
368 |
|
369 |
} |
370 |
|
371 |
std::vector<TreeNode*> NameFinder::getAllChildren(TreeNode* node) { |
372 |
std::vector<TreeNode*> childNodes; |
373 |
std::map<std::string, TreeNode*>::iterator i; |
374 |
for (i = node->children.begin(); i != node->children.end(); ++i) { |
375 |
childNodes.push_back(i->second); |
376 |
} |
377 |
return childNodes; |
378 |
} |
379 |
|
380 |
std::vector<TreeNode*> NameFinder::getMatchedChildren(TreeNode* node, const std::string& name) { |
381 |
std::vector<TreeNode*> matchedNodes; |
382 |
std::map<std::string, TreeNode*>::iterator i; |
383 |
for (i = node->children.begin(); i != node->children.end(); ++i) { |
384 |
if (isMatched( i->first, name)) { |
385 |
matchedNodes.push_back(i->second); |
386 |
} |
387 |
} |
388 |
|
389 |
return matchedNodes; |
390 |
} |
391 |
|
392 |
bool NameFinder::isMatched(const std::string& str, const std::string& wildcard) { |
393 |
return Wildcard::wildcardfit(wildcard.c_str(), str.c_str()) > 0 ? true : false; |
394 |
} |
395 |
|
396 |
|
397 |
void NameFinder::matchInternalIndex(const std::string& name, int internalIndex, SelectionSet& bs){ |
398 |
|
399 |
SimInfo::MoleculeIterator mi; |
400 |
Molecule* mol; |
401 |
|
402 |
for (mol = info_->beginMolecule(mi); mol != NULL; |
403 |
mol = info_->nextMolecule(mi)) { |
404 |
|
405 |
if (isMatched(mol->getMoleculeName(), name) ) { |
406 |
int natoms = mol->getNAtoms(); |
407 |
int nrigidbodies = mol->getNRigidBodies(); |
408 |
if (internalIndex >= natoms + nrigidbodies) { |
409 |
continue; |
410 |
} else if (internalIndex < natoms) { |
411 |
bs.bitsets_[STUNTDOUBLE].setBitOn(mol->getAtomAt(internalIndex)->getGlobalIndex()); |
412 |
continue; |
413 |
} else if ( internalIndex < natoms + nrigidbodies) { |
414 |
bs.bitsets_[STUNTDOUBLE].setBitOn(mol->getRigidBodyAt(internalIndex - natoms)->getGlobalIndex()); |
415 |
} |
416 |
} |
417 |
} |
418 |
} |
419 |
|
420 |
bool NameFinder::isInteger(const std::string &str) { |
421 |
for(unsigned int i = 0; i < str.size(); ++i){ |
422 |
if (!std::isdigit(str[i])) { |
423 |
return false; |
424 |
} |
425 |
} |
426 |
return true; |
427 |
} |
428 |
} |