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 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 56 | Line 56
56   #include "primitives/GhostTorsion.hpp"
57   #include "types/AtomType.hpp"
58   #include "types/FixedBondType.hpp"
59 + #include "types/BondTypeParser.hpp"
60 + #include "types/BendTypeParser.hpp"
61 + #include "types/TorsionTypeParser.hpp"
62 + #include "types/InversionTypeParser.hpp"
63   #include "utils/simError.h"
64   #include "utils/StringUtils.hpp"
65  
# Line 65 | Line 69 | namespace OpenMD {
69                                              MoleculeStamp *molStamp,
70                                              int stampId, int globalIndex,
71                                              LocalIndexManager* localIndexMan) {
72 <    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName());
72 >    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName(),
73 >                                 molStamp->getRegion() );
74  
75      //create atoms
76      Atom* atom;
# Line 95 | Line 100 | namespace OpenMD {
100      int nBonds = molStamp->getNBonds();
101  
102      for (int i = 0; i < nBonds; ++i) {
103 <      currentBondStamp = molStamp->getBondStamp(i);
104 <      bond = createBond(ff, mol, currentBondStamp);
103 >      currentBondStamp = molStamp->getBondStamp(i);        
104 >      bond = createBond(ff, mol, currentBondStamp, localIndexMan);
105        mol->addBond(bond);
106      }
107  
# Line 106 | Line 111 | namespace OpenMD {
111      int nBends = molStamp->getNBends();
112      for (int i = 0; i < nBends; ++i) {
113        currentBendStamp = molStamp->getBendStamp(i);
114 <      bend = createBend(ff, mol, currentBendStamp);
114 >      bend = createBend(ff, mol, currentBendStamp, localIndexMan);
115        mol->addBend(bend);
116      }
117  
# Line 116 | Line 121 | namespace OpenMD {
121      int nTorsions = molStamp->getNTorsions();
122      for (int i = 0; i < nTorsions; ++i) {
123        currentTorsionStamp = molStamp->getTorsionStamp(i);
124 <      torsion = createTorsion(ff, mol, currentTorsionStamp);
124 >      torsion = createTorsion(ff, mol, currentTorsionStamp, localIndexMan);
125        mol->addTorsion(torsion);
126      }
127  
# Line 126 | Line 131 | namespace OpenMD {
131      int nInversions = molStamp->getNInversions();
132      for (int i = 0; i < nInversions; ++i) {
133        currentInversionStamp = molStamp->getInversionStamp(i);
134 <      inversion = createInversion(ff, mol, currentInversionStamp);
134 >      inversion = createInversion(ff, mol, currentInversionStamp,
135 >                                  localIndexMan);
136        if (inversion != NULL ) {
137          mol->addInversion(inversion);
138        }
# Line 138 | Line 144 | namespace OpenMD {
144      int nCutoffGroups = molStamp->getNCutoffGroups();
145      for (int i = 0; i < nCutoffGroups; ++i) {
146        currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
147 <      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp, localIndexMan);
147 >      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp,
148 >                                      localIndexMan);
149        mol->addCutoffGroup(cutoffGroup);
150      }
151  
# Line 172 | Line 179 | namespace OpenMD {
179        cutoffGroup = createCutoffGroup(mol, *fai, localIndexMan);
180        mol->addCutoffGroup(cutoffGroup);
181      }
182 <    //create constraints
182 >
183 >    //create bonded constraintPairs:
184      createConstraintPair(mol);
185 +
186 +    //create non-bonded constraintPairs
187 +    for (std::size_t i = 0; i < molStamp->getNConstraints(); ++i) {
188 +      ConstraintStamp* cStamp = molStamp->getConstraintStamp(i);
189 +      Atom* atomA;
190 +      Atom* atomB;
191 +      
192 +      atomA = mol->getAtomAt(cStamp->getA());
193 +      atomB = mol->getAtomAt(cStamp->getB());
194 +      assert( atomA && atomB );
195 +      
196 +      bool printConstraintForce = false;
197 +
198 +      if (cStamp->havePrintConstraintForce()) {
199 +        printConstraintForce = cStamp->getPrintConstraintForce();
200 +      }
201 +
202 +      if (!cStamp->haveConstrainedDistance()) {
203 +        sprintf(painCave.errMsg,
204 +                "Constraint Error: A non-bond constraint was specified\n"
205 +                "\twithout providing a value for the constrainedDistance.\n");
206 +        painCave.isFatal = 1;
207 +        simError();      
208 +      } else {
209 +        RealType distance = cStamp->getConstrainedDistance();
210 +        ConstraintElem* consElemA = new ConstraintElem(atomA);
211 +        ConstraintElem* consElemB = new ConstraintElem(atomB);
212 +        ConstraintPair* cPair = new ConstraintPair(consElemA, consElemB,
213 +                                                   distance,
214 +                                                   printConstraintForce);
215 +        mol->addConstraintPair(cPair);
216 +      }  
217 +    }
218 +
219 +    // now create the constraint elements:
220 +
221      createConstraintElem(mol);
222      
223      // Does this molecule stamp define a total constrained charge value?
# Line 183 | Line 227 | namespace OpenMD {
227        mol->setConstrainTotalCharge( molStamp->getConstrainTotalCharge() );
228      }
229  
230 <    //the construction of this molecule is finished
230 >    // The construction of this molecule is finished:
231      mol->complete();
232      
233      return mol;
# Line 225 | Line 269 | namespace OpenMD {
269    RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp,
270                                                Molecule* mol,
271                                                RigidBodyStamp* rbStamp,
272 <                                              LocalIndexManager* localIndexMan) {
272 >                                              LocalIndexManager* localIndexMan){
273      Atom* atom;
274      int nAtoms;
275      Vector3d refCoor;
# Line 249 | Line 293 | namespace OpenMD {
293      //set the local index of this rigid body, global index will be set later
294      rb->setLocalIndex(localIndexMan->getNextRigidBodyIndex());
295  
296 <    //the rule for naming rigidbody MoleculeName_RB_Integer
297 <    //The first part is the name of the molecule
298 <    //The second part is alway fixed as "RB"
299 <    //The third part is the index of the rigidbody defined in meta-data file
300 <    //For example, Butane_RB_0 is a valid rigid body name of butane molecule
301 <    /**@todo replace itoa by lexi_cast */
296 >    // The rule for naming a rigidbody is: MoleculeName_RB_Integer
297 >    // The first part is the name of the molecule
298 >    // The second part is always fixed as "RB"
299 >    // The third part is the index of the rigidbody defined in meta-data file
300 >    // For example, Butane_RB_0 is a valid rigid body name of butane molecule
301 >
302      std::string s = OpenMD_itoa(mol->getNRigidBodies(), 10);
303      rb->setType(mol->getType() + "_RB_" + s.c_str());
260 –
304      return rb;
305    }    
306  
307    Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol,
308 <                                    BondStamp* stamp) {
309 <    BondType* bondType;
308 >                                    BondStamp* stamp,
309 >                                    LocalIndexManager* localIndexMan) {
310 >    BondTypeParser btParser;        
311 >    BondType* bondType = NULL;
312      Atom* atomA;
313      Atom* atomB;
314      
# Line 271 | Line 316 | namespace OpenMD {
316      atomB = mol->getAtomAt(stamp->getB());
317      
318      assert( atomA && atomB);
274 –    
275 –    bondType = ff->getBondType(atomA->getType(), atomB->getType());
319  
320 <    if (bondType == NULL) {
321 <      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
322 <              atomA->getType().c_str(),
323 <              atomB->getType().c_str());
324 <      
325 <      painCave.isFatal = 1;
326 <      simError();
320 >    if (stamp->hasOverride()) {
321 >
322 >      try {
323 >        bondType = btParser.parseTypeAndPars(stamp->getOverrideType(),
324 >                                             stamp->getOverridePars() );
325 >      }
326 >      catch( OpenMDException e) {
327 >        sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
328 >                "for molecule %s\n",
329 >                e.what(), mol->getType().c_str() );
330 >        painCave.isFatal = 1;
331 >        simError();
332 >      }
333 >
334 >    } else {
335 >      bondType = ff->getBondType(atomA->getType(), atomB->getType());
336 >
337 >      if (bondType == NULL) {
338 >        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
339 >                atomA->getType().c_str(),
340 >                atomB->getType().c_str());
341 >        
342 >        painCave.isFatal = 1;
343 >        simError();
344 >      }
345      }
346 <    return new Bond(atomA, atomB, bondType);    
346 >    
347 >    Bond* bond = new Bond(atomA, atomB, bondType);
348 >
349 >    //set the local index of this bond, the global index will be set later
350 >    bond->setLocalIndex(localIndexMan->getNextBondIndex());
351 >
352 >    // The rule for naming a bond is: MoleculeName_Bond_Integer
353 >    // The first part is the name of the molecule
354 >    // The second part is always fixed as "Bond"
355 >    // The third part is the index of the bond defined in meta-data file
356 >    // For example, Butane_bond_0 is a valid Bond name in a butane molecule
357 >
358 >    std::string s = OpenMD_itoa(mol->getNBonds(), 10);
359 >    bond->setName(mol->getType() + "_Bond_" + s.c_str());
360 >    return bond;    
361    }    
362    
363    Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol,
364 <                                    BendStamp* stamp) {
365 <    Bend* bend = NULL;
366 <    std::vector<int> bendAtoms = stamp->getMembers();
364 >                                    BendStamp* stamp,
365 >                                    LocalIndexManager* localIndexMan) {
366 >    BendTypeParser btParser;
367 >    BendType* bendType = NULL;
368 >    Bend* bend = NULL;
369 >    
370 >    std::vector<int> bendAtoms = stamp->getMembers();
371      if (bendAtoms.size() == 3) {
372        Atom* atomA = mol->getAtomAt(bendAtoms[0]);
373        Atom* atomB = mol->getAtomAt(bendAtoms[1]);
374        Atom* atomC = mol->getAtomAt(bendAtoms[2]);
375        
376 <      assert( atomA && atomB && atomC);
377 <      
378 <      BendType* bendType = ff->getBendType(atomA->getType().c_str(),
300 <                                           atomB->getType().c_str(),
301 <                                           atomC->getType().c_str());
302 <      
303 <      if (bendType == NULL) {
304 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
305 <                atomA->getType().c_str(),
306 <                atomB->getType().c_str(),
307 <                atomC->getType().c_str());
376 >      assert( atomA && atomB && atomC );
377 >
378 >      if (stamp->hasOverride()) {
379          
380 <        painCave.isFatal = 1;
381 <        simError();
380 >        try {
381 >          bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
382 >                                               stamp->getOverridePars() );
383 >        }
384 >        catch( OpenMDException e) {
385 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
386 >                  "for molecule %s\n",
387 >                  e.what(), mol->getType().c_str() );
388 >          painCave.isFatal = 1;
389 >          simError();
390 >        }
391 >      } else {
392 >        
393 >        bendType = ff->getBendType(atomA->getType().c_str(),
394 >                                             atomB->getType().c_str(),
395 >                                             atomC->getType().c_str());
396 >      
397 >        if (bendType == NULL) {
398 >          sprintf(painCave.errMsg,
399 >                  "Can not find Matching Bend Type for[%s, %s, %s]",
400 >                  atomA->getType().c_str(),
401 >                  atomB->getType().c_str(),
402 >                  atomC->getType().c_str());
403 >          
404 >          painCave.isFatal = 1;
405 >          simError();
406 >        }
407        }
408        
409        bend = new Bend(atomA, atomB, atomC, bendType);
410 +      
411      } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
412        int ghostIndex = stamp->getGhostVectorSource();
413 <      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1];
413 >      int normalIndex = ghostIndex != bendAtoms[0] ?
414 >        bendAtoms[0] : bendAtoms[1];
415        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
416        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
417        if (ghostAtom == NULL) {
# Line 321 | Line 419 | namespace OpenMD {
419          painCave.isFatal = 1;
420          simError();
421        }
324 –                
325 –      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
422  
423 <      if (bendType == NULL) {
424 <        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
425 <                normalAtom->getType().c_str(),
426 <                ghostAtom->getType().c_str(),
427 <                "GHOST");
428 <
429 <        painCave.isFatal = 1;
430 <        simError();
423 >      if (stamp->hasOverride()) {
424 >        
425 >        try {
426 >          bendType = btParser.parseTypeAndPars(stamp->getOverrideType(),
427 >                                               stamp->getOverridePars() );
428 >        }
429 >        catch( OpenMDException e) {
430 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
431 >                  "for molecule %s\n",
432 >                  e.what(), mol->getType().c_str() );
433 >          painCave.isFatal = 1;
434 >          simError();
435 >        }
436 >      } else {
437 >      
438 >        bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(),
439 >                                   "GHOST");
440 >        
441 >        if (bendType == NULL) {
442 >          sprintf(painCave.errMsg,
443 >                  "Can not find Matching Bend Type for[%s, %s, %s]",
444 >                  normalAtom->getType().c_str(),
445 >                  ghostAtom->getType().c_str(),
446 >                  "GHOST");
447 >          
448 >          painCave.isFatal = 1;
449 >          simError();
450 >        }
451        }
452        
453        bend = new GhostBend(normalAtom, ghostAtom, bendType);      
454        
455      }
456      
457 +    //set the local index of this bend, the global index will be set later
458 +    bend->setLocalIndex(localIndexMan->getNextBendIndex());
459 +    
460 +    // The rule for naming a bend is: MoleculeName_Bend_Integer
461 +    // The first part is the name of the molecule
462 +    // The second part is always fixed as "Bend"
463 +    // The third part is the index of the bend defined in meta-data file
464 +    // For example, Butane_Bend_0 is a valid Bend name in a butane molecule
465 +    
466 +    std::string s = OpenMD_itoa(mol->getNBends(), 10);
467 +    bend->setName(mol->getType() + "_Bend_" + s.c_str());    
468      return bend;
469    }    
470  
471    Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol,
472 <                                          TorsionStamp* stamp) {
472 >                                          TorsionStamp* stamp,
473 >                                          LocalIndexManager* localIndexMan) {
474  
475 +    TorsionTypeParser ttParser;
476 +    TorsionType* torsionType = NULL;
477      Torsion* torsion = NULL;
478 +
479      std::vector<int> torsionAtoms = stamp->getMembers();
480      if (torsionAtoms.size() < 3) {
481          return torsion;
# Line 357 | Line 488 | namespace OpenMD {
488      if (torsionAtoms.size() == 4) {
489        Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
490  
491 <      assert(atomA && atomB && atomC && atomD);
491 >      assert(atomA && atomB && atomC && atomD );
492 >
493 >      if (stamp->hasOverride()) {
494          
495 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(),
496 <                                                    atomB->getType(),
497 <                                                    atomC->getType(),
498 <                                                    atomD->getType());
499 <      if (torsionType == NULL) {
500 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
501 <                atomA->getType().c_str(),
502 <                atomB->getType().c_str(),
503 <                atomC->getType().c_str(),
504 <                atomD->getType().c_str());
495 >        try {
496 >          torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
497 >                                                  stamp->getOverridePars() );
498 >        }
499 >        catch( OpenMDException e) {
500 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
501 >                  "for molecule %s\n",
502 >                  e.what(), mol->getType().c_str() );
503 >          painCave.isFatal = 1;
504 >          simError();
505 >        }
506 >      } else {
507 >
508          
509 <        painCave.isFatal = 1;
510 <        simError();
509 >        torsionType = ff->getTorsionType(atomA->getType(),
510 >                                         atomB->getType(),
511 >                                         atomC->getType(),
512 >                                         atomD->getType());
513 >        if (torsionType == NULL) {
514 >          sprintf(painCave.errMsg,
515 >                  "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
516 >                  atomA->getType().c_str(),
517 >                  atomB->getType().c_str(),
518 >                  atomC->getType().c_str(),
519 >                  atomD->getType().c_str());
520 >          
521 >          painCave.isFatal = 1;
522 >          simError();
523 >        }
524        }
525        
526        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);      
527 <    }
379 <    else {
527 >    } else {
528        
529        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
530        if (dAtom == NULL) {
# Line 384 | Line 532 | namespace OpenMD {
532          painCave.isFatal = 1;
533          simError();
534        }        
535 <      
536 <      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
389 <                                                    atomC->getType(), "GHOST");
390 <      
391 <      if (torsionType == NULL) {
392 <        sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
393 <                atomA->getType().c_str(),
394 <                atomB->getType().c_str(),
395 <                atomC->getType().c_str(),
396 <                "GHOST");
535 >
536 >      if (stamp->hasOverride()) {
537          
538 <        painCave.isFatal = 1;
539 <        simError();
540 <      }
538 >        try {
539 >          torsionType = ttParser.parseTypeAndPars(stamp->getOverrideType(),
540 >                                                  stamp->getOverridePars() );
541 >        }
542 >        catch( OpenMDException e) {
543 >          sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
544 >                  "for molecule %s\n",
545 >                  e.what(), mol->getType().c_str() );
546 >          painCave.isFatal = 1;
547 >          simError();
548 >        }
549 >      } else {              
550 >        torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(),
551 >                                         atomC->getType(), "GHOST");
552        
553 +        if (torsionType == NULL) {
554 +          sprintf(painCave.errMsg,
555 +                  "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
556 +                  atomA->getType().c_str(),
557 +                  atomB->getType().c_str(),
558 +                  atomC->getType().c_str(),
559 +                  "GHOST");
560 +          
561 +          painCave.isFatal = 1;
562 +          simError();
563 +        }
564 +      }
565 +
566        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);              
567      }
568 +
569 +    //set the local index of this torsion, the global index will be set later
570 +    torsion->setLocalIndex(localIndexMan->getNextTorsionIndex());
571      
572 +    // The rule for naming a torsion is: MoleculeName_Torsion_Integer
573 +    // The first part is the name of the molecule
574 +    // The second part is always fixed as "Torsion"
575 +    // The third part is the index of the torsion defined in meta-data file
576 +    // For example, Butane_Torsion_0 is a valid Torsion name in a
577 +    // butane molecule
578 +
579 +    std::string s = OpenMD_itoa(mol->getNTorsions(), 10);
580 +    torsion->setName(mol->getType() + "_Torsion_" + s.c_str());
581      return torsion;
582    }    
583  
584    Inversion* MoleculeCreator::createInversion(ForceField* ff, Molecule* mol,
585 <                                              InversionStamp* stamp) {
586 <    
585 >                                              InversionStamp* stamp,
586 >                                              LocalIndexManager* localIndexMan) {
587 >
588 >    InversionTypeParser itParser;
589 >    InversionType* inversionType = NULL;
590      Inversion* inversion = NULL;
591 +    
592      int center = stamp->getCenter();
593      std::vector<int> satellites = stamp->getSatellites();
594      if (satellites.size() != 3) {
# Line 421 | Line 601 | namespace OpenMD {
601      Atom* atomD = mol->getAtomAt(satellites[2]);
602        
603      assert(atomA && atomB && atomC && atomD);
424 –    
425 –    InversionType* inversionType = ff->getInversionType(atomA->getType(),
426 –                                                        atomB->getType(),
427 –                                                        atomC->getType(),
428 –                                                        atomD->getType());
604  
605 <    if (inversionType == NULL) {
431 <      sprintf(painCave.errMsg, "No Matching Inversion Type for[%s, %s, %s, %s]\n"
432 <              "\t(May not be a problem: not all inversions are parametrized)\n",
433 <              atomA->getType().c_str(),
434 <              atomB->getType().c_str(),
435 <              atomC->getType().c_str(),
436 <              atomD->getType().c_str());
605 >    if (stamp->hasOverride()) {
606        
607 <      painCave.isFatal = 0;
608 <      painCave.severity = OPENMD_INFO;
609 <      simError();
610 <      return NULL;
607 >      try {
608 >        inversionType = itParser.parseTypeAndPars(stamp->getOverrideType(),
609 >                                                  stamp->getOverridePars() );
610 >      }
611 >      catch( OpenMDException e) {
612 >        sprintf(painCave.errMsg, "MoleculeCreator Error: %s "
613 >                "for molecule %s\n",
614 >                e.what(), mol->getType().c_str() );
615 >        painCave.isFatal = 1;
616 >        simError();
617 >      }
618      } else {
619        
620 +      inversionType = ff->getInversionType(atomA->getType(),
621 +                                           atomB->getType(),
622 +                                           atomC->getType(),
623 +                                           atomD->getType());
624 +      
625 +      if (inversionType == NULL) {
626 +        sprintf(painCave.errMsg,
627 +                "No Matching Inversion Type for[%s, %s, %s, %s]\n"
628 +                "\t(May not be a problem: not all inversions are parametrized)\n",
629 +                atomA->getType().c_str(),
630 +                atomB->getType().c_str(),
631 +                atomC->getType().c_str(),
632 +                atomD->getType().c_str());
633 +        
634 +        painCave.isFatal = 0;
635 +        painCave.severity = OPENMD_INFO;
636 +        simError();
637 +      }
638 +    }
639 +    if (inversionType != NULL) {
640 +      
641        inversion = new Inversion(atomA, atomB, atomC, atomD, inversionType);
642 +      
643 +      // set the local index of this inversion, the global index will
644 +      // be set later
645 +      inversion->setLocalIndex(localIndexMan->getNextInversionIndex());
646 +      
647 +      // The rule for naming an inversion is: MoleculeName_Inversion_Integer
648 +      // The first part is the name of the molecule
649 +      // The second part is always fixed as "Inversion"
650 +      // The third part is the index of the inversion defined in meta-data file
651 +      // For example, Benzene_Inversion_0 is a valid Inversion name in a
652 +      // Benzene molecule
653 +
654 +      std::string s = OpenMD_itoa(mol->getNInversions(), 10);
655 +      inversion->setName(mol->getType() + "_Inversion_" + s.c_str());
656        return inversion;
657 +    } else {
658 +      return NULL;
659      }
660    }
661    
# Line 485 | Line 698 | namespace OpenMD {
698      //add bond constraints
699      Molecule::BondIterator bi;
700      Bond* bond;
701 +    ConstraintPair* cPair;
702 +
703      for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
704          
705        BondType* bt = bond->getBondType();
706  
492 –      //class Parent1 {};
493 –      //class Child1 : public Parent {};
494 –      //class Child2 : public Parent {};
495 –      //Child1* ch1 = new Child1();
496 –      //Child2* ch2 = dynamic_cast<Child2*>(ch1);
497 –      //the dynamic_cast is succeed in above line. A compiler bug?        
498 –
707        if (typeid(FixedBondType) == typeid(*bt)) {
708          FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
709  
710          ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
711          ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
712 <        ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
713 <        mol->addConstraintPair(consPair);
712 >        cPair = new ConstraintPair(consElemA, consElemB,
713 >                                   fbt->getEquilibriumBondLength(), false);
714 >        mol->addConstraintPair(cPair);
715        }
716      }
717  
# Line 514 | Line 723 | namespace OpenMD {
723      ConstraintPair* consPair;
724      Molecule::ConstraintPairIterator cpi;
725      std::set<StuntDouble*> sdSet;
726 <    for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
726 >    for (consPair = mol->beginConstraintPair(cpi); consPair != NULL;
727 >         consPair = mol->nextConstraintPair(cpi)) {
728  
729        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
730        if (sdSet.find(sdA) == sdSet.end()){
731          sdSet.insert(sdA);
732          mol->addConstraintElem(new ConstraintElem(sdA));
733        }
734 <
734 >      
735        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
736        if (sdSet.find(sdB) == sdSet.end()){
737          sdSet.insert(sdB);
738          mol->addConstraintElem(new ConstraintElem(sdB));
739 <      }
530 <        
739 >      }      
740      }
532 –
741    }
534 –    
742   }

Diff Legend

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