ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/types/MoleculeStamp.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/types/MoleculeStamp.cpp (file contents):
Revision 2469 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 2515 by tim, Fri Dec 16 18:26:41 2005 UTC

# Line 38 | Line 38
38   * University of Notre Dame has been advised of the possibility of
39   * such damages.
40   */
41 –
42 – #include <stdlib.h>
43 – #include <stdio.h>
44 – #include <string.h>
45 – #include <iostream>
41  
42 + #include <functional>
43 + #include <iostream>
44 + #include <sstream>
45   #include "types/MoleculeStamp.hpp"
46 <
46 > #include "utils/Tuple.hpp"
47 > #include "utils/MemoryUtils.hpp"
48   namespace oopse {
49 +
50 + template<class ContainerType>
51 + std::string containerToString(ContainerType& cont) {
52 +    std::ostringstream oss;
53 +    oss << "(";
54 +    typename ContainerType::iterator i = cont.begin();
55 +    if (i != cont.end()) {
56 +        oss << *i;
57 +        ++i;
58 +    }
59 +    for (; i != cont.end();++i) {
60 +        oss << ", ";
61 +        oss << *i;
62 +    }
63 +    oss << ")";
64 +    return oss.str();
65 + }
66 +
67   MoleculeStamp::MoleculeStamp() {
68      DefineParameter(Name, "name");
69      
# Line 60 | Line 77 | MoleculeStamp::~MoleculeStamp() {
77   }
78  
79   MoleculeStamp::~MoleculeStamp() {
80 <
80 >    MemoryUtils::deletePointers(atomStamps_);
81 >    MemoryUtils::deletePointers(bondStamps_);
82 >    MemoryUtils::deletePointers(bendStamps_);
83 >    MemoryUtils::deletePointers(torsionStamps_);
84 >    MemoryUtils::deletePointers(rigidBodyStamps_);
85 >    MemoryUtils::deletePointers(cutoffGroupStamps_);
86 >    MemoryUtils::deletePointers(fragmentStamps_);    
87   }
88  
89   bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
90      bool ret = addIndexSensitiveStamp(atomStamps_, atom);
91      if (!ret) {
92 <        std::cout << "multiple atoms have the same index: " << atom->getIndex() <<" in " << getName()  << " Molecule\n";
92 >         std::cout<< "Error in Molecule " << getName()  << ": multiple atoms have the same indices"<< atom->getIndex() <<"\n";
93      }
94      return ret;
95      
# Line 90 | Line 113 | bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp*
113   bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
114      bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
115      if (!ret) {
116 <        std::cout << "multiple rigidbodies have the same index: " << rigidbody->getIndex() <<" in " << getName()  << " Molecule\n";
116 >        std::cout<< "Error in Molecule " << getName()  << ": multiple rigidbodies have the same indices: " << rigidbody->getIndex() <<"\n";
117      }
118      return ret;
119   }
# Line 107 | Line 130 | void MoleculeStamp::validate() {
130   void MoleculeStamp::validate() {
131      DataHolder::validate();
132  
133 +    atom2Rigidbody.resize(getNAtoms());
134 +    // negative number means atom is a free atom, does not belong to rigidbody
135 +    //every element in atom2Rigidbody has unique negative number at the very beginning
136 +    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
137 +        atom2Rigidbody[i] = -1 - i;
138 +    }
139 +    for (int i = 0; i < getNRigidBodies(); ++i) {
140 +        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
141 +        std::vector<int> members = rbStamp->getMembers();
142 +        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
143 +            atom2Rigidbody[*j] = i;                
144 +        }
145 +    }
146 +
147 +    checkAtoms();
148 +    checkBonds();
149 +    fillBondInfo();
150 +    checkBends();
151 +    checkTorsions();
152 +    checkRigidBodies();
153 +    checkCutoffGroups();
154 +    checkFragments();
155 +
156 +    int nrigidAtoms = 0;
157 +    for (int i = 0; i < getNRigidBodies(); ++i) {
158 +        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
159 +        nrigidAtoms += rbStamp->getNMembers();
160 +    }
161 +    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
162 +
163 + }
164 +
165 + void MoleculeStamp::checkAtoms() {
166      std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
167      if (ai != atomStamps_.end()) {
168          std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
169      }
170  
171 <     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
116 <     if (ri != rigidBodyStamps_.end()) {
117 <         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
118 <     }
119 <    
120 <    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
121 <    if (fi != fragmentStamps_.end()) {
122 <        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
123 <    }
171 > }
172  
173 + void MoleculeStamp::checkBonds() {
174      //make sure index is not out of range
175      int natoms = getNAtoms();
176      for(int i = 0; i < getNBonds(); ++i) {
177          BondStamp* bondStamp = getBondStamp(i);
178          if (bondStamp->getA() >=  natoms && bondStamp->getB() >= natoms) {
179 <            std::cout << "Error in Molecule " << getName() <<  ": bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " is invalid\n";
179 >            std::cout << "Error in Molecule " << getName() <<  ": bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") is invalid\n";
180          }
181      }
182 +    
183 +    //make sure bonds are unique
184 +    std::set<std::pair<int, int> > allBonds;
185 +    for(int i = 0; i < getNBonds(); ++i) {
186 +        BondStamp* bondStamp= getBondStamp(i);        
187 +        std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
188 +        //make sure bondTuple.first is always less than or equal to bondTuple.third
189 +        if (bondPair.first > bondPair.second) {
190 +            std::swap(bondPair.first, bondPair.second);
191 +        }
192 +        
193 +        std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
194 +        if ( iter != allBonds.end()) {
195 +            std::cout << "Error in Molecule " << getName() << ": " << "bond(" <<iter->first << ", "<< iter->second << ")appears multiple times\n";
196 +        } else {
197 +            allBonds.insert(bondPair);
198 +        }
199 +    }
200 +    
201 +    //make sure atoms belong to same rigidbody do not bond to each other
202 +    for(int i = 0; i < getNBonds(); ++i) {
203 +        BondStamp* bondStamp = getBondStamp(i);
204 +        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
205 +            std::cout << "Error in Molecule " << getName() << ": "<<"bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
206 +        }
207 +    }
208 +    
209 + }
210 +
211 + void MoleculeStamp::checkBends() {
212      for(int i = 0; i < getNBends(); ++i) {
213          BendStamp* bendStamp = getBendStamp(i);
214          std::vector<int> bendAtoms =  bendStamp->getMembers();
215 <        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
215 >        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
216          if (j != bendAtoms.end()) {
217 <            std::cout << "Error in Molecule " << getName();
217 >            std::cout << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << "have invalid indices\n";
218          }
219  
220 <        if (bendAtoms.size() == 2 && !bendStamp->haveGhostVectorSource()) {
221 <            std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing";
220 >        if (bendAtoms.size() == 2 ) {
221 >            if (!bendStamp->haveGhostVectorSource()) {
222 >                std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
223 >            }else{
224 >                int ghostIndex = bendStamp->getGhostVectorSource();
225 >                if (ghostIndex < getNAtoms()) {
226 >                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
227 >                      std::cout <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
228 >                    }
229 >                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
230 >                        std::cout <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
231 >                    }
232 >                }else {
233 >                    std::cout << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
234 >                }
235 >            }
236 >        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
237 >            std::cout <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
238          }
144 –    }    
145 –    for(int i = 0; i < getNBends(); ++i) {
146 –        TorsionStamp* torsionStamp = getTorsionStamp(i);
147 –        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
148 –        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
149 –        if (j != torsionAtoms.end()) {
150 –            std::cout << "Error in Molecule " << getName();
151 –        }
239      }
153 –    for(int i = 0; i < getNCutoffGroups(); ++i) {
154 –        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
155 –        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
156 –        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
157 –        if (j != cutoffGroupAtoms.end()) {
158 –            std::cout << "Error in Molecule " << getName();
159 –        }
160 –    }
161 –        
162 –    atom2Rigidbody.resize(natoms);
163 –    // negative number means atom is a free atom, does not belong to rigidbody
164 –    //every element in atom2Rigidbody has unique negative number at the very beginning
165 –    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
166 –        atom2Rigidbody[i] = -1 - i;
167 –    }
240  
169 –    for (int i = 0; i < getNRigidBodies(); ++i) {
170 –        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
171 –        std::vector<int> members = rbStamp->getMembers();
172 –        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
173 –            atom2Rigidbody[*j] = i;                
174 –        }
175 –    }
176 –    //make sure atoms belong to same rigidbody do not bond to each other
177 –    for(int i = 0; i < getNBonds(); ++i) {
178 –        BondStamp* bondStamp = getBondStamp(i);
179 –        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()])
180 –            std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << "belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
181 –        }
182 –        
241      for(int i = 0; i < getNBends(); ++i) {
242          BendStamp* bendStamp = getBendStamp(i);
243          std::vector<int> bendAtoms =  bendStamp->getMembers();
# Line 190 | Line 248 | void MoleculeStamp::validate() {
248              if (rigidbodyIndex >= 0) {
249                  ++rigidSet[rigidbodyIndex];
250                  if (rigidSet[rigidbodyIndex] > 1) {
251 <                    std::cout << "Error in Molecule " << getName() << ": ";
194 <                    //std::cout << "atoms of bend " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
251 >                    std::cout << "Error in Molecule " << getName() << ": bend" << containerToString(bendAtoms) << " belong to same rigidbody " << rigidbodyIndex << "\n";                    
252                  }
253              }
254          }
255 <    }      
255 >    }
256 >    
257 >    
258 >    std::set<IntTuple3> allBends;
259 >    std::set<IntTuple3>::iterator iter;
260 >    for(int i = 0; i < getNBends(); ++i) {
261 >        BendStamp* bendStamp= getBendStamp(i);
262 >        std::vector<int> bend = bendStamp->getMembers();
263 >        if (bend.size() == 2) {
264 >        // in case we have two ghost bend. For example,
265 >        // bend {
266 >        // members (0, 1);
267 >        //   ghostVectorSource = 0;
268 >        // }
269 >        // and
270 >        // bend {
271 >        //   members (0, 1);
272 >        // ghostVectorSource = 0;
273 >        // }
274 >        // In order to distinguish them. we expand them to Tuple3.
275 >        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
276 >             int ghostIndex = bendStamp->getGhostVectorSource();
277 >             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
278 >             if (j != bend.end()) {
279 >                bend.insert(j, ghostIndex);
280 >             }
281 >        }
282 >        
283 >        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
284 >        //make sure bendTuple.first is always less than or equal to bendTuple.third
285 >        if (bendTuple.first > bendTuple.third) {
286 >            std::swap(bendTuple.first, bendTuple.third);
287 >        }
288 >        
289 >        iter = allBends.find(bendTuple);
290 >        if ( iter != allBends.end()) {
291 >            std::cout << "Error in Molecule " << getName() << ": " << "Bend appears multiple times\n";
292 >        } else {
293 >            allBends.insert(bendTuple);
294 >        }
295 >    }
296 >
297 >    for (int i = 0; i < getNBonds(); ++i) {
298 >        BondStamp* bondStamp = getBondStamp(i);
299 >        int a = bondStamp->getA();
300 >        int b = bondStamp->getB();
301 >
302 >        AtomStamp* atomA = getAtomStamp(a);
303 >        AtomStamp* atomB = getAtomStamp(b);
304 >
305 >        //find bend c--a--b
306 >        AtomStamp::AtomIter ai;
307 >        for(int c= atomA->getFirstBonedAtom(ai);c != -1;c = atomA->getNextBonedAtom(ai))
308 >        {
309 >            if(b == c)
310 >                continue;          
311 >            
312 >            IntTuple3 newBend(c, a, b);
313 >            if (newBend.first > newBend.third) {
314 >                std::swap(newBend.first, newBend.third);
315 >            }
316 >
317 >            if (allBends.find(newBend) == allBends.end() ) {                
318 >                allBends.insert(newBend);
319 >                BendStamp * newBendStamp = new BendStamp();
320 >                newBendStamp->setMembers(newBend);
321 >                addBendStamp(newBendStamp);
322 >            }
323 >        }        
324 >
325 >        //find bend a--b--c
326 >        for(int c= atomB->getFirstBonedAtom(ai);c != -1;c = atomB->getNextBonedAtom(ai))
327 >        {
328 >            if(a == c)
329 >                continue;          
330 >
331 >            IntTuple3 newBend( a, b, c);
332 >            if (newBend.first > newBend.third) {
333 >                std::swap(newBend.first, newBend.third);
334 >            }            
335 >            if (allBends.find(newBend) == allBends.end() ) {                
336 >                allBends.insert(newBend);
337 >                BendStamp * newBendStamp = new BendStamp();
338 >                newBendStamp->setMembers(newBend);
339 >                addBendStamp(newBendStamp);
340 >            }
341 >        }        
342 >    }
343 >
344 > }
345 >
346 > void MoleculeStamp::checkTorsions() {
347      for(int i = 0; i < getNTorsions(); ++i) {
348          TorsionStamp* torsionStamp = getTorsionStamp(i);
349 +        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
350 +        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
351 +        if (j != torsionAtoms.end()) {
352 +            std::cout << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n";
353 +        }
354 +    }
355 +    
356 +    for(int i = 0; i < getNTorsions(); ++i) {
357 +        TorsionStamp* torsionStamp = getTorsionStamp(i);
358          std::vector<int> torsionAtoms =  torsionStamp->getMembers();
359          std::vector<int> rigidSet(getNRigidBodies(), 0);
360          std::vector<int>::iterator j;
# Line 206 | Line 363 | void MoleculeStamp::validate() {
363              if (rigidbodyIndex >= 0) {
364                  ++rigidSet[rigidbodyIndex];
365                  if (rigidSet[rigidbodyIndex] > 1) {
366 <                    std::cout << "Error in Molecule " << getName() << ": ";
210 <                    //std::cout << "atoms of torsion " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
366 >                    std::cout << "Error in Molecule " << getName() << ": torsion" << containerToString(torsionAtoms) << "is invalid\n";                  
367                  }
368              }
369          }
370 <    }
370 >    }    
371  
372 +    std::set<IntTuple4> allTorsions;
373 +    std::set<IntTuple4>::iterator iter;
374 +     for(int i = 0; i < getNTorsions(); ++i) {
375 +         TorsionStamp* torsionStamp= getTorsionStamp(i);
376 +         std::vector<int> torsion = torsionStamp->getMembers();
377 +         if (torsion.size() == 3) {
378 +             int ghostIndex = torsionStamp->getGhostVectorSource();
379 +             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
380 +             if (j != torsion.end()) {
381 +                torsion.insert(j, ghostIndex);
382 +             }
383 +         }
384  
385 <    //fill in bond information into atom
386 <    fillBondInfo();
387 <    findBends();
388 <    findTorsions();
385 >        IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
386 >        if (torsionTuple.first > torsionTuple.fourth) {
387 >            std::swap(torsionTuple.first, torsionTuple.fourth);
388 >            std::swap(torsionTuple.second, torsionTuple.third);                    
389 >        }                
390  
391 <    int nrigidAtoms = 0;
391 >         iter = allTorsions.find(torsionTuple);
392 >         if ( iter == allTorsions.end()) {
393 >            allTorsions.insert(torsionTuple);
394 >         } else {
395 >            std::cout << "Error in Molecule " << getName() << ": " << "Torsion appears multiple times\n";
396 >         }
397 >     }
398 >
399 >    for (int i = 0; i < getNBonds(); ++i) {
400 >        BondStamp* bondStamp = getBondStamp(i);
401 >        int b = bondStamp->getA();
402 >        int c = bondStamp->getB();
403 >
404 >        AtomStamp* atomB = getAtomStamp(b);
405 >        AtomStamp* atomC = getAtomStamp(c);
406 >
407 >        AtomStamp::AtomIter ai2;
408 >        AtomStamp::AtomIter ai3;
409 >
410 >        for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
411 >        {
412 >            if(a == c)
413 >                continue;
414 >
415 >            for(int d = atomC->getFirstBonedAtom(ai3);d != -1;d = atomC->getNextBonedAtom(ai3))
416 >            {
417 >                if(d == b)
418 >                    continue;
419 >                
420 >                IntTuple4 newTorsion(a, b, c, d);
421 >                //make sure the first element is always less than or equal to the fourth element in IntTuple4
422 >                if (newTorsion.first > newTorsion.fourth) {
423 >                    std::swap(newTorsion.first, newTorsion.fourth);
424 >                    std::swap(newTorsion.second, newTorsion.third);                    
425 >                }                
426 >                if (allTorsions.find(newTorsion) == allTorsions.end() ) {                
427 >                    allTorsions.insert(newTorsion);
428 >                    TorsionStamp * newTorsionStamp = new TorsionStamp();
429 >                    newTorsionStamp->setMembers(newTorsion);
430 >                    addTorsionStamp(newTorsionStamp);                    
431 >                }            
432 >            }
433 >        }    
434 >    }
435 >
436 > }
437 >
438 > void MoleculeStamp::checkRigidBodies() {
439 >     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
440 >     if (ri != rigidBodyStamps_.end()) {
441 >         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
442 >     }
443 >
444      for (int i = 0; i < getNRigidBodies(); ++i) {
445          RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
446 <        nrigidAtoms += rbStamp->getNMembers();
447 <    }
448 <    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
446 >        std::vector<int> rigidAtoms =  rbStamp ->getMembers();
447 >        std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
448 >        if (j != rigidAtoms.end()) {
449 >            std::cout << "Error in Molecule " << getName();
450 >        }
451 >        
452 >    }    
453 > }
454  
455 < }
455 > void MoleculeStamp::checkCutoffGroups() {
456  
457 < void MoleculeStamp::fillBondInfo() {
457 >    for(int i = 0; i < getNCutoffGroups(); ++i) {
458 >        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
459 >        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
460 >        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
461 >        if (j != cutoffGroupAtoms.end()) {
462 >            std::cout << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
463 >        }
464 >    }    
465   }
466  
467 < void MoleculeStamp::findBends() {
467 > void MoleculeStamp::checkFragments() {
468  
469 +    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
470 +    if (fi != fragmentStamps_.end()) {
471 +        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
472 +    }
473 +    
474   }
475  
476 < void MoleculeStamp::findTorsions() {
476 > void MoleculeStamp::fillBondInfo() {
477  
478 +    for (int i = 0; i < getNBonds(); ++i) {
479 +        BondStamp* bondStamp = getBondStamp(i);
480 +        int a = bondStamp->getA();
481 +        int b = bondStamp->getB();
482 +        AtomStamp* atomA = getAtomStamp(a);
483 +        AtomStamp* atomB = getAtomStamp(b);
484 +        atomA->addBond(i);
485 +        atomA->addBondedAtom(b);
486 +        atomB->addBond(i);        
487 +        atomB->addBondedAtom(a);
488 +
489 +    }
490   }
491  
492 +
493 +
494   //Function Name: isBondInSameRigidBody
495   //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
496   bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
# Line 262 | Line 514 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
514   // Function Name: isAtomInRigidBody
515   //return false if atom does not belong to a rigid body, otherwise return true
516   bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
517 <  int whichRigidBody;
266 <  int consAtomIndex;
267 <
268 <  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
517 >  return atom2Rigidbody[atomIndex] >=0 ;
518    
519   }
520  
# Line 276 | Line 525 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex, i
525   //whichRigidBody: the index of rigidbody in component
526   //consAtomIndex:  the position of joint atom apears in  rigidbody's definition
527   bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
279 –  RigidBodyStamp* rbStamp;
280 –  int numRb;
281 –  int numAtom;
528  
529 +  
530 +
531    whichRigidBody = -1;
532    consAtomIndex = -1;
533  
534 <  numRb = this->getNRigidBodies();
535 <  
536 <  for(int i = 0 ; i < numRb; i++){
537 <    rbStamp = this->getRigidBodyStamp(i);
538 <    numAtom = rbStamp->getNMembers();
291 <    for(int j = 0; j < numAtom; j++)
534 >  if (atom2Rigidbody[atomIndex] >=0) {
535 >    whichRigidBody = atom2Rigidbody[atomIndex];
536 >    RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
537 >    int numAtom = rbStamp->getNMembers();
538 >    for(int j = 0; j < numAtom; j++) {
539        if (rbStamp->getMemberAt(j) == atomIndex){
293 –        whichRigidBody = i;
540          consAtomIndex = j;
541          return true;
542        }
543 +    }
544    }
545  
546    return false;

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines