diff --git a/CMakeLists.txt b/CMakeLists.txt
index 9d0d32bc1..abea4d695 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -171,6 +171,7 @@ set(DCA_LIBS
   signals
   coarsegraining
   ${DCA_CONCURRENCY_LIB}
+  parallel_no_concurrency
   ${DCA_THREADING_LIBS}
   lapack
   models
@@ -195,12 +196,12 @@ if (DCA_HAVE_GPU)
 )  
 endif()
 
-if (DCA_WITH_ADIOS2 AND DCA_HAVE_ADIOS2)
-  list(APPEND DCA_LIBS
-    dca_adios2 adios2::adios2
-    )
-  message("-- Add dca_adios2 to targets")
-endif()
+# if (DCA_WITH_ADIOS2 AND DCA_HAVE_ADIOS2)
+#   list(APPEND DCA_LIBS
+#     dca_adios2 adios2::adios2
+#     )
+#   message("-- Add dca_adios2 to targets")
+# endif()
 
 # The BLAS and LAPACK libraries in DCA_EXTERNAL_LIBS should be linked after MAGMA.
 list(APPEND DCA_LIBS lapack ${DCA_EXTERNAL_LIBS})
diff --git a/applications/analysis/chi_q_omega.cpp b/applications/analysis/chi_q_omega.cpp
index dbd655783..8c31370d9 100644
--- a/applications/analysis/chi_q_omega.cpp
+++ b/applications/analysis/chi_q_omega.cpp
@@ -77,7 +77,7 @@ int main(int argc, char** argv) {
   if (dca::io::stringToIOType(parameters.get_output_format()) == dca::io::IOType::ADIOS2) {
     int rank = concurrency.id();
     std::cout << "\nProcessor " << concurrency.id() << " is writing data." << std::endl;
-    dca::io::Writer writer(adios, concurrency, parameters.get_output_format(), true);
+    dca::io::Writer writer(concurrency, parameters.get_output_format(), true);
     std::string filename_bse(parameters.get_directory() + parameters.getAppropriateFilenameAnalysis());
     writer.open_file(filename_bse);
 
diff --git a/applications/analysis/main_analysis.cpp b/applications/analysis/main_analysis.cpp
index e3cbf939a..617ff415c 100644
--- a/applications/analysis/main_analysis.cpp
+++ b/applications/analysis/main_analysis.cpp
@@ -60,16 +60,13 @@ int main(int argc, char** argv) {
   // Create and initialize the DCA data object and read the output of the DCA(+) calculation.
   DcaDataType dca_data(parameters);
   dca_data.initialize();
-#ifdef DCA_HAVE_ADIOS2
-  adios2::ADIOS adios;
-
   if (dca::io::stringToIOType(parameters.get_output_format()) == dca::io::IOType::ADIOS2) {
     std::cout << "\nProcessor " << concurrency.id() << " is writing data." << std::endl;
-    dca::io::Writer writer(adios, concurrency, parameters.get_output_format(), true);
+    dca::io::Writer writer(concurrency, parameters.get_output_format(), true);
     std::string filename_bse(parameters.get_directory() + parameters.getAppropriateFilenameAnalysis());
     writer.open_file(filename_bse);
 
-    dca_data.read(adios, parameters.get_directory() + parameters.get_filename_dca());
+    dca_data.read(parameters.get_directory() + parameters.get_filename_dca());
     BseSolverType bse_solver(parameters, dca_data);
     bse_solver.calculateSusceptibilities();
 
@@ -82,7 +79,6 @@ int main(int argc, char** argv) {
     }
   }
   else
-#endif
   {
     dca_data.read(parameters.get_directory() + parameters.get_filename_dca());
     BseSolverType bse_solver(parameters, dca_data);
diff --git a/applications/dca/CMakeLists.txt b/applications/dca/CMakeLists.txt
index 6bc3d37e6..aac9a9cf6 100644
--- a/applications/dca/CMakeLists.txt
+++ b/applications/dca/CMakeLists.txt
@@ -5,8 +5,8 @@ if (DCA_BUILD_DCA)
   target_include_directories(main_dca PRIVATE ${DCA_INCLUDE_DIRS})
 
   if (DCA_HAVE_GPU)
-    target_link_libraries(main_dca PRIVATE ${DCA_GPU_LIBS} g0_interpolation ${DCA_KERNEL_LIBS})
+    target_link_libraries(main_dca PRIVATE ${DCA_KERNEL_LIBS})
   endif()
   
-  target_link_libraries(main_dca PUBLIC FFTW::Double signals  ${DCA_LIBS})
+  target_link_libraries(main_dca PUBLIC FFTW::Double signals ${DCA_LIBS} dca_io)
 endif()
diff --git a/applications/dca/main_dca.cpp b/applications/dca/main_dca.cpp
index dcb5fbe23..835561e18 100644
--- a/applications/dca/main_dca.cpp
+++ b/applications/dca/main_dca.cpp
@@ -17,6 +17,8 @@
 #include "dca/config/dca.hpp"
 #include "dca/application/dca_loop_dispatch.hpp"
 #include "dca/config/cmake_options.hpp"
+#include "dca/config/haves_defines.hpp"
+
 // Defines Concurrency, Threading, ParametersType, DcaData, DcaLoop, and Profiler.
 #include "dca/io/json/json_reader.hpp"
 #include "dca/util/git_version.hpp"
diff --git a/build-aux/frontier_rocm6_build.sh b/build-aux/frontier_rocm6_build.sh
new file mode 100644
index 000000000..cf68705ed
--- /dev/null
+++ b/build-aux/frontier_rocm6_build.sh
@@ -0,0 +1,40 @@
+
+cmake -DDCA_WITH_CUDA=off -DDCA_WITH_HIP=ON \
+      -DFFTW_ROOT=$FFTW_PATH \
+      -DDCA_FIX_BROKEN_MPICH=ON \
+      -DROCM_ROOT=${ROCM_PATH} \
+      -DMAGMA_ROOT=${MAGMA_ROOT} \
+      -DLAPACK_ROOT=${OPENBLAS_ROOT} \
+      -DBLAS_ROOT=${OPENBLAS_ROOT} \
+      -DDCA_WITH_TESTS_FAST=ON \
+      -DTEST_RUNNER="srun" \
+      -DGPU_TARGETS=gfx90a \
+      -DAMDGPU_TARGETS=gfx90a \
+      -DCMAKE_C_COMPILER=mpicc \
+      -DCMAKE_CXX_COMPILER=mpic++ \
+      -DCMAKE_HIP_COMPILER=/opt/rocm-6.0.0/llvm/bin/clang++ \
+      -DCMAKE_INSTALL_PREFIX=$INST \
+      -DCMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}" \
+      -GNinja \
+      ..
+
+#cmake -DDCA_WITH_CUDA=off -DDCA_WITH_HIP=ON \
+      -DFFTW_ROOT=$FFTW_PATH \
+      -DDCA_FIX_BROKEN_MPICH=ON \
+      -DROCM_ROOT=${ROCM_PATH} \
+      -DMAGMA_ROOT=${MAGMA_ROOT} \
+      -DLAPACK_ROOT=${OPENBLAS_ROOT} \
+      -DBLAS_ROOT=${OPENBLAS_ROOT} \
+      -DDCA_WITH_TESTS_FAST=ON \
+      -DTEST_RUNNER="srun" \
+      -DGPU_TARGETS=gfx90a \
+      -DAMDGPU_TARGETS=gfx90a \
+      -DCMAKE_C_COMPILER=mpicc \
+      -DCMAKE_CXX_COMPILER=mpic++ \
+      -DCMAKE_HIP_COMPILER=/opt/rocm-6.0.0/llvm/bin/clang++ \
+      -DCMAKE_INSTALL_PREFIX=$INST \
+      -DCMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}" \
+      -GNinja \
+      ..
+# cmake -DDCA_WITH_CUDA=off -DDCA_WITH_HIP=ON -DFFTW_ROOT=$FFTW_PATH -DDCA_FIX_BROKEN_MPICH=ON -DROCM_ROOT=${ROCM_PATH} -DMAGMA_ROOT=${MAGMA_ROOT} -DLAPACK_ROOT=${OPENBLAS_ROOT} -DBLAS_ROOT=${OPENBLAS_ROOT} -DDCA_WITH_TESTS_FAST=ON -DTEST_RUNNER="srun" -DGPU_TARGETS=gfx90a -DAMDGPU_TARGETS=gfx90a -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpic++ -DCMAKE_HIP_COMPILER=/opt/rocm-6.0.0/llvm/bin/clang++ -DCMAKE_INSTALL_PREFIX=$INST -DCMAKE_PREFIX_PATH="${CMAKE_PREFIX_PATH}" -GNinja ..
+      ..
diff --git a/build-aux/frontier_rocm6_just_exports.sh b/build-aux/frontier_rocm6_just_exports.sh
new file mode 100644
index 000000000..059c296f3
--- /dev/null
+++ b/build-aux/frontier_rocm6_just_exports.sh
@@ -0,0 +1,7 @@
+export CC=mpicc
+export CXX=mpicxx
+
+export OPENBLAS_ROOT=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/gcc-11.2.0/openblas-0.3.25-scaywvuh5zsm5u7smg54plj2oyf7nekv
+export HDF5_ROOT=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/rocmcc-6.0.0/hdf5-1.12.1-ajskwiaabdvgc36ozb6hzqnrwu2becha
+export MAGMA_ROOT=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/rocmcc-6.0.0/magma-master-rizw3ajkhfcq5cjutoykgkkv5hexftoz
+export FFTW_PATH=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/rocmcc-6.0.0/fftw-3.3.10-2mykijticsr5rfbyunax4zrwhhzcb7qm
diff --git a/build-aux/frontier_rocm6_load_modules.sh b/build-aux/frontier_rocm6_load_modules.sh
new file mode 100644
index 000000000..44951217c
--- /dev/null
+++ b/build-aux/frontier_rocm6_load_modules.sh
@@ -0,0 +1,26 @@
+#!/bin/bash
+#
+# Loads all modules that are required to build DCA++ on ORNL's Frontier.
+# A reset is done at the beginning to restore to the default programming environment on Frontier.
+# This is for development only at this point.
+#
+# Usage: source frontier_load_modules.sh
+
+
+module reset
+module load amd-mixed/6.0.0
+spack load cmake%gcc@11.2.0
+spack load ninja%gcc@11.2.0
+spack load magma@master amdgpu_target=gfx90a
+spack load hdf5@1.12.1 +cxx ~mpi api=v112 %rocmcc@6.0.0
+spack load fftw ~mpi %rocmcc@6.0.0
+spack load openblas@0.3.25 %gcc@11.2.0
+
+export CC=mpicc
+export CXX=mpicxx
+
+export OPENBLAS_ROOT=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/gcc-11.2.0/openblas-0.3.25-scaywvuh5zsm5u7smg54plj2oyf7nekv
+export HDF5_ROOT=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/rocmcc-6.0.0/hdf5-1.12.1-ajskwiaabdvgc36ozb6hzqnrwu2becha
+export MAGMA_ROOT=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/rocmcc-6.0.0/magma-master-rizw3ajkhfcq5cjutoykgkkv5hexftoz
+export FFTW_PATH=/lustre/orion/cph102/proj-shared/epd/spack/opt/spack/linux-sles15-zen3/rocmcc-6.0.0/fftw-3.3.10-2mykijticsr5rfbyunax4zrwhhzcb7qm
+#export LD_PRELOAD=/opt/cray/pe/lib64/cce/libtcmalloc_minimal.so.1
diff --git a/cmake/dca_cuda.cmake b/cmake/dca_cuda.cmake
index 54a2295b2..746eace27 100644
--- a/cmake/dca_cuda.cmake
+++ b/cmake/dca_cuda.cmake
@@ -31,7 +31,7 @@ if (CMAKE_CUDA_COMPILER)
   dca_add_haves_define(DCA_HAVE_CUDA)
   dca_add_haves_define(DCA_HAVE_GPU)
 
-  list(APPEND DCA_GPU_LIBS CUDA::cudart CUDA::cublas)
+  list(APPEND DCA_GPU_LIBS CUDA::cudart CUDA::cublas CUDA::cusparse)
   set(DCA_CUDA_PROPERTIES "CMAKE_CUDA_ARCHITECTURES 70")
   list(APPEND CUDAFLAGS "--expt-relaxed-constexpr" ${DCA_CUDA_OPTIONS})
   set(CMAKE_CUDA_STANDARD 17)
diff --git a/cmake/dca_external_libs.cmake b/cmake/dca_external_libs.cmake
index f2b01f519..29af1285f 100644
--- a/cmake/dca_external_libs.cmake
+++ b/cmake/dca_external_libs.cmake
@@ -41,7 +41,7 @@ list(APPEND DCA_EXTERNAL_LIBS ${BLAS_LIBRARIES})
 # HDF5
 
 if (NOT HDF5_LIBRARIES)
-set(HDF5_NO_FIND_PACKAGE_CONFIG_FILE true)
+set(HDF5_NO_FIND_PACKAGE_CONFIG_FILE false)
 set(HDF5_PREFER_PARALLEL false)
 find_package(HDF5 REQUIRED COMPONENTS C CXX)
 message("HDF5: ${HDF5_FOUND} ${HDF5_LIBRARIES}")
diff --git a/cmake/dca_hip.cmake b/cmake/dca_hip.cmake
index bc19fc9ec..2d3e0ebb3 100644
--- a/cmake/dca_hip.cmake
+++ b/cmake/dca_hip.cmake
@@ -39,14 +39,19 @@ if(DCA_WITH_HIP)
 #-------------------------------------------------------------------
 #  set up HIP compiler options
 #-------------------------------------------------------------------
-  set(CMAKE_MODULE_PATH "${ROCM_ROOT}/hip/cmake" ${CMAKE_MODULE_PATH})
+  set(CMAKE_MODULE_PATH "${ROCM_ROOT}/hip/cmake" "${ROCM_ROOT}/lib/cmake/hip" "${ROCM_ROOT}/lib/cmake/hipblas" "${ROCM_ROOT}/lib/cmake/rocthrust" ${CMAKE_MODULE_PATH})
   find_package(HIP REQUIRED)
   find_package(hipblas REQUIRED)
   find_package(hipsparse REQUIRED)
   find_package(rocsolver REQUIRED)
+  find_package(rocthrust REQUIRED)
 
 endif(DCA_WITH_HIP)
 
+get_property(hipblas_include_dirs TARGET roc::hipblas PROPERTY INTERFACE_INCLUDE_DIRECTORIES)
+message("hipblas includes: ${hipblas_include_dirs}")
+
+
 #set(CUDA_ARCHITECTURES "sm_60" CACHE STRING "Name of the real architecture to build for.")
 set(MAGMA_ROOT "" CACHE PATH "Path to the MAGMA installation directory. Hint for CMake to find MAGMA.")
 
@@ -66,11 +71,12 @@ if (CMAKE_HIP_COMPILER)
   dca_add_haves_define(DCA_HAVE_GPU)
   dca_add_haves_define(__HIP_PLATFORM_AMD__)
   list(APPEND DCA_GPU_LIBS hip::host roc::hipblas roc::hipsparse)
-  set(DCA_HIP_PROPERTIES "CMAKE_HIP_ARCHITECTURES gfx906,gfx908")
+  set(DCA_HIP_PROPERTIES "CMAKE_HIP_ARCHITECTURES gfx908,gfx90a")
   set(CMAKE_HIP_STANDARD 17)
   list(APPEND HIP_HIPCC_FLAGS "-fPIC")
   # doesn't appear to work
   set(CMAKE_HIP_SOURCE_FILE_EXTENSIONS cu)
+  message("Enabled HIP as a language")
   # NOTE: this is solved by dca_linking.cmake: dca_gpu_device_link()
   # alternative method (same issue)
   #file(GLOB_RECURSE CUDA_KERNELS_SRC ${PROJECT_SOURCE_DIR} *.cu)
@@ -106,6 +112,7 @@ if (MAGMA_LIBRARY AND MAGMA_INCLUDE_DIR)
   # I have built magma without openmp for
   # CI. But if you naively use a random systems
   # magma expect to have a link error.
-  target_link_libraries(magma::sparse INTERFACE magma::magma roc::hipblas roc::hipsparse)
-  list(APPEND DCA_GPU_LIBS ${MAGMA_LIBRARY} roc::hipsparse)
+  target_link_libraries(magma::magma INTERFACE roc::hipblas roc::hipsparse LAPACK::LAPACK BLAS::BLAS)
+  target_link_libraries(magma::sparse INTERFACE magma::magma)
+  list(APPEND DCA_GPU_LIBS ${MAGMA_LIBRARY} roc::hipsparse roc::hipblas)
 endif()
diff --git a/cmake/dca_linking.cmake b/cmake/dca_linking.cmake
index 040e0bada..9fd39f2b7 100644
--- a/cmake/dca_linking.cmake
+++ b/cmake/dca_linking.cmake
@@ -2,7 +2,7 @@
 #link the correct gpu runtime library
 function(dca_gpu_runtime_link target_name)
   if(DCA_HAVE_HIP)
-    target_link_libraries(${target_name} PUBLIC hip::host)
+    target_link_libraries(${target_name} PUBLIC hip::host roc::hipblas roc::hipsparse)
     message("linking target ${target_name} to hip::host")
   elseif(DCA_HAVE_CUDA)
     target_link_libraries(${target_name} PUBLIC CUDA::cudart)
@@ -12,7 +12,7 @@ endfunction()
 #link the correct gpu runtime library
 function(dca_gpu_blas_link target_name)
   if(DCA_HAVE_HIP)
-    target_link_libraries(${target_name} PUBLIC roc::hipblas)
+    target_link_libraries(${target_name} PUBLIC roc::hipblas roc::hipsparse)
     message("linking target ${target_name} to roc::hipblas")
   elseif(DCA_HAVE_CUDA)
     target_link_libraries(${target_name} PUBLIC CUDA::cublas)
@@ -26,7 +26,7 @@ function(dca_gpu_device_link target_name)
       PROPERTIES HIP_SEPARABLE_COMPILATION ON)
     set_target_properties( ${target_name}
       PROPERTIES HIP_RESOLVE_DEVICE_SYMBOLS ON)
-    target_link_libraries(${target_name} PRIVATE hip::device)
+    target_link_libraries(${target_name} PRIVATE hip::device roc::hipblas roc::hipsparse roc::rocthrust)
     get_target_property(_srcs ${target_name} SOURCES)
     get_target_property(_src_dir ${target_name} SOURCE_DIR)
     #
diff --git a/cmake/dca_testing.cmake b/cmake/dca_testing.cmake
index 9fa851207..cc597e4e4 100644
--- a/cmake/dca_testing.cmake
+++ b/cmake/dca_testing.cmake
@@ -34,7 +34,7 @@ endif()
 # MPI or CUDA may be given to indicate that the test requires these libraries. MPI_NUMPROC is the
 # number of MPI processes to use for a test with MPI, the default value is 1.
 function(dca_add_gtest name)
-  set(options FAST EXTENSIVE STOCHASTIC PERFORMANCE GTEST_MAIN THREADED MPI CUDA CUDA_MPI)
+  set(options FAST EXTENSIVE STOCHASTIC PERFORMANCE GTEST_MAIN GTEST_MPI_MAIN THREADED MPI CUDA CUDA_MPI)
   set(oneValueArgs MPI_NUMPROC)
   set(multiValueArgs INCLUDE_DIRS SOURCES LIBS)
   cmake_parse_arguments(DCA_ADD_GTEST "${options}" "${oneValueArgs}" "${multiValueArgs}" ${ARGN})
@@ -108,6 +108,10 @@ function(dca_add_gtest name)
     set(DCA_ADD_GTEST_SOURCES ${PROJECT_SOURCE_DIR}/test/dca_gtest_main.cpp ${DCA_ADD_GTEST_SOURCES})
   endif()
 
+  if (DCA_ADD_GTEST_GTEST_MPI_MAIN)
+    set(DCA_ADD_GTEST_SOURCES ${PROJECT_SOURCE_DIR}/test/dca_gtest_main_mpi.cpp ${DCA_ADD_GTEST_SOURCES})
+  endif()
+
   add_executable(${name} ${name}.cpp ${DCA_ADD_GTEST_SOURCES})
 
   target_compile_definitions(${name} PRIVATE DCA_SOURCE_DIR=\"${PROJECT_SOURCE_DIR}\")
diff --git a/include/dca/function/function.hpp b/include/dca/function/function.hpp
index 487786a3a..359402f1e 100644
--- a/include/dca/function/function.hpp
+++ b/include/dca/function/function.hpp
@@ -436,7 +436,7 @@ function<scalartype, domain, DT>::function(const std::string& name)
       std::cerr << "large functions need names give yourself a chance.\n";
   }
   // will zero real or complex values
-  fnc_values_.resize(dmn.get_size(), {});
+  fnc_values_.resize(dmn.get_size());
 }
 
 /** copy constructor
@@ -512,7 +512,7 @@ function<scalartype, domain, DT>::function(std::initializer_list<scalartype> ini
       Nb_sbdms(dmn.get_leaf_domain_sizes().size()) {
   start_ = 0;
   end_ = dmn.get_size();
-  fnc_values_.resize(dmn.get_size(), {});
+  fnc_values_.resize(dmn.get_size());
   std::copy_n(init_list.begin(), init_list.size(), fnc_values_.begin());
 }
 
@@ -657,7 +657,7 @@ function<scalartype, domain, DT>& function<scalartype, domain, DT>::operator=(
       Nb_sbdms = other.dmn.get_leaf_domain_sizes().size();
       start_ = other.start_;
       end_ = other.end_;
-      fnc_values_.resize(other.size(), {});
+      fnc_values_.resize(other.size());
     }
     fnc_values_ = other.fnc_values_;
   }
@@ -679,7 +679,7 @@ inline function<Scalar, domain, DT>& function<Scalar, domain, DT>::operator=(
         Nb_sbdms = other.dmn.get_leaf_domain_sizes().size();
         start_ = other.start_;
         end_ = other.end_;
-        fnc_values_.resize(other.size(), {});
+        fnc_values_.resize(other.size());
       }
       fnc_values_ = other.fnc_values_;
     }
@@ -694,7 +694,7 @@ inline function<Scalar, domain, DT>& function<Scalar, domain, DT>::operator=(
       Nb_sbdms = other.dmn.get_leaf_domain_sizes().size();
       start_ = other.start_;
       end_ = other.end_;
-      fnc_values_.resize(other.size(), {});
+      fnc_values_.resize(other.size());
     }
     auto kConvert = [](auto& kvec) -> std::vector<Scalar> {
       std::vector<Scalar> k_converted(kvec.size());
diff --git a/include/dca/io/adios2/adios2_global.hpp b/include/dca/io/adios2/adios2_global.hpp
index 0e2772596..6f06144fc 100644
--- a/include/dca/io/adios2/adios2_global.hpp
+++ b/include/dca/io/adios2/adios2_global.hpp
@@ -10,7 +10,7 @@
 // ADIOS2 global object
 
 
-/** For testing only, in main_dca adios2::ADIOS is a member of dca_loop
+/** For testing only, in main_dca adios2::ADIOS is a member of concurrency the concurrency context owned by main.
  */
 class GlobalAdios
 {
diff --git a/include/dca/io/adios2/adios2_reader.hpp b/include/dca/io/adios2/adios2_reader.hpp
index b6b7e63fa..8e5e08a96 100644
--- a/include/dca/io/adios2/adios2_reader.hpp
+++ b/include/dca/io/adios2/adios2_reader.hpp
@@ -43,12 +43,8 @@ class ADIOS2Reader {
   typedef adios2::ADIOS file_type;
 
 public:
-  // In: verbose. If true, the reader outputs a short log whenever it is executed.
-  ADIOS2Reader(adios2::ADIOS& adios, const CT* concurrency,
-               bool verbose = false);
-
-  /** maybe this one uses a singleton adios2::ADIOS accessor */
-  ADIOS2Reader(const CT* concurrency,
+  ADIOS2Reader() = delete;
+  ADIOS2Reader(CT& concurrency,
                bool verbose = false);
 
   ~ADIOS2Reader();
@@ -150,9 +146,10 @@ class ADIOS2Reader {
   template <class T>
   std::string VectorToString(const std::vector<T>& v);
 
+  // concurrency_ must be before the adios reference
+  const CT& concurrency_;
   adios2::ADIOS& adios_;
   const bool verbose_;
-  const CT* concurrency_;
 
   adios2::IO io_;
   std::string io_name_;
@@ -394,7 +391,7 @@ bool ADIOS2Reader<CT>::execute(const std::string& name, func::function<Scalartyp
 
   if (verbose_) {
     std::cout << "\t ADIOS2Reader: Read function : " << name
-              << " in linear distributed manner, rank = " << concurrency_->id()
+              << " in linear distributed manner, rank = " << concurrency_.id()
               << " start = " << std::to_string(start) << " end = " << std::to_string(end) << "\n";
   }
 
@@ -428,7 +425,7 @@ bool ADIOS2Reader<CT>::execute(const std::string& name, func::function<Scalartyp
 
     if (verbose_) {
       std::cout << "\t ADIOS2Reader: Read function : " << name
-                << " in linear distributed manner, rank = " << concurrency_->id()
+                << " in linear distributed manner, rank = " << concurrency_.id()
                 << " shape = " << VectorToString(var.Shape()) << " start = " << VectorToString(s)
                 << " count = " << VectorToString(c) << "\n";
     }
@@ -476,7 +473,7 @@ bool ADIOS2Reader<CT>::execute(const std::string& name, func::function<Scalartyp
 
   if (verbose_) {
     std::cout << "\t ADIOS2Reader: Read function : " << name
-              << " in distributed manner, rank = " << concurrency_->id()
+              << " in distributed manner, rank = " << concurrency_.id()
               << " start = " << VectorToString(start) << " end = " << VectorToString(end) << "\n";
   }
 
@@ -521,7 +518,7 @@ bool ADIOS2Reader<CT>::execute(const std::string& name, func::function<Scalartyp
 
     if (verbose_) {
       std::cout << "\t ADIOS2Reader: Read function : " << name
-                << " in distributed manner, rank = " << concurrency_->id()
+                << " in distributed manner, rank = " << concurrency_.id()
                 << " shape = " << VectorToString(var.Shape()) << " start = " << VectorToString(s)
                 << " count = " << VectorToString(c) << "\n";
     }
@@ -614,10 +611,14 @@ std::string ADIOS2Reader<CT>::VectorToString(const std::vector<T>& v) {
 }
 
 extern template class ADIOS2Reader<dca::parallel::NoConcurrency>;
+extern template class ADIOS2Reader<dca::parallel::NoConcurrency const>;
 #ifdef DCA_HAVE_MPI
 extern template class ADIOS2Reader<dca::parallel::MPIConcurrency>;
+//this has become required for some reason. Const correctness for the ADIOS2Reader has gone wrong
+extern template class ADIOS2Reader<dca::parallel::MPIConcurrency const>;
 #endif
 }  // namespace io
 }  // namespace dca
 
 #endif  // DCA_IO_ADIOS2_ADIOS2_READER_HPP
+
diff --git a/include/dca/io/adios2/adios2_writer.hpp b/include/dca/io/adios2/adios2_writer.hpp
index 9c3f9ed77..e1ff3481a 100644
--- a/include/dca/io/adios2/adios2_writer.hpp
+++ b/include/dca/io/adios2/adios2_writer.hpp
@@ -48,6 +48,8 @@ class ADIOS2Writer {
 public:
   ADIOS2Writer() = delete;
   // In: verbose. If true, the writer outputs a short log whenever it is executed.
+  // Concurrency objects now carry an adios pointer if DCA_HAVE_ADIOS2.
+  ADIOS2Writer(const CT* concurrency, bool verbose = false);
   ADIOS2Writer(adios2::ADIOS& adios, const CT* concurrency, bool verbose = false);
   ~ADIOS2Writer();
 
@@ -171,9 +173,9 @@ class ADIOS2Writer {
   }
 
 private:
+  const CT* concurrency_;
   adios2::ADIOS& adios_;
   bool verbose_;
-  const CT* concurrency_;
 
   /** For the case of sample it is necessary to explicitly write a rank "local" variable
    *  of this type.
diff --git a/include/dca/io/hdf5/hdf5_reader.hpp b/include/dca/io/hdf5/hdf5_reader.hpp
index d8147f351..efe60c060 100644
--- a/include/dca/io/hdf5/hdf5_reader.hpp
+++ b/include/dca/io/hdf5/hdf5_reader.hpp
@@ -185,7 +185,7 @@ bool HDF5Reader::execute(const std::string& name, std::vector<std::vector<Scalar
 
   std::vector<hvl_t> data(size);
 
-  H5::DataSet dataset = file_->openDataSet(name.c_str());
+  H5::DataSet dataset = file_->openDataSet(full_name.c_str());
   dataset.read(data.data(), type);
 
   value.resize(size);
diff --git a/include/dca/io/reader.hpp b/include/dca/io/reader.hpp
index ba9251641..930cb4f6d 100644
--- a/include/dca/io/reader.hpp
+++ b/include/dca/io/reader.hpp
@@ -34,18 +34,20 @@ namespace dca::io {
 template <class Concurrency>
 class Reader {
 public:
-  using DCAReaderVariant = std::variant<io::HDF5Reader, io::JSONReader
+  using DCAReaderVariant = std::variant<std::monostate,
 #ifdef DCA_HAVE_ADIOS2
-                                        ,
-                                        io::ADIOS2Reader<Concurrency>
+                                        io::ADIOS2Reader<Concurrency>,
 #endif
-                                        >;
+                                        io::HDF5Reader, io::JSONReader>;
   /**
    * \param[in] concureency   reference to current concurrency env
    * \param[in] format        format as IOType
    * \param[in] verbose       If true, reader does some logging
    */
-  Reader(const Concurrency& concurrency, const IOType format, bool verbose = true)
+  Reader(Concurrency& concurrency, const std::string& format_tag, bool verbose = true)
+    : Reader(concurrency, stringToIOType(format_tag), verbose) {}
+  
+  Reader(Concurrency& concurrency, IOType format, bool verbose = true)
       : concurrency_(concurrency) {
     switch (format) {
       case IOType::HDF5:
@@ -56,81 +58,126 @@ class Reader {
         break;
       case IOType::ADIOS2:
 #ifdef DCA_HAVE_ADIOS2
-        reader_.template emplace<io::ADIOS2Reader<Concurrency>>(&concurrency, verbose);
+        reader_.template emplace<io::ADIOS2Reader<Concurrency>>(concurrency, verbose);
 
 #endif
         break;
     }
   }
 
-  /** DEPRECATED -- Support for format type as string
-   * \param[in] format        string representation of format, since parameters still use string
-   *
-   * \todo remove need for this constructor store IO format as enum class type.
-   */
-  Reader(const Concurrency& concurrency, const std::string& format, bool verbose = true)
-      : Reader(concurrency, stringToIOType(format), verbose) {}
-
-#ifdef DCA_HAVE_ADIOS2
-  Reader(adios2::ADIOS& adios, const Concurrency& concurrency, const std::string& format,
-         bool verbose = true)
-      : concurrency_(concurrency) {
-    if (format == "HDF5") {
-      throw(std::logic_error("ADIOS2 reference not an argument for hdf5 reader"));
-    }
-    else if (format == "JSON") {
-      throw(std::logic_error("ADIOS2 reference not an argument for json reader"));
-    }
-    else if (format == "ADIOS2") {
-      reader_.template emplace<io::ADIOS2Reader<Concurrency>>(adios, &concurrency, verbose);
-    }
-    else {
-      throw(std::logic_error("Invalid input format"));
-    }
-  }
-#endif
-
+  
   constexpr static bool is_reader = true;
   constexpr static bool is_writer = false;
 
   void open_file(const std::string& file_name) {
-    std::visit([&](auto& var) { var.open_file(file_name); }, reader_);
+    std::visit(
+        [&](auto& var) {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            var.open_file(file_name);
+        },
+        reader_);
   }
   void close_file() {
-    std::visit([&](auto& var) { var.close_file(); }, reader_);
+    std::visit(
+        [&](auto& var) {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            var.close_file();
+        },
+        reader_);
   }
 
   /** For reading input there is great utility in knowing if a group is present.
    *  It isn't an exceptional circumstance if a group is not present.
    */
   bool open_group(const std::string& new_path) {
-    return std::visit([&](auto& var) -> bool { return var.open_group(new_path); }, reader_);
+    return std::visit(
+        [&](auto& var) -> bool {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            return var.open_group(new_path);
+        },
+        reader_);
   }
 
   void close_group() {
-    std::visit([&](auto& var) { var.close_group(); }, reader_);
+    std::visit(
+        [&](auto& var) {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            var.close_group();
+        },
+        reader_);
   }
 
   long getStepCount() {
-    return std::visit([&](auto& var) ->std::size_t { return var.getStepCount(); }, reader_);
+    return std::visit(
+        [&](auto& var) -> std::size_t {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            return var.getStepCount();
+        },
+        reader_);
   }
 
   void begin_step() {
-    std::visit([&](auto& var) { var.begin_step(); }, reader_);
+    std::visit(
+        [&](auto& var) {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            var.begin_step();
+        },
+        reader_);
   }
 
   void end_step() {
-    std::visit([&](auto& var) { var.end_step(); }, reader_);
+    std::visit(
+        [&](auto& var) {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            var.end_step();
+        },
+        reader_);
   }
 
   std::string get_path() {
-    return std::visit([&](auto& var) -> std::string { return var.get_path(); }, reader_);
+    return std::visit(
+        [&](auto& var) -> std::string {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            return var.get_path();
+        },
+        reader_);
   }
 
   template <class... Args>
   bool execute(Args&&... args) noexcept {
-    return std::visit([&](auto& var) -> bool { return var.execute(std::forward<Args>(args)...); },
-                      reader_);
+    return std::visit(
+        [&](auto& var) -> bool {
+          if constexpr (std::is_same_v<std::monostate&, decltype(var)>)
+            throw std::runtime_error(
+                "No operations should ever occur on monostate in reader variant");
+          else
+            return var.execute(std::forward<Args>(args)...);
+        },
+        reader_);
   }
 
   DCAReaderVariant& getUnderlying() {
@@ -139,7 +186,7 @@ class Reader {
 
 private:
   DCAReaderVariant reader_;
-  const Concurrency& concurrency_;
+  Concurrency& concurrency_;
 };
 
 extern template class Reader<dca::parallel::NoConcurrency>;
diff --git a/include/dca/io/writer.hpp b/include/dca/io/writer.hpp
index ec79f8484..3a1767f06 100644
--- a/include/dca/io/writer.hpp
+++ b/include/dca/io/writer.hpp
@@ -38,7 +38,6 @@ class Writer {
 #endif
                                         >;
   /** Constructor for writer, pretty tortured due to optional Adios2.
-   *  \param [in] adios        if build with ADIOS2 support adios2 env reference.
    *  \param [in] concurrency  reference to the concurrency environment
    *  \param [in] format       output format: ADIOS2, HDF5 or JSON.
    *  \param [in] verbose      If true, the writer outputs a short log whenever it is executed.
@@ -47,13 +46,10 @@ class Writer {
    *  used by the writer.
    */
   Writer(
-#ifdef DCA_HAVE_ADIOS2
-      adios2::ADIOS& adios,
-#endif
       Concurrency& concurrency, const std::string& format, bool verbose = true)
       :
 #ifdef DCA_HAVE_ADIOS2
-        adios_(adios),
+    adios_(concurrency.get_adios()),
 #endif
         concurrency_(concurrency) {
     if (format == "HDF5") {
@@ -64,7 +60,7 @@ class Writer {
     }
 #ifdef DCA_HAVE_ADIOS2
     else if (format == "ADIOS2") {
-      writer_.template emplace<io::ADIOS2Writer<Concurrency>>(adios_, &concurrency_, verbose);
+      writer_.template emplace<io::ADIOS2Writer<Concurrency>>(&concurrency_, verbose);
     }
 #endif
     else {
diff --git a/include/dca/linalg/matrix.hpp b/include/dca/linalg/matrix.hpp
index 138eff146..e55181138 100644
--- a/include/dca/linalg/matrix.hpp
+++ b/include/dca/linalg/matrix.hpp
@@ -306,7 +306,7 @@ Matrix<ScalarType, device_name,  ALLOC>::Matrix(const Matrix<ScalarRhs, rhs_devi
     : name_(name), size_(rhs.size_), capacity_(rhs.capacity_) {
   if (sizeof(ScalarType) != sizeof(ScalarRhs))
     throw std::runtime_error("conversion of both type and location of Matrix not currently possible!");
-  data_ = Allocator::allocate(nrElements(capacity_));
+  data_ = ALLOC::allocate(nrElements(capacity_));
   util::memoryCopy(data_, leadingDimension(), rhs.data_, rhs.leadingDimension(), size_);
 }
 
@@ -319,7 +319,7 @@ Matrix<ScalarType, device_name, ALLOC>::Matrix(const std::string& name, std::pai
   assert(capacity.first >= size_.first && capacity.second >= size_.second);
   assert(capacity_.first >= capacity.first && capacity_.second >= capacity.second);
 
-  data_ = Allocator::allocate(nrElements(capacity_));
+  data_ = ALLOC::allocate(nrElements(capacity_));
   util::Memory<device_name>::setToZero(data_, nrElements(capacity_));
 }
 
diff --git a/include/dca/linalg/matrixop.hpp b/include/dca/linalg/matrixop.hpp
index 2fb4a4feb..87c650d25 100644
--- a/include/dca/linalg/matrixop.hpp
+++ b/include/dca/linalg/matrixop.hpp
@@ -68,7 +68,8 @@ inline void copyMatrixToArray(const Matrix<Scalar, CPU, ALLOC>& mat, Scalar* a,
 // Copies the m by n matrix stored in a to the matrix mat.
 // Preconditions: lda >= m.
 template <typename Scalar, class ALLOC>
-inline void copyArrayToMatrix(int m, int n, const Scalar* a, int lda, Matrix<Scalar, CPU, ALLOC>& mat) {
+inline void copyArrayToMatrix(int m, int n, const Scalar* a, int lda,
+                              Matrix<Scalar, CPU, ALLOC>& mat) {
   assert(lda >= m);
   mat.resizeNoCopy(std::make_pair(m, n));
   lapack::lacpy("A", mat.nrRows(), mat.nrCols(), a, lda, mat.ptr(), mat.leadingDimension());
@@ -95,8 +96,9 @@ inline void copyCol(const Matrix<Scalar, device_name>& mat_x, int jx,
 //                0 <= j_x[i] < mat_x.nrCols() for 0 <= i < j_x.size(),
 //                0 <= j_y[i] < mat_y.nrCols() for 0 <= i < j_x.size().
 template <typename Scalar, class Vec, class ALLOC>
-inline void copyCols(const Matrix<Scalar, CPU, ALLOC>& mat_x, const Vec& j_x, Matrix<Scalar, CPU, ALLOC>& mat_y,
-                     const Vec& j_y, int /*thread_id*/ = 0, int /*stream_id*/ = 0) {
+inline void copyCols(const Matrix<Scalar, CPU, ALLOC>& mat_x, const Vec& j_x,
+                     Matrix<Scalar, CPU, ALLOC>& mat_y, const Vec& j_y, int /*thread_id*/ = 0,
+                     int /*stream_id*/ = 0) {
   assert(j_x.size() <= j_y.size());
 
   for (int ind_j = 0; ind_j < j_x.size(); ++ind_j)
@@ -127,7 +129,6 @@ inline void copyCols(const Matrix<Scalar, GPU>& mat_x, const Vector<int, GPU>& j
   blas::copyCols(mat_x.nrRows(), j_x.size(), j_x.ptr(), mat_x.ptr(), mat_x.leadingDimension(),
                  mat_y.ptr(), mat_y.leadingDimension(), thread_id, stream_id);
   checkErrorsCudaDebug();
-
 }
 #endif  // DCA_HAVE_GPU
 
@@ -153,8 +154,9 @@ inline void copyRow(const Matrix<Scalar, device_name>& mat_x, int ix,
 //                0 <= i_x[i] < mat_x.nrRows() for 0 <= i < i_x.size(),
 //                0 <= i_y[i] < mat_y.nrRows() for 0 <= i < i_x.size().
 template <typename Scalar, class Vec, class ALLOC>
-inline void copyRows(const Matrix<Scalar, CPU, ALLOC>& mat_x, const Vec& i_x, Matrix<Scalar, CPU, ALLOC>& mat_y,
-                     const Vec& i_y, const int /*thread_id*/ = 0, const int /*stream_id*/ = 0) {
+inline void copyRows(const Matrix<Scalar, CPU, ALLOC>& mat_x, const Vec& i_x,
+                     Matrix<Scalar, CPU, ALLOC>& mat_y, const Vec& i_y, const int /*thread_id*/ = 0,
+                     const int /*stream_id*/ = 0) {
   assert(i_x.size() <= i_y.size());
   assert(mat_x.nrCols() == mat_y.nrCols());
 
@@ -200,7 +202,6 @@ auto difference(const Matrix<Scalar, CPU, ALLOC>& a, const Matrix<Scalar, CPU, A
   auto max_diff = std::abs(Scalar(0));
   assert(a.size() == b.size());
 
-
   for (int j = 0; j < a.nrCols(); ++j) {
     for (int i = 0; i < a.nrRows(); ++i) {
       max_diff = std::max(max_diff, std::abs(a(i, j) - b(i, j)));
@@ -228,7 +229,7 @@ auto difference(const Matrix<Scalar, CPU, ALLOC>& a, const Matrix<Scalar, CPU, A
 
   return max_diff;
 }
-  template <typename Scalar, DeviceType device_name, class ALLOC>
+template <typename Scalar, DeviceType device_name, class ALLOC>
 auto difference(const Matrix<Scalar, device_name>& a, const Matrix<Scalar, CPU, ALLOC>& b,
                 double diff_threshold = 1e-3) {
   Matrix<Scalar, CPU, ALLOC> cp_a(a);
@@ -295,8 +296,9 @@ void insertRow(Matrix<Scalar, CPU, ALLOC>& mat, int i) {
 // Preconditions: mat is a square matrix.
 // Postconditions: ipiv and work are resized to the needed dimension.
 // \todo consider doing inverse at full precision reguardless of incoming Scalar precision
-  template <typename Scalar, DeviceType device_name, class ALLOC, template <typename, DeviceType, class> class MatrixType>
-  void inverse(MatrixType<Scalar, device_name, ALLOC>& mat, Vector<int, CPU>& ipiv,
+template <typename Scalar, DeviceType device_name, template <typename> class ALLOC,
+          template <typename, DeviceType, class> class MatrixType>
+void inverse(MatrixType<Scalar, device_name, ALLOC<Scalar>>& mat, Vector<int, CPU>& ipiv,
              Vector<Scalar, device_name>& work) {
   assert(mat.is_square());
 
@@ -312,15 +314,17 @@ void insertRow(Matrix<Scalar, CPU, ALLOC>& mat, int i) {
                                         work.ptr(), lwork);
 }
 
-  template <typename Scalar, DeviceType device_name, class ALLOC, template <typename, DeviceType, class> class MatrixType>
-  void inverse(MatrixType<Scalar, device_name, ALLOC>& mat) {
+  template <typename Scalar, DeviceType device_name, template<typename> class ALLOC,
+          template <typename, DeviceType, class> class MatrixType>
+void inverse(MatrixType<Scalar, device_name, ALLOC<Scalar>>& mat) {
   Vector<int, CPU> ipiv;
   Vector<Scalar, device_name> work;
   inverse(mat, ipiv, work);
 }
 
 template <typename Scalar, class ALLOC>
-void smallInverse(Matrix<Scalar, CPU, ALLOC>& m_inv, Vector<int, CPU>& ipiv, Vector<Scalar, CPU>& work) {
+void smallInverse(Matrix<Scalar, CPU, ALLOC>& m_inv, Vector<int, CPU>& ipiv,
+                  Vector<Scalar, CPU>& work) {
   assert(m_inv.is_square());
   switch (m_inv.nrCols()) {
     case 1:
@@ -631,8 +635,8 @@ inline void swapRowAndCol(Matrix<Scalar, device_name>& mat, int i1, int i2, int
 //                a.nrRows() == y.size() if transa == 'N', a.nrCols() == y.size() otherwise,
 //                a.nrCols() == x.size() if transa == 'N', a.nrRows() == x.size() otherwise.
 template <typename Scalar, class ALLOC>
-void gemv(char transa, Scalar alpha, const Matrix<Scalar, CPU, ALLOC>& a, const Vector<Scalar, CPU>& x,
-          Scalar beta, Vector<Scalar, CPU>& y) {
+void gemv(char transa, Scalar alpha, const Matrix<Scalar, CPU, ALLOC>& a,
+          const Vector<Scalar, CPU>& x, Scalar beta, Vector<Scalar, CPU>& y) {
   if (transa == 'N') {
     assert(a.nrRows() == y.size());
     assert(a.nrCols() == x.size());
@@ -669,11 +673,12 @@ void gemv(char transa, const Matrix<Scalar, CPU, ALLOC>& a, const Vector<Scalar,
 //                b.nrCols() == c.nrCols() if transb == 'N', b.nrRows() == c.nrCols() otherwise,
 //                ka == kb, where ka = a.nrCols() if transa == 'N', ka = a.nrRows() otherwise and
 //                          kb = b.nrRows() if transb == 'N', kb = b.nrCols() otherwise.
-  template <typename Scalar, DeviceType device_name, template <typename, DeviceType, class> class MatrixA,
-	    template <typename, DeviceType, class> class MatrixB, template <typename, DeviceType, class> class MatrixC, class ALLOC1, class ALLOC2, class ALLOC3>
-  void gemm(char transa, char transb, Scalar alpha, const MatrixA<Scalar, device_name, ALLOC1>& a,
-	    const MatrixB<Scalar, device_name, ALLOC2>& b, Scalar beta, MatrixC<Scalar, device_name, ALLOC3>& c,
-          int thread_id = 0, int stream_id = 0) {
+template <typename Scalar, DeviceType device_name, template <typename, DeviceType, class> class MatrixA,
+          template <typename, DeviceType, class> class MatrixB,
+          template <typename, DeviceType, class> class MatrixC, class ALLOC1, class ALLOC2, class ALLOC3>
+void gemm(char transa, char transb, Scalar alpha, const MatrixA<Scalar, device_name, ALLOC1>& a,
+          const MatrixB<Scalar, device_name, ALLOC2>& b, Scalar beta,
+          MatrixC<Scalar, device_name, ALLOC3>& c, int thread_id = 0, int stream_id = 0) {
   int m = c.nrRows();
   int n = c.nrCols();
   int k;
@@ -707,21 +712,26 @@ void gemv(char transa, const Matrix<Scalar, CPU, ALLOC>& a, const Vector<Scalar,
 // Performs the matrix-matrix multiplication c <- a * b
 // Out: c
 // Preconditions: a.nrRows() == c.nrRows(), b.nrCols() == c.nrCols() and a.nrCols() == b.nrRows()
-  template <typename Scalar, DeviceType device_name, class ALLOC, template <typename, DeviceType, class> class MatrixA,
-	    template <typename, DeviceType, class> class MatrixB, template <typename, DeviceType, class> class MatrixC>
-  inline void gemm(const MatrixA<Scalar, device_name, ALLOC>& a, const MatrixB<Scalar, device_name, ALLOC>& b,
-		   MatrixC<Scalar, device_name, ALLOC>& c, int thread_id = 0, int stream_id = 0) {
+template <typename Scalar, DeviceType device_name, class ALLOC,
+          template <typename, DeviceType, class> class MatrixA,
+          template <typename, DeviceType, class> class MatrixB,
+          template <typename, DeviceType, class> class MatrixC>
+inline void gemm(const MatrixA<Scalar, device_name, ALLOC>& a,
+                 const MatrixB<Scalar, device_name, ALLOC>& b,
+                 MatrixC<Scalar, device_name, ALLOC>& c, int thread_id = 0, int stream_id = 0) {
   gemm<Scalar, device_name>('N', 'N', 1., a, b, 0., c, thread_id, stream_id);
 }
 
 // Performs the matrix-matrix multiplication c <- alpha * a * b + beta * c,
 // In/Out: c ('In' only if beta != 0)
 // Preconditions: a.nrRows() == c.nrRows(), b.nrCols() == c.nrCols() and a.nrCols() == b.nrRows()
-  template <typename Scalar, DeviceType device_name, class ALLOC, template <typename, DeviceType, class> class MatrixA,
-	    template <typename, DeviceType,class> class MatrixB, template <typename, DeviceType, class> class MatrixC>
-  inline void gemm(Scalar alpha, const MatrixA<Scalar, device_name, ALLOC>& a,
-		   const MatrixB<Scalar, device_name, ALLOC>& b, Scalar beta,
-		   MatrixC<Scalar, device_name, ALLOC>& c, int thread_id = 0, int stream_id = 0) {
+template <typename Scalar, DeviceType device_name, class ALLOC,
+          template <typename, DeviceType, class> class MatrixA,
+          template <typename, DeviceType, class> class MatrixB,
+          template <typename, DeviceType, class> class MatrixC>
+inline void gemm(Scalar alpha, const MatrixA<Scalar, device_name, ALLOC>& a,
+                 const MatrixB<Scalar, device_name, ALLOC>& b, Scalar beta,
+                 MatrixC<Scalar, device_name, ALLOC>& c, int thread_id = 0, int stream_id = 0) {
   gemm<Scalar, device_name>('N', 'N', alpha, a, b, beta, c, thread_id, stream_id);
 }
 
@@ -734,7 +744,7 @@ void gemv(char transa, const Matrix<Scalar, CPU, ALLOC>& a, const Vector<Scalar,
 //                b.nrCols() == c.nrCols() if transb == 'N', b.nrRows() == c.nrCols() otherwise,
 //                ka == kb, where ka = a.nrCols() if transa == 'N', ka = a.nrRows() otherwise and
 //                          kb = b.nrRows() if transb == 'N', kb = b.nrCols() otherwise.
-  template <typename Scalar, DeviceType device_name, class ALLOC>
+template <typename Scalar, DeviceType device_name, class ALLOC>
 inline void gemm(char transa, char transb, const Matrix<Scalar, device_name, ALLOC>& a,
                  const Matrix<Scalar, device_name, ALLOC>& b, Matrix<Scalar, device_name, ALLOC>& c,
                  int thread_id = 0, int stream_id = 0) {
@@ -774,8 +784,8 @@ void trsm(char uplo, char diag, const Matrix<Scalar, device_name>& a,
 //                ka == kb, where ka = a.nrCols() if transa == 'N', ka = a.nrRows() otherwise and
 //                          kb = b.nrRows() if transb == 'N', kb = b.nrCols() otherwise.
 template <typename Scalar, class ALLOC>
-void gemm(char transa, char transb, Matrix<Scalar, CPU, ALLOC>& a, Matrix<std::complex<Scalar>, CPU, ALLOC>& b,
-          Matrix<std::complex<Scalar>, CPU, ALLOC>& c) {
+void gemm(char transa, char transb, Matrix<Scalar, CPU, ALLOC>& a,
+          Matrix<std::complex<Scalar>, CPU, ALLOC>& b, Matrix<std::complex<Scalar>, CPU, ALLOC>& c) {
   Matrix<Scalar, CPU, ALLOC> b_part(b.size());
   Matrix<Scalar, CPU, ALLOC> c_re(c.size());
   Matrix<Scalar, CPU, ALLOC> c_im(c.size());
@@ -858,7 +868,8 @@ static void gemm(char transa, char transb, Matrix<std::complex<Scalar>, CPU, ALL
 //                and kb = b[0].nrRows() if transb == 'N', kb = b[0].nrCols() otherwise.
 template <typename Scalar, class ALLOC>
 void multiply(char transa, char transb, const std::array<Matrix<Scalar, CPU, ALLOC>, 2>& a,
-              const std::array<Matrix<Scalar, CPU, ALLOC>, 2>& b, std::array<Matrix<Scalar, CPU, ALLOC>, 2>& c,
+              const std::array<Matrix<Scalar, CPU, ALLOC>, 2>& b,
+              std::array<Matrix<Scalar, CPU, ALLOC>, 2>& c,
               std::array<Matrix<Scalar, CPU, ALLOC>, 5>& work) {
   assert(a[0].size() == a[1].size());
   assert(b[0].size() == b[1].size());
@@ -895,7 +906,8 @@ void multiply(char transa, char transb, const std::array<Matrix<Scalar, CPU, ALL
 
 template <typename Scalar, class ALLOC>
 void multiply(const std::array<Matrix<Scalar, CPU, ALLOC>, 2>& a,
-              const std::array<Matrix<Scalar, CPU, ALLOC>, 2>& b, std::array<Matrix<Scalar, CPU, ALLOC>, 2>& c,
+              const std::array<Matrix<Scalar, CPU, ALLOC>, 2>& b,
+              std::array<Matrix<Scalar, CPU, ALLOC>, 2>& c,
               std::array<Matrix<Scalar, CPU, ALLOC>, 5>& work) {
   multiply('N', 'N', a, b, c, work);
 }
@@ -1018,8 +1030,9 @@ inline void multiplyDiagonalRight(const Matrix<Scalar, GPU>& a, const Vector<Sca
 // Postcondition: lambda_re, lambda_i, are resized, vl if jobvl == 'V', vr if jobvr == 'V' are
 // resized.
 template <typename Scalar, class ALLOC>
-void eigensolver(char jobvl, char jobvr, const Matrix<Scalar, CPU, ALLOC>& a, Vector<Scalar, CPU>& lambda_re,
-                 Vector<Scalar, CPU>& lambda_im, Matrix<Scalar, CPU, ALLOC>& vl, Matrix<Scalar, CPU, ALLOC>& vr) {
+void eigensolver(char jobvl, char jobvr, const Matrix<Scalar, CPU, ALLOC>& a,
+                 Vector<Scalar, CPU>& lambda_re, Vector<Scalar, CPU>& lambda_im,
+                 Matrix<Scalar, CPU, ALLOC>& vl, Matrix<Scalar, CPU, ALLOC>& vr) {
   assert(a.is_square());
 
   Matrix<Scalar, CPU, ALLOC> a_copy(a);
@@ -1055,8 +1068,9 @@ void eigensolver(char jobvl, char jobvr, const Matrix<Scalar, CPU, ALLOC>& a, Ve
 //               a is a square matrix.
 // Postcondition: lambda, is resized, vl if jobvl == 'V', vr if jobvr == 'V' are resized.
 template <typename Scalar, class ALLOC>
-void eigensolver(char jobvl, char jobvr, const Matrix<std::complex<Scalar>, CPU>& a,
-                 Vector<std::complex<Scalar>, CPU>& lambda, Matrix<std::complex<Scalar>, CPU, ALLOC>& vl,
+void eigensolver(char jobvl, char jobvr, const Matrix<std::complex<Scalar>, CPU, ALLOC>& a,
+                 Vector<std::complex<Scalar>, CPU>& lambda,
+                 Matrix<std::complex<Scalar>, CPU, ALLOC>& vl,
                  Matrix<std::complex<Scalar>, CPU, ALLOC>& vr) {
   assert(a.is_square());
 
@@ -1145,7 +1159,8 @@ void eigensolverHermitian(char jobv, char uplo, const Matrix<std::complex<Scalar
 }
 
 template <typename Scalar, class ALLOC>
-void eigensolverGreensFunctionMatrix(char jobv, char uplo, const Matrix<std::complex<Scalar>, CPU, ALLOC>& a,
+void eigensolverGreensFunctionMatrix(char jobv, char uplo,
+                                     const Matrix<std::complex<Scalar>, CPU, ALLOC>& a,
                                      Vector<Scalar, CPU>& lambda,
                                      Matrix<std::complex<Scalar>, CPU, ALLOC>& v) {
   assert(a.is_square());
@@ -1174,7 +1189,8 @@ void eigensolverGreensFunctionMatrix(char jobv, char uplo, const Matrix<std::com
 // Out: a_inv
 // Postconditions: a_inv is resized to the needed dimension.
 template <typename Scalar, class ALLOC>
-void pseudoInverse(const Matrix<Scalar, CPU, ALLOC>& a, Matrix<Scalar, CPU, ALLOC>& a_inv, double eps = 1.e-6) {
+void pseudoInverse(const Matrix<Scalar, CPU, ALLOC>& a, Matrix<Scalar, CPU, ALLOC>& a_inv,
+                   double eps = 1.e-6) {
   int m = a.nrRows();
   int n = a.nrCols();
   a_inv.resizeNoCopy(std::make_pair(n, m));
@@ -1236,7 +1252,7 @@ void pseudoInverse(const Matrix<Scalar, CPU, ALLOC>& a, Matrix<Scalar, CPU, ALLO
 // Computes (in place) the determinant of the matrix.
 // Returns: determinant.
 // Postcondition: M is its LU decomposition.
-  template <template <typename, DeviceType, class> class MatrixType, typename Scalar, class ALLOC>
+template <template <typename, DeviceType, class> class MatrixType, typename Scalar, class ALLOC>
 Scalar determinantIP(MatrixType<Scalar, CPU, ALLOC>& M) {
   assert(M.nrCols() == M.nrRows());
   const int n = M.nrCols();
@@ -1263,7 +1279,7 @@ Scalar determinantIP(MatrixType<Scalar, CPU, ALLOC>& M) {
 
 // Copy and computes the determinant of the matrix.
 // Returns: determinant.
-  template <typename Scalar, DeviceType device, class ALLOC>
+template <typename Scalar, DeviceType device, class ALLOC>
 Scalar determinant(const Matrix<Scalar, device>& M) {
   Matrix<Scalar, CPU, ALLOC> M_copy(M);
   return determinantIP(M_copy);
@@ -1305,7 +1321,8 @@ auto logDeterminantIP(MatrixType<Scalar, CPU, ALLOC>& M, std::vector<int>& ipiv)
   return std::make_pair(log_det, phase);
 }
 
-template <template <typename, DeviceType, class> class MatrixType, typename Scalar, DeviceType device, class ALLOC>
+template <template <typename, DeviceType, class> class MatrixType, typename Scalar,
+          DeviceType device, class ALLOC>
 auto logDeterminant(const MatrixType<Scalar, device, ALLOC>& m) {
   Matrix<Scalar, CPU, ALLOC> m_copy(m);
   std::vector<int> ipiv;
@@ -1330,8 +1347,8 @@ auto inverseAndLogDeterminant(MatrixType<Scalar, CPU, ALLOC>& mat) {
 }
 
 template <typename Scalar, class ALLOC>
-bool areNear(const Matrix<Scalar, CPU, ALLOC>& A, const Matrix<Scalar, CPU, ALLOC>& B, const double err = 1e-16) {
-
+bool areNear(const Matrix<Scalar, CPU, ALLOC>& A, const Matrix<Scalar, CPU, ALLOC>& B,
+             const double err = 1e-16) {
   if (A.size() != B.size())
     return false;
 
diff --git a/include/dca/linalg/util/allocators/pinned_allocator.hpp b/include/dca/linalg/util/allocators/pinned_allocator.hpp
index bdf0a108d..3e5353d54 100644
--- a/include/dca/linalg/util/allocators/pinned_allocator.hpp
+++ b/include/dca/linalg/util/allocators/pinned_allocator.hpp
@@ -17,7 +17,7 @@
 #include "dca/platform/dca_gpu.h"
 #else
 #error "This file requires GPU support."
-#endif 
+#endif
 
 namespace dca {
 namespace linalg {
@@ -34,11 +34,22 @@ class PinnedAllocator {
   using const_pointer = const T*;
   using value_type = T;
 
+  /** Allocates pinned memory through the GPU runtime.
+   *  It does this in the crudest way possible on HIP which is to just always allocate as many 4k
+   * blocks as the allocation can minimally fit, the rest of the memory is wasted.  We have no way
+   * of knowing if this was what the runtime always did or not.
+   */
   T* allocate(std::size_t n, const void* /*hint*/ = nullptr) {
     if (n == 0)
       return nullptr;
+#ifdef DCA_HAVE_HIP
+    std::size_t actual_size = ((n * sizeof(T)) / 4096 + 1) * 4096;
+#else
+    std::size_t actual_size = n * sizeof(T);
+#endif
     void* ptr;
-    cudaError_t ret = cudaHostAlloc(&ptr, n * sizeof(T), cudaHostAllocDefault);
+    cudaError_t ret =
+        cudaHostAlloc(&ptr, actual_size, cudaHostAllocPortable);  // cudaHostAllocDefault
     if (ret != cudaSuccess) {
       printErrorMessage(ret, __FUNCTION__, __FILE__, __LINE__,
                         "\t HOST size requested : " + std::to_string(n * sizeof(T)));
@@ -48,14 +59,19 @@ class PinnedAllocator {
   }
 
   void deallocate(T*& ptr, std::size_t /*n*/ = 0) noexcept {
+    //    assert(ptr == host_ptr_);
     cudaError_t ret = cudaFreeHost(ptr);
     if (ret != cudaSuccess) {
       printErrorMessage(ret, __FUNCTION__, __FILE__, __LINE__);
       std::terminate();
     }
+    //host_ptr_ = nullptr;
     ptr = nullptr;
   }
 
+private:
+  //std::size_t actual_size_ = 0;
+  //T* host_ptr_{nullptr};
 };
 
 // These are part of the requirements for a C++ Allocator however these are insufficient
@@ -63,12 +79,16 @@ class PinnedAllocator {
 // They are part of what's needed to do a std::move on a PinnedAllocator backed object
 // and avoid deallocation and reallocation.
 template <class T, class U>
-bool operator==(const PinnedAllocator<T>&, const PinnedAllocator<U>&) { return true; }
+bool operator==(const PinnedAllocator<T>&, const PinnedAllocator<U>&) {
+  return true;
+}
 template <class T, class U>
-bool operator!=(const PinnedAllocator<T>&, const PinnedAllocator<U>&) { return false; }
+bool operator!=(const PinnedAllocator<T>&, const PinnedAllocator<U>&) {
+  return false;
+}
 
-}  // util
-}  // linalg
-}  // dca
+}  // namespace util
+}  // namespace linalg
+}  // namespace dca
 
 #endif  // DCA_LINALG_UTIL_ALLOCATORS_PINNED_ALLOCATOR_HPP
diff --git a/include/dca/linalg/util/error_gpuBLAS.hpp b/include/dca/linalg/util/error_gpuBLAS.hpp
index 267be6b89..37b7c2e94 100644
--- a/include/dca/linalg/util/error_gpuBLAS.hpp
+++ b/include/dca/linalg/util/error_gpuBLAS.hpp
@@ -22,7 +22,7 @@
 #include "dca/linalg/util/error_cuda.hpp"
 #elif defined(DCA_HAVE_HIP)
 #include <hip/hip_runtime.h>
-#include <hipblas.h>
+#include <hipblas/hipblas.h>
 #include <hip/hip_complex.h>
 #include "dca/util/cuda2hip.h"
 #include "dca/linalg/util/error_hip.hpp"
diff --git a/include/dca/linalg/util/gpuBLAS_handles.hpp b/include/dca/linalg/util/gpuBLAS_handles.hpp
index eba9a41b0..dba2682aa 100644
--- a/include/dca/linalg/util/gpuBLAS_handles.hpp
+++ b/include/dca/linalg/util/gpuBLAS_handles.hpp
@@ -21,7 +21,7 @@
 #include <cublas_v2.h>
 #elif defined(DCA_HAVE_HIP)
 #include <hip/hip_runtime.h>
-#include <hipblas.h>
+#include <hipblas/hipblas.h>
 #include <hip/hip_complex.h>
 #include "dca/util/cuda2hip.h"
 #endif
diff --git a/include/dca/linalg/util/util_matrixop.hpp b/include/dca/linalg/util/util_matrixop.hpp
index 180f8f5e1..bd8e26400 100644
--- a/include/dca/linalg/util/util_matrixop.hpp
+++ b/include/dca/linalg/util/util_matrixop.hpp
@@ -33,8 +33,8 @@ namespace util {
 
 // Returns optimal lwork for inverse.
 // In: mat
-  template <typename ScalarType, class ALLOC, template <typename, DeviceType, class> class MatrixType>
-int getInverseWorkSize(MatrixType<ScalarType, CPU, ALLOC>& mat) {
+  template <typename ScalarType, template <typename> class ALLOC, template <typename, DeviceType, class> class MatrixType>
+  int getInverseWorkSize(MatrixType<ScalarType, CPU, ALLOC<ScalarType>>& mat) {
   assert(mat.is_square());
 
   ScalarType tmp;
@@ -42,8 +42,8 @@ int getInverseWorkSize(MatrixType<ScalarType, CPU, ALLOC>& mat) {
   return lapack::util::getWorkSize(tmp);
 }
 #ifdef DCA_HAVE_GPU
-template <typename ScalarType,  template <typename, DeviceType> class MatrixType>
-int getInverseWorkSize(const MatrixType<ScalarType, GPU>& mat) {
+  template <typename ScalarType,  template<typename> class ALLOC, template <typename, DeviceType, class> class MatrixType>
+  int getInverseWorkSize(const MatrixType<ScalarType, GPU, ALLOC<ScalarType>>& mat) {
   assert(mat.is_square());
 
   return mat.nrRows() * magma::get_getri_nb<ScalarType>(mat.nrRows());
diff --git a/include/dca/math/function_transform/domainwise_function_transform.hpp b/include/dca/math/function_transform/domainwise_function_transform.hpp
index a6310903a..06ce4cb11 100644
--- a/include/dca/math/function_transform/domainwise_function_transform.hpp
+++ b/include/dca/math/function_transform/domainwise_function_transform.hpp
@@ -112,10 +112,10 @@ struct DomainwiseFunctionTransform {
     }
   }
 
-  template <typename scalartype_input, typename scalartype_output, typename scalartype_T>
+  template <typename scalartype_input, typename scalartype_output, typename scalartype_T, class ALLOC>
   static void execute_on_all(const func::function<scalartype_input, domain_input>& f_input,
                              func::function<scalartype_output, domain_output>& f_output,
-                             const linalg::Matrix<scalartype_T, linalg::CPU>& T) {
+                             const linalg::Matrix<scalartype_T, linalg::CPU, ALLOC>& T) {
     if (VERBOSE) {
       std::cout << "\n\n\t" << __FUNCTION__ << "\t" << f_input.get_name() << " --> "
                 << f_output.get_name() << "\n\n";
diff --git a/include/dca/math/function_transform/function_transform.hpp b/include/dca/math/function_transform/function_transform.hpp
index 97dd0c8a4..38a836967 100644
--- a/include/dca/math/function_transform/function_transform.hpp
+++ b/include/dca/math/function_transform/function_transform.hpp
@@ -85,10 +85,10 @@ class FunctionTransform<func::dmn_0<type_input>, func::dmn_0<type_output>> {
   }
 
   template <typename scalartype_input, class domain_input, typename scalartype_output,
-            class domain_output, typename scalartype_T>
+            class domain_output, typename scalartype_T, class ALLOC>
   static void execute_on_all(const func::function<scalartype_input, domain_input>& f_input,
                              func::function<scalartype_output, domain_output>& f_output,
-                             const linalg::Matrix<scalartype_T, linalg::CPU>& T) {
+                             const linalg::Matrix<scalartype_T, linalg::CPU, ALLOC>& T) {
     if (VERBOSE)
       print_types(f_input, f_output);
 
diff --git a/include/dca/math/function_transform/transform_domain_procedure.hpp b/include/dca/math/function_transform/transform_domain_procedure.hpp
index e78ade316..10729d5e7 100644
--- a/include/dca/math/function_transform/transform_domain_procedure.hpp
+++ b/include/dca/math/function_transform/transform_domain_procedure.hpp
@@ -30,31 +30,30 @@ class TRANSFORM_DOMAIN_PROCEDURE {
   static void characterize_transformation(const f_input_t& f_input, const f_output_t& f_output,
                                           int& M, int& K, int& N, int& P);
 
-  template <typename scalartype, class domain_input, class domain_output>
-  static void transform(const func::function<scalartype, domain_input>& f_input,
-                        func::function<scalartype, domain_output>& f_output,
-                        const linalg::Matrix<scalartype, linalg::CPU>& T);
+  template <typename Scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
+  static void transform(const func::function<Scalartype, domain_input>& f_input,
+                        func::function<Scalartype, domain_output>& f_output,
+                        const linalg::Matrix<Scalartype, linalg::CPU, ALLOC<Scalartype>>& tm);
 
-  template <typename scalartype, class domain_input, class domain_output>
+  template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
   static void transform(const func::function<std::complex<scalartype>, domain_input>& f_input,
                         func::function<std::complex<scalartype>, domain_output>& f_output,
-                        const linalg::Matrix<scalartype, linalg::CPU>& T);
+                        const linalg::Matrix<scalartype, linalg::CPU, ALLOC<scalartype>>& tm);
 
-  template <typename Scalartype, typename Scalartype2, class DomainInput, class DomainOutput,
-            class = std::enable_if_t<std::is_scalar<Scalartype>::value>>
+  template <typename Scalartype, class DomainInput, class DomainOutput, template<typename> class ALLOC>
   static void transform(const func::function<Scalartype, DomainInput>& f_input,
                         func::function<Scalartype, DomainOutput>& f_output,
-                        const linalg::Matrix<std::complex<Scalartype2>, linalg::CPU>& T);
+                        const linalg::Matrix<std::complex<Scalartype>, linalg::CPU, ALLOC<std::complex<Scalartype>>>& tm);
 
-  template <typename scalartype, class domain_input, class domain_output>
+  template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
   static void transform(const func::function<scalartype, domain_output>& f_input,
                         func::function<std::complex<scalartype>, domain_input>& f_output,
-                        const linalg::Matrix<std::complex<scalartype>, linalg::CPU>& T);
+                        const linalg::Matrix<std::complex<scalartype>, linalg::CPU, ALLOC<std::complex<scalartype>>>& tm);
 
-  template <typename scalartype, class domain_input, class domain_output>
+  template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
   static void transform(const func::function<std::complex<scalartype>, domain_input>& f_input,
                         func::function<scalartype, domain_output>& f_output,
-                        const linalg::Matrix<std::complex<scalartype>, linalg::CPU>& T);
+                        const linalg::Matrix<std::complex<scalartype>, linalg::CPU, ALLOC<std::complex<scalartype>>>& T);
 };
 
 template <int DMN_INDEX>
@@ -76,11 +75,11 @@ void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::characterize_transformation(const f_
 }
 
 template <int DMN_INDEX>
-template <typename scalartype, class domain_input, class domain_output>
+template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
 void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
     const func::function<scalartype, domain_input>& f_input,
     func::function<scalartype, domain_output>& f_output,
-    const linalg::Matrix<scalartype, linalg::CPU>& T) {
+    const linalg::Matrix<scalartype, linalg::CPU, ALLOC<scalartype>>& T) {
   int M, K, N, P;
   characterize_transformation(f_input, f_output, M, K, N, P);
 
@@ -104,11 +103,11 @@ void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
 }
 
 template <int DMN_INDEX>
-template <typename scalartype, class domain_input, class domain_output>
+template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
 void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
     const func::function<std::complex<scalartype>, domain_input>& f_input,
     func::function<std::complex<scalartype>, domain_output>& f_output,
-    const linalg::Matrix<scalartype, linalg::CPU>& T) {
+    const linalg::Matrix<scalartype, linalg::CPU, ALLOC<scalartype>>& T) {
   int M, K, N, P;
   characterize_transformation(f_input, f_output, M, K, N, P);
 
@@ -125,12 +124,13 @@ void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
 }
 
 template <int DMN_INDEX>
-template <typename Scalartype, typename Scalartype2, class DomainInput, class DomainOutput, class>
+template <typename Scalartype,  class DomainInput, class DomainOutput, template<typename> class ALLOC>
 void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
     const func::function<Scalartype, DomainInput>& f_input,
     func::function<Scalartype, DomainOutput>& f_output,
-    const linalg::Matrix<std::complex<Scalartype2>, linalg::CPU>& T) {
-  linalg::Matrix<Scalartype, linalg::CPU> T_re("T_re", T.size());
+    const linalg::Matrix<std::complex<Scalartype>, linalg::CPU, ALLOC<std::complex<Scalartype>>>& T) {
+  auto matrix_size = T.size();
+  linalg::Matrix<dca::util::RealAlias<Scalartype>, linalg::CPU, ALLOC<dca::util::RealAlias<Scalartype>>> T_re("T_re", matrix_size);
 
   for (int j = 0; j < T.size().second; j++)
     for (int i = 0; i < T.size().first; i++)
@@ -140,28 +140,28 @@ void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
 }
 
 template <int DMN_INDEX>
-template <typename scalartype, class domain_input, class domain_output>
+template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
 void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
     const func::function<scalartype, domain_output>& f_input,
     func::function<std::complex<scalartype>, domain_input>& f_output,
-    const linalg::Matrix<std::complex<scalartype>, linalg::CPU>& T) {
+    const linalg::Matrix<std::complex<scalartype>, linalg::CPU, ALLOC<std::complex<scalartype>>>& T) {
   const auto f_cmplx = func::util::complex(f_input);
 
   transform(f_cmplx, f_output, T);
 }
 
 template <int DMN_INDEX>
-template <typename scalartype, class domain_input, class domain_output>
+template <typename scalartype, class domain_input, class domain_output, template<typename> class ALLOC>
 void TRANSFORM_DOMAIN_PROCEDURE<DMN_INDEX>::transform(
     const func::function<std::complex<scalartype>, domain_input>& f_input,
     func::function<scalartype, domain_output>& f_output,
-    const linalg::Matrix<std::complex<scalartype>, linalg::CPU>& T) {
+    const linalg::Matrix<std::complex<scalartype>, linalg::CPU, ALLOC<std::complex<scalartype>>>& T) {
   f_output = 0.;
 
   func::function<scalartype, domain_input> f_in("f_in");
   func::function<scalartype, domain_output> f_out("f_out");
 
-  linalg::Matrix<scalartype, linalg::CPU> T_tmp("T_tmp", T.size());
+  linalg::Matrix<scalartype, linalg::CPU, ALLOC<scalartype>> T_tmp("T_tmp", T.size());
 
   {
     for (int i = 0; i < f_input.size(); i++)
diff --git a/include/dca/parallel/mpi_concurrency/mpi_concurrency.hpp b/include/dca/parallel/mpi_concurrency/mpi_concurrency.hpp
index 529207e23..4408dd528 100644
--- a/include/dca/parallel/mpi_concurrency/mpi_concurrency.hpp
+++ b/include/dca/parallel/mpi_concurrency/mpi_concurrency.hpp
@@ -87,15 +87,18 @@ class MPIConcurrency final : public virtual MPIInitializer,
 
 #ifdef DCA_HAVE_ADIOS2
   adios2::ADIOS& get_adios() {
-    return adios_;
+    return *adios_;
   }
 #endif
 
+  //friend class adios2_reader<MPIConcurrency>;
+  
 private:
   constexpr static char parallel_type_str_[] = "MPIConcurrency";
 
+public:
 #ifdef DCA_HAVE_ADIOS2
-  adios2::ADIOS adios_;
+  std::unique_ptr<adios2::ADIOS> adios_;
 #endif
 };
 
diff --git a/include/dca/parallel/no_concurrency/no_concurrency.hpp b/include/dca/parallel/no_concurrency/no_concurrency.hpp
index 03c3fb2f4..0b24b152c 100644
--- a/include/dca/parallel/no_concurrency/no_concurrency.hpp
+++ b/include/dca/parallel/no_concurrency/no_concurrency.hpp
@@ -29,7 +29,7 @@ class NoConcurrency : public SerialCollectiveSum, public SerialProcessorGrouping
 public:
   using Gang = SerialGang;
 
-  NoConcurrency(int /*argc*/, char** /*argv*/){};
+  NoConcurrency(int /*argc*/, char** /*argv*/);
 
   void abort() const;
 
@@ -64,7 +64,7 @@ class NoConcurrency : public SerialCollectiveSum, public SerialProcessorGrouping
 
 #ifdef DCA_HAVE_ADIOS2
   adios2::ADIOS& get_adios() {
-    return adios_;
+    return *adios_;
   }
 #endif
 
@@ -74,8 +74,9 @@ class NoConcurrency : public SerialCollectiveSum, public SerialProcessorGrouping
 
 private:
   constexpr static char parallel_type_str_[] = "NoConcurrency";
+public:
 #ifdef DCA_HAVE_ADIOS2
-  adios2::ADIOS adios_;
+  std::unique_ptr<adios2::ADIOS> adios_;
 #endif
 
 };
diff --git a/include/dca/phys/dca_data/dca_data.hpp b/include/dca/phys/dca_data/dca_data.hpp
index 490c76019..8e5b05ccd 100644
--- a/include/dca/phys/dca_data/dca_data.hpp
+++ b/include/dca/phys/dca_data/dca_data.hpp
@@ -116,9 +116,6 @@ class DcaData {
   /** These reads are used by analysis programs only for now.
    */
   void read(const std::string& filename);
-#ifdef DCA_HAVE_ADIOS2
-  void read(adios2::ADIOS& adios, std::string filename);
-#endif
 
   /** prefer this as it allows for more sensible handling of appendable files like bp4.
    *  currently this is used by chi_q_omega only.
@@ -135,11 +132,6 @@ class DcaData {
   void initialize();
   void initializeH0_and_H_i();
   void initialize_G0();
-#ifdef DCA_HAVE_ADIOS2
-  /** read initializeSigma from adios file, it is probably already open.
-   */
-  void initializeSigma(adios2::ADIOS& adios, const std::string& filename);
-#endif
   void initializeSigma(const std::string& filename);
   void readSigmaFile(io::Reader<Concurrency>& reader);
 
@@ -150,7 +142,7 @@ class DcaData {
 
 private:
   Parameters& parameters_;
-  const Concurrency& concurrency_;
+  Concurrency& concurrency_;
 
 public:
   func::function<int, NuNuDmn> H_symmetry;
@@ -367,31 +359,6 @@ void DcaData<Parameters, DT>::read(const std::string& filename) {
   }
 }
 
-#ifdef DCA_HAVE_ADIOS2
-template <class Parameters, DistType DT>
-void DcaData<Parameters, DT>::read(adios2::ADIOS& adios, std::string filename) {
-  if (concurrency_.id() == concurrency_.first()) {
-    std::cout << "\n\n\t starts reading \n\n";
-
-    dca::io::Reader<typename Parameters::concurrency_type> reader(adios, concurrency_,
-                                                                  parameters_.get_output_format());
-    reader.open_file(filename);
-    read(reader);
-    reader.close_file();
-  }
-  // realize that if you need anyting else from the prior data that you need to add its broadcast here.
-  concurrency_.broadcast(parameters_.get_chemical_potential());
-  concurrency_.broadcast_object(Sigma);
-
-  if (parameters_.isAccumulatingG4()) {
-    concurrency_.broadcast_object(G_k_w);
-
-    for (auto& G4_channel : G4_)
-      concurrency_.broadcast_object(G4_channel);
-  }
-}
-#endif
-
 template <class Parameters, DistType DT>
 void DcaData<Parameters, DT>::read(dca::io::Reader<typename Parameters::concurrency_type>& reader) {
   std::size_t step_count = reader.getStepCount();
@@ -616,65 +583,56 @@ void DcaData<Parameters, DT>::initialize_G0() {
   }
 }
 
-#ifdef DCA_HAVE_ADIOS2
-template <class Parameters, DistType DT>
-void DcaData<Parameters, DT>::initializeSigma(adios2::ADIOS& adios [[maybe_unused]],
-                                              const std::string& filename) {
-  if (concurrency_.id() == concurrency_.first()) {
-    std::cout << "reading Sigma File\n";
-    io::IOType sigma_file_io = io::extensionToIOType(filename);
-    io::Reader reader(concurrency_, sigma_file_io);
-    reader.open_file(filename);
-    std::size_t step_count = reader.getStepCount();
-    for (std::size_t i = 0; i < step_count; ++i) {
-      reader.begin_step();
-      reader.end_step();
-    }
-    readSigmaFile(reader);
-    reader.close_file();
-  }
-  concurrency_.broadcast(parameters_.get_chemical_potential());
-  concurrency_.broadcast(Sigma);
-}
-#endif
-
-// Strong assumption here that this only gets called for HDF5 input.
 template <class Parameters, DistType DT>
 void DcaData<Parameters, DT>::initializeSigma(const std::string& filename) {
   if (concurrency_.id() == concurrency_.first()) {
     std::cout << "reading Sigma File\n";
     io::IOType sigma_file_io = io::extensionToIOType(filename);
     io::Reader reader(concurrency_, sigma_file_io);
-    int hdf5_last_iteration = -1;
-    reader.open_file(filename);
-    long step_count = reader.getStepCount();
-    // Work around odd way hdf5 steps get written
-    int completed_iteration = 0;
-    std::cout << step_count << " steps found.\n";
-    if (step_count >= 0) {
-    find_step:
-      for (std::size_t i = 0; i < step_count; ++i) {
-        reader.begin_step();
-        std::cerr << "current step " << i << '\n';
-        bool has_iteration =
-            reader.execute("DCA-loop-functions/completed-iteration", completed_iteration);
-        std::cerr << "completed_iteration " << completed_iteration << '\n';
-        if (has_iteration && (i > completed_iteration)) {
-          std::cerr << "past complete iterations " << completed_iteration << "at step " << i << '\n';
-          hdf5_last_iteration = i;
-          reader.close_file();
-          reader.open_file(filename);
-          step_count = hdf5_last_iteration;
-          goto find_step;
+    if (sigma_file_io == io::IOType::HDF5 || sigma_file_io == io::IOType::ADIOS2) {
+      io::Reader reader(concurrency_, sigma_file_io);
+      int hdf5_last_iteration = -1;
+      reader.open_file(filename);
+      long step_count = reader.getStepCount();
+      // Work around odd way hdf5 steps get written
+      int completed_iteration = 0;
+      std::cout << step_count << " steps found.\n";
+      if (step_count >= 0) {
+      find_step:
+        for (std::size_t i = 0; i < step_count; ++i) {
+          reader.begin_step();
+          std::cerr << "current step " << i << '\n';
+          bool has_iteration =
+              reader.execute("DCA-loop-functions/completed-iteration", completed_iteration);
+          std::cerr << "completed_iteration " << completed_iteration << '\n';
+          if (has_iteration && (i > completed_iteration)) {
+            std::cerr << "past complete iterations " << completed_iteration << "at step " << i
+                      << '\n';
+            hdf5_last_iteration = i;
+            reader.close_file();
+            reader.open_file(filename);
+            step_count = hdf5_last_iteration;
+            goto find_step;
+          }
+          if (i < step_count - 1)
+            reader.end_step();
         }
-        if (i < step_count - 1)
-          reader.end_step();
+        readSigmaFile(reader);
+        reader.end_step();
+      }
+      else {
+        readSigmaFile(reader);
       }
-      readSigmaFile(reader);
-      reader.end_step();
     }
     else {
+      reader.open_file(filename);
+      std::size_t step_count = reader.getStepCount();
+      for (std::size_t i = 0; i < step_count; ++i) {
+        reader.begin_step();
+        reader.end_step();
+      }
       readSigmaFile(reader);
+      reader.close_file();
     }
   }
   concurrency_.broadcast(parameters_.get_chemical_potential());
diff --git a/include/dca/phys/dca_loop/dca_loop.hpp b/include/dca/phys/dca_loop/dca_loop.hpp
index 195b98c64..40d23505c 100644
--- a/include/dca/phys/dca_loop/dca_loop.hpp
+++ b/include/dca/phys/dca_loop/dca_loop.hpp
@@ -56,7 +56,6 @@ class DcaLoop {
   using Lattice = typename ParametersType::lattice_type;
   using Real = typename ParametersType::Real;
 
-  
   using b = func::dmn_0<domains::electron_band_domain>;
   using s = func::dmn_0<domains::electron_spin_domain>;
   using k_DCA =
@@ -143,18 +142,8 @@ DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::DcaLoop(
       cluster_mapping_obj(parameters),
       lattice_mapping_obj(parameters),
       update_chemical_potential_obj(parameters, MOMS, cluster_mapping_obj),
-#ifdef DCA_HAVE_ADIOS2
-      output_file_(
-          parameters.get_output_format() == "ADIOS2"
-              ? std::make_shared<io::Writer<concurrency_type>>(
-                    concurrency.get_adios(), concurrency_ref, parameters.get_output_format(), false)
-              : std::make_shared<io::Writer<concurrency_type>>(
-                           concurrency.get_adios(), concurrency_ref, parameters.get_output_format(),
-                           false)),
-#else
       output_file_(std::make_shared<io::Writer<concurrency_type>>(
-								  concurrency_ref, parameters.get_output_format(), false)),
-#endif
+          concurrency_ref, parameters.get_output_format(), false)),
       monte_carlo_integrator_(parameters_ref, MOMS_ref, output_file_) {
   file_name_ = parameters.get_directory() + parameters.get_filename_dca();
 
@@ -169,7 +158,7 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::write() {
   // We assume DCALoop write is called once at the end of the run and we leave a step open for it to
   // write into.
 
-  //output_file_->begin_step();
+  // output_file_->begin_step();
 
   if (concurrency.id() == concurrency.first()) {
     // This should probably happen first not at the end
@@ -208,15 +197,8 @@ void DcaLoop<ParametersType, DDT, MCIntegratorType, DIST>::initialize() {
   int last_completed = -1;
   auto& autoresume_filename = parameters.get_autoresume_filename();
   if (parameters.autoresume()) {
-#ifdef DCA_HAVE_ADIOS2
-    io::IOType iotype = io::extensionToIOType(autoresume_filename);
-    if (iotype == io::IOType::ADIOS2)
-      last_completed = DCA_info_struct.readData(autoresume_filename, parameters.get_output_format(),
-                                                concurrency, concurrency.get_adios());
-    else
-#endif
-      last_completed =
-          DCA_info_struct.readData(autoresume_filename, parameters.get_output_format(), concurrency);
+    last_completed =
+        DCA_info_struct.readData(autoresume_filename, parameters.get_output_format(), concurrency);
 
     if (last_completed >= 0) {
       if (concurrency.id() == concurrency.first())
@@ -225,22 +207,12 @@ void DcaLoop<ParametersType, DDT, MCIntegratorType, DIST>::initialize() {
                   << std::endl;
 
       dca_iteration_ = std::min(last_completed + 1, parameters.get_dca_iterations() - 1);
-#ifdef DCA_HAVE_ADIOS2
-      if (iotype == io::IOType::ADIOS2)
-        MOMS.initializeSigma(concurrency.get_adios(), autoresume_filename);
-      else
-#endif
-        MOMS.initializeSigma(autoresume_filename);
+      MOMS.initializeSigma(autoresume_filename);
       perform_lattice_mapping();
     }
   }
   else if (parameters.get_initial_self_energy() != "zero") {
-#ifdef DCA_HAVE_ADIOS2
-    if (io::extensionToIOType(parameters.get_initial_self_energy()) == io::IOType::ADIOS2)
-      MOMS.initializeSigma(concurrency.get_adios(), parameters.get_initial_self_energy());
-    else
-#endif
-      MOMS.initializeSigma(parameters.get_initial_self_energy());
+    MOMS.initializeSigma(parameters.get_initial_self_energy());
     perform_lattice_mapping();
   }
 
@@ -302,16 +274,10 @@ void DcaLoop<ParametersType, DcaDataType, MCIntegratorType, DIST>::execute() {
     }
 
     if (parameters.do_not_update_sigma()) {
-	if (parameters.get_initial_self_energy() == "zero")
-	  throw std::runtime_error("If there is no initial self energy it must be updated!");
-#ifdef DCA_HAVE_ADIOS2
-	if (io::extensionToIOType(parameters.get_initial_self_energy()) == io::IOType::ADIOS2)
-	  MOMS.initializeSigma(concurrency.get_adios(), parameters.get_initial_self_energy());
-	else
-#endif
-	  MOMS.initializeSigma(parameters.get_initial_self_energy());
+      if (parameters.get_initial_self_energy() == "zero")
+        throw std::runtime_error("If there is no initial self energy it must be updated!");
+      MOMS.initializeSigma(parameters.get_initial_self_energy());
     }
-
   }
 }
 
diff --git a/include/dca/phys/dca_loop/dca_loop_data.hpp b/include/dca/phys/dca_loop/dca_loop_data.hpp
index f61a38e4a..76e128f23 100644
--- a/include/dca/phys/dca_loop/dca_loop_data.hpp
+++ b/include/dca/phys/dca_loop/dca_loop_data.hpp
@@ -39,7 +39,7 @@ class DcaLoopData {
 
   using Real = typename Parameters::Real;
   using Scalar = typename Parameters::Scalar;
-  
+
   using b = func::dmn_0<domains::electron_band_domain>;
   using s = func::dmn_0<domains::electron_spin_domain>;
   using nu = func::dmn_variadic<b, s>;  // orbital-spin index
@@ -56,10 +56,6 @@ class DcaLoopData {
   // the last completed iteration from the input file, otherwise it returns -1.
   int readData(const std::string& filename, const std::string& format,
                const Concurrency& concurrency);
-#ifdef DCA_HAVE_ADIOS2
-  int readData(const std::string& filename, const std::string& format,
-               const Concurrency& concurrency, adios2::ADIOS& adios);
-#endif
   template <class READER>
   void readLoopDataCommon(READER& reader, const std::string& filename, const std::string& format);
   func::function<double, DCA_iteration_domain_type> Gflop_per_mpi_task;
@@ -160,39 +156,24 @@ void DcaLoopData<Parameters>::write(Writer& writer, Concurrency& concurrency) {
 
 template <typename Parameters>
 int DcaLoopData<Parameters>::readData(const std::string& filename, const std::string& format,
-                                          const Concurrency& concurrency) {
+                                      const Concurrency& concurrency) {
   if (concurrency.id() == concurrency.first() && filesystem::exists(filename)) {
     io::Reader reader(concurrency, format, false);
-    readLoopDataCommon(reader, filename, format);
-  }
-  concurrency.broadcast(last_completed_iteration);
-  return last_completed_iteration;
-}
-
-#ifdef DCA_HAVE_ADIOS2
-template <typename Parameters>
-int DcaLoopData<Parameters>::readData(const std::string& filename, const std::string& format,
-                                          const Concurrency& concurrency, adios2::ADIOS& adios) {
-  std::cout << "Reading dca_loop_data with ADIOS2\n";
-  if (concurrency.id() == concurrency.first() && filesystem::exists(filename)) {
-    io::Reader reader(adios, concurrency, format, false);
-    auto& adios2_reader = std::get<io::ADIOS2Reader<Concurrency>>(reader.getUnderlying());
-    std::size_t step_count = adios2_reader.getStepCount();
+    std::size_t step_count = reader.getStepCount();
     for (std::size_t i = 0; i < step_count; ++i) {
-      adios2_reader.begin_step();
-      adios2_reader.end_step();
+      reader.begin_step();
+      reader.end_step();
     }
     readLoopDataCommon(reader, filename, format);
   }
   concurrency.broadcast(last_completed_iteration);
   return last_completed_iteration;
 }
-#endif
 
 template <typename Parameters>
 template <class READER>
 void DcaLoopData<Parameters>::readLoopDataCommon(READER& reader, const std::string& filename,
-                                                     const std::string& format [[maybe_unused]]) {
+                                                 const std::string& format [[maybe_unused]]) {
   reader.open_file(filename);
   reader.open_group("DCA-loop-functions");
 
diff --git a/include/dca/phys/parameters/output_parameters.hpp b/include/dca/phys/parameters/output_parameters.hpp
index 2be6d4915..864a52421 100644
--- a/include/dca/phys/parameters/output_parameters.hpp
+++ b/include/dca/phys/parameters/output_parameters.hpp
@@ -29,10 +29,9 @@ class OutputParameters {
       : directory_("./"),
 	autoresume_(false),
 #ifdef DCA_HAVE_ADIOS2
-        output_format_("ADIOS2"),
-        g4_output_format_("ADIOS2"),
-        filename_g4_("dca_g4.bp"),
-        filename_dca_("dca.bp"),
+        output_format_("HDF5"),
+        g4_output_format_(""),
+        filename_dca_("dca.hdf5"),
         filename_analysis_("sofqomega.bp"),
 #else
         output_format_("HDF5"),
diff --git a/include/dca/testing/gtest_h_w_warning_blocking.h b/include/dca/testing/gtest_h_w_warning_blocking.h
index 14e670b0b..6bbd9615a 100644
--- a/include/dca/testing/gtest_h_w_warning_blocking.h
+++ b/include/dca/testing/gtest_h_w_warning_blocking.h
@@ -8,14 +8,25 @@
 // Author: Peter W. Doak  doakpw@ornl.gov  Oak Ridge National Lab
 //
 // This file implements a work around for at least older version of google test
-// giving little effort to removing warning trips
+// tripping newer compiler warnings. Eventually we will upgrade our bundled version
+// of google test or migrate to a different test framework.
 
-#ifndef DCA_GTEST_H_W_WARNING_BLOCKING_
+#ifndef DCA_GTEST_H_W_WARNING_BLOCKING_H
 #define DCA_GTEST_H_W_WARNING_BLOCKING_H
 
+#if defined(__clang__)
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunknown-warning-option"
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#pragma GCC diagnostic ignored "-Wunknown-pragmas"
+#else
 #pragma GCC diagnostic push
 #pragma GCC diagnostic ignored "-Wsuggest-override"
+#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
+#endif
+
 #include "gtest/gtest.h"
+
 #pragma GCC diagnostic_pop
 
 #endif
diff --git a/include/dca/util/containers/random_access_map_iterator.hpp b/include/dca/util/containers/random_access_map_iterator.hpp
index 041972573..32388d7ce 100644
--- a/include/dca/util/containers/random_access_map_iterator.hpp
+++ b/include/dca/util/containers/random_access_map_iterator.hpp
@@ -12,6 +12,8 @@
 #ifndef DCA_UTIL_TYPE_RANDOM_ACCESS_MAP_ITERATOR_HPP
 #define DCA_UTIL_TYPE_RANDOM_ACCESS_MAP_ITERATOR_HPP
 
+#include <cstdint>
+
 namespace dca {
 namespace util {
 namespace details {
diff --git a/include/dca/util/cuda2hip.h b/include/dca/util/cuda2hip.h
index 4a87e6a59..4de6d8caa 100644
--- a/include/dca/util/cuda2hip.h
+++ b/include/dca/util/cuda2hip.h
@@ -136,6 +136,7 @@
 #define cudaHostAlloc                   hipHostMalloc
 #define cudaHostAllocMapped             hipHostMallocMapped
 #define cudaHostAllocDefault            hipHostMallocDefault
+#define cudaHostAllocPortable           hipHostMallocPortable
 #define cudaIpcGetMemHandle             hipIpcGetMemHandle
 #define cudaIpcMemHandle_t              hipIpcMemHandle_t
 #define cudaIpcMemLazyEnablePeerAccess  hipIpcMemLazyEnablePeerAccess
diff --git a/include/dca/util/type_help.hpp b/include/dca/util/type_help.hpp
index bff88cb98..b9a095b5a 100644
--- a/include/dca/util/type_help.hpp
+++ b/include/dca/util/type_help.hpp
@@ -13,6 +13,7 @@
 #include <type_traits>
 #include <complex>
 #include <tuple>
+#include <cstdint>
 
 // namespace dca {
 // namespace util {
diff --git a/libs/googletest-1.8.0/include/gtest/gtest.h b/libs/googletest-1.8.0/include/gtest/gtest.h
index c97fc99aa..051efad5a 100644
--- a/libs/googletest-1.8.0/include/gtest/gtest.h
+++ b/libs/googletest-1.8.0/include/gtest/gtest.h
@@ -55,9 +55,6 @@
 #include <ostream>
 #include <vector>
 
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wsuggest-override"
-#pragma GCC diagnostic ignored "-Wmaybe-uninitialized"
 
 #include "gtest/internal/gtest-internal.h"
 #include "gtest/internal/gtest-string.h"
@@ -69,8 +66,6 @@
 #include "gtest/gtest-test-part.h"
 #include "gtest/gtest-typed-test.h"
 
-#pragma GCC diagnostic_pop
-
 // Depending on the platform, different string classes are available.
 // On Linux, in addition to ::std::string, Google also makes use of
 // class ::string, which has the same interface as ::std::string, but
diff --git a/src/io/CMakeLists.txt b/src/io/CMakeLists.txt
index bba962353..52e0337ee 100644
--- a/src/io/CMakeLists.txt
+++ b/src/io/CMakeLists.txt
@@ -4,12 +4,15 @@ add_subdirectory(json)
 #don't use this variable use the dca_io target
 set(DCA_IO_LIBS "dca_hdf5;json")
 
+add_library(dca_io reader.cpp io_types.cpp)
+
 if (DCA_HAVE_ADIOS2)
   add_subdirectory(adios2)
   set(DCA_IO_LIBS "${DCA_IO_LIBS};dca_adios2")
+  target_link_libraries(dca_io PUBLIC dca_adios2)
 endif()
 
-add_library(dca_io reader.cpp io_types.cpp)
+
 
 dca_gpu_runtime_link(dca_io)
 dca_gpu_blas_link(dca_io)
diff --git a/src/io/adios2/adios2_reader.cpp b/src/io/adios2/adios2_reader.cpp
index bb4e22b20..f7156165a 100644
--- a/src/io/adios2/adios2_reader.cpp
+++ b/src/io/adios2/adios2_reader.cpp
@@ -22,12 +22,8 @@ namespace io {
 // dca::io::
 
 template <class CT>
-ADIOS2Reader<CT>::ADIOS2Reader(adios2::ADIOS& adios, const CT* concurrency, bool verbose)
-    : adios_(adios), verbose_(verbose), concurrency_(concurrency) {}
-
-template <class CT>
-ADIOS2Reader<CT>::ADIOS2Reader(const CT* concurrency, bool verbose)
-    : adios_(GlobalAdios::getAdios()), verbose_(verbose), concurrency_(concurrency) {}
+ADIOS2Reader<CT>::ADIOS2Reader(CT& concurrency, bool verbose)
+  : concurrency_(concurrency), adios_(*concurrency.adios_), verbose_(verbose)  {}
 
 template <class CT>
 ADIOS2Reader<CT>::~ADIOS2Reader() {
@@ -43,6 +39,7 @@ void ADIOS2Reader<CT>::open_file(const std::string& file_name) {
   io_name_ = file_name;
   file_name_ = file_name;
   io_ = adios_.DeclareIO(io_name_);
+  io_.SetEngine("BP4");
   file_ = io_.Open(file_name_, adios2::Mode::Read);
 }
 
@@ -154,8 +151,10 @@ void ADIOS2Reader<CT>::dumpAvailableVars() {
 }
 
 template class ADIOS2Reader<dca::parallel::NoConcurrency>;
+template class ADIOS2Reader<dca::parallel::NoConcurrency const>;
 #ifdef DCA_HAVE_MPI
 template class ADIOS2Reader<dca::parallel::MPIConcurrency>;
+template class ADIOS2Reader<dca::parallel::MPIConcurrency const>;
 #endif
 
 }  // namespace io
diff --git a/src/io/adios2/adios2_writer.cpp b/src/io/adios2/adios2_writer.cpp
index c3f636c35..04d625a22 100644
--- a/src/io/adios2/adios2_writer.cpp
+++ b/src/io/adios2/adios2_writer.cpp
@@ -30,15 +30,21 @@ ADIOS2Writer<Concurrency>::ADIOS2Writer(adios2::ADIOS& adios, const Concurrency*
   }
 }
 
-// template <class CT>
-// ADIOS2Writer<CT>::ADIOS2Writer(const CT* concurrency, bool verbose)
-//     : adios_(GlobalAdios::getAdios()), verbose_(verbose), concurrency_(concurrency) {}
+template <class CT>
+ADIOS2Writer<CT>::ADIOS2Writer(const CT* concurrency, bool verbose)
+    : concurrency_(concurrency), adios_(*(concurrency->adios_)), verbose_(verbose) {
+  if constexpr (std::is_same<decltype(concurrency_), dca::parallel::MPIConcurrency>::value) {
+    std::cout << "AdiosWriter on MPI world:" << concurrency_->get_world_id()
+              << " of size: " << concurrency_->get_world_size() << '\n';
+    std::cout << "comm size: " << concurrency_->size() << " id: " << concurrency_->id() << '\n';
+  }
+}
 
 template <class Concurrency>
 ADIOS2Writer<Concurrency>::~ADIOS2Writer() {
   // close_file won't close an invalid engine.
   // but this could still be a problem if the adios2::ADIOS object has been destroyed.
-  if(file_)
+  if (file_)
     close_file();
 }
 
@@ -62,9 +68,10 @@ void ADIOS2Writer<Concurrency>::open_file(const std::string& file_name_ref, bool
   }
   io_name_ = file_name_ref;
   file_name_ = file_name_ref;
-  if(!io_)
-    io_ = adios_.DeclareIO(io_name_);
+  io_ = adios_.DeclareIO(io_name_);
+  io_.SetEngine("BP4");
 
+  
   file_ = io_.Open(file_name_, mode, concurrency_->get());
   // This is true if m_isClosed is false, that doesn't mean the "file" is open.
   if (!file_) {
@@ -77,12 +84,13 @@ void ADIOS2Writer<Concurrency>::open_file(const std::string& file_name_ref, bool
 template <class Concurrency>
 void ADIOS2Writer<Concurrency>::close_file() {
   if (static_cast<bool>(file_)) {
+    file_.Flush();
     file_.Close();
     adios_.RemoveIO(io_.Name());
-    //adios_.RemoveIO(io_name_);
-    //file_.Close();
-    // This combined with overwrite seems to create issues.
     // adios_.RemoveIO(io_name_);
+    // file_.Close();
+    //  This combined with overwrite seems to create issues.
+    //  adios_.RemoveIO(io_name_);
   }
 }
 
@@ -166,8 +174,8 @@ void ADIOS2Writer<Concurrency>::addAttribute(const std::string& set, const std::
 }
 
 template class ADIOS2Writer<dca::parallel::NoConcurrency>;
-
+#ifdef DCA_HAVE_MPI
 template class ADIOS2Writer<dca::parallel::MPIConcurrency>;
-
+#endif
 }  // namespace io
 }  // namespace dca
diff --git a/src/io/json/CMakeLists.txt b/src/io/json/CMakeLists.txt
index 337d7e7d2..8792c4fa0 100644
--- a/src/io/json/CMakeLists.txt
+++ b/src/io/json/CMakeLists.txt
@@ -5,6 +5,6 @@ if(DCA_HAVE_CUDA)
   target_link_libraries(json PUBLIC CUDA::cudart)
 endif()
 if(DCA_HAVE_HIP)
-  target_link_libraries(json PUBLIC hip::host)
+  target_link_libraries(json PUBLIC hip::host roc::hipblas roc::hipsparse)
 endif()
 		      
diff --git a/src/linalg/util/CMakeLists.txt b/src/linalg/util/CMakeLists.txt
index a30127299..519d30b16 100644
--- a/src/linalg/util/CMakeLists.txt
+++ b/src/linalg/util/CMakeLists.txt
@@ -5,6 +5,11 @@ if(DCA_HAVE_GPU)
   add_library(gpu_utils STATIC ${srcs})
   dca_gpu_runtime_link(gpu_utils)
   dca_gpu_blas_link(gpu_utils)
-  target_link_libraries(gpu_utils PUBLIC magma::magma) #"$<TARGET_PROPERTY:magma::magma,INTERFACE_LINK_LIBRARIES>")
+  target_link_libraries(gpu_utils PUBLIC magma::magma)
+  get_property(util_include_dirs TARGET gpu_utils PROPERTY INCLUDE_DIRECTORIES)
+  message("gpu_utils includes: ${util_include_dirs}")
+  get_property(interface_util_include_dirs TARGET gpu_utils PROPERTY INTERFACE_INCLUDE_DIRECTORIES)
+  message("gpu_utils interface includes: ${interface_util_include_dirs}")
+
 endif()
 
diff --git a/src/linalg/util/util_gpublas.cpp b/src/linalg/util/util_gpublas.cpp
index 8af51ff25..ea7a7b3b7 100644
--- a/src/linalg/util/util_gpublas.cpp
+++ b/src/linalg/util/util_gpublas.cpp
@@ -10,11 +10,11 @@
 //
 // This file implements cublas related utilities.
 
-#include "dca/linalg/util/handle_functions.hpp"
 #include "dca/linalg/util/util_gpublas.hpp"
+#include "dca/linalg/util/handle_functions.hpp"
 
-#include <mutex>
 #include <magma_v2.h>
+#include <mutex>
 
 
 namespace dca {
diff --git a/src/math/CMakeLists.txt b/src/math/CMakeLists.txt
index 4030f11b0..4dcf046c3 100644
--- a/src/math/CMakeLists.txt
+++ b/src/math/CMakeLists.txt
@@ -3,3 +3,12 @@ add_subdirectory(geometry)
 add_subdirectory(nfft)
 add_subdirectory(random)
 add_subdirectory(statistical_testing)
+
+add_library(dca_math INTERFACE)
+target_link_libraries(dca_math INTERFACE
+  function_transform
+  gaussian_quadrature
+  tetrahedron_mesh
+  nfft
+  random
+)
diff --git a/src/parallel/mpi_concurrency/mpi_concurrency.cpp b/src/parallel/mpi_concurrency/mpi_concurrency.cpp
index b3a98f2e9..d86b1713b 100644
--- a/src/parallel/mpi_concurrency/mpi_concurrency.cpp
+++ b/src/parallel/mpi_concurrency/mpi_concurrency.cpp
@@ -17,7 +17,12 @@
 namespace dca {
 namespace parallel {
 
-MPIConcurrency::MPIConcurrency(int argc, char** argv) : MPIInitializer(argc, argv) {
+  MPIConcurrency::MPIConcurrency(int argc, char** argv) : MPIInitializer(argc, argv)
+#ifdef DCA_HAVE_ADIOS2
+							,adios_(std::make_unique<adios2::ADIOS>(MPIProcessorGrouping::get()))
+#endif
+
+  {
   if (!MPIProcessorGrouping::isValid()) {
     std::string error_string;
     if (MPIProcessorGrouping::get_size() == MPIProcessorGrouping::get_world_size())
diff --git a/src/parallel/no_concurrency/no_concurrency.cpp b/src/parallel/no_concurrency/no_concurrency.cpp
index 3af8320cb..9ecd492b5 100644
--- a/src/parallel/no_concurrency/no_concurrency.cpp
+++ b/src/parallel/no_concurrency/no_concurrency.cpp
@@ -17,6 +17,12 @@
 namespace dca {
 namespace parallel {
 
+#ifdef DCA_HAVE_ADIOS2
+  NoConcurrency::NoConcurrency(int argc, char** argv) : adios_(std::make_unique<adios2::ADIOS>("", MPI_COMM_SELF)) {}
+#else
+  NoConcurrency::NoConcurrency(int argc, char** argv) {}
+#endif
+
 void NoConcurrency::abort() const {
   std::cout << "\nAborting process.\n";
   std::terminate();
diff --git a/test/dca_gtest_main_mpi.cpp b/test/dca_gtest_main_mpi.cpp
new file mode 100644
index 000000000..16ca98424
--- /dev/null
+++ b/test/dca_gtest_main_mpi.cpp
@@ -0,0 +1,31 @@
+
+
+#include "dca/testing/gtest_h_w_warning_blocking.h"
+#include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
+
+#ifdef DCA_HAVE_HPX
+#include <hpx/hpx_start.hpp>
+
+int hpx_main(int argc, char* argv[])
+{
+  dca::parallel::MPIConcurrency concurrency(argc, argv);
+  ::testing::InitGoogleTest(&argc, argv);
+  auto ret_gtest = RUN_ALL_TESTS();
+  hpx::finalize();
+  return ret_gtest;
+}
+
+int main(int argc, char** argv) {
+  std::cout << "doing HPX gtest main\n";
+  hpx::start(argc, argv);
+  // returns the value that hpx::main returns
+  return hpx::stop();
+}
+
+#else
+int main(int argc, char** argv) {
+  dca::parallel::MPIConcurrency concurrency(argc, argv);
+  ::testing::InitGoogleTest(&argc, argv);
+  return RUN_ALL_TESTS();
+}
+#endif
diff --git a/test/integration/cluster_solver/stdthread_qmci/stdthread_ctint_cluster_solver_io_test.cpp b/test/integration/cluster_solver/stdthread_qmci/stdthread_ctint_cluster_solver_io_test.cpp
index 7e61db898..ffe5db075 100644
--- a/test/integration/cluster_solver/stdthread_qmci/stdthread_ctint_cluster_solver_io_test.cpp
+++ b/test/integration/cluster_solver/stdthread_qmci/stdthread_ctint_cluster_solver_io_test.cpp
@@ -88,7 +88,7 @@ TEST(PosixCtintClusterSolverTest, PerMeasurementIO) {
   data.initialize();
 
   auto writer = std::make_shared<dca::io::Writer<TestConcurrency>>(
-      concurrency.get_adios(), std::ref(concurrency), parameters.get_output_format(), false);
+      std::ref(concurrency), parameters.get_output_format(), false);
   writer->open_file(parameters.get_filename_dca(), true);
   // Do one integration step.
   QmcSolver qmc_solver(parameters, data, writer);
diff --git a/test/integration/exact_diagonalization_advanced/CMakeLists.txt b/test/integration/exact_diagonalization_advanced/CMakeLists.txt
index fba8e75d9..f6d5269b9 100644
--- a/test/integration/exact_diagonalization_advanced/CMakeLists.txt
+++ b/test/integration/exact_diagonalization_advanced/CMakeLists.txt
@@ -1,21 +1,21 @@
 # Advanced exact diagonalization integration tests
 
 dca_add_gtest(hamiltonian_test
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${FFTW_INCLUDE_DIR};${PROJECT_SOURCE_DIR}
   LIBS function function_transform dca_io json time_and_frequency_domains cluster_domains enumerations quantum_domains ${LAPACK_LIBRARIES}
-       ${DCA_GPU_LIBS} ${UTIL_LIBS})
+       ${DCA_GPU_LIBS} ${UTIL_LIBS} parallel_no_concurrency)
 
 dca_add_gtest(ed_cluster_solver_single_site_test
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   CUDA
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
   LIBS function function_transform json time_and_frequency_domains cluster_domains enumerations quantum_domains dca_io timer
-       dca_algorithms ${LAPACK_LIBRARIES} ${DCA_THREADING_LIBS} ${DCA_GPU_LIBS} dca_hdf5 gpu_utils)
+       dca_algorithms ${LAPACK_LIBRARIES} ${DCA_THREADING_LIBS} ${DCA_GPU_LIBS} dca_hdf5 gpu_utils parallel_no_concurrency)
 
 dca_add_gtest(ed_cluster_solver_four_site_test
   EXTENSIVE
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
   LIBS function function_transform json time_and_frequency_domains cluster_domains enumerations quantum_domains dca_io timer
-       dca_algorithms ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS} ${DCA_THREADING_LIBS})
+       dca_algorithms ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS} ${DCA_THREADING_LIBS} parallel_no_concurrency)
diff --git a/test/integration/phys/lattice_mapping/CMakeLists.txt b/test/integration/phys/lattice_mapping/CMakeLists.txt
index 8a1abc1a8..df6f0c1a7 100644
--- a/test/integration/phys/lattice_mapping/CMakeLists.txt
+++ b/test/integration/phys/lattice_mapping/CMakeLists.txt
@@ -8,4 +8,4 @@ dca_add_gtest(lattice_mapping_sp_test EXTENSIVE
        time_and_frequency_domains
        quantum_domains gaussian_quadrature
        tetrahedron_mesh coarsegraining enumerations dca_hdf5 ${LAPACK_LIBRARIES} ${HDF5_LIBRARIES} ${DCA_GPU_LIBS}
-       lapack gpu_utils)
+       lapack gpu_utils parallel_no_concurrency)
diff --git a/test/system-level/dca/CMakeLists.txt b/test/system-level/dca/CMakeLists.txt
index b8647a241..26ee2a01f 100644
--- a/test/system-level/dca/CMakeLists.txt
+++ b/test/system-level/dca/CMakeLists.txt
@@ -2,26 +2,26 @@ dca_add_gtest(dca_sp_DCA+_mpi_test
   EXTENSIVE
   MPI MPI_NUMPROC 8
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS ${DCA_LIBS} ${DCA_KERNEL_LIBS})
+  LIBS dca_io ${DCA_LIBS} ${DCA_KERNEL_LIBS} dca_io)
 
 dca_add_gtest(dca_sp_DCA+_mpi_cuda_test
   EXTENSIVE
   MPI MPI_NUMPROC 8
   CUDA
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS ${DCA_LIBS} ${DCA_GPU_LIBS} ${DCA_KERNEL_LIBS})
+  LIBS ${DCA_LIBS} ${DCA_GPU_LIBS} ${DCA_KERNEL_LIBS} dca_io)
 
 dca_add_gtest(dca_sp_DCA+_thread_test
   GTEST_MAIN
   EXTENSIVE
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS ${DCA_LIBS} ${DCA_KERNEL_LIBS})
+  LIBS dca_io ${DCA_LIBS} ${DCA_KERNEL_LIBS} dca_io)
 
 dca_add_gtest(ctint_sp_DCA_mpi_thread_test
   MPI MPI_NUMPROC 2
   EXTENSIVE
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS ${DCA_LIBS} ${DCA_KERNEL_LIBS})
+  LIBS ${DCA_LIBS} ${DCA_KERNEL_LIBS} dca_io)
 
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/initial_sigma.hdf5
                ${CMAKE_CURRENT_BINARY_DIR}/initial_sigma.hdf5 COPYONLY)
diff --git a/test/unit/io/CMakeLists.txt b/test/unit/io/CMakeLists.txt
index 0a282dbdb..8783cdb58 100644
--- a/test/unit/io/CMakeLists.txt
+++ b/test/unit/io/CMakeLists.txt
@@ -4,11 +4,12 @@ dca_add_gtest(buffer_test
   GTEST_MAIN)
 
 dca_add_gtest(reader_test
-  LIBS function dca_io)
+  MPI
+  LIBS function ${DCA_LIBS} dca_io)
 
 dca_add_gtest(writer_test
   MPI
-  LIBS function dca_io parallel_mpi_concurrency)
+  LIBS function dca_io parallel_mpi_concurrency parallel_no_concurrency)
 
 dca_add_gtest(json_reader_test
   GTEST_MAIN
diff --git a/test/unit/io/adios2_parallelIO_test.cpp b/test/unit/io/adios2_parallelIO_test.cpp
index 243256870..00ec6bb5f 100644
--- a/test/unit/io/adios2_parallelIO_test.cpp
+++ b/test/unit/io/adios2_parallelIO_test.cpp
@@ -27,8 +27,23 @@
 #include "gtest/gtest.h"
 
 int rank, comm_size;
-dca::parallel::MPIConcurrency* concurrency_ptr;
-adios2::ADIOS* adios_ptr;
+
+#ifdef DCA_HAVE_MPI
+#include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
+using Concurrency = dca::parallel::MPIConcurrency;
+#else
+#error "This file needs GPU support."
+#endif
+
+Concurrency* concurrency_ptr = nullptr;
+
+#ifdef DCA_HAVE_ADIOS2
+#include "dca/io/adios2/adios2_writer.hpp"
+#include "dca/io/adios2/adios2_global.hpp"
+adios2::ADIOS* adios2_ptr = nullptr;
+#else
+#error "This file needs ADIOS2 support."
+#endif
 
 template <typename Scalar>
 class ADIOS2ParallelIOTest : public ::testing::Test {};
@@ -101,16 +116,17 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWrite) {
 
   const std::string fname("ADIOS2ParallelIOTest_" + typeStr + ".bp");
   {
-    dca::io::ADIOS2Writer<dca::parallel::MPIConcurrency> writer(*adios_ptr, concurrency_ptr, false);
+    dca::io::ADIOS2Writer<dca::parallel::MPIConcurrency> writer(concurrency_ptr);
     writer.open_file(fname, true);
-
+    writer.begin_step();
     // Because the caller needs to know if its function is distributed or not we will assume this is
     // so for the API as well. in the future I think something more sophisticated needs to be done
     // and the function will need to know its distribution, but for now we distribute only over the
     // fastest index
 
     writer.execute(f1, subind_start, subind_end);
-    adios_ptr->FlushAll();
+    //adios_ptr->FlushAll();
+    writer.end_step();
     MPI_Barrier(concurrency_ptr->get());
     writer.close_file();
   }
@@ -120,11 +136,10 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWrite) {
       std::cout << " Read back data with 3D selection " << std::endl;
     }
 
-    dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, false);
+    dca::io::ADIOS2Reader reader(*concurrency_ptr);
     reader.open_file(fname);
-
+    reader.begin_step();
     dca::func::function<Scalar, Dmn, dca::DistType::BLOCKED> f2("parallelFunc", *concurrency_ptr);
-
     EXPECT_TRUE(reader.execute(f2, subind_start, subind_end));
 
     /* TODO: This should be working on every rank */
@@ -146,7 +161,7 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWrite) {
       if (f1(i) != f2(i))
         break;
     }
-
+    reader.end_step();
     reader.close_file();
   }
 }
@@ -193,11 +208,11 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionImplicitBlockingReadWrite) {
 
   const std::string fname("ADIOS2ParallelIOTest_ImplicitSubIndices_" + typeStr + ".bp");
   {
-    dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr, true);
+    dca::io::ADIOS2Writer writer(concurrency_ptr, true);
     writer.open_file(fname, true);
-
+    writer.begin_step();
     writer.execute(f1);
-
+    writer.end_step();
     MPI_Barrier(concurrency_ptr->get());
     writer.close_file();
   }
@@ -207,11 +222,11 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionImplicitBlockingReadWrite) {
       std::cout << " Read back data with 3D selection from " << fname << std::endl;
     }
 
-    dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, true);
+    dca::io::ADIOS2Reader reader(*concurrency_ptr, true);
     reader.open_file(fname);
-
+    reader.begin_step();
     dca::func::function<Scalar, Dmn, dca::DistType::BLOCKED> f2("parallelFunc", *concurrency_ptr);
-
+    MPI_Barrier(concurrency_ptr->get());
     EXPECT_TRUE(reader.execute(f2));
 
     for (int i = 0; i < f2.get_end() - f2.get_start() + 1; ++i) {
@@ -219,7 +234,7 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionImplicitBlockingReadWrite) {
       if (f1(i) != f2(i))
         break;
     }
-
+    reader.end_step();
     reader.close_file();
   }
 }
@@ -256,15 +271,16 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWriteLinear) {
 
   const std::string fname("ADIOS2ParallelIOTest_linear_" + typeStr + ".bp");
   {
-    dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr, true);
+    dca::io::ADIOS2Writer writer(concurrency_ptr, true);
     writer.open_file(fname, true);
-
+    writer.begin_step();
     // Because the caller needs to know if its function is distributed or not we will assume this is
     // so for the API as well. in the future I think something more sophisticated needs to be done
     // and the function will need to know its distribution, but for now we distribute only over the
     // fastest index
 
     writer.execute(f1, start, end);
+    writer.end_step();
     MPI_Barrier(concurrency_ptr->get());
     writer.close_file();
   }
@@ -273,11 +289,11 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWriteLinear) {
     if (!rank) {
       std::cout << " Read back data with linear 1D selection " << std::endl;
     }
-    dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, true);
+    dca::io::ADIOS2Reader reader(*concurrency_ptr, true);
     reader.open_file(fname);
-
+    reader.begin_step();
     dca::func::function<Scalar, Dmn, dca::DistType::LINEAR> f2("parallelFuncLin", *concurrency_ptr);
-
+    MPI_Barrier(concurrency_ptr->get());
     EXPECT_TRUE(reader.execute(f2, start, end));
 
     /* TODO: This should be working on every rank */
@@ -299,7 +315,7 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWriteLinear) {
 
     /* 1D array on disk cannot be read back with 3D selection */
     EXPECT_FALSE(reader.execute(f3, subind_start, subind_end));
-
+    reader.end_step();
     reader.close_file();
   }
 }
@@ -332,15 +348,16 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWriteLinearIrregular) {
 
   const std::string fname("ADIOS2ParallelIOTest_linear_" + typeStr + ".bp");
   {
-    dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr, true);
+    dca::io::ADIOS2Writer writer(concurrency_ptr, true);
     writer.open_file(fname, true);
-
+    writer.begin_step();
     // Because the caller needs to know if its function is distributed or not we will assume this is
     // so for the API as well. in the future I think something more sophisticated needs to be done
     // and the function will need to know its distribution, but for now we distribute only over the
     // fastest index
 
     writer.execute(f1, start, end);
+    writer.end_step();
     MPI_Barrier(concurrency_ptr->get());
     writer.close_file();
   }
@@ -349,11 +366,11 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWriteLinearIrregular) {
     if (!rank) {
       std::cout << " Read back data with linear 1D selection " << std::endl;
     }
-    dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, true);
+    dca::io::ADIOS2Reader reader(*concurrency_ptr, true);
     reader.open_file(fname);
-
+    reader.begin_step();
     dca::func::function<Scalar, Dmn, dca::DistType::LINEAR> f2("parallelFuncLin", *concurrency_ptr);
-
+    MPI_Barrier(concurrency_ptr->get());
     EXPECT_TRUE(reader.execute(f2, start, end));
 
     /* TODO: This should be working on every rank */
@@ -382,24 +399,26 @@ TYPED_TEST(ADIOS2ParallelIOTest, FunctionReadWriteLinearIrregular) {
 }
 
 int main(int argc, char** argv) {
-  int result = 0;
+  // This results in a copy constructor beging called at somepoint,  resulting in an MPI_INIT after
+  // the finalize. concurrency = std::make_unique<dca::parallel::MPIConcurrency>(argc, argv);
+  concurrency_ptr = new Concurrency(argc, argv);
 
-  dca::parallel::MPIConcurrency concurrency(argc, argv);
-  rank = concurrency.id();
-  comm_size = concurrency.number_of_processors();
-  concurrency_ptr = &concurrency;
-  ::testing::InitGoogleTest(&argc, argv);
+  rank = concurrency_ptr->id();
+  std::cout << "rank: " << rank << '\n';
+  comm_size = concurrency_ptr->number_of_processors();
 
+  ::testing::InitGoogleTest(&argc, argv);
   ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
-  if (rank != 0) {
+
+  if (concurrency_ptr->id() != 0) {
     delete listeners.Release(listeners.default_result_printer());
     listeners.Append(new dca::testing::MinimalistPrinter);
   }
 
-  adios2::ADIOS adios("", concurrency_ptr->get());
-  adios_ptr = &adios;
+  int result = RUN_ALL_TESTS();
 
-  result = RUN_ALL_TESTS();
+  delete concurrency_ptr;
 
   return result;
-}
+
+  }
diff --git a/test/unit/io/adios2_reader_writer_test.cpp b/test/unit/io/adios2_reader_writer_test.cpp
index 3eb08d98c..0f2766e45 100644
--- a/test/unit/io/adios2_reader_writer_test.cpp
+++ b/test/unit/io/adios2_reader_writer_test.cpp
@@ -41,7 +41,7 @@ TEST(ADIOS2ReaderWriterTest, ReaderDestructorCleanUp) {
   std::string object_name = "forty-two";
 
   // Create test file.
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   const int i = 42;
 
   writer.open_file(test_file_name);
@@ -51,7 +51,7 @@ TEST(ADIOS2ReaderWriterTest, ReaderDestructorCleanUp) {
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   int j;
 
   reader.open_file(test_file_name);
@@ -70,7 +70,7 @@ TEST(ADIOS2ReaderWriterTest, WriterDestructorCleanUp) {
   std::string group_name_2 = "magic-numbers";
   std::string object_name = "forty-two";
 
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);;
+  dca::io::ADIOS2Writer writer(concurrency_ptr);;
   const int i = 42;
 
   writer.open_file(test_file_name);
@@ -91,15 +91,18 @@ TEST(ADIOS2ReaderWriterTest, VectorReadWrite) {
       std::complex<double>(1., 0.), std::complex<double>(0., 1.), std::complex<double>(23.4, -1.5)};
 
   // Create test file.
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   writer.open_file(file_name);
+  writer.begin_step();
   writer.execute(object_name, a_vector);
+  writer.end_step();
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   std::vector<std::complex<double>> vector_read;
   reader.open_file(file_name);
+  reader.begin_step();
   EXPECT_TRUE(reader.execute(object_name, vector_read));
 
   ASSERT_EQ(a_vector.size(), vector_read.size());
@@ -107,28 +110,31 @@ TEST(ADIOS2ReaderWriterTest, VectorReadWrite) {
     EXPECT_DOUBLE_EQ(std::real(a_vector[i]), std::real(vector_read[i]));
     EXPECT_DOUBLE_EQ(std::imag(a_vector[i]), std::imag(vector_read[i]));
   }
-
+  reader.end_step();
   reader.close_file();
 }
 
 TEST(ADIOS2ReaderWriterTest, VectorOfVectorsReadWrite) {
   const std::string object_name = "a_vector";
-  const std::string file_name = "test_vector_vector.bp";
+  const std::string file_name = "adios2_test_vector_vector.bp";
   const std::vector<std::vector<double>> data_unequal_size{{1, 0, 2}, {1}, {1, 0}, {}, {2, 2}};
 
   // Create test file.
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr, true);
+  dca::io::ADIOS2Writer writer(concurrency_ptr, true);
   writer.open_file(file_name);
+  writer.begin_step();
   writer.execute(object_name, data_unequal_size);
+  writer.end_step();
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, true);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr, true);
   std::vector<std::vector<double>> data_read;
   reader.open_file(file_name);
+  reader.begin_step();
   EXPECT_TRUE(reader.execute(object_name, data_read));
-
   EXPECT_EQ(data_unequal_size, data_read);
+  reader.end_step();
   reader.close_file();
 }
 
@@ -139,19 +145,22 @@ TEST(ADIOS2ReaderWriterTest, VectorOfArraysReadWrite) {
   std::vector<std::array<int, 3>> data{{-1, 2, 3}, {5, -7, 0}, {9, 2, -1}};
 
   // Create test file.
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr, true);
+  dca::io::ADIOS2Writer writer(concurrency_ptr, true);
   writer.open_file(file_name);
+  writer.begin_step();
   writer.execute(object_name, data);
+  writer.end_step();
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, true);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr, true);
   std::vector<std::array<int, 3>> data_read;
   reader.open_file(file_name);
+  reader.begin_step();
   EXPECT_TRUE(reader.execute(object_name, data_read));
 
   EXPECT_EQ(data, data_read);
-
+  reader.end_step();
   reader.close_file();
 }
 
@@ -161,22 +170,25 @@ TEST(ADIOS2ReaderWriterTest, StringAndVectorOfStringsReadWrite) {
   const std::string filename = "test_vec_of_strings.bp";
 
   // Create test file.
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   writer.open_file(filename);
+  writer.begin_step();
   writer.execute("single-string", s1);
   writer.execute("strings", s_vec1);
+  writer.end_step();
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   reader.open_file(filename);
-  //
+  reader.begin_step();
   std::vector<std::string> s_vec2;
   std::string s2;
   EXPECT_TRUE(reader.execute("strings", s_vec2));
   EXPECT_TRUE(reader.execute("single-string", s2));
   EXPECT_EQ(s_vec1, s_vec2);
   EXPECT_EQ(s1, s2);
+  reader.end_step();
   reader.close_file();
 }
 
@@ -197,24 +209,24 @@ TYPED_TEST(ADIOS2ReaderWriterTestTyped, FunctionReadWrite) {
   for (auto& x : f1)
     x = ++val;
 
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   writer.open_file("test_func.bp", true);
-
+  writer.begin_step();
   writer.execute(f1);
-
+  writer.end_step();
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   reader.open_file("test_func.bp");
-
+  reader.begin_step();
   dca::func::function<Scalar, Dmn> f2("myfunc");
 
   EXPECT_TRUE(reader.execute(f2));
 
   for (int i = 0; i < f1.size(); ++i)
     EXPECT_EQ(f1(i), f2(i));
-
+  reader.end_step();
   reader.close_file();
 }
 
@@ -228,13 +240,13 @@ TYPED_TEST(ADIOS2ReaderWriterTestTyped, MatrixReadWrite) {
     for (int i = 0; i < m1.nrRows(); ++i)
       m1(i, j) = ++val;
 
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   writer.open_file("test_mat.bp", true);
   writer.execute(m1);
   writer.close_file();
 
   // Read test file.
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   reader.open_file("test_mat.bp");
 
   dca::linalg::Matrix<Scalar, dca::linalg::CPU> m2;
@@ -254,7 +266,7 @@ TEST(ADIOS2ReaderWriterTest, NonAccessibleFile) {
   //dca::io::ADIOS2Writer writer;
   //EXPECT_THROW(writer.open_file("/not_accessible_directory/file.txt"), std::ios_base::failure);
 
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   // Exception type doesn't seem portable
   EXPECT_ANY_THROW(reader.open_file("not_existing_file.txt"));
 }
@@ -265,7 +277,7 @@ TEST(ADIOS2ReaderWriterTest, FunctionNotPresent) {
   dca::func::function<int, Dmn> present("present");
   present = 1;
 
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   writer.open_file("adios2_missing_func.bp");
   writer.execute(present);
   writer.close_file();
@@ -273,7 +285,7 @@ TEST(ADIOS2ReaderWriterTest, FunctionNotPresent) {
   dca::func::function<int, Dmn> not_present("not_present");
   present = 0;
 
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   reader.open_file("adios2_missing_func.bp");
   EXPECT_FALSE(reader.execute(not_present));
   EXPECT_TRUE(reader.execute(present));
@@ -285,7 +297,7 @@ TEST(ADIOS2ReaderWriterTest, FunctionNotPresent) {
 }
 
 TEST(ADIOS2ReaderWriterTest, GroupOpenclose) {
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
   writer.open_file("group_open_close.bp");
 
   writer.open_group("foo");
@@ -302,7 +314,7 @@ TEST(ADIOS2ReaderWriterTest, GroupOpenclose) {
 
   writer.close_file();
 
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);
   reader.open_file("group_open_close.bp");
 
   int i_val;
@@ -324,21 +336,21 @@ TEST(ADIOS2ReaderWriterTest, GroupOpenclose) {
 TEST(ADIOS2ReaderWriterTest, VectorOfVector) {
   dca::io::ADIOS2Writer writer(GlobalAdios::getAdios(), concurrency_ptr, true);
   writer.open_file("vector.bp");
-
+  writer.begin_step();
   // Simple 3D vector
   std::vector<double> vec_1{1., 2., 3.};
   writer.execute("simple-vector", vec_1);
 
   std::vector<std::vector<float>> vec_2{{1.2, 3.4}, {5.6, 7.8, 4.4}, {1.0}};
   writer.execute("vector-of-vectors", vec_2);
-
+  writer.end_step();
   writer.close_file();
 }
 
 /** Semantics of adios_writer have changed to allow overwrite
  */
 // TEST(ADIOS2ReaderWriterTest, Overwrite) {
-//   dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
+//   dca::io::ADIOS2Writer writer(concurrency_ptr);
 //   writer.open_file("test.bp", true);
 
 //   writer.open_group("foo");
@@ -350,7 +362,7 @@ TEST(ADIOS2ReaderWriterTest, VectorOfVector) {
 
 //   writer.close_file();
 
-//   dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);
+//   dca::io::ADIOS2Reader reader(concurrency_ptr);
 //   reader.open_file("test.bp");
 
 //   int i_val;
diff --git a/test/unit/io/reader_test.cpp b/test/unit/io/reader_test.cpp
index 2815c9a7a..0532ed656 100644
--- a/test/unit/io/reader_test.cpp
+++ b/test/unit/io/reader_test.cpp
@@ -13,9 +13,9 @@
 #include "dca/io/csv/csv_reader.hpp"
 #include "dca/io/hdf5/hdf5_reader.hpp"
 #include "dca/io/json/json_reader.hpp"
+#include <functional>
 #ifdef DCA_HAVE_ADIOS2
 #include "dca/io/adios2/adios2_reader.hpp"
-adios2::ADIOS* adios_ptr;
 #endif
 #include <vector>
 #include "gtest/gtest.h"
@@ -26,7 +26,14 @@ template <typename T>
 class ReaderTest : public ::testing::Test {};
 
 
+#ifdef DCA_HAVE_MPI
+using Concurrency = dca::parallel::MPIConcurrency;
+#else
 using Concurrency = dca::parallel::NoConcurrency;
+#endif
+
+
+
 
 using ReaderTypes = ::testing::Types<dca::io::CSVReader, dca::io::HDF5Reader, dca::io::JSONReader
 #ifdef DCA_HAVE_ADIOS2
@@ -36,7 +43,7 @@ using ReaderTypes = ::testing::Types<dca::io::CSVReader, dca::io::HDF5Reader, dc
 
 TYPED_TEST_CASE(ReaderTest, ReaderTypes);
 
-dca::parallel::NoConcurrency* concurrency_ptr;
+std::unique_ptr<Concurrency> concurrency_ptr;
 
 TEST(JSONReaderTest, Vector) {
   dca::io::JSONReader reader;
@@ -81,7 +88,7 @@ TEST(HDF5ReaderTest, Vector) {
 #ifdef DCA_HAVE_ADIOS2
 // more extensive testing in adios2_reader_writer_test
 TEST(ADIOS2ReaderTest, Vector) {
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr, true);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr, true);
   reader.open_file(DCA_SOURCE_DIR "/test/unit/io/vector.bp");
 
   // Simple 3D vector
@@ -103,15 +110,7 @@ TEST(ADIOS2ReaderTest, Vector) {
 #endif
 
 int main(int argc, char** argv) {
-  dca::parallel::NoConcurrency concurrency(argc, argv);
-  concurrency_ptr = &concurrency;
-
-#ifdef DCA_HAVE_ADIOS2
-  //ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", false);
-  adios_ptr = &adios;
-#endif
-
+  concurrency_ptr = std::make_unique<Concurrency>(argc, argv);
   ::testing::InitGoogleTest(&argc, argv);
 
   // ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
diff --git a/test/unit/io/reader_writer_test.cpp b/test/unit/io/reader_writer_test.cpp
index 7df32bde5..2f2c3edb6 100644
--- a/test/unit/io/reader_writer_test.cpp
+++ b/test/unit/io/reader_writer_test.cpp
@@ -16,12 +16,14 @@
 
 #ifdef DCA_HAVE_MPI
 #include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
-dca::parallel::MPIConcurrency* concurrency_ptr = nullptr;
+using Concurrency = dca::parallel::MPIConcurrency;
 #else
 #include "dca/parallel/no_concurrency/no_concurrency.hpp"
-dca::parallel::NoConcurrency* concurrency_ptr;
+using Concurrency = parallel::NoConcurrency;
 #endif
 
+Concurrency* concurrency_ptr = nullptr;
+
 #ifdef DCA_HAVE_ADIOS2
 #include "dca/io/adios2/adios2_writer.hpp"
 #include "dca/io/adios2/adios2_global.hpp"
@@ -63,11 +65,7 @@ TEST(ReaderWriterTest, RAIIDestructor) {
 
     const int i = 42;
     {
-#ifdef DCA_HAVE_ADIOS2
-      dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type, true);
-#else
       dca::io::Writer writer(*concurrency_ptr, type);
-#endif
       writer.open_file(test_file_name);
       writer.open_group(group_name_1);
       writer.open_group(group_name_2);
@@ -93,19 +91,18 @@ TEST(ReaderWriterTest, VectorReadWrite) {
         std::complex<double>(1., 0.), std::complex<double>(0., 1.), std::complex<double>(23.4, -1.5)};
 
     // Create test file.
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file(file_name);
+    writer.begin_step();
     writer.execute(object_name, a_vector);
+    writer.end_step();
     writer.close_file();
 
     // Read test file.
     dca::io::Reader reader(*concurrency_ptr, type);
     std::vector<std::complex<double>> vector_read;
     reader.open_file(file_name);
+    reader.begin_step();
     EXPECT_TRUE(reader.execute(object_name, vector_read));
 
     ASSERT_EQ(a_vector.size(), vector_read.size());
@@ -125,49 +122,49 @@ TEST(ReaderWriterTest, VectorOfVectorsReadWrite) {
     const std::vector<std::vector<double>> data_unequal_size{{0, 0, 2}, {1}, {1, 0}, {}, {0, 0}};
 
     // Create test file.
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file(file_name);
+    writer.begin_step();
     writer.execute(object_name, data_unequal_size);
+    writer.end_step();
     writer.close_file();
 
     // Read test file.
     dca::io::Reader reader(*concurrency_ptr, type);
     std::vector<std::vector<double>> data_read;
     reader.open_file(file_name);
+    reader.begin_step();
     EXPECT_TRUE(reader.execute(object_name, data_read));
 
     EXPECT_EQ(data_unequal_size, data_read);
+    reader.end_step();
     reader.close_file();
   }
 }
 
 TEST(ReaderWriterTest, VectorOfArraysReadWrite) {
   for (auto type : types) {
-    const std::string object_name = "obj";
+    const std::string object_name = "obj_11";
     const std::string file_name = "test_vec_of_arr." + toLower(type);
 
     std::vector<std::array<int, 3>> data{{-1, 2, 3}, {5, -7, 0}};
 
-// Create test file.
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
+    // Create test file.
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file(file_name);
-    writer.execute(object_name, data);
+    writer.begin_step();
+    writer.execute("test_ovj", data);
+    writer.end_step();
     writer.close_file();
 
     // Read test file.
     dca::io::Reader reader(*concurrency_ptr, type);
-    std::vector<std::array<int, 3>> data_read;
+    std::vector<std::array<int, 3>> data_read(2);
     reader.open_file(file_name);
-    EXPECT_TRUE(reader.execute(object_name, data_read));
+    reader.begin_step();
+    EXPECT_TRUE(reader.execute("test_ovj", data_read));
     EXPECT_EQ(data, data_read);
+    reader.end_step();
     reader.close_file();
   }
 }
@@ -178,12 +175,8 @@ TEST(ReaderWriterTest, StringAndVectorOfStringsReadWrite) {
     std::string s1{"bazinga"};
     const std::string filename = "test_vec_of_strings" + toLower(type);
 
-// Create test file.
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
+    // Create test file.
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file(filename);
     writer.execute("single-string", s1);
     writer.execute("strings", s_vec1);
@@ -221,11 +214,7 @@ TYPED_TEST(ReaderWriterTest, FunctionReadWrite) {
     for (auto& x : f1)
       x = ++val;
     {
-#ifdef DCA_HAVE_ADIOS2
-      dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
       dca::io::Writer writer(*concurrency_ptr, type);
-#endif
       writer.open_file("test_func." + toLower(type), true);
       writer.execute(f1);
     }
@@ -258,11 +247,7 @@ TYPED_TEST(ReaderWriterTest, MatrixReadWrite) {
         m1(i, j) = ++val;
 
     {
-#ifdef DCA_HAVE_ADIOS2
-      dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
       dca::io::Writer writer(*concurrency_ptr, type);
-#endif
       writer.open_file(filename, true);
       EXPECT_TRUE(writer.execute(m1));
       writer.close_file();
@@ -282,21 +267,13 @@ TYPED_TEST(ReaderWriterTest, MatrixReadWrite) {
 
 TEST(ReaderWriterTest, NonAccessibleFile) {
   for (auto type : types) {
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     H5::Exception::dontPrint();
     // Exception type doesn't seem portable
     if (type != "ADIOS2")
       EXPECT_ANY_THROW(writer.open_file("not_existing_directory/file.txt"));
 
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Reader reader(*concurrency_ptr, type);
-#else
     dca::io::Reader reader(*concurrency_ptr, type);
-#endif
     // Exception type doesn't seem portable
     EXPECT_ANY_THROW(reader.open_file("not_existing_file.txt"));
   }
@@ -308,11 +285,7 @@ TEST(ReaderWriterTest, FunctionNotPresent) {
     dca::func::function<int, Dmn> present("present");
     present = 1;
 
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file("missing_func." + toLower(type));
     writer.execute(present);
     writer.close_file();
@@ -320,11 +293,7 @@ TEST(ReaderWriterTest, FunctionNotPresent) {
     dca::func::function<int, Dmn> not_present("not_present");
     present = 0;
 
-#ifdef DCA_HAVE_ADIOS2
     dca::io::Reader reader(*concurrency_ptr, type);
-#else
-    dca::io::Reader reader(*concurrency_ptr, type);
-#endif
     reader.open_file("missing_func." + toLower(type));
     EXPECT_FALSE(reader.execute(not_present));
     EXPECT_TRUE(reader.execute(present));
@@ -338,11 +307,7 @@ TEST(ReaderWriterTest, FunctionNotPresent) {
 
 TEST(ReaderWriterTest, GroupOpenclose) {
   for (auto type : types) {
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file("group_open_close" + toLower(type));
 
     writer.open_group("foo");
@@ -359,11 +324,7 @@ TEST(ReaderWriterTest, GroupOpenclose) {
 
     writer.close_file();
 
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Reader reader(*concurrency_ptr, type);
-#else
     dca::io::Reader reader(*concurrency_ptr, type);
-#endif
     reader.open_file("group_open_close" + toLower(type));
 
     int i_val;
@@ -385,11 +346,7 @@ TEST(ReaderWriterTest, GroupOpenclose) {
 
 TEST(ReaderWriterTest, SteppedGroupOpenclose) {
   for (auto type : stepped_types) {
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     std::string test_filename{"group_open_close_stepped." + toLower(type)};
     writer.open_file(test_filename);
 
@@ -423,11 +380,7 @@ TEST(ReaderWriterTest, SteppedGroupOpenclose) {
 
     writer.close_file();
 
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Reader reader(*concurrency_ptr, type);
-#else
     dca::io::Reader reader(*concurrency_ptr, type);
-#endif
     reader.open_file(test_filename);
 
     int i_val;
@@ -487,17 +440,12 @@ TEST(ReaderWriterTest, SteppedGroupOpenclose) {
     reader.open_group("bar");
     reader.execute("b2", d_val);
     EXPECT_EQ(1.5, d_val);
-
   }
 }
 
 TEST(ReaderWriterTest, Overwrite) {
   for (auto type : types) {
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
     dca::io::Writer writer(*concurrency_ptr, type);
-#endif
     writer.open_file("test" + toLower(type), true);
 
     writer.open_group("foo");
@@ -509,11 +457,7 @@ TEST(ReaderWriterTest, Overwrite) {
 
     writer.close_file();
 
-#ifdef DCA_HAVE_ADIOS2
-    dca::io::Reader reader(*concurrency_ptr, type);
-#else
     dca::io::Reader reader(*concurrency_ptr, type);
-#endif
     reader.open_file("test" + toLower(type));
 
     int i_val;
@@ -527,11 +471,7 @@ TEST(ReaderWriterTest, Overwrite) {
 TEST(ReaderWriterTest, CompoundName) {
   std::string type{"HDF5"};
 
-#ifdef DCA_HAVE_ADIOS2
-  dca::io::Writer writer(*adios2_ptr, *concurrency_ptr, type);
-#else
   dca::io::Writer writer(*concurrency_ptr, type);
-#endif
   writer.open_file("compound_names." + toLower(type), true);
 
   // compound names are actually only supported by the reader for hdf5
@@ -547,11 +487,7 @@ TEST(ReaderWriterTest, CompoundName) {
 
   writer.close_file();
 
-#ifdef DCA_HAVE_ADIOS2
-  dca::io::Reader reader(*concurrency_ptr, type);
-#else
   dca::io::Reader reader(*concurrency_ptr, type);
-#endif
   reader.open_file("compound_names." + toLower(type));
 
   int i_val;
@@ -567,13 +503,7 @@ TEST(ReaderWriterTest, CompoundName) {
 int main(int argc, char** argv) {
   // This results in a copy constructor beging called at somepoint,  resulting in an MPI_INIT after
   // the finalize. concurrency = std::make_unique<dca::parallel::MPIConcurrency>(argc, argv);
-#ifdef DCA_HAVE_MPI
-  concurrency_ptr = new dca::parallel::MPIConcurrency(argc, argv);
-#ifdef DCA_HAVE_ADIOS2
-  adios2::ADIOS adios("", concurrency_ptr->get());
-  adios2_ptr = &(concurrency_ptr->get_adios());
-#endif
-#endif
+  concurrency_ptr = new Concurrency(argc, argv);
 
   ::testing::InitGoogleTest(&argc, argv);
   ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
diff --git a/test/unit/io/writer_test.cpp b/test/unit/io/writer_test.cpp
index 1851b36d2..0be5b1c8b 100644
--- a/test/unit/io/writer_test.cpp
+++ b/test/unit/io/writer_test.cpp
@@ -1,5 +1,5 @@
-// Copyright (C) 2022 ETH Zurich
-// Copyright (C) 2022 UT-Battelle, LLC
+// Copyright (C) 2024 ETH Zurich
+// Copyright (C) 2024 UT-Battelle, LLC
 // All rights reserved.
 //
 // See LICENSE for terms of usage.
@@ -22,25 +22,26 @@
 #include "dca/util/ignore.hpp"
 #ifdef DCA_HAVE_MPI
 #include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
-dca::parallel::MPIConcurrency* mpi_concurrency_ptr = nullptr;
 #ifdef DCA_HAVE_ADIOS2
 #include "dca/io/adios2/adios2_writer.hpp"
 adios2::ADIOS* adios_ptr = nullptr;
 #endif
-#endif
+using Concurrency = dca::parallel::MPIConcurrency;
+#else
 #include "dca/parallel/no_concurrency/no_concurrency.hpp"
-dca::parallel::NoConcurrency* concurrency_ptr = nullptr;
+using Concurrency = dca::parallel::NoCurrency;
+#endif
+
+Concurrency* concurrency_ptr = nullptr;
 
 template <typename T>
 class WriterTest : public ::testing::Test {};
 
-using Concurrency = dca::parallel::NoConcurrency;
-using MPIConcurrency = dca::parallel::MPIConcurrency;
 
 using WriterTypes = ::testing::Types<dca::io::HDF5Writer, dca::io::JSONWriter
 #ifdef DCA_HAVE_ADIOS2
                                      ,
-                                     dca::io::ADIOS2Writer<MPIConcurrency>
+                                     dca::io::ADIOS2Writer<Concurrency>
 #endif
                                      >;  //, dca::io::CSVWriter>;
 TYPED_TEST_CASE(WriterTest, WriterTypes);
@@ -48,8 +49,8 @@ TYPED_TEST_CASE(WriterTest, WriterTypes);
 TYPED_TEST(WriterTest, Constructor) {
   std::unique_ptr<TypeParam> writer_ptr;
 #ifdef DCA_HAVE_ADIOS2
-  if constexpr (std::is_same<TypeParam, dca::io::ADIOS2Writer<MPIConcurrency>>::value)
-    writer_ptr = std::make_unique<TypeParam>(*adios_ptr, mpi_concurrency_ptr);
+  if constexpr (std::is_same<TypeParam, dca::io::ADIOS2Writer<Concurrency>>::value)
+    writer_ptr = std::make_unique<TypeParam>(concurrency_ptr);
   else
 #endif
     writer_ptr = std::make_unique<TypeParam>();
@@ -59,8 +60,8 @@ TYPED_TEST(WriterTest, Constructor) {
 TYPED_TEST(WriterTest, Unique_Ptr) {
   std::unique_ptr<TypeParam> writer_ptr;
 #ifdef DCA_HAVE_ADIOS2
-  if constexpr (std::is_same<TypeParam, dca::io::ADIOS2Writer<MPIConcurrency>>::value)
-    writer_ptr = std::make_unique<TypeParam>(*adios_ptr, mpi_concurrency_ptr);
+  if constexpr (std::is_same<TypeParam, dca::io::ADIOS2Writer<Concurrency>>::value)
+    writer_ptr = std::make_unique<TypeParam>(concurrency_ptr);
   else
 #endif
     writer_ptr = std::make_unique<TypeParam>();
@@ -92,17 +93,11 @@ TYPED_TEST(WriterTest, Unique_Ptr) {
 }
 
 int main(int argc, char** argv) {
-  dca::parallel::NoConcurrency concurrency(argc, argv);
-
+  Concurrency concurrency(argc, argv);
   concurrency_ptr = &concurrency;
 
-#ifdef DCA_HAVE_MPI
-  dca::parallel::MPIConcurrency mpi_concurrency(argc, argv);
-  mpi_concurrency_ptr = &mpi_concurrency;
 #ifdef DCA_HAVE_ADIOS2
-  adios2::ADIOS adios("", mpi_concurrency_ptr->get());
-  adios_ptr = &adios;
-#endif
+  adios_ptr = &(concurrency_ptr->get_adios());
 #endif
 
   ::testing::InitGoogleTest(&argc, argv);
diff --git a/test/unit/linalg/CMakeLists.txt b/test/unit/linalg/CMakeLists.txt
index 81cfdf190..4d13dae2d 100644
--- a/test/unit/linalg/CMakeLists.txt
+++ b/test/unit/linalg/CMakeLists.txt
@@ -22,9 +22,8 @@ dca_add_gtest(matrix_cpu_gpu_test
               LIBS ${DCA_LIBS})
 
 dca_add_gtest(matrix_gpu_test
-              GTEST_MAIN
               CUDA
-              LIBS ${DCA_LIBS})
+              LIBS ${DCA_LIBS} ${DCA_KERNEL_LIBS})
 
 # dca_add_gtest(matrixop_complex_cpu_test
 #               GTEST_MAIN
@@ -50,11 +49,11 @@ dca_add_gtest(matrixop_complex_gpu_test
               CUDA
               LIBS ${DCA_LIBS} lapack lapack_kernels blas_kernels gpu_utils ${MAGMA_LIBARY} lapack_kernels)
 
-dca_add_gtest(reshapable_matrix_cpu_test GTEST_MAIN LIBS ${UTIL_LIBS})
+dca_add_gtest(reshapable_matrix_cpu_test GTEST_MAIN LIBS ${UTIL_LIBS} ${DCA_LIBS} magma::magma)
 
-dca_add_gtest(reshapable_matrix_gpu_test CUDA GTEST_MAIN LIBS gpu_utils)
+dca_add_gtest(reshapable_matrix_gpu_test CUDA GTEST_MAIN LIBS ${DCA_LIBS} gpu_utils magma::magma)
 
-dca_add_gtest(multi_vector_test GTEST_MAIN LIBS ${UTIL_LIBS})
+dca_add_gtest(multi_vector_test GTEST_MAIN LIBS ${DCA_LIBS} ${UTIL_LIBS})
 
 add_subdirectory(blas)
 add_subdirectory(util)
diff --git a/test/unit/linalg/gpu_test_util.hpp b/test/unit/linalg/gpu_test_util.hpp
index ce9f6a524..bcc3be0ae 100644
--- a/test/unit/linalg/gpu_test_util.hpp
+++ b/test/unit/linalg/gpu_test_util.hpp
@@ -27,7 +27,11 @@ cudaMemoryType PointerType(const ScalarType* ptr) {
 #if defined(DCA_HAVE_CUDA)
   return attributes.type;
 #elif defined(DCA_HAVE_HIP)
+  #if HIP_VERSION_MAJOR >= 6
+  return attributes.type;
+  #else
   return attributes.memoryType;
+  #endif
 #endif
 }
 
diff --git a/test/unit/linalg/matrix_gpu_test.cpp b/test/unit/linalg/matrix_gpu_test.cpp
index f9c66f9b7..126236404 100644
--- a/test/unit/linalg/matrix_gpu_test.cpp
+++ b/test/unit/linalg/matrix_gpu_test.cpp
@@ -15,6 +15,9 @@
 #include <utility>
 #include "gtest/gtest.h"
 #include "gpu_test_util.hpp"
+#include "dca/linalg/util/util_gpublas.hpp"
+#include "dca/linalg/util/info_gpu.hpp"
+
 
 TEST(MatrixGPUTest, Constructors) {
   int size = 3;
@@ -649,3 +652,17 @@ TEST(MatrixGPUTest, setTo) {
     for (int i = 0; i < mat_host.nrRows(); ++i)
       EXPECT_EQ(mat_host_ret(i,j), mat_host(i, j)); 
 }
+
+int main(int argc, char** argv) {
+
+      dca::linalg::util::printInfoDevices();
+  dca::linalg::util::initializeMagma();
+  ::testing::InitGoogleTest(&argc, argv);
+
+  // ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
+  // delete listeners.Release(listeners.default_result_printer());
+  // listeners.Append(new dca::testing::MinimalistPrinter);
+
+  int result = RUN_ALL_TESTS();
+  return result;
+}
diff --git a/test/unit/linalg/util/CMakeLists.txt b/test/unit/linalg/util/CMakeLists.txt
index 94c03d215..d5d54352e 100644
--- a/test/unit/linalg/util/CMakeLists.txt
+++ b/test/unit/linalg/util/CMakeLists.txt
@@ -8,7 +8,7 @@ dca_add_gtest(check_RC_test
 dca_add_gtest(magma_batched_gemm_test
 	      GTEST_MAIN
               CUDA
-              LIBS lapack gpu_utils magma::magma magma::sparse ${DCA_EXTERNAL_LIBS})
+              LIBS lapack gpu_utils magma::magma magma::sparse ${DCA_LIBS} ${DCA_EXTERNAL_LIBS})
 
 dca_add_gtest(stream_handle_test
               GTEST_MAIN
diff --git a/test/unit/math/function_transform/CMakeLists.txt b/test/unit/math/function_transform/CMakeLists.txt
index 0ac2b7041..f457065ec 100644
--- a/test/unit/math/function_transform/CMakeLists.txt
+++ b/test/unit/math/function_transform/CMakeLists.txt
@@ -1,12 +1,18 @@
 # unit/math/function_transform
+
+set(gtest_main_flag GTEST_MAIN)
+if(DCA_HAVE_ADIOS2)
+  set(gtest_main_flag GTEST_MPI_MAIN)
+endif(DCA_HAVE_ADIOS2)
+
 dca_add_gtest(function_transform_test
+        ${gtest_main_flag}
         INCLUDE_DIRS ${DCA_INCLUDES};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
-        LIBS FFTW::Double ${DCA_LIBS})
+        LIBS FFTW::Double ${DCA_LIBS} parallel_no_concurrency)
 
 dca_add_gtest(space_transform_2D_test
         INCLUDE_DIRS ${DCA_INCLUDES};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
-        LIBS FFTW::Double ${DCA_LIBS})
-
+        LIBS FFTW::Double ${DCA_LIBS} parallel_no_concurrency)
 
 dca_add_gtest(space_transform_2D_gpu_test
     CUDA
diff --git a/test/unit/math/function_transform/function_transform_test.cpp b/test/unit/math/function_transform/function_transform_test.cpp
index a742316ae..6a1f68dba 100644
--- a/test/unit/math/function_transform/function_transform_test.cpp
+++ b/test/unit/math/function_transform/function_transform_test.cpp
@@ -31,16 +31,6 @@ using McOptions = MockMcOptions<double>;
 #include "dca/parallel/no_concurrency/no_concurrency.hpp"
 #include "dca/parallel/no_threading/no_threading.hpp"
 #include "dca/profiling/null_profiler.hpp"
-#ifdef DCA_HAVE_ADIOS2
-adios2::ADIOS* adios_ptr;
-#endif
-#ifdef DCA_HAVE_MPI
-#include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
-dca::parallel::MPIConcurrency* concurrency_ptr;
-#else
-#include "dca/parallel/no_concurrency/no_concurrency.hpp"
-dca::parallel::NoConcurrency* concurrency_ptr;
-#endif
 
 using Model =
     dca::phys::models::TightBindingModel<dca::phys::models::RashbaHubbard<dca::phys::domains::D4>>;
@@ -152,35 +142,11 @@ void spTestImplementation(const bool direct) {
 }
 
 TEST(FunctionTransformTest, SpaceToMomentumCmplx) {
+  initialize();
   spTestImplementation<RDmn, KDmn>(true);
 }
 
 TEST(FunctionTransformTest, MomentumToSpaceCmplx) {
-  spTestImplementation<KDmn, RDmn>(false);
-}
-
-int main(int argc, char** argv) {
-#ifdef DCA_HAVE_MPI
-  dca::parallel::MPIConcurrency concurrency(argc, argv);
-  concurrency_ptr = &concurrency;
-#else
-  dca::parallel::NoConcurrency concurrency(argc, argv);
-  concurrency_ptr = &concurrency;
-#endif
-
-#ifdef DCA_HAVE_ADIOS2
-  //ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), "C++");
-  adios_ptr = &adios;
-#endif
-  ::testing::InitGoogleTest(&argc, argv);
-
-  // ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
-  // delete listeners.Release(listeners.default_result_printer());
-  // listeners.Append(new dca::testing::MinimalistPrinter);
-
   initialize();
-
-  int result = RUN_ALL_TESTS();
-  return result;
+  spTestImplementation<KDmn, RDmn>(false);
 }
diff --git a/test/unit/math/function_transform/space_transform_2D_gpu_test.cpp b/test/unit/math/function_transform/space_transform_2D_gpu_test.cpp
index dcec13579..319d0099c 100644
--- a/test/unit/math/function_transform/space_transform_2D_gpu_test.cpp
+++ b/test/unit/math/function_transform/space_transform_2D_gpu_test.cpp
@@ -142,8 +142,8 @@ TYPED_TEST(SpaceTransform2DGpuTest, Execute) {
 
 #ifdef DCA_HAVE_ADIOS2
   if (write_transforms) {
-    dca::io::Writer writer(*adios_ptr, *concurrency_ptr, "ADIOS2", true);
-    dca::io::Writer writer_h5(*adios_ptr, *concurrency_ptr, "HDF5", true);
+    dca::io::Writer writer(*concurrency_ptr, "ADIOS2", true);
+    dca::io::Writer writer_h5(*concurrency_ptr, "HDF5", true);
 
     //writer.open_file("tp_single_band_gpu_test_G4.bp");
     writer_h5.open_file("space_transform_2D_gpu_test.hdf5");
diff --git a/test/unit/math/function_transform/space_transform_2D_test.cpp b/test/unit/math/function_transform/space_transform_2D_test.cpp
index a43c99827..1dee64bb6 100644
--- a/test/unit/math/function_transform/space_transform_2D_test.cpp
+++ b/test/unit/math/function_transform/space_transform_2D_test.cpp
@@ -26,7 +26,9 @@
 #include "dca/parallel/no_concurrency/no_concurrency.hpp"
 #include "dca/parallel/no_threading/no_threading.hpp"
 #include "dca/profiling/null_profiler.hpp"
-
+#ifdef DCA_HAVE_ADIOS2
+#include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
+#endif
 dca::parallel::NoConcurrency* concurrency_ptr;
 
 using Model = dca::phys::models::TightBindingModel<
@@ -121,6 +123,11 @@ TYPED_TEST(SpaceTransform2DTest, Execute) {
 }
 
 int main(int argc, char** argv) {
+#ifdef DCA_HAVE_ADIOS2
+  // Until I figure out how to prevent globals inside the adios library from calling
+  // MPI_COMM_dup why're stuck with this if it is linked.
+  dca::parallel::MPIConcurrency mpi_concurrency(argc, argv);
+#endif
   dca::parallel::NoConcurrency concurrency(argc, argv);
   concurrency_ptr = &concurrency;
 
diff --git a/test/unit/math/inference/CMakeLists.txt b/test/unit/math/inference/CMakeLists.txt
index e888a009b..16383a26e 100644
--- a/test/unit/math/inference/CMakeLists.txt
+++ b/test/unit/math/inference/CMakeLists.txt
@@ -2,4 +2,4 @@
 
 dca_add_gtest(richardson_lucy_deconvolution_test
   GTEST_MAIN
-  LIBS function ${LAPACK_LIBRARIES} ${UTIL_LIBS})
+  LIBS function ${DCA_LIBS} ${LAPACK_LIBRARIES} ${UTIL_LIBS})
diff --git a/test/unit/math/random/std_random_wrapper_custom_seed_test.cpp b/test/unit/math/random/std_random_wrapper_custom_seed_test.cpp
index 81511d7da..781088f3d 100644
--- a/test/unit/math/random/std_random_wrapper_custom_seed_test.cpp
+++ b/test/unit/math/random/std_random_wrapper_custom_seed_test.cpp
@@ -34,6 +34,19 @@ inline std::vector<double>& getFirstRandomNumbersCustomSeed<std::mt19937_64>() {
                                  0.41722249930558986, 0.095007484445949017};
   return vec;
 }
+
+#if defined(__clang__)
+template <>
+inline std::vector<double>& getFirstRandomNumbersCustomSeed<std::ranlux48_base>() {
+  static std::vector<double> vec{0.59121028582286772,0.1686617322410229, 0.076703898677588533,0.73282393484576336,0.3330563551366435};
+  return vec;
+}
+template <>
+inline std::vector<double>& getFirstRandomNumbersCustomSeed<std::ranlux48>() {
+  static std::vector<double> vec{0.59121028582286772, 0.1686617322410229, 0.076703898677588533, 0.73282393484576336, 0.3330563551366435};
+  return vec;
+}
+#else
 template <>
 inline std::vector<double>& getFirstRandomNumbersCustomSeed<std::ranlux48_base>() {
   static std::vector<double> vec{0.0068332189419341447, 0.86802894579254541, 0.25453437079295799,
@@ -46,7 +59,7 @@ inline std::vector<double>& getFirstRandomNumbersCustomSeed<std::ranlux48>() {
                                  0.24682368736319102, 0.88335484565944544};
   return vec;
 }
-
+#endif  
 }  // testing
 }  // dca
 
diff --git a/test/unit/math/random/std_random_wrapper_default_seed_test.cpp b/test/unit/math/random/std_random_wrapper_default_seed_test.cpp
index 453397264..4ff99bb51 100644
--- a/test/unit/math/random/std_random_wrapper_default_seed_test.cpp
+++ b/test/unit/math/random/std_random_wrapper_default_seed_test.cpp
@@ -35,6 +35,18 @@ inline std::vector<double>& getFirstRandomNumbersDefaultSeed<std::mt19937_64>()
                                  0.79004886762619775, 0.73670732232143421};
   return vec;
 }
+#if defined (__clang__)
+template <>
+inline std::vector<double>& getFirstRandomNumbersDefaultSeed<std::ranlux48_base>() {
+  static std::vector<double> vec{0.67052872832964927, 0.89987293212401465, 0.48383422722711483,0.28699620000273923,0.37575706648963081};
+  return vec;
+}
+template <>
+inline std::vector<double>& getFirstRandomNumbersDefaultSeed<std::ranlux48>() {
+  static std::vector<double> vec{0.67052872832964927, 0.89987293212401465, 0.48383422722711483,0.28699620000273923,0.37575706648963081};
+  return vec;
+}
+#else
 template <>
 inline std::vector<double>& getFirstRandomNumbersDefaultSeed<std::ranlux48_base>() {
   static std::vector<double> vec{0.72777987290424495, 0.6717218115467265, 0.67916691101138849,
@@ -47,7 +59,7 @@ inline std::vector<double>& getFirstRandomNumbersDefaultSeed<std::ranlux48>() {
                                  0.040260746353792579, 0.94799739988482223};
   return vec;
 }
-
+#endif
 }  // testing
 }  // dca
 
diff --git a/test/unit/math/util/CMakeLists.txt b/test/unit/math/util/CMakeLists.txt
index b540f0f4f..82f43a7d6 100644
--- a/test/unit/math/util/CMakeLists.txt
+++ b/test/unit/math/util/CMakeLists.txt
@@ -4,7 +4,7 @@ dca_add_gtest(phase_test
 
 dca_add_gtest(vector_operations_test
   GTEST_MAIN
-  LIBS ${MAGMA_LIBRARY} ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS})
+  LIBS ${MAGMA_LIBRARY} ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS} ${DCA_LIBS})
 
 dca_add_gtest(comparison_methods_test
   GTEST_MAIN
diff --git a/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp b/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp
index 86b045d90..8f4e76978 100644
--- a/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp
+++ b/test/unit/parallel/mpi_concurrency/mpi_collective_sum_test.cpp
@@ -46,6 +46,8 @@ TEST_F(MPICollectiveSumTest, SumScalar) {
 
   scalar_expected = size_ * (size_ - 1) / 2;
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_EQ(scalar_expected, scalar_test);
 }
 
@@ -59,10 +61,13 @@ TEST_F(MPICollectiveSumTest, SumFunction) {
     function_test(i) = i * rank_;
 
   sum_interface_.sum(function_test);
-
+  
+  
   for (int i = 0; i < function_test.size(); i++)
     function_expected(i) = i * size_ * (size_ - 1) / 2;
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   for (int i = 0; i < function_test.size(); i++)
     EXPECT_EQ(function_expected(i), function_test(i));
 }
@@ -74,7 +79,7 @@ TEST_F(MPICollectiveSumTest, LeaveOneOutAvgAndSum) {
     values[i] = 3.14 + i;
     sum += values[i];
   }
-
+  std::cout << "looaas:" << sum << "\n";
   // Expected result
   const double expected = (sum - values[rank_]) / double(size_ - 1);
 
@@ -83,6 +88,9 @@ TEST_F(MPICollectiveSumTest, LeaveOneOutAvgAndSum) {
   double sum_one_out = scalar;
   sum_interface_.leaveOneOutAvg(scalar);
   sum_interface_.leaveOneOutSum(sum_one_out);
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_DOUBLE_EQ(expected, scalar);
   EXPECT_DOUBLE_EQ(expected * (size_ - 1), sum_one_out);
 
@@ -96,6 +104,8 @@ TEST_F(MPICollectiveSumTest, LeaveOneOutAvgAndSum) {
   sum_interface_.leaveOneOutAvg(f);
   sum_interface_.leaveOneOutSum(f_sum);
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_DOUBLE_EQ(expected, f(0));
   EXPECT_DOUBLE_EQ(expected * (size_ - 1), f_sum(0));
 
@@ -116,6 +126,8 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorReal) {
 
   auto err_trivial = sum_interface_.jackknifeError(f);
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_NEAR(0., err_trivial(0), epsilon_);
   EXPECT_NEAR(0., err_trivial(1), epsilon_);
 
@@ -160,6 +172,8 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorReal) {
   // Do not overwrite the jackknife estimates with their average.
   auto err_no_overwriting = sum_interface_.jackknifeError(f, false);
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_DOUBLE_EQ(f_copy(0), f(0));
   EXPECT_DOUBLE_EQ(f_copy(1), f(1));
 
@@ -190,6 +204,9 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorComplex) {
   f(1) = std::complex<double>(2.72, 3.4);
 
   auto err_trivial = sum_interface_.jackknifeError(f);
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_NEAR(0., err_trivial(0).real(), epsilon_);
   EXPECT_NEAR(0., err_trivial(1).real(), epsilon_);
   EXPECT_NEAR(0., err_trivial(0).imag(), epsilon_);
@@ -238,6 +255,8 @@ TEST_F(MPICollectiveSumTest, JackknifeErrorComplex) {
   // Do not overwrite the jackknife estimates with their average.
   auto err_no_overwriting = sum_interface_.jackknifeError(f, false);
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   EXPECT_DOUBLE_EQ(f_copy(0).real(), f(0).real());
   EXPECT_DOUBLE_EQ(f_copy(1).real(), f(1).real());
 
@@ -311,6 +330,8 @@ TEST_F(MPICollectiveSumTest, ComputeCovarianceScalar) {
   // Compute covariance matrix with respect to the mean of f.
   sum_interface_.computeCovarianceAndAvg(f, covariance_from_computeCovarianceAndAvg);
 
+  MPI_Barrier(MPI_COMM_WORLD);
+
   for (int i = 0; i < f_mean.size(); ++i)
     EXPECT_DOUBLE_EQ(f_mean(i), f(i));
 
@@ -335,6 +356,9 @@ TEST_F(MPICollectiveSumTest, ComputeCovarianceComplex) {
     f_mean(i) = f(i);
   }
   sum_interface_.sum(f_mean);
+
+  MPI_Barrier(MPI_COMM_WORLD);
+
   f_mean /= size_;
 
   /// Calculate the equivalent covariance matrix without using function's etc
@@ -390,8 +414,11 @@ TEST_F(MPICollectiveSumTest, ComputeCovarianceComplex) {
 
   // Compute covariance matrix with respect to the *precomputed* mean of f.
   sum_interface_.computeCovariance(f, f_mean, covariance_from_computeCovariance);
+  MPI_Barrier(MPI_COMM_WORLD);
+
   // Compute covariance matrix with respect to the mean of f.
   sum_interface_.computeCovarianceAndAvg(f, covariance_from_computeCovarianceAndAvg);
+  MPI_Barrier(MPI_COMM_WORLD);
 
   for (int i = 0; i < f_mean.size(); ++i) {
     EXPECT_DOUBLE_EQ(f_mean(i).real(), f(i).real());
@@ -420,11 +447,12 @@ TEST_F(MPICollectiveSumTest, AvgNormalizedMomenta) {
 
   std::vector<double> momenta(orders.size());
 
+  int rank = rank_ + 1;
+  
   for (int i = 0; i < f.size(); i++) {
     f(i) = rank_ * (i + 1);
   }
 
-
   std::vector<std::vector<double>> m_momenta;
   for (int i = 0; i < orders.size(); ++i)
     m_momenta.emplace_back(f.size(), 0.0);
@@ -445,6 +473,7 @@ TEST_F(MPICollectiveSumTest, AvgNormalizedMomenta) {
       f_mean[i] += ir * (i + 1);
     }
   }
+
   for (int i = 0; i < f.size(); ++i)
     f_mean[i] /= size_;
 
@@ -455,6 +484,7 @@ TEST_F(MPICollectiveSumTest, AvgNormalizedMomenta) {
 
     unormalizedMomenta(samples, f_mean);
   }
+
   std::vector<double> expected_momenta_avg(orders.size(), 0.0);
   for (int i = 0; i < f.size(); ++i) {
     const double var = std::sqrt(var2[i] / size_);
@@ -466,6 +496,8 @@ TEST_F(MPICollectiveSumTest, AvgNormalizedMomenta) {
     expected_momenta_avg[ir] /= f.size();
 
   momenta = sum_interface_.avgNormalizedMomenta(f, orders);
+  MPI_Barrier(MPI_COMM_WORLD);
+
   // Expected values obtained with python
   EXPECT_NEAR(0., momenta[0], 1e-8);
   EXPECT_NEAR(expected_momenta_avg[1], momenta[1], 1e-8);
@@ -496,6 +528,7 @@ TEST_F(MPICollectiveSumTest, DelayedSum) {
   sum_interface_.delayedSum(function_cmplx_test);
 
   sum_interface_.resolveSums();
+  MPI_Barrier(MPI_COMM_WORLD);
 
   EXPECT_EQ(scalar_expected, scalar_test);
   EXPECT_EQ(function_expected, function_test);
diff --git a/test/unit/parallel/no_concurrency/CMakeLists.txt b/test/unit/parallel/no_concurrency/CMakeLists.txt
index 2ffeb6638..f75738b3c 100644
--- a/test/unit/parallel/no_concurrency/CMakeLists.txt
+++ b/test/unit/parallel/no_concurrency/CMakeLists.txt
@@ -4,14 +4,26 @@ set(optional_io_lib "")
 
 if (DCA_HAVE_ADIOS2)
   set(optional_io_lib dca_adios2)
-endif()
 
-dca_add_gtest(no_concurrency_test
-  GTEST_MAIN
-  LIBS function
-  ${optional_io_lib}
+  dca_add_gtest(no_concurrency_test
+    GTEST_MPI_MAIN
+    LIBS function
+    ${optional_io_lib}
+    ${DCA_LIBS}
+    parallel_no_concurrency
+  )
+
+else()
+  dca_add_gtest(no_concurrency_test
+    GTEST_MAIN
+    LIBS function
+    ${optional_io_lib}
+    ${DCA_LIBS}
+    parallel_no_concurrency
   )
 
+endif()
+
 dca_add_gtest(serial_collective_sum_test
   GTEST_MAIN
   LIBS function)
diff --git a/test/unit/phys/dca_algorithms/CMakeLists.txt b/test/unit/phys/dca_algorithms/CMakeLists.txt
index 5ad1ae048..5531b7c4f 100644
--- a/test/unit/phys/dca_algorithms/CMakeLists.txt
+++ b/test/unit/phys/dca_algorithms/CMakeLists.txt
@@ -1,10 +1,15 @@
 # DCA algorithms unit tests
 
+set(gtest_main_flag GTEST_MAIN)
+if(DCA_HAVE_ADIOS2)
+  set(gtest_main_flag GTEST_MPI_MAIN)
+endif(DCA_HAVE_ADIOS2)
+
 dca_add_gtest(compute_band_structure_test
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${PROJECT_SOURCE_DIR}
   LIBS FFTW::Double dca_io json function cluster_domains enumerations time_and_frequency_domains quantum_domains ${UTIL_LIBS}
-  ${LAPACK_LIBRARIES})
+  ${LAPACK_LIBRARIES} parallel_no_concurrency)
 
 if (DCA_HAVE_ADIOS2)
   set(optional_io_lib dca_io)
@@ -12,10 +17,10 @@ endif()
 
 # \todo I think this should be able to run without GPU
 dca_add_gtest(compute_free_greens_function_test
-  GTEST_MAIN
+  ${gtest_main_flag}
   CUDA
-  LIBS function FFTW::Double ${DCA_THREADING_LIBS} ${LAPACK_LIBRARIES} ${UTIL_LIBS} ${DCA_GPU_LIBS} ${optional_io_lib} parallel_util)
+  LIBS function FFTW::Double ${DCA_THREADING_LIBS} ${LAPACK_LIBRARIES} ${UTIL_LIBS} ${DCA_GPU_LIBS} ${optional_io_lib} parallel_util parallel_no_concurrency)
 
 dca_add_gtest(compute_greens_function_test
-  GTEST_MAIN
-  LIBS function FFTW::Double ${DCA_THREADING_LIBS} ${LAPACK_LIBRARIES} ${UTIL_LIBS} ${DCA_GPU_LIBS} ${optional_io_lib} parallel_util)
+  ${gtest_main_flag}
+  LIBS function FFTW::Double ${DCA_THREADING_LIBS} ${LAPACK_LIBRARIES} ${UTIL_LIBS} ${DCA_GPU_LIBS} ${optional_io_lib} parallel_util parallel_no_concurrency)
diff --git a/test/unit/phys/dca_data/CMakeLists.txt b/test/unit/phys/dca_data/CMakeLists.txt
index 10ca03917..e27e5959c 100644
--- a/test/unit/phys/dca_data/CMakeLists.txt
+++ b/test/unit/phys/dca_data/CMakeLists.txt
@@ -3,4 +3,4 @@
 dca_add_gtest(dca_data_test
   FAST
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS function ${DCA_LIBS})
+  LIBS function ${DCA_LIBS} parallel_no_concurrency)
diff --git a/test/unit/phys/dca_data/dca_data_test.cpp b/test/unit/phys/dca_data/dca_data_test.cpp
index 69bd85bb5..9323c5391 100644
--- a/test/unit/phys/dca_data/dca_data_test.cpp
+++ b/test/unit/phys/dca_data/dca_data_test.cpp
@@ -44,13 +44,15 @@ constexpr char dca_data_test_input[] = DCA_SOURCE_DIR "/test/unit/phys/dca_data/
 
 template <typename Scalar, class Concurrency, class Lattice = LatticeSquare,
           ClusterSolverId solver_name = ClusterSolverId::CT_AUX,
-          const char* input_name = dca_data_test_input, DistType DT = DistType::NONE, class NUMTRAITS = dca::NumericalTraits<dca::util::RealAlias<Scalar>, Scalar>>
+          const char* input_name = dca_data_test_input, DistType DT = DistType::NONE,
+          class NUMTRAITS = dca::NumericalTraits<dca::util::RealAlias<Scalar>, Scalar>>
 struct DCADataSetup {
   using LatticeType = Lattice;
   using Model = phys::models::TightBindingModel<Lattice>;
   using RngType = testing::StubRng;
-  using Parameters = phys::params::Parameters<Concurrency, parallel::NoThreading,
-                                              profiling::NullProfiler, Model, RngType, solver_name, NUMTRAITS>;
+  using Parameters =
+      phys::params::Parameters<Concurrency, parallel::NoThreading, profiling::NullProfiler, Model,
+                               RngType, solver_name, NUMTRAITS>;
   using Data = phys::DcaData<Parameters, DT>;
   using LoopData = phys::DcaLoopData<Parameters>;
 
@@ -159,8 +161,6 @@ TEST_F(DCADataTest, ReadSigmaAdios2AnotherPot) {
 
 #endif
 
-
-
 int main(int argc, char** argv) {
 #ifdef DCA_HAVE_MPI
   dca::parallel::MPIConcurrency concurrency(argc, argv);
@@ -172,7 +172,7 @@ int main(int argc, char** argv) {
 
 #ifdef DCA_HAVE_ADIOS2
   // ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), false);
+  adios2::ADIOS adios("", concurrency_ptr->get());
   adios_ptr = &adios;
 #endif
   ::testing::InitGoogleTest(&argc, argv);
@@ -186,11 +186,7 @@ int main(int argc, char** argv) {
     dca_data.loop_data_->last_completed_iteration = 0;
     dca_data.data_->Sigma = 1.0;
     std::string format{"HDF5"};
-    dca::io::Writer<Concurrency> writer(
-#ifdef DCA_HAVE_ADIOS2
-        *adios_ptr,
-#endif
-        *concurrency_ptr, format);
+    dca::io::Writer<Concurrency> writer(*concurrency_ptr, format);
     writer.open_file("dca_data_test.hdf5", true);
     writer.begin_step();
     dca_data.data_->write(writer);
@@ -213,11 +209,7 @@ int main(int argc, char** argv) {
     writer.close_file();
 
     // Fake legacy hdf5 file by omitting steps
-    dca::io::Writer<Concurrency> writer_legacy(
-#ifdef DCA_HAVE_ADIOS2
-        *adios_ptr,
-#endif
-        *concurrency_ptr, format);
+    dca::io::Writer<Concurrency> writer_legacy(*concurrency_ptr, format);
     writer_legacy.open_file("dca_data_test_legacy.hdf5", true);
     dca_data.loop_data_->chemical_potential(3) = 4.0;
     dca_data.loop_data_->last_completed_iteration = 3;
@@ -239,45 +231,43 @@ int main(int argc, char** argv) {
     dca_data.loop_data_->last_completed_iteration = 0;
     dca_data.data_->Sigma = 1.0;
     std::string format{"ADIOS2"};
-    dca::io::Writer<Concurrency> writer(*adios_ptr, *concurrency_ptr, format);
+    dca::io::Writer<Concurrency> writer(*concurrency_ptr, format);
     writer.open_file("dca_data_test.bp", true);
-    //insure chemical potentials past the last complete iteration are ignored.
-    dca::io::Writer<Concurrency> writer_extra_cpot(
-        *adios_ptr,
-        *concurrency_ptr, format);
+    // insure chemical potentials past the last complete iteration are ignored.
+    dca::io::Writer<Concurrency> writer_extra_cpot(*concurrency_ptr, format);
     writer_extra_cpot.open_file("dca_data_test_extra_cpot.bp", true);
-    //insure chemical potentials past the last complete iteration are ignored.
-    dca::io::Writer<Concurrency> writer_another_cpot(
-        *adios_ptr,
-        *concurrency_ptr, format);
+    // insure chemical potentials past the last complete iteration are ignored.
+    dca::io::Writer<Concurrency> writer_another_cpot(*concurrency_ptr, format);
     writer_another_cpot.open_file("dca_data_test_another_cpot.bp", true);
 
-    std::vector<std::reference_wrapper<dca::io::Writer<Concurrency>>> writers{writer, writer_extra_cpot, writer_another_cpot};
+    std::vector<std::reference_wrapper<dca::io::Writer<Concurrency>>> writers{
+        writer, writer_extra_cpot, writer_another_cpot};
 
     auto begin_steps = [](dca::io::Writer<Concurrency>& writer) { writer.begin_step(); };
-    auto writes = [&dca_data](dca::io::Writer<Concurrency>& writer) { dca_data.data_->write(writer);
-      dca_data.loop_data_->write(writer, *concurrency_ptr);};
+    auto writes = [&dca_data](dca::io::Writer<Concurrency>& writer) {
+      dca_data.data_->write(writer);
+      dca_data.loop_data_->write(writer, *concurrency_ptr);
+    };
     auto end_steps = [](dca::io::Writer<Concurrency>& writer) { writer.end_step(); };
 
     std::for_each(writers.begin(), writers.end(), begin_steps);
     std::for_each(writers.begin(), writers.end(), writes);
-    std::for_each(writers.begin(), writers.end(), end_steps);    
+    std::for_each(writers.begin(), writers.end(), end_steps);
 
     dca_data.loop_data_->chemical_potential(1) = 2.0;
     dca_data.loop_data_->last_completed_iteration = 1;
     dca_data.data_->Sigma = 2.0;
     std::for_each(writers.begin(), writers.end(), begin_steps);
     std::for_each(writers.begin(), writers.end(), writes);
-    std::for_each(writers.begin(), writers.end(), end_steps);    
+    std::for_each(writers.begin(), writers.end(), end_steps);
 
     dca_data.loop_data_->chemical_potential(2) = 3.0;
     dca_data.loop_data_->last_completed_iteration = 2;
     dca_data.data_->Sigma = 3.0;
     std::for_each(writers.begin(), writers.end(), begin_steps);
     std::for_each(writers.begin(), writers.end(), writes);
-    std::for_each(writers.begin(), writers.end(), end_steps);    
-    
-    
+    std::for_each(writers.begin(), writers.end(), end_steps);
+
     dca_data.loop_data_->chemical_potential(3) = 4.0;
     dca_data.data_->Sigma = 4.0;
     writer_extra_cpot.begin_step();
@@ -286,7 +276,6 @@ int main(int argc, char** argv) {
     writer_extra_cpot.end_step();
     writer_extra_cpot.close_file();
 
-    
     dca_data.loop_data_->last_completed_iteration = 3;
     writer_another_cpot.begin_step();
     dca_data.data_->write(writer_another_cpot);
diff --git a/test/unit/phys/dca_step/cluster_solver/ctaux/accumulator/tp/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctaux/accumulator/tp/CMakeLists.txt
index 9cbce1d09..ea801d402 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctaux/accumulator/tp/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctaux/accumulator/tp/CMakeLists.txt
@@ -1,6 +1,6 @@
 dca_add_gtest(tp_equal_time_accumulator_test
   FAST
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS FFTW::Double ${DCA_LIBS} ${KERNELS_LIB}
+  LIBS FFTW::Double ${DCA_LIBS} ${KERNELS_LIB} parallel_no_concurrency
   )
diff --git a/test/unit/phys/dca_step/cluster_solver/ctaux/structs/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctaux/structs/CMakeLists.txt
index 8d2d86bc6..33414c352 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctaux/structs/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctaux/structs/CMakeLists.txt
@@ -2,7 +2,7 @@
 
 dca_add_gtest(read_write_config_test
         FAST
-        GTEST_MAIN
+        GTEST_MPI_MAIN
         INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
         LIBS     FFTW::Double ${DCA_LIBS}
         )
diff --git a/test/unit/phys/dca_step/cluster_solver/ctaux/walker/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctaux/walker/CMakeLists.txt
index 5409d9e93..3c07dfcb8 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctaux/walker/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctaux/walker/CMakeLists.txt
@@ -7,7 +7,7 @@ message("DCA_KERNEL_LIBS: ${DCA_KERNEL_LIBS}")
   
 dca_add_gtest(ct_aux_walker_test
     FAST
-    GTEST_MAIN
+    GTEST_MPI_MAIN
     CUDA
     INCLUDE_DIRS ${TEST_INCLUDES}
     LIBS     FFTW::Double ${TEST_LIBS} g0_interpolation ${DCA_KERNEL_LIBS}
diff --git a/test/unit/phys/dca_step/cluster_solver/ctint/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctint/CMakeLists.txt
index 033d66b84..70d824688 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctint/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctint/CMakeLists.txt
@@ -4,7 +4,7 @@ set(TEST_LIBS     ${DCA_LIBS})
 dca_add_gtest(ct_int_solver_methods_test
   FAST
   CUDA
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${TEST_INCLUDES}
   LIBS     FFTW::Double ${TEST_LIBS} ctint g0_interpolation ${KERNELS_LIB}
 )
diff --git a/test/unit/phys/dca_step/cluster_solver/ctint/structs/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctint/structs/CMakeLists.txt
index ba453b3ca..d2933ef52 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctint/structs/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctint/structs/CMakeLists.txt
@@ -3,7 +3,7 @@ set(TEST_INCLUDES ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR})
 
 dca_add_gtest(ct_int_interaction_vertices_test
     FAST
-    GTEST_MAIN
+    GTEST_MPI_MAIN
     INCLUDE_DIRS ${TEST_INCLUDES}
     CUDA
     LIBS FFTW::Double ${DCA_LIBS}
diff --git a/test/unit/phys/dca_step/cluster_solver/ctint/walker/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctint/walker/CMakeLists.txt
index 2ccf46731..d98942217 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctint/walker/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctint/walker/CMakeLists.txt
@@ -5,21 +5,21 @@ set(TEST_LIBS     ${DCA_LIBS})
 
 dca_add_gtest(ct_int_walker_test
     FAST
-    GTEST_MAIN
+    GTEST_MPI_MAIN
     INCLUDE_DIRS ${TEST_INCLUDES}
     LIBS     FFTW::Double ${TEST_LIBS} g0_interpolation ${DCA_KERNEL_LIBS}
     )
 
 dca_add_gtest(ct_int_walker_submatrix_test
         FAST
-        GTEST_MAIN
+        GTEST_MPI_MAIN
         INCLUDE_DIRS ${TEST_INCLUDES}
         LIBS     ${TEST_LIBS} g0_interpolation ${DCA_GPU_LIBS}
         )
 
 dca_add_gtest(ct_int_walker_multinsertion_test
         FAST
-        GTEST_MAIN
+        GTEST_MPI_MAIN
         INCLUDE_DIRS ${TEST_INCLUDES}
         LIBS     ${TEST_LIBS} g0_interpolation
         )
@@ -27,7 +27,7 @@ dca_add_gtest(ct_int_walker_multinsertion_test
 dca_add_gtest(ct_int_walker_submatrix_gpu_test
         CUDA
         FAST
-        GTEST_MAIN
+        GTEST_MPI_MAIN
         INCLUDE_DIRS ${TEST_INCLUDES}
         LIBS     ${TEST_LIBS} g0_interpolation ${DCA_GPU_LIBS} ctint ${DCA_KERNEL_LIBS}
       )
@@ -35,7 +35,7 @@ dca_add_gtest(ct_int_walker_submatrix_gpu_test
 dca_add_gtest(ct_int_walker_submatrix_gpu_complex_test
         CUDA
         FAST
-        GTEST_MAIN
+        GTEST_MPI_MAIN
         INCLUDE_DIRS ${TEST_INCLUDES}
         LIBS     ${TEST_LIBS} g0_interpolation ${DCA_GPU_LIBS} ctint ${DCA_KERNEL_LIBS}
         )
diff --git a/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/CMakeLists.txt
index db10ac4ec..934ca8827 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/CMakeLists.txt
@@ -2,14 +2,10 @@ set(TEST_INCLUDES ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR})
 set(TEST_LIBS     ${DCA_LIBS})
 
 set(THIS_TEST_LIBS ${DCA_LIBS} ctint)
-if (DCA_HAVE_ADIOS2)
-  set(THIS_TEST_LIBS ${THIS_TEST_LIBS} dca_adios2)
-endif()
 
 dca_add_gtest(d_matrix_builder_gpu_test
     FAST
     CUDA
-    GTEST_MAIN
     INCLUDE_DIRS ${TEST_INCLUDES}
     LIBS     ${THIS_TEST_LIBS} g0_interpolation ctint ${DCA_KERNEL_LIBS}
     )
diff --git a/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu_test.cpp b/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu_test.cpp
index bd2e1339d..220e6e00d 100644
--- a/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/ctint/walker/tools/d_matrix_builder_gpu_test.cpp
@@ -117,3 +117,14 @@ TYPED_TEST(DMatrixBuilderGpuTest, RemoveAndInstertVertex) {
     }
   }
 }
+
+int main(int argc, char** argv) {
+  int result = 0;
+
+  dca::parallel::MPIConcurrency concurrency(argc, argv);
+  ::testing::InitGoogleTest(&argc, argv);
+
+  result = RUN_ALL_TESTS();
+
+  return result;
+}
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/CMakeLists.txt b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/CMakeLists.txt
index b40cb0f9a..7f99ed90f 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/CMakeLists.txt
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/CMakeLists.txt
@@ -2,7 +2,7 @@
 
 dca_add_gtest(hdf5_g4_file_io_test
   FAST
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
   LIBS     ${DCA_LIBS}
   )
@@ -12,26 +12,28 @@ if (DCA_HAVE_ADIOS2)
     FAST
     INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
     LIBS     ${DCA_LIBS}
-    )
-dca_add_gtest(tp_accumulator_test
-  FAST
-  INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
-  LIBS     ${DCA_LIBS} ${DCA_GPU_LIBS}
   )
 
-  endif()
+endif()
 
+# Baseline write is broken
+# dca_add_gtest(tp_accumulator_test
+#   FAST
+#   GTEST_MPI_MAIN
+#   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
+#   LIBS enumerations function timer parallel_stdthread parallel_util dca_domains dca_math modern_string_utils dca_io parallel_no_concurrency
+#   )
 
 dca_add_gtest(tp_accumulator_test_complex
   FAST
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
   LIBS     ${DCA_LIBS}
   )
 
 dca_add_gtest(tp_accumulator_singleband_test
   FAST
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR}
   LIBS     ${DCA_LIBS}
   )
@@ -46,10 +48,10 @@ endif()
 dca_add_gtest(tp_accumulator_particle_hole_test
   FAST
   CUDA
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   INCLUDE_DIRS ${PROJECT_SOURCE_DIR} ${FFTW_INCLUDE_DIR} ${HDF5_INCLUDE_DIR}
   LIBS function FFTW::Double cluster_domains time_and_frequency_domains quantum_domains timer json random
-       enumerations dca_algorithms ${DCA_THREADING_LIBS} dca_io parallel_util ${LAPACK_LIBRARIES} gpu_utils)
+       enumerations dca_algorithms ${DCA_THREADING_LIBS} dca_io parallel_util ${LAPACK_LIBRARIES} gpu_utils parallel_no_concurrency)
 
 dca_add_gtest(tp_accumulator_gpu_test
   FAST
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/adios2_g4_file_io_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/adios2_g4_file_io_test.cpp
index f0fc3e133..c734ab7e7 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/adios2_g4_file_io_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/adios2_g4_file_io_test.cpp
@@ -67,8 +67,8 @@ TEST_F(G4FileIoTest, ReadWrite) {
       G4(i) = rng();
   };
 
-  dca::io::ADIOS2Writer writer(*adios_ptr, concurrency_ptr);
-  dca::io::ADIOS2Reader reader(*adios_ptr, concurrency_ptr);;
+  dca::io::ADIOS2Writer writer(concurrency_ptr);
+  dca::io::ADIOS2Reader reader(*concurrency_ptr);;
 
   std::map<dca::phys::FourPointType, std::string> func_names;
   func_names[dca::phys::FourPointType::PARTICLE_HOLE_CHARGE] = "G4_ph_charge";
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_complex_g0_gpu_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_complex_g0_gpu_test.cpp
index 504170f66..909eee6a7 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_complex_g0_gpu_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_complex_g0_gpu_test.cpp
@@ -141,12 +141,12 @@ TYPED_TEST(TpAccumulatorComplexG0GpuTest, Accumulate) {
 
 #ifdef DCA_HAVE_ADIOS2
   if (write_G4s) {
-    dca::io::Writer writer(*adios_ptr, *concurrency_ptr, "ADIOS2", true);
-    dca::io::Writer writer_h5(*adios_ptr, *concurrency_ptr, "HDF5", true);
+    dca::io::Writer writer(*concurrency_ptr, "ADIOS2", true);
+    dca::io::Writer writer_h5(*concurrency_ptr, "HDF5", true);
 
     writer.open_file("tp_gpu_test_complex_G0_G4.bp");
     writer_h5.open_file("tp_gpu_test_complex_G0_G4.hdf5");
-
+    writer.begin_step();
     this->host_setup.parameters_.write(writer);
     this->host_setup.parameters_.write(writer_h5);
     this->host_setup.data_->write(writer);
@@ -160,7 +160,7 @@ TYPED_TEST(TpAccumulatorComplexG0GpuTest, Accumulate) {
       writer_h5.execute("accumulatorHOST_" + channel_str, accumulatorHost.get_G4()[channel]);
       writer_h5.execute("accumulatorDevice_" + channel_str, accumulatorDevice.get_G4()[channel]);
     }
-
+    
 #ifndef NDEBUG
     writer_h5.execute("accumulatorDevice_G_0", accumulatorDevice.get_G_Debug()[0]);
     writer_h5.execute("accumulatorDevice_G_1", accumulatorDevice.get_G_Debug()[1]);
@@ -206,7 +206,7 @@ TYPED_TEST(TpAccumulatorComplexG0GpuTest, Accumulate) {
             << "( " << i << ", " << j << " )";
       }
 #endif
-
+    writer.end_step();
     writer.close_file();
     writer_h5.close_file();
   }
@@ -262,7 +262,7 @@ int main(int argc, char** argv) {
 
 #ifdef DCA_HAVE_ADIOS2
   // ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), false);
+  adios2::ADIOS adios("", concurrency_ptr->get());
   adios_ptr = &adios;
 #endif
   ::testing::InitGoogleTest(&argc, argv);
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_test.cpp
index a0c1cf5f0..1eb09a529 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_gpu_test.cpp
@@ -113,12 +113,13 @@ TYPED_TEST(TpAccumulatorGpuTest, Accumulate) {
 
 #ifdef DCA_HAVE_ADIOS2
   if (write_G4s) {
-    dca::io::Writer writer(*adios_ptr, *concurrency_ptr, "ADIOS2", true);
-    dca::io::Writer writer_h5(*adios_ptr, *concurrency_ptr, "HDF5", true);
+    dca::io::Writer writer(*concurrency_ptr, "ADIOS2", true);
+    dca::io::Writer writer_h5(*concurrency_ptr, "HDF5", true);
 
     writer.open_file("tp_gpu_test_G4.bp");
     writer_h5.open_file("tp_gpu_test_G4.hdf5");
 
+    writer.begin_step();
     this->host_setup.parameters_.write(writer);
     this->host_setup.parameters_.write(writer_h5);
     this->host_setup.data_->write(writer);
@@ -131,6 +132,7 @@ TYPED_TEST(TpAccumulatorGpuTest, Accumulate) {
       writer_h5.execute("accumulatorHOST_" + channel_str, accumulatorHost.get_G4()[channel]);
       writer_h5.execute("accumulatorDevice_" + channel_str, accumulatorDevice.get_G4()[channel]);
     }
+    writer.end_step();
     writer.close_file();
     writer_h5.close_file();
   }
@@ -254,7 +256,7 @@ int main(int argc, char** argv) {
 
 #ifdef DCA_HAVE_ADIOS2
   // ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), false);
+  adios2::ADIOS adios("", concurrency_ptr->get());
   adios_ptr = &adios;
 #endif
   ::testing::InitGoogleTest(&argc, argv);
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_multiband_gpu_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_multiband_gpu_test.cpp
index 30bbe45b6..0ca1ec6ab 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_multiband_gpu_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_multiband_gpu_test.cpp
@@ -123,12 +123,12 @@ TYPED_TEST(TpAccumulatorGpuTest, Accumulate) {
 
 #ifdef DCA_HAVE_ADIOS2
   if (write_G4s) {
-    dca::io::Writer writer(*adios_ptr, *concurrency_ptr, "ADIOS2", true);
-    dca::io::Writer writer_h5(*adios_ptr, *concurrency_ptr, "HDF5", true);
+    dca::io::Writer writer(*concurrency_ptr, "ADIOS2", true);
+    dca::io::Writer writer_h5(*concurrency_ptr, "HDF5", true);
 
     writer.open_file("tp_gpu_test_G4.bp");
     writer_h5.open_file("tp_gpu_test_G4.hdf5");
-
+    writer.begin_step();
     this->host_setup.parameters_.write(writer);
     this->host_setup.parameters_.write(writer_h5);
     this->host_setup.data_->write(writer);
@@ -184,7 +184,7 @@ TYPED_TEST(TpAccumulatorGpuTest, Accumulate) {
 	EXPECT_NEAR(G_down(i, j).imag(), G_spin_separated[1](i, j).imag(), 1E-12) << "( " << i << ", " << j << " )";
       }
 #endif
-
+    writer.end_step();
     writer.close_file();
     writer_h5.close_file();
   }
@@ -407,7 +407,7 @@ int main(int argc, char** argv) {
 
 #ifdef DCA_HAVE_ADIOS2
   // ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), false);
+  adios2::ADIOS adios("", concurrency_ptr->get());
   adios_ptr = &adios;
 #endif
   ::testing::InitGoogleTest(&argc, argv);
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_singleband_gpu_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_singleband_gpu_test.cpp
index 410be1a5d..3123f8404 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_singleband_gpu_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_singleband_gpu_test.cpp
@@ -100,12 +100,12 @@ TEST_F(TpAccumulatorGpuSinglebandTest, Accumulate) {
 
 #ifdef DCA_HAVE_ADIOS2
   if (write_G4s) {
-    dca::io::Writer writer(*adios_ptr, *concurrency_ptr, "ADIOS2", true);
-    dca::io::Writer writer_h5(*adios_ptr, *concurrency_ptr, "HDF5", true);
+    dca::io::Writer writer(*concurrency_ptr, "ADIOS2", true);
+    dca::io::Writer writer_h5(*concurrency_ptr, "HDF5", true);
 
     writer.open_file("tp_single_band_gpu_test_G4.bp");
     writer_h5.open_file("tp_single_band_gpu_test_G4.hdf5");
-
+    writer.begin_step();
     parameters_.write(writer);
     parameters_.write(writer_h5);
     data_->write(writer);
@@ -121,6 +121,7 @@ TEST_F(TpAccumulatorGpuSinglebandTest, Accumulate) {
       writer_h5.execute("accumulatorHOST_" + channel_str, accumulatorHost.get_G4()[channel]);
       writer_h5.execute("accumulatorDevice_" + channel_str, accumulatorDevice.get_G4()[channel]);
     }
+    writer.end_step();
     writer.close_file();
     writer_h5.close_file();
   }
@@ -143,7 +144,7 @@ int main(int argc, char** argv) {
 
 #ifdef DCA_HAVE_ADIOS2
   // ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), false);
+  adios2::ADIOS adios("", concurrency_ptr->get());
   adios_ptr = &adios;
 #endif
 
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test.cpp b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test.cpp
index 4f0fb5a47..77747acc2 100644
--- a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test.cpp
+++ b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test.cpp
@@ -35,21 +35,9 @@ using McOptions = MockMcOptions<Scalar>;
 #include "dca/phys/four_point_type.hpp"
 #include "test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/accumulation_test.hpp"
 
-#ifdef DCA_HAVE_ADIOS2
-adios2::ADIOS* adios_ptr;
-#endif
-
-#ifdef DCA_HAVE_MPI
-#include "dca/parallel/mpi_concurrency/mpi_concurrency.hpp"
-dca::parallel::MPIConcurrency* concurrency_ptr;
-#else
-#include "dca/parallel/no_concurrency/no_concurrency.hpp"
-dca::parallel::NoConcurrency* concurrency_ptr;
-#endif
-
 constexpr bool update_baseline = true;
 
-constexpr bool write_G4s = true;
+constexpr bool write_G4s = false;
 
 #define INPUT_DIR \
   DCA_SOURCE_DIR "/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/"
@@ -109,9 +97,11 @@ TEST_F(TpAccumulatorTest, Accumulate) {
 
   const auto& G4 = accumulator.get_G4();
 
+  auto& concurrency = parameters_.get_concurrency();
+
   if (write_G4s) {
-    dca::io::Writer writer(*adios_ptr, *concurrency_ptr, "ADIOS2", true);
-    dca::io::Writer writer_h5(*adios_ptr, *concurrency_ptr, "HDF5", true);
+    dca::io::Writer writer(concurrency, "ADIOS2", true);
+    dca::io::Writer writer_h5(concurrency, "HDF5", true);
     writer.open_file("tp_accumulator_test_G4.bp");
     writer_h5.open_file("tp_accumulator_test_G4.hdf5");
     parameters_.write(writer);
@@ -128,8 +118,8 @@ TEST_F(TpAccumulatorTest, Accumulate) {
   }
 
   if (update_baseline) {
-    baseline_writer.end_step();
     data_->write(baseline_writer);
+    baseline_writer.end_step();
     baseline_writer.close_file();
   }
   else {
@@ -138,33 +128,7 @@ TEST_F(TpAccumulatorTest, Accumulate) {
       reader.execute(G4_check);
       const auto diff = dca::func::util::difference(G4[0], G4_check);
       EXPECT_GT(1e-8, diff.l_inf);
-
       reader.close_file();
     }
   }
 }
-
-int main(int argc, char** argv) {
-#ifdef DCA_HAVE_MPI
-  dca::parallel::MPIConcurrency concurrency(argc, argv);
-  concurrency_ptr = &concurrency;
-#else
-  dca::parallel::NoConcurrency concurrency(argc, argv);
-  concurrency_ptr = &concurrency;
-#endif
-
-#ifdef DCA_HAVE_ADIOS2
-  // ADIOS expects MPI_COMM pointer or nullptr
-  adios2::ADIOS adios("", concurrency_ptr->get(), false);
-  adios_ptr = &adios;
-#endif
-
-  ::testing::InitGoogleTest(&argc, argv);
-
-  // ::testing::TestEventListeners& listeners = ::testing::UnitTest::GetInstance()->listeners();
-  // delete listeners.Release(listeners.default_result_printer());
-  // listeners.Append(new dca::testing::MinimalistPrinter);
-
-  int result = RUN_ALL_TESTS();
-  return result;
-}
diff --git a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_baseline.hdf5 b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_baseline.hdf5
index ea94fcc16..8f96dea2f 100644
Binary files a/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_baseline.hdf5 and b/test/unit/phys/dca_step/cluster_solver/shared_tools/accumulation/tp/tp_accumulator_test_baseline.hdf5 differ
diff --git a/test/unit/phys/dca_step/lattice_mapping/deconvolution/CMakeLists.txt b/test/unit/phys/dca_step/lattice_mapping/deconvolution/CMakeLists.txt
index 60ac2491f..ebb551b98 100644
--- a/test/unit/phys/dca_step/lattice_mapping/deconvolution/CMakeLists.txt
+++ b/test/unit/phys/dca_step/lattice_mapping/deconvolution/CMakeLists.txt
@@ -1,9 +1,9 @@
 # (Phys) deconvolution unit tests
 
 dca_add_gtest(deconvolution_routines_test
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   CUDA
   INCLUDE_DIRS ${SIMPLEX_GM_RULE_INCLUDE_DIR} ${FFTW_INCLUDE_DIR} ${PROJECT_SOURCE_DIR}
-  LIBS dca_io function function_transform cluster_domains  ${DCA_THREADING_LIBS} parallel_util time_and_frequency_domains
+  LIBS dca_io function function_transform cluster_domains  ${DCA_THREADING_LIBS} parallel_util parallel_no_concurrency time_and_frequency_domains
        quantum_domains gaussian_quadrature tetrahedron_mesh coarsegraining enumerations
-       ${LAPACK_LIBRARIES} lapack ${UTIL_LIBS})
+       ${LAPACK_LIBRARIES} lapack ${UTIL_LIBS} )
diff --git a/test/unit/phys/dca_step/symmetrization/CMakeLists.txt b/test/unit/phys/dca_step/symmetrization/CMakeLists.txt
index 666b36272..709a60918 100644
--- a/test/unit/phys/dca_step/symmetrization/CMakeLists.txt
+++ b/test/unit/phys/dca_step/symmetrization/CMakeLists.txt
@@ -2,7 +2,7 @@
 
 dca_add_gtest(symmetrize_test
         FAST
-        GTEST_MAIN
+        GTEST_MPI_MAIN
         INCLUDE_DIRS ${DCA_INCLUDES};${PROJECT_SOURCE_DIR};${FFTW_INCLUDE_DIR}
         LIBS  ${DCA_LIBS}
         )
diff --git a/test/unit/phys/domains/CMakeLists.txt b/test/unit/phys/domains/CMakeLists.txt
index 75e41607d..1fb2c30c6 100644
--- a/test/unit/phys/domains/CMakeLists.txt
+++ b/test/unit/phys/domains/CMakeLists.txt
@@ -2,3 +2,10 @@
 
 add_subdirectory(cluster)
 add_subdirectory(time_and_frequency)
+
+add_library(dca_domains INTERFACE)
+target_link_libraries(dca_domains INTERFACE  convert_domains
+  cluster_domains
+  quantum_domains
+  time_and_frequency_domains
+)
diff --git a/test/unit/phys/domains/cluster/CMakeLists.txt b/test/unit/phys/domains/cluster/CMakeLists.txt
index 0a7078295..821b36155 100644
--- a/test/unit/phys/domains/cluster/CMakeLists.txt
+++ b/test/unit/phys/domains/cluster/CMakeLists.txt
@@ -5,7 +5,7 @@ dca_add_gtest(cluster_operations_test
   LIBS ${LAPACK_LIBRARIES} ${DCA_GPU_LIBS})
 
 dca_add_gtest(momentum_exchage_domain_test
-  GTEST_MAIN
+  GTEST_MPI_MAIN
   LIBS ${DCA_LIBS}
   INCLUDE_DIRS ${DCA_INCLUDE_DIRS};${PROJECT_SOURCE_DIR})
 
@@ -13,6 +13,5 @@ dca_add_gtest(cluster_domain_test
   GTEST_MAIN
   LIBS ${LAPACK_LIBRARIES} function enumerations quantum_domains cluster_domains dca_io FFTW::Double ${DCA_GPU_LIBS})
 
-
 # deprecated (requires NFFT)
 # add_subdirectory(interpolation/wannier_interpolation)
diff --git a/test/unit/phys/parameters/four_point_parameters/CMakeLists.txt b/test/unit/phys/parameters/four_point_parameters/CMakeLists.txt
index c94cab018..d29861b2c 100644
--- a/test/unit/phys/parameters/four_point_parameters/CMakeLists.txt
+++ b/test/unit/phys/parameters/four_point_parameters/CMakeLists.txt
@@ -1,5 +1,6 @@
 # Four-point parameters' unit tests
 
+# why this test pulls in a bunch of magma references is a mystery
 dca_add_gtest(four_point_parameters_test
   GTEST_MAIN
-  LIBS json enumerations ${DCA_GPU_LIBS})
+  LIBS json enumerations magma::magma ${DCA_GPU_LIBS})
diff --git a/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp b/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp
index 2db2fd372..a8f19a75b 100644
--- a/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp
+++ b/test/unit/phys/parameters/output_parameters/output_parameters_test.cpp
@@ -22,15 +22,9 @@ TEST(OutputParametersTest, DefaultValues) {
 
 
   EXPECT_EQ("./", pars.get_directory());
-#ifdef DCA_HAVE_ADIOS2
-  EXPECT_EQ("ADIOS2", pars.get_output_format());
-  EXPECT_EQ("ADIOS2", pars.get_g4_output_format());
-  EXPECT_EQ("dca.bp", pars.get_filename_dca());
-#else
   EXPECT_EQ("HDF5", pars.get_output_format());
   EXPECT_EQ("", pars.get_g4_output_format());
   EXPECT_EQ("dca.hdf5", pars.get_filename_dca());
-#endif
   EXPECT_FALSE(pars.autoresume());
   EXPECT_EQ("", pars.get_directory_config_read());
   EXPECT_EQ("", pars.get_directory_config_write());