ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/MoleculeCreator.cpp
(Generate patch)

Comparing trunk/src/brains/MoleculeCreator.cpp (file contents):
Revision 273 by tim, Tue Jan 25 17:45:23 2005 UTC vs.
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 39 | Line 39
39   * such damages.
40   */
41    
42 < /**
43 <  * @file MoleculeCreator.cpp
44 <  * @author tlin
45 <  * @date 11/04/2004
46 <  * @time 13:44am
47 <  * @version 1.0
48 <  */
42 > /**
43 > * @file MoleculeCreator.cpp
44 > * @author tlin
45 > * @date 11/04/2004
46 > * @time 13:44am
47 > * @version 1.0
48 > */
49  
50   #include <cassert>
51   #include <set>
# Line 60 | Line 60 | namespace oopse {
60  
61   namespace oopse {
62  
63 < Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
64 <    int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
63 >  Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
64 >                                            int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
65  
66      Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
67      
# Line 70 | Line 70 | Molecule* MoleculeCreator::createMolecule(ForceField*
70      AtomStamp* currentAtomStamp;
71      int nAtom = molStamp->getNAtoms();
72      for (int i = 0; i < nAtom; ++i) {
73 <        currentAtomStamp = molStamp->getAtom(i);
74 <        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
75 <        mol->addAtom(atom);
73 >      currentAtomStamp = molStamp->getAtom(i);
74 >      atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
75 >      mol->addAtom(atom);
76      }
77  
78      //create rigidbodies
# Line 81 | Line 81 | Molecule* MoleculeCreator::createMolecule(ForceField*
81      int nRigidbodies = molStamp->getNRigidBodies();
82  
83      for (int i = 0; i < nRigidbodies; ++i) {
84 <        currentRigidBodyStamp = molStamp->getRigidBody(i);
85 <        rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
86 <        mol->addRigidBody(rb);
84 >      currentRigidBodyStamp = molStamp->getRigidBody(i);
85 >      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
86 >      mol->addRigidBody(rb);
87      }
88  
89      //create bonds
# Line 92 | Line 92 | Molecule* MoleculeCreator::createMolecule(ForceField*
92      int nBonds = molStamp->getNBonds();
93  
94      for (int i = 0; i < nBonds; ++i) {
95 <        currentBondStamp = molStamp->getBond(i);
96 <        bond = createBond(ff, mol, currentBondStamp);
97 <        mol->addBond(bond);
95 >      currentBondStamp = molStamp->getBond(i);
96 >      bond = createBond(ff, mol, currentBondStamp);
97 >      mol->addBond(bond);
98      }
99  
100      //create bends
# Line 102 | Line 102 | Molecule* MoleculeCreator::createMolecule(ForceField*
102      BendStamp* currentBendStamp;
103      int nBends = molStamp->getNBends();
104      for (int i = 0; i < nBends; ++i) {
105 <        currentBendStamp = molStamp->getBend(i);
106 <        bend = createBend(ff, mol, currentBendStamp);
107 <        mol->addBend(bend);
105 >      currentBendStamp = molStamp->getBend(i);
106 >      bend = createBend(ff, mol, currentBendStamp);
107 >      mol->addBend(bend);
108      }
109  
110      //create torsions
# Line 112 | Line 112 | Molecule* MoleculeCreator::createMolecule(ForceField*
112      TorsionStamp* currentTorsionStamp;
113      int nTorsions = molStamp->getNTorsions();
114      for (int i = 0; i < nTorsions; ++i) {
115 <        currentTorsionStamp = molStamp->getTorsion(i);
116 <        torsion = createTorsion(ff, mol, currentTorsionStamp);
117 <        mol->addTorsion(torsion);
115 >      currentTorsionStamp = molStamp->getTorsion(i);
116 >      torsion = createTorsion(ff, mol, currentTorsionStamp);
117 >      mol->addTorsion(torsion);
118      }
119  
120      //create cutoffGroups
# Line 122 | Line 122 | Molecule* MoleculeCreator::createMolecule(ForceField*
122      CutoffGroupStamp* currentCutoffGroupStamp;
123      int nCutoffGroups = molStamp->getNCutoffGroups();
124      for (int i = 0; i < nCutoffGroups; ++i) {
125 <        currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
126 <        cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
127 <        mol->addCutoffGroup(cutoffGroup);
125 >      currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
126 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
127 >      mol->addCutoffGroup(cutoffGroup);
128      }
129  
130      //every free atom is a cutoff group    
131      std::set<Atom*> allAtoms;
132 <     Molecule::AtomIterator ai;
132 >    Molecule::AtomIterator ai;
133  
134      //add all atoms into allAtoms set
135      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
136 <        allAtoms.insert(atom);
136 >      allAtoms.insert(atom);
137      }
138  
139      Molecule::CutoffGroupIterator ci;
# Line 143 | Line 143 | Molecule* MoleculeCreator::createMolecule(ForceField*
143      //add all of the atoms belong to cutoff groups into cutoffAtoms set
144      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
145  
146 <        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
147 <            cutoffAtoms.insert(atom);
148 <        }
146 >      for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
147 >        cutoffAtoms.insert(atom);
148 >      }
149  
150      }      
151      
# Line 155 | Line 155 | Molecule* MoleculeCreator::createMolecule(ForceField*
155      //[cutoffAtoms.begin(), cutoffAtoms.end()).
156      std::vector<Atom*> freeAtoms;    
157      std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
158 <                            std::back_inserter(freeAtoms));
158 >                        std::back_inserter(freeAtoms));
159  
160      if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
161 <        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
162 <        sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
161 >      //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
162 >      sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
163  
164 <        painCave.isFatal = 1;
165 <        simError();        
164 >      painCave.isFatal = 1;
165 >      simError();        
166      }
167  
168      //loop over the free atoms and then create one cutoff group for every single free atom
169      std::vector<Atom*>::iterator fai;
170  
171      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
172 <        cutoffGroup = createCutoffGroup(mol, *fai);
173 <        mol->addCutoffGroup(cutoffGroup);
172 >      cutoffGroup = createCutoffGroup(mol, *fai);
173 >      mol->addCutoffGroup(cutoffGroup);
174      }
175      //create constraints
176      createConstraintPair(mol);
# Line 180 | Line 180 | Molecule* MoleculeCreator::createMolecule(ForceField*
180      mol->complete();
181  
182      return mol;
183 < }    
183 >  }    
184  
185  
186 < Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
187 <                                                                  LocalIndexManager* localIndexMan) {
186 >  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
187 >                                    LocalIndexManager* localIndexMan) {
188      AtomType * atomType;
189      Atom* atom;
190  
191      atomType =  ff->getAtomType(stamp->getType());
192  
193      if (atomType == NULL) {
194 <        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
195 <                   stamp->getType());
194 >      sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
195 >              stamp->getType());
196  
197 <        painCave.isFatal = 1;
198 <        simError();
197 >      painCave.isFatal = 1;
198 >      simError();
199      }
200      
201      //below code still have some kind of hard-coding smell
202      if (atomType->isDirectional()){
203      
204 <        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
204 >      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
205          
206 <        if (dAtomType == NULL) {
207 <            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
206 >      if (dAtomType == NULL) {
207 >        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
208  
209 <            painCave.isFatal = 1;
210 <            simError();
211 <        }
209 >        painCave.isFatal = 1;
210 >        simError();
211 >      }
212  
213 <        DirectionalAtom* dAtom;
214 <        dAtom = new DirectionalAtom(dAtomType);
215 <        atom = dAtom;    
213 >      DirectionalAtom* dAtom;
214 >      dAtom = new DirectionalAtom(dAtomType);
215 >      atom = dAtom;    
216      }
217      else{
218 <        atom = new Atom(atomType);
218 >      atom = new Atom(atomType);
219      }
220  
221      atom->setLocalIndex(localIndexMan->getNextAtomIndex());
222  
223      return atom;
224 < }
224 >  }
225  
226 < RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
227 <                                                                                    RigidBodyStamp* rbStamp,
228 <                                                                                    LocalIndexManager* localIndexMan) {
226 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
227 >                                              RigidBodyStamp* rbStamp,
228 >                                              LocalIndexManager* localIndexMan) {
229      Atom* atom;
230      int nAtoms;
231      Vector3d refCoor;
# Line 234 | Line 234 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
234      RigidBody* rb = new RigidBody();
235      nAtoms = rbStamp->getNMembers();    
236      for (int i = 0; i < nAtoms; ++i) {
237 <        //rbStamp->getMember(i) return the local index of current atom inside the molecule.
238 <        //It is not the same as local index of atom which is the index of atom at DataStorage class
239 <        atom = mol->getAtomAt(rbStamp->getMember(i));
240 <        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
241 <        rb->addAtom(atom, atomStamp);
237 >      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
238 >      //It is not the same as local index of atom which is the index of atom at DataStorage class
239 >      atom = mol->getAtomAt(rbStamp->getMember(i));
240 >      atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
241 >      rb->addAtom(atom, atomStamp);
242      }
243  
244      //after all of the atoms are added, we need to calculate the reference coordinates
# Line 253 | Line 253 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
253      //The third part is the index of the rigidbody defined in meta-data file
254      //For example, Butane_RB_0 is a valid rigid body name of butane molecule
255      /**@todo replace itoa by lexi_cast */
256 <    rb->setType(mol->getType() + "_RB_" + toString(mol->getNRigidBodies()));
257 <    
256 >    std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10);
257 >    rb->setType(mol->getType() + "_RB_" + s.c_str());
258 >
259      return rb;
260 < }    
260 >  }    
261  
262 < Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
262 >  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
263      BondType* bondType;
264      Atom* atomA;
265      Atom* atomB;
# Line 271 | Line 272 | Bond* MoleculeCreator::createBond(ForceField* ff, Mole
272      bondType = ff->getBondType(atomA->getType(), atomB->getType());
273  
274      if (bondType == NULL) {
275 <        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
276 <                   atomA->getType().c_str(),
277 <                   atomB->getType().c_str());
275 >      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
276 >              atomA->getType().c_str(),
277 >              atomB->getType().c_str());
278  
279 <        painCave.isFatal = 1;
280 <        simError();
279 >      painCave.isFatal = 1;
280 >      simError();
281      }
282      return new Bond(atomA, atomB, bondType);    
283 < }    
283 >  }    
284  
285 < Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
285 >  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
286      bool isGhostBend = false;
287      int ghostIndex;
288  
289      
290      //
291      if (stamp->haveExtras()){
292 <        LinkedAssign* extras = stamp->getExtras();
293 <        LinkedAssign* currentExtra = extras;
292 >      LinkedAssign* extras = stamp->getExtras();
293 >      LinkedAssign* currentExtra = extras;
294  
295 <        while (currentExtra != NULL){
296 <            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
297 <                switch (currentExtra->getType()){
298 <                case 0:
299 <                    ghostIndex = currentExtra->getInt();
300 <                    isGhostBend = true;
301 <                    break;
295 >      while (currentExtra != NULL){
296 >        if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
297 >          switch (currentExtra->getType()){
298 >          case 0:
299 >            ghostIndex = currentExtra->getInt();
300 >            isGhostBend = true;
301 >            break;
302  
303 <                default:
304 <                sprintf(painCave.errMsg,
305 <                "SimSetup Error: ghostVectorSource must be an int.\n");
306 <                painCave.isFatal = 1;
307 <                simError();
308 <                }
309 <            } else{
310 <                sprintf(painCave.errMsg,
311 <                "SimSetup Error: unhandled bend assignment:\n");
312 <                painCave.isFatal = 1;
313 <                simError();
314 <            }
315 <            currentExtra = currentExtra->getNext();
316 <        }
303 >          default:
304 >            sprintf(painCave.errMsg,
305 >                    "SimSetup Error: ghostVectorSource must be an int.\n");
306 >            painCave.isFatal = 1;
307 >            simError();
308 >          }
309 >        } else{
310 >          sprintf(painCave.errMsg,
311 >                  "SimSetup Error: unhandled bend assignment:\n");
312 >          painCave.isFatal = 1;
313 >          simError();
314 >        }
315 >        currentExtra = currentExtra->getNext();
316 >      }
317          
318      }
319  
320      if (isGhostBend) {
321  
322 <        int indexA = stamp->getA();
323 <        int indexB= stamp->getB();
322 >      int indexA = stamp->getA();
323 >      int indexB= stamp->getB();
324  
325 <        assert(indexA != indexB);
325 >      assert(indexA != indexB);
326  
327 <        int normalIndex;
328 <        if (indexA == ghostIndex) {
329 <            normalIndex = indexB;
330 <        } else if (indexB == ghostIndex) {
331 <            normalIndex = indexA;
332 <        }
327 >      int normalIndex;
328 >      if (indexA == ghostIndex) {
329 >        normalIndex = indexB;
330 >      } else if (indexB == ghostIndex) {
331 >        normalIndex = indexA;
332 >      }
333          
334 <        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
335 <        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
336 <        if (ghostAtom == NULL) {
337 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
338 <            painCave.isFatal = 1;
339 <            simError();
340 <        }
334 >      Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
335 >      DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
336 >      if (ghostAtom == NULL) {
337 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
338 >        painCave.isFatal = 1;
339 >        simError();
340 >      }
341                  
342 <        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
342 >      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
343  
344 <        if (bendType == NULL) {
345 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
346 <                       normalAtom->getType().c_str(),
347 <                       ghostAtom->getType().c_str(),
348 <                       "GHOST");
344 >      if (bendType == NULL) {
345 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
346 >                normalAtom->getType().c_str(),
347 >                ghostAtom->getType().c_str(),
348 >                "GHOST");
349  
350 <            painCave.isFatal = 1;
351 <            simError();
352 <        }
350 >        painCave.isFatal = 1;
351 >        simError();
352 >      }
353          
354 <        return new GhostBend(normalAtom, ghostAtom, bendType);      
354 >      return new GhostBend(normalAtom, ghostAtom, bendType);      
355  
356      } else {
357              
358 <        Atom* atomA = mol->getAtomAt(stamp->getA());
359 <        Atom* atomB = mol->getAtomAt(stamp->getB());
360 <        Atom* atomC = mol->getAtomAt(stamp->getC());
358 >      Atom* atomA = mol->getAtomAt(stamp->getA());
359 >      Atom* atomB = mol->getAtomAt(stamp->getB());
360 >      Atom* atomC = mol->getAtomAt(stamp->getC());
361  
362 <        assert( atomA && atomB && atomC);
362 >      assert( atomA && atomB && atomC);
363          
364 <        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
364 >      BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
365  
366 <        if (bendType == NULL) {
367 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
368 <                       atomA->getType().c_str(),
369 <                       atomB->getType().c_str(),
370 <                       atomC->getType().c_str());
366 >      if (bendType == NULL) {
367 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
368 >                atomA->getType().c_str(),
369 >                atomB->getType().c_str(),
370 >                atomC->getType().c_str());
371  
372 <            painCave.isFatal = 1;
373 <            simError();
374 <        }
372 >        painCave.isFatal = 1;
373 >        simError();
374 >      }
375  
376 <        return new Bend(atomA, atomB, atomC, bendType);      
376 >      return new Bend(atomA, atomB, atomC, bendType);      
377      }
378 < }    
378 >  }    
379  
380 < Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
380 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
381  
382      Atom* atomA = mol->getAtomAt(stamp->getA());
383      Atom* atomB = mol->getAtomAt(stamp->getB());
# Line 384 | Line 385 | Torsion* MoleculeCreator::createTorsion(ForceField* ff
385      Torsion* torsion;
386  
387      if (stamp->getD() != -1) {
388 <        Atom* atomD = mol->getAtomAt(stamp->getD());
388 >      Atom* atomD = mol->getAtomAt(stamp->getD());
389  
390 <        assert(atomA && atomB && atomC && atomD);
390 >      assert(atomA && atomB && atomC && atomD);
391          
392 <        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
393 <                                                           atomC->getType(), atomD->getType());
392 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
393 >                                                    atomC->getType(), atomD->getType());
394  
395 <        if (torsionType == NULL) {
396 <            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
397 <                       atomA->getType().c_str(),
398 <                       atomB->getType().c_str(),
399 <                       atomC->getType().c_str(),
400 <                       atomD->getType().c_str());
395 >      if (torsionType == NULL) {
396 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
397 >                atomA->getType().c_str(),
398 >                atomB->getType().c_str(),
399 >                atomC->getType().c_str(),
400 >                atomD->getType().c_str());
401  
402 <            painCave.isFatal = 1;
403 <            simError();
404 <        }
402 >        painCave.isFatal = 1;
403 >        simError();
404 >      }
405          
406 <        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
406 >      torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
407      }
408      else {
409  
410 <        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
411 <        if (dAtom == NULL) {
412 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
413 <            painCave.isFatal = 1;
414 <            simError();
415 <        }        
410 >      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
411 >      if (dAtom == NULL) {
412 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
413 >        painCave.isFatal = 1;
414 >        simError();
415 >      }        
416  
417 <        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
418 <                                                           atomC->getType(), "GHOST");
417 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
418 >                                                    atomC->getType(), "GHOST");
419  
420 <        if (torsionType == NULL) {
421 <            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
422 <                       atomA->getType().c_str(),
423 <                       atomB->getType().c_str(),
424 <                       atomC->getType().c_str(),
425 <                       "GHOST");
420 >      if (torsionType == NULL) {
421 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
422 >                atomA->getType().c_str(),
423 >                atomB->getType().c_str(),
424 >                atomC->getType().c_str(),
425 >                "GHOST");
426  
427 <            painCave.isFatal = 1;
428 <            simError();
429 <        }
427 >        painCave.isFatal = 1;
428 >        simError();
429 >      }
430          
431 <        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
431 >      torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
432      }
433  
434      return torsion;
435 < }    
435 >  }    
436  
437 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
437 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
438      int nAtoms;
439      CutoffGroup* cg;
440      Atom* atom;
# Line 441 | Line 442 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
442      
443      nAtoms = stamp->getNMembers();
444      for (int i =0; i < nAtoms; ++i) {
445 <        atom = mol->getAtomAt(stamp->getMember(i));
446 <        assert(atom);
447 <        cg->addAtom(atom);
445 >      atom = mol->getAtomAt(stamp->getMember(i));
446 >      assert(atom);
447 >      cg->addAtom(atom);
448      }
449  
450      return cg;
451 < }    
451 >  }    
452  
453 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
453 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
454      CutoffGroup* cg;
455      cg  = new CutoffGroup();
456      cg->addAtom(atom);
457      return cg;
458 < }
458 >  }
459  
460 < void MoleculeCreator::createConstraintPair(Molecule* mol) {
460 >  void MoleculeCreator::createConstraintPair(Molecule* mol) {
461  
462      //add bond constraints
463      Molecule::BondIterator bi;
464      Bond* bond;
465      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
466          
467 <        BondType* bt = bond->getBondType();
467 >      BondType* bt = bond->getBondType();
468  
469 <        //class Parent1 {};
470 <        //class Child1 : public Parent {};
471 <        //class Child2 : public Parent {};
472 <        //Child1* ch1 = new Child1();
473 <        //Child2* ch2 = dynamic_cast<Child2*>(ch1);
474 <        //the dynamic_cast is succeed in above line. A compiler bug?        
469 >      //class Parent1 {};
470 >      //class Child1 : public Parent {};
471 >      //class Child2 : public Parent {};
472 >      //Child1* ch1 = new Child1();
473 >      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
474 >      //the dynamic_cast is succeed in above line. A compiler bug?        
475  
476 <        if (typeid(FixedBondType) == typeid(*bt)) {
477 <            FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
476 >      if (typeid(FixedBondType) == typeid(*bt)) {
477 >        FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
478  
479 <            ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
480 <            ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
481 <            ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
482 <            mol->addConstraintPair(consPair);
483 <        }
479 >        ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
480 >        ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
481 >        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
482 >        mol->addConstraintPair(consPair);
483 >      }
484      }
485  
486      //rigidbody -- rigidbody constraint is not support yet
487 < }
487 >  }
488  
489 < void MoleculeCreator::createConstraintElem(Molecule* mol) {
489 >  void MoleculeCreator::createConstraintElem(Molecule* mol) {
490  
491      ConstraintPair* consPair;
492      Molecule::ConstraintPairIterator cpi;
493      std::set<StuntDouble*> sdSet;
494      for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
495  
496 <        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
497 <        if (sdSet.find(sdA) == sdSet.end()){
498 <            sdSet.insert(sdA);
499 <            mol->addConstraintElem(new ConstraintElem(sdA));
500 <        }
496 >      StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
497 >      if (sdSet.find(sdA) == sdSet.end()){
498 >        sdSet.insert(sdA);
499 >        mol->addConstraintElem(new ConstraintElem(sdA));
500 >      }
501  
502 <        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
503 <        if (sdSet.find(sdB) == sdSet.end()){
504 <            sdSet.insert(sdB);
505 <            mol->addConstraintElem(new ConstraintElem(sdB));
506 <        }
502 >      StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
503 >      if (sdSet.find(sdB) == sdSet.end()){
504 >        sdSet.insert(sdB);
505 >        mol->addConstraintElem(new ConstraintElem(sdB));
506 >      }
507          
508      }
509  
510 < }
510 >  }
511      
512   }

Diff Legend

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