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 647 by tim, Wed Oct 5 19:34:01 2005 UTC vs.
Revision 1290 by cli2, Wed Sep 10 19:51:45 2008 UTC

# Line 48 | Line 48
48   */
49  
50   #include <cassert>
51 + #include <typeinfo>
52   #include <set>
53  
54   #include "brains/MoleculeCreator.hpp"
# Line 60 | Line 61 | namespace oopse {
61  
62   namespace oopse {
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());
64 >  Molecule* MoleculeCreator::createMolecule(ForceField* ff,
65 >                                            MoleculeStamp *molStamp,
66 >                                            int stampId, int globalIndex,
67 >                                            LocalIndexManager* localIndexMan) {
68 >    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName());
69      
70      //create atoms
71      Atom* atom;
72      AtomStamp* currentAtomStamp;
73      int nAtom = molStamp->getNAtoms();
74      for (int i = 0; i < nAtom; ++i) {
75 <      currentAtomStamp = molStamp->getAtom(i);
75 >      currentAtomStamp = molStamp->getAtomStamp(i);
76        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
77        mol->addAtom(atom);
78      }
# Line 81 | Line 83 | namespace oopse {
83      int nRigidbodies = molStamp->getNRigidBodies();
84  
85      for (int i = 0; i < nRigidbodies; ++i) {
86 <      currentRigidBodyStamp = molStamp->getRigidBody(i);
87 <      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
86 >      currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
87 >      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp,
88 >                           localIndexMan);
89        mol->addRigidBody(rb);
90      }
91 <
91 >    
92      //create bonds
93      Bond* bond;
94      BondStamp* currentBondStamp;
95      int nBonds = molStamp->getNBonds();
96  
97      for (int i = 0; i < nBonds; ++i) {
98 <      currentBondStamp = molStamp->getBond(i);
98 >      currentBondStamp = molStamp->getBondStamp(i);
99        bond = createBond(ff, mol, currentBondStamp);
100        mol->addBond(bond);
101      }
# Line 102 | Line 105 | namespace oopse {
105      BendStamp* currentBendStamp;
106      int nBends = molStamp->getNBends();
107      for (int i = 0; i < nBends; ++i) {
108 <      currentBendStamp = molStamp->getBend(i);
108 >      currentBendStamp = molStamp->getBendStamp(i);
109        bend = createBend(ff, mol, currentBendStamp);
110        mol->addBend(bend);
111      }
# Line 112 | Line 115 | namespace oopse {
115      TorsionStamp* currentTorsionStamp;
116      int nTorsions = molStamp->getNTorsions();
117      for (int i = 0; i < nTorsions; ++i) {
118 <      currentTorsionStamp = molStamp->getTorsion(i);
118 >      currentTorsionStamp = molStamp->getTorsionStamp(i);
119        torsion = createTorsion(ff, mol, currentTorsionStamp);
120        mol->addTorsion(torsion);
121      }
122  
123 +    //create inversions
124 +    Inversion* inversion;
125 +    InversionStamp* currentInversionStamp;
126 +    int nInversions = molStamp->getNInversions();
127 +    for (int i = 0; i < nInversions; ++i) {
128 +      currentInversionStamp = molStamp->getInversionStamp(i);
129 +      inversion = createInversion(ff, mol, currentInversionStamp);
130 +      if (inversion != NULL ) {
131 +        mol->addInversion(inversion);
132 +      }
133 +    }
134 +
135      //create cutoffGroups
136      CutoffGroup* cutoffGroup;
137      CutoffGroupStamp* currentCutoffGroupStamp;
138      int nCutoffGroups = molStamp->getNCutoffGroups();
139      for (int i = 0; i < nCutoffGroups; ++i) {
140 <      currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
140 >      currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
141        cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
142        mol->addCutoffGroup(cutoffGroup);
143      }
# Line 140 | Line 155 | namespace oopse {
155      Molecule::CutoffGroupIterator ci;
156      CutoffGroup* cg;
157      
158 <    for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
159 <
158 >    for (cg = mol->beginCutoffGroup(ci); cg != NULL;
159 >         cg = mol->nextCutoffGroup(ci)) {
160 >      
161        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
162 <           //erase the atoms belong to cutoff groups from freeAtoms vector
163 <           freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom), freeAtoms.end());
164 <      }
165 <
162 >        //erase the atoms belong to cutoff groups from freeAtoms vector
163 >        freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom),
164 >                        freeAtoms.end());
165 >      }      
166      }      
167      
168 <    //loop over the free atoms and then create one cutoff group for every single free atom
168 >    // loop over the free atoms and then create one cutoff group for
169 >    // every single free atom
170      
171      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
172        cutoffGroup = createCutoffGroup(mol, *fai);
# Line 161 | Line 178 | namespace oopse {
178      
179      //the construction of this molecule is finished
180      mol->complete();
181 <
181 >    
182      return mol;
183    }    
184  
185  
186 <  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp,
187 <                                    LocalIndexManager* localIndexMan) {
186 >  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol,
187 >                                    AtomStamp* stamp,
188 >                                    LocalIndexManager* localIndexMan) {
189      AtomType * atomType;
190      Atom* atom;
191  
192      atomType =  ff->getAtomType(stamp->getType());
193 <
193 >    
194      if (atomType == NULL) {
195        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
196 <              stamp->getType());
196 >              stamp->getType().c_str());
197  
198        painCave.isFatal = 1;
199        simError();
# Line 205 | Line 223 | namespace oopse {
223  
224      return atom;
225    }
226 <
227 <  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
226 >  
227 >  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
228 >                                              Molecule* mol,
229                                                RigidBodyStamp* rbStamp,
230                                                LocalIndexManager* localIndexMan) {
231      Atom* atom;
# Line 217 | Line 236 | namespace oopse {
236      RigidBody* rb = new RigidBody();
237      nAtoms = rbStamp->getNMembers();    
238      for (int i = 0; i < nAtoms; ++i) {
239 <      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
240 <      //It is not the same as local index of atom which is the index of atom at DataStorage class
241 <      atom = mol->getAtomAt(rbStamp->getMember(i));
242 <      atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
239 >      //rbStamp->getMember(i) return the local index of current atom
240 >      //inside the molecule.  It is not the same as local index of
241 >      //atom which is the index of atom at DataStorage class
242 >      atom = mol->getAtomAt(rbStamp->getMemberAt(i));
243 >      atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));    
244        rb->addAtom(atom, atomStamp);
245      }
246  
247 <    //after all of the atoms are added, we need to calculate the reference coordinates
247 >    //after all of the atoms are added, we need to calculate the
248 >    //reference coordinates
249      rb->calcRefCoords();
250  
251      //set the local index of this rigid body, global index will be set later
# Line 242 | Line 263 | namespace oopse {
263      return rb;
264    }    
265  
266 <  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
266 >  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
267 >                                    BondStamp* stamp) {
268      BondType* bondType;
269      Atom* atomA;
270      Atom* atomB;
271 <
271 >    
272      atomA = mol->getAtomAt(stamp->getA());
273      atomB = mol->getAtomAt(stamp->getB());
274 <
274 >    
275      assert( atomA && atomB);
276      
277      bondType = ff->getBondType(atomA->getType(), atomB->getType());
# Line 258 | Line 280 | namespace oopse {
280        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
281                atomA->getType().c_str(),
282                atomB->getType().c_str());
283 <
283 >      
284        painCave.isFatal = 1;
285        simError();
286      }
287      return new Bond(atomA, atomB, bondType);    
288    }    
289 <
290 <  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
291 <    bool isGhostBend = false;
292 <    int ghostIndex;
293 <
294 <    
295 <    //
296 <    if (stamp->haveExtras()){
297 <      LinkedAssign* extras = stamp->getExtras();
298 <      LinkedAssign* currentExtra = extras;
299 <
300 <      while (currentExtra != NULL){
301 <        if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
302 <          switch (currentExtra->getType()){
303 <          case 0:
304 <            ghostIndex = currentExtra->getInt();
305 <            isGhostBend = true;
306 <            break;
307 <
308 <          default:
309 <            sprintf(painCave.errMsg,
288 <                    "SimSetup Error: ghostVectorSource must be an int.\n");
289 <            painCave.isFatal = 1;
290 <            simError();
291 <          }
292 <        } else{
293 <          sprintf(painCave.errMsg,
294 <                  "SimSetup Error: unhandled bend assignment:\n");
295 <          painCave.isFatal = 1;
296 <          simError();
297 <        }
298 <        currentExtra = currentExtra->getNext();
299 <      }
300 <        
301 <    }
302 <
303 <    if (isGhostBend) {
304 <
305 <      int indexA = stamp->getA();
306 <      int indexB= stamp->getB();
307 <
308 <      assert(indexA != indexB);
309 <
310 <      int normalIndex;
311 <      if (indexA == ghostIndex) {
312 <        normalIndex = indexB;
313 <      } else if (indexB == ghostIndex) {
314 <        normalIndex = indexA;
315 <      }
289 >  
290 >  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
291 >                                    BendStamp* stamp) {
292 >    Bend* bend = NULL;
293 >    std::vector<int> bendAtoms = stamp->getMembers();
294 >    if (bendAtoms.size() == 3) {
295 >      Atom* atomA = mol->getAtomAt(bendAtoms[0]);
296 >      Atom* atomB = mol->getAtomAt(bendAtoms[1]);
297 >      Atom* atomC = mol->getAtomAt(bendAtoms[2]);
298 >      
299 >      assert( atomA && atomB && atomC);
300 >      
301 >      BendType* bendType = ff->getBendType(atomA->getType().c_str(),
302 >                                           atomB->getType().c_str(),
303 >                                           atomC->getType().c_str());
304 >      
305 >      if (bendType == NULL) {
306 >        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
307 >                atomA->getType().c_str(),
308 >                atomB->getType().c_str(),
309 >                atomC->getType().c_str());
310          
311 +        painCave.isFatal = 1;
312 +        simError();
313 +      }
314 +      
315 +      bend = new Bend(atomA, atomB, atomC, bendType);
316 +    } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
317 +      int ghostIndex = stamp->getGhostVectorSource();
318 +      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
319        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
320        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
321        if (ghostAtom == NULL) {
# Line 333 | Line 335 | namespace oopse {
335          painCave.isFatal = 1;
336          simError();
337        }
338 <        
339 <      return new GhostBend(normalAtom, ghostAtom, bendType);      
338 >      
339 >      bend = new GhostBend(normalAtom, ghostAtom, bendType);      
340 >      
341 >    }
342 >    
343 >    return bend;
344 >  }    
345  
346 <    } else {
347 <            
341 <      Atom* atomA = mol->getAtomAt(stamp->getA());
342 <      Atom* atomB = mol->getAtomAt(stamp->getB());
343 <      Atom* atomC = mol->getAtomAt(stamp->getC());
346 >  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
347 >                                          TorsionStamp* stamp) {
348  
349 <      assert( atomA && atomB && atomC);
350 <        
351 <      BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
352 <
349 <      if (bendType == NULL) {
350 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
351 <                atomA->getType().c_str(),
352 <                atomB->getType().c_str(),
353 <                atomC->getType().c_str());
354 <
355 <        painCave.isFatal = 1;
356 <        simError();
357 <      }
358 <
359 <      return new Bend(atomA, atomB, atomC, bendType);      
349 >    Torsion* torsion = NULL;
350 >    std::vector<int> torsionAtoms = stamp->getMembers();
351 >    if (torsionAtoms.size() < 3) {
352 >        return torsion;
353      }
361 –  }    
354  
355 <  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
355 >    Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
356 >    Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
357 >    Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
358  
359 <    Atom* atomA = mol->getAtomAt(stamp->getA());
360 <    Atom* atomB = mol->getAtomAt(stamp->getB());
367 <    Atom* atomC = mol->getAtomAt(stamp->getC());
368 <    Torsion* torsion;
359 >    if (torsionAtoms.size() == 4) {
360 >      Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
361  
370 –    if (stamp->getD() != -1) {
371 –      Atom* atomD = mol->getAtomAt(stamp->getD());
372 –
362        assert(atomA && atomB && atomC && atomD);
363          
364 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
365 <                                                    atomC->getType(), atomD->getType());
366 <
364 >      TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
365 >                                                    atomB->getType(),
366 >                                                    atomC->getType(),
367 >                                                    atomD->getType());
368        if (torsionType == NULL) {
369          sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
370                  atomA->getType().c_str(),
371                  atomB->getType().c_str(),
372                  atomC->getType().c_str(),
373                  atomD->getType().c_str());
374 <
374 >        
375          painCave.isFatal = 1;
376          simError();
377        }
378 <        
378 >      
379        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
380      }
381      else {
382 <
383 <      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
382 >      
383 >      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
384        if (dAtom == NULL) {
385          sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
386          painCave.isFatal = 1;
387          simError();
388        }        
389 <
389 >      
390        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
391                                                      atomC->getType(), "GHOST");
392 <
392 >      
393        if (torsionType == NULL) {
394          sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
395                  atomA->getType().c_str(),
396                  atomB->getType().c_str(),
397                  atomC->getType().c_str(),
398                  "GHOST");
399 <
399 >        
400          painCave.isFatal = 1;
401          simError();
402        }
403 <        
403 >      
404        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
405      }
406 <
406 >    
407      return torsion;
408    }    
409  
410 +  Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
411 +                                              InversionStamp* stamp) {
412 +    
413 +    Inversion* inversion = NULL;
414 +    int center = stamp->getCenter();
415 +    std::vector<int> satellites = stamp->getSatellites();
416 +    if (satellites.size() != 3) {
417 +        return inversion;
418 +    }
419 +
420 +    Atom* atomA = mol->getAtomAt(center);
421 +    Atom* atomB = mol->getAtomAt(satellites[0]);
422 +    Atom* atomC = mol->getAtomAt(satellites[1]);
423 +    Atom* atomD = mol->getAtomAt(satellites[2]);
424 +      
425 +    assert(atomA && atomB && atomC && atomD);
426 +    
427 +    InversionType* inversionType = ff->getInversionType(atomA->getType(),
428 +                                                        atomB->getType(),
429 +                                                        atomC->getType(),
430 +                                                        atomD->getType());
431 +
432 +    if (inversionType == NULL) {
433 +      sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
434 +              "\t(May not be a problem: not all inversions are parametrized)\n",
435 +              atomA->getType().c_str(),
436 +              atomB->getType().c_str(),
437 +              atomC->getType().c_str(),
438 +              atomD->getType().c_str());
439 +      
440 +      painCave.isFatal = 0;
441 +      painCave.severity = OOPSE_INFO;
442 +      simError();
443 +      return NULL;
444 +    } else {
445 +      
446 +      inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
447 +      return inversion;
448 +    }
449 +  }
450 +  
451 +
452    CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
453      int nAtoms;
454      CutoffGroup* cg;
# Line 425 | Line 457 | namespace oopse {
457      
458      nAtoms = stamp->getNMembers();
459      for (int i =0; i < nAtoms; ++i) {
460 <      atom = mol->getAtomAt(stamp->getMember(i));
460 >      atom = mol->getAtomAt(stamp->getMemberAt(i));
461        assert(atom);
462        cg->addAtom(atom);
463      }

Diff Legend

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