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 1442 by gezelter, Mon May 10 17:28:26 2010 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42    
43   /**
44   * @file MoleculeCreator.cpp
45   * @author tlin
46   * @date 11/04/2004
46 – * @time 13:44am
47   * @version 1.0
48   */
49  
# Line 54 | Line 54
54   #include "brains/MoleculeCreator.hpp"
55   #include "primitives/GhostBend.hpp"
56   #include "primitives/GhostTorsion.hpp"
57 < #include "types/DirectionalAtomType.hpp"
57 > #include "types/AtomType.hpp"
58   #include "types/FixedBondType.hpp"
59   #include "utils/simError.h"
60   #include "utils/StringUtils.hpp"
# Line 65 | Line 65 | namespace OpenMD {
65                                              MoleculeStamp *molStamp,
66                                              int stampId, int globalIndex,
67                                              LocalIndexManager* localIndexMan) {
68 <    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName());
69 <    
68 >    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName(),
69 >                                 molStamp->getRegion() );
70 >
71      //create atoms
72      Atom* atom;
73      AtomStamp* currentAtomStamp;
# Line 96 | Line 97 | namespace OpenMD {
97  
98      for (int i = 0; i < nBonds; ++i) {
99        currentBondStamp = molStamp->getBondStamp(i);
100 <      bond = createBond(ff, mol, currentBondStamp);
100 >      bond = createBond(ff, mol, currentBondStamp, localIndexMan);
101        mol->addBond(bond);
102      }
103  
# Line 106 | Line 107 | namespace OpenMD {
107      int nBends = molStamp->getNBends();
108      for (int i = 0; i < nBends; ++i) {
109        currentBendStamp = molStamp->getBendStamp(i);
110 <      bend = createBend(ff, mol, currentBendStamp);
110 >      bend = createBend(ff, mol, currentBendStamp, localIndexMan);
111        mol->addBend(bend);
112      }
113  
# Line 116 | Line 117 | namespace OpenMD {
117      int nTorsions = molStamp->getNTorsions();
118      for (int i = 0; i < nTorsions; ++i) {
119        currentTorsionStamp = molStamp->getTorsionStamp(i);
120 <      torsion = createTorsion(ff, mol, currentTorsionStamp);
120 >      torsion = createTorsion(ff, mol, currentTorsionStamp, localIndexMan);
121        mol->addTorsion(torsion);
122      }
123  
# Line 126 | Line 127 | namespace OpenMD {
127      int nInversions = molStamp->getNInversions();
128      for (int i = 0; i < nInversions; ++i) {
129        currentInversionStamp = molStamp->getInversionStamp(i);
130 <      inversion = createInversion(ff, mol, currentInversionStamp);
130 >      inversion = createInversion(ff, mol, currentInversionStamp,
131 >                                  localIndexMan);
132        if (inversion != NULL ) {
133          mol->addInversion(inversion);
134        }
# Line 138 | Line 140 | namespace OpenMD {
140      int nCutoffGroups = molStamp->getNCutoffGroups();
141      for (int i = 0; i < nCutoffGroups; ++i) {
142        currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
143 <      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
143 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp,
144 >                                      localIndexMan);
145        mol->addCutoffGroup(cutoffGroup);
146      }
147  
# Line 169 | Line 172 | namespace OpenMD {
172      // every single free atom
173      
174      for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
175 <      cutoffGroup = createCutoffGroup(mol, *fai);
175 >      cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
176        mol->addCutoffGroup(cutoffGroup);
177      }
178      //create constraints
179      createConstraintPair(mol);
180      createConstraintElem(mol);
181      
182 +    // Does this molecule stamp define a total constrained charge value?
183 +    // If so, let the created molecule know about it.
184 +
185 +    if (molStamp->haveConstrainTotalCharge() ) {
186 +      mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() );
187 +    }
188 +
189      //the construction of this molecule is finished
190      mol->complete();
191      
# Line 198 | Line 208 | namespace OpenMD {
208        painCave.isFatal = 1;
209        simError();
210      }
211 <    
211 >
212      //below code still have some kind of hard-coding smell
213      if (atomType->isDirectional()){
204 –    
205 –      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
206 –        
207 –      if (dAtomType == NULL) {
208 –        sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
214  
210 –        painCave.isFatal = 1;
211 –        simError();
212 –      }
213 –
215        DirectionalAtom* dAtom;
216 <      dAtom = new DirectionalAtom(dAtomType);
216 >      dAtom = new DirectionalAtom(atomType);
217        atom = dAtom;    
218      }
219      else{
# Line 227 | Line 228 | namespace OpenMD {
228    RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
229                                                Molecule* mol,
230                                                RigidBodyStamp* rbStamp,
231 <                                              LocalIndexManager* localIndexMan) {
231 >                                              LocalIndexManager* localIndexMan){
232      Atom* atom;
233      int nAtoms;
234      Vector3d refCoor;
# Line 251 | Line 252 | namespace OpenMD {
252      //set the local index of this rigid body, global index will be set later
253      rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
254  
255 <    //the rule for naming rigidbody MoleculeName_RB_Integer
256 <    //The first part is the name of the molecule
257 <    //The second part is alway fixed as "RB"
258 <    //The third part is the index of the rigidbody defined in meta-data file
259 <    //For example, Butane_RB_0 is a valid rigid body name of butane molecule
260 <    /**@todo replace itoa by lexi_cast */
255 >    // The rule for naming a rigidbody is: MoleculeName_RB_Integer
256 >    // The first part is the name of the molecule
257 >    // The second part is always fixed as "RB"
258 >    // The third part is the index of the rigidbody defined in meta-data file
259 >    // For example, Butane_RB_0 is a valid rigid body name of butane molecule
260 >
261      std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
262      rb->setType(mol->getType() + "_RB_" + s.c_str());
262 –
263      return rb;
264    }    
265  
266    Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
267 <                                    BondStamp* stamp) {
267 >                                    BondStamp* stamp,
268 >                                    LocalIndexManager* localIndexMan) {
269      BondType* bondType;
270      Atom* atomA;
271      Atom* atomB;
# Line 284 | Line 285 | namespace OpenMD {
285        painCave.isFatal = 1;
286        simError();
287      }
288 <    return new Bond(atomA, atomB, bondType);    
288 >    Bond* bond = new Bond(atomA, atomB, bondType);
289 >
290 >    //set the local index of this bond, the global index will be set later
291 >    bond->setLocalIndex(localIndexMan->getNextBondIndex());
292 >
293 >    // The rule for naming a bond is: MoleculeName_Bond_Integer
294 >    // The first part is the name of the molecule
295 >    // The second part is always fixed as "Bond"
296 >    // The third part is the index of the bond defined in meta-data file
297 >    // For example, Butane_bond_0 is a valid Bond name in a butane molecule
298 >
299 >    std::string s = OpenMD_itoa(mol->getNBonds(), 10);
300 >    bond->setName(mol->getType() + "_Bond_" + s.c_str());
301 >    return bond;    
302    }    
303    
304    Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
305 <                                    BendStamp* stamp) {
305 >                                    BendStamp* stamp,
306 >                                    LocalIndexManager* localIndexMan) {
307      Bend* bend = NULL;
308      std::vector<int> bendAtoms = stamp->getMembers();
309      if (bendAtoms.size() == 3) {
# Line 303 | Line 318 | namespace OpenMD {
318                                             atomC->getType().c_str());
319        
320        if (bendType == NULL) {
321 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
321 >        sprintf(painCave.errMsg,
322 >                "Can not find Matching Bend Type for[%s, %s, %s]",
323                  atomA->getType().c_str(),
324                  atomB->getType().c_str(),
325                  atomC->getType().c_str());
# Line 327 | Line 343 | namespace OpenMD {
343        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
344  
345        if (bendType == NULL) {
346 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
346 >        sprintf(painCave.errMsg,
347 >                "Can not find Matching Bend Type for[%s, %s, %s]",
348                  normalAtom->getType().c_str(),
349                  ghostAtom->getType().c_str(),
350                  "GHOST");
# Line 339 | Line 356 | namespace OpenMD {
356        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
357        
358      }
359 <    
359 >
360 >    //set the local index of this bend, the global index will be set later
361 >    bend->setLocalIndex(localIndexMan->getNextBendIndex());
362 >
363 >    // The rule for naming a bend is: MoleculeName_Bend_Integer
364 >    // The first part is the name of the molecule
365 >    // The second part is always fixed as "Bend"
366 >    // The third part is the index of the bend defined in meta-data file
367 >    // For example, Butane_Bend_0 is a valid Bend name in a butane molecule
368 >
369 >    std::string s = OpenMD_itoa(mol->getNBends(), 10);
370 >    bend->setName(mol->getType() + "_Bend_" + s.c_str());    
371      return bend;
372    }    
373  
374    Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
375 <                                          TorsionStamp* stamp) {
375 >                                          TorsionStamp* stamp,
376 >                                          LocalIndexManager* localIndexMan) {
377  
378      Torsion* torsion = NULL;
379      std::vector<int> torsionAtoms = stamp->getMembers();
# Line 366 | Line 395 | namespace OpenMD {
395                                                      atomC->getType(),
396                                                      atomD->getType());
397        if (torsionType == NULL) {
398 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
398 >        sprintf(painCave.errMsg,
399 >                "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
400                  atomA->getType().c_str(),
401                  atomB->getType().c_str(),
402                  atomC->getType().c_str(),
# Line 403 | Line 433 | namespace OpenMD {
433        
434        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
435      }
436 +
437 +    //set the local index of this torsion, the global index will be set later
438 +    torsion->setLocalIndex(localIndexMan->getNextTorsionIndex());
439      
440 +    // The rule for naming a torsion is: MoleculeName_Torsion_Integer
441 +    // The first part is the name of the molecule
442 +    // The second part is always fixed as "Torsion"
443 +    // The third part is the index of the torsion defined in meta-data file
444 +    // For example, Butane_Torsion_0 is a valid Torsion name in a
445 +    // butane molecule
446 +
447 +    std::string s = OpenMD_itoa(mol->getNTorsions(), 10);
448 +    torsion->setName(mol->getType() + "_Torsion_" + s.c_str());
449      return torsion;
450    }    
451  
452    Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
453 <                                              InversionStamp* stamp) {
453 >                                              InversionStamp* stamp,
454 >                                              LocalIndexManager* localIndexMan) {
455      
456      Inversion* inversion = NULL;
457      int center = stamp->getCenter();
# Line 444 | Line 487 | namespace OpenMD {
487      } else {
488        
489        inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
490 +
491 +      // set the local index of this inversion, the global index will
492 +      // be set later
493 +      inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
494 +
495 +      // The rule for naming an inversion is: MoleculeName_Inversion_Integer
496 +      // The first part is the name of the molecule
497 +      // The second part is always fixed as "Inversion"
498 +      // The third part is the index of the inversion defined in meta-data file
499 +      // For example, Benzene_Inversion_0 is a valid Inversion name in a
500 +      // Benzene molecule
501 +
502 +      std::string s = OpenMD_itoa(mol->getNInversions(), 10);
503 +      inversion->setName(mol->getType() + "_Inversion_" + s.c_str());
504        return inversion;
505      }
506    }
507    
508  
509 <  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
509 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol,
510 >                                                  CutoffGroupStamp* stamp,
511 >                                                  LocalIndexManager* localIndexMan) {
512      int nAtoms;
513      CutoffGroup* cg;
514      Atom* atom;
# Line 461 | Line 520 | namespace OpenMD {
520        assert(atom);
521        cg->addAtom(atom);
522      }
523 <
523 >    
524 >    //set the local index of this cutoffGroup, global index will be set later
525 >    cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
526 >    
527      return cg;
528    }    
529 <
530 <  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
529 >  
530 >  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom,
531 >                                                  LocalIndexManager* localIndexMan) {
532      CutoffGroup* cg;
533      cg  = new CutoffGroup();
534      cg->addAtom(atom);
535 +
536 +    //set the local index of this cutoffGroup, global index will be set later
537 +    cg->setLocalIndex(localIndexMan->getNextCutoffGroupIndex());
538 +
539      return cg;
540    }
541  
# Line 481 | Line 548 | namespace OpenMD {
548          
549        BondType* bt = bond->getBondType();
550  
484 –      //class Parent1 {};
485 –      //class Child1 : public Parent {};
486 –      //class Child2 : public Parent {};
487 –      //Child1* ch1 = new Child1();
488 –      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
489 –      //the dynamic_cast is succeed in above line. A compiler bug?        
490 –
551        if (typeid(FixedBondType) == typeid(*bt)) {
552          FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
553  

Diff Legend

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