--- trunk/src/brains/MoleculeCreator.cpp	2005/03/08 21:06:49	403
+++ trunk/src/brains/MoleculeCreator.cpp	2005/12/02 15:38:03	770
@@ -1,4 +1,4 @@
- /*
+/*
  * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
  * The University of Notre Dame grants you ("Licensee") a
@@ -39,13 +39,13 @@
  * such damages.
  */
   
- /**
-  * @file MoleculeCreator.cpp
-  * @author tlin
-  * @date 11/04/2004
-  * @time 13:44am
-  * @version 1.0
-  */
+/**
+ * @file MoleculeCreator.cpp
+ * @author tlin
+ * @date 11/04/2004
+ * @time 13:44am
+ * @version 1.0
+ */
 
 #include <cassert>
 #include <set>
@@ -59,20 +59,20 @@ namespace oopse {
 #include "utils/StringUtils.hpp"
 
 namespace oopse {
+  
+  Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
+					    int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
 
-Molecule* MoleculeCreator::createMolecule(ForceField* ff, MoleculeStamp *molStamp,
-    int stampId, int globalIndex, LocalIndexManager* localIndexMan) {
-
-    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getID());
+    Molecule* mol = new Molecule(stampId, globalIndex, molStamp->getName());
     
     //create atoms
     Atom* atom;
     AtomStamp* currentAtomStamp;
     int nAtom = molStamp->getNAtoms();
     for (int i = 0; i < nAtom; ++i) {
-        currentAtomStamp = molStamp->getAtom(i);
-        atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
-        mol->addAtom(atom);
+      currentAtomStamp = molStamp->getAtomStamp(i);
+      atom = createAtom(ff, mol, currentAtomStamp, localIndexMan);
+      mol->addAtom(atom);
     }
 
     //create rigidbodies
@@ -81,9 +81,9 @@ Molecule* MoleculeCreator::createMolecule(ForceField* 
     int nRigidbodies = molStamp->getNRigidBodies();
 
     for (int i = 0; i < nRigidbodies; ++i) {
-        currentRigidBodyStamp = molStamp->getRigidBody(i);
-        rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
-        mol->addRigidBody(rb);
+      currentRigidBodyStamp = molStamp->getRigidBodyStamp(i);
+      rb = createRigidBody(molStamp, mol, currentRigidBodyStamp, localIndexMan);
+      mol->addRigidBody(rb);
     }
 
     //create bonds
@@ -92,9 +92,9 @@ Molecule* MoleculeCreator::createMolecule(ForceField* 
     int nBonds = molStamp->getNBonds();
 
     for (int i = 0; i < nBonds; ++i) {
-        currentBondStamp = molStamp->getBond(i);
-        bond = createBond(ff, mol, currentBondStamp);
-        mol->addBond(bond);
+      currentBondStamp = molStamp->getBondStamp(i);
+      bond = createBond(ff, mol, currentBondStamp);
+      mol->addBond(bond);
     }
 
     //create bends
@@ -102,9 +102,9 @@ Molecule* MoleculeCreator::createMolecule(ForceField* 
     BendStamp* currentBendStamp;
     int nBends = molStamp->getNBends();
     for (int i = 0; i < nBends; ++i) {
-        currentBendStamp = molStamp->getBend(i);
-        bend = createBend(ff, mol, currentBendStamp);
-        mol->addBend(bend);
+      currentBendStamp = molStamp->getBendStamp(i);
+      bend = createBend(ff, mol, currentBendStamp);
+      mol->addBend(bend);
     }
 
     //create torsions
@@ -112,9 +112,9 @@ Molecule* MoleculeCreator::createMolecule(ForceField* 
     TorsionStamp* currentTorsionStamp;
     int nTorsions = molStamp->getNTorsions();
     for (int i = 0; i < nTorsions; ++i) {
-        currentTorsionStamp = molStamp->getTorsion(i);
-        torsion = createTorsion(ff, mol, currentTorsionStamp);
-        mol->addTorsion(torsion);
+      currentTorsionStamp = molStamp->getTorsionStamp(i);
+      torsion = createTorsion(ff, mol, currentTorsionStamp);
+      mol->addTorsion(torsion);
     }
 
     //create cutoffGroups
@@ -122,55 +122,38 @@ Molecule* MoleculeCreator::createMolecule(ForceField* 
     CutoffGroupStamp* currentCutoffGroupStamp;
     int nCutoffGroups = molStamp->getNCutoffGroups();
     for (int i = 0; i < nCutoffGroups; ++i) {
-        currentCutoffGroupStamp = molStamp->getCutoffGroup(i);
-        cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
-        mol->addCutoffGroup(cutoffGroup);
+      currentCutoffGroupStamp = molStamp->getCutoffGroupStamp(i);
+      cutoffGroup = createCutoffGroup(mol, currentCutoffGroupStamp);
+      mol->addCutoffGroup(cutoffGroup);
     }
 
     //every free atom is a cutoff group    
-    std::set<Atom*> allAtoms;
-     Molecule::AtomIterator ai;
+    std::vector<Atom*> freeAtoms;
+    std::vector<Atom*>::iterator ai;
+    std::vector<Atom*>::iterator fai;
 
     //add all atoms into allAtoms set
-    for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
-        allAtoms.insert(atom);
+    for(atom = mol->beginAtom(fai); atom != NULL; atom = mol->nextAtom(fai)) {
+      freeAtoms.push_back(atom);
     }
 
     Molecule::CutoffGroupIterator ci;
     CutoffGroup* cg;
-    std::set<Atom*> cutoffAtoms;    
     
-    //add all of the atoms belong to cutoff groups into cutoffAtoms set
     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
 
-        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
-            cutoffAtoms.insert(atom);
-        }
+      for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
+           //erase the atoms belong to cutoff groups from freeAtoms vector
+           freeAtoms.erase(std::remove(freeAtoms.begin(), freeAtoms.end(), atom), freeAtoms.end());
+      }
 
     }       
     
-    //find all free atoms (which do not belong to cutoff groups)  
-    //performs the "difference" operation from set theory,  the output range contains a copy of every
-    //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in 
-    //[cutoffAtoms.begin(), cutoffAtoms.end()).
-    std::vector<Atom*> freeAtoms;    
-    std::set_difference(allAtoms.begin(), allAtoms.end(), cutoffAtoms.begin(), cutoffAtoms.end(),
-                            std::back_inserter(freeAtoms));
-
-    if (freeAtoms.size() != allAtoms.size() - cutoffAtoms.size()) {
-        //Some atoms in rigidAtoms are not in allAtoms, something must be wrong
-        sprintf(painCave.errMsg, "Atoms in cutoff groups are not in the atom list of the same molecule");
-
-        painCave.isFatal = 1;
-        simError();        
-    }
-
     //loop over the free atoms and then create one cutoff group for every single free atom
-    std::vector<Atom*>::iterator fai;
-
+    
     for (fai = freeAtoms.begin(); fai != freeAtoms.end(); ++fai) {
-        cutoffGroup = createCutoffGroup(mol, *fai);
-        mol->addCutoffGroup(cutoffGroup);
+      cutoffGroup = createCutoffGroup(mol, *fai);
+      mol->addCutoffGroup(cutoffGroup);
     }
     //create constraints
     createConstraintPair(mol);
@@ -180,52 +163,52 @@ Molecule* MoleculeCreator::createMolecule(ForceField* 
     mol->complete();
 
     return mol;
-}    
+  }    
 
 
-Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp, 
-                                                                  LocalIndexManager* localIndexMan) {
+  Atom* MoleculeCreator::createAtom(ForceField* ff, Molecule* mol, AtomStamp* stamp, 
+				    LocalIndexManager* localIndexMan) {
     AtomType * atomType;
     Atom* atom;
 
     atomType =  ff->getAtomType(stamp->getType());
 
     if (atomType == NULL) {
-        sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
-                   stamp->getType());
+      sprintf(painCave.errMsg, "Can not find Matching Atom Type for[%s]",
+	      stamp->getType().c_str());
 
-        painCave.isFatal = 1;
-        simError();
+      painCave.isFatal = 1;
+      simError();
     }
     
     //below code still have some kind of hard-coding smell
     if (atomType->isDirectional()){
      
-        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
+      DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
         
-        if (dAtomType == NULL) {
-            sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
+      if (dAtomType == NULL) {
+	sprintf(painCave.errMsg, "Can not cast AtomType to DirectionalAtomType");
 
-            painCave.isFatal = 1;
-            simError();
-        }
+	painCave.isFatal = 1;
+	simError();
+      }
 
-        DirectionalAtom* dAtom;
-        dAtom = new DirectionalAtom(dAtomType);
-        atom = dAtom;    
+      DirectionalAtom* dAtom;
+      dAtom = new DirectionalAtom(dAtomType);
+      atom = dAtom;    
     }
     else{
-        atom = new Atom(atomType);
+      atom = new Atom(atomType);
     }
 
     atom->setLocalIndex(localIndexMan->getNextAtomIndex());
 
     return atom;
-}
+  }
 
-RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
-                                                                                    RigidBodyStamp* rbStamp, 
-                                                                                    LocalIndexManager* localIndexMan) {
+  RigidBody* MoleculeCreator::createRigidBody(MoleculeStamp *molStamp, Molecule* mol,
+					      RigidBodyStamp* rbStamp, 
+					      LocalIndexManager* localIndexMan) {
     Atom* atom;
     int nAtoms;
     Vector3d refCoor;
@@ -234,11 +217,11 @@ RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
     RigidBody* rb = new RigidBody();
     nAtoms = rbStamp->getNMembers();    
     for (int i = 0; i < nAtoms; ++i) {
-        //rbStamp->getMember(i) return the local index of current atom inside the molecule.
-        //It is not the same as local index of atom which is the index of atom at DataStorage class
-        atom = mol->getAtomAt(rbStamp->getMember(i));
-        atomStamp= molStamp->getAtom(rbStamp->getMember(i));    
-        rb->addAtom(atom, atomStamp);
+      //rbStamp->getMember(i) return the local index of current atom inside the molecule.
+      //It is not the same as local index of atom which is the index of atom at DataStorage class
+      atom = mol->getAtomAt(rbStamp->getMemberAt(i));
+      atomStamp= molStamp->getAtomStamp(rbStamp->getMemberAt(i));    
+      rb->addAtom(atom, atomStamp);
     }
 
     //after all of the atoms are added, we need to calculate the reference coordinates
@@ -257,9 +240,9 @@ RigidBody* MoleculeCreator::createRigidBody(MoleculeSt
     rb->setType(mol->getType() + "_RB_" + s.c_str());
 
     return rb;
-}    
+  }    
 
-Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
+  Bond* MoleculeCreator::createBond(ForceField* ff, Molecule* mol, BondStamp* stamp) {
     BondType* bondType;
     Atom* atomA;
     Atom* atomB;
@@ -272,169 +255,132 @@ Bond* MoleculeCreator::createBond(ForceField* ff, Mole
     bondType = ff->getBondType(atomA->getType(), atomB->getType());
 
     if (bondType == NULL) {
-        sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
-                   atomA->getType().c_str(),
-                   atomB->getType().c_str());
+      sprintf(painCave.errMsg, "Can not find Matching Bond Type for[%s, %s]",
+	      atomA->getType().c_str(),
+	      atomB->getType().c_str());
 
-        painCave.isFatal = 1;
-        simError();
+      painCave.isFatal = 1;
+      simError();
     }
     return new Bond(atomA, atomB, bondType);    
-}    
+  }    
 
-Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
-    bool isGhostBend = false;
-    int ghostIndex;
+  Bend* MoleculeCreator::createBend(ForceField* ff, Molecule* mol, BendStamp* stamp) {
+    Bend* bend = NULL; 
+    std::vector<int> bendAtoms = stamp->getMembers(); 
+    if (bendAtoms.size() == 3) {
+      Atom* atomA = mol->getAtomAt(bendAtoms[0]);
+      Atom* atomB = mol->getAtomAt(bendAtoms[1]);
+      Atom* atomC = mol->getAtomAt(bendAtoms[2]);
 
-    
-    //
-    if (stamp->haveExtras()){
-        LinkedAssign* extras = stamp->getExtras();
-        LinkedAssign* currentExtra = extras;
+      assert( atomA && atomB && atomC);
 
-        while (currentExtra != NULL){
-            if (!strcmp(currentExtra->getlhs(), "ghostVectorSource")){
-                switch (currentExtra->getType()){
-                case 0:
-                    ghostIndex = currentExtra->getInt();
-                    isGhostBend = true;
-                    break;
+      BendType* bendType = ff->getBendType(atomA->getType().c_str(), atomB->getType().c_str(), atomC->getType().c_str());
 
-                default:
-                sprintf(painCave.errMsg,
-                "SimSetup Error: ghostVectorSource must be an int.\n");
-                painCave.isFatal = 1;
-                simError();
-                }
-            } else{
-                sprintf(painCave.errMsg,
-                "SimSetup Error: unhandled bend assignment:\n");
-                painCave.isFatal = 1;
-                simError();
-            }
-            currentExtra = currentExtra->getNext();
-        }
-        
-    }
+      if (bendType == NULL) {
+        sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
+                atomA->getType().c_str(),
+                atomB->getType().c_str(),
+                atomC->getType().c_str());
 
-    if (isGhostBend) {
+        painCave.isFatal = 1;
+        simError();
+      }
 
-        int indexA = stamp->getA();
-        int indexB= stamp->getB();
+      bend = new Bend(atomA, atomB, atomC, bendType);
+    } else if ( bendAtoms.size() == 2 && stamp->haveGhostVectorSource()) {
+      int ghostIndex = stamp->getGhostVectorSource();
+      int normalIndex = ghostIndex != bendAtoms[0] ? bendAtoms[0] : bendAtoms[1]; 
+      Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
+      DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
+      if (ghostAtom == NULL) {
+	sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
+	painCave.isFatal = 1;
+	simError();
+      }
+                
+      BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
 
-        assert(indexA != indexB);
+      if (bendType == NULL) {
+	sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
+		normalAtom->getType().c_str(),
+		ghostAtom->getType().c_str(),
+		"GHOST");
 
-        int normalIndex;
-        if (indexA == ghostIndex) {
-            normalIndex = indexB;
-        } else if (indexB == ghostIndex) {
-            normalIndex = indexA;
-        }
+	painCave.isFatal = 1;
+	simError();
+      }
         
-        Atom* normalAtom = mol->getAtomAt(normalIndex) ;        
-        DirectionalAtom* ghostAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(ghostIndex));
-        if (ghostAtom == NULL) {
-            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
-            painCave.isFatal = 1;
-            simError();
-        }
-                
-        BendType* bendType = ff->getBendType(normalAtom->getType(), ghostAtom->getType(), "GHOST");
+      bend = new GhostBend(normalAtom, ghostAtom, bendType);       
 
-        if (bendType == NULL) {
-            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
-                       normalAtom->getType().c_str(),
-                       ghostAtom->getType().c_str(),
-                       "GHOST");
+    } 
+    
+    return bend;
+  }    
 
-            painCave.isFatal = 1;
-            simError();
-        }
-        
-        return new GhostBend(normalAtom, ghostAtom, bendType);       
+  Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
 
-    } else {
-            
-        Atom* atomA = mol->getAtomAt(stamp->getA());
-        Atom* atomB = mol->getAtomAt(stamp->getB());
-        Atom* atomC = mol->getAtomAt(stamp->getC());
-
-        assert( atomA && atomB && atomC);
-        
-        BendType* bendType = ff->getBendType(atomA->getType(), atomB->getType(), atomC->getType());
-
-        if (bendType == NULL) {
-            sprintf(painCave.errMsg, "Can not find Matching Bend Type for[%s, %s, %s]",
-                       atomA->getType().c_str(),
-                       atomB->getType().c_str(),
-                       atomC->getType().c_str());
-
-            painCave.isFatal = 1;
-            simError();
-        }
-
-        return new Bend(atomA, atomB, atomC, bendType);       
+    Torsion* torsion = NULL;
+    std::vector<int> torsionAtoms = stamp->getMembers();
+    if (torsionAtoms.size() < 3) {
+	return torsion;
     }
-}    
 
-Torsion* MoleculeCreator::createTorsion(ForceField* ff, Molecule* mol, TorsionStamp* stamp) {
+    Atom* atomA = mol->getAtomAt(torsionAtoms[0]);
+    Atom* atomB = mol->getAtomAt(torsionAtoms[1]);
+    Atom* atomC = mol->getAtomAt(torsionAtoms[2]);
 
-    Atom* atomA = mol->getAtomAt(stamp->getA());
-    Atom* atomB = mol->getAtomAt(stamp->getB());
-    Atom* atomC = mol->getAtomAt(stamp->getC());
-    Torsion* torsion;
+    if (torsionAtoms.size() == 4) {
+      Atom* atomD = mol->getAtomAt(torsionAtoms[3]);
 
-    if (stamp->getD() != -1) {
-        Atom* atomD = mol->getAtomAt(stamp->getD());
-
-        assert(atomA && atomB && atomC && atomD);
+      assert(atomA && atomB && atomC && atomD);
         
-        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), 
-                                                           atomC->getType(), atomD->getType());
+      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), 
+						    atomC->getType(), atomD->getType());
 
-        if (torsionType == NULL) {
-            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
-                       atomA->getType().c_str(),
-                       atomB->getType().c_str(),
-                       atomC->getType().c_str(),
-                       atomD->getType().c_str());
+      if (torsionType == NULL) {
+	sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
+		atomA->getType().c_str(),
+		atomB->getType().c_str(),
+		atomC->getType().c_str(),
+		atomD->getType().c_str());
 
-            painCave.isFatal = 1;
-            simError();
-        }
+	painCave.isFatal = 1;
+	simError();
+      }
         
-        torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);       
+      torsion = new Torsion(atomA, atomB, atomC, atomD, torsionType);       
     }
     else {
 
-        DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(atomC);
-        if (dAtom == NULL) {
-            sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
-            painCave.isFatal = 1;
-            simError();
-        }        
+      DirectionalAtom* dAtom = dynamic_cast<DirectionalAtom*>(mol->getAtomAt(stamp->getGhostVectorSource()));
+      if (dAtom == NULL) {
+	sprintf(painCave.errMsg, "Can not cast Atom to DirectionalAtom");
+	painCave.isFatal = 1;
+	simError();
+      }        
 
-        TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), 
-                                                           atomC->getType(), "GHOST");
+      TorsionType* torsionType = ff->getTorsionType(atomA->getType(), atomB->getType(), 
+						    atomC->getType(), "GHOST");
 
-        if (torsionType == NULL) {
-            sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
-                       atomA->getType().c_str(),
-                       atomB->getType().c_str(),
-                       atomC->getType().c_str(),
-                       "GHOST");
+      if (torsionType == NULL) {
+	sprintf(painCave.errMsg, "Can not find Matching Torsion Type for[%s, %s, %s, %s]",
+		atomA->getType().c_str(),
+		atomB->getType().c_str(),
+		atomC->getType().c_str(),
+		"GHOST");
 
-            painCave.isFatal = 1;
-            simError();
-        }
+	painCave.isFatal = 1;
+	simError();
+      }
         
-        torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);               
+      torsion = new GhostTorsion(atomA, atomB, dAtom, torsionType);               
     }
 
     return torsion;
-}    
+  }    
 
-CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
+  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule* mol, CutoffGroupStamp* stamp) {
     int nAtoms;
     CutoffGroup* cg;
     Atom* atom;
@@ -442,71 +388,71 @@ CutoffGroup* MoleculeCreator::createCutoffGroup(Molecu
     
     nAtoms = stamp->getNMembers();
     for (int i =0; i < nAtoms; ++i) {
-        atom = mol->getAtomAt(stamp->getMember(i));
-        assert(atom);
-        cg->addAtom(atom);
+      atom = mol->getAtomAt(stamp->getMemberAt(i));
+      assert(atom);
+      cg->addAtom(atom);
     }
 
     return cg;
-}    
+  }    
 
-CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
+  CutoffGroup* MoleculeCreator::createCutoffGroup(Molecule * mol, Atom* atom) {
     CutoffGroup* cg;
     cg  = new CutoffGroup();
     cg->addAtom(atom);
     return cg;
-}
+  }
 
-void MoleculeCreator::createConstraintPair(Molecule* mol) {
+  void MoleculeCreator::createConstraintPair(Molecule* mol) {
 
     //add bond constraints
     Molecule::BondIterator bi;
     Bond* bond;
     for (bond = mol->beginBond(bi); bond != NULL; bond = mol->nextBond(bi)) {
         
-        BondType* bt = bond->getBondType();
+      BondType* bt = bond->getBondType();
 
-        //class Parent1 {};
-        //class Child1 : public Parent {};
-        //class Child2 : public Parent {};
-        //Child1* ch1 = new Child1();
-        //Child2* ch2 = dynamic_cast<Child2*>(ch1); 
-        //the dynamic_cast is succeed in above line. A compiler bug?        
+      //class Parent1 {};
+      //class Child1 : public Parent {};
+      //class Child2 : public Parent {};
+      //Child1* ch1 = new Child1();
+      //Child2* ch2 = dynamic_cast<Child2*>(ch1); 
+      //the dynamic_cast is succeed in above line. A compiler bug?        
 
-        if (typeid(FixedBondType) == typeid(*bt)) {
-            FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
+      if (typeid(FixedBondType) == typeid(*bt)) {
+	FixedBondType* fbt = dynamic_cast<FixedBondType*>(bt);
 
-            ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
-            ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
-            ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
-            mol->addConstraintPair(consPair);
-        }
+	ConstraintElem* consElemA = new ConstraintElem(bond->getAtomA());
+	ConstraintElem* consElemB = new ConstraintElem(bond->getAtomB());            
+	ConstraintPair* consPair = new ConstraintPair(consElemA, consElemB, fbt->getEquilibriumBondLength());
+	mol->addConstraintPair(consPair);
+      }
     }
 
     //rigidbody -- rigidbody constraint is not support yet
-}
+  }
 
-void MoleculeCreator::createConstraintElem(Molecule* mol) {
+  void MoleculeCreator::createConstraintElem(Molecule* mol) {
 
     ConstraintPair* consPair;
     Molecule::ConstraintPairIterator cpi;
     std::set<StuntDouble*> sdSet;
     for (consPair = mol->beginConstraintPair(cpi); consPair != NULL; consPair = mol->nextConstraintPair(cpi)) {
 
-        StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
-        if (sdSet.find(sdA) == sdSet.end()){
-            sdSet.insert(sdA);
-            mol->addConstraintElem(new ConstraintElem(sdA));
-        }
+      StuntDouble* sdA = consPair->getConsElem1()->getStuntDouble();            
+      if (sdSet.find(sdA) == sdSet.end()){
+	sdSet.insert(sdA);
+	mol->addConstraintElem(new ConstraintElem(sdA));
+      }
 
-        StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
-        if (sdSet.find(sdB) == sdSet.end()){
-            sdSet.insert(sdB);
-            mol->addConstraintElem(new ConstraintElem(sdB));
-        }
+      StuntDouble* sdB = consPair->getConsElem2()->getStuntDouble();            
+      if (sdSet.find(sdB) == sdSet.end()){
+	sdSet.insert(sdB);
+	mol->addConstraintElem(new ConstraintElem(sdB));
+      }
         
     }
 
-}
+  }
     
 }