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 1719 by tim, Fri Nov 5 23:38:27 2004 UTC vs.
Revision 1805 by tim, Tue Nov 30 20:50:47 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) {
44 <    Molecule* mol = new Molecule(stampId, 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
49      Atom* atom;
# 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 +     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 +    Molecule::CutoffGroupIterator ci;
120 +    CutoffGroup* cg;
121 +    std::set<Atom*> cutoffAtoms;    
122 +    
123 +    //add all of the atoms belong to cutoff groups into cutoffAtoms set
124 +    for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
125 +
126 +        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
127 +            cutoffAtoms.insert(atom);
128 +        }
129 +
130 +    }      
131 +    
132 +    //find all free atoms (which do not belong to cutoff groups)  
133 +    //performs the "difference" operation from set theory,  the output range contains a copy of every
134 +    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
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));
139 +
140 +    if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
141 +        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
142 +        sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
143 +
144 +        painCave.isFatal = 1;
145 +        simError();        
146 +    }
147 +
148 +    //loop over the free atoms and then create one cutoff group for every single free atom
149 +    std::vector<Atom*>::iterator fai;
150 +
151 +    for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
152 +        createCutoffGroup(mol, *fai);
153 +        mol->addCutoffGroup(cutoffGroup);
154 +    }
155      //create constraints
156  
157      //the construction of this molecule is finished
# Line 117 | Line 168 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
168  
169      atomType =  ff->getAtomType(stamp->getType());
170  
171 +    if (atomType == NULL) {
172 +        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
173 +                   stamp->getType());
174 +
175 +        painCave.isFatal = 1;
176 +        simError();
177 +    }
178 +    
179      //below code still have some kind of hard-coding smell
180      if (stamp->haveOrientation()){
181          DirectionalAtom* dAtom;
182          double phi;
183          double theta;
184          double psi;
185 +
186 +        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
187 +        if (dAtomType == NULL) {
188 +            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
189 +
190 +            painCave.isFatal = 1;
191 +            simError();
192 +        }
193          
194 <        dAtom = new DirectionalAtom(atomType);
194 >        dAtom = new DirectionalAtom(dAtomType);
195  
196          // Directional Atoms have standard unit vectors which are oriented
197          // in space using the three Euler angles.  We assume the standard
# Line 142 | Line 209 | Atom* MoleculeCreator::createAtom(ForceField* ff, Mole
209      }
210  
211      atom->setLocalIndex(localIndexMan->getNextAtomIndex());
212 +
213 +    return atom;
214   }
215  
216   RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
# Line 152 | Line 221 | RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
221      Vector3d refCoor;
222      AtomStamp* atomStamp;
223      
155 –    nAtoms = stamp->getNMembers();
156 –
224      RigidBody* rb = new RigidBody();
225 <    
225 >    nAtoms = rbStamp->getNMembers();    
226      for (int i = 0; i < nAtoms; ++i) {
227          //rbStamp->getMember(i) return the local index of current atom inside the molecule.
228          //It is not the same as local index of atom which is the index of atom at DataStorage class
229          atom = mol->getAtomAt(rbStamp->getMember(i));
230 <        atomStamp= molStamp->molStamp->getAtom(rbStamp->getMember(i));    
230 >        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
231          rb->addAtom(atom, atomStamp);
232      }
233 <    
233 >
234 >    //after all of the atoms are added, we need to calculate the reference coordinates
235      rb->calcRefCoords();
236 +
237 +    //set the local index of this rigid body, global index will be set later
238      rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
239 <    rb->setType(mol->getType() + "_" + itoa(mol.getNRigidBodies()));
239 >
240 >    //the rule for naming rigidbody MoleculeName_RB_Integer
241 >    //The first part is the name of the molecule
242 >    //The second part is alway fixed as "RB"
243 >    //The third part is the index of the rigidbody defined in meta-data file
244 >    //For example, Butane_RB_0 is a valid rigid body name of butane molecule
245 >    /**@todo replace itoa by lexi_cast */
246 >    rb->setType(mol->getType() + "_RB_" + toString(mol->getNRigidBodies()));
247      
248      return rb;
249   }    
# Line 181 | Line 258 | Bond* MoleculeCreator::createBond(ForceField* ff, Mole
258  
259      assert( atomA && atomB);
260      
261 <    bondType = ff->getBondType(atomA, atomB);
261 >    bondType = ff->getBondType(atomA->getType(), atomB->getType());
262  
263 +    if (bondType == NULL) {
264 +        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
265 +                   atomA->getType().c_str(),
266 +                   atomB->getType().c_str());
267 +
268 +        painCave.isFatal = 1;
269 +        simError();
270 +    }
271      return new Bond(atomA, atomB, bondType);    
272   }    
273  
274   Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
275 <    BendType* benType;
276 <    Atom* atomA;
277 <    Atom* atomB;
278 <    Atom* atomC;
275 >    bool isGhostBend = false;
276 >    int ghostIndex;
277 >
278 >    
279 >    //
280 >    if (stamp->haveExtras()){
281 >        LinkedAssign* extras = stamp->getExtras();
282 >        LinkedAssign* currentExtra = extras;
283  
284 <    //need to consider the ghost bend
285 <    atomA = mol->getAtomAt(stamp->getA());
286 <    atomB = mol->getAtomAt(stamp->getB());
287 <    atomC = mol->getAtomAt(stamp->getC());
284 >        while (currentExtra != NULL){
285 >            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
286 >                switch (currentExtra->getType()){
287 >                case 0:
288 >                    ghostIndex = currentExtra->getInt();
289 >                    isGhostBend = true;
290 >                    break;
291  
292 <    assert( atomA && atomB && atomC);
293 <    
294 <    benType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
292 >                default:
293 >                sprintf(painCave.errMsg,
294 >                "SimSetup Error: ghostVectorSource must be an int.\n");
295 >                painCave.isFatal = 1;
296 >                simError();
297 >                }
298 >            } else{
299 >                sprintf(painCave.errMsg,
300 >                "SimSetup Error: unhandled bend assignment:\n");
301 >                painCave.isFatal = 1;
302 >                simError();
303 >            }
304 >            currentExtra = currentExtra->getNext();
305 >        }
306 >        
307 >    }
308  
309 <    return new Bond(atomA, atomB, benType);      
309 >    if (isGhostBend) {
310  
311 +        int indexA = stamp->getA();
312 +        int indexB= stamp->getB();
313 +
314 +        assert(indexA != indexB);
315 +
316 +        int normalIndex;
317 +        if (indexA == ghostIndex) {
318 +            normalIndex = indexB;
319 +        } else if (indexB == ghostIndex) {
320 +            normalIndex = indexA;
321 +        }
322 +        
323 +        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
324 +        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
325 +        if (ghostAtom == NULL) {
326 +            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
327 +            painCave.isFatal = 1;
328 +            simError();
329 +        }
330 +                
331 +        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
332 +
333 +        if (bendType == NULL) {
334 +            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
335 +                       normalAtom->getType().c_str(),
336 +                       ghostAtom->getType().c_str(),
337 +                       "GHOST");
338 +
339 +            painCave.isFatal = 1;
340 +            simError();
341 +        }
342 +        
343 +        return new GhostBend(normalAtom, ghostAtom, bendType);      
344 +
345 +    } else {
346 +            
347 +        Atom* atomA = mol->getAtomAt(stamp->getA());
348 +        Atom* atomB = mol->getAtomAt(stamp->getB());
349 +        Atom* atomC = mol->getAtomAt(stamp->getC());
350 +
351 +        assert( atomA && atomB && atomC);
352 +        
353 +        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
354 +
355 +        if (bendType == NULL) {
356 +            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
357 +                       atomA->getType().c_str(),
358 +                       atomB->getType().c_str(),
359 +                       atomC->getType().c_str());
360 +
361 +            painCave.isFatal = 1;
362 +            simError();
363 +        }
364 +
365 +        return new Bend(atomA, atomB, atomC, bendType);      
366 +    }
367   }    
368  
369   Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
# Line 219 | Line 380 | Torsion* MoleculeCreator::createTorsion(ForceField* ff
380  
381      assert(atomA && atomB && atomC && atomD);
382      
383 <    torsionType = ff->getTosionType(atomA->getType(), atomB->getType(),
383 >    torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
384                                                         atomC->getType(), atomD->getType());
385  
386 <    return new Torsion(atomA, atomB, torsionType);      
386 >    if (torsionType == NULL) {
387 >        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
388 >                   atomA->getType().c_str(),
389 >                   atomB->getType().c_str(),
390 >                   atomC->getType().c_str(),
391 >                   atomD->getType().c_str());
392 >
393 >        painCave.isFatal = 1;
394 >        simError();
395 >    }
396 >    
397 >    return new Torsion(atomA, atomB, atomC, atomD, torsionType);      
398   }    
399  
400   CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
# Line 241 | Line 413 | CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
413      return cg;
414   }    
415  
416 + CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
417 +    CutoffGroup* cg;
418 +    cg  = new CutoffGroup();
419 +    cg->addAtom(atom);
420 +    return cg;
421 + }
422   //Constraint* MoleculeCreator::createConstraint() {
423  
424   //}

Diff Legend

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