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

Comparing trunk/OOPSE-4/src/types/MoleculeStamp.cpp (file contents):
Revision 3172 by tim, Wed Jan 11 19:01:20 2006 UTC vs.
Revision 3173 by gezelter, Fri Jul 13 18:10:52 2007 UTC

# Line 46 | Line 46 | namespace oopse {
46   #include "utils/Tuple.hpp"
47   #include "utils/MemoryUtils.hpp"
48   namespace oopse {
49 <
50 < template<class ContainerType>
51 < bool hasDuplicateElement(const ContainerType& cont) {
49 >  
50 >  template<class ContainerType>
51 >  bool hasDuplicateElement(const ContainerType& cont) {
52      ContainerType tmp = cont;
53      std::sort(tmp.begin(), tmp.end());
54      tmp.erase(std::unique(tmp.begin(), tmp.end()), tmp.end());
55      return tmp.size() != cont.size();
56 < }
57 <
58 < MoleculeStamp::MoleculeStamp() {
56 >  }
57 >  
58 >  MoleculeStamp::MoleculeStamp() {
59      DefineParameter(Name, "name");
60      
61      deprecatedKeywords_.insert("nAtoms");
# Line 65 | Line 65 | MoleculeStamp::MoleculeStamp() {
65      deprecatedKeywords_.insert("nRigidBodies");
66      deprecatedKeywords_.insert("nCutoffGroups");
67      
68 < }
69 <
70 < MoleculeStamp::~MoleculeStamp() {
68 >  }
69 >  
70 >  MoleculeStamp::~MoleculeStamp() {
71      MemoryUtils::deletePointers(atomStamps_);
72      MemoryUtils::deletePointers(bondStamps_);
73      MemoryUtils::deletePointers(bendStamps_);
# Line 75 | Line 75 | MoleculeStamp::~MoleculeStamp() {
75      MemoryUtils::deletePointers(rigidBodyStamps_);
76      MemoryUtils::deletePointers(cutoffGroupStamps_);
77      MemoryUtils::deletePointers(fragmentStamps_);    
78 < }
79 <
80 < bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
78 >  }
79 >  
80 >  bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
81      bool ret = addIndexSensitiveStamp(atomStamps_, atom);
82      if (!ret) {
83 <         std::ostringstream oss;
84 <         oss<< "Error in Molecule " << getName()  << ": multiple atoms have the same indices"<< atom->getIndex() <<"\n";
85 <         throw OOPSEException(oss.str());
83 >      std::ostringstream oss;
84 >      oss<< "Error in Molecule " << getName()  <<
85 >        ": multiple atoms have the same indices"<< atom->getIndex() <<"\n";
86 >      throw OOPSEException(oss.str());
87      }
88      return ret;
89      
90 < }
91 <
92 < bool MoleculeStamp::addBondStamp( BondStamp* bond) {
90 >  }
91 >  
92 >  bool MoleculeStamp::addBondStamp( BondStamp* bond) {
93      bondStamps_.push_back(bond);
94      return true;
95 < }
96 <
97 < bool MoleculeStamp::addBendStamp( BendStamp* bend) {
95 >  }
96 >  
97 >  bool MoleculeStamp::addBendStamp( BendStamp* bend) {
98      bendStamps_.push_back(bend);
99      return true;
100 < }
101 <
102 < bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
100 >  }
101 >  
102 >  bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
103      torsionStamps_.push_back(torsion);
104      return true;
105 < }
106 <
107 < bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
105 >  }
106 >  
107 >  bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
108      bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
109      if (!ret) {
110 <        std::ostringstream oss;
111 <        oss<< "Error in Molecule " << getName()  << ": multiple rigidbodies have the same indices: " << rigidbody->getIndex() <<"\n";
112 <        throw OOPSEException(oss.str());
110 >      std::ostringstream oss;
111 >      oss << "Error in Molecule " << getName()  <<
112 >        ": multiple rigidbodies have the same indices: " <<
113 >        rigidbody->getIndex() <<"\n";
114 >      throw OOPSEException(oss.str());
115      }
116      return ret;
117 < }
118 <
119 < bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
117 >  }
118 >  
119 >  bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
120      cutoffGroupStamps_.push_back(cutoffgroup);
121      return true;
122 < }
123 <
124 < bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
122 >  }
123 >  
124 >  bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
125      return addIndexSensitiveStamp(fragmentStamps_, fragment);
126 < }
127 <    
128 < void MoleculeStamp::validate() {
126 >  }
127 >  
128 >  void MoleculeStamp::validate() {
129      DataHolder::validate();
130 <
130 >    
131      atom2Rigidbody.resize(getNAtoms());
132 <    // negative number means atom is a free atom, does not belong to rigidbody
133 <    //every element in atom2Rigidbody has unique negative number at the very beginning
132 >
133 >    // A negative number means the atom is a free atom, and does not
134 >    // belong to rigidbody. Every element in atom2Rigidbody has unique
135 >    // negative number at the very beginning
136 >
137      for(int i = 0; i < atom2Rigidbody.size(); ++i) {
138 <        atom2Rigidbody[i] = -1 - i;
138 >      atom2Rigidbody[i] = -1 - i;
139      }
140      for (int i = 0; i < getNRigidBodies(); ++i) {
141 <        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
142 <        std::vector<int> members = rbStamp->getMembers();
143 <        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
144 <            atom2Rigidbody[*j] = i;                
145 <        }
141 >      RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
142 >      std::vector<int> members = rbStamp->getMembers();
143 >      for(std::vector<int>::iterator j = members.begin();
144 >          j != members.end(); ++j) {
145 >        atom2Rigidbody[*j] = i;                
146 >      }
147      }
148 <
148 >    
149      checkAtoms();
150      checkBonds();
151      fillBondInfo();
# Line 147 | Line 154 | void MoleculeStamp::validate() {
154      checkRigidBodies();
155      checkCutoffGroups();
156      checkFragments();
157 <
157 >    
158      int nrigidAtoms = 0;
159      for (int i = 0; i < getNRigidBodies(); ++i) {
160 <        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
161 <        nrigidAtoms += rbStamp->getNMembers();
160 >      RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
161 >      nrigidAtoms += rbStamp->getNMembers();
162      }
163      nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
164 <
165 < }
166 <
167 < void MoleculeStamp::checkAtoms() {
168 <    std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
164 >    
165 >  }
166 >  
167 >  void MoleculeStamp::checkAtoms() {
168 >    std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(),
169 >                                                     atomStamps_.end(),
170 >                                                     static_cast<AtomStamp*>(NULL));
171      if (ai != atomStamps_.end()) {
172 <        std::ostringstream oss;
173 <        oss << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
174 <        throw OOPSEException(oss.str());
172 >      std::ostringstream oss;
173 >      oss << "Error in Molecule " << getName() << ": atom[" <<
174 >        ai - atomStamps_.begin()<< "] is missing\n";
175 >      throw OOPSEException(oss.str());
176      }
177 +    
178 +  }
179  
180 < }
169 <
170 < void MoleculeStamp::checkBonds() {
180 >  void MoleculeStamp::checkBonds() {
181      std::ostringstream oss;
182      //make sure index is not out of range
183      int natoms = getNAtoms();
184      for(int i = 0; i < getNBonds(); ++i) {
185 <        BondStamp* bondStamp = getBondStamp(i);
186 <        if (bondStamp->getA() > natoms-1 ||  bondStamp->getA() < 0 || bondStamp->getB() > natoms-1 || bondStamp->getB() < 0 || bondStamp->getA() == bondStamp->getB()) {
187 <            
188 <            oss << "Error in Molecule " << getName() <<  ": bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") is invalid\n";
189 <            throw OOPSEException(oss.str());
190 <        }
185 >      BondStamp* bondStamp = getBondStamp(i);
186 >      if (bondStamp->getA() > natoms-1 ||  bondStamp->getA() < 0 ||
187 >          bondStamp->getB() > natoms-1 || bondStamp->getB() < 0 ||
188 >          bondStamp->getA() == bondStamp->getB()) {
189 >        
190 >        oss << "Error in Molecule " << getName() <<  ": bond(" <<
191 >          bondStamp->getA() << ", " << bondStamp->getB() << ") is invalid\n";
192 >        throw OOPSEException(oss.str());
193 >      }
194      }
195      
196      //make sure bonds are unique
197      std::set<std::pair<int, int> > allBonds;
198      for(int i = 0; i < getNBonds(); ++i) {
199 <        BondStamp* bondStamp= getBondStamp(i);        
200 <        std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
201 <        //make sure bondTuple.first is always less than or equal to bondTuple.third
202 <        if (bondPair.first > bondPair.second) {
203 <            std::swap(bondPair.first, bondPair.second);
204 <        }
199 >      BondStamp* bondStamp= getBondStamp(i);        
200 >      std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
201 >      //make sure bondTuple.first is always less than or equal to
202 >      //bondTuple.third
203 >      if (bondPair.first > bondPair.second) {
204 >        std::swap(bondPair.first, bondPair.second);
205 >      }
206 >      
207 >      std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
208 >      if ( iter != allBonds.end()) {
209          
210 <        std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
211 <        if ( iter != allBonds.end()) {
212 <            
213 <            oss << "Error in Molecule " << getName() << ": " << "bond(" <<iter->first << ", "<< iter->second << ") appears multiple times\n";
214 <            throw OOPSEException(oss.str());
215 <        } else {
199 <            allBonds.insert(bondPair);
200 <        }
210 >        oss << "Error in Molecule " << getName() << ": " << "bond(" <<
211 >          iter->first << ", "<< iter->second << ") appears multiple times\n";
212 >        throw OOPSEException(oss.str());
213 >      } else {
214 >        allBonds.insert(bondPair);
215 >      }
216      }
217      
218      //make sure atoms belong to same rigidbody do not bond to each other
219      for(int i = 0; i < getNBonds(); ++i) {
220 <        BondStamp* bondStamp = getBondStamp(i);
221 <        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
222 <            
223 <            oss << "Error in Molecule " << getName() << ": "<<"bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
224 <            throw OOPSEException(oss.str());
225 <        }
226 <    }
227 <    
228 < }
229 <
230 < void MoleculeStamp::checkBends() {
220 >      BondStamp* bondStamp = getBondStamp(i);
221 >      if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
222 >        
223 >        oss << "Error in Molecule " << getName() << ": "<<"bond(" <<
224 >          bondStamp->getA() << ", " << bondStamp->getB() <<
225 >          ") belong to same rigidbody " <<
226 >          atom2Rigidbody[bondStamp->getA()] << "\n";
227 >        throw OOPSEException(oss.str());
228 >      }
229 >    }    
230 >  }
231 >  
232 >  void MoleculeStamp::checkBends() {
233      std::ostringstream oss;
234      for(int i = 0; i < getNBends(); ++i) {
235 <        BendStamp* bendStamp = getBendStamp(i);
236 <        std::vector<int> bendAtoms =  bendStamp->getMembers();
237 <        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
238 <        std::vector<int>::iterator k =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::less<int>(), 0));
239 <
240 <        if (j != bendAtoms.end() || k != bendAtoms.end()) {
241 <            
242 <            oss << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << " have invalid indices\n";
243 <            throw OOPSEException(oss.str());
244 <        }
235 >      BendStamp* bendStamp = getBendStamp(i);
236 >      std::vector<int> bendAtoms =  bendStamp->getMembers();
237 >      std::vector<int>::iterator j =std::find_if(bendAtoms.begin(),
238 >                                                 bendAtoms.end(),
239 >                                                 std::bind2nd(std::greater<int>(),
240 >                                                              getNAtoms()-1));
241 >      std::vector<int>::iterator k =std::find_if(bendAtoms.begin(),
242 >                                                 bendAtoms.end(),
243 >                                                 std::bind2nd(std::less<int>(),
244 >                                                              0));
245 >      
246 >      if (j != bendAtoms.end() || k != bendAtoms.end()) {
247          
248 <        if (hasDuplicateElement(bendAtoms)) {
249 <            oss << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << " have duplicated indices\n";    
250 <            throw OOPSEException(oss.str());            
251 <        }
252 <            
253 <        if (bendAtoms.size() == 2 ) {
254 <            if (!bendStamp->haveGhostVectorSource()) {
255 <                
256 <                oss << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
257 <                throw OOPSEException(oss.str());
258 <            }else{
259 <                int ghostIndex = bendStamp->getGhostVectorSource();
260 <                if (ghostIndex < getNAtoms()) {
261 <                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
262 <                      
263 <                      oss <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
264 <                      throw OOPSEException(oss.str());
265 <                    }
266 <                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
267 <                        
268 <                        oss <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
269 <                        throw OOPSEException(oss.str());
270 <                    }
271 <                }else {
272 <                    oss << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
273 <                    throw OOPSEException(oss.str());
255 <                }
248 >        oss << "Error in Molecule " << getName() << " : atoms of bend" <<
249 >          containerToString(bendAtoms) << " have invalid indices\n";
250 >        throw OOPSEException(oss.str());
251 >      }
252 >      
253 >      if (hasDuplicateElement(bendAtoms)) {
254 >        oss << "Error in Molecule " << getName() << " : atoms of bend" <<
255 >          containerToString(bendAtoms) << " have duplicated indices\n";    
256 >        throw OOPSEException(oss.str());            
257 >      }
258 >      
259 >      if (bendAtoms.size() == 2 ) {
260 >        if (!bendStamp->haveGhostVectorSource()) {
261 >          
262 >          oss << "Error in Molecule " << getName() <<
263 >            ": ghostVectorSouce is missing\n";
264 >          throw OOPSEException(oss.str());
265 >        }else{
266 >          int ghostIndex = bendStamp->getGhostVectorSource();
267 >          if (ghostIndex < getNAtoms()) {
268 >            if (std::find(bendAtoms.begin(), bendAtoms.end(),
269 >                          ghostIndex) == bendAtoms.end()) {
270 >              
271 >              oss <<  "Error in Molecule " << getName() <<
272 >                ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
273 >              throw OOPSEException(oss.str());
274              }
275 <        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
276 <            oss <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
275 >            if (!getAtomStamp(ghostIndex)->haveOrientation()) {
276 >              
277 >              oss <<  "Error in Molecule " << getName() <<
278 >                ": ghost atom must be a directioanl atom\n";
279 >              throw OOPSEException(oss.str());
280 >            }
281 >          } else {
282 >            oss << "Error in Molecule " << getName() <<  
283 >              ": ghostVectorSource " << ghostIndex<< "  is invalid\n";
284              throw OOPSEException(oss.str());
285 +          }
286          }
287 +      } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
288 +        oss <<  "Error in Molecule " << getName() <<
289 +          ": normal bend should not have ghostVectorSouce\n";
290 +        throw OOPSEException(oss.str());
291 +      }
292      }
293 <
293 >    
294      for(int i = 0; i < getNBends(); ++i) {
295 <        BendStamp* bendStamp = getBendStamp(i);
296 <        std::vector<int> bendAtoms =  bendStamp->getMembers();
297 <        std::vector<int> rigidSet(getNRigidBodies(), 0);
298 <        std::vector<int>::iterator j;
299 <        for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
300 <            int rigidbodyIndex = atom2Rigidbody[*j];
301 <            if (rigidbodyIndex >= 0) {
302 <                ++rigidSet[rigidbodyIndex];
303 <                if (rigidSet[rigidbodyIndex] > 1) {
304 <                    oss << "Error in Molecule " << getName() << ": bend" << containerToString(bendAtoms) << " belong to same rigidbody " << rigidbodyIndex << "\n";  
305 <                    throw OOPSEException(oss.str());
306 <                }
307 <            }
295 >      BendStamp* bendStamp = getBendStamp(i);
296 >      std::vector<int> bendAtoms =  bendStamp->getMembers();
297 >      std::vector<int> rigidSet(getNRigidBodies(), 0);
298 >      std::vector<int>::iterator j;
299 >      for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
300 >        int rigidbodyIndex = atom2Rigidbody[*j];
301 >        if (rigidbodyIndex >= 0) {
302 >          ++rigidSet[rigidbodyIndex];
303 >          if (rigidSet[rigidbodyIndex] > 1) {
304 >            oss << "Error in Molecule " << getName() << ": bend" <<
305 >              containerToString(bendAtoms) << " belong to same rigidbody " <<
306 >              rigidbodyIndex << "\n";  
307 >            throw OOPSEException(oss.str());
308 >          }
309          }
310 +      }
311      }
312      
313      
314      std::set<IntTuple3> allBends;
315      std::set<IntTuple3>::iterator iter;
316      for(int i = 0; i < getNBends(); ++i) {
317 <        BendStamp* bendStamp= getBendStamp(i);
318 <        std::vector<int> bend = bendStamp->getMembers();
319 <        if (bend.size() == 2) {
317 >      BendStamp* bendStamp= getBendStamp(i);
318 >      std::vector<int> bend = bendStamp->getMembers();
319 >      if (bend.size() == 2) {
320          // in case we have two ghost bend. For example,
321          // bend {
322          // members (0, 1);
# Line 295 | Line 328 | void MoleculeStamp::checkBends() {
328          // ghostVectorSource = 0;
329          // }
330          // In order to distinguish them. we expand them to Tuple3.
331 <        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
332 <             int ghostIndex = bendStamp->getGhostVectorSource();
333 <             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
334 <             if (j != bend.end()) {
335 <                bend.insert(j, ghostIndex);
336 <             }
331 >        // the first one is expanded to (0, 0, 1) while the second one
332 >        // is expaned to (0, 1, 1)
333 >        int ghostIndex = bendStamp->getGhostVectorSource();
334 >        std::vector<int>::iterator j = std::find(bend.begin(), bend.end(),
335 >                                                 ghostIndex);
336 >        if (j != bend.end()) {
337 >          bend.insert(j, ghostIndex);
338          }
339 +      }
340 +      
341 +      IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
342 +
343 +      // make sure bendTuple.first is always less than or equal to
344 +      // bendTuple.third
345 +      if (bendTuple.first > bendTuple.third) {
346 +        std::swap(bendTuple.first, bendTuple.third);
347 +      }
348 +      
349 +      iter = allBends.find(bendTuple);
350 +      if ( iter != allBends.end()) {
351 +        oss << "Error in Molecule " << getName() << ": " << "Bend" <<
352 +          containerToString(bend)<< " appears multiple times\n";
353 +        throw OOPSEException(oss.str());
354 +      } else {
355 +        allBends.insert(bendTuple);
356 +      }
357 +    }
358 +    
359 +    for (int i = 0; i < getNBonds(); ++i) {
360 +      BondStamp* bondStamp = getBondStamp(i);
361 +      int a = bondStamp->getA();
362 +      int b = bondStamp->getB();
363 +      
364 +      AtomStamp* atomA = getAtomStamp(a);
365 +      AtomStamp* atomB = getAtomStamp(b);
366 +      
367 +      //find bend c--a--b
368 +      AtomStamp::AtomIter ai;
369 +      for(int c= atomA->getFirstBondedAtom(ai);c != -1;
370 +          c = atomA->getNextBondedAtom(ai)) {
371 +        if(b == c)
372 +          continue;          
373          
374 <        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
375 <        //make sure bendTuple.first is always less than or equal to bendTuple.third
376 <        if (bendTuple.first > bendTuple.third) {
309 <            std::swap(bendTuple.first, bendTuple.third);
374 >        IntTuple3 newBend(c, a, b);
375 >        if (newBend.first > newBend.third) {
376 >          std::swap(newBend.first, newBend.third);
377          }
378          
379 <        iter = allBends.find(bendTuple);
380 <        if ( iter != allBends.end()) {
381 <            oss << "Error in Molecule " << getName() << ": " << "Bend" << containerToString(bend)<< " appears multiple times\n";
382 <            throw OOPSEException(oss.str());
383 <        } else {
317 <            allBends.insert(bendTuple);
379 >        if (allBends.find(newBend) == allBends.end() ) {                
380 >          allBends.insert(newBend);
381 >          BendStamp * newBendStamp = new BendStamp();
382 >          newBendStamp->setMembers(newBend);
383 >          addBendStamp(newBendStamp);
384          }
385 <    }
386 <
387 <    for (int i = 0; i < getNBonds(); ++i) {
388 <        BondStamp* bondStamp = getBondStamp(i);
389 <        int a = bondStamp->getA();
390 <        int b = bondStamp->getB();
391 <
392 <        AtomStamp* atomA = getAtomStamp(a);
393 <        AtomStamp* atomB = getAtomStamp(b);
394 <
395 <        //find bend c--a--b
396 <        AtomStamp::AtomIter ai;
397 <        for(int c= atomA->getFirstBondedAtom(ai);c != -1;c = atomA->getNextBondedAtom(ai))
398 <        {
399 <            if(b == c)
400 <                continue;          
401 <            
402 <            IntTuple3 newBend(c, a, b);
403 <            if (newBend.first > newBend.third) {
404 <                std::swap(newBend.first, newBend.third);
405 <            }
406 <
407 <            if (allBends.find(newBend) == allBends.end() ) {                
342 <                allBends.insert(newBend);
343 <                BendStamp * newBendStamp = new BendStamp();
344 <                newBendStamp->setMembers(newBend);
345 <                addBendStamp(newBendStamp);
346 <            }
347 <        }        
348 <
349 <        //find bend a--b--c
350 <        for(int c= atomB->getFirstBondedAtom(ai);c != -1;c = atomB->getNextBondedAtom(ai))
351 <        {
352 <            if(a == c)
353 <                continue;          
354 <
355 <            IntTuple3 newBend( a, b, c);
356 <            if (newBend.first > newBend.third) {
357 <                std::swap(newBend.first, newBend.third);
358 <            }            
359 <            if (allBends.find(newBend) == allBends.end() ) {                
360 <                allBends.insert(newBend);
361 <                BendStamp * newBendStamp = new BendStamp();
362 <                newBendStamp->setMembers(newBend);
363 <                addBendStamp(newBendStamp);
364 <            }
365 <        }        
366 <    }
367 <
368 < }
369 <
370 < void MoleculeStamp::checkTorsions() {
385 >      }        
386 >      
387 >      //find bend a--b--c
388 >      for(int c= atomB->getFirstBondedAtom(ai);c != -1;
389 >          c = atomB->getNextBondedAtom(ai)) {
390 >        if(a == c)
391 >          continue;          
392 >        
393 >        IntTuple3 newBend( a, b, c);
394 >        if (newBend.first > newBend.third) {
395 >          std::swap(newBend.first, newBend.third);
396 >        }            
397 >        if (allBends.find(newBend) == allBends.end() ) {                
398 >          allBends.insert(newBend);
399 >          BendStamp * newBendStamp = new BendStamp();
400 >          newBendStamp->setMembers(newBend);
401 >          addBendStamp(newBendStamp);
402 >        }
403 >      }        
404 >    }    
405 >  }
406 >  
407 >  void MoleculeStamp::checkTorsions() {
408      std::ostringstream oss;
409      for(int i = 0; i < getNTorsions(); ++i) {
410 <        TorsionStamp* torsionStamp = getTorsionStamp(i);
411 <        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
412 <        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
413 <        std::vector<int>::iterator k =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::less<int>(), 0));
414 <
415 <        if (j != torsionAtoms.end() || k != torsionAtoms.end()) {
416 <            oss << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n";
417 <            throw OOPSEException(oss.str());
418 <        }
419 <        if (hasDuplicateElement(torsionAtoms)) {
420 <            oss << "Error in Molecule " << getName() << " : atoms of torsion" << containerToString(torsionAtoms) << " have duplicated indices\n";    
421 <            throw OOPSEException(oss.str());            
422 <        }        
410 >      TorsionStamp* torsionStamp = getTorsionStamp(i);
411 >      std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
412 >      std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(),
413 >                                                 torsionAtoms.end(),
414 >                                                 std::bind2nd(std::greater<int>(),
415 >                                                              getNAtoms()-1));
416 >      std::vector<int>::iterator k =std::find_if(torsionAtoms.begin(),
417 >                                                 torsionAtoms.end(),
418 >                                                 std::bind2nd(std::less<int>(), 0));
419 >      
420 >      if (j != torsionAtoms.end() || k != torsionAtoms.end()) {
421 >        oss << "Error in Molecule " << getName() << ": atoms of torsion" <<
422 >          containerToString(torsionAtoms) << " have invalid indices\n";
423 >        throw OOPSEException(oss.str());
424 >      }
425 >      if (hasDuplicateElement(torsionAtoms)) {
426 >        oss << "Error in Molecule " << getName() << " : atoms of torsion" <<
427 >          containerToString(torsionAtoms) << " have duplicated indices\n";    
428 >        throw OOPSEException(oss.str());            
429 >      }        
430      }
431      
432      for(int i = 0; i < getNTorsions(); ++i) {
433 <        TorsionStamp* torsionStamp = getTorsionStamp(i);
434 <        std::vector<int> torsionAtoms =  torsionStamp->getMembers();
435 <        std::vector<int> rigidSet(getNRigidBodies(), 0);
436 <        std::vector<int>::iterator j;
437 <        for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
438 <            int rigidbodyIndex = atom2Rigidbody[*j];
439 <            if (rigidbodyIndex >= 0) {
440 <                ++rigidSet[rigidbodyIndex];
441 <                if (rigidSet[rigidbodyIndex] > 1) {
442 <                    oss << "Error in Molecule " << getName() << ": torsion" << containerToString(torsionAtoms) << "is invalid\n";          
443 <                    throw OOPSEException(oss.str());
444 <                }
445 <            }
433 >      TorsionStamp* torsionStamp = getTorsionStamp(i);
434 >      std::vector<int> torsionAtoms =  torsionStamp->getMembers();
435 >      std::vector<int> rigidSet(getNRigidBodies(), 0);
436 >      std::vector<int>::iterator j;
437 >      for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
438 >        int rigidbodyIndex = atom2Rigidbody[*j];
439 >        if (rigidbodyIndex >= 0) {
440 >          ++rigidSet[rigidbodyIndex];
441 >          if (rigidSet[rigidbodyIndex] > 1) {
442 >            oss << "Error in Molecule " << getName() << ": torsion" <<
443 >              containerToString(torsionAtoms) << "is invalid\n";          
444 >            throw OOPSEException(oss.str());
445 >          }
446          }
447 +      }
448      }    
449 <
449 >    
450      std::set<IntTuple4> allTorsions;
451      std::set<IntTuple4>::iterator iter;
452 <     for(int i = 0; i < getNTorsions(); ++i) {
453 <         TorsionStamp* torsionStamp= getTorsionStamp(i);
454 <         std::vector<int> torsion = torsionStamp->getMembers();
455 <         if (torsion.size() == 3) {
456 <             int ghostIndex = torsionStamp->getGhostVectorSource();
457 <             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
458 <             if (j != torsion.end()) {
459 <                torsion.insert(j, ghostIndex);
460 <             }
416 <         }
417 <
418 <        IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
419 <        if (torsionTuple.first > torsionTuple.fourth) {
420 <            std::swap(torsionTuple.first, torsionTuple.fourth);
421 <            std::swap(torsionTuple.second, torsionTuple.third);                    
422 <        }                
423 <
424 <         iter = allTorsions.find(torsionTuple);
425 <         if ( iter == allTorsions.end()) {
426 <            allTorsions.insert(torsionTuple);
427 <         } else {
428 <            oss << "Error in Molecule " << getName() << ": " << "Torsion" << containerToString(torsion)<< " appears multiple times\n";
429 <            throw OOPSEException(oss.str());
430 <         }
431 <     }
432 <
433 <    for (int i = 0; i < getNBonds(); ++i) {
434 <        BondStamp* bondStamp = getBondStamp(i);
435 <        int b = bondStamp->getA();
436 <        int c = bondStamp->getB();
437 <
438 <        AtomStamp* atomB = getAtomStamp(b);
439 <        AtomStamp* atomC = getAtomStamp(c);
440 <
441 <        AtomStamp::AtomIter ai2;
442 <        AtomStamp::AtomIter ai3;
443 <
444 <        for(int a = atomB->getFirstBondedAtom(ai2);a != -1;a = atomB->getNextBondedAtom(ai2))
445 <        {
446 <            if(a == c)
447 <                continue;
448 <
449 <            for(int d = atomC->getFirstBondedAtom(ai3);d != -1;d = atomC->getNextBondedAtom(ai3))
450 <            {
451 <                if(d == b)
452 <                    continue;
453 <                
454 <                IntTuple4 newTorsion(a, b, c, d);
455 <                //make sure the first element is always less than or equal to the fourth element in IntTuple4
456 <                if (newTorsion.first > newTorsion.fourth) {
457 <                    std::swap(newTorsion.first, newTorsion.fourth);
458 <                    std::swap(newTorsion.second, newTorsion.third);                    
459 <                }                
460 <                if (allTorsions.find(newTorsion) == allTorsions.end() ) {                
461 <                    allTorsions.insert(newTorsion);
462 <                    TorsionStamp * newTorsionStamp = new TorsionStamp();
463 <                    newTorsionStamp->setMembers(newTorsion);
464 <                    addTorsionStamp(newTorsionStamp);                    
465 <                }            
466 <            }
467 <        }    
468 <    }
469 <
470 < }
471 <
472 < void MoleculeStamp::checkRigidBodies() {
473 <     std::ostringstream oss;
474 <     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
475 <     if (ri != rigidBodyStamps_.end()) {
476 <         oss << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
477 <         throw OOPSEException(oss.str());
478 <     }
479 <
480 <    for (int i = 0; i < getNRigidBodies(); ++i) {
481 <        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
482 <        std::vector<int> rigidAtoms =  rbStamp ->getMembers();
483 <        std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
484 <        if (j != rigidAtoms.end()) {
485 <            oss << "Error in Molecule " << getName();
486 <            throw OOPSEException(oss.str());
452 >    for(int i = 0; i < getNTorsions(); ++i) {
453 >      TorsionStamp* torsionStamp= getTorsionStamp(i);
454 >      std::vector<int> torsion = torsionStamp->getMembers();
455 >      if (torsion.size() == 3) {
456 >        int ghostIndex = torsionStamp->getGhostVectorSource();
457 >        std::vector<int>::iterator j = std::find(torsion.begin(),
458 >                                                 torsion.end(), ghostIndex);
459 >        if (j != torsion.end()) {
460 >          torsion.insert(j, ghostIndex);
461          }
462 <        
463 <    }    
464 < }
465 <
466 < void MoleculeStamp::checkCutoffGroups() {
467 <    for(int i = 0; i < getNCutoffGroups(); ++i) {
468 <        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
469 <        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
470 <        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
471 <        if (j != cutoffGroupAtoms.end()) {
472 <            std::ostringstream oss;
473 <            oss << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
474 <            throw OOPSEException(oss.str());
475 <        }
502 <    }    
503 < }
504 <
505 < void MoleculeStamp::checkFragments() {
506 <
507 <    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
508 <    if (fi != fragmentStamps_.end()) {
509 <        std::ostringstream oss;
510 <        oss << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
462 >      }
463 >      
464 >      IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
465 >      if (torsionTuple.first > torsionTuple.fourth) {
466 >        std::swap(torsionTuple.first, torsionTuple.fourth);
467 >        std::swap(torsionTuple.second, torsionTuple.third);                    
468 >      }                
469 >      
470 >      iter = allTorsions.find(torsionTuple);
471 >      if ( iter == allTorsions.end()) {
472 >        allTorsions.insert(torsionTuple);
473 >      } else {
474 >        oss << "Error in Molecule " << getName() << ": " << "Torsion" <<
475 >          containerToString(torsion)<< " appears multiple times\n";
476          throw OOPSEException(oss.str());
477 +      }
478      }
479      
514 – }
515 –
516 – void MoleculeStamp::fillBondInfo() {
517 –
480      for (int i = 0; i < getNBonds(); ++i) {
481 <        BondStamp* bondStamp = getBondStamp(i);
482 <        int a = bondStamp->getA();
483 <        int b = bondStamp->getB();
484 <        AtomStamp* atomA = getAtomStamp(a);
485 <        AtomStamp* atomB = getAtomStamp(b);
486 <        atomA->addBond(i);
487 <        atomA->addBondedAtom(b);
488 <        atomB->addBond(i);        
489 <        atomB->addBondedAtom(a);
490 <
481 >      BondStamp* bondStamp = getBondStamp(i);
482 >      int b = bondStamp->getA();
483 >      int c = bondStamp->getB();
484 >      
485 >      AtomStamp* atomB = getAtomStamp(b);
486 >      AtomStamp* atomC = getAtomStamp(c);
487 >      
488 >      AtomStamp::AtomIter ai2;
489 >      AtomStamp::AtomIter ai3;
490 >      
491 >      for(int a = atomB->getFirstBondedAtom(ai2);a != -1;
492 >          a = atomB->getNextBondedAtom(ai2)) {
493 >        if(a == c)
494 >          continue;
495 >        
496 >        for(int d = atomC->getFirstBondedAtom(ai3);d != -1;
497 >            d = atomC->getNextBondedAtom(ai3)) {          
498 >          if(d == b)
499 >            continue;
500 >          
501 >          IntTuple4 newTorsion(a, b, c, d);
502 >          //make sure the first element is always less than or equal
503 >          //to the fourth element in IntTuple4
504 >          if (newTorsion.first > newTorsion.fourth) {
505 >            std::swap(newTorsion.first, newTorsion.fourth);
506 >            std::swap(newTorsion.second, newTorsion.third);                    
507 >          }                
508 >          if (allTorsions.find(newTorsion) == allTorsions.end() ) {
509 >            allTorsions.insert(newTorsion);
510 >            TorsionStamp * newTorsionStamp = new TorsionStamp();
511 >            newTorsionStamp->setMembers(newTorsion);
512 >            addTorsionStamp(newTorsionStamp);                    
513 >          }            
514 >        }
515 >      }    
516      }
517 < }
517 >    
518 >  }
519 >  
520 >  void MoleculeStamp::checkRigidBodies() {
521 >    std::ostringstream oss;
522 >    std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(),
523 >                                                          rigidBodyStamps_.end(),
524 >                                                          static_cast<RigidBodyStamp*>(NULL));
525 >    if (ri != rigidBodyStamps_.end()) {
526 >      oss << "Error in Molecule " << getName() << ":rigidBody[" <<  
527 >        ri - rigidBodyStamps_.begin()<< "] is missing\n";
528 >      throw OOPSEException(oss.str());
529 >    }
530 >    
531 >    for (int i = 0; i < getNRigidBodies(); ++i) {
532 >      RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
533 >      std::vector<int> rigidAtoms =  rbStamp ->getMembers();
534 >      std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(),
535 >                                                 rigidAtoms.end(),
536 >                                                 std::bind2nd(std::greater<int>(),
537 >                                                              getNAtoms()-1));
538 >      if (j != rigidAtoms.end()) {
539 >        oss << "Error in Molecule " << getName();
540 >        throw OOPSEException(oss.str());
541 >      }      
542 >    }    
543 >  }
544 >  
545 >  void MoleculeStamp::checkCutoffGroups() {
546 >    for(int i = 0; i < getNCutoffGroups(); ++i) {
547 >      CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
548 >      std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
549 >      std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(),
550 >                                                 cutoffGroupAtoms.end(),
551 >                                                 std::bind2nd(std::greater<int>(),
552 >                                                              getNAtoms()-1));
553 >      if (j != cutoffGroupAtoms.end()) {
554 >        std::ostringstream oss;
555 >        oss << "Error in Molecule " << getName() << ": cutoffGroup" <<
556 >          " is out of range\n";
557 >        throw OOPSEException(oss.str());
558 >      }
559 >    }    
560 >  }
561  
562 +  void MoleculeStamp::checkFragments() {
563  
564 +    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(),
565 +                                                         fragmentStamps_.end(),
566 +                                                         static_cast<FragmentStamp*>(NULL));
567 +    if (fi != fragmentStamps_.end()) {
568 +      std::ostringstream oss;
569 +      oss << "Error in Molecule " << getName() << ":fragment[" <<  
570 +        fi - fragmentStamps_.begin()<< "] is missing\n";
571 +      throw OOPSEException(oss.str());
572 +    }
573 +    
574 +  }
575 +  
576 +  void MoleculeStamp::fillBondInfo() {
577 +    
578 +    for (int i = 0; i < getNBonds(); ++i) {
579 +      BondStamp* bondStamp = getBondStamp(i);
580 +      int a = bondStamp->getA();
581 +      int b = bondStamp->getB();
582 +      AtomStamp* atomA = getAtomStamp(a);
583 +      AtomStamp* atomB = getAtomStamp(b);
584 +      atomA->addBond(i);
585 +      atomA->addBondedAtom(b);
586 +      atomB->addBond(i);        
587 +      atomB->addBondedAtom(a);
588  
589 < //Function Name: isBondInSameRigidBody
590 < //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
536 < bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
537 <  int rbA;
538 <  int rbB;
539 <  int consAtomA;
540 <  int consAtomB;
589 >    }
590 >  }
591  
542 –  if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
543 –    return false;
592  
545 –  if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
546 –    return false;
593  
594 <  if(rbB == rbA)
595 <    return true;
596 <  else
597 <    return false;
598 < }
594 >  // Function Name: isBondInSameRigidBody
595 >  // Returns true is both atoms of the bond belong to the same rigid
596 >  // body, otherwise return false
597 >  bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
598 >    int rbA;
599 >    int rbB;
600 >    int consAtomA;
601 >    int consAtomB;
602  
603 < // Function Name: isAtomInRigidBody
604 < //return false if atom does not belong to a rigid body, otherwise return true
556 < bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
557 <  return atom2Rigidbody[atomIndex] >=0 ;
558 <  
559 < }
603 >    if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
604 >      return false;
605  
606 < // Function Name: isAtomInRigidBody
607 < //return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
563 < //and consAtomIndex
564 < //atomIndex : the index of atom in component
565 < //whichRigidBody: the index of rigidbody in component
566 < //consAtomIndex:  the position of joint atom apears in  rigidbody's definition
567 < bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
606 >    if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
607 >      return false;
608  
609 <  
609 >    if(rbB == rbA)
610 >      return true;
611 >    else
612 >      return false;
613 >  }
614  
615 <  whichRigidBody = -1;
616 <  consAtomIndex = -1;
617 <
618 <  if (atom2Rigidbody[atomIndex] >=0) {
619 <    whichRigidBody = atom2Rigidbody[atomIndex];
620 <    RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
621 <    int numAtom = rbStamp->getNMembers();
622 <    for(int j = 0; j < numAtom; j++) {
623 <      if (rbStamp->getMemberAt(j) == atomIndex){
624 <        consAtomIndex = j;
625 <        return true;
615 >  // Function Name: isAtomInRigidBody
616 >  // Returns false if atom does not belong to a rigid body, otherwise
617 >  // returns true
618 >  bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
619 >    return atom2Rigidbody[atomIndex] >=0 ;    
620 >  }
621 >  
622 >  // Function Name: isAtomInRigidBody
623 >  // Returns false if atom does not belong to a rigid body otherwise
624 >  // returns true and sets whichRigidBody and consAtomIndex
625 >  // atomIndex : the index of atom in component
626 >  // whichRigidBody: the index of the rigidbody in the component
627 >  // consAtomIndex:  the position the joint atom appears in the rigidbody's
628 >  //                 definition
629 >  bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody,
630 >                                        int& consAtomIndex){
631 >    whichRigidBody = -1;
632 >    consAtomIndex = -1;
633 >    
634 >    if (atom2Rigidbody[atomIndex] >=0) {
635 >      whichRigidBody = atom2Rigidbody[atomIndex];
636 >      RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
637 >      int numAtom = rbStamp->getNMembers();
638 >      for(int j = 0; j < numAtom; j++) {
639 >        if (rbStamp->getMemberAt(j) == atomIndex){
640 >          consAtomIndex = j;
641 >          return true;
642 >        }
643        }
644      }
645 +    
646 +    return false;    
647    }
585 –
586 –  return false;
587 –  
588 – }
589 –
590 – //return the position of joint atom apears in  rigidbody's definition
591 – //for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
592 – //actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
593 – std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
594 –  RigidBodyStamp* rbStamp1;
595 –  RigidBodyStamp* rbStamp2;
596 –  int natomInRb1;
597 –  int natomInRb2;
598 –  int atomIndex1;
599 –  int atomIndex2;
600 –  std::vector<std::pair<int, int> > jointAtomIndexPair;
648    
649 <  rbStamp1 = this->getRigidBodyStamp(rb1);
650 <  natomInRb1 =rbStamp1->getNMembers();
649 >  // Returns the position of joint atoms apearing in a rigidbody's definition
650 >  // For the time being, we will use the most inefficient algorithm,
651 >  // the complexity is O(N^2).  We could improve the
652 >  // complexity to O(NlogN) by sorting the atom index in rigid body
653 >  // first
654 >  std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1,
655 >                                                                 int rb2){
656 >    RigidBodyStamp* rbStamp1;
657 >    RigidBodyStamp* rbStamp2;
658 >    int natomInRb1;
659 >    int natomInRb2;
660 >    int atomIndex1;
661 >    int atomIndex2;
662 >    std::vector<std::pair<int, int> > jointAtomIndexPair;
663 >    
664 >    rbStamp1 = this->getRigidBodyStamp(rb1);
665 >    natomInRb1 =rbStamp1->getNMembers();
666 >    
667 >    rbStamp2 = this->getRigidBodyStamp(rb2);
668 >    natomInRb2 =rbStamp2->getNMembers();
669  
670 <  rbStamp2 = this->getRigidBodyStamp(rb2);
671 <  natomInRb2 =rbStamp2->getNMembers();
607 <
608 <  for(int i = 0; i < natomInRb1; i++){
609 <    atomIndex1 = rbStamp1->getMemberAt(i);
670 >    for(int i = 0; i < natomInRb1; i++){
671 >      atomIndex1 = rbStamp1->getMemberAt(i);
672        
673 <    for(int j= 0; j < natomInRb1; j++){
674 <      atomIndex2 = rbStamp2->getMemberAt(j);
673 >      for(int j= 0; j < natomInRb1; j++){
674 >        atomIndex2 = rbStamp2->getMemberAt(j);
675  
676 <      if(atomIndex1 == atomIndex2){
677 <        jointAtomIndexPair.push_back(std::make_pair(i, j));
678 <        break;
679 <      }
676 >        if(atomIndex1 == atomIndex2){
677 >          jointAtomIndexPair.push_back(std::make_pair(i, j));
678 >          break;
679 >        }
680        
681 +      }
682 +
683      }
684  
685 +    return jointAtomIndexPair;
686    }
687  
623 –  return jointAtomIndexPair;
688   }
625 –
626 – }

Diff Legend

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