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

Comparing trunk/src/brains/MoleculeCreator.cpp (file contents):
Revision 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 809 by gezelter, Mon Dec 12 19:32:50 2005 UTC

# Line 48 | Line 48
48   */
49  
50   #include <cassert>
51 + #include <typeinfo>
52   #include <set>
53  
54   #include "brains/MoleculeCreator.hpp"
# Line 59 | Line 60 | namespace oopse {
60   #include "utils/StringUtils.hpp"
61  
62   namespace oopse {
63 <
63 >  
64    Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
65                                              int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
66  
67 <    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);
74 >      currentAtomStamp = molStamp->getAtomStamp(i);
75        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
76        mol->addAtom(atom);
77      }
# Line 81 | Line 82 | namespace oopse {
82      int nRigidbodies = molStamp->getNRigidBodies();
83  
84      for (int i = 0; i < nRigidbodies; ++i) {
85 <      currentRigidBodyStamp = molStamp->getRigidBody(i);
85 >      currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
86        rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
87        mol->addRigidBody(rb);
88      }
# Line 92 | Line 93 | namespace oopse {
93      int nBonds = molStamp->getNBonds();
94  
95      for (int i = 0; i < nBonds; ++i) {
96 <      currentBondStamp = molStamp->getBond(i);
96 >      currentBondStamp = molStamp->getBondStamp(i);
97        bond = createBond(ff, mol, currentBondStamp);
98        mol->addBond(bond);
99      }
# Line 102 | Line 103 | namespace oopse {
103      BendStamp* currentBendStamp;
104      int nBends = molStamp->getNBends();
105      for (int i = 0; i < nBends; ++i) {
106 <      currentBendStamp = molStamp->getBend(i);
106 >      currentBendStamp = molStamp->getBendStamp(i);
107        bend = createBend(ff, mol, currentBendStamp);
108        mol->addBend(bend);
109      }
# Line 112 | Line 113 | namespace oopse {
113      TorsionStamp* currentTorsionStamp;
114      int nTorsions = molStamp->getNTorsions();
115      for (int i = 0; i < nTorsions; ++i) {
116 <      currentTorsionStamp = molStamp->getTorsion(i);
116 >      currentTorsionStamp = molStamp->getTorsionStamp(i);
117        torsion = createTorsion(ff, mol, currentTorsionStamp);
118        mol->addTorsion(torsion);
119      }
# Line 122 | Line 123 | namespace oopse {
123      CutoffGroupStamp* currentCutoffGroupStamp;
124      int nCutoffGroups = molStamp->getNCutoffGroups();
125      for (int i = 0; i < nCutoffGroups; ++i) {
126 <      currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
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;
141 –    std::set<Atom*> cutoffAtoms;    
143      
143 –    //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);
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      
152 –    //find all free atoms (which do not belong to cutoff groups)  
153 –    //performs the "difference" operation from set theory,  the output range contains a copy of every
154 –    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in
155 –    //[cutoffAtoms.begin(), cutoffAtoms.end()).
156 –    std::vector<Atom*> freeAtoms;    
157 –    std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
158 –                        std::back_inserter(freeAtoms));
159 –
160 –    if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
161 –      //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
162 –      sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
163 –
164 –      painCave.isFatal = 1;
165 –      simError();        
166 –    }
167 –
153      //loop over the free atoms and then create one cutoff group for every single free atom
154 <    std::vector<Atom*>::iterator fai;
170 <
154 >    
155      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
156        cutoffGroup = createCutoffGroup(mol, *fai);
157        mol->addCutoffGroup(cutoffGroup);
# Line 192 | Line 176 | namespace oopse {
176  
177      if (atomType == NULL) {
178        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
179 <              stamp->getType());
179 >              stamp->getType().c_str());
180  
181        painCave.isFatal = 1;
182        simError();
# Line 236 | Line 220 | namespace oopse {
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));    
223 >      atom = mol->getAtomAt(rbStamp->getMemberAt(i));
224 >      atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));    
225        rb->addAtom(atom, atomStamp);
226      }
227  
# Line 283 | Line 267 | namespace oopse {
267    }    
268  
269    Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
270 <    bool isGhostBend = false;
271 <    int ghostIndex;
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 <    
290 <    //
291 <    if (stamp->haveExtras()){
292 <      LinkedAssign* extras = stamp->getExtras();
293 <      LinkedAssign* currentExtra = extras;
277 >      assert( atomA && atomB && atomC);
278  
279 <      while (currentExtra != NULL){
296 <        if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
297 <          switch (currentExtra->getType()){
298 <          case 0:
299 <            ghostIndex = currentExtra->getInt();
300 <            isGhostBend = true;
301 <            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();
308 <          }
309 <        } else{
310 <          sprintf(painCave.errMsg,
311 <                  "SimSetup Error: unhandled bend assignment:\n");
312 <          painCave.isFatal = 1;
313 <          simError();
314 <        }
315 <        currentExtra = currentExtra->getNext();
316 <      }
317 <        
318 <    }
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 <    if (isGhostBend) {
288 <
322 <      int indexA = stamp->getA();
323 <      int indexB= stamp->getB();
324 <
325 <      assert(indexA != indexB);
326 <
327 <      int normalIndex;
328 <      if (indexA == ghostIndex) {
329 <        normalIndex = indexB;
330 <      } else if (indexB == ghostIndex) {
331 <        normalIndex = indexA;
287 >        painCave.isFatal = 1;
288 >        simError();
289        }
290 <        
290 >
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) {
# Line 351 | Line 312 | namespace oopse {
312          simError();
313        }
314          
315 <      return new GhostBend(normalAtom, ghostAtom, bendType);      
315 >      bend = new GhostBend(normalAtom, ghostAtom, bendType);      
316  
317 <    } else {
318 <            
319 <      Atom* atomA = mol->getAtomAt(stamp->getA());
359 <      Atom* atomB = mol->getAtomAt(stamp->getB());
360 <      Atom* atomC = mol->getAtomAt(stamp->getC());
361 <
362 <      assert( atomA && atomB && atomC);
363 <        
364 <      BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
365 <
366 <      if (bendType == NULL) {
367 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
368 <                atomA->getType().c_str(),
369 <                atomB->getType().c_str(),
370 <                atomC->getType().c_str());
371 <
372 <        painCave.isFatal = 1;
373 <        simError();
374 <      }
375 <
376 <      return new Bend(atomA, atomB, atomC, bendType);      
377 <    }
317 >    }
318 >    
319 >    return bend;
320    }    
321  
322    Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
323  
324 <    Atom* atomA = mol->getAtomAt(stamp->getA());
325 <    Atom* atomB = mol->getAtomAt(stamp->getB());
326 <    Atom* atomC = mol->getAtomAt(stamp->getC());
327 <    Torsion* torsion;
324 >    Torsion* torsion = NULL;
325 >    std::vector<int> torsionAtoms = stamp->getMembers();
326 >    if (torsionAtoms.size() < 3) {
327 >        return torsion;
328 >    }
329  
330 <    if (stamp->getD() != -1) {
331 <      Atom* atomD = mol->getAtomAt(stamp->getD());
330 >    Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
331 >    Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
332 >    Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
333  
334 +    if (torsionAtoms.size() == 4) {
335 +      Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
336 +
337        assert(atomA && atomB && atomC && atomD);
338          
339        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
# Line 407 | Line 354 | namespace oopse {
354      }
355      else {
356  
357 <      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
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;
# Line 442 | Line 389 | namespace oopse {
389      
390      nAtoms = stamp->getNMembers();
391      for (int i =0; i < nAtoms; ++i) {
392 <      atom = mol->getAtomAt(stamp->getMember(i));
392 >      atom = mol->getAtomAt(stamp->getMemberAt(i));
393        assert(atom);
394        cg->addAtom(atom);
395      }

Diff Legend

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