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

Comparing branches/new_design/OOPSE-4/src/brains/MoleculeCreator.cpp (file contents):
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC vs.
Revision 1813 by tim, Wed Dec 1 17:38:32 2004 UTC

# Line 31 | Line 31
31    * @version 1.0
32    */
33  
34 – #include "brains/MoleculeCreator.hpp"
34   #include <cassert>
35  
36 + #include "brains/MoleculeCreator.hpp"
37 + #include "primitives/GhostBend.hpp"
38 + #include "utils/simError.h"
39 + #include "utils/StringUtils.hpp"
40 +
41   namespace oopse {
42  
43 < Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp, int stampId, int globalIndex) {
43 > Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
44 >    int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
45 >
46      Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
47      
48      //create atoms
# Line 45 | Line 51 | Molecule* MoleculeCreator::createMolecule(ForceField*
51      int nAtom = molStamp->getNAtoms();
52      for (int i = 0; i < nAtom; ++i) {
53          currentAtomStamp = molStamp->getAtom(i);
54 <        atom = createAtom(ff, currentAtomStamp);
54 >        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
55          mol->addAtom(atom);
56      }
57  
# Line 56 | Line 62 | Molecule* MoleculeCreator::createMolecule(ForceField*
62  
63      for (int i = 0; i < nRigidbodies; ++i) {
64          currentRigidBodyStamp = molStamp->getRigidBody(i);
65 <        rb = createRigidBody(ff, mol, currentRigidBodyStamp);
65 >        rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
66          mol->addRigidBody(rb);
67      }
68  
# Line 66 | Line 72 | Molecule* MoleculeCreator::createMolecule(ForceField*
72      int nBonds = molStamp->getNBends();
73  
74      for (int i = 0; i < nBonds; ++i) {
75 <        currentBondStamp = molStamp->getBend(i);
75 >        currentBondStamp = molStamp->getBond(i);
76          bond = createBond(ff, mol, currentBondStamp);
77          mol->addBond(bond);
78      }
# Line 97 | Line 103 | Molecule* MoleculeCreator::createMolecule(ForceField*
103      int nCutoffGroups = molStamp->getNCutoffGroups();
104      for (int i = 0; i < nCutoffGroups; ++i) {
105          currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
106 <        cutoffGroup = createCutoffGroup(ff, mol, currentCutoffGroupStamp);
106 >        cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
107          mol->addCutoffGroup(cutoffGroup);
108      }
109  
110      //every free atom is a cutoff group    
111      std::set<Atom*> allAtoms;
112 <    typename  Molecule::AtomIterator ai;
112 >     Molecule::AtomIterator ai;
113  
114      //add all atoms into allAtoms set
115      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
116          allAtoms.insert(atom);
117      }
118  
119 <    typename Molecule::CutoffGroupIterator ci;
119 >    Molecule::CutoffGroupIterator ci;
120      CutoffGroup* cg;
121      std::set<Atom*> cutoffAtoms;    
122      
# Line 129 | Line 135 | Molecule* MoleculeCreator::createMolecule(ForceField*
135      //[cutoffAtoms.begin(), cutoffAtoms.end()).
136      std::vector<Atom*> freeAtoms;    
137      std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
138 <                            std::back_inserter(freeAtoms.end()));
138 >                            std::back_inserter(freeAtoms));
139  
140      if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
141          //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
# Line 162 | Line 168 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
168  
169      atomType =  ff->getAtomType(stamp->getType());
170  
171 <    if (bondType == NULL) {
171 >    if (atomType == NULL) {
172          sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
173                     stamp->getType());
174  
# Line 171 | Line 177 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
177      }
178      
179      //below code still have some kind of hard-coding smell
180 <    if (stamp->haveOrientation()){
181 <        DirectionalAtom* dAtom;
182 <        double phi;
177 <        double theta;
178 <        double psi;
180 >    if (atomType->isDirectional()){
181 >    
182 >        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
183          
184 <        dAtom = new DirectionalAtom(atomType);
184 >        if (dAtomType == NULL) {
185 >            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
186  
187 <        // Directional Atoms have standard unit vectors which are oriented
188 <        // in space using the three Euler angles.  We assume the standard
189 <        // unit vector was originally along the z axis below.
187 >            painCave.isFatal = 1;
188 >            simError();
189 >        }
190  
191 <        phi = stamp->getEulerPhi() * M_PI / 180.0;
192 <        theta = stamp->getEulerTheta() * M_PI / 180.0;
188 <        psi = stamp->getEulerPsi()* M_PI / 180.0;
189 <
190 <        dAtom->setUnitFrameFromEuler(phi, theta, psi);
191 >        DirectionalAtom* dAtom;
192 >        dAtom = new DirectionalAtom(dAtomType);
193          atom = dAtom;    
194      }
195      else{
# Line 195 | Line 197 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
197      }
198  
199      atom->setLocalIndex(localIndexMan->getNextAtomIndex());
200 +
201 +    return atom;
202   }
203  
204   RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
# Line 205 | Line 209 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
209      Vector3d refCoor;
210      AtomStamp* atomStamp;
211      
208 –    nAtoms = molStamp->getNMembers();
209 –
212      RigidBody* rb = new RigidBody();
213 <    
213 >    nAtoms = rbStamp->getNMembers();    
214      for (int i = 0; i < nAtoms; ++i) {
215          //rbStamp->getMember(i) return the local index of current atom inside the molecule.
216          //It is not the same as local index of atom which is the index of atom at DataStorage class
217          atom = mol->getAtomAt(rbStamp->getMember(i));
218 <        atomStamp= molStamp->molStamp->getAtom(rbStamp->getMember(i));    
218 >        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
219          rb->addAtom(atom, atomStamp);
220      }
221  
# Line 228 | Line 230 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
230      //The second part is alway fixed as "RB"
231      //The third part is the index of the rigidbody defined in meta-data file
232      //For example, Butane_RB_0 is a valid rigid body name of butane molecule
233 <    rb->setType(mol->getType() + "_RB_" + itoa(mol.getNRigidBodies()));
233 >    /**@todo replace itoa by lexi_cast */
234 >    rb->setType(mol->getType() + "_RB_" + toString(mol->getNRigidBodies()));
235      
236      return rb;
237   }    
# Line 243 | Line 246 | Bond* MoleculeCreator::createBond(ForceField* ff, Mole
246  
247      assert( atomA && atomB);
248      
249 <    bondType = ff->getBondType(atomA, atomB);
249 >    bondType = ff->getBondType(atomA->getType(), atomB->getType());
250  
251      if (bondType == NULL) {
252          sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
# Line 263 | Line 266 | Bend* MoleculeCreator::createBend(ForceField* ff, Mole
266      
267      //
268      if (stamp->haveExtras()){
269 <        LinkedAssign* extras = currentBend->getExtras();
269 >        LinkedAssign* extras = stamp->getExtras();
270          LinkedAssign* currentExtra = extras;
271  
272          while (currentExtra != NULL){
# Line 305 | Line 308 | Bend* MoleculeCreator::createBend(ForceField* ff, Mole
308              normalIndex = indexA;
309          }
310          
311 <        Atom* normalAtom = mol->getAtomAt(normalIndex) ;
312 <        Atom* ghostAtom = mol->getAtomAt(ghostIndex);
313 <        
311 >        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
312 >        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
313 >        if (ghostAtom == NULL) {
314 >            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
315 >            painCave.isFatal = 1;
316 >            simError();
317 >        }
318 >                
319          BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
320  
321          if (bendType == NULL) {
# Line 319 | Line 327 | Bend* MoleculeCreator::createBend(ForceField* ff, Mole
327              painCave.isFatal = 1;
328              simError();
329          }
330 <
330 >        
331          return new GhostBend(normalAtom, ghostAtom, bendType);      
332  
333      } else {
# Line 360 | Line 368 | Torsion* MoleculeCreator::createTorsion(ForceField* ff
368  
369      assert(atomA && atomB && atomC && atomD);
370      
371 <    torsionType = ff->getTosionType(atomA->getType(), atomB->getType(),
371 >    torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
372                                                         atomC->getType(), atomD->getType());
373  
374      if (torsionType == NULL) {
# Line 374 | Line 382 | Torsion* MoleculeCreator::createTorsion(ForceField* ff
382          simError();
383      }
384      
385 <    return new Torsion(atomA, atomB, torsionType);      
385 >    return new Torsion(atomA, atomB, atomC, atomD, torsionType);      
386   }    
387  
388   CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
# Line 396 | Line 404 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
404   CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
405      CutoffGroup* cg;
406      cg  = new CutoffGroup();
407 <    cg->addAtom();
407 >    cg->addAtom(atom);
408      return cg;
409   }
410   //Constraint* MoleculeCreator::createConstraint() {

Diff Legend

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