ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/selection/NameFinder.cpp
Revision: 1953
Committed: Thu Dec 5 18:19:26 2013 UTC (11 years, 5 months ago) by gezelter
File size: 15798 byte(s)
Log Message:
Rewrote much of selection module, added a bond correlation function

File Contents

# User Rev Content
1 tim 283 /*
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 gezelter 1390 * 1. Redistributions of source code must retain the above copyright
10 tim 283 * notice, this list of conditions and the following disclaimer.
11     *
12 gezelter 1390 * 2. Redistributions in binary form must reproduce the above copyright
13 tim 283 * 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 gezelter 1390 *
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 gezelter 1879 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 gezelter 1782 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40     * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 tim 283 */
42     #include "selection/NameFinder.hpp"
43 tim 284 #include "utils/wildcards.hpp"
44 tim 287 #include "utils/StringTokenizer.hpp"
45     #include "primitives/Molecule.hpp"
46 tim 303 #include "utils/StringUtils.hpp"
47 gezelter 1390 namespace OpenMD {
48 tim 283
49 gezelter 507 TreeNode::~TreeNode(){
50 tim 284 std::map<std::string, TreeNode*>::iterator i;
51     for ( i = children.begin(); i != children.end(); ++i) {
52 gezelter 507 i->second->~TreeNode();
53 tim 284 }
54     children.clear();
55 gezelter 507 }
56 tim 284
57 gezelter 507 NameFinder::NameFinder(SimInfo* info) : info_(info), root_(NULL){
58 gezelter 1953 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 tim 283 loadNames();
64 gezelter 507 }
65 tim 283
66 gezelter 507 NameFinder::~NameFinder(){
67 tim 284 delete root_;
68 gezelter 507 }
69 tim 283
70 gezelter 507 void NameFinder::loadNames() {
71 gezelter 1953 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 tim 283
79     Molecule* mol;
80     Atom* atom;
81     RigidBody* rb;
82 gezelter 1953 Bond* bond;
83     Bend* bend;
84     Torsion* torsion;
85     Inversion* inversion;
86 tim 284
87     root_ = new TreeNode;
88 gezelter 1953 root_->bs.resize(nObjects_);
89 tim 284 root_->bs.setAll(); //
90    
91 gezelter 1801 for (mol = info_->beginMolecule(mi); mol != NULL;
92     mol = info_->nextMolecule(mi)) {
93 tim 288
94 gezelter 507 std::string molName = mol->getMoleculeName();
95 gezelter 1953 TreeNode* molNode = createNode(root_, molName);
96 tim 283
97 gezelter 507 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
98     std::string atomName = atom->getType();
99 gezelter 1953 TreeNode* atomNode = createNode(molNode, atomName);
100 tim 288
101 gezelter 1953 molNode->bs.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
102     atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
103 gezelter 507 }
104 tim 284
105 gezelter 1801 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
106     rb = mol->nextRigidBody(rbIter)) {
107 gezelter 507 std::string rbName = rb->getType();
108 gezelter 1953 TreeNode* rbNode = createNode(molNode, rbName);
109 tim 284
110 gezelter 1953 molNode->bs.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
111     rbNode->bs.bitsets_[STUNTDOUBLE].setBitOn(rb->getGlobalIndex());
112 tim 284
113 gezelter 507 //create nodes for atoms belong to this rigidbody
114     for(atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
115     std::string rbAtomName = atom->getType();
116 gezelter 1953 TreeNode* rbAtomNode = createNode(rbNode, rbAtomName);
117 tim 288
118 gezelter 1953 rbAtomNode->bs.bitsets_[STUNTDOUBLE].setBitOn(atom->getGlobalIndex());
119 gezelter 507 }
120     }
121 gezelter 1953
122     for (bond = mol->beginBond(bondIter); bond != NULL;
123     bond = mol->nextBond(bondIter)) {
124    
125     std::string bondName = bond->getName();
126     TreeNode* bondNode = createNode(molNode, bondName);
127    
128     molNode->bs.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
129     bondNode->bs.bitsets_[BOND].setBitOn(bond->getGlobalIndex());
130    
131     std::vector<Atom*> atoms = bond->getAtoms();
132     std::vector<Atom*>::iterator ai;
133    
134     for (ai = atoms.begin(); ai != atoms.end(); ++ai) {
135     std::string atomName = (*ai)->getType();
136     TreeNode* atomNode = createNode(bondNode, atomName);
137     atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex());
138     }
139     }
140     for (bend = mol->beginBend(bendIter); bend != NULL;
141     bend = mol->nextBend(bendIter)) {
142    
143     std::string bendName = bend->getName();
144     TreeNode* bendNode = createNode(molNode, bendName);
145    
146     molNode->bs.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
147     bendNode->bs.bitsets_[BEND].setBitOn(bend->getGlobalIndex());
148    
149     std::vector<Atom*> atoms = bend->getAtoms();
150     std::vector<Atom*>::iterator ai;
151    
152     for (ai = atoms.begin(); ai != atoms.end(); ++ai) {
153     std::string atomName = (*ai)->getType();
154     TreeNode* atomNode = createNode(bendNode, atomName);
155     atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex());
156     }
157    
158     }
159     for (torsion = mol->beginTorsion(torsionIter); torsion != NULL;
160     torsion = mol->nextTorsion(torsionIter)) {
161    
162     std::string torsionName = torsion->getName();
163     TreeNode* torsionNode = createNode(molNode, torsionName);
164    
165     molNode->bs.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
166     torsionNode->bs.bitsets_[TORSION].setBitOn(torsion->getGlobalIndex());
167    
168     std::vector<Atom*> atoms = torsion->getAtoms();
169     std::vector<Atom*>::iterator ai;
170    
171     for (ai = atoms.begin(); ai != atoms.end(); ++ai) {
172     std::string atomName = (*ai)->getType();
173     TreeNode* atomNode = createNode(torsionNode, atomName);
174     atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex());
175     }
176    
177     }
178     for (inversion = mol->beginInversion(inversionIter); inversion != NULL;
179     inversion = mol->nextInversion(inversionIter)) {
180    
181     std::string inversionName = inversion->getName();
182     TreeNode* inversionNode = createNode(molNode, inversionName);
183    
184     molNode->bs.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
185     inversionNode->bs.bitsets_[INVERSION].setBitOn(inversion->getGlobalIndex());
186     std::vector<Atom*> atoms = inversion->getAtoms();
187     std::vector<Atom*>::iterator ai;
188    
189     for (ai = atoms.begin(); ai != atoms.end(); ++ai) {
190     std::string atomName = (*ai)->getType();
191     TreeNode* atomNode = createNode(inversionNode, atomName);
192     atomNode->bs.bitsets_[STUNTDOUBLE].setBitOn((*ai)->getGlobalIndex());
193     }
194     }
195 gezelter 1801 }
196 gezelter 507 }
197 tim 283
198 gezelter 507 TreeNode* NameFinder::createNode(TreeNode* parent, const std::string& name) {
199 tim 288 TreeNode* node;
200     std::map<std::string, TreeNode*>::iterator foundIter;
201     foundIter = parent->children.find(name);
202     if ( foundIter == parent->children.end()) {
203 gezelter 507 node = new TreeNode;
204     node->name = name;
205 gezelter 1953 node->bs.resize(nObjects_);
206 gezelter 507 parent->children.insert(std::make_pair(name, node));
207 tim 288 }else {
208 gezelter 507 node = foundIter->second;
209 tim 288 }
210     return node;
211 gezelter 507 }
212 tim 283
213 gezelter 1953 SelectionSet NameFinder::match(const std::string& name){
214     SelectionSet bs(nObjects_);
215 tim 288
216 tim 283 StringTokenizer tokenizer(name, ".");
217    
218     std::vector<std::string> names;
219     while(tokenizer.hasMoreTokens()) {
220 gezelter 507 names.push_back(tokenizer.nextToken());
221 tim 283 }
222    
223     int size = names.size();
224 gezelter 1953
225 tim 283 switch(size) {
226 gezelter 507 case 1 :
227     //could be molecule name, atom name and rigidbody name
228     matchMolecule(names[0], bs);
229     matchStuntDouble("*", names[0], bs);
230 gezelter 1953 matchBond("*", names[0], bs);
231     matchBend("*", names[0], bs);
232     matchTorsion("*", names[0], bs);
233     matchInversion("*", names[0], bs);
234 tim 284
235 gezelter 507 break;
236     case 2:
237     //could be molecule.*(include atoms and rigidbodies) or rigidbody.*(atoms belong to rigidbody)
238 tim 303
239 gezelter 507 if (!isInteger(names[1])){
240     matchRigidAtoms("*", names[0], names[1], bs);
241     matchStuntDouble(names[0], names[1], bs);
242     } else {
243     int internalIndex = lexi_cast<int>(names[1]);
244     if (internalIndex < 0) {
245     std::cerr << names[0] << ". " << names[1] << " is an invalid name" << std::endl;
246     } else {
247     matchInternalIndex(names[0], internalIndex, bs);
248     }
249     }
250 tim 284
251 gezelter 507 break;
252     case 3:
253     //must be molecule.rigidbody.*
254     matchRigidAtoms(names[0], names[1], names[2], bs);
255     break;
256     default:
257     std::cerr << "invalid name: " << name << std::endl;
258     break;
259 tim 283 }
260    
261 tim 288 return bs;
262 gezelter 507 }
263 tim 283
264 gezelter 1953 void NameFinder::matchMolecule(const std::string& molName, SelectionSet& bs) {
265 tim 284 std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
266     std::vector<TreeNode*>::iterator i;
267     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
268 gezelter 507 bs |= (*i)->bs;
269 tim 284 }
270 gezelter 507 }
271 tim 284
272 gezelter 1953 void NameFinder::matchStuntDouble(const std::string& molName, const std::string& sdName, SelectionSet& bs){
273 tim 284 std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
274     std::vector<TreeNode*>::iterator i;
275     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
276 gezelter 507 std::vector<TreeNode*> sdNodes = getMatchedChildren(*i, sdName);
277     std::vector<TreeNode*>::iterator j;
278     for (j = sdNodes.begin(); j != sdNodes.end(); ++j) {
279     bs |= (*j)->bs;
280     }
281 tim 284 }
282    
283 gezelter 507 }
284 tim 284
285 gezelter 1953 void NameFinder::matchBond(const std::string& molName,
286     const std::string& bondName, SelectionSet& bs){
287 tim 284 std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
288     std::vector<TreeNode*>::iterator i;
289     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
290 gezelter 1953 std::vector<TreeNode*> bondNodes = getMatchedChildren(*i, bondName);
291     std::vector<TreeNode*>::iterator j;
292     for (j = bondNodes.begin(); j != bondNodes.end(); ++j) {
293     bs |= (*j)->bs;
294     std::vector<TreeNode*> bondAtomNodes = getAllChildren(*j);
295     std::vector<TreeNode*>::iterator k;
296     for(k = bondAtomNodes.begin(); k != bondAtomNodes.end(); ++k){
297     bs |= (*k)->bs;
298     }
299     }
300     }
301     }
302    
303     void NameFinder::matchBend(const std::string& molName, const std::string& bendName, SelectionSet& bs){
304     std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
305     std::vector<TreeNode*>::iterator i;
306     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
307     std::vector<TreeNode*> bendNodes = getMatchedChildren(*i, bendName);
308     std::vector<TreeNode*>::iterator j;
309     for (j = bendNodes.begin(); j != bendNodes.end(); ++j) {
310     std::vector<TreeNode*> bendAtomNodes = getAllChildren(*j);
311     std::vector<TreeNode*>::iterator k;
312     for(k = bendAtomNodes.begin(); k != bendAtomNodes.end(); ++k){
313     bs |= (*k)->bs;
314     }
315     }
316     }
317     }
318     void NameFinder::matchTorsion(const std::string& molName, const std::string& torsionName, SelectionSet& bs){
319     std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
320     std::vector<TreeNode*>::iterator i;
321     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
322     std::vector<TreeNode*> torsionNodes = getMatchedChildren(*i, torsionName);
323     std::vector<TreeNode*>::iterator j;
324     for (j = torsionNodes.begin(); j != torsionNodes.end(); ++j) {
325     std::vector<TreeNode*> torsionAtomNodes = getAllChildren(*j);
326     std::vector<TreeNode*>::iterator k;
327     for(k = torsionAtomNodes.begin(); k != torsionAtomNodes.end(); ++k){
328     bs |= (*k)->bs;
329     }
330     }
331     }
332     }
333     void NameFinder::matchInversion(const std::string& molName, const std::string& inversionName, SelectionSet& bs){
334     std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
335     std::vector<TreeNode*>::iterator i;
336     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
337     std::vector<TreeNode*> inversionNodes = getMatchedChildren(*i, inversionName);
338     std::vector<TreeNode*>::iterator j;
339     for (j = inversionNodes.begin(); j != inversionNodes.end(); ++j) {
340     std::vector<TreeNode*> inversionAtomNodes = getAllChildren(*j);
341     std::vector<TreeNode*>::iterator k;
342     for(k = inversionAtomNodes.begin(); k != inversionAtomNodes.end(); ++k){
343     bs |= (*k)->bs;
344     }
345     }
346     }
347     }
348    
349     void NameFinder::matchRigidAtoms(const std::string& molName, const std::string& rbName, const std::string& rbAtomName, SelectionSet& bs){
350     std::vector<TreeNode*> molNodes = getMatchedChildren(root_, molName);
351     std::vector<TreeNode*>::iterator i;
352     for( i = molNodes.begin(); i != molNodes.end(); ++i ) {
353 gezelter 507 std::vector<TreeNode*> rbNodes = getMatchedChildren(*i, rbName);
354     std::vector<TreeNode*>::iterator j;
355     for (j = rbNodes.begin(); j != rbNodes.end(); ++j) {
356     std::vector<TreeNode*> rbAtomNodes = getMatchedChildren(*j, rbAtomName);
357     std::vector<TreeNode*>::iterator k;
358     for(k = rbAtomNodes.begin(); k != rbAtomNodes.end(); ++k){
359     bs |= (*k)->bs;
360     }
361     }
362 tim 284 }
363    
364 gezelter 507 }
365 tim 284
366 gezelter 1953 std::vector<TreeNode*> NameFinder::getAllChildren(TreeNode* node) {
367     std::vector<TreeNode*> childNodes;
368     std::map<std::string, TreeNode*>::iterator i;
369     for (i = node->children.begin(); i != node->children.end(); ++i) {
370     childNodes.push_back(i->second);
371     }
372     return childNodes;
373     }
374 tim 284
375 gezelter 507 std::vector<TreeNode*> NameFinder::getMatchedChildren(TreeNode* node, const std::string& name) {
376 tim 284 std::vector<TreeNode*> matchedNodes;
377     std::map<std::string, TreeNode*>::iterator i;
378     for (i = node->children.begin(); i != node->children.end(); ++i) {
379 gezelter 507 if (isMatched( i->first, name)) {
380     matchedNodes.push_back(i->second);
381     }
382 tim 284 }
383    
384     return matchedNodes;
385 gezelter 507 }
386 tim 284
387 gezelter 507 bool NameFinder::isMatched(const std::string& str, const std::string& wildcard) {
388 gezelter 1782 return Wildcard::wildcardfit(wildcard.c_str(), str.c_str()) > 0 ? true : false;
389 gezelter 507 }
390 tim 284
391 tim 303
392 gezelter 1953 void NameFinder::matchInternalIndex(const std::string& name, int internalIndex, SelectionSet& bs){
393 tim 303
394     SimInfo::MoleculeIterator mi;
395     Molecule* mol;
396    
397 gezelter 1801 for (mol = info_->beginMolecule(mi); mol != NULL;
398     mol = info_->nextMolecule(mi)) {
399 tim 303
400 gezelter 507 if (isMatched(mol->getMoleculeName(), name) ) {
401     int natoms = mol->getNAtoms();
402     int nrigidbodies = mol->getNRigidBodies();
403     if (internalIndex >= natoms + nrigidbodies) {
404     continue;
405     } else if (internalIndex < natoms) {
406 gezelter 1953 bs.bitsets_[STUNTDOUBLE].setBitOn(mol->getAtomAt(internalIndex)->getGlobalIndex());
407 gezelter 507 continue;
408     } else if ( internalIndex < natoms + nrigidbodies) {
409 gezelter 1953 bs.bitsets_[STUNTDOUBLE].setBitOn(mol->getRigidBodyAt(internalIndex - natoms)->getGlobalIndex());
410 gezelter 507 }
411     }
412 gezelter 1801 }
413 gezelter 507 }
414 tim 303
415 gezelter 1879 bool NameFinder::isInteger(const std::string &str) {
416 gezelter 1782 for(unsigned int i = 0; i < str.size(); ++i){
417 gezelter 507 if (!std::isdigit(str[i])) {
418     return false;
419     }
420 tim 303 }
421     return true;
422 gezelter 507 }
423 tim 303 }

Properties

Name Value
svn:executable *
svn:keywords Author Id Revision Date