diff --git a/modules/mol/mm/examples/ethanol_example_using_topology.py b/modules/mol/mm/examples/ethanol_example_using_topology.py
index 8cca033a9eba699b4558c013f20822d97aef721a..41104e2a2ecbc72e177eb4fa6706e626fb3ec287 100644
--- a/modules/mol/mm/examples/ethanol_example_using_topology.py
+++ b/modules/mol/mm/examples/ethanol_example_using_topology.py
@@ -39,7 +39,7 @@ prof = io.IOProfile(dialect='PDB', fault_tolerant=False,
 
 ent = io.LoadPDB('ethanol.pdb', profile=prof)
 masses = [12.011,12.011,15.999,1.008,1.008,1.008,1.008,1.008,1.008]
-top = Topology(ent,masses)
+top = Topology(masses)
 
 #all per atom parameters have to be set at once...
 charges = [-0.18,0.145,-0.683,0.06,0.06,0.06,0.06,0.06,0.418]
@@ -144,7 +144,7 @@ settings.integrator = VerletIntegrator(0.0001)
 cc_bond_index = top.GetHarmonicBondIndices(0,1)[0]
 
 #create the simulation
-sim = Simulation(top,settings)
+sim = Simulation(top,ent,settings)
 go = gfx.Entity("ethanol",sim.GetEntity())
 go.SetRenderMode(gfx.CUSTOM)
 scene.Add(go)
diff --git a/modules/mol/mm/examples/gb_example.py b/modules/mol/mm/examples/gb_example.py
index 60a5ddf8934ecbbf29b01dc8d3623be88453d0c4..cad4e5e1dddf2d68e1fd7bbbbf8865ac5fbc95bc 100644
--- a/modules/mol/mm/examples/gb_example.py
+++ b/modules/mol/mm/examples/gb_example.py
@@ -40,7 +40,7 @@ settings.nonbonded_method = mm.NonbondedMethod.CutoffNonPeriodic
 
 sim = mm.Simulation(prot,settings)
 sim.MinimizeEnergy(type = "lbfgs",tolerance = 1.0, max_iterations = 200)
-sim.UpdateTopologyPositions()
+sim.UpdatePositions()
 
 ent = sim.GetEntity()
 go = gfx.Entity("crambin",ent)
diff --git a/modules/mol/mm/examples/gb_example_traj.dcd b/modules/mol/mm/examples/gb_example_traj.dcd
index 913d8d33f5f4f0a55f44002966ec332007c7dffe..f2dda449862fad41d41373d4ec76e52b72617d76 100644
Binary files a/modules/mol/mm/examples/gb_example_traj.dcd and b/modules/mol/mm/examples/gb_example_traj.dcd differ
diff --git a/modules/mol/mm/examples/gb_example_view_trajectory.py b/modules/mol/mm/examples/gb_example_view_trajectory.py
index 2075c4c4535593ae2d2c21f033a9a194eed0c36d..ba544247935af808cff70e3869a0cdd4492e6c17 100644
--- a/modules/mol/mm/examples/gb_example_view_trajectory.py
+++ b/modules/mol/mm/examples/gb_example_view_trajectory.py
@@ -1,5 +1,7 @@
-from ost.gui import traj
+from ost.gui import trajectory_viewer
 
 t = io.LoadCHARMMTraj("gb_example_traj.pdb","gb_example_traj.dcd")
-w = traj.TrajWidget(t)
+go = gfx.Entity("gb",t.GetEntity())
+scene.Add(go)
+w = trajectory_viewer.TrajWidget([t],[go])
 w.show()
diff --git a/modules/mol/mm/pymod/export_simulation.cc b/modules/mol/mm/pymod/export_simulation.cc
index 91f5e25b95362210e39788654b4cf7dee2fbbc6d..1b358eed00b9a2517cc115f197bb5f14a4a97b44 100644
--- a/modules/mol/mm/pymod/export_simulation.cc
+++ b/modules/mol/mm/pymod/export_simulation.cc
@@ -48,8 +48,8 @@ void export_Simulation()
 {
 
   class_<ost::mol::mm::Simulation>("Simulation",no_init)
-    .def(init<ost::mol::EntityHandle, ost::mol::mm::MMSettingsPtr>())
-    .def(init<ost::mol::mm::TopologyPtr,ost::mol::mm::MMSettingsPtr>())
+    .def(init<const ost::mol::EntityHandle, const ost::mol::mm::MMSettingsPtr>())
+    .def(init<const ost::mol::mm::TopologyPtr,const ost::mol::EntityHandle&,const ost::mol::mm::MMSettingsPtr>())
     .def("Save",&ost::mol::mm::Simulation::Save)
     .def("Load",&ost::mol::mm::Simulation::Load).staticmethod("Load")
     .def("Steps",&ost::mol::mm::Simulation::Steps)
@@ -64,7 +64,7 @@ void export_Simulation()
     .def("GetPotentialEnergy",&ost::mol::mm::Simulation::GetPotentialEnergy)                                                                             
     .def("GetKineticEnergy",&ost::mol::mm::Simulation::GetKineticEnergy) 
     .def("Register",&ost::mol::mm::Simulation::Register) 
-    .def("UpdateTopologyPositions",&ost::mol::mm::Simulation::UpdateTopologyPositions, (arg("enforce_periodic_box")=false))
+    .def("UpdatePositions",&ost::mol::mm::Simulation::UpdatePositions, (arg("enforce_periodic_box")=false))
     .def("GetTopology",&ost::mol::mm::Simulation::GetTopology)
     .def("ResetHarmonicBond",&ost::mol::mm::Simulation::ResetHarmonicBond) 
     .def("ResetHarmonicAngle",&ost::mol::mm::Simulation::ResetHarmonicAngle)
diff --git a/modules/mol/mm/pymod/export_topology.cc b/modules/mol/mm/pymod/export_topology.cc
index 1307ad6816cd02437c0c232c957ecb37643e4f37..eb7a16fcb111afb55b9736a2fb420c564dba274d 100644
--- a/modules/mol/mm/pymod/export_topology.cc
+++ b/modules/mol/mm/pymod/export_topology.cc
@@ -36,10 +36,9 @@ namespace{
     return l;
   }
 
-  ost::mol::mm::TopologyPtr WrapTopologyConstructor(ost::mol::EntityHandle& handle, 
-                                     const boost::python::list& masses_list){
+  ost::mol::mm::TopologyPtr WrapTopologyConstructor(const boost::python::list& masses_list){
     std::vector<Real> masses_vector = ListToVec<Real>(masses_list);
-    ost::mol::mm::TopologyPtr p(new ost::mol::mm::Topology(handle,masses_vector));
+    ost::mol::mm::TopologyPtr p(new ost::mol::mm::Topology(masses_vector));
     return p;
   }
 
@@ -296,16 +295,12 @@ namespace{
   }
 
   boost::python::list GetHarmonicPositionRestraintIndicesSingleIndex(ost::mol::mm::TopologyPtr p, uint i){
-    std::vector<uint> return_vec = p->GetHarmonicPositionRestraintIndices(i);
+  std::vector<uint> return_vec = p->GetHarmonicPositionRestraintIndices(i);
     return VecToList<uint>(return_vec);
   }
 
 }
 
-long GetIndexA(ost::mol::mm::TopologyPtr top,ost::mol::AtomHandle& at) { return top->GetAtomIndex(at); }
-long GetIndexB(ost::mol::mm::TopologyPtr top,int r_idx,const String& aname) { return top->GetAtomIndex(r_idx,aname); }
-
-
 void export_Topology()
 {
 
@@ -318,7 +313,6 @@ void export_Topology()
     .def("__init__",make_constructor(&WrapTopologyConstructor))
     .def("Save",&ost::mol::mm::Topology::Save)
     .def("Load",&ost::mol::mm::Topology::Load).staticmethod("Load")
-    .def("GetEntity",&ost::mol::mm::Topology::GetEntity)
     //interaction adding functions
     .def("AddHarmonicBond",&ost::mol::mm::Topology::AddHarmonicBond)
     .def("AddHarmonicAngle",&ost::mol::mm::Topology::AddHarmonicAngle)
@@ -335,7 +329,6 @@ void export_Topology()
     .def("ResetExclusions",&ost::mol::mm::Topology::ResetExclusions)
     .def("AddHarmonicPositionRestraint",&ost::mol::mm::Topology::AddHarmonicPositionRestraint,(arg("index"),arg("ref_position"),arg("k"),arg("x_scale")=1.0,arg("y_scale")=1.0,arg("z_scale")=1.0))
     .def("AddHarmonicDistanceRestraint",&ost::mol::mm::Topology::AddHarmonicDistanceRestraint)
-    .def("Merge",&ost::mol::mm::Topology::Merge)
     //single atom parameter getter and setter functions
     .def("SetSigmas",&WrapSetSigmas)
     .def("SetSigma",&ost::mol::mm::Topology::SetSigma)
@@ -412,7 +405,7 @@ void export_Topology()
     .def("GetHarmonicPositionRestraintIndices",&GetHarmonicPositionRestraintIndicesSingleIndex)
 
     //functions to get amount of data in topology
-    .def("GetNumAtoms",&ost::mol::mm::Topology::GetNumAtoms)
+    .def("GetNumParticles",&ost::mol::mm::Topology::GetNumParticles)
     .def("GetNumHarmonicBonds",&ost::mol::mm::Topology::GetNumHarmonicBonds)
     .def("GetNumHarmonicAngles",&ost::mol::mm::Topology::GetNumHarmonicAngles)
     .def("GetNumPeriodicDihedrals",&ost::mol::mm::Topology::GetNumPeriodicDihedrals)
@@ -423,10 +416,9 @@ void export_Topology()
     .def("GetNumDistanceConstraints",&ost::mol::mm::Topology::GetNumDistanceConstraints)
     .def("GetNumExclusions",&ost::mol::mm::Topology::GetNumExclusions)
     .def("GetNumHarmonicPositionRestraints",&ost::mol::mm::Topology::GetNumHarmonicPositionRestraints)
-    .def("GetNumHarmonicDistanceRestraints",&ost::mol::mm::Topology::GetNumHarmonicDistanceRestraints)
-    //get internal indices
-    .def("GetAtomIndex",&GetIndexA)
-    .def("GetAtomIndex",&GetIndexB)                                                               
+    .def("GetNumHarmonicDistanceRestraints",&ost::mol::mm::Topology::GetNumHarmonicDistanceRestraints)  
+
+    .def("Merge",&ost::mol::mm::Topology::Merge)                                                           
   ;
 
   boost::python::register_ptr_to_python<ost::mol::mm::TopologyPtr>();
diff --git a/modules/mol/mm/src/mm_observer.cc b/modules/mol/mm/src/mm_observer.cc
index ed1f14b7a9396f3765ede452b2d26d7c37bec101..a8c4da9ee2448f8a7c4a2a380046bf3005c74589 100644
--- a/modules/mol/mm/src/mm_observer.cc
+++ b/modules/mol/mm/src/mm_observer.cc
@@ -3,7 +3,8 @@
 namespace ost { namespace mol{ namespace mm{
 
 void TrajObserver::Init(boost::shared_ptr<OpenMM::Context> c, 
-                        TopologyPtr top){
+                        TopologyPtr top,
+                        ost::mol::EntityHandle& ent){
 
   if(registered_){
     throw ost::Error("Can register observer to only one simulation!");
@@ -12,8 +13,6 @@ void TrajObserver::Init(boost::shared_ptr<OpenMM::Context> c,
   registered_ = true;
 
   context_ = c;
-  ost::mol::EntityHandle ent = top->GetEntity().Copy();
-  MMModeller::AssignPDBNaming(ent);
   c_group_ = ost::mol::CreateCoordGroup(ent.GetAtomList());
 }
 
@@ -26,7 +25,8 @@ void TrajObserver::Notify(){
 
 
 void TrajWriter::Init(boost::shared_ptr<OpenMM::Context> c, 
-                      TopologyPtr top){
+                      TopologyPtr top,
+                      ost::mol::EntityHandle& ent){
 
   if(registered_){
     throw ost::Error("Can register observer to only one simulation!");
@@ -37,13 +37,13 @@ void TrajWriter::Init(boost::shared_ptr<OpenMM::Context> c,
   context_ = c;
   ost::io::IOProfile profile("CHARMM",false,false,false,false,false);
   ost::io::PDBWriter writer(pdb_filename_, profile);
-  writer.Write(top->GetEntity().GetAtomList());
+  writer.Write(ent.GetAtomList());
 
   stream_.open(dcd_filename_.c_str(), std::ios::binary);
 
-  x.resize(top->GetNumAtoms());
-  y.resize(top->GetNumAtoms());
-  z.resize(top->GetNumAtoms());
+  x.resize(top->GetNumParticles());
+  y.resize(top->GetNumParticles());
+  z.resize(top->GetNumParticles());
 
 
   // size of first header block in bytes
@@ -101,7 +101,7 @@ void TrajWriter::Init(boost::shared_ptr<OpenMM::Context> c,
 
   // atom count block
   stream_.write(reinterpret_cast<char*>(&four), 4);
-  int32_t atom_count=top->GetNumAtoms();  
+  int32_t atom_count=top->GetNumParticles();  
   stream_.write(reinterpret_cast<char*>(&atom_count), 4);
   stream_.write(reinterpret_cast<char*>(&four), 4);
 }
diff --git a/modules/mol/mm/src/mm_observer.hh b/modules/mol/mm/src/mm_observer.hh
index 92b52971a4c5b60f3fef5d3e60620e8b9c067588..1e95f092ce08f16ded321c53ee182d124ad2fae0 100644
--- a/modules/mol/mm/src/mm_observer.hh
+++ b/modules/mol/mm/src/mm_observer.hh
@@ -31,7 +31,8 @@ public:
   virtual void Notify() = 0;
 
   virtual void Init(boost::shared_ptr<OpenMM::Context> c,
-                    TopologyPtr top) = 0;
+                    TopologyPtr top,
+                    ost::mol::EntityHandle& ent) = 0;
 
   virtual int Rythm() = 0;
 
@@ -45,7 +46,8 @@ public:
   TrajObserver(int rythm): rythm_(rythm), registered_(false) { }
 
   void Init(boost::shared_ptr<OpenMM::Context> c, 
-            TopologyPtr top);
+            TopologyPtr top,
+            ost::mol::EntityHandle& ent);
 
   void Notify();
 
@@ -75,7 +77,8 @@ public:
                                                                                  frames_(0) { }
 
   void Init(boost::shared_ptr<OpenMM::Context> c, 
-            TopologyPtr top);
+            TopologyPtr top,
+            ost::mol::EntityHandle& ent);
 
   void Notify();
 
diff --git a/modules/mol/mm/src/simulation.cc b/modules/mol/mm/src/simulation.cc
index fe15d7329a86b8210e387e1921f08416c2311756..26a03e933eec7cb352f213dae661813dbaf592b6 100644
--- a/modules/mol/mm/src/simulation.cc
+++ b/modules/mol/mm/src/simulation.cc
@@ -8,13 +8,19 @@ Simulation::Simulation(const ost::mol::EntityHandle& handle,
   //note, that ent_ will be "completed" inside this function!
   //(hydrogens and shit)
 
-  TopologyPtr top = TopologyCreator::Create(handle,settings);
+  ent_ = handle.Copy();
+  TopologyPtr top = TopologyCreator::Create(ent_,settings);
   this->Init(top, settings);
 }
 
 Simulation::Simulation(const TopologyPtr top,
+                       const ost::mol::EntityHandle& handle,
                        const MMSettingsPtr settings){
 
+  if(static_cast<uint>(handle.GetAtomCount()) != top->GetNumParticles()){
+    throw ost::Error("Number of atoms in entity must be consistent with number of particles in topology!");
+  }
+  ent_ = handle.Copy();
   this->Init(top, settings);
 }
 
@@ -29,6 +35,86 @@ void Simulation::Save(const String& filename){
     ds & (*i)[1];
     ds & (*i)[2];
   }
+
+  uint num_chains;
+  uint num_residues;
+  uint num_atoms;
+  uint num_bonded_atoms;
+  Real bfac;
+  Real occ;
+  bool is_hetatm;
+  String chain_name;
+  String res_name;
+  int resnum_num;
+  char resnum_code;
+  String atom_name;
+  String atom_element;
+  uint atom_index;
+
+  num_chains = ent_.GetChainCount();
+  ds & num_chains;
+  ost::mol::ChainHandleList chain_list = ent_.GetChainList();
+  for(ost::mol::ChainHandleList::iterator i = chain_list.begin();
+      i != chain_list.end(); ++i){
+    chain_name = i->GetName();
+    num_residues = i->GetResidueCount();
+    ds & chain_name;
+    ds & num_residues;
+    ost::mol::ResidueHandleList res_list = i->GetResidueList();
+    for(ost::mol::ResidueHandleList::iterator j = res_list.begin();
+        j != res_list.end(); ++j){
+      res_name = j->GetKey();
+      resnum_num = j->GetNumber().GetNum();
+      resnum_code = j->GetNumber().GetInsCode(); 
+      num_atoms = j->GetAtomCount();
+      ds & res_name;
+      ds & resnum_num;
+      ds & resnum_code;
+      ds & num_atoms;
+      ost::mol::AtomHandleList atom_list = j->GetAtomList();
+      for(ost::mol::AtomHandleList::iterator k = atom_list.begin();
+          k != atom_list.end(); ++k){
+        atom_name = k->GetName();
+        atom_element = k->GetElement();
+        geom::Vec3 pos = k->GetPos();
+        bfac = k->GetBFactor();
+        occ = k->GetOccupancy();
+        is_hetatm = k->IsHetAtom();
+        ds & atom_name;
+        ds & atom_element;
+        ds & pos[0];
+        ds & pos[1];
+        ds & pos[2];
+        ds & bfac;
+        ds & occ;
+        ds & is_hetatm;
+      }
+    }
+  }
+
+  ost::mol::AtomHandleList atom_list = ent_.GetAtomList();
+  ost::mol::AtomHandleList bonded_atoms;
+
+  std::map<long,int> atom_indices;
+  int actual_index = 0;
+  for(ost::mol::AtomHandleList::const_iterator i = atom_list.begin(), e = atom_list.end(); 
+      i != e; ++i){
+    atom_indices[i->GetHashCode()] = actual_index;
+    ++actual_index;
+  }
+
+  for(ost::mol::AtomHandleList::iterator i = atom_list.begin();
+      i != atom_list.end(); ++i){
+    bonded_atoms = i->GetBondPartners();
+    num_bonded_atoms = bonded_atoms.size();
+    ds & num_bonded_atoms;
+    for(ost::mol::AtomHandleList::iterator j = bonded_atoms.begin();
+        j != bonded_atoms.end(); ++j){
+      atom_index = atom_indices[j->GetHashCode()];
+      ds & atom_index;
+    }
+  }
+
   context_->createCheckpoint(stream);  
 }
 
@@ -94,6 +180,63 @@ SimulationPtr Simulation::Load(const String& filename, MMSettingsPtr settings){
   }
   sim_ptr->context_->setPositions(positions);
 
+  uint num_chains;
+  uint num_residues;
+  uint num_atoms;
+  uint num_bonded_atoms;
+  Real x_pos;
+  Real y_pos;
+  Real z_pos;
+  Real bfac;
+  Real occ;
+  bool is_hetatm;
+  String chain_name;
+  String res_name;
+  int resnum_num;
+  char resnum_code;
+  String atom_name;
+  String atom_element;
+  uint atom_index;
+
+  ost::mol::EntityHandle ent = ost::mol::CreateEntity();
+  ost::mol::XCSEditor ed = ent.EditXCS();
+  ds & num_chains;
+  for(uint i = 0; i < num_chains; ++i){
+    ds & chain_name;
+    ds & num_residues;
+    ost::mol::ChainHandle chain = ed.InsertChain(chain_name);
+    for(uint j = 0; j < num_residues; ++j){
+      ds & res_name;
+      ds & resnum_num;
+      ds & resnum_code;
+      ds & num_atoms;
+      ost::mol::ResNum num(resnum_num,resnum_code);
+      ost::mol::ResidueHandle res = ed.AppendResidue(chain,res_name,num);
+      for(uint k = 0; k < num_atoms; ++k){
+        ds & atom_name;
+        ds & atom_element;
+        ds & x_pos;
+        ds & y_pos;
+        ds & z_pos;
+        ds & bfac;
+        ds & occ;
+        ds & is_hetatm;
+        geom::Vec3 pos(x_pos,y_pos,z_pos);
+        ed.InsertAtom(res,atom_name,pos,atom_element,occ,bfac,is_hetatm);
+      }
+    }
+  }
+  ost::mol::AtomHandleList atom_list = ent.GetAtomList();
+  for(uint i = 0; i < atom_list.size(); ++i){
+    ds & num_bonded_atoms;
+    for(uint j = 0; j < num_bonded_atoms; ++j){
+      ds & atom_index;
+      ed.Connect(atom_list[i],atom_list[atom_index]);
+    }
+  }
+
+  sim_ptr->ent_ = ent;
+
   sim_ptr->context_->loadCheckpoint(stream);
 
   return sim_ptr;
@@ -111,7 +254,6 @@ void Simulation::Init(const TopologyPtr top,
     throw ost::Error("Settings must have a valid integrator attached to set up a simulation!");
   }
   system_ = SystemCreator::Create(top_,settings,system_force_mapper_); 
-  ost::mol::EntityHandle ent = top_->GetEntity();
 
   //setting up the context, which combines the system with an integrator
   //to proceed in time, but first we have to load the proper platform
@@ -140,7 +282,7 @@ void Simulation::Init(const TopologyPtr top,
 
   context_ = ContextPtr(new OpenMM::Context(*system_,*integrator_,*platform));
 
-  ost::mol::AtomHandleList atom_list = ent.GetAtomList();
+  ost::mol::AtomHandleList atom_list = ent_.GetAtomList();
   std::vector<OpenMM::Vec3> positions;
   geom::Vec3 ost_vec;
   OpenMM::Vec3 open_mm_vec;
@@ -152,6 +294,7 @@ void Simulation::Init(const TopologyPtr top,
     open_mm_vec[2] = ost_vec[2]/10;
     positions.push_back(open_mm_vec);
   }
+
   context_->setPositions(positions);
 
   //make sure the context satisfies the distance constraints
@@ -181,51 +324,61 @@ geom::Vec3List Simulation::GetForces(){
 }
 
 void Simulation::SetPositions(const geom::Vec3List& positions, bool in_angstrom){
-  if(top_->GetNumAtoms() != positions.size()){
-    throw ost::Error("Number of positions does not correspond to number of atoms in topology!");
+  if(top_->GetNumParticles() != positions.size()){
+    throw ost::Error("Number of positions does not correspond to number of particles in topology!");
   }
-  std::vector<OpenMM::Vec3> openmm_positions;
-  OpenMM::Vec3 open_mm_vec;
+  std::vector<OpenMM::Vec3> openmm_positions(top_->GetNumParticles());
+  OpenMM::Vec3* actual_pos = &openmm_positions[0];
   if(in_angstrom){
     for(geom::Vec3List::const_iterator i = positions.begin();
         i != positions.end(); ++i){
-      open_mm_vec[0] = (*i)[0]*0.1;
-      open_mm_vec[1] = (*i)[1]*0.1;
-      open_mm_vec[2] = (*i)[2]*0.1;
-      openmm_positions.push_back(open_mm_vec);
+      (*actual_pos)[0] = (*i)[0]*0.1;
+      (*actual_pos)[1] = (*i)[1]*0.1;
+      (*actual_pos)[2] = (*i)[2]*0.1;
+      ++actual_pos;
     }
   }
   else{
     for(geom::Vec3List::const_iterator i = positions.begin();
         i != positions.end(); ++i){
-      open_mm_vec[0] = (*i)[0];
-      open_mm_vec[1] = (*i)[1];
-      open_mm_vec[2] = (*i)[2];
-      openmm_positions.push_back(open_mm_vec);
+      (*actual_pos)[0] = (*i)[0];
+      (*actual_pos)[1] = (*i)[1];
+      (*actual_pos)[2] = (*i)[2];
+      ++actual_pos;
     }
   }
   context_->setPositions(openmm_positions);
 }
 
 void Simulation::SetVelocities(geom::Vec3List& velocities){
-  if(top_->GetNumAtoms() != velocities.size()){
-    throw ost::Error("Number of velocities does not correspond to number of atoms in topology!");
+  if(top_->GetNumParticles() != velocities.size()){
+    throw ost::Error("Number of velocities does not correspond to number of particles in topology!");
   }
-  std::vector<OpenMM::Vec3> openmm_velocities;
-  OpenMM::Vec3 open_mm_vec;
+  std::vector<OpenMM::Vec3> openmm_velocities(top_->GetNumParticles());
+  OpenMM::Vec3* actual_vel = &openmm_velocities[0];
   for(geom::Vec3List::iterator i = velocities.begin();
       i != velocities.end(); ++i){
-    open_mm_vec[0] = (*i)[0];
-    open_mm_vec[1] = (*i)[1];
-    open_mm_vec[2] = (*i)[2];
-    openmm_velocities.push_back(open_mm_vec);
+    (*actual_vel)[0] = (*i)[0];
+    (*actual_vel)[1] = (*i)[1];
+    (*actual_vel)[2] = (*i)[2];
+    ++actual_vel;
   }
   context_->setVelocities(openmm_velocities);
 }
 
-void Simulation::UpdateTopologyPositions(bool enforce_periodic_box){
-  geom::Vec3List positions = this->GetPositions(enforce_periodic_box);
-  top_->SetEntityPositions(positions);
+void Simulation::UpdatePositions(bool enforce_periodic_box){
+  if(top_->GetNumParticles() != static_cast<uint>(ent_.GetAtomCount())){
+    throw ost::Error("Num particles in topology and num atoms in entity are not consistent!");
+  }
+  geom::Vec3List positions = this->GetPositions(enforce_periodic_box, true);
+  ost::mol::XCSEditor ed = ent_.EditXCS();
+  ost::mol::AtomHandleList atom_list = ent_.GetAtomList();
+  ost::mol::AtomHandleList::iterator a = atom_list.begin();
+  ost::mol::AtomHandleList::iterator ae = atom_list.end();
+  geom::Vec3List::iterator v = positions.begin();
+  for(; a != ae; ++a, ++v){
+    ed.SetAtomPos(*a,*v);
+  }
 }
 
 void Simulation::MinimizeEnergy(const String& type, Real tolerance, int max_iterations){
@@ -291,7 +444,7 @@ Real Simulation::GetPotentialEnergy(){
 void Simulation::Register(MMObserverPtr o){
   observers_.push_back(o);
   time_to_notify_.push_back(o->Rythm());
-  o->Init(context_,top_);
+  o->Init(context_,top_,ent_);
 }
 
 int Simulation::TimeToNextNotification(){
@@ -472,7 +625,7 @@ void Simulation::ResetDistanceConstraint(uint index, Real constraint_length){
 }
 
 void Simulation::AddPositionConstraint(uint index){
-  if(index >= top_->GetNumAtoms()){
+  if(index >= top_->GetNumParticles()){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   system_->setParticleMass(index,0.0);
@@ -484,7 +637,7 @@ void Simulation::AddPositionConstraints(const std::vector<uint>& index){
 
   for(std::vector<uint>::const_iterator i = index.begin(); 
       i != index.end(); ++i){
-    if(*i >= top_->GetNumAtoms()){
+    if(*i >= top_->GetNumParticles()){
       throw ost::Error("Provided index exceeds number of atoms!");
     }
     system_->setParticleMass(*i,0.0);
@@ -553,7 +706,7 @@ void Simulation::ResetHarmonicDistanceRestraint(uint index, Real length, Real fo
 }
 
 void Simulation::ResetLJ(uint index, Real sigma, Real epsilon){
-  if(index >= top_->GetNumAtoms()){
+  if(index >= top_->GetNumParticles()){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   if(system_force_mapper_.find(LJ) == system_force_mapper_.end()){
@@ -573,7 +726,7 @@ void Simulation::ResetLJ(uint index, Real sigma, Real epsilon){
 }
 
 void Simulation::ResetGBSA(uint index, Real radius, Real scaling){
-  if(index >= top_->GetNumAtoms()){
+  if(index >= top_->GetNumParticles()){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   if(system_force_mapper_.find(GBSA) == system_force_mapper_.end()){
@@ -593,7 +746,7 @@ void Simulation::ResetGBSA(uint index, Real radius, Real scaling){
 }
 
 void Simulation::ResetCharge(uint index, Real charge){
-  if(index >= top_->GetNumAtoms()){
+  if(index >= top_->GetNumParticles()){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   if(system_force_mapper_.find(LJ) == system_force_mapper_.end()){
@@ -641,7 +794,7 @@ void Simulation::ResetCharge(uint index, Real charge){
 }
 
 void Simulation::ResetMass(uint index, Real mass){
-  if(index >= top_->GetNumAtoms()){
+  if(index >= top_->GetNumParticles()){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   system_->setParticleMass(index,mass);
diff --git a/modules/mol/mm/src/simulation.hh b/modules/mol/mm/src/simulation.hh
index 3b918286830fbc5008b91ee0f65fa2f19484689f..14c835d95f8a0dd755cc7d17994e5e08a897d0fc 100644
--- a/modules/mol/mm/src/simulation.hh
+++ b/modules/mol/mm/src/simulation.hh
@@ -45,13 +45,14 @@ public:
              const MMSettingsPtr settings);
 
   Simulation(const ost::mol::mm::TopologyPtr top,
+             const ost::mol::EntityHandle& handle,
              const MMSettingsPtr settings);
 
   void Save(const String& filename);
 
   static SimulationPtr Load(const String& filename, MMSettingsPtr settings);
 
-  ost::mol::EntityHandle GetEntity() { return top_->GetEntity(); }
+  ost::mol::EntityHandle GetEntity() { return ent_; }
 
   geom::Vec3List GetPositions(bool enforce_periodic_box = false, bool in_angstrom = true);
 
@@ -63,7 +64,7 @@ public:
 
   void SetVelocities(geom::Vec3List& velocities);
 
-  void UpdateTopologyPositions(bool enforce_periodic_box = false);
+  void UpdatePositions(bool enforce_periodic_box = false);
 
   void MinimizeEnergy(const String& type = "steep", Real tolerance = 1.0, int max_iterations = 1000);
 
@@ -134,6 +135,7 @@ private:
   std::vector<MMObserverPtr> observers_;
   std::vector<int> time_to_notify_;
   std::map<FuncType,uint> system_force_mapper_;
+  ost::mol::EntityHandle ent_;
 };
 
 }}} //ns
diff --git a/modules/mol/mm/src/topology.cc b/modules/mol/mm/src/topology.cc
index 2ed6edad2866c580ab0246aba779438389105ddb..d5b238b842a710fb011b56a3c3599bf23c8a0448 100644
--- a/modules/mol/mm/src/topology.cc
+++ b/modules/mol/mm/src/topology.cc
@@ -6,21 +6,11 @@ namespace ost{ namespace mol{ namespace mm{
 
 
 
-Topology::Topology(const ost::mol::EntityHandle& ent, const std::vector<Real>& masses){
-
-  num_atoms_ = ent.GetAtomCount();
-  if(num_atoms_ != masses.size()){
-    throw ost::Error("Number of atoms in entity and number of masses do not match when creating topology!");
-  }
-  ent_ = ent.Copy();
-  num_residues_ = ent_.GetResidueCount();
-  atom_list_ = ent_.GetAtomList();
-  res_list_ = ent_.GetResidueList();
+Topology::Topology(const std::vector<Real>& masses){
+  num_particles_ = masses.size();
   atom_masses_ = masses;
   fudge_lj_ = 1.0;
   fudge_qq_ = 1.0;
-
-  this->InitMappers();
 }
 
 TopologyPtr Topology::Load(const String& filename){
@@ -50,8 +40,8 @@ uint Topology::AddHarmonicBond(uint index_one,
                                Real bond_length,
                                Real force_constant){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_){
-    throw ost::Error("Index of bond atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_){
+    throw ost::Error("Index of bond particle exceeds number of particles present in topology!");
   }
   Index<2> index(index_one,index_two);
   std::vector<Real> parameters;
@@ -67,9 +57,9 @@ uint Topology::AddHarmonicAngle(uint index_one,
                                 Real angle,
                                 Real force_constant){
 
-  if(uint(index_one) >= num_atoms_ || uint(index_two) >= num_atoms_ || 
-     uint(index_three) >= num_atoms_){
-    throw ost::Error("Index of angle atom exceeds number of atoms present in topology!");
+  if(uint(index_one) >= num_particles_ || uint(index_two) >= num_particles_ || 
+     uint(index_three) >= num_particles_){
+    throw ost::Error("Index of angle particle exceeds number of particles present in topology!");
   }
   Index<3> index(index_one,index_two,index_three); 
   std::vector<Real> parameters;
@@ -87,9 +77,9 @@ uint Topology::AddUreyBradleyAngle(uint index_one,
                                    Real bond_length,
                                    Real bond_force_constant){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_ || 
-     index_three >= num_atoms_){
-    throw ost::Error("Index of angle atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_ || 
+     index_three >= num_particles_){
+    throw ost::Error("Index of angle particle exceeds number of particles present in topology!");
   }
   Index<3> index(index_one,index_two,index_three); 
   std::vector<Real> parameters;
@@ -109,9 +99,9 @@ uint Topology::AddPeriodicDihedral(uint index_one,
                                    Real phase,
                                    Real force_constant){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_||
-     index_three >= num_atoms_ || index_four >= num_atoms_){
-    throw ost::Error("Index of dihedral atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_||
+     index_three >= num_particles_ || index_four >= num_particles_){
+    throw ost::Error("Index of dihedral particle exceeds number of particles present in topology!");
   }
   Index<4> index(index_one,index_two,index_three,index_four);
   std::vector<Real> parameters;
@@ -130,9 +120,9 @@ uint Topology::AddPeriodicImproper(uint index_one,
                                    Real phase,
                                    Real force_constant){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_||
-     index_three >= num_atoms_ || index_four >= num_atoms_){
-    throw ost::Error("Index of improper atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_||
+     index_three >= num_particles_ || index_four >= num_particles_){
+    throw ost::Error("Index of improper particle exceeds number of particles present in topology!");
   }
   Index<4> index(index_one,index_two,index_three,index_four);
   std::vector<Real> parameters;
@@ -150,9 +140,9 @@ uint Topology::AddHarmonicImproper(uint index_one,
                                    Real angle,
                                    Real force_constant){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_||
-     index_three >= num_atoms_ || index_four >= num_atoms_){
-    throw ost::Error("Index of improper atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_||
+     index_three >= num_particles_ || index_four >= num_particles_){
+    throw ost::Error("Index of improper particle exceeds number of particles present in topology!");
   }
   Index<4> index(index_one,index_two,index_three,index_four);
   std::vector<Real> parameters;
@@ -170,10 +160,10 @@ uint Topology::AddCMap(uint index_one,
                        int dimension,
                        std::vector<Real> values){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_ ||
-     index_three >= num_atoms_ || index_four >= num_atoms_ ||
-     index_five >= num_atoms_){
-    throw ost::Error("Index of cmap atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_ ||
+     index_three >= num_particles_ || index_four >= num_particles_ ||
+     index_five >= num_particles_){
+    throw ost::Error("Index of cmap particle exceeds number of particles present in topology!");
   }
   if(uint(dimension * dimension) != values.size()){
     throw ost::Error("Dimension of CMap is inconsistent with its values!");
@@ -188,8 +178,8 @@ uint Topology::AddLJPair(uint index_one,
                          Real sigma,
                          Real epsilon){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_){
-    throw ost::Error("Index of lj pair atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_){
+    throw ost::Error("Index of lj pair particle exceeds number of particles present in topology!");
   }
   Index<2> index(index_one,index_two);
   Index<2> reverse_index(index_two,index_one);
@@ -208,8 +198,8 @@ uint Topology::AddExclusion(uint index_one,
                             uint index_two){
   Index<2> index(index_one,index_two);
   Index<2> reverse_index(index_two,index_one);
-  if(index_one >= num_atoms_ || index_two >= num_atoms_){
-    throw ost::Error("Index of exclusion atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_){
+    throw ost::Error("Index of exclusion particle exceeds number of particles present in topology!");
   }
   if(added_exclusions_.find(index) != added_exclusions_.end() || added_exclusions_.find(reverse_index) != added_exclusions_.end()){
     throw ost::Error("This particular exclusion has already been set!");
@@ -223,8 +213,8 @@ uint Topology::AddDistanceConstraint(uint index_one,
                                      uint index_two,
                                      Real distance){
 
-  if(index_one >= num_atoms_ || index_two >= num_atoms_){
-    throw ost::Error("Index of distance constraint atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_){
+    throw ost::Error("Index of distance constraint atom exceeds number of particles present in topology!");
   }
   Index<2> index(index_one,index_two);
   Index<2> reverse_index(index_two,index_one);
@@ -239,8 +229,8 @@ uint Topology::AddDistanceConstraint(uint index_one,
 }
 
 void Topology::AddPositionConstraint(uint index){
-  if(index >= num_atoms_){
-    throw ost::Error("Index of position constraint exceeds number of atoms present in topology!");
+  if(index >= num_particles_){
+    throw ost::Error("Index of position constraint exceeds number of particles present in topology!");
   }
   std::vector<uint>::iterator it = std::find(position_constraints_.begin(),position_constraints_.end(),index);
   if(it != position_constraints_.end()) return; //already present
@@ -249,8 +239,8 @@ void Topology::AddPositionConstraint(uint index){
 
 uint Topology::AddHarmonicPositionRestraint(uint ind, const geom::Vec3& ref_position, Real k, 
                                             Real x_scale, Real y_scale, Real z_scale){
-  if(ind >= num_atoms_){
-    throw ost::Error("Index of harmonic distance constraint exceeds number of atoms present in topology!");
+  if(ind >= num_particles_){
+    throw ost::Error("Index of harmonic distance constraint exceeds number of particles present in topology!");
   }
   Index<1> index(ind);
   std::vector<Real> parameters;
@@ -268,8 +258,8 @@ uint Topology::AddHarmonicPositionRestraint(uint ind, const geom::Vec3& ref_posi
 
 uint Topology::AddHarmonicDistanceRestraint(uint index_one, uint index_two,
                                             Real length, Real force_constant){
-  if(index_one >= num_atoms_ || index_two >= num_atoms_){
-    throw ost::Error("Index of distance restraint atom exceeds number of atoms present in topology!");
+  if(index_one >= num_particles_ || index_two >= num_particles_){
+    throw ost::Error("Index of distance restraint atom exceeds number of particles present in topology!");
   }
   Index<2> index(index_one,index_two);
   std::vector<Real> parameters;
@@ -280,8 +270,8 @@ uint Topology::AddHarmonicDistanceRestraint(uint index_one, uint index_two,
 }
 
 void Topology::SetSigmas(const std::vector<Real>& sigmas){
-  if(sigmas.size() != num_atoms_){
-    throw ost::Error("Expect the same number of sigma parameters than atoms!");
+  if(sigmas.size() != num_particles_){
+    throw ost::Error("Expect the same number of sigma parameters than particles!");
   }
   sigmas_ = sigmas;
 }
@@ -290,15 +280,15 @@ void Topology::SetSigma(uint index, Real sigma){
   if(sigmas_.empty()){
     throw ost::Error("Modification of sigma parameter requires initial assignment of sigma parameters in globo!");
   }
-  if(index >= num_atoms_){
-    throw ost::Error("Given atom index exceeds number off available atoms!");
+  if(index >= num_particles_){
+    throw ost::Error("Given particle index exceeds number of available particles!");
   }
   sigmas_[index] = sigma;
 }
 
 void Topology::SetEpsilons(const std::vector<Real>& epsilons){
-  if(epsilons.size() != num_atoms_){
-    throw ost::Error("Expect the same number of epsilon parameters than atoms!");
+  if(epsilons.size() != num_particles_){
+    throw ost::Error("Expect the same number of epsilon parameters than particles!");
   }
   epsilons_ = epsilons;
 }
@@ -307,15 +297,15 @@ void Topology::SetEpsilon(uint index, Real epsilon){
   if(epsilons_.empty()){
     throw ost::Error("Modification of epsilon parameter requires initial assignment of epsilon parameters in globo!");
   }
-  if(index >= num_atoms_){
-    throw ost::Error("Given atom index exceeds number off available atoms!");
+  if(index >= num_particles_){
+    throw ost::Error("Given particle index exceeds number of available particles!");
   }
   epsilons_[index] = epsilon;
 }
 
 void Topology::SetGBSARadii(const std::vector<Real>& gbsa_radii){
-  if(gbsa_radii.size() != num_atoms_){
-    throw ost::Error("Expect the same number of gbsa parameters than atoms!");
+  if(gbsa_radii.size() != num_particles_){
+    throw ost::Error("Expect the same number of gbsa parameters than particles!");
   }
   gbsa_radii_ = gbsa_radii;
 }
@@ -324,15 +314,15 @@ void Topology::SetGBSARadius(uint index, Real radius){
   if(gbsa_radii_.empty()){
     throw ost::Error("Modification of gbsa radius requires initial assignment of gbsa radii in globo!");
   }
-  if(index >= num_atoms_){
-    throw ost::Error("Given atom index exceeds number off available atoms!");
+  if(index >= num_particles_){
+    throw ost::Error("Given particle index exceeds number of available particles!");
   }
   gbsa_radii_[index] = radius;
 }
 
 void Topology::SetOBCScalings(const std::vector<Real>& obc_scaling){
-  if(obc_scaling.size() != num_atoms_){
-    throw ost::Error("Expect the same number of gbsa parameters than atoms!");
+  if(obc_scaling.size() != num_particles_){
+    throw ost::Error("Expect the same number of gbsa parameters than particles!");
   }
   obc_scaling_ = obc_scaling;
 }
@@ -341,15 +331,15 @@ void Topology::SetOBCScaling(uint index, Real scaling){
   if(obc_scaling_.empty()){
     throw ost::Error("Modification of obc scaling parameter requires initial assignment of obc scaling parameters in globo!");
   }
-  if(index >= num_atoms_){
-    throw ost::Error("Given atom index exceeds number off available atoms!");
+  if(index >= num_particles_){
+    throw ost::Error("Given particle index exceeds number of available particles!");
   }
   obc_scaling_[index] = scaling;
 }
 
 void Topology::SetCharges(const std::vector<Real>& charges){
-  if(charges.size() != num_atoms_){
-    throw ost::Error("Expect the same number of charges than atoms!");
+  if(charges.size() != num_particles_){
+    throw ost::Error("Expect the same number of charges than particles!");
   }
   charges_ = charges;
 }
@@ -358,15 +348,15 @@ void Topology::SetCharge(uint index, Real charge){
   if(charges_.empty()){
     throw ost::Error("Modification of charge parameter requires initial assignment of charges in globo!");
   }
-  if(index >= num_atoms_){
-    throw ost::Error("Given atom index exceeds number off available atoms!");
+  if(index >= num_particles_){
+    throw ost::Error("Given particle index exceeds number of available particles!");
   }
   charges_[index] = charge;  
 }
 
 void Topology::SetMasses(const std::vector<Real>& masses){
-  if(masses.size() != num_atoms_){
-    throw ost::Error("Expect masses vector to have the same number of elements than atoms in the topologies entity!");
+  if(masses.size() != num_particles_){
+    throw ost::Error("Expect masses vector to have the same number of elements than particles in the topologies entity!");
   }
   atom_masses_ = masses;
 }
@@ -375,8 +365,8 @@ void Topology::SetMass(uint index, Real mass){
   if(atom_masses_.empty()){
     throw ost::Error("Modification of mass parameter requires initial assignment of masses in globo!");
   }
-  if(index >= num_atoms_){
-    throw ost::Error("Given atom index exceeds number off available atoms!");
+  if(index >= num_particles_){
+    throw ost::Error("Given particle index exceeds number of available particles!");
   }
   atom_masses_[index] = mass;
 
@@ -625,114 +615,42 @@ void Topology::SetHarmonicDistanceRestraintParameters(uint index, Real length,
 
 }
 
-uint Topology::GetAtomIndex(const ost::mol::AtomHandle& at) const{
-  std::map<long,uint>::const_iterator i = atom_index_mapper_.find(at.GetHashCode());
-  if(i == atom_index_mapper_.end()){
-    throw ost::Error("Could not find atom in attached entity!");
-  }
-  return i->second;
-}
-
-
-uint Topology::GetResidueIndex(const ost::mol::ResidueHandle& res) const{
-  std::map<long,uint>::const_iterator i = residue_index_mapper_.find(res.GetHashCode());
-  if(i == residue_index_mapper_.end()){
-    throw ost::Error("Could not find residue in attached entity!");
-  }
-  return i->second;
-}
-
-uint Topology::GetAtomIndex(uint residue_index, const String& atom_name) const{
-  if(residue_index < 0 || residue_index >= num_residues_){
-    throw ost::Error("Residue index out of bounds");
-  }
-  if(atom_name[0] == '+'){
-    ost::mol::ResidueHandle actual_res = res_list_[residue_index];
-    ost::mol::ResidueHandle next_residue = actual_res.GetNext();
-    if(!next_residue){
-      throw ost::Error("Could not get next residue!");
-    }
-    String next_atom_name = atom_name.substr(1);
-    std::map<String,uint>::const_iterator it = atom_name_mapper_[this->GetResidueIndex(next_residue)].find(next_atom_name);
-    if(it == atom_name_mapper_[this->GetResidueIndex(next_residue)].end()){
-      std::stringstream ss;
-      ss << "Could not find required atom " << next_atom_name;
-      ss << " from residue " << next_residue.GetQualifiedName();
-    }
-    return it->second;
-  }
-  if(atom_name[0] == '-'){
-    ost::mol::ResidueHandle actual_res = res_list_[residue_index];
-    ost::mol::ResidueHandle prev_residue = actual_res.GetPrev();
-    if(!prev_residue){
-      throw ost::Error("Could not get prev residue!");
-    }
-    String prev_atom_name = atom_name.substr(1);
-    std::map<String,uint>::const_iterator it = atom_name_mapper_[this->GetResidueIndex(prev_residue)].find(prev_atom_name);
-    if(it == atom_name_mapper_[this->GetResidueIndex(prev_residue)].end()){
-      std::stringstream ss;
-      ss << "Could not find required atom " << prev_atom_name;
-      ss << " from residue " << prev_residue.GetQualifiedName();
-    }
-    return it->second;
-  }
-  std::map<String,uint>::const_iterator it = atom_name_mapper_[residue_index].find(atom_name);
-  if(it == atom_name_mapper_[residue_index].end()){
-    std::stringstream ss;
-    ss << "Could not find required atom " << atom_name;
-    ss << " from residue " << res_list_[residue_index].GetQualifiedName();
-  }
-  return it->second;
-}
-
 Real Topology::GetMass(uint index) const{
-  if(index >= num_atoms_) throw ost::Error("Provided index exceeds number of atoms present in topology!");
+  if(index >= num_particles_) throw ost::Error("Provided index exceeds number of particles present in topology!");
   return atom_masses_[index];
 }
 
 Real Topology::GetSigma(uint index) const{
-  if(index >= num_atoms_) throw ost::Error("Provided index exceeds number of atoms present in topology!");
+  if(index >= num_particles_) throw ost::Error("Provided index exceeds number of particles present in topology!");
   if(sigmas_.empty()) throw ost::Error("No sigmas set!");
   return sigmas_[index];
 }
 
 Real Topology::GetEpsilon(uint index) const{
-  if(index >= num_atoms_) throw ost::Error("Provided index exceeds number of atoms present in topology!");
+  if(index >= num_particles_) throw ost::Error("Provided index exceeds number of particles present in topology!");
   if(epsilons_.empty()) throw ost::Error("No epsilons set!");
   return epsilons_[index];
 }
 
 Real Topology::GetGBSARadius(uint index) const{
-  if(index >= num_atoms_) throw ost::Error("Provided index exceeds number of atoms present in topology!");
+  if(index >= num_particles_) throw ost::Error("Provided index exceeds number of particles present in topology!");
   if(gbsa_radii_.empty()) throw ost::Error("No gbsa radii set!");
   return gbsa_radii_[index];
 }
 
 Real Topology::GetOBCScaling(uint index) const{
-  if(index >= num_atoms_) throw ost::Error("Provided index exceeds number of atoms present in topology!");
+  if(index >= num_particles_) throw ost::Error("Provided index exceeds number of particles present in topology!");
   if(obc_scaling_.empty()) throw ost::Error("No obc scalings set!");
   return obc_scaling_[index];
 }
 
 
 Real Topology::GetCharge(uint index) const{
-  if(index >= num_atoms_) throw ost::Error("Provided index exceeds number of atoms present in topology!");
+  if(index >= num_particles_) throw ost::Error("Provided index exceeds number of particles present in topology!");
   if(charges_.empty()) throw ost::Error("No charges set!");
   return charges_[index];
 }
 
-
-void Topology::SetEntityPositions(const geom::Vec3List& positions){
-  if(positions.size() != num_atoms_){
-    throw ost::Error("Expect positions vector to have the same number of elements than atoms in the topologies entity!");
-  }
-  ost::mol::XCSEditor ed = ent_.EditXCS(ost::mol::BUFFERED_EDIT);
-  for(uint i = 0; i < num_atoms_; ++i){
-    ed.SetAtomPos(atom_list_[i],positions[i]);
-  }
-  ed.UpdateICS();
-}
-
 std::vector<uint> Topology::GetHarmonicBondIndices(uint index_one,
                                                     uint index_two) const{
   Index<2> index(index_one,index_two);
@@ -990,33 +908,77 @@ std::vector<uint> Topology::GetHarmonicDistanceRestraintIndices(uint atom_index)
   return return_indices;
 }
 
-void Topology::Merge(TopologyPtr p){  
+void Topology::Merge(ost::mol::EntityHandle& ent, TopologyPtr other, 
+                     const ost::mol::EntityHandle& other_ent){
+
+  //check whether the particle numbers match
+  if(num_particles_ != static_cast<uint>(ent.GetAtomCount())){
+    throw ost::Error("Num Particles is not consistent with num atoms of provided ent!");
+  }
+
+  if(other->GetNumParticles() != static_cast<uint>(other_ent.GetAtomCount())){
+    throw ost::Error("Num Particles is not consistent with num atoms of provided other_ent!");
+  } 
+
   //check whether there is chain from the new entity is already present in the actual entity 
-  ost::mol::EntityHandle source_ent = p->GetEntity();
-  ost::mol::ChainHandleList source_chains = source_ent.GetChainList();
-  for(ost::mol::ChainHandleList::iterator i = source_chains.begin();
-      i != source_chains.end(); ++i){
-    if(ent_.FindChain(i->GetName()).IsValid()){
+  ost::mol::ChainHandleList other_chains = other_ent.GetChainList();
+  for(ost::mol::ChainHandleList::iterator i = other_chains.begin(), 
+      e = other_chains.end(); i != e; ++i){
+    if(ent.FindChain(i->GetName()).IsValid()){
       std::stringstream ss;
-      ss << "Chain with name \"" << i->GetName() << "\" from source topology";
+      ss << "Chain with name \"" << i->GetName() << "\" from other topology";
       ss << "is already present in destination topology!";
       throw ost::Error(ss.str());
     }
   }
   //check whether the fudge parameters are consistent
-  if(p->GetFudgeLJ() != fudge_lj_ || p->GetFudgeQQ() != fudge_qq_){
+  if(other->fudge_lj_ != fudge_lj_ || other->fudge_qq_ != fudge_qq_){
     throw ost::Error("Expect the fudge parameters of source topology to consistent with the fudge parameters from the destination topology!");
   }
 
+  if(!charges_.empty()){
+    if(other->charges_.empty()){
+      throw ost::Error("Cannot merge topology without charges into a topology with defined charges!");
+    }
+  }
+
+  if(!sigmas_.empty()){
+    if(other->sigmas_.empty()){
+      throw ost::Error("Cannot merge topology without lj sigmas into a topology with defined sigmas!");
+    }
+  }
+
+  if(!epsilons_.empty()){
+    if(other->epsilons_.empty()){
+      throw ost::Error("Cannot merge topology without lj epsilons into a topology with defined epsilons!");
+    }
+  }
+
+  if(!gbsa_radii_.empty()){
+    if(other->gbsa_radii_.empty()){
+      throw ost::Error("Cannot merge topology without gbsa radii into a topology with defined radii!");
+    }
+  }
+
+  if(!obc_scaling_.empty()){
+    if(other->obc_scaling_.empty()){
+      throw ost::Error("Cannot merge topology without obc scaling into a topology with defined scaling!");
+    }
+  }
+
+  uint old_num_particles = num_particles_;
+  num_particles_ = num_particles_ + other->GetNumParticles();
+
   //mapper of hashcode from source atom to index in added_atom list
   std::map<long,int> index_mapper;
   ost::mol::AtomHandleList added_atoms;
 
+
   //let's create an editor to copy over all chains, residues and atoms
-  //from the source topologies entity
-  ost::mol::XCSEditor ed = ent_.EditXCS(ost::mol::BUFFERED_EDIT);
-  for(ost::mol::ChainHandleList::iterator i = source_chains.begin();
-      i != source_chains.end(); ++i){
+  //from the other topologies entity
+  ost::mol::XCSEditor ed = ent.EditXCS(ost::mol::BUFFERED_EDIT);
+  for(ost::mol::ChainHandleList::iterator i = other_chains.begin(),
+      e = other_chains.end(); i != e; ++i){
     ost::mol::ResidueHandleList res_list = i->GetResidueList();
     ost::mol::ChainHandle added_chain = ed.InsertChain(i->GetName());
 
@@ -1034,7 +996,7 @@ void Topology::Merge(TopologyPtr p){
   }
 
   //let's rebuild the connectivity
-  ost::mol::BondHandleList bond_list = source_ent.GetBondList();
+  ost::mol::BondHandleList bond_list = other_ent.GetBondList();
   for(ost::mol::BondHandleList::iterator i = bond_list.begin();
       i != bond_list.end(); ++i){
     ed.Connect(added_atoms[index_mapper[i->GetFirst().GetHashCode()]],
@@ -1042,253 +1004,157 @@ void Topology::Merge(TopologyPtr p){
   }
   ed.UpdateICS();
 
-  //let's update the remaining structural data
-  atom_list_ = ent_.GetAtomList();
-  res_list_ = ent_.GetResidueList();
-  num_atoms_ = atom_list_.size();
-  num_residues_ = res_list_.size();
-
-  //let's update the index mappers
-  for(uint i = 0; i < num_atoms_; ++i){
-    if(atom_index_mapper_.find(atom_list_[i].GetHashCode()) == atom_index_mapper_.end()){
-      atom_index_mapper_[atom_list_[i].GetHashCode()] = i;
-      continue;
-    }
-    //indices of atoms that have already been there must not be changed
-    if(atom_index_mapper_[atom_list_[i].GetHashCode()] != i){
-      throw ost::Error("Updating indices during merging process went horribly wrong!");
-    }
-  }
-
-  for(uint i = 0; i != num_residues_; ++i){
-    if(residue_index_mapper_.find(res_list_[i].GetHashCode()) == residue_index_mapper_.end()){
-      residue_index_mapper_[res_list_[i].GetHashCode()] = i;
-      continue;
-    }
-    //indices of residues that have already been there must not be changed
-    if(residue_index_mapper_[res_list_[i].GetHashCode()] != i){
-      throw ost::Error("Updating indices during merging process went horribly wrong!");
-    }
-  }
-
-  atom_name_mapper_.clear(); //let's rebuild from scratch
-  std::map<String,uint> per_residue_mapper;
-  ost::mol::AtomHandleList residue_atom_list;
-  for(uint i = 0; i < num_residues_; ++i){
-    per_residue_mapper.clear();
-    residue_atom_list = res_list_[i].GetAtomList();
-    for(ost::mol::AtomHandleList::iterator j = residue_atom_list.begin();
-        j != residue_atom_list.end(); ++j){
-      per_residue_mapper[j->GetName()] = atom_index_mapper_[j->GetHashCode()];
-    }
-    atom_name_mapper_.push_back(per_residue_mapper);
-  }
- 
-  //let's first generate an index mapper from source atom indices to the indices in
-  //the destination entity
-  ost::mol::AtomHandleList source_atoms = source_ent.GetAtomList();
-  std::map<int,int> final_index_mapper;
-  for(uint i = 0; i < source_atoms.size(); ++i){
-    final_index_mapper[i] = atom_index_mapper_[added_atoms[index_mapper[source_atoms[i].GetHashCode()]].GetHashCode()];
-  }
-
   //let's map over masses
-  std::vector<Real> source_masses = p->GetMasses();
-
-  for(uint i = 0; i < source_atoms.size(); ++i){
-    atom_masses_.push_back(0.0);
-  }
-  for(uint i = 0; i < source_atoms.size(); ++i){
-    atom_masses_[final_index_mapper[i]] = source_masses[i];
-  }
+  atom_masses_.resize(old_num_particles + other->GetNumParticles());
+  memcpy(&atom_masses_[old_num_particles],&other->atom_masses_[0],other->GetNumParticles() * sizeof(Real));
 
   //let's map over the other per atom parameters if they're defined in the destination topology
   if(!charges_.empty()){
-    std::vector<Real> source_charges = p->GetCharges();
-    if(source_charges.empty()){
-      throw ost::Error("Cannot merge topology without charges into a topology with defined charges!");
-    }
-    for(uint i = 0; i < source_charges.size(); ++i){
-      charges_.push_back(0.0);
-    }
-    for(uint i = 0; i < source_charges.size(); ++i){
-      charges_[final_index_mapper[i]] = source_charges[i];
-    }
+    charges_.resize(old_num_particles + other->GetNumParticles());
+    memcpy(&charges_[old_num_particles],&other->charges_[0],other->GetNumParticles() * sizeof(Real));
   }
 
   if(!sigmas_.empty()){
-    std::vector<Real> source_sigmas = p->GetSigmas();
-    if(source_sigmas.empty()){
-      throw ost::Error("Cannot merge topology without lj sigmas into a topology with defined sigmas!");
-    }
-    for(uint i = 0; i < source_sigmas.size(); ++i){
-      sigmas_.push_back(0.0);
-    }
-    for(uint i = 0; i < source_sigmas.size(); ++i){
-      sigmas_[final_index_mapper[i]] = source_sigmas[i];
-    }
+    sigmas_.resize(old_num_particles + other->GetNumParticles());
+    memcpy(&sigmas_[old_num_particles],&other->sigmas_[0],other->GetNumParticles() * sizeof(Real));
   }
 
   if(!epsilons_.empty()){
-    std::vector<Real> source_epsilons= p->GetEpsilons();
-    if(source_epsilons.empty()){
-      throw ost::Error("Cannot merge topology without lj epsilons into a topology with defined epsilons!");
-    }
-    for(uint i = 0; i < source_epsilons.size(); ++i){
-      epsilons_.push_back(0.0);
-    }
-    for(uint i = 0; i < source_epsilons.size(); ++i){
-      epsilons_[final_index_mapper[i]] = source_epsilons[i];
-    }
+    epsilons_.resize(old_num_particles + other->GetNumParticles());
+    memcpy(&epsilons_[old_num_particles],&other->epsilons_[0],other->GetNumParticles() * sizeof(Real));
   }
 
   if(!gbsa_radii_.empty()){
-    std::vector<Real> source_radii= p->GetGBSARadii();
-    if(source_radii.empty()){
-      throw ost::Error("Cannot merge topology without gbsa radii into a topology with defined radii!");
-    }
-    for(uint i = 0; i < source_radii.size(); ++i){
-      gbsa_radii_.push_back(0.0);
-    }
-    for(uint i = 0; i < source_radii.size(); ++i){
-      gbsa_radii_[final_index_mapper[i]] = source_radii[i];
-    }
+    gbsa_radii_.resize(old_num_particles + other->GetNumParticles());
+    memcpy(&gbsa_radii_[old_num_particles],&other->gbsa_radii_[0],other->GetNumParticles() * sizeof(Real));
   }
 
   if(!obc_scaling_.empty()){
-    std::vector<Real> source_scaling= p->GetOBCScalings();
-    if(source_scaling.empty()){
-      throw ost::Error("Cannot merge topology without obc scaling into a topology with defined scaling!");
-    }
-    for(uint i = 0; i < source_scaling.size(); ++i){
-      obc_scaling_.push_back(0.0);
-    }
-    for(uint i = 0; i < source_scaling.size(); ++i){
-      obc_scaling_[final_index_mapper[i]] = source_scaling[i];
-    }
+    obc_scaling_.resize(old_num_particles + other->GetNumParticles());
+    memcpy(&obc_scaling_[old_num_particles],&other->obc_scaling_[0],other->GetNumParticles() * sizeof(Real));
   }
 
   //finally, all the interactions get mapped over
 
-  const std::vector<std::pair<Index<2>, std::vector<Real> > >& harmonic_bonds = p->GetHarmonicBonds();
+  const std::vector<std::pair<Index<2>, std::vector<Real> > >& harmonic_bonds = other->GetHarmonicBonds();
   if(!harmonic_bonds.empty()){
     for(std::vector<std::pair<Index<2>, std::vector<Real> > >::const_iterator i = harmonic_bonds.begin();
         i != harmonic_bonds.end(); ++i){
-      this->AddHarmonicBond(final_index_mapper[i->first[0]],
-                            final_index_mapper[i->first[1]],
+      this->AddHarmonicBond(old_num_particles + i->first[0],
+                            old_num_particles + i->first[1],
                             i->second[0],i->second[1]);
     }
   }
 
-  const std::vector<std::pair<Index<3>, std::vector<Real> > >& harmonic_angles = p->GetHarmonicAngles();
+  const std::vector<std::pair<Index<3>, std::vector<Real> > >& harmonic_angles = other->GetHarmonicAngles();
   if(!harmonic_angles.empty()){
     for(std::vector<std::pair<Index<3>, std::vector<Real> > >::const_iterator i = harmonic_angles.begin();
         i != harmonic_angles.end(); ++i){
-      this->AddHarmonicAngle(final_index_mapper[i->first[0]],
-                             final_index_mapper[i->first[1]],
-                             final_index_mapper[i->first[2]],
+      this->AddHarmonicAngle(old_num_particles + i->first[0],
+                             old_num_particles + i->first[1],
+                             old_num_particles + i->first[2],
                              i->second[0],i->second[1]);
     }
   }
 
-  const std::vector<std::pair<Index<3>, std::vector<Real> > >& urey_bradley_angles = p->GetUreyBradleyAngles();
+  const std::vector<std::pair<Index<3>, std::vector<Real> > >& urey_bradley_angles = other->GetUreyBradleyAngles();
   if(!urey_bradley_angles.empty()){
     for(std::vector<std::pair<Index<3>, std::vector<Real> > >::const_iterator i = urey_bradley_angles.begin();
         i != urey_bradley_angles.end(); ++i){
-      this->AddUreyBradleyAngle(final_index_mapper[i->first[0]],
-                                final_index_mapper[i->first[1]],
-                                final_index_mapper[i->first[2]],
+      this->AddUreyBradleyAngle(old_num_particles + i->first[0],
+                                old_num_particles + i->first[1],
+                                old_num_particles + i->first[2],
                                 i->second[0],i->second[1],
                                 i->second[2],i->second[3]);
     }
   }
 
-  const std::vector<std::pair<Index<4>, std::vector<Real> > >& periodic_dihedrals = p->GetPeriodicDihedrals();
+  const std::vector<std::pair<Index<4>, std::vector<Real> > >& periodic_dihedrals = other->GetPeriodicDihedrals();
   if(!periodic_dihedrals.empty()){
     for(std::vector<std::pair<Index<4>, std::vector<Real> > >::const_iterator i = periodic_dihedrals.begin();
         i != periodic_dihedrals.end(); ++i){
-      this->AddPeriodicDihedral(final_index_mapper[i->first[0]],
-                                final_index_mapper[i->first[1]],
-                                final_index_mapper[i->first[2]],
-                                final_index_mapper[i->first[3]],
+      this->AddPeriodicDihedral(old_num_particles + i->first[0],
+                                old_num_particles + i->first[1],
+                                old_num_particles + i->first[2],
+                                old_num_particles + i->first[3],
                                 int(i->second[0]),i->second[1],i->second[2]);
     }
   }
 
-  const std::vector<std::pair<Index<4>, std::vector<Real> > >& periodic_impropers = p->GetPeriodicImpropers();
+  const std::vector<std::pair<Index<4>, std::vector<Real> > >& periodic_impropers = other->GetPeriodicImpropers();
   if(!periodic_impropers.empty()){
     for(std::vector<std::pair<Index<4>, std::vector<Real> > >::const_iterator i = periodic_impropers.begin();
         i != periodic_impropers.end(); ++i){
-      this->AddPeriodicImproper(final_index_mapper[i->first[0]],
-                                final_index_mapper[i->first[1]],
-                                final_index_mapper[i->first[2]],
-                                final_index_mapper[i->first[3]],
+      this->AddPeriodicImproper(old_num_particles + i->first[0],
+                                old_num_particles + i->first[1],
+                                old_num_particles + i->first[2],
+                                old_num_particles + i->first[3],
                                 int(i->second[0]),i->second[1],i->second[2]);
     }
   }
 
-  const std::vector<std::pair<Index<4>, std::vector<Real> > >& harmonic_impropers = p->GetHarmonicImpropers();
+  const std::vector<std::pair<Index<4>, std::vector<Real> > >& harmonic_impropers = other->GetHarmonicImpropers();
   if(!harmonic_impropers.empty()){
     for(std::vector<std::pair<Index<4>, std::vector<Real> > >::const_iterator i = harmonic_impropers.begin();
         i != harmonic_impropers.end(); ++i){
-      this->AddHarmonicImproper(final_index_mapper[i->first[0]],
-                                final_index_mapper[i->first[1]],
-                                final_index_mapper[i->first[2]],
-                                final_index_mapper[i->first[3]],
+      this->AddHarmonicImproper(old_num_particles + i->first[0],
+                                old_num_particles + i->first[1],
+                                old_num_particles + i->first[2],
+                                old_num_particles + i->first[3],
                                 i->second[0],i->second[1]);
     }
   }
 
-  const std::vector<std::pair<Index<5>, std::vector<Real> > >& cmaps = p->GetCMaps();
+  const std::vector<std::pair<Index<5>, std::vector<Real> > >& cmaps = other->GetCMaps();
   if(!cmaps.empty()){
     for(std::vector<std::pair<Index<5>, std::vector<Real> > >::const_iterator i = cmaps.begin();
         i != cmaps.end(); ++i){
       uint dimension = sqrt(i->second.size());
-      this->AddCMap(final_index_mapper[i->first[0]],
-                    final_index_mapper[i->first[1]],
-                    final_index_mapper[i->first[2]],
-                    final_index_mapper[i->first[3]],
-                    final_index_mapper[i->first[4]],
+      this->AddCMap(old_num_particles + i->first[0],
+                    old_num_particles + i->first[1],
+                    old_num_particles + i->first[2],
+                    old_num_particles + i->first[3],
+                    old_num_particles + i->first[4],
                     dimension,i->second);
     }
   }
 
-  const std::vector<std::pair<Index<2>, std::vector<Real> > >& ljpairs = p->GetLJPairs();
+  const std::vector<std::pair<Index<2>, std::vector<Real> > >& ljpairs = other->GetLJPairs();
   if(!ljpairs.empty()){
     for(std::vector<std::pair<Index<2>, std::vector<Real> > >::const_iterator i = ljpairs.begin();
         i != ljpairs.end(); ++i){
-      this->AddLJPair(final_index_mapper[i->first[0]],
-                      final_index_mapper[i->first[1]],
+      this->AddLJPair(old_num_particles + i->first[0],
+                      old_num_particles + i->first[1],
                       i->second[0],i->second[1]);
     }
   }
 
 
-  const std::vector<std::pair<Index<2>, std::vector<Real> > >& distance_constraints = p->GetDistanceConstraints();
+  const std::vector<std::pair<Index<2>, std::vector<Real> > >& distance_constraints = other->GetDistanceConstraints();
   if(!distance_constraints.empty()){
     for(std::vector<std::pair<Index<2>, std::vector<Real> > >::const_iterator i = distance_constraints.begin();
         i != distance_constraints.end(); ++i){
-      this->AddDistanceConstraint(final_index_mapper[i->first[0]],
-                                  final_index_mapper[i->first[1]],
+      this->AddDistanceConstraint(old_num_particles + i->first[0],
+                                  old_num_particles + i->first[1],
                                   i->second[0]);
     }
   }
 
-  const std::vector<Index<2> >& exclusions = p->GetExclusions();
+  const std::vector<Index<2> >& exclusions = other->GetExclusions();
   if(!exclusions.empty()){
     for(std::vector<Index<2> >::const_iterator i = exclusions.begin();
         i != exclusions.end(); ++i){
-      this->AddExclusion(final_index_mapper[(*i)[0]],
-                          final_index_mapper[(*i)[1]]);
+      this->AddExclusion(old_num_particles + (*i)[0],
+                         old_num_particles + (*i)[1]);
     }
   }
 
-  const std::vector<std::pair<Index<1>,std::vector<Real> > >& harmonic_position_restraints = p->GetHarmonicPositionRestraints();
+  const std::vector<std::pair<Index<1>,std::vector<Real> > >& harmonic_position_restraints = other->GetHarmonicPositionRestraints();
   if(!harmonic_position_restraints.empty()){
     for(std::vector<std::pair<Index<1>,std::vector<Real> > >::const_iterator i = harmonic_position_restraints.begin();
         i != harmonic_position_restraints.end(); ++i){
       geom::Vec3 ref_pos(i->second[0],i->second[1],i->second[2]);
-      this->AddHarmonicPositionRestraint(final_index_mapper[i->first[0]],
+      this->AddHarmonicPositionRestraint(old_num_particles + i->first[0],
                                          ref_pos,
                                          i->second[3],
                                          i->second[4],
@@ -1297,66 +1163,43 @@ void Topology::Merge(TopologyPtr p){
     }
   }
 
-  const std::vector<std::pair<Index<2>,std::vector<Real> > >& harmonic_distance_restraints = p->GetHarmonicDistanceRestraints();
+  const std::vector<std::pair<Index<2>,std::vector<Real> > >& harmonic_distance_restraints = other->GetHarmonicDistanceRestraints();
   if(!harmonic_distance_restraints.empty()){
     for(std::vector<std::pair<Index<2>,std::vector<Real> > >::const_iterator i = harmonic_distance_restraints.begin();
         i != harmonic_distance_restraints.end(); ++i){
-      this->AddHarmonicDistanceRestraint(final_index_mapper[i->first[0]],
-                                         final_index_mapper[i->first[1]],
+      this->AddHarmonicDistanceRestraint(old_num_particles + i->first[0],
+                                         old_num_particles + i->first[1],
                                          i->second[0],
                                          i->second[1]);
     }
   }
 
-  const std::vector<uint>& position_constraints = p->GetPositionConstraints();
+  const std::vector<uint>& position_constraints = other->GetPositionConstraints();
   if(!position_constraints.empty()){
     for(std::vector<uint>::const_iterator i = position_constraints.begin();
         i != position_constraints.end(); ++i){
-      this->AddPositionConstraint(final_index_mapper[*i]);
+      this->AddPositionConstraint(old_num_particles + (*i));
     }
   }
 
-  for(std::set<Index<2> >::iterator i = p->added_lj_pairs_.begin();
-      i != p->added_lj_pairs_.end(); ++i){
-    added_lj_pairs_.insert(Index<2>(final_index_mapper[(*i)[0]],
-                                    final_index_mapper[(*i)[1]]));
-  }
-
-  for(std::set<Index<2> >::iterator i = p->added_distance_constraints_.begin();
-      i != p->added_distance_constraints_.end(); ++i){
-    added_distance_constraints_.insert(Index<2>(final_index_mapper[(*i)[0]],
-                                                final_index_mapper[(*i)[1]]));
+  for(std::set<Index<2> >::iterator i = other->added_lj_pairs_.begin();
+      i != other->added_lj_pairs_.end(); ++i){
+    added_lj_pairs_.insert(Index<2>(old_num_particles + (*i)[0],
+                                    old_num_particles + (*i)[1]));
   }
 
-  for(std::set<Index<2> >::iterator i = p->added_exclusions_.begin();
-      i != p->added_exclusions_.end(); ++i){
-    added_exclusions_.insert(Index<2>(final_index_mapper[(*i)[0]],
-                                      final_index_mapper[(*i)[1]]));
+  for(std::set<Index<2> >::iterator i = other->added_distance_constraints_.begin();
+      i != other->added_distance_constraints_.end(); ++i){
+    added_distance_constraints_.insert(Index<2>(old_num_particles + (*i)[0],
+                                                old_num_particles + (*i)[1]));
   }
 
-}
-
-void Topology::InitMappers(){
-  for(uint i = 0; i < atom_list_.size(); ++i){
-    atom_index_mapper_[atom_list_[i].GetHashCode()] = i;
-  }  
-
-  for(uint i = 0; i < res_list_.size(); ++i){
-    residue_index_mapper_[res_list_[i].GetHashCode()] = i;
+  for(std::set<Index<2> >::iterator i = other->added_exclusions_.begin();
+      i != other->added_exclusions_.end(); ++i){
+    added_exclusions_.insert(Index<2>(old_num_particles + (*i)[0],
+                                      old_num_particles + (*i)[1]));
   }
 
-  std::map<String,uint> per_residue_mapper;
-  ost::mol::AtomHandleList residue_atom_list;
-  for(uint i = 0; i < num_residues_; ++i){
-    per_residue_mapper.clear();
-    residue_atom_list = res_list_[i].GetAtomList();
-    for(ost::mol::AtomHandleList::iterator j = residue_atom_list.begin();
-        j != residue_atom_list.end(); ++j){
-      per_residue_mapper[j->GetName()] = atom_index_mapper_[j->GetHashCode()];
-    }
-    atom_name_mapper_.push_back(per_residue_mapper);
-  }
 }
 
-
 }}}//ns
diff --git a/modules/mol/mm/src/topology.hh b/modules/mol/mm/src/topology.hh
index 93284524c65bf24e522a4ea0df1a71cdc622f9b7..5298a60071478ade883b17c3086cff13da3ef8af 100644
--- a/modules/mol/mm/src/topology.hh
+++ b/modules/mol/mm/src/topology.hh
@@ -36,7 +36,7 @@ class Topology{
 
 public:
 
-  Topology(const ost::mol::EntityHandle& ent, const std::vector<Real>& masses);
+  Topology(const std::vector<Real>& masses);
 
   Topology() { } //should not be accessible from Python to avoid messing around
                  //with empty topology
@@ -206,14 +206,6 @@ public:
 
   void SetHarmonicDistanceRestraintParameters(uint index, Real length, Real force_constant);
 
-  ost::mol::EntityHandle GetEntity() const { return ent_; };
-
-  uint GetAtomIndex(const ost::mol::AtomHandle& at) const;
-
-  uint GetAtomIndex(uint residue_index, const String& atom_name) const;
-
-  uint GetResidueIndex(const ost::mol::ResidueHandle& res) const;
-
   const std::vector<std::pair<Index<2>, std::vector<Real> > >& GetHarmonicBonds() const { return harmonic_bonds_; }
  
   const std::vector<std::pair<Index<3>, std::vector<Real> > >& GetHarmonicAngles() const { return harmonic_angles_; }
@@ -334,9 +326,7 @@ public:
 
   std::vector<uint> GetHarmonicDistanceRestraintIndices(uint atom_index) const;
 
-  uint GetNumAtoms() { return num_atoms_; }
-
-  uint GetNumResidues() { return num_residues_; }
+  uint GetNumParticles() { return num_particles_; }
 
   uint GetNumHarmonicBonds() { return harmonic_bonds_.size(); }
 
@@ -364,132 +354,16 @@ public:
 
   uint GetNumExclusions() { return exclusions_.size(); }
 
-  void Merge(TopologyPtr p);
-
-  void SetEntityPositions(const geom::Vec3List& positions);
+  void Merge(ost::mol::EntityHandle& ent, TopologyPtr other, const ost::mol::EntityHandle& other_ent);
 
   template <typename DS>
   void Serialize(DS& ds){
-    
-    ds & num_atoms_;
-    ds & num_residues_;
-
-    uint num_chains = 0;
-    uint num_residues = 0;
-    uint num_atoms = 0;
-    uint num_bonded_atoms = 0;
-    Real x_pos = 0.0;
-    Real y_pos = 0.0;
-    Real z_pos = 0.0;
-    Real bfac = 0.0;
-    Real occ = 0.0;
-    bool is_hetatm = false;
-    String chain_name = "X";
-    String res_name = "XXX";
-    int resnum_num = 0;
-    char resnum_code = '\0';
-    String atom_name = "X";
-    String atom_element = "X";
-    uint atom_index = 0;
+
+
     uint num_items = 0;
     Index<2> actual_index;
 
-    if(ds.IsSource()){
-      ost::mol::EntityHandle ent = ost::mol::CreateEntity();
-      ost::mol::XCSEditor ed = ent.EditXCS();
-      ds & num_chains;
-      for(uint i = 0; i < num_chains; ++i){
-        ds & chain_name;
-        ds & num_residues;
-        ost::mol::ChainHandle chain = ed.InsertChain(chain_name);
-        for(uint j = 0; j < num_residues; ++j){
-          ds & res_name;
-          ds & resnum_num;
-          ds & resnum_code;
-          ds & num_atoms;
-          ost::mol::ResNum num(resnum_num,resnum_code);
-          ost::mol::ResidueHandle res = ed.AppendResidue(chain,res_name,num);
-          for(uint k = 0; k < num_atoms; ++k){
-            ds & atom_name;
-            ds & atom_element;
-            ds & x_pos;
-            ds & y_pos;
-            ds & z_pos;
-            ds & bfac;
-            ds & occ;
-            ds & is_hetatm;
-            geom::Vec3 pos(x_pos,y_pos,z_pos);
-            ed.InsertAtom(res,atom_name,pos,atom_element,occ,bfac,is_hetatm);
-          }
-        }
-      }
-      ent_ = ent;
-      atom_list_ = ent_.GetAtomList();
-      res_list_ = ent_.GetResidueList();
-      for(uint i = 0; i < atom_list_.size(); ++i){
-        ds & num_bonded_atoms;
-        for(uint j = 0; j < num_bonded_atoms; ++j){
-          ds & atom_index;
-          ed.Connect(atom_list_[i],atom_list_[atom_index]);
-        }
-      }
-      InitMappers();
-    }
-    else{
-      num_chains = ent_.GetChainCount();
-      ds & num_chains;
-      ost::mol::ChainHandleList chain_list = ent_.GetChainList();
-      for(ost::mol::ChainHandleList::iterator i = chain_list.begin();
-          i != chain_list.end(); ++i){
-        chain_name = i->GetName();
-        num_residues = i->GetResidueCount();
-        ds & chain_name;
-        ds & num_residues;
-        ost::mol::ResidueHandleList res_list = i->GetResidueList();
-        for(ost::mol::ResidueHandleList::iterator j = res_list.begin();
-            j != res_list.end(); ++j){
-          res_name = j->GetKey();
-          resnum_num = j->GetNumber().GetNum();
-          resnum_code = j->GetNumber().GetInsCode(); 
-          num_atoms = j->GetAtomCount();
-          ds & res_name;
-          ds & resnum_num;
-          ds & resnum_code;
-          ds & num_atoms;
-          ost::mol::AtomHandleList atom_list = j->GetAtomList();
-          for(ost::mol::AtomHandleList::iterator k = atom_list.begin();
-              k != atom_list.end(); ++k){
-            atom_name = k->GetName();
-            atom_element = k->GetElement();
-            geom::Vec3 pos = k->GetPos();
-            bfac = k->GetBFactor();
-            occ = k->GetOccupancy();
-            is_hetatm = k->IsHetAtom();
-            ds & atom_name;
-            ds & atom_element;
-            ds & pos[0];
-            ds & pos[1];
-            ds & pos[2];
-            ds & bfac;
-            ds & occ;
-            ds & is_hetatm;
-          }
-        }
-      }
-      ost::mol::AtomHandleList bonded_atoms;
-      for(ost::mol::AtomHandleList::iterator i = atom_list_.begin();
-          i != atom_list_.end(); ++i){
-        bonded_atoms = i->GetBondPartners();
-        num_bonded_atoms = bonded_atoms.size();
-        ds & num_bonded_atoms;
-        for(ost::mol::AtomHandleList::iterator j = bonded_atoms.begin();
-            j != bonded_atoms.end(); ++j){
-          atom_index = this->GetAtomIndex(*j);
-          ds & atom_index;
-        }
-      }
-    }
-
+    ds & num_particles_;
     ds & fudge_qq_;
     ds & fudge_lj_;
 
@@ -791,18 +665,7 @@ public:
 
 private:
 
-  void InitMappers();
-
-  ost::mol::EntityHandle ent_;
-  ost::mol::AtomHandleList atom_list_;
-  ost::mol::ResidueHandleList res_list_;
-
-  uint num_atoms_;
-  uint num_residues_;
-
-  std::map<long,uint> atom_index_mapper_;
-  std::map<long,uint> residue_index_mapper_;
-  std::vector<std::map<String,uint> > atom_name_mapper_;
+  uint num_particles_;
 
   //fudge parameters for lj 1,4 pairs
   Real fudge_qq_;
diff --git a/modules/mol/mm/src/topology_creator.cc b/modules/mol/mm/src/topology_creator.cc
index d21526c2428ae47b01c28e658d1d32bd7a85462a..0c8742d9cf3022c088a3ee99910c7c36ceb9e4f3 100644
--- a/modules/mol/mm/src/topology_creator.cc
+++ b/modules/mol/mm/src/topology_creator.cc
@@ -1,25 +1,70 @@
 #include <ost/mol/mm/topology_creator.hh>
 #include <ost/mol/mm/mm_modeller.hh>
 
+namespace{
+
+int GetAtomIndex(uint residue_index, 
+                 ost::mol::ResidueHandleList& res_list, 
+                 std::map<long,int>& atom_indices,
+                 const String& atom_name){
+
+  if(atom_name[0] == '+'){
+    ost::mol::ResidueHandle actual_res = res_list[residue_index];
+    ost::mol::ResidueHandle next_residue = actual_res.GetNext();
+    if(!next_residue){
+      throw ost::Error("Could not find required atom when settings up topology!");
+    }
+    String next_atom_name = atom_name.substr(1);
+    ost::mol::AtomHandle at = next_residue.FindAtom(next_atom_name);
+    if(!at.IsValid()){
+      throw ost::Error("Could not find required atom when settings up topology!");      
+    }
+    return atom_indices[at.GetHashCode()];
+  }
+
+  if(atom_name[0] == '-'){
+    ost::mol::ResidueHandle actual_res = res_list[residue_index];
+    ost::mol::ResidueHandle prev_residue = actual_res.GetPrev();
+    if(!prev_residue){
+      throw ost::Error("Could not find required atom when settings up topology!");
+    }
+    String prev_atom_name = atom_name.substr(1);
+    ost::mol::AtomHandle at = prev_residue.FindAtom(prev_atom_name);
+    if(!at.IsValid()){
+      throw ost::Error("Could not find required atom when settings up topology!");      
+    }
+    return atom_indices[at.GetHashCode()];
+  }
+
+  ost::mol::AtomHandle at = res_list[residue_index].FindAtom(atom_name);
+  if(!at.IsValid()){
+    throw ost::Error("Could not find required atom when settings up topology!");      
+  }
+  return atom_indices[at.GetHashCode()];
+}
+
+
+
+
+
+
+}
+
 namespace ost{ namespace mol{ namespace mm{
 
-TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle, 
+TopologyPtr TopologyCreator::Create(ost::mol::EntityHandle& ent, 
                                     const MMSettingsPtr settings){
 
   if(settings->constrain_hangles == true && settings->constrain_hbonds == false){
     throw ost::Error("If hangles is true, hbonds must also be true in settings object!");
   }
 
-  //make sure the input gets not modified
-  ost::mol::EntityHandle ent = handle.Copy();
   ost::mol::ResidueHandleList res_list = ent.GetResidueList();
   ost::mol::XCSEditor ed = ent.EditXCS(ost::mol::BUFFERED_EDIT);
 
-
   //rename all the stuff to the gromacs naming...
   MMModeller::AssignGromacsNaming(ent);
 
-
   //even if not reconnected yet, it gets assumed, that
   //peptide or nucleotide bonds are correctly set in the input entity
   //to allow the identification and tagging of terminal residues
@@ -137,6 +182,7 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
   ed.UpdateICS();
 
 
+  //check, whether all building blocks match now with the previously applied modifications
   String match_fail_info;
   for(unsigned int i = 0; i < building_blocks.size(); ++i){
     if(!building_blocks[i]->Match(res_list[i], false, match_fail_info)){
@@ -149,14 +195,31 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
 
   std::vector<Real> initial_masses(ent.GetAtomCount());
 
-  TopologyPtr top = TopologyPtr(new Topology(ent,initial_masses));
-  ent = top->GetEntity();
+  TopologyPtr top = TopologyPtr(new Topology(initial_masses));
 
   //note, that we have to get the residue list again, since there is a new entity handle
   //created when initializing the topology
   res_list = ent.GetResidueList();
   ost::mol::AtomHandleList atom_list = ent.GetAtomList();
 
+  std::map<long,int> atom_indices;
+  std::map<long,int> residue_indices;
+  int actual_index = 0;
+
+  for(ost::mol::AtomHandleList::const_iterator i = atom_list.begin(), e = atom_list.end(); 
+      i != e; ++i){
+    atom_indices[i->GetHashCode()] = actual_index;
+    ++actual_index;
+  }
+
+  actual_index = 0;
+  for(ost::mol::ResidueHandleList::const_iterator i = res_list.begin(), e = res_list.end(); 
+      i != e; ++i){
+    residue_indices[i->GetHashCode()] = actual_index;
+    ++actual_index;
+  }
+
+
   //let's extract atom specific informations
   //extract information about single atoms
   std::vector<std::vector<uint> > bound_to;
@@ -168,27 +231,23 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
   std::vector<String> residue_names_of_atoms;
   ost::mol::AtomHandleList bonded_atoms;
   std::vector<uint> temp;
-  //extract masses, types, charges and bondpartners
+  int residue_index;
 
-  ost::mol::ResidueHandle temp_res;
+  //extract masses, types, charges and bondpartners
 
   for(ost::mol::AtomHandleList::iterator i = atom_list.begin();
       i!=atom_list.end();++i){
-    temp_res = i->GetResidue();
-    residue_names_of_atoms.push_back(temp_res.GetName());
+    residue_index = residue_indices[i->GetResidue().GetHashCode()];
+    residue_names_of_atoms.push_back(i->GetResidue().GetName());
     atom_elements.push_back(i->GetElement());
-    atom_types.push_back(building_blocks[top->GetResidueIndex(temp_res)]->GetType(i->GetName()));
-    atom_charges.push_back(building_blocks[top->GetResidueIndex(temp_res)]->GetCharge(i->GetName()));
-    //atom_masses.push_back(building_blocks[top->GetResidueIndex(temp_res)]->GetMass(i->GetName()));
-    //if( atom_masses.back() != atom_masses.back() ){
-    //  atom_masses.back() = ff->GetMass(atom_types.back());  
-    //}
+    atom_types.push_back(building_blocks[residue_index]->GetType(i->GetName()));
+    atom_charges.push_back(building_blocks[residue_index]->GetCharge(i->GetName()));
     atom_masses.push_back(ff->GetMass(atom_types.back()));
     bonded_atoms = i->GetBondPartners();
     temp.clear();
     for(ost::mol::AtomHandleList::iterator j = bonded_atoms.begin();
         j!=bonded_atoms.end();++j){
-      temp.push_back(top->GetAtomIndex(*j));
+      temp.push_back(atom_indices[j->GetHashCode()]);
     }
     bound_to.push_back(temp);
   }
@@ -297,7 +356,7 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
         interaction_atom_names = (*j)->GetNames();
         Index<4> improper_index;
         for(uint k = 0; k < 4; ++k){
-          improper_index[k] = top->GetAtomIndex(i,interaction_atom_names[k]);
+          improper_index[k] = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[k]);
         }
         impropers.insert(improper_index);
       }
@@ -312,7 +371,7 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
         interaction_atom_names = (*j)->GetNames();
         Index<5> cmap_index;
         for(uint k = 0; k < 5; ++k){
-          cmap_index[k] = top->GetAtomIndex(i,interaction_atom_names[k]);
+          cmap_index[k] = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[k]);
         }
         cmaps.insert(cmap_index);
       }
@@ -325,8 +384,8 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
       for(std::vector<MMInteractionPtr>::iterator j = interaction_list.begin();
           j != interaction_list.end(); ++j){
         interaction_atom_names = (*j)->GetNames();
-        uint one = top->GetAtomIndex(i,interaction_atom_names[0]);
-        uint two = top->GetAtomIndex(i,interaction_atom_names[1]);
+        uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+        uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]);
         exclusions.insert(Index<2>(std::min(one,two),std::max(one,two)));
       }
     }
@@ -372,8 +431,8 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
     for(std::vector<MMInteractionPtr>::iterator j = interaction_list.begin();
         j != interaction_list.end(); ++j){
       interaction_atom_names = (*j)->GetNames();
-      uint one = top->GetAtomIndex(i, interaction_atom_names[0]);
-      uint two = top->GetAtomIndex(i, interaction_atom_names[1]);
+      uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+      uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]);
       Index<2> index(std::min(one,two),std::max(one,two));
       //we don't want contradicting constraints!
       if(distance_constraints.find(index) != distance_constraints.end()) continue;
@@ -528,14 +587,14 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
   std::vector<String> types;
   if(settings->add_bonds){
     //handle bonds
-    for(uint i = 0; i < top->GetNumResidues(); ++i){
+    for(uint i = 0; i < res_list.size(); ++i){
       interaction_list = building_blocks[i]->GetBonds();
       for(std::vector<MMInteractionPtr>::iterator j = interaction_list.begin();
           j != interaction_list.end(); ++j){
         if((*j)->IsParametrized()){
           interaction_atom_names = (*j)->GetNames();
-          uint one = top->GetAtomIndex(i, interaction_atom_names[0]);
-          uint two = top->GetAtomIndex(i, interaction_atom_names[1]); 
+          uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+          uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]); 
           Index<2> bond_index(std::min(one,two),std::max(one,two));
           bonds.erase(bond_index);
           //There are only harmonic bonds supported
@@ -571,15 +630,15 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
 
   if(settings->add_angles){
     //handle angles
-    for(uint i = 0; i < top->GetNumResidues(); ++i){
+    for(uint i = 0; i < res_list.size(); ++i){
       interaction_list = building_blocks[i]->GetAngles();
       for(std::vector<MMInteractionPtr>::iterator j = interaction_list.begin();
           j != interaction_list.end(); ++j){
         if((*j)->IsParametrized()){
           interaction_atom_names = (*j)->GetNames();
-          uint one = top->GetAtomIndex(i, interaction_atom_names[0]);
-          uint two = top->GetAtomIndex(i, interaction_atom_names[1]); 
-          uint three = top->GetAtomIndex(i, interaction_atom_names[2]);
+          uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+          uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]); 
+          uint three = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[2]);
           Index<3> angle_index(std::min(one,three),two,std::max(one,three));
           if(angles.find(angle_index) == angles.end()){
             std::stringstream ss;
@@ -653,16 +712,16 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
   if(settings->add_dihedrals){
     //handle dihedrals
     std::set<Index<4> > dihedrals_to_delete;
-    for(uint i = 0; i < top->GetNumResidues(); ++i){
+    for(uint i = 0; i < res_list.size(); ++i){
       interaction_list = building_blocks[i]->GetDihedrals();
       for(std::vector<MMInteractionPtr>::iterator j = interaction_list.begin();
           j != interaction_list.end(); ++j){
         if((*j)->IsParametrized()){
           interaction_atom_names = (*j)->GetNames();
-          uint one = top->GetAtomIndex(i, interaction_atom_names[0]);
-          uint two = top->GetAtomIndex(i, interaction_atom_names[1]); 
-          uint three = top->GetAtomIndex(i, interaction_atom_names[2]);
-          uint four = top->GetAtomIndex(i, interaction_atom_names[3]);
+          uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+          uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]); 
+          uint three = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[2]);
+          uint four = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[3]);
           Index<4> dihedral_index;
           if(one<four){
             dihedral_index[0] = one;
@@ -743,10 +802,10 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
           //impropers from the building block => they will have the same ordering 
           //anyway  (or at least hopefully ;) 
           interaction_atom_names = (*j)->GetNames();
-          uint one = top->GetAtomIndex(i, interaction_atom_names[0]);
-          uint two = top->GetAtomIndex(i, interaction_atom_names[1]); 
-          uint three = top->GetAtomIndex(i, interaction_atom_names[2]);
-          uint four = top->GetAtomIndex(i, interaction_atom_names[3]);
+          uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+          uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]); 
+          uint three = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[2]);
+          uint four = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[3]);
           Index<4> improper_index(one,two,three,four);
           impropers_to_delete.insert(improper_index);
           parameters = (*j)->GetParam();
@@ -824,7 +883,7 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
 
   if(settings->add_cmaps){
     //handle cmaps
-    for(uint i = 0; i < top->GetNumResidues(); ++i){
+    for(uint i = 0; i < res_list.size(); ++i){
       interaction_list = building_blocks[i]->GetCMaps();
       for(std::vector<MMInteractionPtr>::iterator j = interaction_list.begin();
           j != interaction_list.end(); ++j){
@@ -833,11 +892,11 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
           //cmaps from the building block => they will have the same ordering 
           //anyway
           interaction_atom_names = (*j)->GetNames();
-          uint one = top->GetAtomIndex(i, interaction_atom_names[0]);
-          uint two = top->GetAtomIndex(i, interaction_atom_names[1]); 
-          uint three = top->GetAtomIndex(i, interaction_atom_names[2]);
-          uint four = top->GetAtomIndex(i, interaction_atom_names[3]);
-          uint five = top->GetAtomIndex(i, interaction_atom_names[4]);
+          uint one = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[0]);
+          uint two = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[1]); 
+          uint three = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[2]);
+          uint four = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[3]);
+          uint five = GetAtomIndex(i,res_list,atom_indices,interaction_atom_names[4]);
           Index<5> cmap_index(one,two,three,four,five);
           cmaps.erase(cmap_index);
           parameters = (*j)->GetParam();
@@ -924,7 +983,7 @@ TopologyPtr TopologyCreator::Create(const ost::mol::EntityHandle& handle,
   if(settings->add_gbsa){
     std::vector<Real> radii;
     std::vector<Real> scaling;
-    for(uint i = 0; i < top->GetNumAtoms(); ++i){
+    for(uint i = 0; i < atom_list.size(); ++i){
       MMInteractionPtr gbsa_ptr = ff->GetImplicitGenborn(atom_types[i]);
       if(!gbsa_ptr){
         std::stringstream ss;
diff --git a/modules/mol/mm/src/topology_creator.hh b/modules/mol/mm/src/topology_creator.hh
index 1d719a4144b7a3571e8bb2e2764bf4327a2d06e2..52cd37938c25537418e6561897c929f95101c20e 100644
--- a/modules/mol/mm/src/topology_creator.hh
+++ b/modules/mol/mm/src/topology_creator.hh
@@ -29,7 +29,7 @@ typedef boost::shared_ptr<TopologyCreator> TopologyCreatorPtr;
 
 class TopologyCreator {
 public:
-  static TopologyPtr Create(const ost::mol::EntityHandle& ent, 
+  static TopologyPtr Create(ost::mol::EntityHandle& ent, 
                             const MMSettingsPtr settings);
 
 };
diff --git a/modules/mol/mm/tests/test_simulation.cc b/modules/mol/mm/tests/test_simulation.cc
index 6d1351972639d9042dc8cf5e0b2c85fa9f600432..6a10eee3b4c1c257f672e30a79530f4bc4840ade 100644
--- a/modules/mol/mm/tests/test_simulation.cc
+++ b/modules/mol/mm/tests/test_simulation.cc
@@ -46,7 +46,7 @@ BOOST_AUTO_TEST_CASE(test_simulation_basics){
 
   settings->integrator = IntegratorPtr(new OpenMM::VerletIntegrator(0.002));
 
-  Simulation sim(top, settings);
+  Simulation sim(top, test_ent, settings);
 
 
   //we check, wether the reset functions have the desired effect on the topology
@@ -134,10 +134,10 @@ BOOST_AUTO_TEST_CASE(test_simulation_energy_calculations){
 
   settings->integrator = IntegratorPtr(new OpenMM::VerletIntegrator(0.002));
 
-  Simulation sim(top, settings);
+  Simulation sim(top, test_ent,settings);
 
   sim.MinimizeEnergy("steep",1.0,200);
-  sim.UpdateTopologyPositions();
+  sim.UpdatePositions();
 
   ost::mol::EntityHandle minimized_ent = sim.GetEntity();
 
diff --git a/modules/mol/mm/tests/test_topology.cc b/modules/mol/mm/tests/test_topology.cc
index 52744ed6cdb261fac98c91a3fc5b0c1f28b64079..59de803c02a02aabe90b635a881c6b4cc7be6461 100644
--- a/modules/mol/mm/tests/test_topology.cc
+++ b/modules/mol/mm/tests/test_topology.cc
@@ -30,9 +30,7 @@ BOOST_AUTO_TEST_CASE(test_topology_basics){
   std::vector<Real> zero_real(0);
   std::vector<Real> lot_of_reals(test_ent.GetAtomCount());
 
-  BOOST_CHECK_THROW(Topology(test_ent,zero_real),ost::Error);
-
-  TopologyPtr top(new Topology(test_ent,lot_of_reals));
+  TopologyPtr top(new Topology(lot_of_reals));
 
   uint ui1(0), ui2(0), ui3(0), ui4(0), ui5(0);
   int i1(0), i5(0);  
@@ -283,56 +281,19 @@ BOOST_AUTO_TEST_CASE(test_topology_basics){
   BOOST_CHECK(top->GetHarmonicDistanceRestraintIndices(1).size() == 10);
   BOOST_CHECK(top->GetHarmonicDistanceRestraintIndices(2).size() == 0);
 
-  ost::mol::EntityHandle top_test_ent = top->GetEntity();
-  ost::mol::AtomHandleList atom_list = top_test_ent.GetAtomList();
-  top->AddPositionConstraint(atom_list[0]);
+  top->AddPositionConstraint(5);
   BOOST_CHECK(top->GetNumPositionConstraints() == 1);
   top->ResetPositionConstraints();
   BOOST_CHECK(top->GetNumPositionConstraints() == 0);
 }
 
-BOOST_AUTO_TEST_CASE(test_topology_index_getters){
-
-  String pdb_name = "1CRN.pdb";
-  ost::io::PDBReader reader(pdb_name, ost::io::IOProfile());
-  ost::mol::EntityHandle test_ent = ost::mol::CreateEntity();
-  reader.Import(test_ent);
-  ost::conop::ProcessorPtr processor(new ost::conop::HeuristicProcessor);
-  processor->Process(test_ent);
-
-  //check initialisation without settings
-  std::vector<Real> fake_masses(test_ent.GetAtomCount());
-  TopologyPtr top(new Topology(test_ent,fake_masses));
-  test_ent = top->GetEntity(); //we do a copy in the constructor of the topology
-
-  ost::mol::AtomHandleList atom_list = test_ent.GetAtomList();
-  ost::mol::ResidueHandleList res_list = test_ent.GetResidueList();
-
-  for(uint i = 0; i < atom_list.size(); ++i){
-    BOOST_CHECK(top->GetAtomIndex(atom_list[i]) == i);
-  }
-  for(uint i = 0; i < res_list.size(); ++i){
-    BOOST_CHECK(top->GetResidueIndex(res_list[i]) == i);
-  }
-
-  ost::mol::AtomHandle test_atom_one = res_list[10].FindAtom("CA");
-  ost::mol::AtomHandle test_atom_two = res_list[9].FindAtom("CA");
-  ost::mol::AtomHandle test_atom_three = res_list[11].FindAtom("CA");
-
-  uint index_one = top->GetAtomIndex(test_atom_one);
-  uint index_two = top->GetAtomIndex(test_atom_two);
-  uint index_three = top->GetAtomIndex(test_atom_three);
-
-  BOOST_CHECK(top->GetAtomIndex(10,"CA") == index_one);
-  BOOST_CHECK(top->GetAtomIndex(10,"-CA") == index_two);
-  BOOST_CHECK(top->GetAtomIndex(10,"+CA") == index_three);
-}
-
 BOOST_AUTO_TEST_CASE(test_topology_merge){
 
   String pdb_name = "1CRN.pdb";
   ost::io::PDBReader reader(pdb_name, ost::io::IOProfile());
   ost::mol::EntityHandle test_ent = ost::mol::CreateEntity();
+  ost::mol::EntityHandle temp;
+
   reader.Import(test_ent);
 
   ForcefieldPtr ff = Forcefield::Load("CHARMM27.dat");
@@ -354,7 +315,7 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
 
   //check whether error gets thrown when chain with same name is already present
   std::vector<Real> real_vec(6);
-  TopologyPtr merge_top_one(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_one(new Topology(real_vec));
   merge_top_one->SetCharges(real_vec);
   merge_top_one->SetOBCScalings(real_vec);
   merge_top_one->SetGBSARadii(real_vec);
@@ -363,10 +324,11 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_one->SetFudgeQQ(1.0);
   merge_top_one->SetFudgeLJ(1.0);
   TopologyPtr top_one = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_one->Merge(merge_top_one),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_one->Merge(temp,merge_top_one,new_ent),ost::Error); 
   ost::mol::XCSEditor ed = new_ent.EditXCS();
   ed.RenameChain(new_ent.GetChainList()[0],"B");
-  TopologyPtr merge_top_one_two(new Topology(new_ent,real_vec)); // I know...
+  TopologyPtr merge_top_one_two(new Topology(real_vec)); // I know...
   merge_top_one_two->SetCharges(real_vec);
   merge_top_one_two->SetOBCScalings(real_vec);
   merge_top_one_two->SetGBSARadii(real_vec);
@@ -374,11 +336,11 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_one_two->SetEpsilons(real_vec);
   merge_top_one_two->SetFudgeQQ(1.0);
   merge_top_one_two->SetFudgeLJ(1.0);
-  BOOST_CHECK_NO_THROW(top_one->Merge(merge_top_one_two));
+  BOOST_CHECK_NO_THROW(top_one->Merge(temp,merge_top_one_two,new_ent));
 
 
   //check whether error gets thrown, when fudge parameters are inconsistent
-  TopologyPtr merge_top_two(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_two(new Topology(real_vec));
   merge_top_two->SetCharges(real_vec);
   merge_top_two->SetOBCScalings(real_vec);
   merge_top_two->SetGBSARadii(real_vec);
@@ -387,14 +349,15 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_two->SetFudgeQQ(42.0);
   merge_top_two->SetFudgeLJ(42.0);
   TopologyPtr top_two = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_two->Merge(merge_top_two),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_two->Merge(temp,merge_top_two,new_ent),ost::Error); 
   merge_top_two->SetFudgeQQ(1.0);
   merge_top_two->SetFudgeLJ(1.0);
-  BOOST_CHECK_NO_THROW(top_two->Merge(merge_top_two));
+  BOOST_CHECK_NO_THROW(top_two->Merge(temp,merge_top_two,new_ent));
 
 
   //check whether error gets thrown when charges are not set
-  TopologyPtr merge_top_three(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_three(new Topology(real_vec));
   merge_top_three->SetOBCScalings(real_vec);
   merge_top_three->SetGBSARadii(real_vec);
   merge_top_three->SetSigmas(real_vec);
@@ -402,13 +365,14 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_three->SetFudgeQQ(1.0);
   merge_top_three->SetFudgeLJ(1.0);
   TopologyPtr top_three = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_three->Merge(merge_top_three),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_three->Merge(temp,merge_top_three,new_ent),ost::Error); 
   TopologyPtr top_three_two = TopologyCreator::Create(test_ent,settings);
   merge_top_three->SetCharges(real_vec);
-  BOOST_CHECK_NO_THROW(top_three_two->Merge(merge_top_three));
+  BOOST_CHECK_NO_THROW(top_three_two->Merge(temp,merge_top_three,new_ent));
 
   //check whether error gets thrown when obc scaling factors are not set
-  TopologyPtr merge_top_four(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_four(new Topology(real_vec));
   merge_top_four->SetCharges(real_vec);
   merge_top_four->SetGBSARadii(real_vec);
   merge_top_four->SetSigmas(real_vec);
@@ -416,13 +380,14 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_four->SetFudgeQQ(1.0);
   merge_top_four->SetFudgeLJ(1.0);
   TopologyPtr top_four = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_four->Merge(merge_top_four),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_four->Merge(temp,merge_top_four,new_ent),ost::Error); 
   TopologyPtr top_four_two = TopologyCreator::Create(test_ent,settings);
   merge_top_four->SetOBCScalings(real_vec);
-  BOOST_CHECK_NO_THROW(top_four_two->Merge(merge_top_four));
+  BOOST_CHECK_NO_THROW(top_four_two->Merge(temp,merge_top_four,new_ent));
 
   //check whether error gets thrown when gbsa radii are not set
-  TopologyPtr merge_top_five(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_five(new Topology(real_vec));
   merge_top_five->SetCharges(real_vec);
   merge_top_five->SetOBCScalings(real_vec);
   merge_top_five->SetSigmas(real_vec);
@@ -430,13 +395,14 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_five->SetFudgeQQ(1.0);
   merge_top_five->SetFudgeLJ(1.0);
   TopologyPtr top_five = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_five->Merge(merge_top_five),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_five->Merge(temp,merge_top_five,new_ent),ost::Error); 
   TopologyPtr top_five_two = TopologyCreator::Create(test_ent,settings);
   merge_top_five->SetGBSARadii(real_vec);
-  BOOST_CHECK_NO_THROW(top_five_two->Merge(merge_top_five));
+  BOOST_CHECK_NO_THROW(top_five_two->Merge(temp,merge_top_five,new_ent));
 
   //check whether error gets thrown when sigmas are not set
-  TopologyPtr merge_top_six(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_six(new Topology(real_vec));
   merge_top_six->SetCharges(real_vec);
   merge_top_six->SetOBCScalings(real_vec);
   merge_top_six->SetGBSARadii(real_vec);  
@@ -444,13 +410,14 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_six->SetFudgeQQ(1.0);
   merge_top_six->SetFudgeLJ(1.0);
   TopologyPtr top_six = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_six->Merge(merge_top_six),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_six->Merge(temp,merge_top_six,new_ent),ost::Error); 
   TopologyPtr top_six_two = TopologyCreator::Create(test_ent,settings);
   merge_top_six->SetSigmas(real_vec);
-  BOOST_CHECK_NO_THROW(top_six_two->Merge(merge_top_six));
+  BOOST_CHECK_NO_THROW(top_six_two->Merge(temp,merge_top_six,new_ent));
 
   //check whether error gets thrown when epsilons are not set
-  TopologyPtr merge_top_seven(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_seven(new Topology(real_vec));
   merge_top_seven->SetCharges(real_vec);
   merge_top_seven->SetOBCScalings(real_vec);
   merge_top_seven->SetGBSARadii(real_vec); 
@@ -458,12 +425,13 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_seven->SetFudgeQQ(1.0);
   merge_top_seven->SetFudgeLJ(1.0);
   TopologyPtr top_seven = TopologyCreator::Create(test_ent,settings);
-  BOOST_CHECK_THROW(top_seven->Merge(merge_top_seven),ost::Error); 
+  temp = test_ent.Copy();
+  BOOST_CHECK_THROW(top_seven->Merge(temp,merge_top_seven,new_ent),ost::Error); 
   TopologyPtr top_seven_two = TopologyCreator::Create(test_ent,settings);
   merge_top_seven->SetEpsilons(real_vec);
-  BOOST_CHECK_NO_THROW(top_seven_two->Merge(merge_top_seven));
+  BOOST_CHECK_NO_THROW(top_seven_two->Merge(temp,merge_top_seven,new_ent));
 
-  TopologyPtr merge_top_eight(new Topology(new_ent,real_vec));
+  TopologyPtr merge_top_eight(new Topology(real_vec));
   merge_top_eight->SetCharges(real_vec);
   merge_top_eight->SetOBCScalings(real_vec);
   merge_top_eight->SetGBSARadii(real_vec); 
@@ -473,7 +441,6 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_eight->SetEpsilons(real_vec);
 
   //let's add every possible interaction
-
   merge_top_eight->AddHarmonicBond(0,1,42.0,42000.0);
   merge_top_eight->AddHarmonicAngle(0,1,2,24.0,24000.0);
   merge_top_eight->AddUreyBradleyAngle(1,2,3,1.0,2.0,3.0,4.0);
@@ -494,7 +461,7 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   merge_top_eight->AddHarmonicPositionRestraint(2,pos,10.0,5.0,6.0,7.0);
   merge_top_eight->AddHarmonicDistanceRestraint(1,2,10.0,1000.0);
 
-  TopologyPtr top_eight = TopologyCreator::Create(test_ent,settings);\
+  TopologyPtr top_eight = TopologyCreator::Create(test_ent,settings);
 
   uint num_harmonic_bonds = top_eight->GetNumHarmonicBonds();
   uint num_harmonic_angles = top_eight->GetNumHarmonicAngles();
@@ -510,7 +477,8 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   uint num_harmonic_position_restraints = top_eight->GetNumHarmonicPositionRestraints();
   uint num_harmonic_distance_restraints = top_eight->GetNumHarmonicDistanceRestraints();
 
-  top_eight->Merge(merge_top_eight);
+  temp = test_ent.Copy();
+  top_eight->Merge(temp, merge_top_eight, new_ent);
 
   BOOST_CHECK(top_eight->GetNumHarmonicBonds() == num_harmonic_bonds+1);
   BOOST_CHECK(top_eight->GetNumHarmonicAngles() == num_harmonic_angles+1);
@@ -526,28 +494,12 @@ BOOST_AUTO_TEST_CASE(test_topology_merge){
   BOOST_CHECK(top_eight->GetNumHarmonicPositionRestraints() == num_harmonic_position_restraints+1);
   BOOST_CHECK(top_eight->GetNumHarmonicDistanceRestraints() == num_harmonic_distance_restraints+1);
 
-
-  //check wether the new indices are right...
-  ost::mol::EntityHandle top_ent = top_eight->GetEntity();
-  ost::mol::ChainHandle top_chain = top_ent.FindChain("B");
-  ost::mol::ResidueHandleList res_list = top_chain.GetResidueList();
-  ost::mol::AtomHandleList atom_list = top_chain.GetAtomList();
-  uint res_index_one = top_eight->GetResidueIndex(res_list[0]);
-  uint res_index_two = top_eight->GetResidueIndex(res_list[1]);
-  uint res_index_three = top_eight->GetResidueIndex(res_list[2]);
-  uint atom_index_zero = top_eight->GetAtomIndex(res_index_one,"A");
-  uint atom_index_one = top_eight->GetAtomIndex(res_index_one,"B");
-  uint atom_index_two = top_eight->GetAtomIndex(res_index_two,"A");
-  uint atom_index_three = top_eight->GetAtomIndex(res_index_two,"B");
-  uint atom_index_four = top_eight->GetAtomIndex(res_index_three,"A");
-  uint atom_index_five = top_eight->GetAtomIndex(res_index_three,"B");
-
-  BOOST_CHECK(atom_index_zero == top_eight->GetAtomIndex(atom_list[0]));
-  BOOST_CHECK(atom_index_one == top_eight->GetAtomIndex(atom_list[1]));
-  BOOST_CHECK(atom_index_two == top_eight->GetAtomIndex(atom_list[2]));
-  BOOST_CHECK(atom_index_three == top_eight->GetAtomIndex(atom_list[3]));
-  BOOST_CHECK(atom_index_four == top_eight->GetAtomIndex(atom_list[4]));
-  BOOST_CHECK(atom_index_five == top_eight->GetAtomIndex(atom_list[5]));
+  uint atom_index_zero = 0 + test_ent.GetAtomCount();
+  uint atom_index_one = 1 + test_ent.GetAtomCount();
+  uint atom_index_two = 2 + test_ent.GetAtomCount();
+  uint atom_index_three = 3 + test_ent.GetAtomCount();
+  uint atom_index_four = 4 + test_ent.GetAtomCount();
+  //uint atom_index_five = 5 + test_ent.GetAtomCount();
 
   //check whether the unique interactions are properly mapped
   BOOST_CHECK_THROW(top_eight->AddExclusion(atom_index_three,atom_index_four),ost::Error);