diff --git a/check/CMakeLists.txt b/check/CMakeLists.txt
index 645e28c456293d0b57fceefb8fe5dc975996b008..6ac9d220719b010f964b2f267bd6e8bfc1b3b403 100644
--- a/check/CMakeLists.txt
+++ b/check/CMakeLists.txt
@@ -24,7 +24,6 @@ set(source_check_files
     check_2d.cpp
 
     # Diwina
-    test-rand-3.cpp
     test-rand-3_view.cpp
 )
 
diff --git a/check/check_1d.cpp b/check/check_1d.cpp
index 25bd5eed0322ed7fbbc23b3f72eefd42f3fa5a9f..3bde180add1aaaa7a6bd2531a2af8ba063e5771b 100644
--- a/check/check_1d.cpp
+++ b/check/check_1d.cpp
@@ -289,7 +289,7 @@ auto run(const int& tree_height, const int& group_size, const std::size_t order,
                                           [&box, tree_height](auto& leaf)
                                           {
                                               const auto nb_elt = leaf.size();
-                                              scalfmm::utils::print_leaf(leaf, box, tree_height - 1);
+                                              scalfmm::io::print_leaf(leaf, box, tree_height - 1);
                                               std::cout << "-----\n\n";
                                           });
     }
diff --git a/check/check_2d.cpp b/check/check_2d.cpp
index 5a6fd01edb8f81e303be5f9e46ac99381d4aeb70..5cc05c80e21ea43f8ca69c107ff7a7fb0ad1eea7 100644
--- a/check/check_2d.cpp
+++ b/check/check_2d.cpp
@@ -292,7 +292,7 @@ auto run(const int &tree_height, const int &group_size, const std::size_t order,
             scalfmm::component::for_each_leaf(std::begin(tree), std::end(tree),
                                               [&box, tree_height](auto& leaf)
                                               {
-                                                  scalfmm::utils::print_leaf(leaf, box, tree_height - 1);
+                                                  scalfmm::io::print_leaf(leaf, box, tree_height - 1);
                                                   std::cout << "-----\n\n";
                                               });
         }
diff --git a/check/check_interaction_lists.cpp b/check/check_interaction_lists.cpp
index d1a6174255d1dae153790febfd731bf2c72a8f01..02888e23f2e2aae0534a06f97eb843d38948ce0b 100644
--- a/check/check_interaction_lists.cpp
+++ b/check/check_interaction_lists.cpp
@@ -116,7 +116,7 @@ template<typename Tree>
 auto print_leaves(Tree const& tree) -> void
 {
     scalfmm::component::for_each_leaf(std::cbegin(tree), std::cend(tree),
-                                      [&tree](auto& leaf) { scalfmm::utils::print_leaf(leaf); });
+                                      [&tree](auto& leaf) { scalfmm::io::print_leaf(leaf); });
 }
 
 template<typename TREE, typename ARRAY1, typename ARRAY2>
diff --git a/check/test-rand-3.cpp b/check/test-rand-3.cpp
deleted file mode 100644
index 04137049788e9ed3ff51fd48c86af30d7b7e8b6b..0000000000000000000000000000000000000000
--- a/check/test-rand-3.cpp
+++ /dev/null
@@ -1,428 +0,0 @@
-// @FUSE_OMP
-/**
- * @file test-rand-3.cpp
- * @author your name (you@domain.com)
- * @brief
- * @version 0.1
- * @date 2023-02-10
- *
- * @copyright Copyright (c) 2023
- *
- */
-
-/**
- *\code
- * check/Release/test-rand-3 -th 3 -gs 10  --order 6 --scale 1.0,1.0,1.0
- *\endcode
- */
-#include "scalfmm/container/particle.hpp"
-#include "scalfmm/container/particle_container.hpp"
-#include "scalfmm/container/point.hpp"
-//
-//
-// #include "scalfmm/meta/type_pack.hpp"
-//
-#include "scalfmm/algorithms/fmm.hpp"
-#include "scalfmm/algorithms/full_direct.hpp"
-#include "scalfmm/algorithms/sequential/sequential.hpp"
-#include "scalfmm/interpolation/interpolation.hpp"
-#include "scalfmm/matrix_kernels/laplace.hpp"
-#include "scalfmm/operators/fmm_operators.hpp"
-// #include "scalfmm/options/options.hpp"
-//
-// Tree
-#include "scalfmm/interpolation/grid_storage.hpp"
-#include "scalfmm/tree/box.hpp"
-#include "scalfmm/tree/cell.hpp"
-#include "scalfmm/tree/group_tree.hpp"
-#include "scalfmm/tree/leaf.hpp"
-#include "scalfmm/tree/utils.hpp"
-
-//
-#include "scalfmm/meta/utils.hpp"
-#include "scalfmm/tools/fma_loader.hpp"
-
-// #include "scalfmm/utils/io_helpers.hpp"
-#include "scalfmm/utils/source_target.hpp"
-//
-#include "parameters.hpp"
-
-#include <array>
-#include <chrono>
-#include <iostream>
-#include <random>
-#include <string>
-#include <tuple>
-
-#include <cpp_tools/cl_parser/help_descriptor.hpp>
-#include <cpp_tools/cl_parser/tcli.hpp>
-#include <cpp_tools/colors/colorized.hpp>
-#include <cpp_tools/timers/simple_timer.hpp>
-
-// using namespace scalfmm::io;
-
-//
-namespace local_args
-{
-    struct output_file
-    {
-        cpp_tools::cl_parser::str_vec flags = {"--output-file", "-fout"};
-        const char* description = "Output particle file (with extension .fma (ascii) or bfma (binary).";
-        using type = std::string;
-    };
-    struct ns
-    {
-        cpp_tools::cl_parser::str_vec flags = {"--number_sources", "-ns"};
-        const char* description = "Number of source particles.";
-        using type = int;
-        type def = 2000000;
-    };
-    struct nt
-    {
-        cpp_tools::cl_parser::str_vec flags = {"--number_targets", "-nt"};
-        const char* description = "Number of target particles.";
-        using type = int;
-        type def = 200000;
-    };
-    struct seed
-    {
-        cpp_tools::cl_parser::str_vec flags = {"-seed"};
-        const char* description = "Seed to generate random particles.";
-        using type = int;
-        type def = 1323;
-    };
-    struct scale
-    {
-        cpp_tools::cl_parser::str_vec flags = {"--scale", "-s"}; /*!< The flags */
-        std::string description = "Scale parameter sx,s,sz";     /*!< The description */
-        std::string input_hint = "value,value,value";            /*!< The description */
-        using type = std::vector<double>;
-    };
-    struct check
-    {
-        cpp_tools::cl_parser::str_vec flags = {"--check"};
-        const char* description = "Check with p2p ";
-        using type = bool;
-        /// The parameter is a flag, it doesn't expect a following value
-        enum
-        {
-            flagged
-        };
-    };
-}   // namespace local_args
-
-using timer_type = cpp_tools::timers::timer<std::chrono::seconds>;
-
-template<typename Tree, typename Value_type>
-auto sum_output(Tree const& tree, const Value_type& Q) -> void
-{
-    Value_type potential = Value_type(0.);
-    Value_type fx = 0.0, fy = 0.0, fz = 0.0;
-    int idx{0};
-    scalfmm::component::for_each_leaf(std::cbegin(tree), std::cend(tree),
-                                      [&potential, &fx, &fy, &fz, &idx](auto& leaf)
-                                      {
-                                          const auto& container = leaf.cparticles();
-                                          const auto nb_elt = std::distance(container.begin(), container.end());
-                                          for(std::size_t i = 0; i < nb_elt; ++i)
-                                          {
-                                              const auto& p = container.particle(i);
-                                              auto output = p.outputs();
-                                              potential += output.at(0);
-                                              fx += output.at(1);
-                                              fy += output.at(2);
-                                              fz += output.at(3);
-                                              std::cout << idx++ << " " << p << std::endl;
-                                          }
-                                      });
-    std::cout << "Forces Sum / Q  x = " << fx / Q << " y = " << fy / Q << " z = " << fz / Q << std::endl;
-    std::cout << "Potential / Q    = " << potential / Q << std::endl;
-    std::cout << std::endl;
-}
-
-template<int dimension, typename value_type, class fmm_operators_type>
-auto fmm_run(const int& tree_height, const int& group_size, const int& order, const int& seed,
-             const std::vector<value_type>& Scale_v, const int& NbSources, const int NbTargets, bool check_direct,
-             const std::string& output_file) -> int
-{
-    bool display_container = false;
-    bool display_tree = false;
-
-    //  The matrix kernel
-    using near_matrix_kernel_type = typename fmm_operators_type::near_field_type::matrix_kernel_type;
-    using far_field_type = typename fmm_operators_type::far_field_type;
-    using interpolator_type = typename far_field_type::approximation_type;
-
-    using far_matrix_kernel_type = typename interpolator_type::matrix_kernel_type;
-    static constexpr std::size_t nb_inputs{near_matrix_kernel_type::km};
-    static constexpr std::size_t nb_outputs{near_matrix_kernel_type::kn};
-    //
-    near_matrix_kernel_type mk{};
-    //
-    // Open particles files
-
-    timer_type time{};
-
-    constexpr int zeros{nb_inputs};   // should be zero
-    using point_type = scalfmm::container::point<value_type, dimension>;
-    using particle_source_type =
-      scalfmm::container::particle<value_type, dimension, value_type, nb_inputs, value_type, nb_outputs, std::size_t>;
-    using particle_target_type =
-      scalfmm::container::particle<value_type, dimension, value_type, zeros, value_type, nb_outputs, std::size_t>;
-    // Construct the container of particles
-    using container_source_type = scalfmm::container::particle_container<particle_source_type>;
-    using container_target_type = scalfmm::container::particle_container<particle_target_type>;
-    using box_type = scalfmm::component::box<point_type>;
-    //
-    using leaf_source_type = scalfmm::component::leaf<particle_source_type>;
-    using leaf_target_type = scalfmm::component::leaf<particle_target_type>;
-    //
-    using cell_type = scalfmm::component::cell<typename interpolator_type::storage_type>;
-    using tree_source_type = scalfmm::component::group_tree<cell_type, leaf_source_type, box_type>;
-    using tree_target_type = scalfmm::component::group_tree<cell_type, leaf_target_type, box_type>;
-    ///////////////////////////////////////////////////////////////////////////////////////////////////
-    std::cout << "Working on particles ..." << std::endl;
-    std::mt19937_64 gen(seed);
-    std::uniform_real_distribution<value_type> dist(-1, 1);
-    const auto ScaleX = Scale_v[0];
-    const auto ScaleY = Scale_v[1];
-    const auto ScaleZ = Scale_v[2];
-    auto Scale = ScaleX;
-    Scale = (ScaleY > Scale ? ScaleY : Scale);
-    Scale = (ScaleZ > Scale ? ScaleZ : Scale);
-    std::cout << cpp_tools::colors::green << "Creating & Inserting ...\n" << cpp_tools::colors::reset;
-    //
-    std::cout << cpp_tools::colors::blue << "Entering tree test...\n" << cpp_tools::colors::reset;
-
-    time.tic();
-    point_type box_center_target{0., 0., 0.};
-    value_type box_width_target{2.01 * Scale};
-    container_target_type container_target(NbTargets);
-    box_type box_target(box_width_target, box_center_target);
-
-    particle_target_type p;
-    for(int i = 0; i < NbTargets; ++i)
-    {   // cibles
-        // double xTarget = (value_type(2.) * drand48() - value_type(1.)) * ScaleX;
-        // double yTarget = (value_type(2.) * drand48() - value_type(1.)) * ScaleY;
-        double xTarget = dist(gen) * ScaleX;
-        double yTarget = dist(gen) * ScaleY;
-        if constexpr(dimension == 3)
-        {
-            // double zTarget = (value_type(2.) * drand48() - value_type(1.)) * ScaleZ;
-            double zTarget = dist(gen) * ScaleZ;
-            p.position() = {xTarget, yTarget, zTarget};
-        }
-        else
-        {
-            p.position() = {xTarget, yTarget};
-        }
-        p.inputs(0) = 1.0;
-        p.variables(i);
-        // std::cout << "target " << i << " " << p << std::endl;
-        container_target.insert_particle(i, p);
-    }
-
-    point_type box_center_source{0., 0., 0.};
-    value_type box_width_source{2.01 * Scale};
-    container_source_type container_source(NbSources);
-    box_type box_source(box_width_source, box_center_source);
-    //
-    for(int i = 0; i < NbSources; ++i)
-    {   // sources
-        particle_source_type p;
-        // double xSource = 2. * drand48() - 1.;
-        // double ySource = 2. * drand48() - 1.;
-        double xSource = dist(gen);
-        double ySource = dist(gen);
-        if constexpr(dimension == 3)
-        {
-            // double zSource = 2. * drand48() - 1.;
-            double zSource = dist(gen);
-            p.position() = {xSource, ySource, zSource};
-        }
-        else
-        {
-            p.position() = {xSource, ySource};
-        }
-        p.inputs(0) = 1.0;
-        p.variables(i);
-        // std::cout << "source " << i << " " << p << std::endl;
-        container_source.insert_particle(i, p);
-    }
-
-    time.tac();
-    std::cout << cpp_tools::colors::yellow << "Container loaded in " << time.elapsed() << " s\n"
-              << cpp_tools::colors::reset;
-
-    auto box = scalfmm::utils::bounding_box(box_source, box_target);
-    // auto container_all = scalfmm::utils::merge(container_source, container_target) ;
-
-    // build trees
-    bool sorted = false;
-    time.tic();
-
-    tree_source_type tree_source(tree_height, order, box, group_size, group_size, container_source, sorted);
-
-    tree_target_type tree_target(tree_height, order, box, group_size, group_size, container_target, sorted);
-    time.tac();
-    std::cout << "Done  "
-              << "(@Building trees = " << time.elapsed() << " s)." << std::endl;
-
-    std::cout << "bounding_box " << box << std::endl;
-    std::cout << cpp_tools::colors::green << "... Done.\n" << cpp_tools::colors::reset;
-
-    /////////////////////////////////////////////////////////////////////////////////////
-    //
-    //              Compute source-target interaction though FMM
-    //
-    /////////////////////////////////////////////////////////////////////////////////////
-    auto box_width = box.width(0);
-    // Far field
-    far_matrix_kernel_type mk_far{};
-    interpolator_type interpolator(mk_far, order, tree_height, box_width);
-    typename fmm_operators_type::far_field_type far_field(interpolator);
-    // Near field
-    near_matrix_kernel_type mk_near{};
-    typename fmm_operators_type::near_field_type near_field(mk_near, false /* no mutual for source!= target*/);
-    //
-    std::cout << cpp_tools::colors::blue << "Fmm with kernels: " << std::endl
-              << "       near " << mk_near.name() << std::endl
-              << "       far  " << mk_far.name() << std::endl;
-
-    fmm_operators_type fmm_operator(near_field, far_field);
-    auto neighbour_separation = fmm_operator.near_field().separation_criterion();
-    //
-    tree_target.build_interaction_lists(tree_source, neighbour_separation, false);
-    auto operator_to_proceed = scalfmm::algorithms::all;
-    value_type Q;
-    for(int iter = 0; iter < 1; ++iter)
-    {
-        Q = pow(2.0, double(iter));
-        std::cout << "iter = " << iter << " Q = " << Q << std::endl;
-        // reset outputs
-        tree_target.reset_outputs();
-        tree_target.reset_locals();
-        tree_source.reset_multipoles();
-        //
-        scalfmm::component::for_each_leaf(std::cbegin(tree_source), std::cend(tree_source),
-                                          [&Q](auto& leaf)
-                                          {
-                                              auto& container = leaf.particles();
-                                              const auto nb_elt = std::distance(container.begin(), container.end());
-                                              auto it = std::begin(container);
-                                              using proxy_type = typename std::decay_t<decltype(container)>::proxy_type;
-                                              for(auto i = 0; i < nb_elt; ++it, ++i)
-                                              {
-                                                  auto proxy = proxy_type(*it);
-                                                  proxy.inputs(0) = Q;
-                                              }
-                                          });
-        //
-        std::cout << "Start algo\n";
-        time.tic();
-        scalfmm::algorithms::fmm[scalfmm::options::_s(scalfmm::options::timit)](tree_source, tree_target, fmm_operator,
-                                                                                operator_to_proceed);
-
-        // scalfmm::algorithms::sequential::sequential(tree_source, tree_target, fmm_operator, operator_to_proceed);
-        time.tac();
-        std::cout << "Done  " << cpp_tools::colors::blue << "(@Algorithm = " << time.elapsed() << " s)." << std::endl
-                  << cpp_tools::colors::reset;
-        ;
-        sum_output(tree_target, Q);
-    }
-    if(check_direct)
-    {
-        Q = 1;
-        std::cout << cpp_tools::colors::green << "full interaction computation  with kernel: " << mk_near.name()
-                  << std::endl
-                  << cpp_tools::colors::reset;
-        time.tic();
-        scalfmm::algorithms::full_direct(container_source, container_target, mk_near);
-        time.tac();
-        std::cout << cpp_tools::colors::green << "... Done.\n" << cpp_tools::colors::reset;
-        std::cout << cpp_tools::colors::yellow << "Computation done in " << time.elapsed() << " ms\n"
-                  << cpp_tools::colors::reset;
-        //
-        value_type potential = value_type(0.);
-        value_type fx = value_type(0.), fy = value_type(0.), fz = value_type(0.);
-        //
-        for(int i{0}; i < container_target.size(); ++i)
-        {
-            std::cout << i << " " << container_target.particle(i) << std::endl;
-
-            auto output = container_target.outputs(i);
-            potential += output.at(0);
-            fx += output.at(1);
-            fy += output.at(2);
-            if constexpr(dimension == 3)
-            {
-                fz += output.at(3);
-            }
-        }
-        std::cout << "Forces Sum / Q  x = " << fx / Q << " y = " << fy / Q << " z = " << fz / Q << std::endl;
-        std::cout << "Potential / Q    = " << potential / Q << std::endl;
-        std::cout << std::endl;
-    }
-
-    std::string outputFile{"random3_target_old.fma"};
-    std::cout << "Write targets in " << outputFile << std::endl;
-    scalfmm::tools::FFmaGenericWriter<double> writer_t(outputFile);
-    writer_t.writeDataFromTree(tree_target, NbTargets);
-    return 0;
-}
-
-auto main([[maybe_unused]] int argc, [[maybe_unused]] char* argv[]) -> int
-{
-    using value_type = double;
-    constexpr int dimension = 3;
-
-    //
-    // Parameter handling
-    auto parser = cpp_tools::cl_parser::make_parser(
-      cpp_tools::cl_parser::help{}, args::tree_height{}, args::block_size{}, args::order{}, local_args::output_file(),
-      local_args::check{}, local_args::ns{}, local_args::nt{}, local_args::scale{}, local_args::seed{});
-    parser.parse(argc, argv);
-
-    const int seed{parser.get<local_args::seed>()};
-    const int NbSources{parser.get<local_args::ns>()};
-    const int NbTargets{parser.get<local_args::nt>()};
-    std::cout << cpp_tools::colors::blue << "<params> seed:      " << seed << cpp_tools::colors::reset << '\n';
-    std::cout << cpp_tools::colors::blue << "<params> NbSources: " << NbSources << cpp_tools::colors::reset << '\n';
-    std::cout << cpp_tools::colors::blue << "<params> NbTargets: " << NbTargets << cpp_tools::colors::reset << '\n';
-    std::vector<value_type> scale(3, 1.0);
-    scale = parser.get<local_args::scale>();
-
-    // const unsigned int NbThreads = FParameters::getValue(argc, argv, "-t", omp_get_max_threads());
-
-    const int NbLevels{parser.get<args::tree_height>()};
-    std::cout << cpp_tools::colors::blue << "<params> Tree height : " << NbLevels << cpp_tools::colors::reset << '\n';
-
-    const int group_size{parser.get<args::block_size>()};
-    const auto order{parser.get<args::order>()};
-
-    std::cout << cpp_tools::colors::blue << "<params> Group Size : " << group_size << cpp_tools::colors::reset << '\n';
-    const auto output_file{parser.get<local_args::output_file>()};
-    std::cout << cpp_tools::colors::blue << "<params> Runtime order : " << order << cpp_tools::colors::reset << '\n';
-
-    bool check_direct{parser.exists<local_args::check>()};
-    if(!output_file.empty())
-    {
-        std::cout << cpp_tools::colors::blue << "<params> Output file : " << output_file << cpp_tools::colors::reset
-                  << '\n';
-    }
-    //
-    //
-    using near_matrix_kernel_type = scalfmm::matrix_kernels::laplace::val_grad_one_over_r<dimension>;
-    using near_field_type = scalfmm::operators::near_field_operator<near_matrix_kernel_type>;
-    using far_matrix_kernel_type = near_matrix_kernel_type;   // scalfmm::matrix_kernels::laplace::one_over_r;
-    //
-    using options = scalfmm::options::uniform_<scalfmm::options::fft_>;
-
-    using interpolation_type = scalfmm::interpolation::interpolator<value_type, 3, far_matrix_kernel_type, options>;
-    using far_field_type = scalfmm::operators::far_field_operator<interpolation_type>;
-
-    using fmm_operators_type = scalfmm::operators::fmm_operators<near_field_type, far_field_type>;
-    fmm_run<3, value_type, fmm_operators_type>(NbLevels, group_size, order, seed, scale, NbSources, NbTargets,
-                                               check_direct, output_file);
-}
diff --git a/docs/user_guide.org b/docs/user_guide.org
index 9965ed696b96a41c5f5735cefa6e07bdc30d49a1..f825472b60463c0fed331a13cf870827a727b200 100644
--- a/docs/user_guide.org
+++ b/docs/user_guide.org
@@ -311,9 +311,9 @@ This is the most efficient method (x3 compared to the previous one) and can be p
 
 We propose several functions in the namesapce ~scalfmm::io~ to
 
-- Display  a leaf  ~print_leaf(leaf)~, we display symbolic information and the particles inside the leaf. 
-- isplay  a cell ~print_cell(cell)~, we display symbolic information and the multipole(s) and local(s) arrays inside the lcell. 
-- Display a tuple (formatted display) 
+- Display  a leaf  ~print_leaf(leaf)~, we display symbolic information and the particles inside the leaf (in ~tree/utils.hpp~).
+- Display  a cell ~print_cell(cell)~, we display symbolic information and the multipole(s) and local(s) arrays inside the cell (in ~tree/utils.hpp~). 
+- Display a tuple in a formatted display (in ~utils/io_helpers.hpp~)
 #+begin_src C++
   std::tuple<double, int> t(0.5, 3);
   io::print()
diff --git a/examples/fmm_source_target.cpp b/examples/fmm_source_target.cpp
index 7f0af31b0d3100c0f857f27b478e436937aac77c..8c3af0a115826290b1f5a3cb8443a0add1449806 100644
--- a/examples/fmm_source_target.cpp
+++ b/examples/fmm_source_target.cpp
@@ -182,7 +182,7 @@ template<typename Tree>
 auto print_leaves(Tree const& tree) -> void
 {
     scalfmm::component::for_each_leaf(std::cbegin(tree), std::cend(tree),
-                                      [&tree](auto& leaf) { scalfmm::utils::print_leaf(leaf); });
+                                      [&tree](auto& leaf) { scalfmm::io::print_leaf(leaf); });
 }
 template<typename Tree, typename Container>
 auto check_output(Container const& part, Tree const& tree)
diff --git a/include/scalfmm/tree/utils.hpp b/include/scalfmm/tree/utils.hpp
index 50c7055dd2803c7dbf05d2d6321142c8a01bf46d..10181a62692ff1fd4900dc452ebc011a30acf5bf 100644
--- a/include/scalfmm/tree/utils.hpp
+++ b/include/scalfmm/tree/utils.hpp
@@ -811,7 +811,7 @@ namespace scalfmm::index
         return shift;
     }
 }   // namespace scalfmm::index
-namespace scalfmm::utils
+namespace scalfmm::io
 {
 
     template<typename Cell>
@@ -878,27 +878,29 @@ namespace scalfmm::utils
             std::cout << i++ << " morton: " << morton << " part= " << pp << std::endl;
         }
     }
-
-           ///
-        /// \brief trace the index of the cells and leaves in the tree
-        ///
-        /// Depending on the level we print more or less details
-        ///  level_trace = 1 print minimal information (height, order, group size)
-        ///  level_trace = 2 print information of the tree (group interval and index inside)
-        ///  level_trace = 3 print information of the tree (leaf interval and index inside and their p2p interaction
-        ///  list) level_trace = 4 print information of the tree (cell interval and index inside and their m2l
-        ///  interaction list)
-        /// level_trace = 5 print information of the tree (leaf and cell interval and index inside
-        ///  and their p2p and m2l interaction lists)
-        ///
-        /// @warning to have the right p2p list we have to have the group size in the tree equal to the number
-        ///  of leaves otherwise, we only print the index inside the group.
-        ///
-        /// @param[in] level_trace level of the trace
-        ///
-        template< typename TreeType>
-        inline auto trace(std::ostream& os, const TreeType &tree, const std::size_t level_trace = 0) -> void 
-        {
+}   // namespace scalfmm::io
+namespace scalfmm::utils
+{
+    ///
+    /// \brief trace the index of the cells and leaves in the tree
+    ///
+    /// Depending on the level we print more or less details
+    ///  level_trace = 1 print minimal information (height, order, group size)
+    ///  level_trace = 2 print information of the tree (group interval and index inside)
+    ///  level_trace = 3 print information of the tree (leaf interval and index inside and their p2p interaction
+    ///  list) level_trace = 4 print information of the tree (cell interval and index inside and their m2l
+    ///  interaction list)
+    /// level_trace = 5 print information of the tree (leaf and cell interval and index inside
+    ///  and their p2p and m2l interaction lists)
+    ///
+    /// @warning to have the right p2p list we have to have the group size in the tree equal to the number
+    ///  of leaves otherwise, we only print the index inside the group.
+    ///
+    /// @param[in] level_trace level of the trace
+    ///
+    template<typename TreeType>
+    inline auto trace(std::ostream& os, const TreeType& tree, const std::size_t level_trace = 0) -> void
+    {
         std::cout << "Trace of the group tree\n";
 
         auto level_0 = []() {};
@@ -1099,7 +1101,7 @@ namespace scalfmm::utils
 
         default:
             level_0();
-            }
         }
+    }
 }   // namespace scalfmm::utils
 #endif   // SCALFMM_TREE_UTILS_HPP
diff --git a/include/scalfmm/utils/sort.hpp b/include/scalfmm/utils/sort.hpp
index 064aca148db8f23391b4bcbed74056c892fa3439..81dfc823b76d2c617dfe61cb4f6aed4625e2a140 100644
--- a/include/scalfmm/utils/sort.hpp
+++ b/include/scalfmm/utils/sort.hpp
@@ -9,6 +9,7 @@
 #include <utility>
 #include <vector>
 
+#include "scalfmm/container/particle_container.hpp"
 #include "scalfmm/meta/utils.hpp"
 #include "scalfmm/tree/for_each.hpp"
 #include "scalfmm/tree/utils.hpp"
@@ -78,19 +79,51 @@ namespace scalfmm
         /// \param[inout] particle_container the particles container
         ///
         ///
-        template<class BOX_T, class CONTAINER_T, typename Int>
-        auto inline sort_container(const BOX_T& box, const Int& level, CONTAINER_T*& particle_container) -> void
+        template<class BOX_T, class ParticleT, typename Int>
+        auto inline sort_container(const BOX_T& box, const Int& level,
+                                   scalfmm::container::particle_container<ParticleT>*& particle_container) -> void
         {
+            using CONTAINER_T = scalfmm::container::particle_container<ParticleT>;
             auto perm = get_morton_permutation(box, level, *particle_container);
             const auto number_of_particles = particle_container->size();
             CONTAINER_T* ordered_container = new CONTAINER_T(number_of_particles);
             for(std::size_t i = 0; i < number_of_particles; ++i)
             {
-                ordered_container->insert_particle(i, particle_container->particle(std::get<0>(perm[i])));
+                ordered_container->insert_particle(i, particle_container->particle(std::get<1>(perm[i])));
             }
             delete particle_container;
             particle_container = ordered_container;
         }
+        template<class BOX_T, class ParticleT, typename Int>
+        auto inline sort_container(const BOX_T& box, const Int& level,
+                                   scalfmm::container::particle_container<ParticleT>& particle_container) -> void
+        {
+            using CONTAINER_T = scalfmm::container::particle_container<ParticleT>;
+            auto perm = get_morton_permutation(box, level, particle_container);
+            const auto number_of_particles = particle_container.size();
+            CONTAINER_T tmp_container(number_of_particles);
+            // for(std::size_t i = 0; i < number_of_particles; ++i)
+            // {
+            //     particle_container.insert_particle(i, particle_container.particle(std::get<1>(perm[i])));
+            // }
+            std::copy(particle_container.begin(), particle_container.end(), tmp_container.begin());
+            for(std::size_t i = 0; i < number_of_particles; ++i)
+            {
+                particle_container.insert_particle(i, tmp_container.particle(std::get<1>(perm[i])));
+            }
+        }
+        template<class BOX_T, class CONTAINER_T, typename Int>
+        auto inline sort_container(const BOX_T& box, const Int& level, CONTAINER_T& particle_container) -> void
+        {
+            auto perm = get_morton_permutation(box, level, particle_container);
+            const auto number_of_particles = particle_container.size();
+            CONTAINER_T tmp_container(number_of_particles);
+            std::copy(particle_container.begin(), particle_container.end(), tmp_container.begin());
+            for(std::size_t i = 0; i < number_of_particles; ++i)
+            {
+                particle_container[i] = tmp_container[std::get<1>(perm[i])];
+            }
+        }
 
         /// \brief sort a vector particles
         ///