diff --git a/include/IsingModel.hpp b/include/IsingModel.hpp
index 8ec4ae9..aac9df3 100644
--- a/include/IsingModel.hpp
+++ b/include/IsingModel.hpp
@@ -18,6 +18,7 @@
 #include "utils.hpp"
 
 #include <armadillo>
+#include <cstdint>
 #include <random>
 #include <unordered_map>
 
@@ -53,7 +54,8 @@ private:
 
     /** @brief A hash map containing all possible energy changes.
      * */
-    std::unordered_map<int, double> energy_diff;
+    //std::unordered_map<int, double> energy_diff;
+    double energy_diff[17];
 
     /** @brief The temperature of the model.
      * */
@@ -65,12 +67,15 @@ private:
 
     /** @brief The current energy state. unit: \f$ J \f$.
      * */
-    int E;
+    int64_t E;
 
     /** @brief The current magnetic strength. unit: Unitless.
      * */
-    int M;
+    int64_t M;
 
+    std::mt19937 engine;
+
+    void initialize_engine();
     /** @brief Initialize the lattice with a random distribution of 1s and
      *  -1s.
      * */
diff --git a/include/data_type.hpp b/include/data_type.hpp
index 5d3bfd5..7192660 100644
--- a/include/data_type.hpp
+++ b/include/data_type.hpp
@@ -49,7 +49,7 @@ public:
         return res;
     }
 
-    template <class T> data_t& operator/=(T num)
+    template <class T> data_t &operator/=(T num)
     {
         this->E /= (double)num;
         this->E2 /= (double)num;
@@ -72,7 +72,7 @@ public:
         return res;
     }
 
-    template <class T> data_t& operator*=(T num)
+    template <class T> data_t &operator*=(T num)
     {
         this->E *= (double)num;
         this->E2 *= (double)num;
@@ -95,7 +95,7 @@ public:
         return res;
     }
 
-    data_t& operator+=(const data_t &b)
+    data_t &operator+=(const data_t &b)
     {
         this->E += b.E;
         this->E2 += b.E2;
@@ -116,4 +116,6 @@ public:
     }
 };
 
+#pragma omp declare reduction(+ : data_t : omp_out += omp_in)
+
 #endif
diff --git a/include/monte_carlo.hpp b/include/monte_carlo.hpp
index 36e7b05..6fc2544 100644
--- a/include/monte_carlo.hpp
+++ b/include/monte_carlo.hpp
@@ -5,7 +5,7 @@
  *
  *  @version 1.0
  *
- *  @brief Functions for monte carlo simulations.
+ *  @brief Functions for Monte Carlo simulations.
  *
  *  @bug No known bugs
  * */
@@ -17,23 +17,13 @@
 #include "utils.hpp"
 
 #include <functional>
-#include <string>
 #include <omp.h>
+#include <string>
 
-//#define BURN_IN_TIME 12500
+// #define BURN_IN_TIME 12500
 #define BURN_IN_TIME 5000
 
-#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 int
- * */
-int test_2x2_lattice(double tol, int max_cycles);
-
+namespace montecarlo {
 /** @brief Write the expected values for each Monte Carlo cycles to file.
  *
  *  @param T Temperature
@@ -41,8 +31,8 @@ int test_2x2_lattice(double tol, int 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, int L, int cycles,
-                             const std::string filename);
+void progression(double T, int L, int cycles,
+                                    const std::string filename);
 
 /** @brief Write the expected values for each Monte Carlo cycles to file.
  *
@@ -52,8 +42,8 @@ void monte_carlo_progression(double T, int L, int cycles,
  *  @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);
+void progression(double T, int L, int cycles, int value,
+                                    const std::string filename);
 
 /** @brief Estimate the probability distribution for the energy.
  *
@@ -62,7 +52,8 @@ void monte_carlo_progression(double T, int L, int cycles, int value,
  *  @param cycles The amount of Monte Carlo cycles to do
  *  @param filename The file to write to
  * */
-void pd_estimate(double T, int L, int 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
  *  Carlo cycles.
@@ -73,7 +64,7 @@ void pd_estimate(double T, int L, int cycles, const std::string filename);
  *
  *  @return data_t
  * */
-data_t monte_carlo_serial(int L, double T, int cycles);
+data_t mcmc_serial(int L, double T, int cycles, int burn_in_time = BURN_IN_TIME);
 
 /** @brief Execute the Metropolis algorithm for a certain amount of Monte
  *  Carlo cycles in parallel.
@@ -84,7 +75,7 @@ data_t monte_carlo_serial(int L, double T, int cycles);
  *
  *  @return data_t
  * */
-data_t monte_carlo_parallel(int L, double T, int cycles);
+data_t mcmc_parallel(int L, double T, int cycles, int burn_in_time = BURN_IN_TIME);
 
 /** @brief Perform the MCMC algorithm using a range of temperatures.
  *
@@ -95,9 +86,10 @@ data_t monte_carlo_parallel(int L, double T, int cycles);
  *  @param monte_carlo Which Monte Carlo implementation to use
  *  @param outfile The file to write the data to
  * */
-void phase_transition(
-    int L, double start_T, double end_T, int points_T,
-    std::function<data_t(int, double, int)> monte_carlo,
-    std::string outfile);
+void
+phase_transition(int L, double start_T, double end_T, int points_T, int cycles,
+                 std::function<data_t(int, double, int, int)> monte_carlo,
+                 std::string outfile, int burn_in_time = BURN_IN_TIME);
+}; // namespace montecarlo
 
 #endif
diff --git a/slurm_scripts/execute.script b/slurm_scripts/execute.script
new file mode 100755
index 0000000..ae29b33
--- /dev/null
+++ b/slurm_scripts/execute.script
@@ -0,0 +1,64 @@
+#!/bin/bash
+
+usage() {
+>&2 cat << EOF
+Usage: $0
+   [ -h | --help ]
+   [ --start-temp input ]
+   [ --end-temp input ] 
+   [ --points input ]
+   [ --samples input ]
+EOF
+exit 1
+}
+
+# Defaults
+start_temp=2.1
+end_temp=2.4
+points_temp=40
+samples=1000000
+array_arg=20
+time_arg="0-00:30:00"
+
+VALID_ARGS=$(getopt -o h --long help,start-temp:,end-temp:,points:,samples:,array:,time: -- "$@")
+if [[ $? -ne 0 ]]; then
+    usage
+fi
+
+eval set -- ${VALID_ARGS}
+while : 
+do
+    case "$1" in
+        -h | --help)
+            usage
+            shift
+            ;;
+        --start-temp)
+            start_temp=$2
+            shift 2
+            ;;
+        --end-temp)
+            end_temp=$2
+            shift 2
+            ;;
+        --points)
+            points=$2
+            shift 2
+            ;;
+        --samples)
+            samples=$2
+            shift 2
+            ;;
+        --array)
+            array_arg=$(echo "${2// /}")
+            shift 2
+            ;;
+        --time)
+            time_arg=$2
+            shift 2
+            ;;
+        --) shift; break ;;
+    esac
+done
+
+sbatch --array=$array_arg --time=$time_arg ./jobs/pt.script $start_temp $end_temp $points_temp $samples
diff --git a/slurm_scripts/pt.script b/slurm_scripts/pt.script
new file mode 100755
index 0000000..db0d0d6
--- /dev/null
+++ b/slurm_scripts/pt.script
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+#SBATCH --account=ec54
+#SBATCH --job-name=pt
+#SBATCH --ntasks=10
+#SBATCH --mem-per-cpu=1G
+#SBATCH --cpus-per-task=10
+
+set -o errexit # Exit the script on any error
+set -o nounset # Treat any unset variables as an error
+
+module --quiet purge # Reset the modules to the system default
+module load Armadillo/11.4.3-foss-2022b
+module load OpenMPI/4.1.5-GCC-12.3.0
+
+# Args
+start_temp=$1
+end_temp=$2
+points_temp=$3
+samples=$4
+
+srun ./phase_transition_mpi $start_temp $end_temp $points_temp ${SLURM_ARRAY_TASK_ID} $samples 0
diff --git a/slurm_scripts/pt_narrow.script b/slurm_scripts/pt_narrow.script
new file mode 100644
index 0000000..d4aec22
--- /dev/null
+++ b/slurm_scripts/pt_narrow.script
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+#SBATCH --account=ec54
+#SBATCH --job-name=pt_narrow
+#SBATCH --time=0-02:00:00
+#SBATCH --ntasks=8
+#SBATCH --mem-per-cpu=1G
+#SBATCH --cpus-per-task=10
+
+set -o errexit # Exit the script on any error
+set -o nounset # Treat any unset variables as an error
+
+module --quiet purge # Reset the modules to the system default
+module load Armadillo/11.4.3-foss-2022b
+module load OpenMPI/4.1.5-GCC-12.3.0
+
+srun ./phase_transition_mpi 2.25 2.35 40 10000000
diff --git a/src/IsingModel.cpp b/src/IsingModel.cpp
index e0a600e..c94478e 100644
--- a/src/IsingModel.cpp
+++ b/src/IsingModel.cpp
@@ -10,18 +10,18 @@
  *  @bug No known bugs
  * */
 #include "IsingModel.hpp"
-
-#include <cmath>
 #include <random>
 
 IsingModel::IsingModel()
 {
+    this->initialize_engine();
 }
 
 IsingModel::IsingModel(int L, double T)
 {
     this->L = L;
     this->T = T;
+    this->initialize_engine();
     this->initialize_lattice();
     this->initialize_neighbors();
     this->initialize_energy_diff();
@@ -33,6 +33,7 @@ IsingModel::IsingModel(int L, double T, int val)
 {
     this->L = L;
     this->T = T;
+    this->initialize_engine();
     this->lattice.set_size(this->L, this->L);
     this->lattice.fill(val);
     this->initialize_neighbors();
@@ -41,16 +42,20 @@ IsingModel::IsingModel(int L, double T, int val)
     this->initialize_energy();
 }
 
+void IsingModel::initialize_engine()
+{
+    std::random_device rd{};
+    this->engine = std::mt19937{rd()};
+}
+
 void IsingModel::initialize_lattice()
 {
     this->lattice.set_size(this->L, this->L);
-    std::random_device rd{};
-    std::mt19937 engine{rd()};
 
     std::uniform_int_distribution<> coin_flip(0, 1);
 
     for (size_t i = 0; i < this->lattice.n_elem; i++)
-        this->lattice(i) = 2 * coin_flip(engine) - 1;
+        this->lattice(i) = 2 * coin_flip(this->engine) - 1;
 }
 
 void IsingModel::initialize_neighbors()
@@ -67,7 +72,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[i+8] = std::exp(-(double)i / this->T);
     }
 }
 
@@ -95,9 +100,6 @@ void IsingModel::initialize_energy()
 
 data_t IsingModel::Metropolis()
 {
-    std::random_device rd{};
-    std::mt19937_64 engine{rd()};
-
     int ri, rj;
     int dE;
 
@@ -120,7 +122,7 @@ data_t IsingModel::Metropolis()
                 + this->lattice(this->neighbors(ri, DOWN), rj));
 
         // Choose whether or not to accept the new configuration
-        if (random_number(engine) <= this->energy_diff[dE]) {
+        if (random_number(engine) <= this->energy_diff[dE+8]) {
             // Update if the configuration is accepted
             this->lattice(ri, rj) *= -1;
             this->M += 2 * this->lattice(ri, rj);
diff --git a/src/Makefile b/src/Makefile
index 720b48c..540efaa 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -3,9 +3,11 @@ CC=mpic++
 
 LIBSRCS=utils.cpp testlib.cpp data_type.cpp
 LIBOBJS=$(LIBSRCS:.cpp=.o)
+LIBPROFOBJS=$(addprefix prof/, $(LIBOBJS))
 
 CLASSSRCS=IsingModel.cpp monte_carlo.cpp
 CLASSOBJS=$(CLASSSRCS:.cpp=.o)
+CLASSPROFOBJS=$(addprefix prof/, $(CLASSOBJS))
 
 INCLUDE=../include
 
@@ -29,16 +31,14 @@ else
 	PROFFLAG=
 endif
 
-.PHONY: clean
+.PHONY: clean instrument
 
-all: main phase_transition_mpi test_suite 
+all: main phase_transition_mpi test_suite  time
 #all: main 
 
 # Instrumentation using scorep for parallel analysis
-instrument:
-	scorep $(CC) -c utils.cpp -o utils.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
-	scorep $(CC) -c main.cpp -o main.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
-	scorep $(CC) $(LIBOBJS) $(CLASSOBJS) main.o -o main $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
+instrument: prof/phase_transition_mpi.o $(LIBPROFOBJS) $(CLASSPROFOBJS)
+	scorep $(CC) $(LIBPROFOBJS) $(CLASSPROFOBJS) $< -o phase_transition_mpi_prof $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
 	
 # Rules for executables
 main: main.o $(LIBOBJS) $(CLASSOBJS)
@@ -50,10 +50,17 @@ phase_transition_mpi: phase_transition_mpi.o $(LIBOBJS) $(CLASSOBJS)
 test_suite: test_suite.o $(LIBOBJS) $(CLASSOBJS)
 	$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
 
+time: time.o $(LIBOBJS) $(CLASSOBJS)
+	$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
+
 # Rule for object files
 %.o: %.cpp
 	$(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
 
+# Rule for instrumented object files
+prof/%.o: %.cpp
+	scorep $(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
+
 clean:
-	rm *.o
-	rm test_suite main 
+	find . -maxdepth 2 -name "*.o" -type f -delete
+	rm test_suite main phase_transition_mpi phase_transition_mpi_prof time
diff --git a/src/main.cpp b/src/main.cpp
index ce838a6..4d9c1d0 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -13,31 +13,43 @@
 #include "monte_carlo.hpp"
 #include "utils.hpp"
 
-#include <csignal>
-#include <cstdlib>
-#include <iostream>
-#include <omp.h>
-
+/** @brief Create the data for the burn-in time for temperatures 1.0 and 2.4
+ *  for both unordered and ordered initial states.
+ * */
 void create_burn_in_time_data()
 {
     // Test burn-in time
-    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");
+    montecarlo::progression(1.0, 20, 20000,
+                            "../output/burn_in_time/unordered_1_0.txt");
+    montecarlo::progression(1.0, 20, 20000, 1,
+                            "../output/burn_in_time/ordered_1_0.txt");
+    montecarlo::progression(2.4, 20, 20000,
+                            "../output/burn_in_time/unordered_2_4.txt");
+    montecarlo::progression(2.4, 20, 20000, 1,
+                            "../output/burn_in_time/ordered_2_4.txt");
 }
 
+/** @brief Create the data used to estimate the probability distribution
+ *  for tempratures 1.0 anbd 2.4.
+ * */
 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");
+    montecarlo::pd_estimate(1.0, 20, 1000000,
+                            "../output/pd_estimate/estimate_1_0.txt");
+    montecarlo::pd_estimate(2.4, 20, 1000000,
+                            "../output/pd_estimate/estimate_2_4.txt");
 }
 
+void test_burn_in_time()
+{
+    montecarlo::phase_transition(100, 2.1, 2.4, 40, 1e5, montecarlo::mcmc_serial,
+                                 "../output/test_burn_in_time/no_burn_in.txt", 0);
+    montecarlo::phase_transition(100, 2.1, 2.4, 40, 1e5, montecarlo::mcmc_serial,
+                                 "../output/test_burn_in_time/burn_in.txt", 5000);
+}
+
+/** @brief Test how much Openmp speeds up.*/
 void test_parallel_speedup()
 {
     // Test the openmp speedup
@@ -46,10 +58,10 @@ void test_parallel_speedup()
     int tries = 5;
     t0 = omp_get_wtime();
     for (size_t i = 0; i < tries; i++)
-        monte_carlo_serial(20, 1.0, 10000);
+        montecarlo::mcmc_serial(20, 1.0, 10000);
     t1 = omp_get_wtime();
     for (size_t i = 0; i < tries; i++)
-        monte_carlo_parallel(20, 1.0, 10000);
+        montecarlo::mcmc_parallel(20, 1.0, 10000);
     t2 = omp_get_wtime();
 
     std::cout << "Time serial     : " << (t1 - t0) / tries << " seconds"
@@ -59,22 +71,29 @@ void test_parallel_speedup()
     std::cout << "Speedup parallel: " << (t1 - t0) / (t2 - t1) << '\n';
 }
 
+/** @brief Create data for studying phase transition.
+ * */
 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");
+    montecarlo::phase_transition(20, 2.1, 2.4, 40, 1e4,
+                                 montecarlo::mcmc_parallel,
+                                 "../output/phase_transition/size_20.txt");
+    montecarlo::phase_transition(40, 2.1, 2.4, 40, 1e4,
+                                 montecarlo::mcmc_parallel,
+                                 "../output/phase_transition/size_40.txt");
+    montecarlo::phase_transition(60, 2.1, 2.4, 40, 1e4,
+                                 montecarlo::mcmc_parallel,
+                                 "../output/phase_transition/size_60.txt");
+    montecarlo::phase_transition(80, 2.1, 2.4, 40, 1e4,
+                                 montecarlo::mcmc_parallel,
+                                 "../output/phase_transition/size_80.txt");
+    montecarlo::phase_transition(100, 2.1, 2.4, 40, 1e4,
+                                 montecarlo::mcmc_parallel,
+                                 "../output/phase_transition/size_100.txt");
     t1 = omp_get_wtime();
 
     std::cout << "Time: " << t1 - t0 << std::endl;
@@ -104,6 +123,9 @@ int main(int argc, char **argv)
     case 4:
         create_phase_transition_data();
         break;
+    case 5:
+        test_burn_in_time();
+        break;
     default:
         std::cout << "Not a valid option!" << std::endl;
         abort();
diff --git a/src/monte_carlo.cpp b/src/monte_carlo.cpp
index 4f64092..efe5113 100644
--- a/src/monte_carlo.cpp
+++ b/src/monte_carlo.cpp
@@ -11,11 +11,8 @@
  * */
 #include "monte_carlo.hpp"
 
-#include <cmath>
-#include <cstdint>
-
-void monte_carlo_progression(double T, int L, int cycles,
-                             const std::string filename)
+namespace montecarlo {
+void progression(double T, int L, int cycles, const std::string filename)
 {
     // Set some variables
     data_t data, tmp;
@@ -48,8 +45,8 @@ void monte_carlo_progression(double T, int L, int cycles,
     ofile.close();
 }
 
-void monte_carlo_progression(double T, int L, int cycles, int value,
-                             const std::string filename)
+void progression(double T, int L, int cycles, int value,
+                 const std::string filename)
 {
     // Set some variables
     data_t data, tmp;
@@ -108,7 +105,7 @@ void pd_estimate(double T, int L, int cycles, const std::string filename)
 }
 
 // Code for seeing phase transitions.
-data_t monte_carlo_serial(int L, double T, int cycles)
+data_t mcmc_serial(int L, double T, int cycles, int burn_in_time)
 {
     data_t data;
     IsingModel model(L, T);
@@ -122,10 +119,12 @@ data_t monte_carlo_serial(int L, double T, int cycles)
         data += model.Metropolis();
     }
 
-    return data;
+    double norm = 1. / (double)cycles;
+
+    return data * norm;
 }
 
-data_t monte_carlo_parallel(int L, double T, int cycles)
+data_t mcmc_parallel(int L, double T, int cycles, int burn_in_time)
 {
     data_t data;
 #pragma omp parallel
@@ -153,12 +152,11 @@ data_t monte_carlo_parallel(int L, double T, int cycles)
     return data * norm;
 }
 
-void phase_transition(int L, double start, double end, int points,
-                      std::function<data_t(int, double, int)> monte_carlo,
-                      std::string outfile)
+void phase_transition(int L, double start, double end, int points, int cycles,
+                      std::function<data_t(int, double, int, int)> monte_carlo,
+                      std::string outfile, int burn_in_time)
 {
     double dt = (end - start) / (double)points;
-    int cycles = 10000;
     int N = L * L;
     std::ofstream ofile;
 
@@ -172,7 +170,7 @@ void phase_transition(int L, double start, double end, int points,
     using utils::scientific_format;
     for (size_t i = 0; i < points; i++) {
         temp = start + dt * i;
-        data = monte_carlo(L, temp, cycles);
+        data = monte_carlo(L, temp, cycles, burn_in_time);
         E_var = (data.E2 - data.E * data.E) / (double)N;
         M_var = (data.M2 - data.M_abs * data.M_abs) / (double)N;
 
@@ -184,3 +182,4 @@ void phase_transition(int L, double start, double end, int points,
     }
     ofile.close();
 }
+} // namespace montecarlo
diff --git a/src/phase_transition_mpi.cpp b/src/phase_transition_mpi.cpp
index ce197f6..720e672 100644
--- a/src/phase_transition_mpi.cpp
+++ b/src/phase_transition_mpi.cpp
@@ -7,38 +7,42 @@
  *
  *  @brief Sweep over different temperatures and generate data.
  *
+ *  @details This program takes in 4 arguments: the start temperature,
+ *  the end temperature, the amount of temperature points to simulate, and
+ *  the amount of monte carlo samples to collect, in that order.
+ *
  *  @bug No known bugs
  * */
 #include "data_type.hpp"
 #include "monte_carlo.hpp"
 #include "utils.hpp"
 
-#include <algorithm>
-#include <fstream>
-#include <iostream>
-#include <iterator>
 #include <mpi.h>
-#include <sstream>
 
-/** @brief The main function*/
+/** @brief The main function
+ *
+ *  */
 int main(int argc, char **argv)
 {
-    if (argc < 5) {
-        std::cout << "You need at least 4 arguments" << std::endl;
+    // Check that the number of arguments is at least 4.
+    if (argc < 7) {
+        std::cout << "You need at least 6 arguments" << std::endl;
         abort();
     }
+
+    // Timing variables
     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;
-    int cycles = atoi(argv[4]);
-    std::ofstream ofile;
 
+    // Define/initialize variables
+    double start = atof(argv[1]), end = atof(argv[2]);
+    int points = atoi(argv[3]), cycles = atoi(argv[5]), L = atoi(argv[4]),
+        burn_in_time = atoi(argv[6]), N = L * L;
+    double dt = (end - start) / points;
+    std::ofstream ofile;
     data_t data[points];
 
-    // MPI stuff
+    // MPI specific variables
     int rank, cluster_size;
 
     // Initialize MPI
@@ -49,8 +53,9 @@ int main(int argc, char **argv)
     MPI_Comm_rank(MPI_COMM_WORLD, &rank);
 
     int remainder = points % cluster_size;
-    double i_start;
-    int i_points;
+    double i_start; // What temperature to start from
+    int i_points;   // How many points to simulate
+
     // Distribute temperature points
     if (rank < remainder) {
         i_points = points / cluster_size + 1;
@@ -61,61 +66,67 @@ int main(int argc, char **argv)
         i_start = start + dt * (i_points * rank + remainder);
     }
 
+    // Initialize array to contains data for each temperature point
     data_t i_data[i_points];
-    std::cout << "Rank " << rank << ": " << i_points << ',' << i_start << '\n';
 
-    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);
-        }
+    // Simulate and save data to array
+    for (size_t i = 0; i < i_points; i++) {
+        i_data[i] = montecarlo::mcmc_parallel(L, i_start + dt * i, cycles,
+                                              burn_in_time);
+    }
 
-        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);
-                }
+    // Rank 0 collects all the data and copies it to the "master"
+    // data array.
+    if (rank == 0) {
+        // Copy its own i_data to the data array
+        std::copy_n(i_data, i_points, data);
+
+        // Collect i_data from other ranks in order and copy to 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);
             }
-            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';
+            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);
             }
-            ofile.close();
         }
-        else {
-            MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
-                     MPI_COMM_WORLD);
+
+        // Write everything from data to file
+        std::stringstream outfile;
+        outfile << "../output/phase_transition/mpi/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();
+    }
+    // For all other ranks, send the data to rank 0
+    else {
+        MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
+                 MPI_COMM_WORLD);
     }
 
     t1 = MPI_Wtime();
diff --git a/src/test_suite.cpp b/src/test_suite.cpp
index 6305451..a9cd8c5 100644
--- a/src/test_suite.cpp
+++ b/src/test_suite.cpp
@@ -12,8 +12,6 @@
 #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))
@@ -29,7 +27,7 @@
  * */
 class IsingModelTest {
 public:
-    /** @brief Test That initializing works as intended.
+    /** @brief Test that initializing works as intended.
      * */
     void test_init_functions()
     {
@@ -84,15 +82,9 @@ public:
         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)
@@ -102,44 +94,6 @@ public:
                 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;
     }
 };
@@ -150,18 +104,23 @@ int main()
     IsingModelTest test;
 
     test.test_init_functions();
+
     int res = 0;
     int tmp;
-    for (size_t i=0; i < 1000; i++) {
+    int iterations = 10000;
+    int accepted_values = 0;
+
+    // Run through the test multiple times to get a better estimate.
+    for (size_t i=0; i < iterations; i++) {
         tmp = test.test_2x2_lattice(1e-2, 1e5);
         if (tmp == 0) {
-            std::cout << "not enough cycles\n";
-            break;
+            continue;
         }
+        accepted_values++;
         res += tmp;
     }
     
-    std::cout << "Res: " << res / 1000 << std::endl;
+    std::cout << "Res: " << res / accepted_values << std::endl;
 
     return 0;
 }
diff --git a/src/time.cpp b/src/time.cpp
new file mode 100644
index 0000000..00f3327
--- /dev/null
+++ b/src/time.cpp
@@ -0,0 +1,68 @@
+/** @file time.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 0.1
+ *
+ *  @brief Timing various things
+ *
+ *  @bug No known bugs
+ * */
+#include "data_type.hpp"
+#include "monte_carlo.hpp"
+#include "utils.hpp"
+
+#include <omp.h>
+#include <ostream>
+
+void time_lattice_sizes()
+{
+    std::string outfile = "output/timing/lattice_sizes.txt";
+    std::ofstream ofile;
+
+    int lattice_sizes[] = {20, 40, 60, 80, 100};
+
+    utils::mkpath(utils::dirname(outfile));
+    ofile.open(outfile);
+    double t0, t1;
+    for (int L : lattice_sizes) {
+        t0 = omp_get_wtime();
+        montecarlo::phase_transition(L, 2.1, 2.4, 40, 100000,
+                                     montecarlo::mcmc_parallel,
+                                     "output/garbage/null.txt");
+        t1 = omp_get_wtime();
+        ofile << utils::scientific_format(L) << ','
+              << utils::scientific_format(t1 - t0) << '\n';
+    }
+    ofile.close();
+}
+
+void time_sample_sizes()
+{
+    std::string outfile = "output/timing/sample_sizes.txt";
+    std::ofstream ofile;
+
+    int sample_sizes[] = {1000, 10000, 100000};
+
+    utils::mkpath(utils::dirname(outfile));
+    ofile.open(outfile);
+    double t0, t1;
+    for (int samples : sample_sizes) {
+        t0 = omp_get_wtime();
+        montecarlo::phase_transition(20, 2.1, 2.4, 40, samples,
+                                     montecarlo::mcmc_parallel,
+                                     "output/garbage/null.txt");
+        t1 = omp_get_wtime();
+        ofile << utils::scientific_format(samples) << ','
+              << utils::scientific_format(t1 - t0) << '\n';
+    }
+    ofile.close();
+}
+
+int main()
+{
+    time_lattice_sizes();
+    time_sample_sizes();
+    return 0;
+}