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

Diff Legend

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