diff --git a/modules/mol/mm/examples/gb_example.py b/modules/mol/mm/examples/gb_example.py
index 61a563e2e8d7621d147832e94c11a38943266cef..60a5ddf8934ecbbf29b01dc8d3623be88453d0c4 100644
--- a/modules/mol/mm/examples/gb_example.py
+++ b/modules/mol/mm/examples/gb_example.py
@@ -28,7 +28,7 @@ settings.integrator = mm.LangevinIntegrator(310,1,0.002)
 settings.constrain_bonds = True
 settings.constrain_hbonds = True
 settings.constrain_hangles = True
-#settings.add_gbsa = True
+settings.add_gbsa = True
 #settings.add_thermostat = True
 #settings.thermostat_temperature = 310.0
 #settings.thermostat_collision_frequency = 1
diff --git a/modules/mol/mm/pymod/export_simulation.cc b/modules/mol/mm/pymod/export_simulation.cc
index 2a8f4436245dd7015bb928814ad740af8f7fa3a8..91f5e25b95362210e39788654b4cf7dee2fbbc6d 100644
--- a/modules/mol/mm/pymod/export_simulation.cc
+++ b/modules/mol/mm/pymod/export_simulation.cc
@@ -50,6 +50,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("Save",&ost::mol::mm::Simulation::Save)
+    .def("Load",&ost::mol::mm::Simulation::Load).staticmethod("Load")
     .def("Steps",&ost::mol::mm::Simulation::Steps)
     .def("GetPositions",&ost::mol::mm::Simulation::GetPositions,(arg("enforce_periodic_box")=false, arg("in_angstrom")=true))
     .def("GetVelocities",&ost::mol::mm::Simulation::GetVelocities)
diff --git a/modules/mol/mm/src/simulation.cc b/modules/mol/mm/src/simulation.cc
index 213bd4385a5439c0704ab72fcdb82f91945f5790..8433735c6af02e6bb757ffcb66fb0ef119ea0716 100644
--- a/modules/mol/mm/src/simulation.cc
+++ b/modules/mol/mm/src/simulation.cc
@@ -18,6 +18,91 @@ Simulation::Simulation(const TopologyPtr top,
   this->Init(top, settings);
 }
 
+void Simulation::Save(const String& filename){
+  std::ofstream stream(filename.c_str(), std::ios_base::binary);
+  io::BinaryDataSink ds(stream);
+  ds << *top_;
+  geom::Vec3List positions = this->GetPositions(false,false);
+  for(geom::Vec3List::iterator i = positions.begin(); 
+      i != positions.end(); ++i){
+    ds & (*i)[0];
+    ds & (*i)[1];
+    ds & (*i)[2];
+  }
+  context_->createCheckpoint(stream);  
+}
+
+SimulationPtr Simulation::Load(const String& filename, MMSettingsPtr settings){
+  if (!boost::filesystem::exists(filename)) {
+    std::stringstream ss;
+    ss << "Could not open topology. File '"
+       << filename << "' does not exist";
+    throw ost::io::IOException(ss.str());
+  }
+
+  SimulationPtr sim_ptr(new Simulation);
+
+  std::ifstream stream(filename.c_str(), std::ios_base::binary);
+  io::BinaryDataSource ds(stream);
+  TopologyPtr top_p(new Topology);
+  ds >> *top_p;
+
+  sim_ptr->original_masses_ = top_p->GetMasses(); 
+
+  sim_ptr->top_ = top_p;
+
+  sim_ptr->system_ = SystemCreator::Create(sim_ptr->top_,settings,
+                                       sim_ptr->system_force_mapper_);
+
+  sim_ptr->integrator_ = settings->integrator;
+
+  OpenMM::Platform::loadPluginsFromDirectory (settings->openmm_plugin_directory);
+  OpenMM::Platform* platform;
+
+  switch(settings->platform){
+    case Reference:{
+      platform = &OpenMM::Platform::getPlatformByName("Reference");
+      break;
+    }
+    case OpenCL:{
+      platform = &OpenMM::Platform::getPlatformByName("OpenCL");
+      break;
+    }
+    case CUDA:{
+      platform = &OpenMM::Platform::getPlatformByName("CUDA");
+      break;
+    }
+    case CPU:{
+      platform = &OpenMM::Platform::getPlatformByName("CPU");
+      break;
+    }
+  }
+
+  sim_ptr->context_ = ContextPtr(new OpenMM::Context(*(sim_ptr->system_),
+                                                     *(sim_ptr->integrator_),
+                                                     *platform));
+
+  std::vector<OpenMM::Vec3> positions;
+  OpenMM::Vec3 open_mm_vec;
+  Real a,b,c;
+  for(int i = 0; i < sim_ptr->system_->getNumParticles(); ++i){
+    ds & a;
+    ds & b;
+    ds & c;
+    open_mm_vec[0] = a;
+    open_mm_vec[1] = b;
+    open_mm_vec[2] = c;
+    positions.push_back(open_mm_vec);
+  }
+  sim_ptr->context_->setPositions(positions);
+
+  sim_ptr->context_->loadCheckpoint(stream);
+
+  return sim_ptr;
+}
+
+
+
 void Simulation::Init(const TopologyPtr top,
                       const MMSettingsPtr settings){
 
@@ -576,7 +661,5 @@ void Simulation::SetPeriodicBoxExtents(geom::Vec3& vec){
   context_->setPeriodicBoxVectors(ucell_a,ucell_b,ucell_c);
 }
 
-
-
 }}}
 
diff --git a/modules/mol/mm/src/simulation.hh b/modules/mol/mm/src/simulation.hh
index a935c63badc5c6ba2b011cdf311dc3cf4c8d88d8..6b7ba1370621a493f411b4284b8e0aefea64ebd3 100644
--- a/modules/mol/mm/src/simulation.hh
+++ b/modules/mol/mm/src/simulation.hh
@@ -47,6 +47,9 @@ public:
   Simulation(const ost::mol::mm::TopologyPtr top,
              const MMSettingsPtr settings);
 
+  void Save(const String& filename);
+
+  static SimulationPtr Load(const String& filename, MMSettingsPtr settings);
 
   ost::mol::EntityHandle GetEntity() { return top_->GetEntity(); }
 
@@ -116,6 +119,9 @@ public:
   void SetPeriodicBoxExtents(geom::Vec3& vec);
 
 private:
+
+  Simulation() { } //hidden constructor... 
+
   void Init(const ost::mol::mm::TopologyPtr top, 
             const MMSettingsPtr settings);
 
@@ -127,7 +133,7 @@ private:
   TopologyPtr top_;
   std::vector<MMObserverPtr> observers_;
   std::vector<int> time_to_notify_;
-  std::map<FuncType,unsigned int> system_force_mapper_;
+  std::map<FuncType,uint> system_force_mapper_;
   std::vector<Real> original_masses_;
 };
 
diff --git a/modules/mol/mm/src/topology.hh b/modules/mol/mm/src/topology.hh
index 16984c7db7c7c15080137ddaca1c8a5c65bc4637..93284524c65bf24e522a4ea0df1a71cdc622f9b7 100644
--- a/modules/mol/mm/src/topology.hh
+++ b/modules/mol/mm/src/topology.hh
@@ -38,6 +38,9 @@ public:
 
   Topology(const ost::mol::EntityHandle& ent, const std::vector<Real>& masses);
 
+  Topology() { } //should not be accessible from Python to avoid messing around
+                 //with empty topology
+
   static TopologyPtr Load(const String& filename);
 
   void Save(const String& filename);
@@ -788,8 +791,6 @@ public:
 
 private:
 
-  Topology() { } //hidden constructor without parameters
-
   void InitMappers();
 
   ost::mol::EntityHandle ent_;