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 2544 by tim, Wed Jan 11 19:01:20 2006 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>
41 > #include <algorithm>
42 > #include <functional>
43   #include <iostream>
44 <
44 > #include <sstream>
45   #include "types/MoleculeStamp.hpp"
46 <
46 > #include "utils/Tuple.hpp"
47 > #include "utils/MemoryUtils.hpp"
48   namespace oopse {
49 +
50 + template<class ContainerType>
51 + 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() {
59      DefineParameter(Name, "name");
60      
# Line 60 | Line 68 | MoleculeStamp::~MoleculeStamp() {
68   }
69  
70   MoleculeStamp::~MoleculeStamp() {
71 <
71 >    MemoryUtils::deletePointers(atomStamps_);
72 >    MemoryUtils::deletePointers(bondStamps_);
73 >    MemoryUtils::deletePointers(bendStamps_);
74 >    MemoryUtils::deletePointers(torsionStamps_);
75 >    MemoryUtils::deletePointers(rigidBodyStamps_);
76 >    MemoryUtils::deletePointers(cutoffGroupStamps_);
77 >    MemoryUtils::deletePointers(fragmentStamps_);    
78   }
79  
80   bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
81      bool ret = addIndexSensitiveStamp(atomStamps_, atom);
82      if (!ret) {
83 <        std::cout << "multiple atoms have the same index: " << atom->getIndex() <<" in " << getName()  << " Molecule\n";
83 >         std::ostringstream oss;
84 >         oss<< "Error in Molecule " << getName()  << ": multiple atoms have the same indices"<< atom->getIndex() <<"\n";
85 >         throw OOPSEException(oss.str());
86      }
87      return ret;
88      
# Line 90 | Line 106 | bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp*
106   bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
107      bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
108      if (!ret) {
109 <        std::cout << "multiple rigidbodies have the same index: " << rigidbody->getIndex() <<" in " << getName()  << " Molecule\n";
109 >        std::ostringstream oss;
110 >        oss<< "Error in Molecule " << getName()  << ": multiple rigidbodies have the same indices: " << rigidbody->getIndex() <<"\n";
111 >        throw OOPSEException(oss.str());
112      }
113      return ret;
114   }
# Line 107 | Line 125 | void MoleculeStamp::validate() {
125   void MoleculeStamp::validate() {
126      DataHolder::validate();
127  
128 +    atom2Rigidbody.resize(getNAtoms());
129 +    // negative number means atom is a free atom, does not belong to rigidbody
130 +    //every element in atom2Rigidbody has unique negative number at the very beginning
131 +    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
132 +        atom2Rigidbody[i] = -1 - i;
133 +    }
134 +    for (int i = 0; i < getNRigidBodies(); ++i) {
135 +        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
136 +        std::vector<int> members = rbStamp->getMembers();
137 +        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
138 +            atom2Rigidbody[*j] = i;                
139 +        }
140 +    }
141 +
142 +    checkAtoms();
143 +    checkBonds();
144 +    fillBondInfo();
145 +    checkBends();
146 +    checkTorsions();
147 +    checkRigidBodies();
148 +    checkCutoffGroups();
149 +    checkFragments();
150 +
151 +    int nrigidAtoms = 0;
152 +    for (int i = 0; i < getNRigidBodies(); ++i) {
153 +        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
154 +        nrigidAtoms += rbStamp->getNMembers();
155 +    }
156 +    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
157 +
158 + }
159 +
160 + void MoleculeStamp::checkAtoms() {
161      std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
162      if (ai != atomStamps_.end()) {
163 <        std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
163 >        std::ostringstream oss;
164 >        oss << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
165 >        throw OOPSEException(oss.str());
166      }
167  
168 <     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 <    }
168 > }
169  
170 + void MoleculeStamp::checkBonds() {
171 +    std::ostringstream oss;
172      //make sure index is not out of range
173      int natoms = getNAtoms();
174      for(int i = 0; i < getNBonds(); ++i) {
175          BondStamp* bondStamp = getBondStamp(i);
176 <        if (bondStamp->getA() >=  natoms && bondStamp->getB() >= natoms) {
177 <            std::cout << "Error in Molecule " << getName() <<  ": bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " is invalid\n";
176 >        if (bondStamp->getA() > natoms-1 ||  bondStamp->getA() < 0 || bondStamp->getB() > natoms-1 || bondStamp->getB() < 0 || bondStamp->getA() == bondStamp->getB()) {
177 >            
178 >            oss << "Error in Molecule " << getName() <<  ": bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") is invalid\n";
179 >            throw OOPSEException(oss.str());
180          }
181      }
182 <    for(int i = 0; i < getNBends(); ++i) {
183 <        BendStamp* bendStamp = getBendStamp(i);
184 <        std::vector<int> bendAtoms =  bendStamp->getMembers();
185 <        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
186 <        if (j != bendAtoms.end()) {
187 <            std::cout << "Error in Molecule " << getName();
182 >    
183 >    //make sure bonds are unique
184 >    std::set<std::pair<int, int> > allBonds;
185 >    for(int i = 0; i < getNBonds(); ++i) {
186 >        BondStamp* bondStamp= getBondStamp(i);        
187 >        std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
188 >        //make sure bondTuple.first is always less than or equal to bondTuple.third
189 >        if (bondPair.first > bondPair.second) {
190 >            std::swap(bondPair.first, bondPair.second);
191          }
192 <
193 <        if (bendAtoms.size() == 2 && !bendStamp->haveGhostVectorSource()) {
194 <            std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing";
192 >        
193 >        std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
194 >        if ( iter != allBonds.end()) {
195 >            
196 >            oss << "Error in Molecule " << getName() << ": " << "bond(" <<iter->first << ", "<< iter->second << ") appears multiple times\n";
197 >            throw OOPSEException(oss.str());
198 >        } else {
199 >            allBonds.insert(bondPair);
200          }
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 –        }
201      }
202 <    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 <    }
168 <
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 <    }
202 >    
203      //make sure atoms belong to same rigidbody do not bond to each other
204      for(int i = 0; i < getNBonds(); ++i) {
205          BondStamp* bondStamp = getBondStamp(i);
206 <        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()])
207 <            std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << "belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
206 >        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
207 >            
208 >            oss << "Error in Molecule " << getName() << ": "<<"bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
209 >            throw OOPSEException(oss.str());
210          }
211 +    }
212 +    
213 + }
214 +
215 + void MoleculeStamp::checkBends() {
216 +    std::ostringstream oss;
217 +    for(int i = 0; i < getNBends(); ++i) {
218 +        BendStamp* bendStamp = getBendStamp(i);
219 +        std::vector<int> bendAtoms =  bendStamp->getMembers();
220 +        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
221 +        std::vector<int>::iterator k =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::less<int>(), 0));
222 +
223 +        if (j != bendAtoms.end() || k != bendAtoms.end()) {
224 +            
225 +            oss << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << " have invalid indices\n";
226 +            throw OOPSEException(oss.str());
227 +        }
228          
229 +        if (hasDuplicateElement(bendAtoms)) {
230 +            oss << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << " have duplicated indices\n";    
231 +            throw OOPSEException(oss.str());            
232 +        }
233 +            
234 +        if (bendAtoms.size() == 2 ) {
235 +            if (!bendStamp->haveGhostVectorSource()) {
236 +                
237 +                oss << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
238 +                throw OOPSEException(oss.str());
239 +            }else{
240 +                int ghostIndex = bendStamp->getGhostVectorSource();
241 +                if (ghostIndex < getNAtoms()) {
242 +                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
243 +                      
244 +                      oss <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
245 +                      throw OOPSEException(oss.str());
246 +                    }
247 +                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
248 +                        
249 +                        oss <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
250 +                        throw OOPSEException(oss.str());
251 +                    }
252 +                }else {
253 +                    oss << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
254 +                    throw OOPSEException(oss.str());
255 +                }
256 +            }
257 +        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
258 +            oss <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
259 +            throw OOPSEException(oss.str());
260 +        }
261 +    }
262 +
263      for(int i = 0; i < getNBends(); ++i) {
264          BendStamp* bendStamp = getBendStamp(i);
265          std::vector<int> bendAtoms =  bendStamp->getMembers();
# Line 190 | Line 270 | void MoleculeStamp::validate() {
270              if (rigidbodyIndex >= 0) {
271                  ++rigidSet[rigidbodyIndex];
272                  if (rigidSet[rigidbodyIndex] > 1) {
273 <                    std::cout << "Error in Molecule " << getName() << ": ";
274 <                    //std::cout << "atoms of bend " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
273 >                    oss << "Error in Molecule " << getName() << ": bend" << containerToString(bendAtoms) << " belong to same rigidbody " << rigidbodyIndex << "\n";  
274 >                    throw OOPSEException(oss.str());
275                  }
276              }
277          }
278 <    }      
278 >    }
279 >    
280 >    
281 >    std::set<IntTuple3> allBends;
282 >    std::set<IntTuple3>::iterator iter;
283 >    for(int i = 0; i < getNBends(); ++i) {
284 >        BendStamp* bendStamp= getBendStamp(i);
285 >        std::vector<int> bend = bendStamp->getMembers();
286 >        if (bend.size() == 2) {
287 >        // in case we have two ghost bend. For example,
288 >        // bend {
289 >        // members (0, 1);
290 >        //   ghostVectorSource = 0;
291 >        // }
292 >        // and
293 >        // bend {
294 >        //   members (0, 1);
295 >        // ghostVectorSource = 0;
296 >        // }
297 >        // In order to distinguish them. we expand them to Tuple3.
298 >        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
299 >             int ghostIndex = bendStamp->getGhostVectorSource();
300 >             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
301 >             if (j != bend.end()) {
302 >                bend.insert(j, ghostIndex);
303 >             }
304 >        }
305 >        
306 >        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
307 >        //make sure bendTuple.first is always less than or equal to bendTuple.third
308 >        if (bendTuple.first > bendTuple.third) {
309 >            std::swap(bendTuple.first, bendTuple.third);
310 >        }
311 >        
312 >        iter = allBends.find(bendTuple);
313 >        if ( iter != allBends.end()) {
314 >            oss << "Error in Molecule " << getName() << ": " << "Bend" << containerToString(bend)<< " appears multiple times\n";
315 >            throw OOPSEException(oss.str());
316 >        } else {
317 >            allBends.insert(bendTuple);
318 >        }
319 >    }
320 >
321 >    for (int i = 0; i < getNBonds(); ++i) {
322 >        BondStamp* bondStamp = getBondStamp(i);
323 >        int a = bondStamp->getA();
324 >        int b = bondStamp->getB();
325 >
326 >        AtomStamp* atomA = getAtomStamp(a);
327 >        AtomStamp* atomB = getAtomStamp(b);
328 >
329 >        //find bend c--a--b
330 >        AtomStamp::AtomIter ai;
331 >        for(int c= atomA->getFirstBondedAtom(ai);c != -1;c = atomA->getNextBondedAtom(ai))
332 >        {
333 >            if(b == c)
334 >                continue;          
335 >            
336 >            IntTuple3 newBend(c, a, b);
337 >            if (newBend.first > newBend.third) {
338 >                std::swap(newBend.first, newBend.third);
339 >            }
340 >
341 >            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() {
371 >    std::ostringstream oss;
372      for(int i = 0; i < getNTorsions(); ++i) {
373          TorsionStamp* torsionStamp = getTorsionStamp(i);
374 +        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
375 +        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
376 +        std::vector<int>::iterator k =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::less<int>(), 0));
377 +
378 +        if (j != torsionAtoms.end() || k != torsionAtoms.end()) {
379 +            oss << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n";
380 +            throw OOPSEException(oss.str());
381 +        }
382 +        if (hasDuplicateElement(torsionAtoms)) {
383 +            oss << "Error in Molecule " << getName() << " : atoms of torsion" << containerToString(torsionAtoms) << " have duplicated indices\n";    
384 +            throw OOPSEException(oss.str());            
385 +        }        
386 +    }
387 +    
388 +    for(int i = 0; i < getNTorsions(); ++i) {
389 +        TorsionStamp* torsionStamp = getTorsionStamp(i);
390          std::vector<int> torsionAtoms =  torsionStamp->getMembers();
391          std::vector<int> rigidSet(getNRigidBodies(), 0);
392          std::vector<int>::iterator j;
# Line 206 | Line 395 | void MoleculeStamp::validate() {
395              if (rigidbodyIndex >= 0) {
396                  ++rigidSet[rigidbodyIndex];
397                  if (rigidSet[rigidbodyIndex] > 1) {
398 <                    std::cout << "Error in Molecule " << getName() << ": ";
399 <                    //std::cout << "atoms of torsion " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
398 >                    oss << "Error in Molecule " << getName() << ": torsion" << containerToString(torsionAtoms) << "is invalid\n";          
399 >                    throw OOPSEException(oss.str());
400                  }
401              }
402          }
403 <    }
403 >    }    
404  
405 +    std::set<IntTuple4> allTorsions;
406 +    std::set<IntTuple4>::iterator iter;
407 +     for(int i = 0; i < getNTorsions(); ++i) {
408 +         TorsionStamp* torsionStamp= getTorsionStamp(i);
409 +         std::vector<int> torsion = torsionStamp->getMembers();
410 +         if (torsion.size() == 3) {
411 +             int ghostIndex = torsionStamp->getGhostVectorSource();
412 +             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
413 +             if (j != torsion.end()) {
414 +                torsion.insert(j, ghostIndex);
415 +             }
416 +         }
417  
418 <    //fill in bond information into atom
419 <    fillBondInfo();
420 <    findBends();
421 <    findTorsions();
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 <    int nrigidAtoms = 0;
425 <    for (int i = 0; i < getNRigidBodies(); ++i) {
426 <        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
427 <        nrigidAtoms += rbStamp->getNMembers();
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      }
227 –    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
469  
470 < }
470 > }
471  
472 < void MoleculeStamp::fillBondInfo() {
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());
487 >        }
488 >        
489 >    }    
490   }
491  
492 < void MoleculeStamp::findBends() {
492 > void MoleculeStamp::checkCutoffGroups() {
493 >    for(int i = 0; i < getNCutoffGroups(); ++i) {
494 >        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
495 >        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
496 >        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
497 >        if (j != cutoffGroupAtoms.end()) {
498 >            std::ostringstream oss;
499 >            oss << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
500 >            throw OOPSEException(oss.str());
501 >        }
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";
511 +        throw OOPSEException(oss.str());
512 +    }
513 +    
514   }
515  
516 < void MoleculeStamp::findTorsions() {
516 > void MoleculeStamp::fillBondInfo() {
517  
518 +    for (int i = 0; i < getNBonds(); ++i) {
519 +        BondStamp* bondStamp = getBondStamp(i);
520 +        int a = bondStamp->getA();
521 +        int b = bondStamp->getB();
522 +        AtomStamp* atomA = getAtomStamp(a);
523 +        AtomStamp* atomB = getAtomStamp(b);
524 +        atomA->addBond(i);
525 +        atomA->addBondedAtom(b);
526 +        atomB->addBond(i);        
527 +        atomB->addBondedAtom(a);
528 +
529 +    }
530   }
531  
532 +
533 +
534   //Function Name: isBondInSameRigidBody
535   //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
536   bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
# Line 262 | Line 554 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
554   // Function Name: isAtomInRigidBody
555   //return false if atom does not belong to a rigid body, otherwise return true
556   bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
557 <  int whichRigidBody;
266 <  int consAtomIndex;
267 <
268 <  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
557 >  return atom2Rigidbody[atomIndex] >=0 ;
558    
559   }
560  
# Line 276 | Line 565 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex, i
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){
279 –  RigidBodyStamp* rbStamp;
280 –  int numRb;
281 –  int numAtom;
568  
569 +  
570 +
571    whichRigidBody = -1;
572    consAtomIndex = -1;
573  
574 <  numRb = this->getNRigidBodies();
575 <  
576 <  for(int i = 0 ; i < numRb; i++){
577 <    rbStamp = this->getRigidBodyStamp(i);
578 <    numAtom = rbStamp->getNMembers();
291 <    for(int j = 0; j < numAtom; j++)
574 >  if (atom2Rigidbody[atomIndex] >=0) {
575 >    whichRigidBody = atom2Rigidbody[atomIndex];
576 >    RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
577 >    int numAtom = rbStamp->getNMembers();
578 >    for(int j = 0; j < numAtom; j++) {
579        if (rbStamp->getMemberAt(j) == atomIndex){
293 –        whichRigidBody = i;
580          consAtomIndex = j;
581          return true;
582        }
583 +    }
584    }
585  
586    return false;

Diff Legend

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