diff --git a/include/testing.hh b/include/testing.hh
index 53bc571..2395db3 100644
--- a/include/testing.hh
+++ b/include/testing.hh
@@ -4,6 +4,7 @@
 #include <general.hh>
 #include <config_file.hh>
 #include <grid_fft.hh>
+#include <cosmology_calculator.hh>
 
 namespace testing{
     void output_potentials_and_densities( 
@@ -27,6 +28,7 @@ namespace testing{
 
     void output_convergence(
         ConfigFile &the_config,
+        CosmologyCalculator* the_cosmo_calc,
         std::size_t ngrid, real_t boxlen, real_t vfac, real_t dplus,
         Grid_FFT<real_t> &phi,
         Grid_FFT<real_t> &phi2,
diff --git a/src/ic_generator.cc b/src/ic_generator.cc
index a915db9..bc2a956 100644
--- a/src/ic_generator.cc
+++ b/src/ic_generator.cc
@@ -335,7 +335,7 @@ int Run( ConfigFile& the_config )
             } else if(testing == "velocity_displacement_symmetries") {
                 testing::output_velocity_displacement_symmetries(the_config, ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3);
             } else if(testing == "convergence") {
-                testing::output_convergence(the_config, ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3);
+                testing::output_convergence(the_config, the_cosmo_calc.get(), ngrid, boxlen, vfac, Dplus0, phi, phi2, phi3a, phi3b, A3);
             } else {
                 csoca::flog << "unknown test '" << testing << "'" << std::endl;
                 std::abort();
diff --git a/src/testing.cc b/src/testing.cc
index bfd088d..533855a 100644
--- a/src/testing.cc
+++ b/src/testing.cc
@@ -242,6 +242,7 @@ void output_velocity_displacement_symmetries(
 
 void output_convergence(
     ConfigFile &the_config,
+    CosmologyCalculator* the_cosmo_calc,
     std::size_t ngrid, real_t boxlen, real_t vfac, real_t dplus,
     Grid_FFT<real_t> &phi,
     Grid_FFT<real_t> &phi2,
@@ -249,7 +250,6 @@ void output_convergence(
     Grid_FFT<real_t> &phi3b,
     std::array<Grid_FFT<real_t> *, 3> &A3)
 {
-
     // scale all potentials to remove dplus0
     phi /= dplus;
     phi2 /= dplus * dplus;
@@ -259,6 +259,90 @@ void output_convergence(
     (*A3[1]) /= dplus * dplus * dplus;
     (*A3[2]) /= dplus * dplus * dplus;
 
+    ////////////////////// theoretical convergence radius //////////////////////
+
+    // compute phi_code
+    Grid_FFT<real_t> phi_code({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
+    phi_code.FourierTransformForward(false);
+    #pragma omp parallel for collapse(3)
+    for (std::size_t i = 0; i < phi_code.size(0); ++i) {
+        for (std::size_t j = 0; j < phi_code.size(1); ++j) {
+            for (std::size_t k = 0; k < phi_code.size(2); ++k) {
+                std::size_t idx = phi_code.get_idx(i, j, k);
+                phi_code.kelem(idx) = -phi.kelem(idx);
+            }
+        }
+    }
+
+    // initialize norm to 0
+    Grid_FFT<real_t> nabla_vini_norm({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
+    #pragma omp parallel for collapse(3)
+    for (std::size_t i = 0; i < nabla_vini_norm.size(0); ++i) {
+        for (std::size_t j = 0; j < nabla_vini_norm.size(1); ++j) {
+            for (std::size_t k = 0; k < nabla_vini_norm.size(2); ++k) {
+                std::size_t idx = nabla_vini_norm.get_idx(i, j, k);
+                nabla_vini_norm.relem(idx) = 0.0;
+            }
+        }
+    }
+
+    Grid_FFT<real_t> nabla_vini_mn({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
+    for(std::size_t m = 0; m < 3; m++) {
+        for(std::size_t n = m; n < 3; n++) {
+            nabla_vini_mn.FourierTransformForward(false);
+            #pragma omp parallel for collapse(3)
+            for (std::size_t i = 0; i < phi_code.size(0); ++i) {
+                for (std::size_t j = 0; j < phi_code.size(1); ++j) {
+                    for (std::size_t k = 0; k < phi_code.size(2); ++k) {
+                        std::size_t idx = phi_code.get_idx(i, j, k);
+                        auto kk = phi_code.get_k<real_t>(i, j, k);
+                        nabla_vini_mn.kelem(idx) = phi_code.kelem(idx) * (kk[m] * kk[n]);
+                    }
+                }
+            }
+            nabla_vini_mn.FourierTransformBackward();
+            nabla_vini_mn *= (3.2144004915 / the_cosmo_calc->CalcGrowthFactor(1.0));
+            // sum of squares
+            #pragma omp parallel for collapse(3)
+            for (std::size_t i = 0; i < nabla_vini_norm.size(0); ++i) {
+                for (std::size_t j = 0; j < nabla_vini_norm.size(1); ++j) {
+                    for (std::size_t k = 0; k < nabla_vini_norm.size(2); ++k) {
+                        std::size_t idx = nabla_vini_norm.get_idx(i, j, k);
+                        if(m != n) {
+                            nabla_vini_norm.relem(idx) += (2.0 * nabla_vini_mn.relem(idx) * nabla_vini_mn.relem(idx));
+                        } else {
+                            nabla_vini_norm.relem(idx) += (nabla_vini_mn.relem(idx) * nabla_vini_mn.relem(idx));
+                        }
+                    }
+                }
+            }
+        }
+    }
+    // square root
+    #pragma omp parallel for collapse(3)
+    for (std::size_t i = 0; i < nabla_vini_norm.size(0); ++i) {
+        for (std::size_t j = 0; j < nabla_vini_norm.size(1); ++j) {
+            for (std::size_t k = 0; k < nabla_vini_norm.size(2); ++k) {
+                std::size_t idx = nabla_vini_norm.get_idx(i, j, k);
+                nabla_vini_norm.relem(idx) = std::sqrt(nabla_vini_norm.relem(idx));
+            }
+        }
+    }
+
+    // get t_eds
+    Grid_FFT<real_t> t_eds({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
+    #pragma omp parallel for collapse(3)
+    for (std::size_t i = 0; i < t_eds.size(0); ++i) {
+        for (std::size_t j = 0; j < t_eds.size(1); ++j) {
+            for (std::size_t k = 0; k < t_eds.size(2); ++k) {
+                std::size_t idx = t_eds.get_idx(i, j, k);
+                t_eds.relem(idx) = 0.0204 / nabla_vini_norm.relem(idx);
+            }
+        }
+    }
+
+    ////////////////////////// 3lpt convergence test ///////////////////////////
+
     // initialize grids to 0
     Grid_FFT<real_t> psi_1({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
     Grid_FFT<real_t> psi_2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
@@ -351,13 +435,17 @@ void output_convergence(
         }
     }
 
-    // write results
-    unlink("convergence_test.hdf5");
-    inv_convergence_radius.Write_to_HDF5("convergence_test.hdf5", "inv_convergence_radius");
-    psi_1.Write_to_HDF5("convergence_test.hdf5", "psi_1_norm");
-    psi_2.Write_to_HDF5("convergence_test.hdf5", "psi_2_norm");
-    psi_3.Write_to_HDF5("convergence_test.hdf5", "psi_3_norm");
-
+    ////////////////////////////// write results ///////////////////////////////
+    std::string convergence_test_filename("convergence_test.hdf5");
+    unlink(convergence_test_filename.c_str());
+#if defined(USE_MPI)
+    MPI_Barrier(MPI_COMM_WORLD);
+#endif
+    t_eds.Write_to_HDF5(convergence_test_filename, "t_eds");
+    inv_convergence_radius.Write_to_HDF5(convergence_test_filename, "inv_convergence_radius");
+    // psi_1.Write_to_HDF5(convergence_test_filename, "psi_1_norm");
+    // psi_2.Write_to_HDF5(convergence_test_filename, "psi_2_norm");
+    // psi_3.Write_to_HDF5(convergence_test_filename, "psi_3_norm");
 }
 
 } // namespace testing