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

Diff Legend

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