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

Comparing trunk/OOPSE-3.0/src/types/MoleculeStamp.cpp (file contents):
Revision 2469 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 2493 by tim, Tue Dec 6 21:01:08 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 <    for(int i = 0; i < getNBends(); ++i) {
177 <        BendStamp* bendStamp = getBendStamp(i);
178 <        std::vector<int> bendAtoms =  bendStamp->getMembers();
179 <        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
180 <        if (j != bendAtoms.end()) {
181 <            std::cout << "Error in Molecule " << getName();
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 <        if (bendAtoms.size() == 2 && !bendStamp->haveGhostVectorSource()) {
188 <            std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing";
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          }
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 –        }
193      }
194 <    for(int i = 0; i < getNCutoffGroups(); ++i) {
195 <        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
196 <        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
197 <        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
198 <        if (j != cutoffGroupAtoms.end()) {
199 <            std::cout << "Error in Molecule " << getName();
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 <    atom2Rigidbody.resize(natoms);
204 <    // negative number means atom is a free atom, does not belong to rigidbody
205 <    //every element in atom2Rigidbody has unique negative number at the very beginning
206 <    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
207 <        atom2Rigidbody[i] = -1 - i;
202 >    
203 > }
204 >
205 > struct BendLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
206 >    bool operator()(IntTuple3 b1, IntTuple3 b2) {
207 >        return b1.first < b2.first
208 >             || (!(b2.first < b1.first) && b1.second < b2.second)
209 >             || (!(b2.first < b1.first) && !(b2.second < b2.second) && b1.third < b2.third);
210      }
211 + };
212  
213 <    for (int i = 0; i < getNRigidBodies(); ++i) {
214 <        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
215 <        std::vector<int> members = rbStamp->getMembers();
216 <        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
217 <            atom2Rigidbody[*j] = i;                
213 > void MoleculeStamp::checkBends() {
214 >    for(int i = 0; i < getNBends(); ++i) {
215 >        BendStamp* bendStamp = getBendStamp(i);
216 >        std::vector<int> bendAtoms =  bendStamp->getMembers();
217 >        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
218 >        if (j != bendAtoms.end()) {
219 >            std::cout << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << "have invalid indices\n";
220          }
221 <    }
222 <    //make sure atoms belong to same rigidbody do not bond to each other
223 <    for(int i = 0; i < getNBonds(); ++i) {
224 <        BondStamp* bondStamp = getBondStamp(i);
225 <        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()])
226 <            std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << "belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
221 >
222 >        if (bendAtoms.size() == 2 ) {
223 >            if (!bendStamp->haveGhostVectorSource()) {
224 >                std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
225 >            }else{
226 >                int ghostIndex = bendStamp->getGhostVectorSource();
227 >                if (ghostIndex < getNAtoms()) {
228 >                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
229 >                      std::cout <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
230 >                    }
231 >                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
232 >                        std::cout <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
233 >                    }
234 >                }else {
235 >                    std::cout << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
236 >                }
237 >            }
238 >        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
239 >            std::cout <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
240          }
241 <        
241 >    }
242 >
243      for(int i = 0; i < getNBends(); ++i) {
244          BendStamp* bendStamp = getBendStamp(i);
245          std::vector<int> bendAtoms =  bendStamp->getMembers();
# Line 190 | Line 250 | void MoleculeStamp::validate() {
250              if (rigidbodyIndex >= 0) {
251                  ++rigidSet[rigidbodyIndex];
252                  if (rigidSet[rigidbodyIndex] > 1) {
253 <                    std::cout << "Error in Molecule " << getName() << ": ";
194 <                    //std::cout << "atoms of bend " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
253 >                    std::cout << "Error in Molecule " << getName() << ": bend" << containerToString(bendAtoms) << " belong to same rigidbody " << rigidbodyIndex << "\n";                    
254                  }
255              }
256          }
257 <    }      
257 >    }
258 >    
259 >    
260 >    std::set<IntTuple3, BendLessThan> allBends;
261 >    std::set<IntTuple3, BendLessThan>::iterator iter;
262 >    for(int i = 0; i < getNBends(); ++i) {
263 >        BendStamp* bendStamp= getBendStamp(i);
264 >        std::vector<int> bend = bendStamp->getMembers();
265 >        if (bend.size() == 2) {
266 >        // in case we have two ghost bend. For example,
267 >        // bend {
268 >        // members (0, 1);
269 >        //   ghostVectorSource = 0;
270 >        // }
271 >        // and
272 >        // bend {
273 >        //   members (0, 1);
274 >        // ghostVectorSource = 0;
275 >        // }
276 >        // In order to distinguish them. we expand them to Tuple3.
277 >        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
278 >             int ghostIndex = bendStamp->getGhostVectorSource();
279 >             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
280 >             if (j != bend.end()) {
281 >                bend.insert(j, ghostIndex);
282 >             }
283 >        }
284 >        
285 >        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
286 >        //make sure bendTuple.first is always less than or equal to bendTuple.third
287 >        if (bendTuple.first > bendTuple.third) {
288 >            std::swap(bendTuple.first, bendTuple.third);
289 >        }
290 >        
291 >        iter = allBends.find(bendTuple);
292 >        if ( iter != allBends.end()) {
293 >            std::cout << "Error in Molecule " << getName() << ": " << "Bend appears multiple times\n";
294 >        } else {
295 >            allBends.insert(bendTuple);
296 >        }
297 >    }
298 >
299 >    for (int i = 0; i < getNBonds(); ++i) {
300 >        BondStamp* bondStamp = getBondStamp(i);
301 >        int a = bondStamp->getA();
302 >        int b = bondStamp->getB();
303 >
304 >        AtomStamp* atomA = getAtomStamp(a);
305 >        AtomStamp* atomB = getAtomStamp(b);
306 >
307 >        //find bend c--a--b
308 >        AtomStamp::AtomIter ai;
309 >        for(int c= atomA->getFirstBonedAtom(ai);c != -1;c = atomA->getNextBonedAtom(ai))
310 >        {
311 >            if(b == c)
312 >                continue;          
313 >            
314 >            IntTuple3 newBend(c, a, b);
315 >            if (newBend.first > newBend.third) {
316 >                std::swap(newBend.first, newBend.third);
317 >            }
318 >
319 >            if (allBends.find(newBend) == allBends.end() ) {                
320 >                allBends.insert(newBend);
321 >                BendStamp * newBendStamp = new BendStamp();
322 >                newBendStamp->setMembers(newBend);
323 >                addBendStamp(newBendStamp);
324 >            }
325 >        }        
326 >
327 >        //find bend a--b--c
328 >        for(int c= atomB->getFirstBonedAtom(ai);c != -1;c = atomB->getNextBonedAtom(ai))
329 >        {
330 >            if(a == c)
331 >                continue;          
332 >
333 >            IntTuple3 newBend( a, b, c);
334 >            if (newBend.first > newBend.third) {
335 >                std::swap(newBend.first, newBend.third);
336 >            }            
337 >            if (allBends.find(newBend) == allBends.end() ) {                
338 >                allBends.insert(newBend);
339 >                BendStamp * newBendStamp = new BendStamp();
340 >                newBendStamp->setMembers(newBend);
341 >                addBendStamp(newBendStamp);
342 >            }
343 >        }        
344 >    }
345 >
346 > }
347 >
348 > struct TorsionLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
349 >    bool operator()(IntTuple4 t1, IntTuple4 t2) {
350 >
351 >        return t1.first < t2.first
352 >             || (!(t2.first < t1.first) && t1.second < t2.second)
353 >             || (!(t2.first < t1.first) && !(t2.second < t2.second) && t1.third < t2.third)
354 >             ||(!(t2.first < t1.first) && !(t2.second < t2.second) && !(t2.third < t1.third) && t1.fourth < t2.fourth);
355 >    }
356 >
357 >
358 >
359 > };
360 >
361 >
362 > void MoleculeStamp::checkTorsions() {
363      for(int i = 0; i < getNTorsions(); ++i) {
364          TorsionStamp* torsionStamp = getTorsionStamp(i);
365 +        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
366 +        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
367 +        if (j != torsionAtoms.end()) {
368 +            std::cout << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n";
369 +        }
370 +    }
371 +    
372 +    for(int i = 0; i < getNTorsions(); ++i) {
373 +        TorsionStamp* torsionStamp = getTorsionStamp(i);
374          std::vector<int> torsionAtoms =  torsionStamp->getMembers();
375          std::vector<int> rigidSet(getNRigidBodies(), 0);
376          std::vector<int>::iterator j;
# Line 206 | Line 379 | void MoleculeStamp::validate() {
379              if (rigidbodyIndex >= 0) {
380                  ++rigidSet[rigidbodyIndex];
381                  if (rigidSet[rigidbodyIndex] > 1) {
382 <                    std::cout << "Error in Molecule " << getName() << ": ";
210 <                    //std::cout << "atoms of torsion " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
382 >                    std::cout << "Error in Molecule " << getName() << ": torsion" << containerToString(torsionAtoms) << "is invalid\n";                  
383                  }
384              }
385          }
386 <    }
386 >    }    
387  
388 +    std::set<IntTuple4, TorsionLessThan> allTorsions;
389 +    std::set<IntTuple4, TorsionLessThan>::iterator iter;
390 +     for(int i = 0; i < getNTorsions(); ++i) {
391 +         TorsionStamp* torsionStamp= getTorsionStamp(i);
392 +         std::vector<int> torsion = torsionStamp->getMembers();
393 +         if (torsion.size() == 3) {
394 +             int ghostIndex = torsionStamp->getGhostVectorSource();
395 +             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
396 +             if (j != torsion.end()) {
397 +                torsion.insert(j, ghostIndex);
398 +             }
399 +         }
400  
401 <    //fill in bond information into atom
402 <    fillBondInfo();
403 <    findBends();
404 <    findTorsions();
401 >        IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
402 >        if (torsionTuple.first > torsionTuple.fourth) {
403 >            std::swap(torsionTuple.first, torsionTuple.fourth);
404 >            std::swap(torsionTuple.second, torsionTuple.third);                    
405 >        }                
406  
407 <    int nrigidAtoms = 0;
407 >         iter = allTorsions.find(torsionTuple);
408 >         if ( iter == allTorsions.end()) {
409 >            allTorsions.insert(torsionTuple);
410 >         } else {
411 >            std::cout << "Error in Molecule " << getName() << ": " << "Torsion appears multiple times\n";
412 >         }
413 >     }
414 >
415 >    for (int i = 0; i < getNBonds(); ++i) {
416 >        BondStamp* bondStamp = getBondStamp(i);
417 >        int b = bondStamp->getA();
418 >        int c = bondStamp->getB();
419 >
420 >        AtomStamp* atomB = getAtomStamp(b);
421 >        AtomStamp* atomC = getAtomStamp(c);
422 >
423 >        AtomStamp::AtomIter ai2;
424 >        AtomStamp::AtomIter ai3;
425 >
426 >        for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
427 >        {
428 >            if(a == c)
429 >                continue;
430 >
431 >            for(int d = atomC->getFirstBonedAtom(ai3);d != -1;d = atomC->getNextBonedAtom(ai3))
432 >            {
433 >                if(d == b)
434 >                    continue;
435 >                
436 >                IntTuple4 newTorsion(a, b, c, d);
437 >                //make sure the first element is always less than or equal to the fourth element in IntTuple4
438 >                if (newTorsion.first > newTorsion.fourth) {
439 >                    std::swap(newTorsion.first, newTorsion.fourth);
440 >                    std::swap(newTorsion.second, newTorsion.third);                    
441 >                }                
442 >                if (allTorsions.find(newTorsion) == allTorsions.end() ) {                
443 >                    allTorsions.insert(newTorsion);
444 >                    TorsionStamp * newTorsionStamp = new TorsionStamp();
445 >                    newTorsionStamp->setMembers(newTorsion);
446 >                    addTorsionStamp(newTorsionStamp);                    
447 >                }            
448 >            }
449 >        }    
450 >    }
451 >
452 > }
453 >
454 > void MoleculeStamp::checkRigidBodies() {
455 >     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
456 >     if (ri != rigidBodyStamps_.end()) {
457 >         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
458 >     }
459 >
460      for (int i = 0; i < getNRigidBodies(); ++i) {
461          RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
462 <        nrigidAtoms += rbStamp->getNMembers();
463 <    }
464 <    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
462 >        std::vector<int> rigidAtoms =  rbStamp ->getMembers();
463 >        std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
464 >        if (j != rigidAtoms.end()) {
465 >            std::cout << "Error in Molecule " << getName();
466 >        }
467 >        
468 >    }    
469 > }
470  
471 < }
471 > void MoleculeStamp::checkCutoffGroups() {
472  
473 < void MoleculeStamp::fillBondInfo() {
473 >    for(int i = 0; i < getNCutoffGroups(); ++i) {
474 >        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
475 >        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
476 >        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
477 >        if (j != cutoffGroupAtoms.end()) {
478 >            std::cout << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
479 >        }
480 >    }    
481   }
482  
483 < void MoleculeStamp::findBends() {
483 > void MoleculeStamp::checkFragments() {
484  
485 +    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
486 +    if (fi != fragmentStamps_.end()) {
487 +        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
488 +    }
489 +    
490   }
491  
492 < void MoleculeStamp::findTorsions() {
492 > void MoleculeStamp::fillBondInfo() {
493  
494 +    for (int i = 0; i < getNBonds(); ++i) {
495 +        BondStamp* bondStamp = getBondStamp(i);
496 +        int a = bondStamp->getA();
497 +        int b = bondStamp->getB();
498 +        AtomStamp* atomA = getAtomStamp(a);
499 +        AtomStamp* atomB = getAtomStamp(b);
500 +        atomA->addBond(i);
501 +        atomA->addBondedAtom(b);
502 +        atomB->addBond(i);        
503 +        atomB->addBondedAtom(a);
504 +
505 +    }
506   }
507  
508 +
509 +
510   //Function Name: isBondInSameRigidBody
511   //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
512   bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
# Line 262 | Line 530 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
530   // Function Name: isAtomInRigidBody
531   //return false if atom does not belong to a rigid body, otherwise return true
532   bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
533 <  int whichRigidBody;
266 <  int consAtomIndex;
267 <
268 <  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
533 >  return atom2Rigidbody[atomIndex] >=0 ;
534    
535   }
536  
# Line 276 | Line 541 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex, i
541   //whichRigidBody: the index of rigidbody in component
542   //consAtomIndex:  the position of joint atom apears in  rigidbody's definition
543   bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
279 –  RigidBodyStamp* rbStamp;
280 –  int numRb;
281 –  int numAtom;
544  
545 +  
546 +
547    whichRigidBody = -1;
548    consAtomIndex = -1;
549  
550 <  numRb = this->getNRigidBodies();
551 <  
552 <  for(int i = 0 ; i < numRb; i++){
553 <    rbStamp = this->getRigidBodyStamp(i);
554 <    numAtom = rbStamp->getNMembers();
291 <    for(int j = 0; j < numAtom; j++)
550 >  if (atom2Rigidbody[atomIndex] >=0) {
551 >    whichRigidBody = atom2Rigidbody[atomIndex];
552 >    RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
553 >    int numAtom = rbStamp->getNMembers();
554 >    for(int j = 0; j < numAtom; j++) {
555        if (rbStamp->getMemberAt(j) == atomIndex){
293 –        whichRigidBody = i;
556          consAtomIndex = j;
557          return true;
558        }
559 +    }
560    }
561  
562    return false;

Diff Legend

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