diff --git a/include/IsingModel.hpp b/include/IsingModel.hpp
index 9f71bb7..8ec4ae9 100644
--- a/include/IsingModel.hpp
+++ b/include/IsingModel.hpp
@@ -21,7 +21,6 @@
 #include <random>
 #include <unordered_map>
 
-// Faster modulo
 #define INDEX(I, N) (I + N) % N
 
 // Indeces for the neighbor matrix.
@@ -50,7 +49,7 @@ private:
      *  for the left and upper neighbor, and we can use the same column for the
      *  right and lower neighbor.
      * */
-    arma::Mat<uint> neighbors;
+    arma::Mat<int> neighbors;
 
     /** @brief A hash map containing all possible energy changes.
      * */
@@ -62,15 +61,15 @@ private:
 
     /** @brief Size of the lattice.
      * */
-    uint L;
+    int L;
 
     /** @brief The current energy state. unit: \f$ J \f$.
      * */
-    double E;
+    int E;
 
     /** @brief The current magnetic strength. unit: Unitless.
      * */
-    double M;
+    int M;
 
     /** @brief Initialize the lattice with a random distribution of 1s and
      *  -1s.
@@ -103,7 +102,7 @@ public:
      *  @param L The size of the lattice.
      *  @param T The temperature for the system.
      * */
-    IsingModel(uint L, double T);
+    IsingModel(int L, double T);
 
     /** @brief Constructor for the Ising model.
      *
@@ -111,23 +110,23 @@ public:
      *  @param T The temperature for the system.
      *  @param val The value to set for all spins.
      * */
-    IsingModel(uint L, double T, int val);
+    IsingModel(int L, double T, int val);
 
     /** @brief The Metropolis algorithm.
      * */
-    data_t Metropolis(std::mt19937 &engine);
+    data_t Metropolis();
 
     /** @brief Get the current energy.
      *
      *  @return double
      * */
-    double get_E();
+    int get_E();
 
     /** @brief Get the current magnetization.
      *
      *  @return double
      * */
-    double get_M();
+    int get_M();
 };
 
 #endif
diff --git a/include/data_type.hpp b/include/data_type.hpp
index 7382db2..5d3bfd5 100644
--- a/include/data_type.hpp
+++ b/include/data_type.hpp
@@ -15,64 +15,105 @@
 #include <sys/types.h>
 #include <type_traits>
 
-/** @brief Data structure that contains the data needed for the project*/
-struct data_t {
-    double E = 0.; ///< The expected energy
-    double M = 0.; ///< The expected magnetization
-    double E2 = 0.; ///< The expected variance of the energy
-    double M2 = 0.; ///< The expected variance of magnetization
-    double M_abs = 0.; ///< The expected absolute magnetization
+class data_t {
+public:
+    double E, M, E2, M2, M_abs;
+
+    data_t()
+    {
+        this->E = 0.;
+        this->E2 = 0.;
+        this->M = 0.;
+        this->M2 = 0.;
+        this->M_abs = 0.;
+    }
+
+    data_t(double E, double E2, double M, double M2, double M_abs)
+    {
+        this->E = E;
+        this->E2 = E2;
+        this->M = M;
+        this->M2 = M2;
+        this->M_abs = M_abs;
+    }
+
+    template <class T> data_t operator/(T num)
+    {
+        data_t res;
+        res.E = this->E / (double)num;
+        res.E2 = this->E2 / (double)num;
+        res.M = this->M / (double)num;
+        res.M2 = this->M2 / (double)num;
+        res.M_abs = this->M_abs / (double)num;
+
+        return res;
+    }
+
+    template <class T> data_t& operator/=(T num)
+    {
+        this->E /= (double)num;
+        this->E2 /= (double)num;
+        this->M /= (double)num;
+        this->M2 /= (double)num;
+        this->M_abs /= (double)num;
+
+        return *this;
+    }
+
+    template <class T> data_t operator*(T num)
+    {
+        data_t res;
+        res.E = this->E * (double)num;
+        res.E2 = this->E2 * (double)num;
+        res.M = this->M * (double)num;
+        res.M2 = this->M2 * (double)num;
+        res.M_abs = this->M_abs * (double)num;
+
+        return res;
+    }
+
+    template <class T> data_t& operator*=(T num)
+    {
+        this->E *= (double)num;
+        this->E2 *= (double)num;
+        this->M *= (double)num;
+        this->M2 *= (double)num;
+        this->M_abs *= (double)num;
+
+        return *this;
+    }
+
+    data_t operator+(const data_t &b)
+    {
+        data_t res;
+        res.E = this->E + b.E;
+        res.E2 = this->E2 + b.E2;
+        res.M = this->M + b.M;
+        res.M2 = this->M2 + b.M2;
+        res.M_abs = this->M_abs + b.M_abs;
+
+        return res;
+    }
+
+    data_t& operator+=(const data_t &b)
+    {
+        this->E += b.E;
+        this->E2 += b.E2;
+        this->M += b.M;
+        this->M2 += b.M2;
+        this->M_abs += b.M_abs;
+
+        return *this;
+    }
+
+    template <class T> void operator=(T num)
+    {
+        this->E = (double)num;
+        this->E2 = (double)num;
+        this->M = (double)num;
+        this->M2 = (double)num;
+        this->M_abs = (double)num;
+    }
 };
 
-/** @brief Define dividing data_t by a number
- *
- *  @param data The data to divide
- *  @param num The number to divide data by
- *
- *  @return data_t
- * */
-template <class T>
-data_t operator/(const data_t &data, T num);
-
-// Explicit instantiation
-extern template data_t operator/(const data_t &, uint);
-extern template data_t operator/(const data_t &, ulong);
-extern template data_t operator/(const data_t &,int);
-extern template data_t operator/(const data_t &,double);
-
-
-/** @brief Define /= on data_t by a number.
- *
- *  @param data The data to divide
- *  @param num The number to divide data by
- *
- *  @return data_t
- * */
-template <class T>
-data_t& operator/=(data_t &data, T num);
-
-// Explicit instantiation
-extern template data_t& operator/=(data_t &, uint);
-extern template data_t& operator/=(data_t &, ulong);
-extern template data_t& operator/=(data_t &,int);
-extern template data_t& operator/=(data_t &,double);
-
-/** @brief Define + on data_t by a data_t.
- *
- *  @param a The left side
- *  @param b The right side
- *
- *  @return data_t
- * */
-data_t operator+(const data_t &a, const data_t &b);
-
-/** @brief Define += on data_t by a data_t.
- *
- *  @param a The left side
- *  @param b The right side
- *
- *  @return data_t
- * */
-data_t& operator+=(data_t &a, const data_t &b);
-
 #endif
diff --git a/include/monte_carlo.hpp b/include/monte_carlo.hpp
index 7be8f4f..36e7b05 100644
--- a/include/monte_carlo.hpp
+++ b/include/monte_carlo.hpp
@@ -18,31 +18,21 @@
 
 #include <functional>
 #include <string>
+#include <omp.h>
 
-#define BURN_IN_TIME 1000
+//#define BURN_IN_TIME 12500
+#define BURN_IN_TIME 5000
 
-#define EPS_2 (-2 * std::sinh(8.)) / (std::cosh(8.) + 3)
-
-#define MAG_2 (std::exp(8.) + 1) / (2 * cosh(8.) + 3)
-
-#define CV_2                                                                   \
-    16                                                                         \
-        * (3 * std::cosh(8.) + std::cosh(8.) * std::cosh(8.)                   \
-           - std::sinh(8.) * std::sinh(8.))                                    \
-        / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))
-
-#define X_2                                                                    \
-    (3 * std::exp(8.) + std::exp(-8.) + 3)                                     \
-        / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))
+#pragma omp declare reduction(+: data_t: omp_out += omp_in)
 
 /** @brief Test numerical data with analytical data.
  *
  *  @param tol The tolerance between the analytical and numerical solution.
  *  @param max_cycles The max number of Monte Carlo cycles.
  *
- *  return uint
+ *  return int
  * */
-uint test_2x2_lattice(double tol, uint max_cycles);
+int test_2x2_lattice(double tol, int max_cycles);
 
 /** @brief Write the expected values for each Monte Carlo cycles to file.
  *
@@ -51,7 +41,18 @@ uint test_2x2_lattice(double tol, uint max_cycles);
  *  @param cycles The amount of Monte Carlo cycles to do
  *  @param filename The file to write to
  * */
-void monte_carlo_progression(double T, uint L, uint cycles,
+void monte_carlo_progression(double T, int L, int cycles,
+                             const std::string filename);
+
+/** @brief Write the expected values for each Monte Carlo cycles to file.
+ *
+ *  @param T Temperature
+ *  @param L The size of the lattice
+ *  @param cycles The amount of Monte Carlo cycles to do
+ *  @param value The value to set the elements in the lattice
+ *  @param filename The file to write to
+ * */
+void monte_carlo_progression(double T, int L, int cycles, int value,
                              const std::string filename);
 
 /** @brief Estimate the probability distribution for the energy.
@@ -61,26 +62,29 @@ void monte_carlo_progression(double T, uint L, uint cycles,
  *  @param cycles The amount of Monte Carlo cycles to do
  *  @param filename The file to write to
  * */
-void pd_estimate(double T, uint L, uint cycles, const std::string filename);
+void pd_estimate(double T, int L, int cycles, const std::string filename);
 
-/** @brief Execute the Metropolis algorithm for a certain amount of Monte 
+/** @brief Execute the Metropolis algorithm for a certain amount of Monte
  *  Carlo cycles.
  *
- *  @param data The data to store the results
- *  @param L The size of the lattice
- *  @param T The Temperature for the Ising model
- *  @param cycles The amount of Monte Carlo cycles to do*/
-void monte_carlo_serial(data_t &data, uint L, double T, uint cycles);
-
-/** @brief Execute the Metropolis algorithm for a certain amount of Monte 
- *  Carlo cycles in parallel.
- *
- *  @param data The data to store the results
  *  @param L The size of the lattice
  *  @param T The Temperature for the Ising model
  *  @param cycles The amount of Monte Carlo cycles to do
+ *
+ *  @return data_t
  * */
-void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles);
+data_t monte_carlo_serial(int L, double T, int cycles);
+
+/** @brief Execute the Metropolis algorithm for a certain amount of Monte
+ *  Carlo cycles in parallel.
+ *
+ *  @param L The size of the lattice
+ *  @param T The Temperature for the Ising model
+ *  @param cycles The amount of Monte Carlo cycles to do
+ *
+ *  @return data_t
+ * */
+data_t monte_carlo_parallel(int L, double T, int cycles);
 
 /** @brief Perform the MCMC algorithm using a range of temperatures.
  *
@@ -92,8 +96,8 @@ void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles);
  *  @param outfile The file to write the data to
  * */
 void phase_transition(
-    uint L, double start_T, double end_T, uint points_T,
-    std::function<void(data_t &, uint, double, uint)> monte_carlo,
+    int L, double start_T, double end_T, int points_T,
+    std::function<data_t(int, double, int)> monte_carlo,
     std::string outfile);
 
 #endif
diff --git a/include/testlib.hpp b/include/testlib.hpp
index 452088c..853c649 100644
--- a/include/testlib.hpp
+++ b/include/testlib.hpp
@@ -8,7 +8,8 @@
  *  @brief A small test library.
  *
  *  @details This a small testing library that is tailored for the needs of the
- * project.
+ *  project. Anything that is in the details namespace should not be used 
+ *  directly, or else it might cause undefined behavior if not used correctly.
  *
  *  @bug No known bugs
  * */
@@ -90,7 +91,7 @@ template <class T,
           class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
 static bool close_to(T a, T b, double tol = 1e-8)
 {
-    return std::abs(a - b) < tol;
+    return std::fabs(a - b) < tol;
 }
 
 /** @brief Test if two armadillo matrices/vectors are equal.
diff --git a/include/utils.hpp b/include/utils.hpp
index 8cd0554..d9a99fd 100644
--- a/include/utils.hpp
+++ b/include/utils.hpp
@@ -8,7 +8,9 @@
  *  @brief Function prototypes and macros that are useful.
  *
  *  These utility function are mainly for convenience and aren't directly
- *  related to the project.
+ *  related to the project. Anything that is in the details namespace should
+ *  not be used directly, or else it might cause undefined behavior if not used
+ *  correctly.
  *
  *  @bug No known bugs
  * */
@@ -47,7 +49,7 @@ namespace details {
  *
  *  @details This function should only be used for the __METHOD_NAME__ macro,
  *  since it takes the output from __PRETTY_FUNCTION__ and strips the return
- *  type.
+ *  type. 
  *
  *  @param pretty_function The string from __PRETTY_FUNCTION__
  *
diff --git a/src/IsingModel.cpp b/src/IsingModel.cpp
index f61aa91..e0a600e 100644
--- a/src/IsingModel.cpp
+++ b/src/IsingModel.cpp
@@ -11,11 +11,14 @@
  * */
 #include "IsingModel.hpp"
 
+#include <cmath>
+#include <random>
+
 IsingModel::IsingModel()
 {
 }
 
-IsingModel::IsingModel(uint L, double T)
+IsingModel::IsingModel(int L, double T)
 {
     this->L = L;
     this->T = T;
@@ -26,7 +29,7 @@ IsingModel::IsingModel(uint L, double T)
     this->initialize_energy();
 }
 
-IsingModel::IsingModel(uint L, double T, int val)
+IsingModel::IsingModel(int L, double T, int val)
 {
     this->L = L;
     this->T = T;
@@ -41,13 +44,13 @@ IsingModel::IsingModel(uint L, double T, int val)
 void IsingModel::initialize_lattice()
 {
     this->lattice.set_size(this->L, this->L);
-    std::random_device rd;
-    std::mt19937 engine(rd());
+    std::random_device rd{};
+    std::mt19937 engine{rd()};
 
-    std::uniform_int_distribution<int> coin_flip(0, 1);
+    std::uniform_int_distribution<> coin_flip(0, 1);
 
     for (size_t i = 0; i < this->lattice.n_elem; i++)
-        this->lattice(i) = coin_flip(engine) == 1 ? 1 : -1;
+        this->lattice(i) = 2 * coin_flip(engine) - 1;
 }
 
 void IsingModel::initialize_neighbors()
@@ -64,7 +67,7 @@ void IsingModel::initialize_neighbors()
 void IsingModel::initialize_energy_diff()
 {
     for (int i = -8; i <= 8; i += 4) {
-        this->energy_diff.insert({i, std::exp(-((double)i / this->T))});
+        this->energy_diff.insert({i, std::exp(-(double)i / this->T)});
     }
 }
 
@@ -72,7 +75,7 @@ void IsingModel::initialize_magnetization()
 {
     this->M = 0.;
     for (size_t i = 0; i < this->lattice.n_elem; i++) {
-        this->M += (double)this->lattice(i);
+        this->M += this->lattice(i);
     }
 }
 
@@ -83,21 +86,23 @@ void IsingModel::initialize_energy()
     // Loop through the matrix
     for (size_t j = 0; j < this->L; j++) {
         for (size_t i = 0; i < this->L; i++) {
-            this->E -= (double)this->lattice(i, j)
+            this->E -= this->lattice(i, j)
                        * (this->lattice(i, this->neighbors(j, RIGHT))
                           + this->lattice(this->neighbors(i, DOWN), j));
         }
     }
 }
 
-data_t IsingModel::Metropolis(std::mt19937 &engine)
+data_t IsingModel::Metropolis()
 {
-    uint ri, rj;
+    std::random_device rd{};
+    std::mt19937_64 engine{rd()};
+
+    int ri, rj;
     int dE;
-    data_t res;
 
     // Create random distribution for indeces
-    std::uniform_int_distribution<uint> random_index(0, this->L - 1);
+    std::uniform_int_distribution<> random_index(0, this->L - 1);
     // Create random distribution for acceptance
     std::uniform_real_distribution<> random_number(0., 1.);
 
@@ -118,25 +123,21 @@ data_t IsingModel::Metropolis(std::mt19937 &engine)
         if (random_number(engine) <= this->energy_diff[dE]) {
             // Update if the configuration is accepted
             this->lattice(ri, rj) *= -1;
-            this->M += 2. * (double)this->lattice(ri, rj);
-            this->E += (double)dE;
+            this->M += 2 * this->lattice(ri, rj);
+            this->E += dE;
         }
     }
-    res.E = this->E;
-    res.E2 = this->E * this->E;
-    res.M = this->M;
-    res.M2 = this->M * this->M;
-    res.M_abs = std::abs(this->M);
 
-    return res;
+    return data_t((double)this->E, (double)(this->E * this->E), (double)this->M,
+                  (double)(this->M * this->M), std::fabs((double)this->M));
 }
 
-double IsingModel::get_E()
+int IsingModel::get_E()
 {
     return this->E;
 }
 
-double IsingModel::get_M()
+int IsingModel::get_M()
 {
     return this->M;
 }
diff --git a/src/Makefile b/src/Makefile
index dbf94c4..720b48c 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -16,7 +16,7 @@ OPENMP=-fopenmp
 # Add a debug flag when compiling (For the DEBUG macro in utils.hpp)
 DEBUG ?= 0
 ifeq ($(DEBUG), 1)
-	DBGFLAG=-DDBG
+	DBGFLAG=-DDBG -g
 else
 	DBGFLAG=
 endif
diff --git a/src/data_type.cpp b/src/data_type.cpp
index e602fe8..533ec79 100644
--- a/src/data_type.cpp
+++ b/src/data_type.cpp
@@ -10,63 +10,3 @@
  *  @bug No known bugs
  * */
 #include "data_type.hpp"
-
-template <class T>
-data_t operator/(const data_t &data, T num)
-{
-    data_t res = data;
-    res.E /= num;
-    res.E2 /= num;
-    res.M /= num;
-    res.M2 /= num;
-    res.M_abs /= num;
-
-    return res;
-}
-
-// Explicit instantiation
-template data_t operator/(const data_t &, uint);
-template data_t operator/(const data_t &, ulong);
-template data_t operator/(const data_t &,int);
-template data_t operator/(const data_t &,double);
-
-template <class T>
-data_t& operator/=(data_t &data, T num)
-{
-    data.E /= num;
-    data.E2 /= num;
-    data.M /= num;
-    data.M2 /= num;
-    data.M_abs /= num;
-
-    return data;
-}
-
-// Explicit instantiation
-template data_t& operator/=(data_t &, uint);
-template data_t& operator/=(data_t &, ulong);
-template data_t& operator/=(data_t &,int);
-template data_t& operator/=(data_t &,double);
-
-data_t operator+(const data_t &a, const data_t &b)
-{
-    data_t res = a;
-    res.E += b.E;
-    res.E2 += b.E2;
-    res.M += b.M;
-    res.M2 += b.M2;
-    res.M_abs += b.M_abs;
-
-    return res;
-}
-
-data_t& operator+=(data_t &a, const data_t &b)
-{
-    a.E += b.E;
-    a.E2 += b.E2;
-    a.M += b.M;
-    a.M2 += b.M2;
-    a.M_abs += b.M_abs;
-
-    return a;
-}
diff --git a/src/main.cpp b/src/main.cpp
index dc4a72d..ce838a6 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -9,35 +9,105 @@
  *
  *  @bug No known bugs
  * */
+#include "data_type.hpp"
 #include "monte_carlo.hpp"
 #include "utils.hpp"
 
+#include <csignal>
+#include <cstdlib>
 #include <iostream>
 #include <omp.h>
 
-/** @brief The main function.*/
-int main()
+void create_burn_in_time_data()
 {
-    // uint test_cycles = test_2x2_lattice(1e-7, 10000);
-    // monte_carlo(1.0, 2, 10000, "output/2_lattice_test.txt");
-
     // Test burn-in time
-    monte_carlo_progression(1.0, 20, 10000, "output/burn_in_time_1_0.txt");
-    monte_carlo_progression(2.4, 20, 10000, "output/burn_in_time_2_4.txt");
+    monte_carlo_progression(1.0, 20, 20000,
+                            "output/burn_in_time/unordered_1_0.txt");
+    monte_carlo_progression(1.0, 20, 20000, 1,
+                            "output/burn_in_time/ordered_1_0.txt");
+    monte_carlo_progression(2.4, 20, 20000,
+                            "output/burn_in_time/unordered_2_4.txt");
+    monte_carlo_progression(2.4, 20, 20000, 1,
+                            "output/burn_in_time/ordered_2_4.txt");
+}
 
+void create_pd_estimate_data()
+{
+    // Estimate pd
+    pd_estimate(1.0, 20, 1000000, "output/pd_estimate/estimate_1_0.txt");
+    pd_estimate(2.4, 20, 1000000, "output/pd_estimate/estimate_2_4.txt");
+}
+
+void test_parallel_speedup()
+{
     // Test the openmp speedup
+    data_t data;
     double t0, t1, t2;
+    int tries = 5;
     t0 = omp_get_wtime();
-    phase_transition(20, 2.1, 2.4, 1000, monte_carlo_serial,
-                     "output/phase_transition/size_20.txt");
+    for (size_t i = 0; i < tries; i++)
+        monte_carlo_serial(20, 1.0, 10000);
     t1 = omp_get_wtime();
-    phase_transition(20, 2.1, 2.4, 1000, monte_carlo_parallel,
-                     "output/phase_transition/size_20.txt");
+    for (size_t i = 0; i < tries; i++)
+        monte_carlo_parallel(20, 1.0, 10000);
     t2 = omp_get_wtime();
 
-    std::cout << "Time serial     : " << t1 - t0 << " seconds" << '\n';
-    std::cout << "Time parallel   : " << t2 - t1 << " seconds" << '\n';
+    std::cout << "Time serial     : " << (t1 - t0) / tries << " seconds"
+              << '\n';
+    std::cout << "Time parallel   : " << (t2 - t1) / tries << " seconds"
+              << '\n';
     std::cout << "Speedup parallel: " << (t1 - t0) / (t2 - t1) << '\n';
+}
+
+void create_phase_transition_data()
+{
+    double t0, t1;
+
+    t0 = omp_get_wtime();
+    // Phase transition
+    phase_transition(20, 2.1, 2.4, 40, monte_carlo_parallel,
+                     "output/phase_transition/size_20.txt");
+    phase_transition(40, 2.1, 2.4, 40, monte_carlo_parallel,
+                     "output/phase_transition/size_40.txt");
+    phase_transition(60, 2.1, 2.4, 40, monte_carlo_parallel,
+                     "output/phase_transition/size_60.txt");
+    phase_transition(80, 2.1, 2.4, 40, monte_carlo_parallel,
+                     "output/phase_transition/size_80.txt");
+    phase_transition(100, 2.1, 2.4, 40, monte_carlo_parallel,
+                     "output/phase_transition/size_100.txt");
+    t1 = omp_get_wtime();
+
+    std::cout << "Time: " << t1 - t0 << std::endl;
+}
+
+/** @brief The main function.*/
+int main(int argc, char **argv)
+{
+    if (argc < 2) {
+        std::cout << "Need at least 1 argument, got " << argc - 1
+                  << " arguments." << std::endl;
+        abort();
+    }
+
+    int arg = atoi(argv[1]);
+
+    switch (arg) {
+    case 1:
+        create_burn_in_time_data();
+        break;
+    case 2:
+        create_pd_estimate_data();
+        break;
+    case 3:
+        test_parallel_speedup();
+        break;
+    case 4:
+        create_phase_transition_data();
+        break;
+    default:
+        std::cout << "Not a valid option!" << std::endl;
+        abort();
+    }
 
     return 0;
 }
diff --git a/src/monte_carlo.cpp b/src/monte_carlo.cpp
index b3bb867..4f64092 100644
--- a/src/monte_carlo.cpp
+++ b/src/monte_carlo.cpp
@@ -11,50 +11,15 @@
  * */
 #include "monte_carlo.hpp"
 
-uint test_2x2_lattice(double tol, uint max_cycles)
-{
-    data_t data, tmp;
-    size_t L = 2;
-    size_t n_spins = L * L;
-    double T = 1.;
-    size_t cycles = 0;
+#include <cmath>
+#include <cstdint>
 
-    // Create random engine using the mersenne twister
-    std::random_device rd;
-    std::mt19937 engine(rd());
-
-    IsingModel test(L, T);
-
-    double E, M;
-
-    std::function<bool(double, double)> is_close = [tol](double a, double b) {
-        return std::abs(a - b) < tol;
-    };
-
-    // Loop through cycles
-    while (cycles++ < max_cycles) {
-        data += test.Metropolis(engine);
-        tmp = data / (cycles * n_spins);
-        // if (close(EPS_2, tmp.E)
-        //&& close(MAG_2, tmp.M)
-        //&& close(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T))
-        //&& close(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T)) {
-        // return cycles;
-        //}
-        if (is_close(EPS_2, tmp.E) && is_close(MAG_2, tmp.M)) {
-            return cycles;
-        }
-    }
-    std::cout << "hello" << std::endl;
-    return 0;
-}
-
-void monte_carlo_progression(double T, uint L, uint cycles,
+void monte_carlo_progression(double T, int L, int cycles,
                              const std::string filename)
 {
     // Set some variables
     data_t data, tmp;
-    uint n_spins = L * L;
+    int n_spins = L * L;
 
     // File stuff
     std::string directory = utils::dirname(filename);
@@ -62,7 +27,10 @@ void monte_carlo_progression(double T, uint L, uint cycles,
 
     // Create random engine using the mersenne twister
     std::random_device rd;
-    std::mt19937 engine(rd());
+    uint32_t rd_sample = rd();
+    std::mt19937 engine(rd_sample);
+
+    std::cout << "Seed: " << rd_sample << std::endl;
 
     IsingModel ising(L, T);
 
@@ -72,7 +40,7 @@ void monte_carlo_progression(double T, uint L, uint cycles,
 
     // Loop through cycles
     for (size_t i = 1; i <= cycles; i++) {
-        data += ising.Metropolis(engine);
+        data += ising.Metropolis();
         tmp = data / (i * n_spins);
         ofile << i << ',' << tmp.E << ',' << tmp.E2 << ',' << tmp.M << ','
               << tmp.M2 << ',' << tmp.M_abs << '\n';
@@ -80,11 +48,12 @@ void monte_carlo_progression(double T, uint L, uint cycles,
     ofile.close();
 }
 
-void pd_estimate(double T, uint L, uint cycles, const std::string filename)
+void monte_carlo_progression(double T, int L, int cycles, int value,
+                             const std::string filename)
 {
     // Set some variables
-    data_t data;
-    uint n_spins = L * L;
+    data_t data, tmp;
+    int n_spins = L * L;
 
     // File stuff
     std::string directory = utils::dirname(filename);
@@ -92,104 +61,126 @@ void pd_estimate(double T, uint L, uint cycles, const std::string filename)
 
     // Create random engine using the mersenne twister
     std::random_device rd;
+    uint32_t rd_sample = rd();
     std::mt19937 engine(rd());
 
+    std::cout << "Seed: " << rd_sample << std::endl;
+
+    IsingModel ising(L, T, value);
+
+    // Create path and open file
+    utils::mkpath(directory);
+    ofile.open(filename);
+
+    // Loop through cycles
+    for (size_t i = 1; i <= cycles; i++) {
+        data += ising.Metropolis();
+        tmp = data / (i * n_spins);
+        ofile << i << ',' << tmp.E << ',' << tmp.E2 << ',' << tmp.M << ','
+              << tmp.M2 << ',' << tmp.M_abs << '\n';
+    }
+    ofile.close();
+}
+
+void pd_estimate(double T, int L, int cycles, const std::string filename)
+{
+    // Set some variables
+    data_t data, tmp;
+    int n_spins = L * L;
+
+    // File stuff
+    std::string directory = utils::dirname(filename);
+    std::ofstream ofile;
+
     IsingModel ising(L, T);
 
     // Create path and open file
     utils::mkpath(directory);
     ofile.open(filename);
 
-    double E, M;
-
-    // Figure out bin widths and such
+    // Loop through cycles
+    for (size_t i = 1; i <= cycles; i++) {
+        data = ising.Metropolis() / n_spins;
+        ofile << data.E << ',' << data.E2 << ',' << data.M << ',' << data.M2
+              << ',' << data.M_abs << '\n';
+    }
+    ofile.close();
 }
 
 // Code for seeing phase transitions.
-void monte_carlo_serial(data_t &data, uint L, double T, uint cycles)
+data_t monte_carlo_serial(int L, double T, int cycles)
 {
+    data_t data;
     IsingModel model(L, T);
-    // Create random engine using the mersenne twister
-    std::random_device rd;
-    std::mt19937 engine(rd());
 
     for (size_t i = 0; i < BURN_IN_TIME; i++) {
-        model.Metropolis(engine);
+        model.Metropolis();
     }
 
-    for (size_t i = 0; i < cycles; i++) {
-        data += model.Metropolis(engine);
+    double E, M;
+    for (size_t i = 0; i < (uint)cycles; i++) {
+        data += model.Metropolis();
     }
-    data /= cycles;
+
+    return data;
 }
 
-void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles)
+data_t monte_carlo_parallel(int L, double T, int cycles)
 {
+    data_t data;
 #pragma omp parallel
     {
         // Each thread creates an instance of IsingModel.
         IsingModel model(L, T);
-        // Each thread creates an instance of the mersenne twister
-        std::random_device rd;
-        std::mt19937 engine(rd());
-
-        data_t tmp;
 
         // Each thread runs the Metropolis algorithm before starting to collect
         // samples
         for (size_t i = 0; i < BURN_IN_TIME; i++) {
-            model.Metropolis(engine);
+            model.Metropolis();
         }
 
         // Now each thread work on one loop together, but using their own
         // instances of things, but the total of cycles add up.
         // static ensure that each thread gets the same amount of iterations
-#pragma omp for schedule(static)
-        for (size_t i = 0; i < cycles; i++) {
-            tmp = tmp + model.Metropolis(engine);
-        }
-
-        // Combine all the data.
-#pragma omp critical
-        {
-            data += tmp;
+#pragma omp for schedule(static) reduction(+ : data)
+        for (size_t i = 0; i < (uint)cycles; i++) {
+            data += model.Metropolis();
         }
     }
 
-    data /= cycles;
+    double norm = 1. / (double)cycles;
+
+    return data * norm;
 }
 
-void phase_transition(
-    uint L, double start_T, double end_T, uint points_T,
-    std::function<void(data_t &, uint, double, uint)> monte_carlo,
-    std::string outfile)
+void phase_transition(int L, double start, double end, int points,
+                      std::function<data_t(int, double, int)> monte_carlo,
+                      std::string outfile)
 {
-    double dt_T = (end_T - start_T) / points_T;
-    uint cycles = 10000;
-    uint N = L * L;
+    double dt = (end - start) / (double)points;
+    int cycles = 10000;
+    int N = L * L;
     std::ofstream ofile;
 
-    data_t data[points_T];
-
-    for (size_t i = 0; i < points_T; i++) {
-        monte_carlo(data[i], L, start_T + dt_T * i, cycles);
-    }
+    data_t data;
 
     utils::mkpath(utils::dirname(outfile));
     ofile.open(outfile);
 
-    double temp, CV, X;
+    double temp, CV, X, E_var, M_var;
 
     using utils::scientific_format;
-    for (size_t i = 0; i < points_T; i++) {
-        temp = start_T + dt_T * i;
-        CV = (data[i].E2 - data[i].E * data[i].E) / (N * temp * temp);
-        X = (data[i].M2 - data[i].M_abs * data[i].M_abs) / (N * temp);
+    for (size_t i = 0; i < points; i++) {
+        temp = start + dt * i;
+        data = monte_carlo(L, temp, cycles);
+        E_var = (data.E2 - data.E * data.E) / (double)N;
+        M_var = (data.M2 - data.M_abs * data.M_abs) / (double)N;
 
         ofile << scientific_format(temp) << ','
-              << scientific_format(data[i].E / N) << ','
-              << scientific_format(data[i].M_abs / N) << ','
-              << scientific_format(CV) << ',' << scientific_format(X) << '\n';
+              << scientific_format(data.E / (double)N) << ','
+              << scientific_format(data.M_abs / N) << ','
+              << scientific_format(E_var / (temp * temp)) << ','
+              << scientific_format(M_var / temp) << '\n';
     }
     ofile.close();
 }
diff --git a/src/phase_transition_mpi.cpp b/src/phase_transition_mpi.cpp
index ac925a0..ce197f6 100644
--- a/src/phase_transition_mpi.cpp
+++ b/src/phase_transition_mpi.cpp
@@ -18,15 +18,22 @@
 #include <iostream>
 #include <iterator>
 #include <mpi.h>
+#include <sstream>
 
 /** @brief The main function*/
 int main(int argc, char **argv)
 {
-    double start = 1., end = 3.;
-    uint points = 1000, L = 20, N;
+    if (argc < 5) {
+        std::cout << "You need at least 4 arguments" << std::endl;
+        abort();
+    }
+    double t0, t1;
+    t0 = MPI_Wtime();
+    double start = atof(argv[1]), end = atof(argv[2]);
+    int points = atoi(argv[3]), N;
+    int lattice_sizes[] = {20, 40, 60, 80, 100};
     double dt = (end - start) / points;
-    uint cycles = 10000;
-    N = L * L;
+    int cycles = atoi(argv[4]);
     std::ofstream ofile;
 
     data_t data[points];
@@ -41,10 +48,10 @@ int main(int argc, char **argv)
     MPI_Comm_size(MPI_COMM_WORLD, &cluster_size);
     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
-    uint remainder = points % cluster_size;
+    int remainder = points % cluster_size;
     double i_start;
-    uint i_points;
-    // The last
+    int i_points;
+    // Distribute temperature points
     if (rank < remainder) {
         i_points = points / cluster_size + 1;
         i_start = start + dt * i_points * rank;
@@ -57,52 +64,65 @@ int main(int argc, char **argv)
     data_t i_data[i_points];
     std::cout << "Rank " << rank << ": " << i_points << ',' << i_start << '\n';
 
-    for (size_t i = 0; i < i_points; i++) {
-        monte_carlo_serial(i_data[i], L, i_start + dt * i, cycles);
-    }
+    for (int L : lattice_sizes) {
+        N = L * L;
+        for (size_t i = 0; i < i_points; i++) {
+            i_data[i] = monte_carlo_parallel(L, i_start + dt * i, cycles);
+        }
 
-    if (rank == 0) {
-        std::copy_n(i_data, i_points, data);
-        for (size_t i = 1; i < cluster_size; i++) {
-            if (rank < remainder) {
-                MPI_Recv((void *)i_data,
-                         sizeof(data_t) * (points / cluster_size + 1), MPI_CHAR,
-                         i, MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-                std::copy_n(i_data, points / cluster_size + 1,
-                            data + (points / cluster_size) * i);
+        if (rank == 0) {
+            std::copy_n(i_data, i_points, data);
+            for (size_t i = 1; i < cluster_size; i++) {
+                if (rank < remainder) {
+                    MPI_Recv((void *)i_data,
+                             sizeof(data_t) * (points / cluster_size + 1),
+                             MPI_CHAR, i, MPI_ANY_TAG, MPI_COMM_WORLD,
+                             MPI_STATUS_IGNORE);
+                    std::copy_n(i_data, points / cluster_size + 1,
+                                data + (points / cluster_size) * i);
+                }
+                else {
+                    MPI_Recv((void *)i_data,
+                             sizeof(data_t) * (points / cluster_size), MPI_CHAR,
+                             i, MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                    std::copy_n(i_data, points / cluster_size,
+                                data + (points / cluster_size) * i + remainder);
+                }
             }
-            else {
-                MPI_Recv((void *)i_data,
-                         sizeof(data_t) * (points / cluster_size), MPI_CHAR, i,
-                         MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
-                std::copy_n(i_data, points / cluster_size,
-                            data + (points / cluster_size) * i + remainder);
+            std::stringstream outfile;
+            outfile << "output/phase_transition/size_" << L << ".txt";
+            utils::mkpath(utils::dirname(outfile.str()));
+            ofile.open(outfile.str());
+
+            double temp, CV, X;
+
+            using utils::scientific_format;
+            for (size_t i = 0; i < points; i++) {
+                temp = start + dt * i;
+                CV = (data[i].E2 - data[i].E * data[i].E)
+                     / ((double)N * temp * temp);
+                X = (data[i].M2 - data[i].M_abs * data[i].M_abs)
+                    / ((double)N * temp);
+
+                ofile << scientific_format(temp) << ','
+                      << scientific_format(data[i].E / N) << ','
+                      << scientific_format(data[i].M_abs / N) << ','
+                      << scientific_format(CV) << ',' << scientific_format(X)
+                      << '\n';
             }
+            ofile.close();
+        }
+        else {
+            MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
+                     MPI_COMM_WORLD);
         }
     }
-    else {
-        MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
-                 MPI_COMM_WORLD);
+
+    t1 = MPI_Wtime();
+
+    if (rank == 0) {
+        std::cout << "Time: " << t1 - t0 << " seconds\n";
     }
 
     MPI_Finalize();
-
-    std::string outfile = "output/phase_transition/size_20.txt";
-    utils::mkpath(utils::dirname(outfile));
-    ofile.open(outfile);
-
-    double temp, CV, X;
-
-    using utils::scientific_format;
-    for (size_t i = 0; i < points; i++) {
-        temp = start + dt * i;
-        CV = (data[i].E2 - data[i].E * data[i].E) / (N * temp * temp);
-        X = (data[i].M2 - data[i].M_abs * data[i].M_abs) / (N * temp);
-
-        ofile << scientific_format(temp) << ','
-              << scientific_format(data[i].E / N) << ','
-              << scientific_format(data[i].M_abs / N) << ','
-              << scientific_format(CV) << ',' << scientific_format(X) << '\n';
-    }
-    ofile.close();
 }
diff --git a/src/test_suite.cpp b/src/test_suite.cpp
index 59ec82c..6305451 100644
--- a/src/test_suite.cpp
+++ b/src/test_suite.cpp
@@ -12,6 +12,19 @@
 #include "IsingModel.hpp"
 #include "testlib.hpp"
 
+#include <fstream>
+
+#define EPS_2 (-2 * std::sinh(8.)) / (std::cosh(8.) + 3)
+
+#define MAG_2 (std::exp(8.) + 1) / (2 * (cosh(8.) + 3))
+
+#define CV_2                                                                   \
+    16 * (3 * std::cosh(8.) + 1) / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))
+
+#define X_2                                                                    \
+    (3 * std::exp(8.) + std::exp(-8.) + 3)                                     \
+        / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))
+
 /** @brief Test class for the Ising model
  * */
 class IsingModelTest {
@@ -31,7 +44,7 @@ public:
                "Test lattice initialization.");
 
         test.initialize_neighbors();
-        arma::Mat<uint> neighbor_matrix("2, 1 ; 0, 2 ; 1, 0");
+        arma::Mat<int> neighbor_matrix("2, 1 ; 0, 2 ; 1, 0");
         ASSERT(testlib::is_equal(neighbor_matrix, test.neighbors),
                "Test neighbor matrix.");
 
@@ -40,11 +53,94 @@ public:
 
         // Test the initial magnetization.
         test.initialize_magnetization();
-        ASSERT(std::abs(test.M - 9.) < 1e-8, "Test intial magnetization");
+        ASSERT(std::fabs(test.M - 9.) < 1e-8, "Test intial magnetization");
 
         // Test that the initial energy is correct
         test.initialize_energy();
-        ASSERT(std::abs(test.E - (-18)) < 1e-8, "Test initial energy.");
+        ASSERT(std::fabs(test.E - (-18)) < 1e-8, "Test initial energy.");
+    }
+
+    /** @brief Test numerical data with analytical data.
+     *
+     *  @param tol The tolerance between the analytical and numerical solution.
+     *  @param max_cycles The max number of Monte Carlo cycles.
+     *
+     *  return int
+     * */
+    int test_2x2_lattice(double tol, int max_cycles)
+    {
+        data_t data, tmp;
+        size_t L = 2;
+        size_t n_spins = L * L;
+        double T = 1.;
+        size_t cycles = 0;
+
+        // Create random engine using the mersenne twister
+        std::random_device rd;
+        std::mt19937 engine(rd());
+
+        IsingModel test(L, T);
+
+        int arr[]{0, 0, 0, 0};
+
+        // Loop through cycles
+        //std::ofstream ofile;
+        //ofile.open("output/test_2x2.txt");
+        while (cycles++ < max_cycles) {
+            data += test.Metropolis();
+            tmp = data / cycles;
+            //ofile << cycles << ',' << tmp.E / n_spins << ','
+                  //<< tmp.M_abs / n_spins << ','
+                  //<< (tmp.E2 - tmp.E * tmp.E) / (T * T) / n_spins << ','
+                  //<< (tmp.M2 - tmp.M_abs * tmp.M_abs) / T / n_spins << '\n';
+            if (testlib::close_to(EPS_2, tmp.E / n_spins, tol)
+                && testlib::close_to(MAG_2, tmp.M_abs / n_spins, tol)
+                && testlib::close_to(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T)
+                                               / n_spins, tol)
+                && testlib::close_to(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T
+                                              / n_spins, tol)) {
+                return cycles;
+            }
+        }
+        //std::cout << EPS_2 << ',' << MAG_2 << ',' << CV_2 << ',' << X_2
+                  //<< std::endl;
+        //ofile.close();
+        // cycles = 0;
+        // data = 0;
+        // IsingModel test_mag(L, T);
+        // while (cycles++ < max_cycles) {
+        // data += test.Metropolis();
+        // tmp = data / (cycles * n_spins);
+        // if (testlib::close_to(MAG_2, tmp.M, tol)) {
+        // arr[1] = cycles;
+        // break;
+        //}
+        //}
+        // cycles = 0;
+        // data = 0;
+        // IsingModel test_CV(L, T);
+        // while (cycles++ < max_cycles) {
+        // data += test.Metropolis();
+        // tmp = data / (cycles * n_spins);
+        // if (testlib::close_to(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T),
+        // tol)) {
+        // arr[2] = cycles;
+        // break;
+        //}
+        //}
+        // cycles = 0;
+        // data = 0;
+        // IsingModel test_X(L, T);
+        // while (cycles++ < max_cycles) {
+        // data += test.Metropolis();
+        // tmp = data / (cycles * n_spins);
+        // if (testlib::close_to(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T,
+        // tol)) {
+        // arr[3] = cycles;
+        // break;
+        //}
+        //}
+        return 0;
     }
 };
 
@@ -54,4 +150,18 @@ int main()
     IsingModelTest test;
 
     test.test_init_functions();
+    int res = 0;
+    int tmp;
+    for (size_t i=0; i < 1000; i++) {
+        tmp = test.test_2x2_lattice(1e-2, 1e5);
+        if (tmp == 0) {
+            std::cout << "not enough cycles\n";
+            break;
+        }
+        res += tmp;
+    }
+    
+    std::cout << "Res: " << res / 1000 << std::endl;
+
+    return 0;
 }