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 2469 by tim, Fri Dec 2 15:38:03 2005 UTC vs.
Revision 2483 by tim, Mon Dec 5 18:23:30 2005 UTC

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

Diff Legend

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