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 1734 by tim, Fri Nov 12 07:05:43 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 162 | Line 169 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
169  
170      atomType =  ff->getAtomType(stamp->getType());
171  
172 <    if (bondType == NULL) {
172 >    if (atomType == NULL) {
173          sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
174                     stamp->getType());
175  
# 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;
183 <        double phi;
177 <        double theta;
178 <        double psi;
181 >    if (atomType->isDirectional()){
182 >    
183 >        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
184          
185 <        dAtom = new DirectionalAtom(atomType);
185 >        if (dAtomType == NULL) {
186 >            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
187  
188 <        // Directional Atoms have standard unit vectors which are oriented
189 <        // in space using the three Euler angles.  We assume the standard
190 <        // unit vector was originally along the z axis below.
188 >            painCave.isFatal = 1;
189 >            simError();
190 >        }
191  
192 <        phi = stamp->getEulerPhi() * M_PI / 180.0;
193 <        theta = stamp->getEulerTheta() * M_PI / 180.0;
188 <        psi = stamp->getEulerPsi()* M_PI / 180.0;
189 <
190 <        dAtom->setUnitFrameFromEuler(phi, theta, psi);
192 >        DirectionalAtom* dAtom;
193 >        dAtom = new DirectionalAtom(dAtomType);
194          atom = dAtom;    
195      }
196      else{
# Line 195 | 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 205 | Line 210 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
210      Vector3d refCoor;
211      AtomStamp* atomStamp;
212      
208 –    nAtoms = molStamp->getNMembers();
209 –
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 228 | 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 243 | 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 257 | Line 261 | Bend* MoleculeCreator::createBend(ForceField* ff, Mole
261   }    
262  
263   Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
264 <    BendType* bendType;
265 <    Atom* atomA;
262 <    Atom* atomB;
263 <    Atom* atomC;
264 >    bool isGhostBend = false;
265 >    int ghostIndex;
266  
265 –    //need to consider the ghost bend
266 –    atomA = mol->getAtomAt(stamp->getA());
267 –    atomB = mol->getAtomAt(stamp->getB());
268 –    atomC = mol->getAtomAt(stamp->getC());
269 –
270 –    assert( atomA && atomB && atomC);
267      
268 <    bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
268 >    //
269 >    if (stamp->haveExtras()){
270 >        LinkedAssign* extras = stamp->getExtras();
271 >        LinkedAssign* currentExtra = extras;
272  
273 <    if (bendType == NULL) {
274 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
275 <                   atomA->getType().c_str(),
276 <                   atomB->getType().c_str(),
277 <                   atomC->getType().c_str());
278 <
279 <        painCave.isFatal = 1;
280 <        simError();
273 >        while (currentExtra != NULL){
274 >            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
275 >                switch (currentExtra->getType()){
276 >                case 0:
277 >                    ghostIndex = currentExtra->getInt();
278 >                    isGhostBend = true;
279 >                    break;
280 >
281 >                default:
282 >                sprintf(painCave.errMsg,
283 >                "SimSetup Error: ghostVectorSource must be an int.\n");
284 >                painCave.isFatal = 1;
285 >                simError();
286 >                }
287 >            } else{
288 >                sprintf(painCave.errMsg,
289 >                "SimSetup Error: unhandled bend assignment:\n");
290 >                painCave.isFatal = 1;
291 >                simError();
292 >            }
293 >            currentExtra = currentExtra->getNext();
294 >        }
295 >        
296      }
297  
298 <    return new Bond(atomA, atomB, bendType);      
298 >    if (isGhostBend) {
299  
300 +        int indexA = stamp->getA();
301 +        int indexB= stamp->getB();
302 +
303 +        assert(indexA != indexB);
304 +
305 +        int normalIndex;
306 +        if (indexA == ghostIndex) {
307 +            normalIndex = indexB;
308 +        } else if (indexB == ghostIndex) {
309 +            normalIndex = indexA;
310 +        }
311 +        
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) {
323 +            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
324 +                       normalAtom->getType().c_str(),
325 +                       ghostAtom->getType().c_str(),
326 +                       "GHOST");
327 +
328 +            painCave.isFatal = 1;
329 +            simError();
330 +        }
331 +        
332 +        return new GhostBend(normalAtom, ghostAtom, bendType);      
333 +
334 +    } else {
335 +            
336 +        Atom* atomA = mol->getAtomAt(stamp->getA());
337 +        Atom* atomB = mol->getAtomAt(stamp->getB());
338 +        Atom* atomC = mol->getAtomAt(stamp->getC());
339 +
340 +        assert( atomA && atomB && atomC);
341 +        
342 +        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
343 +
344 +        if (bendType == NULL) {
345 +            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
346 +                       atomA->getType().c_str(),
347 +                       atomB->getType().c_str(),
348 +                       atomC->getType().c_str());
349 +
350 +            painCave.isFatal = 1;
351 +            simError();
352 +        }
353 +
354 +        return new Bend(atomA, atomB, atomC, bendType);      
355 +    }
356   }    
357  
358   Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
# Line 299 | 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 313 | 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 335 | 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