diff --git a/CHANGELOG.md b/CHANGELOG.md
index 64ea413..4f0d076 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,7 @@
 
 ## latest
 
+- Add consistent homogenized tangent from the materials' own tangents (NEML2 and native models) [#163](https://github.com/DataAnalyticsEngineering/FANS/pull/163)
 - Add NEML2 material support [#162](https://github.com/DataAnalyticsEngineering/FANS/pull/162)
 - Account for grain orientations in GBDiffusion material model [#161](https://github.com/DataAnalyticsEngineering/FANS/pull/161)
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5b41a00..061fd3e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -156,7 +156,8 @@ target_include_directories(FANS_main PRIVATE "${PROJECT_BINARY_DIR}/include")
 # SOURCES
 # ##############################################################################
 
-target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp)
+target_sources(FANS_FANS PRIVATE src/reader.cpp src/plugin_material.cpp
+                                 src/tangent.cpp)
 target_sources(FANS_main PRIVATE src/main.cpp)
 
 # ##############################################################################
diff --git a/docs/NEML2.md b/docs/NEML2.md
index 90e5974..b4169b1 100644
--- a/docs/NEML2.md
+++ b/docs/NEML2.md
@@ -27,10 +27,10 @@ This gives `test/FANS_neml2`. Without pixi, configure with `-DFANS_NEML2=ON` in
 ## Compile a model
 
 ```bash
-pixi run -e dev-neml2 neml2-compile model.i "cpu cuda"
+pixi run -e dev-neml2 neml2-compile model.i --model model --device cpu cuda --dtype float64 --output-dir compiled_models -d state/S:forces/E
 ```
 
-This compiles the model block named like the file (`[model]`) for both devices into `compiled_models/model`, next to `model.i`. A third argument lists Python files the model needs, e.g. `"a.py b.py"`.
+This compiles the block `[model]` of `model.i` for both devices into `compiled_models/model`. `--load a.py` (repeatable) imports a Python file the model needs. `-d flux:gradient` also compiles the model's tangent, which the `homogenized_tangent` result needs.
 
 ## Input file
 
diff --git a/include/MaterialManager.h b/include/MaterialManager.h
index 0c4ecfe..a927a83 100644
--- a/include/MaterialManager.h
+++ b/include/MaterialManager.h
@@ -138,6 +138,13 @@ class MaterialManager {
         }
     }
 
+    // One linear model for all phases, e.g. the frozen tangent of src/tangent.cpp
+    MaterialManager(Matmodel<howmany, n_str> *model, int n_phases)
+        : phase_to_info(new MaterialInfo<howmany, n_str>[n_phases]), n_phases(n_phases), models{model}, all_stiffness(false)
+    {
+        std::fill_n(phase_to_info, n_phases, MaterialInfo<howmany, n_str>{model, nullptr, 0, true});
+    }
+
     ~MaterialManager()
     {
         for (auto *m : models)
@@ -221,10 +228,6 @@ class MaterialManager {
     {
         return phase_to_info[phase_id];
     }
-    inline size_t get_num_models() const
-    {
-        return models.size();
-    } // Get number of unique material models
     inline int get_num_phases() const
     {
         return n_phases;
diff --git a/include/fans_plugin.h b/include/fans_plugin.h
index 54a344e..5463f49 100644
--- a/include/fans_plugin.h
+++ b/include/fans_plugin.h
@@ -53,10 +53,11 @@ FANSPluginModel *fans_plugin_load(const char *config, char *msg, size_t msglen);
 int fans_plugin_set_table(FANSPluginModel *model, size_t field, int per_phase, size_t n_rows, const double *table,
                           char *err, size_t errlen);
 
-/* The step runs from time t_old to t. */
+/* The step runs from time t_old to t. Unless NULL, `tangent` gets the consistent
+   tangent d flux / d gradient, [n_points][n_str][n_str]. */
 int fans_plugin_evaluate(FANSPluginModel *model, size_t n_points, double t_old, double t, const double *gradient,
                          const int *voxel, const int *phase, const double *history_old, double *flux,
-                         double *history_new, char *err, size_t errlen);
+                         double *history_new, double *tangent, char *err, size_t errlen);
 
 void fans_plugin_free(FANSPluginModel *model);
 
diff --git a/include/material_models/large_strain/CompressibleNeoHookean.h b/include/material_models/large_strain/CompressibleNeoHookean.h
index 824c3b2..f65a9f5 100644
--- a/include/material_models/large_strain/CompressibleNeoHookean.h
+++ b/include/material_models/large_strain/CompressibleNeoHookean.h
@@ -47,6 +47,26 @@ class CompressibleNeoHookean : public LargeStrainMechModel {
         return lambda[mat_index] * logJ * C_inv + mu[mat_index] * (Matrix3d::Identity() - C_inv);
     }
 
+    bool has_tangent() const override
+    {
+        return true;
+    }
+
+    // A = dP/dF of P = mu F + (lambda log(J) - mu) G with G = F^{-T}:
+    //   A_iJkL = mu d_ik d_JL + lambda G_iJ G_kL + (mu - lambda log(J)) G_iL G_kJ
+    void get_tangent(int i_eps, int mat_index, ptrdiff_t element_idx, Tangent A) override
+    {
+        const Matrix3d F = extract_F(i_eps);
+        const Matrix3d G = F.inverse().transpose();
+        const double   c = mu[mat_index] - lambda[mat_index] * log(F.determinant());
+
+        for (int i = 0; i < 3; ++i)
+            for (int J = 0; J < 3; ++J)
+                for (int k = 0; k < 3; ++k)
+                    for (int L = 0; L < 3; ++L)
+                        A(3 * i + J, 3 * k + L) = mu[mat_index] * (i == k) * (J == L) + lambda[mat_index] * G(i, J) * G(k, L) + c * G(i, L) * G(k, J);
+    }
+
     // At F = I the tangent is that of linear isotropic elasticity
     Matrix<double, 9, 9> get_reference_stiffness() override
     {
diff --git a/include/material_models/thermal/GBDiffusion.h b/include/material_models/thermal/GBDiffusion.h
index 2021e84..8f0b824 100644
--- a/include/material_models/thermal/GBDiffusion.h
+++ b/include/material_models/thermal/GBDiffusion.h
@@ -74,8 +74,8 @@ class GBDiffusion : public ThermalModel, public LinearModel<1, 3> {
             const int num_GB       = static_cast<int>(boundary_count);
             n_mat                  = num_crystals + num_GB;
 
-            const vector<double> grain_rot_matrices = Reader::ReadData<double>(reader.ms_filename, reader.MSGroup() + "rotation_matrices");
-            const vector<double> GB_normals         = Reader::ReadData<double>(reader.ms_filename, reader.MSGroup() + "GB_normals");
+            const vector<double> grain_rot_matrices = Reader::ReadData<double>(reader.ms_filename, reader.ms_group + "rotation_matrices");
+            const vector<double> GB_normals         = Reader::ReadData<double>(reader.ms_filename, reader.ms_group + "GB_normals");
             if (grain_rot_matrices.size() < 9 * size_t(num_crystals) || GB_normals.size() < 3 * size_t(n_mat))
                 throw std::runtime_error("rotation_matrices or GB_normals has too few entries");
 
diff --git a/include/matmodel.h b/include/matmodel.h
index ddd7e5c..83f0cf4 100644
--- a/include/matmodel.h
+++ b/include/matmodel.h
@@ -12,13 +12,9 @@ class Matmodel {
   public:
     static constexpr int num_str = n_str; // length of strain and stress
 
-    int    verbosity; //!< output verbosity
-    int    n_mat;     //!< Number of Materials
-    int    n_gp;      //!< Number of Gauss points (computed from FE_type)
-    string FE_type;   //!< Finite element type: "HEX8", "HEX8R", "BBAR"
-
-    double *strain; //!< Gradient
-    double *stress; //!< Flux
+    int    n_mat;   //!< Number of Materials
+    int    n_gp;    //!< Number of Gauss points (computed from FE_type)
+    string FE_type; //!< Finite element type: "HEX8", "HEX8R", "BBAR"
 
     Matmodel(const Reader &reader);
 
@@ -27,6 +23,22 @@ class Matmodel {
     void                                     getStrainStress(double *strain, double *stress, Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx);
     void                                     setGradient(vector<double> _g0);
 
+    // As getStrainStress, for the tangents C = d sigma / d eps of the converged step at the
+    // element's Gauss points, [n_gp][n_str][n_str]
+    using Tangent = Map<Matrix<double, n_str, n_str, RowMajor>>;
+    void getTangent(double *tangent, Matrix<double, howmany * 8, 1> &ue, int mat_index, ptrdiff_t element_idx)
+    {
+        eps.noalias() = B * ue + g0;
+        for (int i = 0; i < n_gp; ++i)
+            get_tangent(n_str * i, mat_index, element_idx, Tangent(tangent + i * n_str * n_str));
+    }
+    // A nonlinear model with a tangent: a batched one gives it in evaluate_batch,
+    // any other writes a get_tangent and says so here
+    virtual bool has_tangent() const
+    {
+        return wants_batch();
+    }
+
     // Accessors for internal Gauss point data (populated after getStrainStress call)
     inline const double *get_eps_data() const
     {
@@ -36,10 +48,6 @@ class Matmodel {
     {
         return sigma.data();
     }
-    inline int get_n_gp() const
-    {
-        return n_gp;
-    }
 
     virtual void postprocess(Solver<howmany, n_str> &solver, Reader &reader, int load_idx, int time_idx) {};
 
@@ -52,8 +60,10 @@ class Matmodel {
     {
         return false;
     }
+    // Unless null, tangent_gp gets the tangent d flux / d gradient of the converged
+    // step, [element][n_gp][n_str][n_str] (src/tangent.cpp).
     virtual void evaluate_batch(const vector<ptrdiff_t> &elems, const unsigned short *phase,
-                                const double *ue, double *sig_gp) {}
+                                const double *ue, double *sig_gp, double *tangent_gp = nullptr) {}
 
     // A model without the element stiffnesses of a LinearModel can still say that its
     // flux is linear in the gradient: then the linear CG applies.
@@ -87,6 +97,16 @@ class Matmodel {
     void                                       Construct_B();
 
     virtual void get_sigma(int i, int mat_index, ptrdiff_t element_idx) = 0;
+    // As get_sigma, for the tangent C = d sigma / d eps at that Gauss point. This one is that
+    // of a linear model, its stresses at unit strains; a nonlinear model writes its own.
+    virtual void get_tangent(int i, int mat_index, ptrdiff_t element_idx, Tangent C)
+    {
+        for (int j = 0; j < n_str; ++j) {
+            eps.segment<n_str>(i) = Matrix<double, n_str, 1>::Unit(j);
+            get_sigma(i, mat_index, element_idx);
+            C.col(j) = sigma.segment<n_str>(i);
+        }
+    }
 };
 
 template <int howmany, int n_str>
diff --git a/include/reader.h b/include/reader.h
index 2295ad5..e2db91b 100644
--- a/include/reader.h
+++ b/include/reader.h
@@ -25,9 +25,9 @@ class Reader {
     ~Reader();
 
     // contents of input file:
-    char             ms_filename[4096]{};    // Name of Micro-structure hdf5 file
-    char             ms_datasetname[4096]{}; // Absolute path of Micro-structure in hdf5 file
-    char             results_prefix[4096]{};
+    string           ms_filename;    // Name of Micro-structure hdf5 file
+    string           ms_datasetname; // Absolute path of Micro-structure in hdf5 file
+    string           ms_group;       // Its group, where the datasets that go with it are: "/sphere/32x32x32/" for "/sphere/32x32x32/ms"
     int              n_mat;
     json             inputJson; // Complete input JSON (for MaterialManager)
     json             materialProperties;
@@ -45,9 +45,8 @@ class Reader {
     string           strain_type{"small"}; // "small" (default) or "large"
     string           FE_type;              // "HEX8" (default), "HEX8R", or "BBAR"
     vector<string>   resultsToWrite;
-    char             results_filename[4096]{}; // Output HDF5 filename
-    char             dataset_name[8192]{};     // Base path for results in HDF5 file
-    hid_t            results_file_id = -1;     // Open HDF5 file handle for results
+    string           dataset_name;         // Base path for results in HDF5 file
+    hid_t            results_file_id = -1; // Open HDF5 file handle for results
 
     // contents of microstructure file:
     vector<int>     dims;
@@ -66,12 +65,10 @@ class Reader {
     ptrdiff_t local_n1;
     ptrdiff_t local_1_start;
 
-    // void Setup(ptrdiff_t howmany);
     void ReadInputFile(const std::string &input_fn);
     void ReadMS(int hm);
     void FreeMS();
     void ComputeVolumeFractions();
-    // void ReadHDF5(char file_name[], char dset_name[]);
     void safe_create_group(hid_t file, const char *const name);
     void OpenResultsFile(const char *output_fn); // Open results file once
     void CloseResultsFile();                     // Explicitly close results file
@@ -94,10 +91,6 @@ class Reader {
     static std::vector<hsize_t> DataShape(const string &file, const string &dset_name);
     template <typename T>
     void ReadSlab(T *data, const std::vector<int> &extra_dims, const string &file, const string &dset_name) const;
-
-    // The microstructure's group, e.g. "/sphere/32x32x32/" for "/sphere/32x32x32/ms",
-    // where the datasets that go with it are
-    string MSGroup() const;
 };
 
 // Opens a dataset for reading; the caller closes it and file_id
@@ -262,9 +255,8 @@ void Reader::writeData(const char *fieldName, int load_idx, int time_idx, const
     if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
         return;
     }
-    char name[5096];
-    snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
-    WriteData(data, name, shape, rank);
+    const string name = dataset_name + "/load" + to_string(load_idx) + "/time_step" + to_string(time_idx) + "/" + fieldName;
+    WriteData(data, name.c_str(), shape, rank);
 }
 
 template <typename T>
@@ -273,9 +265,8 @@ void Reader::writeSlab(const char *fieldName, int load_idx, int time_idx, const
     if (std::find(resultsToWrite.begin(), resultsToWrite.end(), fieldName) == resultsToWrite.end()) {
         return;
     }
-    char name[5096];
-    snprintf(name, sizeof(name), "%s/load%i/time_step%i/%s", dataset_name, load_idx, time_idx, fieldName);
-    WriteSlab(data, extra_dims, name);
+    const string name = dataset_name + "/load" + to_string(load_idx) + "/time_step" + to_string(time_idx) + "/" + fieldName;
+    WriteSlab(data, extra_dims, name.c_str());
 }
 
 #endif
diff --git a/include/solver.h b/include/solver.h
index b6f8789..df2fe36 100644
--- a/include/solver.h
+++ b/include/solver.h
@@ -4,8 +4,6 @@
 #include "matmodel.h"
 #include "MaterialManager.h"
 
-class J2Plasticity;
-
 typedef Map<Array<double, Dynamic, Dynamic>, Unaligned, OuterStride<>> RealArray;
 
 template <int howmany, int n_str>
@@ -51,6 +49,7 @@ class Solver : private MixedBCController<howmany> {
 
     template <int padding, typename F>
     void iterateCubes(F f);
+    void update_ghost_layer(double *u); //!< the next rank's first layer of nodal values, behind this rank's own
 
     void         solve();
     void         extrapolateDisplacement(); //!< Linear extrapolation for next time step
@@ -73,7 +72,7 @@ class Solver : private MixedBCController<howmany> {
     void     evaluate_batched_stress(double *u); //!< in src/plugin_material.cpp
 
     MatrixXd homogenized_tangent;
-    MatrixXd get_homogenized_tangent(double pert_param);
+    MatrixXd get_homogenized_tangent(double pert_param); //!< in src/tangent.cpp
 
     void enableMixedBC(const MixedBC &mbc, size_t step)
     {
@@ -229,6 +228,15 @@ void Solver<howmany, n_str>::CreateFFTWPlans(double *in, fftw_complex *transform
     new (&rhat) Map<VectorXcd>((std::complex<double> *) transformed, local_n1 * n_x * (n_z / 2 + 1) * howmany);
 }
 
+// The elements of a rank's last layer reach into the first layer of nodes of the next rank
+template <int howmany, int n_str>
+void Solver<howmany, n_str>::update_ghost_layer(double *u)
+{
+    const int layer = n_y * n_z * howmany;
+    MPI_Sendrecv(u, layer, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0,
+                 u + local_n0 * layer, layer, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE);
+}
+
 // TODO: possibly circumvent the padding problem by accessing r as a matrix?
 template <int howmany, int n_str>
 template <int padding, typename F>
@@ -243,10 +251,7 @@ void Solver<howmany, n_str>::compute_residual_basic(RealArray &r_matrix, RealArr
         r[i] = 0;
     }
 
-    // int MPI_Sendrecv(void *sendbuf, int sendcount, MPI_Datatype sendtype, int dest, int sendtag, void *recvbuf,
-    //           int recvcount, MPI_Datatype recvtype, int source, int recvtag, MPI_Comm comm, MPI_Status *status)
-    MPI_Sendrecv(u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0,
-                 u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE);
+    update_ghost_layer(u);
 
     Matrix<double, howmany * 8, 1> ue;
 
@@ -481,11 +486,11 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
     // Determine if we need to compute stress/strain at all
     bool need_compute = need_stress || need_stress_gp || need_strain || need_strain_gp || need_global_avg || need_phase_avg;
 
-    // Conditional allocation
-    VectorXd *strain_elem = need_strain ? new VectorXd(local_n0 * n_y * n_z * n_str) : nullptr;
-    VectorXd *stress_elem = need_stress ? new VectorXd(local_n0 * n_y * n_z * n_str) : nullptr;
-    VectorXd *strain_gp   = need_strain_gp ? new VectorXd(local_n0 * n_y * n_z * n_gp * n_str) : nullptr;
-    VectorXd *stress_gp   = need_stress_gp ? new VectorXd(local_n0 * n_y * n_z * n_gp * n_str) : nullptr;
+    // The requested fields; the others stay empty
+    VectorXd strain_elem(need_strain ? local_n0 * n_y * n_z * n_str : 0);
+    VectorXd stress_elem(need_stress ? local_n0 * n_y * n_z * n_str : 0);
+    VectorXd strain_gp(need_strain_gp ? local_n0 * n_y * n_z * n_gp * n_str : 0);
+    VectorXd stress_gp(need_stress_gp ? local_n0 * n_y * n_z * n_gp * n_str : 0);
 
     VectorXd stress_average = VectorXd::Zero(n_str);
     VectorXd strain_average = VectorXd::Zero(n_str);
@@ -496,8 +501,7 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
     vector<VectorXd> phase_strain_average(n_mat, VectorXd::Zero(n_str));
     vector<int>      phase_counts(n_mat, 0);
 
-    MPI_Sendrecv(v_u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0,
-                 v_u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE);
+    update_ghost_layer(v_u);
 
     Matrix<double, howmany * 8, 1> ue;
     int                            phase_id;
@@ -523,23 +527,23 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
             // Store element averages if requested
             if (need_stress) {
                 for (int c = 0; c < n_str; ++c) {
-                    (*stress_elem)(idx[0] * n_str + c) = elem_stress_avg[c];
+                    stress_elem(idx[0] * n_str + c) = elem_stress_avg[c];
                 }
             }
             if (need_strain) {
                 for (int c = 0; c < n_str; ++c) {
-                    (*strain_elem)(idx[0] * n_str + c) = elem_strain_avg[c];
+                    strain_elem(idx[0] * n_str + c) = elem_strain_avg[c];
                 }
             }
 
             // Copy all GP data if requested
             if (need_stress_gp) {
-                memcpy(&(*stress_gp)(idx[0] * n_gp * n_str),
+                memcpy(&stress_gp(idx[0] * n_gp * n_str),
                        info.model->get_sigma_data(),
                        n_gp * n_str * sizeof(double));
             }
             if (need_strain_gp) {
-                memcpy(&(*strain_gp)(idx[0] * n_gp * n_str),
+                memcpy(&strain_gp(idx[0] * n_gp * n_str),
                        info.model->get_eps_data(),
                        n_gp * n_str * sizeof(double));
             }
@@ -585,65 +589,32 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
     homogenized_stress = stress_average;
     homogenized_strain = strain_average;
 
-    /* ====================================================================== *
-     *  u_total = g0·X  +  ũ          (vector or scalar, decided at compile time)
-     * ====================================================================== */
-    const vector<double> &ml  = matmanager->models[0]->macroscale_loading;
-    const double          dx  = reader.l_e[0];
-    const double          dy  = reader.l_e[1];
-    const double          dz  = reader.l_e[2];
-    const double          Lx2 = reader.L[0] / 2.0;
-    const double          Ly2 = reader.L[1] / 2.0;
-    const double          Lz2 = reader.L[2] / 2.0;
-    constexpr double      rs2 = 0.7071067811865475; // 1.0 / std::sqrt(2.0)
-    VectorXd              u_total(local_n0 * n_y * n_z * howmany);
-    /* ---------- single sweep ------------------------------------------------- */
+    // u_total = u + G X, with G the macroscale gradient of the nodal values: the temperature
+    // gradient, the strain (ml holds its Mandel components) or F - I
+    const vector<double>      &ml  = matmanager->models[0]->macroscale_loading;
+    constexpr double           rs2 = 0.7071067811865475; // 1.0 / std::sqrt(2.0)
+    Matrix<double, howmany, 3> G;
+    if constexpr (n_str == 3)
+        G << ml[0], ml[1], ml[2];
+    else if constexpr (n_str == 6)
+        G << ml[0], ml[3] * rs2, ml[4] * rs2,
+            ml[3] * rs2, ml[1], ml[5] * rs2,
+            ml[4] * rs2, ml[5] * rs2, ml[2];
+    else
+        G << ml[0] - 1.0, ml[1], ml[2],
+            ml[3], ml[4] - 1.0, ml[5],
+            ml[6], ml[7], ml[8] - 1.0;
+
+    VectorXd  u_total(local_n0 * n_y * n_z * howmany);
     ptrdiff_t n = 0;
     for (ptrdiff_t ix = 0; ix < local_n0; ++ix) {
-        const double x = (local_0_start + ix) * dx - Lx2;
+        const double x = (local_0_start + ix) * reader.l_e[0] - reader.L[0] / 2.0;
         for (ptrdiff_t iy = 0; iy < n_y; ++iy) {
-            const double y = iy * dy - Ly2;
+            const double y = iy * reader.l_e[1] - reader.L[1] / 2.0;
             for (ptrdiff_t iz = 0; iz < n_z; ++iz, ++n) {
-                const double z = iz * dz - Lz2;
-                if (howmany == 3) { /* ===== mechanics (vector) ===== */
-                    if constexpr (n_str == 6) {
-                        /* Small strain: ml holds 6 Mandel components of the mean strain */
-                        const double    g11 = ml[0];
-                        const double    g22 = ml[1];
-                        const double    g33 = ml[2];
-                        const double    g12 = ml[3] * rs2;
-                        const double    g13 = ml[4] * rs2;
-                        const double    g23 = ml[5] * rs2;
-                        const double    ux  = g11 * x + g12 * y + g13 * z;
-                        const double    uy  = g12 * x + g22 * y + g23 * z;
-                        const double    uz  = g13 * x + g23 * y + g33 * z;
-                        const ptrdiff_t b   = 3 * n;
-                        u_total[b]          = v_u[b] + ux;
-                        u_total[b + 1]      = v_u[b + 1] + uy;
-                        u_total[b + 2]      = v_u[b + 2] + uz;
-                    } else if constexpr (n_str == 9) {
-                        /* Large strain: ml holds 9 components of the mean deformation gradient F */
-                        const double F11 = ml[0];
-                        const double F12 = ml[1];
-                        const double F13 = ml[2];
-                        const double F21 = ml[3];
-                        const double F22 = ml[4];
-                        const double F23 = ml[5];
-                        const double F31 = ml[6];
-                        const double F32 = ml[7];
-                        const double F33 = ml[8];
-                        // u = (F - I) * X
-                        const double    ux = (F11 - 1.0) * x + F12 * y + F13 * z;
-                        const double    uy = F21 * x + (F22 - 1.0) * y + F23 * z;
-                        const double    uz = F31 * x + F32 * y + (F33 - 1.0) * z;
-                        const ptrdiff_t b  = 3 * n;
-                        u_total[b]         = v_u[b] + ux;
-                        u_total[b + 1]     = v_u[b + 1] + uy;
-                        u_total[b + 2]     = v_u[b + 2] + uz;
-                    }
-                } else { /* ===== scalar (howmany==1) ==== */
-                    u_total[n] = v_u[n] + (ml[0] * x + ml[1] * y + ml[2] * z);
-                }
+                const double z = iz * reader.l_e[2] - reader.L[2] / 2.0;
+                for (int i = 0; i < howmany; ++i)
+                    u_total[howmany * n + i] = v_u[howmany * n + i] + (G(i, 0) * x + G(i, 1) * y + G(i, 2) * z);
             }
         }
     }
@@ -670,26 +641,16 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
     reader.writeSlab("residual", load_idx, time_idx, v_r, {howmany});
 
     if (need_strain)
-        reader.writeSlab("strain", load_idx, time_idx, strain_elem->data(), {n_str});
+        reader.writeSlab("strain", load_idx, time_idx, strain_elem.data(), {n_str});
     if (need_stress)
-        reader.writeSlab("stress", load_idx, time_idx, stress_elem->data(), {n_str});
+        reader.writeSlab("stress", load_idx, time_idx, stress_elem.data(), {n_str});
     if (need_strain_gp)
-        reader.writeSlab("strain_gp", load_idx, time_idx, strain_gp->data(), {n_gp, n_str});
+        reader.writeSlab("strain_gp", load_idx, time_idx, strain_gp.data(), {n_gp, n_str});
     if (need_stress_gp)
-        reader.writeSlab("stress_gp", load_idx, time_idx, stress_gp->data(), {n_gp, n_str});
+        reader.writeSlab("stress_gp", load_idx, time_idx, stress_gp.data(), {n_gp, n_str});
 
     matmanager->postprocess(*this, reader, load_idx, time_idx);
 
-    // Cleanup
-    if (strain_elem)
-        delete strain_elem;
-    if (stress_elem)
-        delete stress_elem;
-    if (strain_gp)
-        delete strain_gp;
-    if (stress_gp)
-        delete stress_gp;
-
     // Compute homogenized tangent only if requested
     if (find(reader.resultsToWrite.begin(), reader.resultsToWrite.end(), "homogenized_tangent") != reader.resultsToWrite.end()) {
         homogenized_tangent = get_homogenized_tangent(1e-6);
@@ -707,17 +668,14 @@ void Solver<howmany, n_str>::postprocess(Reader &reader, int load_idx, int time_
 template <int howmany, int n_str>
 VectorXd Solver<howmany, n_str>::get_homogenized_stress()
 {
-
-    VectorXd strain    = VectorXd::Zero(local_n0 * n_y * n_z * n_str);
-    VectorXd stress    = VectorXd::Zero(local_n0 * n_y * n_z * n_str);
     homogenized_stress = VectorXd::Zero(n_str);
 
-    MPI_Sendrecv(v_u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0,
-                 v_u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE);
+    update_ghost_layer(v_u);
     if (matmanager->any_batched)
         evaluate_batched_stress(v_u);
 
     Matrix<double, howmany * 8, 1> ue;
+    Matrix<double, n_str, 1>       strain, stress; // of one element
     int                            phase_id;
     iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *idxPadding) {
         for (int i = 0; i < 8; ++i) {
@@ -728,8 +686,8 @@ VectorXd Solver<howmany, n_str>::get_homogenized_stress()
         phase_id = ms[idx[0]];
 
         const MaterialInfo<howmany, n_str> &info = matmanager->get_info(phase_id);
-        info.model->getStrainStress(strain.segment(n_str * idx[0], n_str).data(), stress.segment(n_str * idx[0], n_str).data(), ue, info.local_mat_id, idx[0]);
-        homogenized_stress += stress.segment(n_str * idx[0], n_str);
+        info.model->getStrainStress(strain.data(), stress.data(), ue, info.local_mat_id, idx[0]);
+        homogenized_stress += stress;
     });
 
     MPI_Allreduce(MPI_IN_PLACE, homogenized_stress.data(), n_str, MPI_DOUBLE, MPI_SUM, communicator);
@@ -738,44 +696,6 @@ VectorXd Solver<howmany, n_str>::get_homogenized_stress()
     return homogenized_stress;
 }
 
-template <int howmany, int n_str>
-MatrixXd Solver<howmany, n_str>::get_homogenized_tangent(double pert_param)
-{
-    homogenized_tangent               = MatrixXd::Zero(n_str, n_str);
-    VectorXd       unperturbed_stress = get_homogenized_stress();
-    VectorXd       perturbed_stress;
-    vector<double> pert_strain(n_str, 0.0);
-    vector<double> g0       = matmanager->get_info(0).model->macroscale_loading;
-    bool           islinear = matmanager->all_linear;
-
-    for (auto *model : matmanager->models) {
-        if (dynamic_cast<J2Plasticity *>(model) != nullptr) {
-            throw std::runtime_error("Homogenized tangent computation not implemented for J2Plasticity models.");
-        }
-    }
-    // TODO: a deep copy of the solver object is needed here to avoid modifying the history of the solver object
-
-    for (int i = 0; i < n_str; i++) {
-        if (islinear) {
-            pert_strain    = vector<double>(n_str, 0.0);
-            pert_strain[i] = 1.0;
-        } else {
-            pert_strain = g0;
-            pert_strain[i] += pert_param;
-        }
-
-        matmanager->set_gradient(pert_strain);
-        disableMixedBC();
-        solve();
-        perturbed_stress = get_homogenized_stress();
-
-        homogenized_tangent.col(i) = islinear ? perturbed_stress : (perturbed_stress - unperturbed_stress) / pert_param;
-    }
-
-    homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval();
-    return homogenized_tangent;
-}
-
 template <int howmany, int n_str>
 Solver<howmany, n_str>::~Solver()
 {
diff --git a/pixi.toml b/pixi.toml
index d5ca5d9..a41fd92 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -70,10 +70,6 @@ neml2 = ">=3.1.0, <4"
 [feature.neml2.tasks]
 build-fans_neml2 = "cmake -S . -B build-neml2 -DFANS_NEML2=ON -DFANS_ENABLE_NATIVE=ON && cmake --build build-neml2 -j"
 
-[feature.neml2.tasks.neml2-compile]
-args = ["file", { arg = "device", default = "cpu" }, { arg = "load", default = "" }]
-cmd  = "cd \"$INIT_CWD\" && neml2-compile {{ file }} --model $(basename {{ file }} .i) --device {{ device }} --dtype float64 --output-dir $(dirname {{ file }})/compiled_models{% for f in load | split %} --load {{ f }}{% endfor %}"
-
 [tasks]
 test = { cmd = "pytest -v -s --from-pixi", cwd = "test/pytest", depends-on = ["test-fans"] }
 test-fans = { cmd = "./run_tests.sh -n {{n}}", args = [
diff --git a/src/neml2_plugin.cpp b/src/neml2_plugin.cpp
index 5ef6482..bdefb40 100644
--- a/src/neml2_plugin.cpp
+++ b/src/neml2_plugin.cpp
@@ -58,6 +58,7 @@ struct FANSPluginModel {
     std::unique_ptr<neml2::aoti::Model> model;
     at::Device                          device{at::kCPU};
     Var                                 gradient, flux;
+    Var                                 tangent; // d flux / d gradient: the flux's shape, then the gradient's
     std::vector<Var>                    history, fields;
     int64_t                             n_history{0};                                     // history doubles per point
     at::Tensor                          perm = at::tensor({0, 1, 2, 5, 4, 3}, at::kLong); // FANS <-> NEML2 Mandel order
@@ -122,6 +123,10 @@ FANSPluginModel *fans_plugin_load(const char *config, char *msg, size_t msglen)
         if (h->flux.size != h->gradient.size)
             throw std::runtime_error("the gradient '" + gradient + "' and the flux '" + flux + "' differ in size");
 
+        std::vector<int64_t> tangent_shape = out_shape[fi];
+        tangent_shape.insert(tangent_shape.end(), in_shape[gi].begin(), in_shape[gi].end());
+        h->tangent = make_var("tangent", tangent_shape);
+
         for (size_t i = 0; i < in.size(); ++i) {
             const std::string &name = in[i];
             if (name == gradient || name == "t" || name == "t~1")
@@ -176,7 +181,7 @@ int fans_plugin_set_table(FANSPluginModel *m, size_t field, int per_phase, size_
 
 int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, double t, const double *gradient,
                          const int *voxel, const int *phase, const double *history_old, double *flux,
-                         double *history_new, char *err, size_t errlen)
+                         double *history_new, double *tangent, char *err, size_t errlen)
 {
     try {
         // A FANS array as an [n][size] tensor, without copying
@@ -199,6 +204,8 @@ int fans_plugin_evaluate(FANSPluginModel *m, size_t n_points, double t_old, doub
         }
 
         const auto outputs = m->model->forward(inputs);
+        if (tangent) // from the Jacobian of the same evaluation
+            m->to_fans(m->model->jacobian(inputs).second.at(m->flux.name).at(m->gradient.name), m->tangent, view(tangent, m->tangent.size));
 
         m->to_fans(outputs.at(m->flux.name), m->flux, view(flux, m->flux.size));
         for (const Var &h : m->history)
diff --git a/src/plugin_material.cpp b/src/plugin_material.cpp
index 18e49fc..6a0b211 100644
--- a/src/plugin_material.cpp
+++ b/src/plugin_material.cpp
@@ -9,7 +9,8 @@
 // History: the plugin's internal variables at every Gauss point: `converged`,
 //   the state a step starts from, and `trial`, the state of the step itself.
 //   Once a step has converged, the solver evaluates it a last time (for the
-//   output) and `trial` becomes `converged`.
+//   output) and `trial` becomes `converged`; `previous` keeps the state that
+//   step started from, for its tangent.
 // Fields: further model inputs, e.g. an orientation, read from the
 //   microstructure file per voxel or per phase (read_field).
 // Linear: "linear": true in the material properties says that the flux is
@@ -49,15 +50,13 @@ Plugin load_plugin(const string &name)
         throw std::runtime_error("Could not load the '" + name + "' material plugin: " + dlerror() +
                                  "\nBuild FANS with -DFANS_" + name + "=ON.");
 
-    auto sym = [&](auto &f, const char *n) {
-        if (!(f = reinterpret_cast<std::remove_reference_t<decltype(f)>>(dlsym(h, n))))
-            throw std::runtime_error(lib + " is missing '" + n + "'");
-    };
     Plugin p;
-    sym(p.load, "fans_plugin_load");
-    sym(p.set_table, "fans_plugin_set_table");
-    sym(p.evaluate, "fans_plugin_evaluate");
-    sym(p.free_model, "fans_plugin_free");
+    p.load       = reinterpret_cast<decltype(p.load)>(dlsym(h, "fans_plugin_load"));
+    p.set_table  = reinterpret_cast<decltype(p.set_table)>(dlsym(h, "fans_plugin_set_table"));
+    p.evaluate   = reinterpret_cast<decltype(p.evaluate)>(dlsym(h, "fans_plugin_evaluate"));
+    p.free_model = reinterpret_cast<decltype(p.free_model)>(dlsym(h, "fans_plugin_free"));
+    if (!p.load || !p.set_table || !p.evaluate || !p.free_model)
+        throw std::runtime_error(lib + " lacks a function of fans_plugin.h");
     return p;
 }
 
@@ -71,7 +70,7 @@ struct Field {
 // absolute path; per voxel if it is [Z][Y][X][...] on the grid, else per phase [n_phase][...]
 Field read_field(const Reader &reader, const string &name, const string &dset, int size)
 {
-    const string          path      = dset[0] == '/' ? dset : reader.MSGroup() + dset;
+    const string          path      = dset[0] == '/' ? dset : reader.ms_group + dset;
     const vector<hsize_t> shape     = Reader::DataShape(reader.ms_filename, path);
     const bool            per_voxel = shape.size() >= 3 && shape[0] == hsize_t(reader.dims[2]) && shape[1] == hsize_t(reader.dims[1]) &&
                                       shape[2] == hsize_t(reader.dims[0]);
@@ -156,7 +155,7 @@ class PluginModel : public Base {
     // The history of every Gauss point, and the arrays of one plugin call
     void initializeInternalVariables(ptrdiff_t num_elements, int num_gauss_points) override
     {
-        trial = converged = vector<double>(num_elements * num_gauss_points * n_history, 0.0);
+        previous = trial = converged = vector<double>(num_elements * num_gauss_points * n_history, 0.0);
         gradient_buf.resize(chunk * n_gp * n_str);
         flux_buf.resize(chunk * n_gp * n_str);
         history_old_buf.resize(chunk * n_gp * n_history);
@@ -168,6 +167,7 @@ class PluginModel : public Base {
     // A step has converged: its state is what the next step starts from
     void updateInternalVariables() override
     {
+        previous.swap(converged);
         converged = trial;
     }
 
@@ -179,11 +179,10 @@ class PluginModel : public Base {
         return {flux_cache + size_t(element_idx) * n_gp * n_str, n_gp * n_str};
     }
 
-    // Called per Gauss point: the first call brings the fluxes of them all
+    // The flux at a Gauss point, from the last evaluate_batch
     void get_sigma(int i, int, ptrdiff_t element_idx) override
     {
-        if (i == 0)
-            sigma = fluxes(element_idx);
+        sigma.segment(i, n_str) = fluxes(element_idx).segment(i, n_str);
     }
 
     // The element's residual straight from its fluxes: the gradient and the loop over
@@ -205,15 +204,19 @@ class PluginModel : public Base {
         return linear;
     }
 
-    // The flux at all Gauss points of `elems`, into flux_all, and their state, into `trial`
+    // The flux at all Gauss points of `elems`, into flux_all, and their state, into `trial`.
+    // With tangent_all, the converged step once more, from the state it started from, for its tangent
     void evaluate_batch(const vector<ptrdiff_t> &elems, const unsigned short *phase, const double *ue_all,
-                        double *flux_all) override
+                        double *flux_all, double *tangent_all) override
     {
         const size_t n_dof     = size_t(howmany) * 8; // nodal values per element
         const size_t per_elem  = size_t(n_gp) * n_str;
         const size_t hist_elem = size_t(n_gp) * n_history;
         flux_cache             = flux_all;
 
+        const vector<double> &old = tangent_all ? previous : converged;
+        vector<double>        tangent_buf(tangent_all ? chunk * per_elem * n_str : 0);
+
         for (size_t b = 0; b < elems.size(); b += chunk) {
             const size_t ne = std::min(chunk, elems.size() - b);
 
@@ -225,19 +228,22 @@ class PluginModel : public Base {
                 std::fill_n(&phase_buf[k * n_gp], n_gp, phase[e]);
                 if (phase[e] >= n_phase_rows)
                     throw std::runtime_error("Phase " + std::to_string(phase[e]) + " has no row in the per-phase fields of its plugin material.");
-                std::copy_n(converged.data() + e * hist_elem, hist_elem, history_old_buf.data() + k * hist_elem);
+                std::copy_n(old.data() + e * hist_elem, hist_elem, history_old_buf.data() + k * hist_elem);
                 std::copy_n(trial.data() + e * hist_elem, hist_elem, history_new_buf.data() + k * hist_elem); // initial guess
             }
 
             char err[FANS_PLUGIN_MSGLEN] = {0};
             if (plugin.evaluate(model, ne * n_gp, this->time_old, this->time, gradient_buf.data(), voxel_buf.data(), phase_buf.data(),
-                                history_old_buf.data(), flux_buf.data(), history_new_buf.data(), err, sizeof(err)) != 0)
+                                history_old_buf.data(), flux_buf.data(), history_new_buf.data(), tangent_all ? tangent_buf.data() : nullptr,
+                                err, sizeof(err)) != 0)
                 throw std::runtime_error(string("Plugin material evaluation failed: ") + err);
 
             for (size_t k = 0; k < ne; ++k) {
                 const ptrdiff_t e = elems[b + k];
                 std::copy_n(&flux_buf[k * per_elem], per_elem, flux_all + e * per_elem);
                 std::copy_n(history_new_buf.data() + k * hist_elem, hist_elem, trial.data() + e * hist_elem);
+                if (tangent_all)
+                    std::copy_n(&tangent_buf[k * per_elem * n_str], per_elem * n_str, tangent_all + e * per_elem * n_str);
             }
         }
     }
@@ -251,7 +257,7 @@ class PluginModel : public Base {
         if (!need && !need_gp)
             return;
         const size_t n_pts = size_t(solver.local_n0 * solver.n_y * solver.n_z) * n_gp;
-        const string dir   = string(reader.dataset_name) + "/load" + std::to_string(load_idx) + "/time_step" + std::to_string(time_idx) + "/";
+        const string dir   = reader.dataset_name + "/load" + std::to_string(load_idx) + "/time_step" + std::to_string(time_idx) + "/";
 
         size_t offset = 0;
         for (auto &[name, size] : history_vars) {
@@ -287,7 +293,7 @@ class PluginModel : public Base {
     vector<std::pair<string, int>> history_vars;                                             // name, doubles per point
     int                            n_history{0};                                             // history doubles per point
     size_t                         n_phase_rows{SIZE_MAX};                                   // phases below have a row in every per-phase field
-    vector<double>                 trial, converged;                                         // history at every Gauss point
+    vector<double>                 trial, converged, previous;                               // history at every Gauss point
     vector<double>                 gradient_buf, flux_buf, history_old_buf, history_new_buf; // one plugin call
     vector<int>                    voxel_buf, phase_buf;
     const double                  *flux_cache{nullptr};
@@ -330,8 +336,7 @@ void Solver<howmany, n_str>::evaluate_batched_stress(double *u)
         batch_gp_stress.resize(size_t(n_elem) * models[0]->n_gp * n_str);
     }
 
-    MPI_Sendrecv(u, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + world_size - 1) % world_size, 0,
-                 u + local_n0 * n_y * n_z * howmany, n_y * n_z * howmany, MPI_DOUBLE, (world_rank + 1) % world_size, 0, communicator, MPI_STATUS_IGNORE);
+    update_ghost_layer(u);
 
     iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *idxPadding) {
         if (matmanager->get_info(ms[idx[0]]).model->wants_batch())
diff --git a/src/reader.cpp b/src/reader.cpp
index 6be704c..1014865 100644
--- a/src/reader.cpp
+++ b/src/reader.cpp
@@ -1,15 +1,8 @@
 #include "general.h"
 #include "reader.h"
 
-#include "H5Cpp.h"
-#include "fftw3-mpi.h"
-#include "hdf5.h"
-#include "mpi.h"
 #include <cstdlib>
 
-#include "H5FDmpi.h"
-#include "H5FDmpio.h"
-
 Reader::Reader(const MPI_Comm &comm)
     : communicator(comm)
 {
@@ -19,47 +12,25 @@ Reader::Reader(const MPI_Comm &comm)
 
 void Reader::ComputeVolumeFractions()
 {
-    unsigned short local_max  = 0;
-    unsigned short local_min  = USHRT_MAX;
-    size_t         local_size = local_n0 * dims[1] * dims[2];
-
-    // Find the local maximum and minimum material indices
-    for (size_t i = 0; i < local_size; i++) {
-        unsigned short val = static_cast<unsigned short>(ms[i]);
-        if (val > local_max) {
-            local_max = val;
-        }
-        if (val < local_min) {
-            local_min = val;
-        }
-    }
-
-    // Find the global maximum and minimum material indices
-    unsigned short global_max, global_min;
-    MPI_Allreduce(&local_max, &global_max, 1, MPI_UNSIGNED_SHORT, MPI_MAX, communicator);
-    MPI_Allreduce(&local_min, &global_min, 1, MPI_UNSIGNED_SHORT, MPI_MIN, communicator);
-
-    if (global_min != 0)
+    const size_t local_size = local_n0 * dims[1] * dims[2];
+
+    // The phases are 0 ... n_mat - 1
+    int highest = *std::max_element(ms, ms + local_size);
+    MPI_Allreduce(MPI_IN_PLACE, &highest, 1, MPI_INT, MPI_MAX, communicator);
+    n_mat = highest + 1;
+
+    // Voxels of each phase, summed over all ranks
+    std::vector<long> voxels(n_mat, 0);
+    for (size_t i = 0; i < local_size; i++)
+        voxels[ms[i]]++;
+    MPI_Allreduce(MPI_IN_PLACE, voxels.data(), n_mat, MPI_LONG, MPI_SUM, communicator);
+    if (voxels[0] == 0)
         throw std::invalid_argument("Microstructure phase IDs must start at 0");
 
-    // Calculate total number of materials
-    n_mat = global_max - global_min + 1;
-
-    Log::logger().info("# Number of materials: {} (from {} to {})", n_mat, global_min, global_max);
+    Log::logger().info("# Number of materials: {} (from 0 to {})", n_mat, highest);
     Log::logger().info("# Volume fractions");
-
-    // Voxels of each phase, summed over all ranks at once
-    std::vector<long> vol_frac(n_mat, 0);
-    for (size_t i = 0; i < local_size; i++) {
-        unsigned short val   = static_cast<unsigned short>(ms[i]);
-        int            index = val - global_min; // Adjust index to start from 0
-        vol_frac[index]++;
-    }
-    MPI_Allreduce(MPI_IN_PLACE, vol_frac.data(), n_mat, MPI_LONG, MPI_SUM, communicator);
-
     for (int i = 0; i < n_mat; i++)
-        Log::logger().info("# material {:4}    vol. frac. {:10.4f}%  ", static_cast<unsigned int>(i) + global_min,
-                           100. * vol_frac[i] / (double(dims[0]) * dims[1] * dims[2]));
+        Log::logger().info("# material {:4}    vol. frac. {:10.4f}%  ", i, 100. * voxels[i] / (double(dims[0]) * dims[1] * dims[2]));
 }
 
 void Reader ::ReadInputFile(const std::string &input_fn)
@@ -71,23 +42,17 @@ void Reader ::ReadInputFile(const std::string &input_fn)
         inputJson = j; // Store complete input JSON for MaterialManager
 
         microstructure = j["microstructure"];
-        std::snprintf(ms_filename, sizeof(ms_filename), "%s", microstructure["filepath"].get<std::string>().c_str());
-        // dataset name handling
-        const auto tmp_str = microstructure["datasetname"].get<std::string>();
-        if (tmp_str.empty())
+        ms_filename    = microstructure["filepath"].get<string>();
+        ms_datasetname = microstructure["datasetname"].get<string>();
+        if (ms_datasetname.empty())
             throw std::invalid_argument("datasetname must not be empty and must refer to a valid HDF5 path");
-        // Ensure absolute HDF5 path, leading slash
-        std::snprintf(ms_datasetname, sizeof(ms_datasetname), "%s%s", tmp_str.front() == '/' ? "" : "/", tmp_str.c_str());
+        if (ms_datasetname.front() != '/')
+            ms_datasetname = "/" + ms_datasetname; // an absolute HDF5 path
         L = microstructure["L"].get<vector<double>>();
 
-        if (j.contains("results_prefix")) {
-            std::snprintf(results_prefix, sizeof(results_prefix), "%s", j["results_prefix"].get<std::string>().c_str());
-        } else {
-            strcpy(results_prefix, "");
-        }
-
-        // Construct dataset_name as "<ms_datasetname>_results/<results_prefix>"
-        std::snprintf(dataset_name, sizeof(dataset_name), "%s_results/%s", ms_datasetname, results_prefix);
+        // The group of the microstructure, and where the results go: "<ms_datasetname>_results/<results_prefix>"
+        ms_group     = ms_datasetname.substr(0, ms_datasetname.find_last_of('/') + 1);
+        dataset_name = ms_datasetname + "_results/" + j.value("results_prefix", "");
 
         errorParameters = j["error_parameters"];
         TOL             = errorParameters["tolerance"].get<double>();
@@ -105,25 +70,13 @@ void Reader ::ReadInputFile(const std::string &input_fn)
         problemType = j["problem_type"].get<string>();
         method      = j["method"].get<string>();
 
-        // Parse strain_type (optional, defaults to "small")
-        if (j.contains("strain_type")) {
-            strain_type = j["strain_type"].get<string>();
-            if (strain_type != "small" && strain_type != "large") {
-                throw std::invalid_argument("strain_type must be either 'small' or 'large'");
-            }
-        } else {
-            strain_type = "small"; // Default to small strain
-        }
+        strain_type = j.value("strain_type", "small");
+        if (strain_type != "small" && strain_type != "large")
+            throw std::invalid_argument("strain_type must be either 'small' or 'large'");
 
-        // Parse FE_type (optional, defaults to "HEX8")
-        if (j.contains("FE_type")) {
-            FE_type = j["FE_type"].get<string>();
-            if (FE_type != "HEX8" && FE_type != "HEX8R" && FE_type != "BBAR") {
-                throw std::invalid_argument("FE_type must be one of: 'HEX8', 'HEX8R', or 'BBAR'");
-            }
-        } else {
-            FE_type = "HEX8"; // Default to full integration
-        }
+        FE_type = j.value("FE_type", "HEX8"); // full integration
+        if (FE_type != "HEX8" && FE_type != "HEX8R" && FE_type != "BBAR")
+            throw std::invalid_argument("FE_type must be one of: 'HEX8', 'HEX8R', or 'BBAR'");
 
         resultsToWrite = j["results"].get<vector<string>>(); // Read the results_to_write field
 
@@ -203,7 +156,7 @@ void Reader ::ReadMS(int hm)
     hid_t   space   = H5Dget_space(dset_id);
     hsize_t file_dims[3];
     if (H5Sget_simple_extent_ndims(space) != 3)
-        throw std::runtime_error("Microstructure dataset " + string(ms_datasetname) + " must be 3D");
+        throw std::runtime_error("Microstructure dataset " + ms_datasetname + " must be 3D");
     H5Sget_simple_extent_dims(space, file_dims, nullptr);
     H5Sclose(space);
 
@@ -282,12 +235,6 @@ bool h5_is_zyx(hid_t dset_id)
     return permute_order.empty() || permute_order[0] == 'z' || permute_order[0] == 'Z';
 }
 
-string Reader::MSGroup() const
-{
-    const string path(ms_datasetname);
-    return path.substr(0, path.find_last_of('/') + 1);
-}
-
 std::vector<hsize_t> Reader::DataShape(const string &file, const string &dset_name)
 {
     hid_t                file_id;
@@ -311,10 +258,9 @@ void Reader::FreeMS()
 
 void Reader::OpenResultsFile(const char *output_fn)
 {
-    std::snprintf(results_filename, sizeof(results_filename), "%s", output_fn);
     hid_t plist_id = H5Pcreate(H5P_FILE_ACCESS);
     H5Pset_fapl_mpio(plist_id, communicator, MPI_INFO_NULL);
-    results_file_id = H5Fcreate(results_filename, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
+    results_file_id = H5Fcreate(output_fn, H5F_ACC_TRUNC, H5P_DEFAULT, plist_id);
     H5Pclose(plist_id);
 
     if (results_file_id < 0) {
@@ -328,7 +274,6 @@ void Reader::CloseResultsFile()
         H5Fclose(results_file_id);
         results_file_id = -1;
     }
-    results_filename[0] = '\0';
 }
 
 Reader::~Reader()
diff --git a/src/tangent.cpp b/src/tangent.cpp
new file mode 100644
index 0000000..7c77bcf
--- /dev/null
+++ b/src/tangent.cpp
@@ -0,0 +1,130 @@
+// The homogenized tangent of a converged step.
+//
+// If every material has a tangent C = d flux / d gradient, the tangents of the
+// step, frozen at every Gauss point, are a linear material: flux = C gradient.
+// Its homogenized stiffness is the consistent homogenized tangent, and FANS
+// finds it as for any linear material: a solve per unit macroscopic gradient,
+// whose mean flux is a column.
+// Otherwise the materials themselves are perturbed, by a finite difference. That
+// is not for materials with internal variables: every perturbed solve advances them.
+
+#include "solver.h"
+#include "LargeStrainMechModel.h"
+
+namespace {
+
+// The model base of the problem: thermal, small or large strain
+template <int n_str>
+using Kinematics = std::conditional_t<n_str == 3, ThermalModel, std::conditional_t<n_str == 6, SmallStrainMechModel, LargeStrainMechModel>>;
+
+// The linear material flux = C gradient, C being [element][n_gp][n_str][n_str]
+template <int n_str>
+class FrozenTangent : public Kinematics<n_str> {
+  public:
+    FrozenTangent(const Reader &reader, size_t n_elem)
+        : Kinematics<n_str>(reader), C(n_elem * this->n_gp * n_str * n_str) {}
+
+    vector<double> C;
+
+    void get_sigma(int i, int, ptrdiff_t element_idx) override
+    {
+        const size_t                                      gp = element_idx * this->n_gp + i / n_str;
+        Map<const Matrix<double, n_str, n_str, RowMajor>> C_gp(&C[gp * n_str * n_str]);
+        this->sigma.template segment<n_str>(i).noalias() = C_gp * this->eps.template segment<n_str>(i);
+    }
+
+    // Not needed: the solver has its reference material already
+    Matrix<double, n_str, n_str> get_reference_stiffness() override
+    {
+        return Matrix<double, n_str, n_str>::Zero();
+    }
+};
+
+// The tangent of every material at the converged step, into C
+template <int howmany, int n_str>
+void collect_tangents(Solver<howmany, n_str> &solver, vector<double> &C)
+{
+    const auto  &materials = *solver.matmanager;
+    const size_t per_elem  = C.size() / (solver.local_n0 * solver.n_y * solver.n_z);
+
+    // The native materials element by element, as Solver::get_homogenized_stress goes through them
+    Matrix<double, howmany * 8, 1> ue;
+    solver.template iterateCubes<0>([&](ptrdiff_t *idx, ptrdiff_t *) {
+        const MaterialInfo<howmany, n_str> &info = materials.get_info(solver.ms[idx[0]]);
+        if (info.model->wants_batch())
+            return;
+        for (int i = 0; i < 8; ++i)
+            for (int j = 0; j < howmany; ++j)
+                ue(howmany * i + j) = solver.v_u[howmany * idx[i] + j];
+        info.model->getTangent(&C[idx[0] * per_elem], ue, info.local_mat_id, idx[0]);
+    });
+
+    // The batched materials in batches
+    for (size_t m = 0; m < materials.models.size(); ++m)
+        if (materials.models[m]->wants_batch())
+            materials.models[m]->evaluate_batch(solver.batch_elems[m], solver.ms, solver.batch_ue.data(), solver.batch_gp_stress.data(), C.data());
+}
+
+} // namespace
+
+template <int howmany, int n_str>
+MatrixXd Solver<howmany, n_str>::get_homogenized_tangent(double pert_param)
+{
+    // The step: its stress (evaluating it also readies the nodal values the tangents are
+    // collected from), and what to come back to
+    const VectorXd                   stress    = get_homogenized_stress();
+    MaterialManager<howmany, n_str> *materials = matmanager;
+    const vector<double>             g0        = matmanager->models[0]->macroscale_loading;
+    const vector<double>             u(v_u, v_u + local_n0 * n_y * n_z * howmany);
+
+    // Does every material have a tangent? Linear models with element stiffnesses are
+    // left as they are: perturbing them is exact, and faster.
+    bool consistent = !matmanager->all_stiffness;
+    for (int phase = 0; phase < matmanager->get_num_phases(); ++phase)
+        if (!matmanager->get_info(phase).is_linear && !matmanager->get_info(phase).model->has_tangent())
+            consistent = false;
+
+    // If so, the frozen tangents take the place of the materials
+    std::unique_ptr<MaterialManager<howmany, n_str>> linearised;
+    if (consistent) {
+        auto *frozen = new FrozenTangent<n_str>(reader, local_n0 * n_y * n_z);
+        collect_tangents(*this, frozen->C);
+        linearised = std::make_unique<MaterialManager<howmany, n_str>>(frozen, matmanager->get_num_phases()); // owns `frozen`
+        matmanager = linearised.get();
+        Log::logger().info("# Homogenized tangent from the materials' tangents");
+    } else if (!matmanager->all_linear) {
+        Log::logger().warn("# Homogenized tangent by perturbation: For materials with internal variables, this is WRONG!");
+    }
+
+    const bool linear = matmanager->all_linear;
+    disableMixedBC();
+
+    homogenized_tangent.resize(n_str, n_str);
+    for (int j = 0; j < n_str; ++j) {
+        vector<double> gradient(n_str, 0.0);
+        if (linear) {
+            gradient[j] = 1.0;
+            std::fill_n(v_u, u.size(), 0.0);
+        } else {
+            gradient = g0;
+            gradient[j] += pert_param;
+        }
+        matmanager->set_gradient(gradient);
+        solve();
+
+        if (linear)
+            homogenized_tangent.col(j) = get_homogenized_stress();
+        else
+            homogenized_tangent.col(j) = (get_homogenized_stress() - stress) / pert_param;
+    }
+    homogenized_tangent = 0.5 * (homogenized_tangent + homogenized_tangent.transpose()).eval();
+
+    // Back to the step
+    matmanager = materials;
+    std::copy(u.begin(), u.end(), v_u);
+    return homogenized_tangent;
+}
+
+template MatrixXd Solver<1, 3>::get_homogenized_tangent(double);
+template MatrixXd Solver<3, 6>::get_homogenized_tangent(double);
+template MatrixXd Solver<3, 9>::get_homogenized_tangent(double);