diff --git a/pyfans/micro.cpp b/pyfans/micro.cpp
index fe9869a..f6b8981 100644
--- a/pyfans/micro.cpp
+++ b/pyfans/micro.cpp
@@ -30,6 +30,8 @@ py::array_t<double> merge_arrays(py::array_t<double> array1, py::array_t<double>
 // Constructor
 MicroSimulation::MicroSimulation(int sim_id)
 {
+    MPI_Init(NULL, NULL);
+
     // initialize fftw mpi
     fftw_mpi_init();
 
@@ -68,75 +70,18 @@ py::dict MicroSimulation::solve(py::dict macro_data, double dt)
 
     vector<double> g0(matmodel->n_str);
 
+    VectorXd homogenized_stress;
+
     for (int i_load = 0; i_load < n_loads; i_load++) {
         for (int i = 0; i < matmodel->n_str; ++i) {
             g0[i] = g0_all[i_load * matmodel->n_str + i];
         }
         matmodel->setGradient(g0);
         solver->solve();
-        solver->postprocess(reader, "result.h5", 0, 0);
         homogenized_stress = solver->get_homogenized_stress();
     }
 
-    Vector3d e1(1, 0, 0);
-    Vector3d e2(0, 1, 0);
-    Vector3d e3(0, 0, 1);
-
-    Matrix<double, 3, 3> delta_strains;
-    Matrix<double, 6, 6> perturbed_strains;
-
-    for (int i = 0; i < matmodel->n_str; i++) {
-        delta_strains.setZero();
-
-        if (i == 0) { // 11
-            delta_strains = (pert_param / 2.0) * (e1 * e1.transpose() + e1 * e1.transpose());
-        } else if (i == 1) { // 22
-            delta_strains = (pert_param / 2.0) * (e2 * e2.transpose() + e2 * e2.transpose());
-        } else if (i == 2) { // 33
-            delta_strains = (pert_param / 2.0) * (e3 * e3.transpose() + e3 * e3.transpose());
-        } else if (i == 3) { // 12
-            delta_strains = (pert_param / 2.0) * (e1 * e2.transpose() + e2 * e1.transpose());
-        } else if (i == 4) { // 13
-            delta_strains = (pert_param / 2.0) * (e1 * e3.transpose() + e3 * e1.transpose());
-        } else if (i == 5) { // 23
-            delta_strains = (pert_param / 2.0) * (e2 * e3.transpose() + e3 * e2.transpose());
-        }
-
-        // Construct perturbed strain matrix according to Mandel notation
-        perturbed_strains(i, 0) = g0[0] + delta_strains(0, 0);
-        perturbed_strains(i, 1) = g0[1] + delta_strains(1, 1);
-        perturbed_strains(i, 2) = g0[2] + delta_strains(2, 2);
-        perturbed_strains(i, 3) = g0[3] + sqrt(2) * delta_strains(0, 1);
-        perturbed_strains(i, 4) = g0[4] + sqrt(2) * delta_strains(0, 2);
-        perturbed_strains(i, 5) = g0[5] + sqrt(2) * delta_strains(1, 2);
-    }
-
-    std::cout << "Perturbed strains: " << perturbed_strains << std::endl;
-
-    vector<double> pert_strain;
-
-    // Calculate the homogenized stiffness matrix C using finite differences
-    for (int i = 0; i < matmodel->n_str; i++) {
-
-        for (int j = 0; j < matmodel->n_str; j++) {
-            pert_strain.push_back(perturbed_strains(i, j));
-        }
-
-        for (int j = 0; j < matmodel->n_str; j++) {
-            std::cout << "Perturbed strain[" << j << "] = " << pert_strain[j] << std::endl;
-        }
-
-        matmodel->setGradient(pert_strain);
-        solver->solve();
-        solver->postprocess(reader, "result.h5", 0, 0);
-        unperturbed_stress = solver->get_homogenized_stress();
-
-        for (int j = 0; j < matmodel->n_str; j++) {
-            C(i, j) = (unperturbed_stress[j] - homogenized_stress.data()[j]) / pert_param;
-        }
-
-        pert_strain.clear();
-    }
+    auto C = solver->get_homogenized_tangent(pert_param);
 
     std::cout << "Homogenized stiffness matrix C: " << C << std::endl;
 
diff --git a/pyfans/micro.hpp b/pyfans/micro.hpp
index 4432533..d2ffb19 100644
--- a/pyfans/micro.hpp
+++ b/pyfans/micro.hpp
@@ -33,9 +33,5 @@ class MicroSimulation {
     Reader              reader;
     Matmodel<3>        *matmodel;
     Solver<3>          *solver;
-
-    double               pert_param = 1e-3; // scalar strain perturbation parameter
-    std::vector<double>  homogenized_stress;
-    std::vector<double>  unperturbed_stress;
-    Matrix<double, 6, 6> C; // Homogenized stiffness matrix C
+    double              pert_param = 1e-3; // scalar strain perturbation parameter
 };
diff --git a/test/micro-manager-config-mech.json b/test/test_pyfans/micro-manager-config-mech.json
similarity index 100%
rename from test/micro-manager-config-mech.json
rename to test/test_pyfans/micro-manager-config-mech.json
diff --git a/test/micro_mech.py b/test/test_pyfans/micro_mech.py
similarity index 87%
rename from test/micro_mech.py
rename to test/test_pyfans/micro_mech.py
index 3ffc66c..31f6c5b 100644
--- a/test/micro_mech.py
+++ b/test/test_pyfans/micro_mech.py
@@ -1,6 +1,8 @@
-import PyFANS as fans  # alternatively PyFANSTHERMAL
+import PyFANS as fans
 import numpy as np
 
+# from mpi4py import MPI
+
 """
 Run FANS from python. For mechanical simulations, import PyFANS. For thermal simulations, use PyFANSTHERMAL.
 To be able to run the code, the python bindings must be compiled with the correct flags to initialize MPI.
@@ -8,6 +10,8 @@
 cmake -DUSE_MPI=ON ..
 """
 
+# MPI.Init()
+
 micro = fans.MicroSimulation(1)
 
 macro_data = dict()
@@ -20,3 +24,5 @@
 output = micro.solve(macro_data, dt)  # solve FANS for one load vector g0
 
 print(output)
+
+# MPI.Finalize()