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

Diff Legend

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