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

Comparing trunk/OOPSE-2.0/src/brains/MoleculeCreator.cpp (file contents):
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2508 by gezelter, Mon Dec 12 19:32:50 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 <typeinfo>
52   #include <set>
53  
54   #include "brains/MoleculeCreator.hpp"
55   #include "primitives/GhostBend.hpp"
56 + #include "primitives/GhostTorsion.hpp"
57   #include "types/DirectionalAtomType.hpp"
58   #include "types/FixedBondType.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
61  
62   namespace oopse {
63 +  
64 +  Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
65 +                                            int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
66  
67 < Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
63 <    int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
64 <
65 <    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
67 >    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName());
68      
69      //create atoms
70      Atom* atom;
71      AtomStamp* currentAtomStamp;
72      int nAtom = molStamp->getNAtoms();
73      for (int i = 0; i < nAtom; ++i) {
74 <        currentAtomStamp = molStamp->getAtom(i);
75 <        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
76 <        mol->addAtom(atom);
74 >      currentAtomStamp = molStamp->getAtomStamp(i);
75 >      atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
76 >      mol->addAtom(atom);
77      }
78  
79      //create rigidbodies
# Line 80 | Line 82 | Molecule* MoleculeCreator::createMolecule(ForceField*
82      int nRigidbodies = molStamp->getNRigidBodies();
83  
84      for (int i = 0; i < nRigidbodies; ++i) {
85 <        currentRigidBodyStamp = molStamp->getRigidBody(i);
86 <        rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
87 <        mol->addRigidBody(rb);
85 >      currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
86 >      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
87 >      mol->addRigidBody(rb);
88      }
89  
90      //create bonds
# Line 91 | Line 93 | Molecule* MoleculeCreator::createMolecule(ForceField*
93      int nBonds = molStamp->getNBonds();
94  
95      for (int i = 0; i < nBonds; ++i) {
96 <        currentBondStamp = molStamp->getBond(i);
97 <        bond = createBond(ff, mol, currentBondStamp);
98 <        mol->addBond(bond);
96 >      currentBondStamp = molStamp->getBondStamp(i);
97 >      bond = createBond(ff, mol, currentBondStamp);
98 >      mol->addBond(bond);
99      }
100  
101      //create bends
# Line 101 | Line 103 | Molecule* MoleculeCreator::createMolecule(ForceField*
103      BendStamp* currentBendStamp;
104      int nBends = molStamp->getNBends();
105      for (int i = 0; i < nBends; ++i) {
106 <        currentBendStamp = molStamp->getBend(i);
107 <        bend = createBend(ff, mol, currentBendStamp);
108 <        mol->addBend(bend);
106 >      currentBendStamp = molStamp->getBendStamp(i);
107 >      bend = createBend(ff, mol, currentBendStamp);
108 >      mol->addBend(bend);
109      }
110  
111      //create torsions
# Line 111 | Line 113 | Molecule* MoleculeCreator::createMolecule(ForceField*
113      TorsionStamp* currentTorsionStamp;
114      int nTorsions = molStamp->getNTorsions();
115      for (int i = 0; i < nTorsions; ++i) {
116 <        currentTorsionStamp = molStamp->getTorsion(i);
117 <        torsion = createTorsion(ff, mol, currentTorsionStamp);
118 <        mol->addTorsion(torsion);
116 >      currentTorsionStamp = molStamp->getTorsionStamp(i);
117 >      torsion = createTorsion(ff, mol, currentTorsionStamp);
118 >      mol->addTorsion(torsion);
119      }
120  
121      //create cutoffGroups
# Line 121 | Line 123 | Molecule* MoleculeCreator::createMolecule(ForceField*
123      CutoffGroupStamp* currentCutoffGroupStamp;
124      int nCutoffGroups = molStamp->getNCutoffGroups();
125      for (int i = 0; i < nCutoffGroups; ++i) {
126 <        currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
127 <        cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
128 <        mol->addCutoffGroup(cutoffGroup);
126 >      currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
127 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
128 >      mol->addCutoffGroup(cutoffGroup);
129      }
130  
131      //every free atom is a cutoff group    
132 <    std::set<Atom*> allAtoms;
133 <     Molecule::AtomIterator ai;
132 >    std::vector<Atom*> freeAtoms;
133 >    std::vector<Atom*>::iterator ai;
134 >    std::vector<Atom*>::iterator fai;
135  
136      //add all atoms into allAtoms set
137 <    for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
138 <        allAtoms.insert(atom);
137 >    for(atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
138 >      freeAtoms.push_back(atom);
139      }
140  
141      Molecule::CutoffGroupIterator ci;
142      CutoffGroup* cg;
140 –    std::set<Atom*> cutoffAtoms;    
143      
142 –    //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 >           //erase the atoms belong to cutoff groups from freeAtoms vector
148 >           freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom), freeAtoms.end());
149 >      }
150  
151      }      
152      
151 –    //find all free atoms (which do not belong to cutoff groups)  
152 –    //performs the "difference" operation from set theory,  the output range contains a copy of every
153 –    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
154 –    //[cutoffAtoms.begin(), cutoffAtoms.end()).
155 –    std::vector<Atom*> freeAtoms;    
156 –    std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
157 –                            std::back_inserter(freeAtoms));
158 –
159 –    if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
160 –        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
161 –        sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
162 –
163 –        painCave.isFatal = 1;
164 –        simError();        
165 –    }
166 –
153      //loop over the free atoms and then create one cutoff group for every single free atom
154 <    std::vector<Atom*>::iterator fai;
169 <
154 >    
155      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
156 <        cutoffGroup = createCutoffGroup(mol, *fai);
157 <        mol->addCutoffGroup(cutoffGroup);
156 >      cutoffGroup = createCutoffGroup(mol, *fai);
157 >      mol->addCutoffGroup(cutoffGroup);
158      }
159      //create constraints
160      createConstraintPair(mol);
# Line 179 | Line 164 | Molecule* MoleculeCreator::createMolecule(ForceField*
164      mol->complete();
165  
166      return mol;
167 < }    
167 >  }    
168  
169  
170 < Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
171 <                                                                  LocalIndexManager* localIndexMan) {
170 >  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
171 >                                    LocalIndexManager* localIndexMan) {
172      AtomType * atomType;
173      Atom* atom;
174  
175      atomType =  ff->getAtomType(stamp->getType());
176  
177      if (atomType == NULL) {
178 <        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
179 <                   stamp->getType());
178 >      sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
179 >              stamp->getType().c_str());
180  
181 <        painCave.isFatal = 1;
182 <        simError();
181 >      painCave.isFatal = 1;
182 >      simError();
183      }
184      
185      //below code still have some kind of hard-coding smell
186      if (atomType->isDirectional()){
187      
188 <        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
188 >      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
189          
190 <        if (dAtomType == NULL) {
191 <            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
190 >      if (dAtomType == NULL) {
191 >        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
192  
193 <            painCave.isFatal = 1;
194 <            simError();
195 <        }
193 >        painCave.isFatal = 1;
194 >        simError();
195 >      }
196  
197 <        DirectionalAtom* dAtom;
198 <        dAtom = new DirectionalAtom(dAtomType);
199 <        atom = dAtom;    
197 >      DirectionalAtom* dAtom;
198 >      dAtom = new DirectionalAtom(dAtomType);
199 >      atom = dAtom;    
200      }
201      else{
202 <        atom = new Atom(atomType);
202 >      atom = new Atom(atomType);
203      }
204  
205      atom->setLocalIndex(localIndexMan->getNextAtomIndex());
206  
207      return atom;
208 < }
208 >  }
209  
210 < RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
211 <                                                                                    RigidBodyStamp* rbStamp,
212 <                                                                                    LocalIndexManager* localIndexMan) {
210 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
211 >                                              RigidBodyStamp* rbStamp,
212 >                                              LocalIndexManager* localIndexMan) {
213      Atom* atom;
214      int nAtoms;
215      Vector3d refCoor;
# Line 233 | Line 218 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
218      RigidBody* rb = new RigidBody();
219      nAtoms = rbStamp->getNMembers();    
220      for (int i = 0; i < nAtoms; ++i) {
221 <        //rbStamp->getMember(i) return the local index of current atom inside the molecule.
222 <        //It is not the same as local index of atom which is the index of atom at DataStorage class
223 <        atom = mol->getAtomAt(rbStamp->getMember(i));
224 <        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
225 <        rb->addAtom(atom, atomStamp);
221 >      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
222 >      //It is not the same as local index of atom which is the index of atom at DataStorage class
223 >      atom = mol->getAtomAt(rbStamp->getMemberAt(i));
224 >      atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));    
225 >      rb->addAtom(atom, atomStamp);
226      }
227  
228      //after all of the atoms are added, we need to calculate the reference coordinates
# Line 252 | Line 237 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
237      //The third part is the index of the rigidbody defined in meta-data file
238      //For example, Butane_RB_0 is a valid rigid body name of butane molecule
239      /**@todo replace itoa by lexi_cast */
240 <    rb->setType(mol->getType() + "_RB_" + toString(mol->getNRigidBodies()));
241 <    
240 >    std::string s = OOPSE_itoa(mol->getNRigidBodies(), 10);
241 >    rb->setType(mol->getType() + "_RB_" + s.c_str());
242 >
243      return rb;
244 < }    
244 >  }    
245  
246 < Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
246 >  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
247      BondType* bondType;
248      Atom* atomA;
249      Atom* atomB;
# Line 270 | Line 256 | Bond* MoleculeCreator::createBond(ForceField* ff, Mole
256      bondType = ff->getBondType(atomA->getType(), atomB->getType());
257  
258      if (bondType == NULL) {
259 <        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
260 <                   atomA->getType().c_str(),
261 <                   atomB->getType().c_str());
259 >      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
260 >              atomA->getType().c_str(),
261 >              atomB->getType().c_str());
262  
263 <        painCave.isFatal = 1;
264 <        simError();
263 >      painCave.isFatal = 1;
264 >      simError();
265      }
266      return new Bond(atomA, atomB, bondType);    
267 < }    
267 >  }    
268  
269 < Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
270 <    bool isGhostBend = false;
271 <    int ghostIndex;
269 >  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
270 >    Bend* bend = NULL;
271 >    std::vector<int> bendAtoms = stamp->getMembers();
272 >    if (bendAtoms.size() == 3) {
273 >      Atom* atomA = mol->getAtomAt(bendAtoms[0]);
274 >      Atom* atomB = mol->getAtomAt(bendAtoms[1]);
275 >      Atom* atomC = mol->getAtomAt(bendAtoms[2]);
276  
277 <    
288 <    //
289 <    if (stamp->haveExtras()){
290 <        LinkedAssign* extras = stamp->getExtras();
291 <        LinkedAssign* currentExtra = extras;
277 >      assert( atomA && atomB && atomC);
278  
279 <        while (currentExtra != NULL){
294 <            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
295 <                switch (currentExtra->getType()){
296 <                case 0:
297 <                    ghostIndex = currentExtra->getInt();
298 <                    isGhostBend = true;
299 <                    break;
279 >      BendType* bendType = ff->getBendType(atomA->getType().c_str(), atomB->getType().c_str(), atomC->getType().c_str());
280  
281 <                default:
282 <                sprintf(painCave.errMsg,
283 <                "SimSetup Error: ghostVectorSource must be an int.\n");
284 <                painCave.isFatal = 1;
285 <                simError();
306 <                }
307 <            } else{
308 <                sprintf(painCave.errMsg,
309 <                "SimSetup Error: unhandled bend assignment:\n");
310 <                painCave.isFatal = 1;
311 <                simError();
312 <            }
313 <            currentExtra = currentExtra->getNext();
314 <        }
315 <        
316 <    }
317 <
318 <    if (isGhostBend) {
319 <
320 <        int indexA = stamp->getA();
321 <        int indexB= stamp->getB();
322 <
323 <        assert(indexA != indexB);
324 <
325 <        int normalIndex;
326 <        if (indexA == ghostIndex) {
327 <            normalIndex = indexB;
328 <        } else if (indexB == ghostIndex) {
329 <            normalIndex = indexA;
330 <        }
331 <        
332 <        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
333 <        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
334 <        if (ghostAtom == NULL) {
335 <            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
336 <            painCave.isFatal = 1;
337 <            simError();
338 <        }
339 <                
340 <        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
341 <
342 <        if (bendType == NULL) {
343 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
344 <                       normalAtom->getType().c_str(),
345 <                       ghostAtom->getType().c_str(),
346 <                       "GHOST");
281 >      if (bendType == NULL) {
282 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
283 >                atomA->getType().c_str(),
284 >                atomB->getType().c_str(),
285 >                atomC->getType().c_str());
286  
287 <            painCave.isFatal = 1;
288 <            simError();
289 <        }
351 <        
352 <        return new GhostBend(normalAtom, ghostAtom, bendType);      
287 >        painCave.isFatal = 1;
288 >        simError();
289 >      }
290  
291 <    } else {
292 <            
293 <        Atom* atomA = mol->getAtomAt(stamp->getA());
294 <        Atom* atomB = mol->getAtomAt(stamp->getB());
295 <        Atom* atomC = mol->getAtomAt(stamp->getC());
291 >      bend = new Bend(atomA, atomB, atomC, bendType);
292 >    } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
293 >      int ghostIndex = stamp->getGhostVectorSource();
294 >      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
295 >      Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
296 >      DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
297 >      if (ghostAtom == NULL) {
298 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
299 >        painCave.isFatal = 1;
300 >        simError();
301 >      }
302 >                
303 >      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
304  
305 <        assert( atomA && atomB && atomC);
305 >      if (bendType == NULL) {
306 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
307 >                normalAtom->getType().c_str(),
308 >                ghostAtom->getType().c_str(),
309 >                "GHOST");
310 >
311 >        painCave.isFatal = 1;
312 >        simError();
313 >      }
314          
315 <        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
315 >      bend = new GhostBend(normalAtom, ghostAtom, bendType);      
316  
317 <        if (bendType == NULL) {
318 <            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
319 <                       atomA->getType().c_str(),
320 <                       atomB->getType().c_str(),
368 <                       atomC->getType().c_str());
317 >    }
318 >    
319 >    return bend;
320 >  }    
321  
322 <            painCave.isFatal = 1;
371 <            simError();
372 <        }
322 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
323  
324 <        return new Bend(atomA, atomB, atomC, bendType);      
324 >    Torsion* torsion = NULL;
325 >    std::vector<int> torsionAtoms = stamp->getMembers();
326 >    if (torsionAtoms.size() < 3) {
327 >        return torsion;
328      }
376 – }    
329  
330 < Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
331 <    TorsionType* torsionType;
332 <    Atom* atomA;
381 <    Atom* atomB;
382 <    Atom* atomC;
383 <    Atom* atomD;
330 >    Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
331 >    Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
332 >    Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
333  
334 <    atomA = mol->getAtomAt(stamp->getA());
335 <    atomB = mol->getAtomAt(stamp->getB());
387 <    atomC = mol->getAtomAt(stamp->getC());
388 <    atomD = mol->getAtomAt(stamp->getD());
334 >    if (torsionAtoms.size() == 4) {
335 >      Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
336  
337 <    assert(atomA && atomB && atomC && atomD);
338 <    
339 <    torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
340 <                                                       atomC->getType(), atomD->getType());
337 >      assert(atomA && atomB && atomC && atomD);
338 >        
339 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
340 >                                                    atomC->getType(), atomD->getType());
341  
342 <    if (torsionType == NULL) {
343 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
344 <                   atomA->getType().c_str(),
345 <                   atomB->getType().c_str(),
346 <                   atomC->getType().c_str(),
347 <                   atomD->getType().c_str());
342 >      if (torsionType == NULL) {
343 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
344 >                atomA->getType().c_str(),
345 >                atomB->getType().c_str(),
346 >                atomC->getType().c_str(),
347 >                atomD->getType().c_str());
348  
349 <        painCave.isFatal = 1;
350 <        simError();
349 >        painCave.isFatal = 1;
350 >        simError();
351 >      }
352 >        
353 >      torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
354      }
355 <    
406 <    return new Torsion(atomA, atomB, atomC, atomD, torsionType);      
407 < }    
355 >    else {
356  
357 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
357 >      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
358 >      if (dAtom == NULL) {
359 >        sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
360 >        painCave.isFatal = 1;
361 >        simError();
362 >      }        
363 >
364 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
365 >                                                    atomC->getType(), "GHOST");
366 >
367 >      if (torsionType == NULL) {
368 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
369 >                atomA->getType().c_str(),
370 >                atomB->getType().c_str(),
371 >                atomC->getType().c_str(),
372 >                "GHOST");
373 >
374 >        painCave.isFatal = 1;
375 >        simError();
376 >      }
377 >        
378 >      torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
379 >    }
380 >
381 >    return torsion;
382 >  }    
383 >
384 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
385      int nAtoms;
386      CutoffGroup* cg;
387      Atom* atom;
# Line 414 | Line 389 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
389      
390      nAtoms = stamp->getNMembers();
391      for (int i =0; i < nAtoms; ++i) {
392 <        atom = mol->getAtomAt(stamp->getMember(i));
393 <        assert(atom);
394 <        cg->addAtom(atom);
392 >      atom = mol->getAtomAt(stamp->getMemberAt(i));
393 >      assert(atom);
394 >      cg->addAtom(atom);
395      }
396  
397      return cg;
398 < }    
398 >  }    
399  
400 < CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
400 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
401      CutoffGroup* cg;
402      cg  = new CutoffGroup();
403      cg->addAtom(atom);
404      return cg;
405 < }
405 >  }
406  
407 < void MoleculeCreator::createConstraintPair(Molecule* mol) {
407 >  void MoleculeCreator::createConstraintPair(Molecule* mol) {
408  
409      //add bond constraints
410      Molecule::BondIterator bi;
411      Bond* bond;
412      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
413          
414 <        BondType* bt = bond->getBondType();
414 >      BondType* bt = bond->getBondType();
415  
416 <        //class Parent1 {};
417 <        //class Child1 : public Parent {};
418 <        //class Child2 : public Parent {};
419 <        //Child1* ch1 = new Child1();
420 <        //Child2* ch2 = dynamic_cast<Child2*>(ch1);
421 <        //the dynamic_cast is succeed in above line. A compiler bug?        
416 >      //class Parent1 {};
417 >      //class Child1 : public Parent {};
418 >      //class Child2 : public Parent {};
419 >      //Child1* ch1 = new Child1();
420 >      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
421 >      //the dynamic_cast is succeed in above line. A compiler bug?        
422  
423 <        if (typeid(FixedBondType) == typeid(*bt)) {
424 <            FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
423 >      if (typeid(FixedBondType) == typeid(*bt)) {
424 >        FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
425  
426 <            ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
427 <            ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
428 <            ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
429 <            mol->addConstraintPair(consPair);
430 <        }
426 >        ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
427 >        ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
428 >        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
429 >        mol->addConstraintPair(consPair);
430 >      }
431      }
432  
433      //rigidbody -- rigidbody constraint is not support yet
434 < }
434 >  }
435  
436 < void MoleculeCreator::createConstraintElem(Molecule* mol) {
436 >  void MoleculeCreator::createConstraintElem(Molecule* mol) {
437  
438      ConstraintPair* consPair;
439      Molecule::ConstraintPairIterator cpi;
440      std::set<StuntDouble*> sdSet;
441      for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
442  
443 <        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
444 <        if (sdSet.find(sdA) == sdSet.end()){
445 <            sdSet.insert(sdA);
446 <            mol->addConstraintElem(new ConstraintElem(sdA));
447 <        }
443 >      StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
444 >      if (sdSet.find(sdA) == sdSet.end()){
445 >        sdSet.insert(sdA);
446 >        mol->addConstraintElem(new ConstraintElem(sdA));
447 >      }
448  
449 <        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
450 <        if (sdSet.find(sdB) == sdSet.end()){
451 <            sdSet.insert(sdB);
452 <            mol->addConstraintElem(new ConstraintElem(sdB));
453 <        }
449 >      StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
450 >      if (sdSet.find(sdB) == sdSet.end()){
451 >        sdSet.insert(sdB);
452 >        mol->addConstraintElem(new ConstraintElem(sdB));
453 >      }
454          
455      }
456  
457 < }
457 >  }
458      
459   }

Diff Legend

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