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 292 by tim, Fri Feb 4 22:44:15 2005 UTC vs.
Revision 647 by tim, Wed Oct 5 19:34:01 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 59 | Line 59 | namespace oopse {
59   #include "utils/StringUtils.hpp"
60  
61   namespace oopse {
62 +  
63 +  Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
64 +                                            int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
65  
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      
68      //create atoms
# 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;
131 >    std::vector<Atom*> freeAtoms;
132 >    std::vector<Atom*>::iterator ai;
133 >    std::vector<Atom*>::iterator fai;
134  
135      //add all atoms into allAtoms set
136 <    for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
137 <        allAtoms.insert(atom);
136 >    for(atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
137 >      freeAtoms.push_back(atom);
138      }
139  
140      Molecule::CutoffGroupIterator ci;
141      CutoffGroup* cg;
141 –    std::set<Atom*> cutoffAtoms;    
142      
143 –    //add all of the atoms belong to cutoff groups into cutoffAtoms set
143      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
144  
145 <        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
146 <            cutoffAtoms.insert(atom);
147 <        }
145 >      for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
146 >           //erase the atoms belong to cutoff groups from freeAtoms vector
147 >           freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom), freeAtoms.end());
148 >      }
149  
150      }      
151      
152 –    //find all free atoms (which do not belong to cutoff groups)  
153 –    //performs the "difference" operation from set theory,  the output range contains a copy of every
154 –    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
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));
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");
163 –
164 –        painCave.isFatal = 1;
165 –        simError();        
166 –    }
167 –
152      //loop over the free atoms and then create one cutoff group for every single free atom
153 <    std::vector<Atom*>::iterator fai;
170 <
153 >    
154      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
155 <        cutoffGroup = createCutoffGroup(mol, *fai);
156 <        mol->addCutoffGroup(cutoffGroup);
155 >      cutoffGroup = createCutoffGroup(mol, *fai);
156 >      mol->addCutoffGroup(cutoffGroup);
157      }
158      //create constraints
159      createConstraintPair(mol);
# Line 180 | Line 163 | Molecule* MoleculeCreator::createMolecule(ForceField*
163      mol->complete();
164  
165      return mol;
166 < }    
166 >  }    
167  
168  
169 < Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
170 <                                                                  LocalIndexManager* localIndexMan) {
169 >  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
170 >                                    LocalIndexManager* localIndexMan) {
171      AtomType * atomType;
172      Atom* atom;
173  
174      atomType =  ff->getAtomType(stamp->getType());
175  
176      if (atomType == NULL) {
177 <        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
178 <                   stamp->getType());
177 >      sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
178 >              stamp->getType());
179  
180 <        painCave.isFatal = 1;
181 <        simError();
180 >      painCave.isFatal = 1;
181 >      simError();
182      }
183      
184      //below code still have some kind of hard-coding smell
185      if (atomType->isDirectional()){
186      
187 <        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
187 >      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
188          
189 <        if (dAtomType == NULL) {
190 <            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
189 >      if (dAtomType == NULL) {
190 >        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
191  
192 <            painCave.isFatal = 1;
193 <            simError();
194 <        }
192 >        painCave.isFatal = 1;
193 >        simError();
194 >      }
195  
196 <        DirectionalAtom* dAtom;
197 <        dAtom = new DirectionalAtom(dAtomType);
198 <        atom = dAtom;    
196 >      DirectionalAtom* dAtom;
197 >      dAtom = new DirectionalAtom(dAtomType);
198 >      atom = dAtom;    
199      }
200      else{
201 <        atom = new Atom(atomType);
201 >      atom = new Atom(atomType);
202      }
203  
204      atom->setLocalIndex(localIndexMan->getNextAtomIndex());
205  
206      return atom;
207 < }
207 >  }
208  
209 < RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
210 <                                                                                    RigidBodyStamp* rbStamp,
211 <                                                                                    LocalIndexManager* localIndexMan) {
209 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
210 >                                              RigidBodyStamp* rbStamp,
211 >                                              LocalIndexManager* localIndexMan) {
212      Atom* atom;
213      int nAtoms;
214      Vector3d refCoor;
# Line 234 | Line 217 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
217      RigidBody* rb = new RigidBody();
218      nAtoms = rbStamp->getNMembers();    
219      for (int i = 0; i < nAtoms; ++i) {
220 <        //rbStamp->getMember(i) return the local index of current atom inside the molecule.
221 <        //It is not the same as local index of atom which is the index of atom at DataStorage class
222 <        atom = mol->getAtomAt(rbStamp->getMember(i));
223 <        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
224 <        rb->addAtom(atom, atomStamp);
220 >      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
221 >      //It is not the same as local index of atom which is the index of atom at DataStorage class
222 >      atom = mol->getAtomAt(rbStamp->getMember(i));
223 >      atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
224 >      rb->addAtom(atom, atomStamp);
225      }
226  
227      //after all of the atoms are added, we need to calculate the reference coordinates
# Line 253 | Line 236 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
236      //The third part is the index of the rigidbody defined in meta-data file
237      //For example, Butane_RB_0 is a valid rigid body name of butane molecule
238      /**@todo replace itoa by lexi_cast */
239 <    char buffer [33];  
240 <    sprintf(buffer, "%d", mol->getNRigidBodies());
258 <    //rb->setType(mol->getType() + "_RB_" + toString<int>(mol->getNRigidBodies()));
259 <    rb->setType(mol->getType() + "_RB_" + buffer);
239 >    std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10);
240 >    rb->setType(mol->getType() + "_RB_" + s.c_str());
241  
242      return rb;
243 < }    
243 >  }    
244  
245 < Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
245 >  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
246      BondType* bondType;
247      Atom* atomA;
248      Atom* atomB;
# Line 274 | Line 255 | Bond* MoleculeCreator::createBond(ForceField* ff, Mole
255      bondType = ff->getBondType(atomA->getType(), atomB->getType());
256  
257      if (bondType == NULL) {
258 <        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
259 <                   atomA->getType().c_str(),
260 <                   atomB->getType().c_str());
258 >      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
259 >              atomA->getType().c_str(),
260 >              atomB->getType().c_str());
261  
262 <        painCave.isFatal = 1;
263 <        simError();
262 >      painCave.isFatal = 1;
263 >      simError();
264      }
265      return new Bond(atomA, atomB, bondType);    
266 < }    
266 >  }    
267  
268 < Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
268 >  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
269      bool isGhostBend = false;
270      int ghostIndex;
271  
272      
273      //
274      if (stamp->haveExtras()){
275 <        LinkedAssign* extras = stamp->getExtras();
276 <        LinkedAssign* currentExtra = extras;
275 >      LinkedAssign* extras = stamp->getExtras();
276 >      LinkedAssign* currentExtra = extras;
277  
278 <        while (currentExtra != NULL){
279 <            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
280 <                switch (currentExtra->getType()){
281 <                case 0:
282 <                    ghostIndex = currentExtra->getInt();
283 <                    isGhostBend = true;
284 <                    break;
278 >      while (currentExtra != NULL){
279 >        if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
280 >          switch (currentExtra->getType()){
281 >          case 0:
282 >            ghostIndex = currentExtra->getInt();
283 >            isGhostBend = true;
284 >            break;
285  
286 <                default:
287 <                sprintf(painCave.errMsg,
288 <                "SimSetup Error: ghostVectorSource must be an int.\n");
289 <                painCave.isFatal = 1;
290 <                simError();
291 <                }
292 <            } else{
293 <                sprintf(painCave.errMsg,
294 <                "SimSetup Error: unhandled bend assignment:\n");
295 <                painCave.isFatal = 1;
296 <                simError();
297 <            }
298 <            currentExtra = currentExtra->getNext();
299 <        }
286 >          default:
287 >            sprintf(painCave.errMsg,
288 >                    "SimSetup Error: ghostVectorSource must be an int.\n");
289 >            painCave.isFatal = 1;
290 >            simError();
291 >          }
292 >        } else{
293 >          sprintf(painCave.errMsg,
294 >                  "SimSetup Error: unhandled bend assignment:\n");
295 >          painCave.isFatal = 1;
296 >          simError();
297 >        }
298 >        currentExtra = currentExtra->getNext();
299 >      }
300          
301      }
302  
303      if (isGhostBend) {
304  
305 <        int indexA = stamp->getA();
306 <        int indexB= stamp->getB();
305 >      int indexA = stamp->getA();
306 >      int indexB= stamp->getB();
307  
308 <        assert(indexA != indexB);
308 >      assert(indexA != indexB);
309  
310 <        int normalIndex;
311 <        if (indexA == ghostIndex) {
312 <            normalIndex = indexB;
313 <        } else if (indexB == ghostIndex) {
314 <            normalIndex = indexA;
315 <        }
310 >      int normalIndex;
311 >      if (indexA == ghostIndex) {
312 >        normalIndex = indexB;
313 >      } else if (indexB == ghostIndex) {
314 >        normalIndex = indexA;
315 >      }
316          
317 <        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
318 <        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
319 <        if (ghostAtom == NULL) {
320 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
321 <            painCave.isFatal = 1;
322 <            simError();
323 <        }
324 <                
325 <        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
317 >      Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
318 >      DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
319 >      if (ghostAtom == NULL) {
320 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
321 >        painCave.isFatal = 1;
322 >        simError();
323 >      }
324 >                
325 >      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
326  
327 <        if (bendType == NULL) {
328 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
329 <                       normalAtom->getType().c_str(),
330 <                       ghostAtom->getType().c_str(),
331 <                       "GHOST");
327 >      if (bendType == NULL) {
328 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
329 >                normalAtom->getType().c_str(),
330 >                ghostAtom->getType().c_str(),
331 >                "GHOST");
332  
333 <            painCave.isFatal = 1;
334 <            simError();
335 <        }
333 >        painCave.isFatal = 1;
334 >        simError();
335 >      }
336          
337 <        return new GhostBend(normalAtom, ghostAtom, bendType);      
337 >      return new GhostBend(normalAtom, ghostAtom, bendType);      
338  
339      } else {
340              
341 <        Atom* atomA = mol->getAtomAt(stamp->getA());
342 <        Atom* atomB = mol->getAtomAt(stamp->getB());
343 <        Atom* atomC = mol->getAtomAt(stamp->getC());
341 >      Atom* atomA = mol->getAtomAt(stamp->getA());
342 >      Atom* atomB = mol->getAtomAt(stamp->getB());
343 >      Atom* atomC = mol->getAtomAt(stamp->getC());
344  
345 <        assert( atomA && atomB && atomC);
345 >      assert( atomA && atomB && atomC);
346          
347 <        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
347 >      BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
348  
349 <        if (bendType == NULL) {
350 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
351 <                       atomA->getType().c_str(),
352 <                       atomB->getType().c_str(),
353 <                       atomC->getType().c_str());
349 >      if (bendType == NULL) {
350 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
351 >                atomA->getType().c_str(),
352 >                atomB->getType().c_str(),
353 >                atomC->getType().c_str());
354  
355 <            painCave.isFatal = 1;
356 <            simError();
357 <        }
355 >        painCave.isFatal = 1;
356 >        simError();
357 >      }
358  
359 <        return new Bend(atomA, atomB, atomC, bendType);      
359 >      return new Bend(atomA, atomB, atomC, bendType);      
360      }
361 < }    
361 >  }    
362  
363 < Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
363 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
364  
365      Atom* atomA = mol->getAtomAt(stamp->getA());
366      Atom* atomB = mol->getAtomAt(stamp->getB());
# Line 387 | Line 368 | Torsion* MoleculeCreator::createTorsion(ForceField* ff
368      Torsion* torsion;
369  
370      if (stamp->getD() != -1) {
371 <        Atom* atomD = mol->getAtomAt(stamp->getD());
371 >      Atom* atomD = mol->getAtomAt(stamp->getD());
372  
373 <        assert(atomA && atomB && atomC && atomD);
373 >      assert(atomA && atomB && atomC && atomD);
374          
375 <        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
376 <                                                           atomC->getType(), atomD->getType());
375 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
376 >                                                    atomC->getType(), atomD->getType());
377  
378 <        if (torsionType == NULL) {
379 <            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
380 <                       atomA->getType().c_str(),
381 <                       atomB->getType().c_str(),
382 <                       atomC->getType().c_str(),
383 <                       atomD->getType().c_str());
378 >      if (torsionType == NULL) {
379 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
380 >                atomA->getType().c_str(),
381 >                atomB->getType().c_str(),
382 >                atomC->getType().c_str(),
383 >                atomD->getType().c_str());
384  
385 <            painCave.isFatal = 1;
386 <            simError();
387 <        }
385 >        painCave.isFatal = 1;
386 >        simError();
387 >      }
388          
389 <        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
389 >      torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
390      }
391      else {
392  
393 <        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
394 <        if (dAtom == NULL) {
395 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
396 <            painCave.isFatal = 1;
397 <            simError();
398 <        }        
393 >      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
394 >      if (dAtom == NULL) {
395 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
396 >        painCave.isFatal = 1;
397 >        simError();
398 >      }        
399  
400 <        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
401 <                                                           atomC->getType(), "GHOST");
400 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
401 >                                                    atomC->getType(), "GHOST");
402  
403 <        if (torsionType == NULL) {
404 <            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
405 <                       atomA->getType().c_str(),
406 <                       atomB->getType().c_str(),
407 <                       atomC->getType().c_str(),
408 <                       "GHOST");
403 >      if (torsionType == NULL) {
404 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
405 >                atomA->getType().c_str(),
406 >                atomB->getType().c_str(),
407 >                atomC->getType().c_str(),
408 >                "GHOST");
409  
410 <            painCave.isFatal = 1;
411 <            simError();
412 <        }
410 >        painCave.isFatal = 1;
411 >        simError();
412 >      }
413          
414 <        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
414 >      torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
415      }
416  
417      return torsion;
418 < }    
418 >  }    
419  
420 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
420 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
421      int nAtoms;
422      CutoffGroup* cg;
423      Atom* atom;
# Line 444 | Line 425 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
425      
426      nAtoms = stamp->getNMembers();
427      for (int i =0; i < nAtoms; ++i) {
428 <        atom = mol->getAtomAt(stamp->getMember(i));
429 <        assert(atom);
430 <        cg->addAtom(atom);
428 >      atom = mol->getAtomAt(stamp->getMember(i));
429 >      assert(atom);
430 >      cg->addAtom(atom);
431      }
432  
433      return cg;
434 < }    
434 >  }    
435  
436 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
436 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
437      CutoffGroup* cg;
438      cg  = new CutoffGroup();
439      cg->addAtom(atom);
440      return cg;
441 < }
441 >  }
442  
443 < void MoleculeCreator::createConstraintPair(Molecule* mol) {
443 >  void MoleculeCreator::createConstraintPair(Molecule* mol) {
444  
445      //add bond constraints
446      Molecule::BondIterator bi;
447      Bond* bond;
448      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
449          
450 <        BondType* bt = bond->getBondType();
450 >      BondType* bt = bond->getBondType();
451  
452 <        //class Parent1 {};
453 <        //class Child1 : public Parent {};
454 <        //class Child2 : public Parent {};
455 <        //Child1* ch1 = new Child1();
456 <        //Child2* ch2 = dynamic_cast<Child2*>(ch1);
457 <        //the dynamic_cast is succeed in above line. A compiler bug?        
452 >      //class Parent1 {};
453 >      //class Child1 : public Parent {};
454 >      //class Child2 : public Parent {};
455 >      //Child1* ch1 = new Child1();
456 >      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
457 >      //the dynamic_cast is succeed in above line. A compiler bug?        
458  
459 <        if (typeid(FixedBondType) == typeid(*bt)) {
460 <            FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
459 >      if (typeid(FixedBondType) == typeid(*bt)) {
460 >        FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
461  
462 <            ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
463 <            ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
464 <            ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
465 <            mol->addConstraintPair(consPair);
466 <        }
462 >        ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
463 >        ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
464 >        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
465 >        mol->addConstraintPair(consPair);
466 >      }
467      }
468  
469      //rigidbody -- rigidbody constraint is not support yet
470 < }
470 >  }
471  
472 < void MoleculeCreator::createConstraintElem(Molecule* mol) {
472 >  void MoleculeCreator::createConstraintElem(Molecule* mol) {
473  
474      ConstraintPair* consPair;
475      Molecule::ConstraintPairIterator cpi;
476      std::set<StuntDouble*> sdSet;
477      for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
478  
479 <        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
480 <        if (sdSet.find(sdA) == sdSet.end()){
481 <            sdSet.insert(sdA);
482 <            mol->addConstraintElem(new ConstraintElem(sdA));
483 <        }
479 >      StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
480 >      if (sdSet.find(sdA) == sdSet.end()){
481 >        sdSet.insert(sdA);
482 >        mol->addConstraintElem(new ConstraintElem(sdA));
483 >      }
484  
485 <        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
486 <        if (sdSet.find(sdB) == sdSet.end()){
487 <            sdSet.insert(sdB);
488 <            mol->addConstraintElem(new ConstraintElem(sdB));
489 <        }
485 >      StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
486 >      if (sdSet.find(sdB) == sdSet.end()){
487 >        sdSet.insert(sdB);
488 >        mol->addConstraintElem(new ConstraintElem(sdB));
489 >      }
490          
491      }
492  
493 < }
493 >  }
494      
495   }

Diff Legend

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