diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b0fdfb2e4def83a7750cb6cf4404c37894bd2cf4..7548629ca3052085eedf5023368683494d7c48d2 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,6 +14,12 @@ target_include_directories(mysimpleheisenbergthreaded PUBLIC ${MYHDF5INCLUDES})
 target_link_libraries(mysimpleheisenbergthreaded libmarqov)
 set_target_properties(mysimpleheisenbergthreaded PROPERTIES CXX_STANDARD 14 CXX_STANDARD_REQUIRED YES)
 
+# configuration for the Klein-Gordon case
+add_executable(kleingordon kleingordon.cpp libmarqov/marqov.cpp)
+target_include_directories(kleingordon PUBLIC ${MYHDF5INCLUDES})
+target_link_libraries(kleingordon libmarqov)
+set_target_properties(kleingordon PROPERTIES CXX_STANDARD 14 CXX_STANDARD_REQUIRED YES)
+
 # configuration of the main parts.
 
 SET(MARQOVTARGETS MARQOV mpiMARQOV)
diff --git a/src/config.h b/src/config.h
new file mode 100644
index 0000000000000000000000000000000000000000..f43a11165def1962f664a3c8f8d42aedcfcd9dcf
--- /dev/null
+++ b/src/config.h
@@ -0,0 +1,3 @@
+#define HAVE_SCANDIR
+#define HAVE_ALPHASORT
+#define HAVE_READDIR_R
diff --git a/src/hamiltonian/MassiveScalarField.h b/src/hamiltonian/MassiveScalarField.h
new file mode 100644
index 0000000000000000000000000000000000000000..a7d957eaf899e7d6ec5100eb22aba3fcd9ebfc25
--- /dev/null
+++ b/src/hamiltonian/MassiveScalarField.h
@@ -0,0 +1,300 @@
+#ifndef MASSIVESCALARFIELD_H
+#define MASSIVESCALARFIELD_H
+#include <array>
+#include <cmath>
+#include <vector>
+#include "util/randomdir.h"
+#include "util/hamparts.h"
+#include "util/termcollection.h" 
+
+
+auto hyperfilter = [](auto p)
+{
+	auto& mp = std::get<1>(p);
+	auto& hp = std::get<2>(p);
+
+	std::string str_repid = std::to_string(mp.repid);
+	std::string str_mass  = "mass"+std::to_string(std::get<3>(hp));
+	mp.outname = str_mass + "_" + str_repid;
+
+	return p;
+};
+
+
+// ------------------------------ OBSERVABLES ---------------------------
+
+#include "util/observables.h"
+
+class BulkMagnetization
+{
+	public:
+		std::string name; ///< The name of the observable
+        std::string desc; ///< A helpful description that will be used in the HDF5 output files.
+
+        /** Construct a bulk magnetization object
+         * 
+         */
+		BulkMagnetization() : name("mb"), desc("bulk magnetization") {}
+
+		template <class StateSpace, class Grid>
+		double measure(const StateSpace& statespace, const Grid& grid)
+		{
+			double mag = 0;
+			int counter = 0;
+
+			for (decltype(grid.size()) i=0; i<grid.size(); i++) 
+			{
+				if (!grid.is_boundary_site(i))
+				{
+					counter++;
+					mag += statespace[i][0]; 
+				}
+			}
+
+			return sqrt(mag*mag)/double(counter);
+		}
+};
+
+
+// ------------------------------ HAMILTONIAN --------------------------
+
+
+// Finite Differences Interaction
+template <class StateSpace, class StateVector> 
+class FinDiffInt
+{
+
+	public:
+		const double k;
+		FinDiffInt(double k) : k(k) {}
+		~FinDiffInt() {}
+
+        double diff (const int rsite,
+						const StateVector& svold,
+						const StateVector& svnew,
+						std::vector<int>& nbrs,
+						StateSpace& s) 
+
+		{
+			double energy_old = 0;
+			double energy_new  = 0;
+			for (std::size_t i=0; i<nbrs.size(); ++i)
+			{                                     
+				auto idx = nbrs[i];                
+				auto nbr = s[idx];                  
+
+				energy_old += pow((svold[0]-nbr[0]),2);
+				energy_new += pow((svnew[0]-nbr[0]),2);
+			}
+			return energy_new-energy_old;
+		}
+
+};
+
+
+
+
+/** Action of a massive scalar field
+ */
+class MassiveScalarField
+{
+	public:
+		//  ----  Parameters  ----
+
+		double beta, k, mass, sqrtg;
+		static constexpr int SymD = 1;
+		const std::string name = "MassiveScalarField";
+
+
+		//  ---- Definitions  -----
+
+		typedef std::array<double, SymD> StateVector;
+		typedef MARQOV::Space<StateVector,HyperbolicRegularFromCSV> StateSpace;
+
+
+
+		//  ----  Hamiltonian terms  ----
+
+		std::array<OnSite<StateVector, double>*,1> onsite 
+			= {new Onsite_Quadratic<StateVector>(0.5*sqrtg*mass)};
+		
+		std::array<FinDiffInt<StateSpace,StateVector>*,1> multisite 
+			= {new FinDiffInt<StateSpace, StateVector>(0.5*k)};
+
+		MassiveScalarField(double k, double sqrtg, double mass) : k(k), 
+																  mass(mass), 
+																  sqrtg(sqrtg),
+																  name("MassiveScalarField") {}
+
+
+
+		//  ----  Observables  ----
+
+		Magnetization obs_m;
+		BulkMagnetization obs_mb;
+        decltype(std::make_tuple(obs_m, obs_mb)) observables = {std::make_tuple(obs_m, obs_mb)};
+
+
+
+		//  ---- Parameter Names ----
+
+		/** Allows to give the Hamiltonian parameter names
+		*
+		* @param i index of the parameter
+		*
+		* @return the parameter name (string)
+		*/
+		std::string paramname(int i)
+		{
+			std::string name;
+			switch (i)
+			{
+				case (0): name = "k"; break;
+				case (1): name = "sqrtg"; break;
+				case (2): name = "mass"; break;
+				default: break;
+			}
+			return name;
+		}
+
+
+		//  ----  Initializer  ----
+
+		template <class StateSpace, class Lattice, class RNG>
+		void initstatespace(StateSpace& statespace, Lattice& grid, RNG& rng) const
+		{
+			double boundary_value = 1000;
+			double central_value = 10;
+
+			for(decltype(grid.size()) i = 0; i < grid.size(); ++i)
+			{
+				// random initialization with mean zero is not a good idea!
+				// it can introduce frustration effects!
+
+				/*
+				// boundary values
+				if (grid.is_boundary_site(i)) statespace[i][0] = boundary_value;
+				// remaining sites
+				else statespace[i][0] = rng.real(0, 2*boundary_value);
+				*/
+
+
+				// central value
+				auto firstlayer = grid.flexnbrs(0,0);
+				if (i==0) statespace[i][0] = central_value;
+				else if (std::find(firstlayer.begin(), firstlayer.end(), i) != firstlayer.end()) 
+				{
+					statespace[i][0] = central_value;
+				}
+				// remaining sites
+				else statespace[i][0] = rng.real(0, 2*central_value);
+
+			}
+		}
+};
+
+
+// ------------------------------ INITIALIZER ---------------------------
+
+template <>
+class Initializer<MassiveScalarField>
+{
+	public:
+		Initializer() {}
+
+		template <class RNGCache>
+		static typename MassiveScalarField::StateVector newsv(const typename MassiveScalarField::StateVector& svold, RNGCache& rng)
+		{
+			double amp = 0.5;
+			double r = rng.real(-1.0, 1.0);
+
+			double oldval = svold[0];
+			double newval = oldval + amp*r;
+			auto nsv = svold;
+			nsv[0] = newval;
+			return nsv;
+		}
+};
+
+
+
+namespace MARQOV
+{
+
+	template <class Lattice>
+	struct Metropolis<MassiveScalarField, Lattice>
+	{
+	    template <class StateSpace, class RNG>
+	    static int move(const MassiveScalarField& ham, const Lattice& grid, StateSpace& statespace, RNG& rng, double beta, int rsite)
+	    {
+
+			// boundary cells are not being updated
+//			if (grid.is_boundary_site(rsite)) return 0;
+			// central cell is not being updated
+			if (rsite == 0) return 0;
+				auto firstlayer = grid.flexnbrs(0,0);
+			if (std::find(firstlayer.begin(), firstlayer.end(), rsite) != firstlayer.end()) 
+				{
+				return 0;
+				}
+
+			// typedefs
+			typedef MassiveScalarField Hamiltonian;
+			typedef typename Hamiltonian::StateVector StateVector;
+
+			// old state vector at index rsite
+			StateVector& svold = statespace[rsite];
+			// propose new configuration
+			StateVector svnew(Initializer<Hamiltonian>::newsv(svold, rng) );
+
+			// no interaction term
+
+			// only one on-site term
+			const int b = 0;
+			auto diffb = ham.onsite[b]->get(svnew) - ham.onsite[b]->get(svold);
+			double onsiteenergydiff = dot(ham.onsite[b]->h, diffb);
+
+			// one flex term
+			const int c = 0; 
+			auto nbrs = getflexnbrs<Lattice>(grid, c, rsite);
+			auto diffc = ham.multisite[c]->diff(rsite, svold, svnew, nbrs, statespace);
+			double flexenergydiff = dot(ham.multisite[c]->k, diffc);
+
+
+
+	    	// sum up energy differences
+	    	double dE 	= onsiteenergydiff + flexenergydiff;
+	
+	    	int retval = 0;
+	    	if ( dE <= 0 )
+	    	{
+	    	    svold = svnew;
+	    	    retval = 1;
+	    	}
+	    	else if (rng.real() < exp(-beta*dE))
+	    	{
+	    	    svold = svnew;
+	    	    retval = 1;
+	    	}
+	    	return retval;
+	    }
+	};
+
+
+
+	// write specialization for Wolff in order to get rid of "interactions" in Hamiltonian
+
+	template <class Lattice>
+	struct Wolff<MassiveScalarField, Lattice>
+	{
+	    template <class StateSpace, class RNG>
+	    static int move(const MassiveScalarField& ham, const Lattice& grid, StateSpace& statespace, RNG& rng, double beta, int rsite)
+	    {
+
+			return 0;
+	    }
+	};
+}
+
+
+#endif
diff --git a/src/hamiltonian/MassiveScalarField2.h b/src/hamiltonian/MassiveScalarField2.h
new file mode 100644
index 0000000000000000000000000000000000000000..ad821dc8ad42f6cd9603f9edb2ce00590e06a844
--- /dev/null
+++ b/src/hamiltonian/MassiveScalarField2.h
@@ -0,0 +1,196 @@
+#ifndef MASSIVESCALARFIELD2_H
+#define MASSIVESCALARFIELD2_H
+#include <array>
+#include <cmath>
+#include <vector>
+#include "util/randomdir.h"
+#include "util/hamparts.h"
+#include "util/termcollection.h" 
+
+
+// ------------------------------ OBSERVABLES ---------------------------
+
+#include "util/observables.h"
+
+
+// ------------------------------ HAMILTONIAN --------------------------
+
+
+/** Action of a massive scalar field
+ */
+class MassiveScalarField2
+{
+	public:
+		//  ----  Parameters  ----
+
+		int q = 7; // hardcoded!!
+		double beta, k, mass, sqrtg;
+		static constexpr int SymD = 1;
+		const std::string name = "MassiveScalarField2";
+
+
+		//  ---- Definitions  -----
+
+		typedef std::array<double, SymD> StateVector;
+		typedef MARQOV::Space<StateVector,HyperbolicRegularFromCSV> StateSpace;
+
+
+
+		//  ----  Hamiltonian terms  ----
+
+		std::array<Standard_Interaction<StateVector>*, 1> interactions 
+			= {new Standard_Interaction<StateVector>(-k)};
+
+		std::array<Onsite_Quadratic<StateVector>*, 1> onsite 
+			= {new Onsite_Quadratic<StateVector>(0.5*(q*k+sqrtg*mass))};
+
+		MassiveScalarField2(double k, double sqrtg, double mass) : k(k), 
+																  mass(mass), 
+																  sqrtg(sqrtg),
+																  name("MassiveScalarField")
+		{
+		}
+
+
+
+		//  ----  Observables  ----
+
+		Magnetization obs_m;
+        decltype(std::make_tuple(obs_m)) observables = {std::make_tuple(obs_m)};
+
+
+
+		//  ---- Parameter Names ----
+
+		/** Allows to give the Hamiltonian parameter names
+		*
+		* @param i index of the parameter
+		*
+		* @return the parameter name (string)
+		*/
+		std::string paramname(int i)
+		{
+			std::string name;
+			switch (i)
+			{
+				case (0): name = "k"; break;
+				case (1): name = "sqrtg"; break;
+				case (2): name = "mass"; break;
+				default: break;
+			}
+			return name;
+		}
+
+
+		//  ----  Initializer  ----
+
+		template <class StateSpace, class Lattice, class RNG>
+		void initstatespace(StateSpace& statespace, Lattice& grid, RNG& rng) const
+		{
+			for(decltype(grid.size()) i = 0; i < grid.size(); ++i)
+			{
+//				auto nnbrs = grid.nbrs(0,i).size();
+//				if (nnbrs < 7) statespace[i][0] = 1;
+//				else statespace[i] = rnddir<RNG, typename StateVector::value_type, SymD>(rng);
+				statespace[i] = rnddir<RNG, typename StateVector::value_type, SymD>(rng);
+			}
+		}
+};
+
+
+// ------------------------------ INITIALIZER ---------------------------
+
+template <>
+class Initializer<MassiveScalarField2>
+{
+	public:
+		Initializer() {}
+
+		template <class RNGCache>
+		// auto?
+		static typename MassiveScalarField2::StateVector newsv(const typename MassiveScalarField2::StateVector& svold, RNGCache& rng)
+		{
+			double amp = 0.5; // Amplitude (TODO: make me a class parameter)
+			double r = rng.real(-1.0, 1.0);
+
+			double oldval = svold[0];
+			double newval = oldval + amp*r;
+			auto nsv = svold;
+			nsv[0] = newval;
+			return nsv;
+		}
+};
+
+
+
+
+namespace MARQOV
+{
+
+	template <class Lattice>
+	struct Metropolis<MassiveScalarField2, Lattice>
+	{
+	    template <class StateSpace, class RNG>
+	    static int move(const MassiveScalarField2& ham, const Lattice& grid, StateSpace& statespace, RNG& rng, double beta, int rsite)
+	    {
+
+			typedef MassiveScalarField2 Hamiltonian;
+			typedef typename Hamiltonian::StateVector StateVector;
+
+
+			// old state vector at index rsite
+			StateVector& svold = statespace[rsite];
+			// propose new configuration
+			StateVector svnew(Initializer<Hamiltonian>::newsv(svold, rng) );
+
+			const int a = 0; // only one interaction term
+			const auto nbrs = getnbrs<Lattice>(grid, a, rsite);
+
+			StateVector neighbourhood = {0};
+
+			// sum over neighbours
+			for (std::size_t i=0; i<nbrs.size(); ++i)
+			{	
+				// index of the neighbour
+				auto idx = nbrs[i];
+
+				// configuration of the neighbour after applying the interaction
+				auto nbr = ham.interactions[a]->get(statespace[idx]);
+
+				// sum up contributions from neighbourbood
+				neighbourhood = neighbourhood + nbr;
+			}
+
+			double interactionenergydiff = ham.interactions[a]->J * (dot(svnew-svold, neighbourhood));
+
+
+
+			const int b = 0; // only one on-site term
+			auto diff = ham.onsite[b]->get(svnew) - ham.onsite[b]->get(svold);
+			double onsiteenergydiff = dot(ham.onsite[b]->h, diff);
+
+
+
+	    	// sum up energy differences
+	    	double dE 	= interactionenergydiff + onsiteenergydiff;
+	
+	    	int retval = 0;
+	    	if ( dE <= 0 )
+	    	{
+	    	    svold = svnew;
+	    	    retval = 1;
+	    	}
+	    	else if (rng.real() < exp(-beta*dE))
+	    	{
+	    	    svold = svnew;
+	    	    retval = 1;
+	    	}
+	    	return retval;
+	    }
+	};
+}
+
+
+
+
+#endif
diff --git a/src/kleingordon.cpp b/src/kleingordon.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ae5bb71c362e8dd0f748e5d412638b05cf01260
--- /dev/null
+++ b/src/kleingordon.cpp
@@ -0,0 +1,116 @@
+/* MARQOV - A modern framework for classical spin models on general topologies
+ * Copyright (C) 2020-2021, The MARQOV Project
+ *
+ * This program is free software: you can redistribute it and/or modify it
+ * under the terms of the GNU General Public License as published by the Free
+ * Software Foundation, either version 3 of the License, or (at your option)
+ * any later version.
+ * 
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+ * See the GNU General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License along with
+ * this program. If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#include <array>
+#include <vector>
+#include <iostream>
+#include <string>
+#include <cstdlib>
+#include <fstream>
+#include <algorithm>
+#include <tuple>
+#include <iomanip>
+
+using std::cout;
+using std::endl;
+using std::flush;
+using std::ofstream;
+
+// MARQOV
+#include "libmarqov/libmarqov.h"
+
+// Lattices
+#include "lattice/graph_from_csv.h"
+
+// Hamiltonians
+#include "hamiltonian/MassiveScalarField.h"
+#include "hamiltonian/MassiveScalarField2.h"
+
+using namespace MARQOV;
+
+
+int main()
+{
+
+    std::cout<<"MARQOV Copyright (C) 2020-2021, The MARQOV Project contributors"<<std::endl;
+    std::cout<<"This program comes with ABSOLUTELY NO WARRANTY."<<std::endl;
+    std::cout<<"This is free software, and you are welcome to redistribute it under certain conditions."<<std::endl;
+    
+
+    
+	// remove old output and prepare new one
+	std::string outbasedir = "../out";
+//	std::string command;
+//	command = "rm -r " + outbasedir;
+//	system(command.c_str());
+//	makeDir(outbasedir);
+
+	
+
+	// Parameters
+	std::vector<double> beta = {1.0};
+	std::vector<double> K = {0.50315223}; // geometric factor
+	std::vector<double> sqrtg = {1.0471975511}; // geometric factor
+	std::vector<double> mass = {-10, -3, -1, -0.3, -0.1, -0.03, -0.01, 0, 1, 3, 10};	// mass squared
+	// note that the geometric factors are only valid for the bulk region!
+
+	auto hp = cart_prod(beta, K, sqrtg, mass);
+	
+	const int nthreads = 4;	
+	const int nreplicas = 1;
+
+	
+	// Typedefs
+	typedef MassiveScalarField Hamiltonian;
+	typedef HyperbolicRegularFromCSV Lattice;
+
+
+	// Lattice
+	int p = 3;
+	int q = 7; 
+	int nlayers = 8; 
+
+	std::string filebase = "/home/schrauth/hyperbolic-";
+	HyperbolicRegularFromCSV lat(filebase, p, q, nlayers);
+
+
+	// Scheduler
+    typedef typename std::tuple<Lattice&, MARQOV::Config, std::tuple<double,double,double,double> > ParameterType;
+	typedef typename GetSchedulerType<Hamiltonian, Lattice, ParameterType>::MarqovScheduler SchedulerType;
+ 	SchedulerType sched(1,nthreads);
+
+
+	// Output
+	std::string outpath = outbasedir+"/fixedcentral-"+std::to_string(nlayers)+"/";
+	makeDir(outpath);
+	
+
+	// Monte Carlo
+	MARQOV::Config mp(outpath);
+	mp.setnmetro(10);
+	mp.setncluster(0);
+	mp.setwarmupsteps(0);
+	mp.setgameloopsteps(1000);
+
+
+	// Schedule and Run
+	auto params  = finalize_parameter(lat, mp, hp);
+    auto rparams = replicator(params, nreplicas);
+ 	for (auto p: rparams) sched.createSimfromParameter(p, hyperfilter);
+	sched.start();
+
+}
diff --git a/src/lattice/graph_from_csv.h b/src/lattice/graph_from_csv.h
index facab010ef336067e828e0a7ba917f654dcfee8d..1f846e2b0fa8badeadaa5c89989d4b0f6cfb0c00 100644
--- a/src/lattice/graph_from_csv.h
+++ b/src/lattice/graph_from_csv.h
@@ -47,6 +47,41 @@ class GraphFromCSV
 };
 
 
+
+class HyperbolicRegularFromCSV : public GraphFromCSV
+{
+	public:
+		int p, q, nlayers;
+		HyperbolicRegularFromCSV(std::string filebasename, int p, int q, int nlayers) : p(p), q(q), nlayers(nlayers), 
+		GraphFromCSV(filebasename+std::to_string(p)+"-"+std::to_string(q)+"-"+std::to_string(nlayers)+".csv") 
+		{
+			// subtract neighbour indices by one
+			for (int i=0; i<neighbours.size(); i++)
+			{
+				for (int j=0; j<neighbours[i].size(); j++)
+				{
+					neighbours[i][j]--;
+				}
+			}
+		}
+
+		// check whether site is boundary site
+		bool is_boundary_site(int idx) const
+		{
+			auto nnbrs = this->neighbours[idx].size();
+			if (nnbrs < q) return true;
+			else return false;
+		}
+
+		// implement flexnbrs
+    	std::vector<int> flexnbrs(const int a, const int i) const
+    	{
+    	    return this->neighbours[i];
+    	}
+
+};
+
+
 /* Read graph with coordinates from csv file.
 * @note for format specifications see general documentation
 */
@@ -58,7 +93,7 @@ class GraphFromCSVwithCoords
     std::vector<std::vector<int>> neighbours;
     std::vector<std::vector<double>> grid;
 
-    GraphFromCSV(std::string filename, ncoords) 
+    GraphFromCSVwithCoords(std::string filename, int ncoords) 
     {
         npoints = import_geometry(filename, grid, neighbours, ncoords);
     }
diff --git a/src/lattice/util/io.h b/src/lattice/util/io.h
index b4d28ae35677fbe57d307a248b8b989ff6aede85..dd49c0c8eb4c7417480fbd2394d719557a6df2a6 100644
--- a/src/lattice/util/io.h
+++ b/src/lattice/util/io.h
@@ -84,55 +84,6 @@ void save_geometry_deluxe(Grid& grid, const std::string path)
 }
 
 
-/**
-* Import CSV lattice file containing coordinates and connections between vertices.
-*
-* @param path the filename
-* @param grid the coordinates will be stored here
-* @param nbrs the bonds will be stored here
-* @param ncoords number of coordinates per point
-*
-* @return number of points
-*/
-int import_geometry(const std::string path, std::vector<std::vector<double>>& grid, std::vector<std::vector<int>>& nbrs, int ncoords)
-{
-	int idxoffs = 1;	
-	// set to 1 if vertex indices in the CSV are counted from 1. Set 0 zero if they are counted from 0
-
-	std::ifstream in(NULL);
-	in.open(path.c_str());
-	std::string row;
-
-	// loop over lines
-    while (!in.eof()) 
-	{
-        std::getline(in, row);
-        if (in.bad() || in.fail()) 
-		{
-            break;
-        }
-		std::istringstream ss(row);
-		std::string substr;
-
-		int counter = 0;
-		std::vector<double> g;
-		std::vector<int> n;
-
-		// loop over "words" in a line
-		while(std::getline(ss, substr, '\t'))
-		{
-			// if coordinate, transform to double
-			if (counter < ncoords) g.push_back(std::stod(substr));	
-			// if vertex, transform to int
-			else n.push_back(std::stoi(substr)-idxoffs);
-			counter++;
-    	}
-		grid.push_back(g);
-		nbrs.push_back(n);
-	}
-	return nbrs.size();
-}
-
 /**
 * @brief Import coordinates and neighbour relations from CSV file
 *
diff --git a/src/libmarqov/core.h b/src/libmarqov/core.h
index 82c0e0e191907a378848f15088bda5d426cb7371..5d6f7a2d740c3720a06234471990da7f8a0ed662 100644
--- a/src/libmarqov/core.h
+++ b/src/libmarqov/core.h
@@ -615,7 +615,7 @@ class Core : public RefType<Grid>
             auto retval = new typename Hamiltonian::StateVector[size];
             if (step > 0)
             {
-                std::cout<<"[MARQOV::Core] Previous data found! Continuing simulation at step "<<step<<std::endl;
+//                std::cout<<"[MARQOV::Core] Previous data found! Continuing simulation at step "<<step<<std::endl;
                 //read in the state space
                 auto prevstepstate = dump.openGroup("/step" + std::to_string(step-1) + "/state");
                 {
diff --git a/src/libmarqov/emcs.h b/src/libmarqov/emcs.h
index 208fc312ccaa4971820413b8556dc2bc542b2875..b92eededb62e37db945cc18eb1b3206536077844 100644
--- a/src/libmarqov/emcs.h
+++ b/src/libmarqov/emcs.h
@@ -64,13 +64,14 @@ namespace MARQOV
 
         // Metropolis sweeps
         mrqvt.switch_clock("metrop");
+		double avgacceptance;
         for (int j=0; j < nmetro; j++)
         {
             // loop sites
             for(decltype(grid.size()) i = 0; i < grid.size(); ++i)
             {
                 const int rsite = rngcache.integer(grid.size());
-                Metropolis<Hamiltonian, Lattice>::move(ham,
+                int avgacceptance = Metropolis<Hamiltonian, Lattice>::move(ham,
                                                 grid, 
                                                 statespace, 
                                                 rngcache, 
@@ -79,7 +80,8 @@ namespace MARQOV
             }
         }
 
-        return avgclustersize/ncluster;
+        return avgacceptance/grid.size()/double(nmetro);
+//        return avgclustersize/ncluster;
     }
     
 template <class Grid, class Hamiltonian, class MutexType, class RNGType, template<class> class RefType>