diff --git a/src/montecarlo.cpp b/src/montecarlo.cpp
index 73e28a6f9..62f83a2ed 100644
--- a/src/montecarlo.cpp
+++ b/src/montecarlo.cpp
@@ -225,6 +225,10 @@ void MetropolisMonteCarlo::State::sync(const State& other, const Change& change)
     pot->sync(&*other.pot, change);
 }
 
+Space& MetropolisMonteCarlo::State::getSpace() { return *spc; }
+
+Energy::Hamiltonian& MetropolisMonteCarlo::State::getHamiltonian() { return *pot; }
+
 void to_json(json& j, const MetropolisMonteCarlo& monte_carlo) {
     j = monte_carlo.state->spc->info();
     j["temperature"] = pc::temperature / 1.0_K;
@@ -238,6 +242,8 @@ void to_json(json& j, const MetropolisMonteCarlo& monte_carlo) {
                            {"last move", monte_carlo.latest_move_name}};
     }
 }
+MetropolisMonteCarlo::State& MetropolisMonteCarlo::getTrialState() { return *state; }
+MetropolisMonteCarlo::State& MetropolisMonteCarlo::getState() { return *trial_state; }
 
 TranslationalEntropy::TranslationalEntropy(Space &trial_space, Space &space) : trial_spc(trial_space), spc(space) {}
 
diff --git a/src/montecarlo.h b/src/montecarlo.h
index 100562961..3742dea31 100644
--- a/src/montecarlo.h
+++ b/src/montecarlo.h
@@ -53,6 +53,8 @@ class MetropolisMonteCarlo {
         std::unique_ptr<Energy::Hamiltonian> pot; //!< Hamiltonian for calc. potential energy
         void sync(const State& other,
                   const Change& change); //!< Sync with another state (the other state is not modified)
+        Space& getSpace();
+        Energy::Hamiltonian& getHamiltonian();
     };
 
   private:
@@ -79,6 +81,8 @@ class MetropolisMonteCarlo {
     void restore(const json& j);                               //!< Restores system from previously store json object
     static bool metropolisCriterion(double energy_change);     //!< Metropolis criterion
     ~MetropolisMonteCarlo();                                   //!< Required due to unique_ptr to incomplete type
+    State& getTrialState();                                    //!< Get trial state where trial moves are performed
+    State& getState();                                         //!< Get current accepted state
 };
 
 void from_json(const json &, MetropolisMonteCarlo::State &); //!< Build state from json object
diff --git a/src/pyfaunus.cpp b/src/pyfaunus.cpp
index 9b0d7f91e..4b6363edb 100644
--- a/src/pyfaunus.cpp
+++ b/src/pyfaunus.cpp
@@ -17,9 +17,7 @@
 namespace py = pybind11;
 using namespace Faunus;
 
-typedef typename Space::GroupType Tgroup;
-typedef Energy::Hamiltonian Thamiltonian;
-typedef MetropolisMonteCarlo Tmcsimulation;
+// @todo The code below should be updated, see https://github.com/pybind/pybind11_json
 
 inline json dict2json(py::dict dict) {
     py::object dumps = py::module::import("json").attr("dumps");
@@ -38,9 +36,7 @@ inline py::dict json2dict(const json &j) {
 
 template<class T>
 std::unique_ptr<T> from_dict(py::dict dict) {
-    auto ptr = new T();
-    *ptr = dict2json(dict);
-    return std::unique_ptr<T>(ptr);
+    return std::make_unique<T>(dict2json(dict));
 } // convert py::dict to T through Faunus::json
 
 PYBIND11_MODULE(pyfaunus, m)
@@ -49,14 +45,12 @@ PYBIND11_MODULE(pyfaunus, m)
 
     // json
     py::class_<json>(m, "json")
-        .def(py::init( [](std::string arg) {
-                    return std::unique_ptr<json>(new json(json::parse(arg)));
-                    } ) )
-    .def(py::init([](py::dict dict) {
-                py::object dumps = py::module::import("json").attr("dumps");
-                std::string s = dumps(dict).cast<std::string>();
-                return std::unique_ptr<json>(new json(json::parse(s)));
-                } ) );
+        .def(py::init([](std::string arg) { return std::make_unique<json>(json(json::parse(arg))); }))
+        .def(py::init([](py::dict dict) {
+            py::object dumps = py::module::import("json").attr("dumps");
+            auto s = dumps(dict).cast<std::string>();
+            return std::make_unique<json>(json(json::parse(s)));
+        }));
 
     // Random
     py::class_<Random>(m, "Random")
@@ -139,25 +133,24 @@ PYBIND11_MODULE(pyfaunus, m)
         .def("end", [](ParticleVector& p) { return p.end(); });
 
     // Group
-    py::class_<Tgroup>(m, "Group")
+    py::class_<Group>(m, "Group")
         .def(py::init<MoleculeData::index_type, ParticleVector::iterator, ParticleVector::iterator>())
-        .def_readwrite("groups", &Tgroup::id, "Molecule id")
-        .def_readwrite("id", &Tgroup::id, "Molecule id")
-        .def_readwrite("cm", &Tgroup::mass_center, "Center of mass")
-        .def("__len__", [](Tgroup& self) { return self.size(); })
+        .def_readwrite("id", &Group::id, "Molecule id")
+        .def_readwrite("mass_center", &Group::mass_center, "Center of mass")
+        .def("__len__", [](Group& self) { return self.size(); })
         .def(
-            "__iter__", [](Tgroup& v) { return py::make_iterator(v.begin(), v.end()); }, py::keep_alive<0, 1>())
-        .def("isAtomic", &Tgroup::isAtomic)
-        .def("isMolecular", &Tgroup::isMolecular)
-        .def("traits", &Tgroup::traits)
-        .def("contains", &Tgroup::contains)
-        .def("capacity", &Tgroup::capacity)
-        .def("deactivate", &Tgroup::deactivate)
-        .def("activate", &Tgroup::activate)
-        .def("begin", (ParticleVector::iterator & (Tgroup::*)()) & Tgroup::begin)
-        .def("end", (ParticleVector::iterator & (Tgroup::*)()) & Tgroup::end);
-
-    py::bind_vector<std::vector<Tgroup>>(m, "GroupVector");
+            "__iter__", [](Group& v) { return py::make_iterator(v.begin(), v.end()); }, py::keep_alive<0, 1>())
+        .def("isAtomic", &Group::isAtomic)
+        .def("isMolecular", &Group::isMolecular)
+        .def("traits", &Group::traits)
+        .def("contains", &Group::contains)
+        .def("capacity", &Group::capacity)
+        .def("deactivate", &Group::deactivate)
+        .def("activate", &Group::activate)
+        .def("begin", (ParticleVector::iterator & (Group::*)()) & Group::begin)
+        .def("end", (ParticleVector::iterator & (Group::*)()) & Group::end);
+
+    py::bind_vector<std::vector<Group>>(m, "GroupVector");
 
     // Region
     py::enum_<Region::RegionBase::RegionType>(m, "RegionType")
@@ -194,6 +187,15 @@ PYBIND11_MODULE(pyfaunus, m)
 
     m.attr("atoms") = &Faunus::atoms; // global instance
 
+    // MoleculeData
+    py::class_<MoleculeData>(m, "MoleculeData")
+        .def(py::init<>())
+        .def_readwrite("name", &MoleculeData::name)
+        .def("numConformations", &MoleculeData::numConformations);
+
+    auto _moleculedatavec = py::bind_vector<decltype(Faunus::molecules)>(m, "MoleculeDataVector");
+    m.attr("molecules") = &Faunus::molecules; // global instance
+
     // Temperature and other globals etc.
     m.def("getTemperature", []() { return pc::temperature; } );
     m.def("setTemperature", [](double T) { pc::temperature = T; } );
@@ -244,31 +246,52 @@ PYBIND11_MODULE(pyfaunus, m)
         // .def("findMolecules", &Space::findMolecules)
         .def("from_dict", [](Space& spc, py::dict dict) { from_json(dict2json(dict), spc); });
 
+    // Energybase
+    py::class_<Energy::Energybase>(m, "Energybase")
+        .def_readonly("name", &Energy::Energybase::name)
+        .def("force", &Energy::Energybase::force)
+        .def("energy", &Energy::Energybase::energy)
+        .def("sync", &Energy::Energybase::sync);
+
     // Hamiltonian
-    py::class_<Thamiltonian>(m, "Hamiltonian")
-        .def(py::init<Space &, const json &>())
-        .def(py::init([](Space &spc, py::list list) {
+    py::class_<Energy::Hamiltonian>(m, "Hamiltonian")
+        .def(py::init<Space&, const json&>())
+        .def(py::init([](Space& spc, py::list list) {
             json j = list2json(list);
-            return std::unique_ptr<Thamiltonian>(new Thamiltonian(spc, j));
+            return std::make_unique<Energy::Hamiltonian>(spc, j);
         }))
-        .def("init", &Thamiltonian::init)
-        .def("energy", &Thamiltonian::energy);
+        .def("init", &Energy::Hamiltonian::init)
+        .def("energy", &Energy::Hamiltonian::energy)
+        .def(
+            "__iter__",
+            [](Energy::Hamiltonian& hamiltonian) { return py::make_iterator(hamiltonian.begin(), hamiltonian.end()); },
+            py::keep_alive<0, 1>());
+
+    // State
+    py::class_<MetropolisMonteCarlo::State>(m, "State")
+        .def("getSpace", &MetropolisMonteCarlo::State::getSpace)
+        .def("getHamiltonian", &MetropolisMonteCarlo::State::getHamiltonian)
+        .def("sync", &MetropolisMonteCarlo::State::sync);
 
     // TranslationalEntropy
     py::class_<TranslationalEntropy>(m, "TranslationalEntropy")
         .def(py::init<Space &, Space &>())
         .def("energy", &TranslationalEntropy::energy);
 
-    // MCSimulation
-    py::class_<Tmcsimulation>(m, "MetropolisMonteCarlo")
+    // MetropolisMonteCarlo
+    py::class_<MetropolisMonteCarlo>(m, "MetropolisMonteCarlo")
         .def(py::init([](py::dict dict) {
             json j = dict2json(dict);
-            return std::unique_ptr<Tmcsimulation>(new Tmcsimulation(j));
+            return std::make_unique<MetropolisMonteCarlo>(j);
         }))
         .def(py::init([](py::dict dict) {
             json j = dict2json(dict);
-            return std::unique_ptr<Tmcsimulation>(new Tmcsimulation(j));
-        }));
+            return std::make_unique<MetropolisMonteCarlo>(j);
+        }))
+        .def("sweep", &MetropolisMonteCarlo::sweep)
+        .def("getState", &MetropolisMonteCarlo::getState)
+        .def("getTrialState", &MetropolisMonteCarlo::getTrialState)
+        .def("relativeEnergyDrift", &MetropolisMonteCarlo::relativeEnergyDrift);
 
     // Analysisbase
     py::class_<Analysis::Analysisbase>(m, "Analysisbase")
@@ -286,16 +309,17 @@ PYBIND11_MODULE(pyfaunus, m)
 
     // CombinedAnalysis
     py::class_<Analysis::CombinedAnalysis>(m, "Analysis")
-        .def(py::init([](Space &spc, Thamiltonian &pot, py::list list) {
+        .def(py::init([](Space& spc, Energy::Hamiltonian& pot, py::list list) {
             json j = list2json(list);
-            return std::unique_ptr<Analysis::CombinedAnalysis>(new Analysis::CombinedAnalysis(j, spc, pot));
+            return std::make_unique<Analysis::CombinedAnalysis>(j, spc, pot);
         }))
-        .def_readwrite("vector", &Analysis::CombinedAnalysis::vec)
+        .def_readwrite("moves", &Analysis::CombinedAnalysis::vec)
         .def("to_dict",
-             [](Analysis::CombinedAnalysis &self) {
+             [](Analysis::CombinedAnalysis& self) {
                  json j;
                  Faunus::to_json(j, self);
                  return json2dict(j);
              })
-        .def("sample", &Analysis::CombinedAnalysis::sample);
+        .def("sample", &Analysis::CombinedAnalysis::sample)
+        .def("to_disk", &Analysis::CombinedAnalysis::to_disk);
 }