diff --git a/.circleci/config.yml b/.circleci/config.yml
index 48cc201f1..84c8625c9 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -1,9 +1,9 @@
 version: 2
 
 variables:
-  ubuntu-2004: &ubuntu-2004
+  ubuntu-2404: &ubuntu-2404
     docker:
-      - image: ghcr.io/mrchemsoft/metamr/circleci_ubuntu-20.04:sha-343e011
+      - image: ghcr.io/mrchemsoft/metamr/circleci_ubuntu-24.04:sha-2b5887e
         name: tsubame
         user: merzbow
     working_directory: ~/mrcpp
@@ -55,14 +55,14 @@ variables:
 
 jobs:
   serial-py3:
-    <<: *ubuntu-2004
+    <<: *ubuntu-2404
     steps:
       - checkout
       - *configure-serial
       - *build
       - *tests
   omp-py3:
-    <<: *ubuntu-2004
+    <<: *ubuntu-2404
     environment:
       - OMP_NUM_THREADS: '2'
     steps:
@@ -71,7 +71,7 @@ jobs:
       - *build
       - *tests
   mpi-py3:
-    <<: *ubuntu-2004
+    <<: *ubuntu-2404
     steps:
       - checkout
       - *configure-mpi
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b9c20d8b8..389a481b6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -2,11 +2,16 @@
 # Copyright (c) 2015-2020 by Radovan Bast, Roberto Di Remigio, Jonas Juselius, and contributors.
 
 # set minimum cmake version
-cmake_minimum_required(VERSION 3.14 FATAL_ERROR)
+cmake_minimum_required(VERSION 3.22 FATAL_ERROR)
 
 # project name
 project(MRCPP LANGUAGES CXX)
 
+# Suppress warning for policy CMP0177
+if(CMP0177)
+    cmake_policy(SET CMP0177 OLD)
+endif()
+
 # do not rebuild if rules (compiler flags) change
 set(CMAKE_SKIP_RULE_DEPENDENCY TRUE)
 
diff --git a/cmake/compiler_flags/CXXFlags.cmake b/cmake/compiler_flags/CXXFlags.cmake
index cbcf32898..7414eac5c 100644
--- a/cmake/compiler_flags/CXXFlags.cmake
+++ b/cmake/compiler_flags/CXXFlags.cmake
@@ -18,24 +18,28 @@ set(CMAKE_CXX_STANDARD_REQUIRED TRUE)
 set(CMAKE_CXX_EXTENSIONS FALSE)
 set(CMAKE_EXPORT_COMPILE_COMMANDS TRUE)
 
-if(ENABLE_ARCH_FLAGS)
-  # iterate over list of flags and use the first one that is compatible with the
-  # compiler in use
-  if(CMAKE_CXX_COMPILER_ID MATCHES GNU)
-    set(_arch_flag "-march=native")
+# Enforce minimum compiler versions for C++17 support
+if(CMAKE_CXX_COMPILER_ID MATCHES GNU)
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "11.2")
+    message(FATAL_ERROR "GCC >= 11.2 required, found ${CMAKE_CXX_COMPILER_VERSION}")
   endif()
-  if(CMAKE_CXX_COMPILER_ID MATCHES Clang)
-    set(_arch_flag "-march=native")
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "14.0")
+    message(FATAL_ERROR "Clang >= 14.0 required, found ${CMAKE_CXX_COMPILER_VERSION}")
   endif()
-  if(CMAKE_CXX_COMPILER_ID MATCHES Intel)
-    set(_arch_flag "-xHost")
+elseif(CMAKE_CXX_COMPILER_ID MATCHES "IntelLLVM")
+  if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "2022.1")
+    message(FATAL_ERROR "Intel oneAPI (icx) >= 2022.1 required, found ${CMAKE_CXX_COMPILER_VERSION}")
   endif()
+endif()
+
+if(ENABLE_ARCH_FLAGS)
+  set(_arch_flag "-march=native")
   message(STATUS "Adding architecture-specific compiler flag: ${_arch_flag}")
   set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${_arch_flag}")
 endif()
 
 include(${CMAKE_CURRENT_LIST_DIR}/GNU.CXX.cmake)
-include(${CMAKE_CURRENT_LIST_DIR}/Intel.CXX.cmake)
 include(${CMAKE_CURRENT_LIST_DIR}/Clang.CXX.cmake)
 
 string(REPLACE " " ";" _cmake_cxx_flags ${CMAKE_CXX_FLAGS})
diff --git a/cmake/compiler_flags/Intel.CXX.cmake b/cmake/compiler_flags/Intel.CXX.cmake
deleted file mode 100644
index 98ff603f2..000000000
--- a/cmake/compiler_flags/Intel.CXX.cmake
+++ /dev/null
@@ -1,8 +0,0 @@
-if(NOT DEFINED ENV{CXXFLAGS})
-    if(CMAKE_CXX_COMPILER_ID MATCHES Intel)
-        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -ftz -fp-speculation=fast -fp-model fast -Wno-unknown-pragmas")
-        set(CMAKE_CXX_FLAGS_RELEASE "-O2 -DNDEBUG")
-        set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O2 -debug -DNDEBUG")
-        set(CMAKE_CXX_FLAGS_DEBUG "-O0 -debug -DDEBUG")
-    endif()
-endif()
diff --git a/docs/install.rst b/docs/install.rst
index 2e37501c3..196a598a3 100644
--- a/docs/install.rst
+++ b/docs/install.rst
@@ -1,42 +1,60 @@
+============
+Installation
+============
+
 ------------------
 Obtaining the code
 ------------------
 
-The latest version of MRCPP is available on `GitHub
-<http://github.com/MRChemSoft/mrcpp>`_::
+The latest development version of MRCPP can be found on the ``master``
+branch on GitHub::
 
-    $ git clone git@github.com:MRChemSoft/mrcpp.git
+    $ git clone https://github.com/MRChemSoft/mrcpp.git
 
+The released versions can be found from Git tags ``vX.Y.Z`` under the
+``release/X.Y`` branches in the same repository, or a zip file can be
+downloaded from `Zenodo <https://doi.org/10.5281/zenodo.3606670>`_.
 
------------------
-Building the code
------------------
+By default, all dependencies will be **fetched** at configure time if they are
+not already available.
 
 
-Prerequisites
--------------
+-------------------
+Build prerequisites
+-------------------
 
-* g++-5.4 or later (``std=c++14``)
-* `CMake <http://cmake.org>`_ version 3.11 or higher.
-* `Eigen <http://eigen.tuxfamily.org>`_ version 3.3 or higher.
-* BLAS (optional)
+- CMake-3.22 (or later)
+- GNU-11.2, Clang-14.0 or IntelLLVM-2022.1 (or later) compilers (C++17 standard)
 
+.. hint::
+    We have collected the recommended modules for the different Norwegian HPC
+    systems under ``tools/<machine>.env``. These files can be sourced in order
+    to get a working environment on the respective machines, and may also serve
+    as a guide for other HPC systems.
+  
 
-Configuration
--------------
+C++ dependencies
+----------------
 
-The configuration and build process is managed through CMake, and a ``setup``
-script is provided for the configuration step. MRCPP's only dependency is
-Eigen3, which will be downloaded at configure time unless it is already
-available on the system. If you have a local version *not* under the system
-path, you can point to it before running ``setup``::
+- Linear algebra: `Eigen-3.4  <https://gitlab.com/libeigen/eigen>`_
+- BLAS (optional)
 
-    $ export EIGEN3_ROOT=/path/to/eigen3
-    $ ./setup [options] [<builddir>]
+Eigen will be downloaded automatically at configure time by CMake,
+but can also be linked manually by setting the variable::
 
-The setup script will create a directory called *<builddir>* (default ``build``)
-and run CMake. There are several options available for the setup, and the most
-important are:
+    EIGEN3_DIR=<path_to_eigen3>/share/eigen3/cmake
+
+
+-----------------
+Building the code
+-----------------
+
+Configure
+---------
+
+The ``setup`` script will create a directory called ``<build-dir>`` and run
+CMake. There are several options available for the setup, the most
+important being:
 
 ``--cxx=<CXX>``
   C++ compiler [default: g++]
@@ -55,31 +73,103 @@ important are:
 ``-h --help``
   List all options
 
+The code can be built with four levels of parallelization:
 
-Compilation
------------
+ - no parallelization
+ - only shared memory (OpenMP)
+ - only distributed memory (MPI)
+ - hybrid OpenMP + MPI
 
-After successful configuration, the code is compiled using the ``make`` command
-in the *<builddir>* directory::
+.. note::
+    In practice we recommend the **shared memory version** for running on your
+    personal laptop/workstation, and the **hybrid version** for running on a
+    HPC cluster. The serial and pure MPI versions are only useful for debugging.
 
-    $ cd <builddir>
-    $ make
+The default build is *without* parallelization and using GNU compilers::
+
+    $ ./setup --prefix=<install-dir> <build-dir>
+
+To use clang compilers you need to specify the ``--cxx`` option::
+
+    $ ./setup --prefix=<install-dir> --cxx=clang++ <build-dir>
+
+To build the code with shared memory (OpenMP) parallelization,
+add the ``--omp`` option::
 
+    $ ./setup --prefix=<install-dir> --omp <build-dir>
 
--------------
-Running tests
--------------
+To build the code with distributed memory (MPI) parallelization, add the
+``--mpi`` option *and* change to the respective MPI compilers (``--cxx=mpicxx``
+for GNU)::
 
-A set of tests is provided with the code to verify that the code compiled
-properly. To compile the test suite, add the ``--enable-tests`` option to
-setup, then run the tests with ``ctest``::
+    $ ./setup --prefix=<install-dir> --omp --mpi --cxx=mpicxx <build-dir>
 
-    $ ./setup --enable-tests build
-    $ cd build
+.. note::
+    If you compile the MRCPP library manually as a separate project, the level
+    of parallelization **must be the same** for MRCPP and MRChem. Similar
+    options apply for the MRCPP setup, see
+    `mrcpp.readthedocs.io <https://mrcpp.readthedocs.io/en/latest/>`_.
+
+
+Build
+-----
+
+If the CMake configuration is successful, the code is compiled with::
+
+    $ cd <build-dir>
     $ make
+
+
+Test
+----
+
+A test suite is provided to make sure that everything compiled properly.
+To run a collection of small tests::
+
+    $ cd <build-dir>
     $ ctest
 
 
+Install
+-------
+
+After the build has been verified with the test suite, it can be installed with
+the following command::
+
+    $ cd <build-dir>
+    $ make install
+
+Now libraries, headers and CMake configuration files can be found under the
+given prefix::
+
+    mrcpp/
+    ├── include/
+    │   └── MRCPP/
+    ├── lib64/
+    │   ├── libmrcpp.a
+    │   ├── libmrcpp.so -> libmrcpp.so.1*
+    │   └── libmrcpp.so.1*
+    └── share/
+        └── cmake/
+
+Please refer to the :ref:`User's Manual` for instructions for how to run the program.
+
+.. hint::
+    We have collected scripts for configure and build of the hybrid OpenMP + MPI
+    version on the different Norwegian HPC systems under ``tools/<machine>.sh``.
+    These scripts will build the current version under ``build-${version}``,
+    run the unit tests and install under ``install-${version}``, e.g. to build
+    version v1.5.0 on Olivia::
+
+        $ cd mrcpp
+        $ git checkout v1.5.0
+        $ tools/olivia.sh
+
+    The configure step requires internet access, so the scripts must be run on
+    the login nodes, and it will run on a single core, so it might take some
+    minutes to complete. 
+
+
 ----------------
 Running examples
 ----------------
@@ -127,6 +217,7 @@ Note that the core of MRCPP is *only* OpenMP parallelized. All MPI data or work
 distribution must be done manually in the application program, using the tools
 provided by MRCPP (see the Parallel section of the API).
 
+
 ----------
 Pilot code
 ----------
@@ -144,32 +235,6 @@ Feel free to do whatever you like in your own pilot code, but please don't add
 this file to git. Also, please don't commit any changes to the existing examples
 (unless you know what you're doing).
 
----------------------
-MRCPP as a dependency
----------------------
-
-Building MRCPP provides CMake configuration files exporting the libraries and
-headers as targets to be consumed by third-party projects also using CMake::
-
-    $ ./setup --prefix=$HOME/Software/mrcpp
-    $ cd build
-    $ make
-    $ ctest
-    $ make install
-
-Now libraries, headers and CMake configuration files can be found under the
-given prefix::
-
-    mrcpp/
-    ├── include/
-    │   └── MRCPP/
-    ├── lib64/
-    │   ├── libmrcpp.a
-    │   ├── libmrcpp.so -> libmrcpp.so.1*
-    │   └── libmrcpp.so.1*
-    └── share/
-        └── cmake/
-
 As an example, the ``pilot`` sample can be built with the following ``CMakeLists.txt``:
 
 .. literalinclude:: snippets/CMakeLists.txt
diff --git a/external/upstream/fetch_eigen3.cmake b/external/upstream/fetch_eigen3.cmake
index 68a6683c8..cd68fcd20 100644
--- a/external/upstream/fetch_eigen3.cmake
+++ b/external/upstream/fetch_eigen3.cmake
@@ -1,11 +1,13 @@
-find_package(Eigen3 CONFIG QUIET
-  NO_CMAKE_PATH
-  NO_CMAKE_PACKAGE_REGISTRY
-  )
-
+if(EIGEN3_FIND_BEHAVIOUR STREQUAL "default" OR EIGEN3_FIND_BEHAVIOUR STREQUAL "onlylocal")
+  find_package(Eigen3 3.4 CONFIG QUIET
+    NO_CMAKE_PATH
+    NO_CMAKE_PACKAGE_REGISTRY
+    NO_CMAKE_SYSTEM_PACKAGE_REGISTRY
+    )
+endif()
 if(TARGET Eigen3::Eigen)
-  message(STATUS "Using Eigen3: ${EIGEN3_ROOT_DIR} (version ${Eigen3_VERSION})")
-else()
+  message(STATUS "Found Eigen3: ${EIGEN3_ROOT_DIR} (version ${Eigen3_VERSION})")
+elseif(EIGEN3_FIND_BEHAVIOUR STREQUAL "default" OR EIGEN3_FIND_BEHAVIOUR STREQUAL "onlyfetch")
   message(STATUS "Suitable Eigen3 could not be located. Fetching and building!")
   include(FetchContent)
 
@@ -13,8 +15,11 @@ else()
     QUIET
     URL
       https://gitlab.com/libeigen/eigen/-/archive/3.4.0/eigen-3.4.0.tar.gz
+    DOWNLOAD_EXTRACT_TIMESTAMP TRUE
     )
 
   set(BUILD_TESTING OFF CACHE BOOL "" FORCE)
   FetchContent_MakeAvailable(eigen3)
+else()
+  message(FATAL_ERROR "No suitable Eigen3 found or fetched. Aborting setup!")
 endif()
diff --git a/setup b/setup
index 259e89b15..7ee86fedc 100755
--- a/setup
+++ b/setup
@@ -5,7 +5,7 @@
 
 import os
 import sys
-assert sys.version_info >= (2, 6), 'Python >= 2.6 is required'
+assert sys.version_info >= (3, 10), 'Python >= 3.10 is required'
 
 sys.path.insert(0, 'cmake')
 from autocmake import configure
@@ -27,6 +27,7 @@ Options:
   --static                               Create only the static library [default: False].
   --enable-tests=<TESTS>                 Enable tests [default: True].
   --enable-examples                      Enable code examples [default: False].
+  --eigen3=<STRING>                      Behaviour for finding/fetching Eigen3 (default, onlylocal, onlyfetch) [default: default].
   --type=<TYPE>                          Set the CMake build type (debug, release, relwithdebinfo, minsizerel) [default: release].
   --generator=<STRING>                   Set the CMake build system generator [default: Unix Makefiles].
   --show                                 Show CMake command and exit.
@@ -52,6 +53,7 @@ def gen_cmake_command(options, arguments):
     command.append('-DBUILD_STATIC_LIBS={0}'.format(arguments['--static']))
     command.append('-DENABLE_TESTS={0}'.format(arguments['--enable-tests']))
     command.append('-DENABLE_EXAMPLES={0}'.format(arguments['--enable-examples']))
+    command.append('-DEIGEN3_FIND_BEHAVIOUR={0}'.format(arguments['--eigen3']))
     command.append('-DCMAKE_BUILD_TYPE={0}'.format(arguments['--type']))
     command.append('-G"{0}"'.format(arguments['--generator']))
     if arguments['--cmake-options'] != "''":
diff --git a/src/core/GaussQuadrature.cpp b/src/core/GaussQuadrature.cpp
index 9c30823f8..94e584453 100644
--- a/src/core/GaussQuadrature.cpp
+++ b/src/core/GaussQuadrature.cpp
@@ -247,6 +247,8 @@ double GaussQuadrature::integrate(RepresentableFunction<3> &func) const {
  * This function has been implemented using a recursive algorithm.
  */
 double GaussQuadrature::integrate_nd(RepresentableFunction<3> &func, int axis) const {
+    (void)func;
+    (void)axis;
     NOT_IMPLEMENTED_ABORT;
     NEEDS_TESTING
 
diff --git a/src/core/ObjectCache.cpp b/src/core/ObjectCache.cpp
index 0df714200..9a56f57c2 100644
--- a/src/core/ObjectCache.cpp
+++ b/src/core/ObjectCache.cpp
@@ -44,6 +44,7 @@ template <class T> void ObjectCache<T>::clear() {
 }
 
 template <class T> void ObjectCache<T>::load(int id) {
+    (void)id;
     MSG_INFO("This routine does nothing in this class.");
 }
 
diff --git a/src/core/ScalingBasis.cpp b/src/core/ScalingBasis.cpp
index 99dc83f24..8bb2fdd94 100644
--- a/src/core/ScalingBasis.cpp
+++ b/src/core/ScalingBasis.cpp
@@ -39,7 +39,7 @@ ScalingBasis::ScalingBasis(int k, int t)
 }
 
 void ScalingBasis::evalf(const double *r, Eigen::MatrixXd &vals) const {
-    if (vals.rows() != this->funcs.size()) MSG_ERROR("Invalid argument");
+    if (static_cast<size_t>(vals.rows()) != this->funcs.size()) MSG_ERROR("Invalid argument");
 
     for (int d = 0; d < vals.cols(); d++) {
         for (int k = 0; k < vals.rows(); k++) { vals(k, d) = getFunc(k).evalf(r[d]); }
diff --git a/src/functions/GaussExp.cpp b/src/functions/GaussExp.cpp
index 21432686d..5874e9e71 100644
--- a/src/functions/GaussExp.cpp
+++ b/src/functions/GaussExp.cpp
@@ -45,16 +45,16 @@ template <int D> GaussExp<D>::GaussExp(int nTerms) {
     for (int i = 0; i < nTerms; i++) { this->funcs.push_back(nullptr); }
 }
 
-template <int D> GaussExp<D>::GaussExp(const GaussExp<D> &gexp) {
+template <int D> GaussExp<D>::GaussExp(const GaussExp<D> &gexp) : RepresentableFunction<D, double>(gexp) {
     screening = gexp.screening;
-    for (unsigned int i = 0; i < gexp.size(); i++) {
+    for (size_t i = 0; i < gexp.size(); i++) {
         Gaussian<D> *gauss = gexp.funcs[i]->copy();
         this->funcs.push_back(gauss);
     }
 }
 
 template <int D> GaussExp<D>::~GaussExp() {
-    for (int i = 0; i < size(); i++) {
+    for (size_t i = 0; i < size(); i++) {
         if (this->funcs[i] != nullptr) {
             delete this->funcs[i];
             this->funcs[i] = nullptr;
@@ -66,7 +66,7 @@ template <int D> GaussExp<D> &GaussExp<D>::operator=(const GaussExp<D> &gexp) {
     if (&gexp == this) return *this;
     // screening = gexp.screening;
     this->funcs.clear();
-    for (unsigned int i = 0; i < gexp.size(); i++) {
+    for (size_t i = 0; i < gexp.size(); i++) {
         if (gexp.funcs[i] == nullptr) {
             this->funcs.push_back(nullptr);
         } else {
@@ -79,26 +79,26 @@ template <int D> GaussExp<D> &GaussExp<D>::operator=(const GaussExp<D> &gexp) {
 
 template <int D> double GaussExp<D>::evalf(const Coord<D> &r) const {
     double val = 0.0;
-    for (int i = 0; i < this->size(); i++) { val += this->getFunc(i).evalf(r); }
+    for (size_t i = 0; i < this->size(); i++) { val += this->getFunc(i).evalf(r); }
     return val;
 }
 
 template <int D> bool GaussExp<D>::isVisibleAtScale(int scale, int nPts) const {
-    for (unsigned int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         if (not this->getFunc(i).isVisibleAtScale(scale, nPts)) { return false; }
     }
     return true;
 }
 
 template <int D> bool GaussExp<D>::isZeroOnInterval(const double *lb, const double *ub) const {
-    for (unsigned int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         if (not this->getFunc(i).isZeroOnInterval(lb, ub)) { return false; }
     }
     return true;
 }
 
 template <int D> void GaussExp<D>::setFunc(int i, const GaussPoly<D> &g, double c) {
-    if (i < 0 or i > (this->size() - 1)) {
+    if (i < 0 or i > (static_cast<int>(this->size()) - 1)) {
         MSG_ERROR("Index out of bounds!");
         return;
     }
@@ -109,7 +109,7 @@ template <int D> void GaussExp<D>::setFunc(int i, const GaussPoly<D> &g, double
 }
 
 template <int D> void GaussExp<D>::setFunc(int i, const GaussFunc<D> &g, double c) {
-    if (i < 0 or i > (this->size() - 1)) {
+    if (i < 0 or i > (static_cast<int>(this->size()) - 1)) {
         MSG_ERROR("Index out of bounds!");
         return;
     }
@@ -125,7 +125,7 @@ template <int D> void GaussExp<D>::append(const Gaussian<D> &g) {
 }
 
 template <int D> void GaussExp<D>::append(const GaussExp<D> &g) {
-    for (int i = 0; i < g.size(); i++) {
+    for (size_t i = 0; i < g.size(); i++) {
         Gaussian<D> *gp = g.getFunc(i).copy();
         this->funcs.push_back(gp);
     }
@@ -134,7 +134,7 @@ template <int D> void GaussExp<D>::append(const GaussExp<D> &g) {
 template <int D> GaussExp<D> GaussExp<D>::differentiate(int dir) const {
     assert(dir >= 0 and dir < D);
     GaussExp<D> result;
-    for (int i = 0; i < this->size(); i++) result.append(this->getFunc(i).differentiate(dir));
+    for (size_t i = 0; i < this->size(); i++) result.append(this->getFunc(i).differentiate(dir));
     return result;
 }
 
@@ -143,11 +143,11 @@ template <int D> GaussExp<D> GaussExp<D>::add(GaussExp<D> &g) {
     GaussExp<D> sum = GaussExp<D>(nsum);
 
     int n = 0;
-    for (int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         sum.funcs[n] = this->funcs[i]->copy();
         n++;
     }
-    for (int i = 0; i < g.size(); i++) {
+    for (size_t i = 0; i < g.size(); i++) {
         sum.funcs[n] = g.funcs[i]->copy();
         n++;
     }
@@ -158,15 +158,15 @@ template <int D> GaussExp<D> GaussExp<D>::add(GaussExp<D> &g) {
 template <int D> GaussExp<D> GaussExp<D>::add(Gaussian<D> &g) {
     int nsum = this->size() + 1;
     GaussExp<D> sum = GaussExp<D>(nsum);
-    for (int n = 0; n < this->size(); n++) { sum.funcs[n] = this->getFunc(n).copy(); }
+    for (size_t n = 0; n < this->size(); n++) { sum.funcs[n] = this->getFunc(n).copy(); }
     sum.funcs[this->size()] = g.copy();
     return sum;
 }
 
 template <int D> GaussExp<D> GaussExp<D>::mult(GaussExp<D> &gexp) {
     GaussExp<D> result;
-    for (int i = 0; i < this->size(); i++) {
-        for (int j = 0; j < gexp.size(); j++) {
+    for (size_t i = 0; i < this->size(); i++) {
+        for (size_t j = 0; j < gexp.size(); j++) {
             if (auto *f = dynamic_cast<GaussFunc<D> *>(this->funcs[i])) {
                 if (auto *g = dynamic_cast<GaussFunc<D> *>(gexp.funcs[j])) {
                     GaussPoly<D> newTerm = (*g) * (*f);
@@ -230,29 +230,29 @@ template <int D> GaussExp<D> GaussExp<D>::mult(GaussPoly<D> &g) {
 
 template <int D> GaussExp<D> GaussExp<D>::mult(double d) {
     GaussExp<D> prod = *this;
-    for (int i = 0; i < this->size(); i++) prod.funcs[i]->multConstInPlace(d);
+    for (size_t i = 0; i < this->size(); i++) prod.funcs[i]->multConstInPlace(d);
     return prod;
 }
 
 template <int D> void GaussExp<D>::multInPlace(double d) {
-    for (int i = 0; i < this->size(); i++) this->funcs[i]->multConstInPlace(d);
+    for (size_t i = 0; i < this->size(); i++) this->funcs[i]->multConstInPlace(d);
 }
 
 template <int D> double GaussExp<D>::calcSquareNorm() const {
     /* computing the squares */
     double norm = 0.0;
-    for (int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         double nc = this->funcs[i]->calcSquareNorm();
         norm += nc;
     }
     /* computing the double products */
-    for (int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         GaussExp<D> funcs_i = getFunc(i).asGaussExp(); // Make sure all entries are GaussFunc
-        for (int fi = 0; fi < funcs_i.size(); fi++) {
+        for (size_t fi = 0; fi < funcs_i.size(); fi++) {
             GaussFunc<D> &func_i = static_cast<GaussFunc<D> &>(funcs_i.getFunc(fi));
-            for (int j = i + 1; j < this->size(); j++) {
+            for (size_t j = i + 1; j < this->size(); j++) {
                 GaussExp<D> funcs_j = getFunc(j).asGaussExp(); // Make sure all entries are GaussFunc
-                for (int fj = 0; fj < funcs_j.size(); fj++) {
+                for (size_t fj = 0; fj < funcs_j.size(); fj++) {
                     GaussFunc<D> &func_j = static_cast<GaussFunc<D> &>(funcs_j.getFunc(fj));
                     double overlap = func_i.calcOverlap(func_j);
                     norm += 2.0 * overlap;
@@ -265,7 +265,7 @@ template <int D> double GaussExp<D>::calcSquareNorm() const {
 
 template <int D> void GaussExp<D>::normalize() {
     double norm = std::sqrt(this->calcSquareNorm());
-    for (int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         double coef = this->funcs[i]->getCoef();
         this->funcs[i]->setCoef(coef / norm);
     }
@@ -273,7 +273,7 @@ template <int D> void GaussExp<D>::normalize() {
 
 template <int D> void GaussExp<D>::calcScreening(double nStdDev) {
     screening = nStdDev;
-    for (int i = 0; i < this->size(); i++) { this->funcs[i]->calcScreening(nStdDev); }
+    for (size_t i = 0; i < this->size(); i++) { this->funcs[i]->calcScreening(nStdDev); }
 }
 
 template <int D> void GaussExp<D>::setScreen(bool screen) {
@@ -282,7 +282,7 @@ template <int D> void GaussExp<D>::setScreen(bool screen) {
     } else {
         this->screening = -std::abs(this->screening);
     }
-    for (int i = 0; i < this->size(); i++) { this->funcs[i]->setScreen(screen); }
+    for (size_t i = 0; i < this->size(); i++) { this->funcs[i]->setScreen(screen); }
 }
 
 // Calculate the scaling and wavelet coefs of all the children, and do the
@@ -326,7 +326,7 @@ template <int D> void GaussExp<D>::setDefaultScreening(double screen) {
 
 template <int D> std::ostream &GaussExp<D>::print(std::ostream &o) const {
     o << "Gaussian expansion: " << size() << " terms" << std::endl;
-    for (int i = 0; i < size(); i++) {
+    for (size_t i = 0; i < size(); i++) {
         o << "Term" << std::setw(3) << i << " :" << std::endl;
         o << getFunc(i) << std::endl << std::endl;
     }
@@ -344,13 +344,13 @@ template <int D> double GaussExp<D>::calcCoulombEnergy() const {
 
 template <> double GaussExp<3>::calcCoulombEnergy() const {
     double energy = 0.0;
-    for (int i = 0; i < this->size(); i++) {
+    for (size_t i = 0; i < this->size(); i++) {
         GaussExp<3> funcs_i = getFunc(i).asGaussExp(); // Make sure all entries are GaussFunc
-        for (int fi = 0; fi < funcs_i.size(); fi++) {
+        for (size_t fi = 0; fi < funcs_i.size(); fi++) {
             GaussFunc<3> &func_i = static_cast<GaussFunc<3> &>(funcs_i.getFunc(fi));
-            for (int j = i; j < this->size(); j++) {
+            for (size_t j = i; j < this->size(); j++) {
                 GaussExp<3> funcs_j = getFunc(j).asGaussExp(); // Make sure all entries are GaussFunc
-                for (int fj = 0; fj < funcs_j.size(); fj++) {
+                for (size_t fj = 0; fj < funcs_j.size(); fj++) {
                     GaussFunc<3> &func_j = static_cast<GaussFunc<3> &>(funcs_j.getFunc(fj));
                     double c = 2.0;
                     if (i == j) c = 1.0;
diff --git a/src/functions/GaussExp.h b/src/functions/GaussExp.h
index 382213214..3e291ea93 100644
--- a/src/functions/GaussExp.h
+++ b/src/functions/GaussExp.h
@@ -95,7 +95,7 @@ template <int D> class GaussExp : public RepresentableFunction<D, double> {
     const std::array<int, D> &getPower(int i) const { return this->funcs[i]->getPower(); }
     const std::array<double, D> &getPos(int i) const { return this->funcs[i]->getPos(); }
 
-    int size() const { return this->funcs.size(); }
+    size_t size() const { return this->funcs.size(); }
     Gaussian<D> &getFunc(int i) { return *this->funcs[i]; }
     const Gaussian<D> &getFunc(int i) const { return *this->funcs[i]; }
 
diff --git a/src/functions/GaussFunc.cpp b/src/functions/GaussFunc.cpp
index 28736be58..57d49008c 100644
--- a/src/functions/GaussFunc.cpp
+++ b/src/functions/GaussFunc.cpp
@@ -204,6 +204,7 @@ template <int D> std::ostream &GaussFunc<D>::print(std::ostream &o) const {
  *  \f$ \alpha = (\beta/\pi)^{D/2} \f$ for this to be correct!
  */
 template <int D> double GaussFunc<D>::calcCoulombEnergy(const GaussFunc<D> &gf) const {
+    (void)gf;
     NOT_IMPLEMENTED_ABORT;
 }
 
diff --git a/src/functions/GaussPoly.cpp b/src/functions/GaussPoly.cpp
index 0dfeaf2cd..145c06e3a 100644
--- a/src/functions/GaussPoly.cpp
+++ b/src/functions/GaussPoly.cpp
@@ -101,9 +101,9 @@ template <int D> Gaussian<D> *GaussPoly<D>::copy() const {
 template<int D> double GaussPoly<D>::calcSquareNorm() const {
     GaussExp<D> this_exp = this->asGaussExp();
     double norm = 0.0;
-    for (int i = 0; i < this_exp.size(); i++) {
+    for (size_t i = 0; i < this_exp.size(); i++) {
         auto func_i = static_cast<GaussFunc<D> &>(this_exp.getFunc(i));
-        for (int j = 0; j < this_exp.size(); j++) {
+        for (size_t j = 0; j < this_exp.size(); j++) {
             auto func_j = static_cast<GaussFunc<D> &>(this_exp.getFunc(j));
             norm += function_utils::calc_overlap(func_i, func_j);
         }
@@ -176,10 +176,12 @@ template <int D> GaussExp<D> GaussPoly<D>::asGaussExp() const {
 }
 
 template <int D> GaussPoly<D> GaussPoly<D>::differentiate(int dir) const {
+    (void)dir;
     NOT_IMPLEMENTED_ABORT;
 }
 
 template <int D> void GaussPoly<D>::multInPlace(const GaussPoly<D> &rhs) {
+    (void)rhs;
     NOT_IMPLEMENTED_ABORT;
 }
 
@@ -229,6 +231,7 @@ void GaussPoly<D>::fillCoefPowVector(std::vector<double> &coefs,
 }
 
 template <int D> GaussPoly<D> GaussPoly<D>::mult(const GaussPoly<D> &rhs) {
+    (void)rhs;
     NOT_IMPLEMENTED_ABORT;
     /*
     GaussPoly<D> &lhs = *this;
diff --git a/src/functions/Gaussian.cpp b/src/functions/Gaussian.cpp
index 0b8e3214b..a5b8e2e5c 100644
--- a/src/functions/Gaussian.cpp
+++ b/src/functions/Gaussian.cpp
@@ -161,9 +161,9 @@ template <int D> double Gaussian<D>::calcOverlap(const Gaussian<D> &inp) const {
     const auto &ket_exp = inp.asGaussExp();   // Make sure all entries are GaussFunc
 
     double S = 0.0;
-    for (int i = 0; i < bra_exp.size(); i++) {
+    for (size_t i = 0; i < bra_exp.size(); i++) {
         const auto &bra_i = static_cast<const GaussFunc<D> &>(bra_exp.getFunc(i));
-        for (int j = 0; j < ket_exp.size(); j++) {
+        for (size_t j = 0; j < ket_exp.size(); j++) {
             const auto &ket_j = static_cast<const GaussFunc<D> &>(ket_exp.getFunc(j));
             S += function_utils::calc_overlap(bra_i, ket_j);
         }
diff --git a/src/functions/JpowerIntegrals.cpp b/src/functions/JpowerIntegrals.cpp
index 179f6fcc6..0fba5c8b8 100644
--- a/src/functions/JpowerIntegrals.cpp
+++ b/src/functions/JpowerIntegrals.cpp
@@ -46,6 +46,8 @@ std::vector<std::complex<double>> &JpowerIntegrals::operator[](int index) {
 std::vector<std::complex<double>> JpowerIntegrals::calculate_J_power_integrals(int l, double a, int M, double threshold) {
     using namespace std::complex_literals;
 
+    (void)threshold;
+
     std::complex<double> J_0 = 0.25 * std::exp(-0.25i * M_PI) / std::sqrt(M_PI * a) * std::exp(0.25i * static_cast<double>(l * l) / a);
     std::complex<double> beta(0, 0.5 / a);
     auto alpha = static_cast<double>(l) * beta;
diff --git a/src/functions/Polynomial.h b/src/functions/Polynomial.h
index 93e3ec77d..c5705845e 100644
--- a/src/functions/Polynomial.h
+++ b/src/functions/Polynomial.h
@@ -76,7 +76,7 @@ class Polynomial : public RepresentableFunction<1, double> {
     void dilate(double n) { this->N *= n; }
     void translate(double l) { this->L += this->N * l; }
 
-    int size() const { return this->coefs.size(); } ///< Length of coefs vector
+    size_t size() const { return this->coefs.size(); } ///< Length of coefs vector
     int getOrder() const;
     void clearCoefs() { this->coefs = Eigen::VectorXd::Zero(1); }
     void setZero() { this->coefs = Eigen::VectorXd::Zero(this->coefs.size()); }
diff --git a/src/functions/RepresentableFunction.cpp b/src/functions/RepresentableFunction.cpp
index 3c55ac92b..765af1708 100644
--- a/src/functions/RepresentableFunction.cpp
+++ b/src/functions/RepresentableFunction.cpp
@@ -75,6 +75,7 @@ template <int D, typename T> RepresentableFunction<D, T>::RepresentableFunction(
 /** Copies function, not bounds. Use copy constructor if you want an
  * identical function. */
 template <int D, typename T> RepresentableFunction<D, T> &RepresentableFunction<D, T>::operator=(const RepresentableFunction<D, T> &func) {
+    (void)func;
     return *this;
 }
 
diff --git a/src/functions/RepresentableFunction.h b/src/functions/RepresentableFunction.h
index 6123e3051..97ec7d509 100644
--- a/src/functions/RepresentableFunction.h
+++ b/src/functions/RepresentableFunction.h
@@ -73,8 +73,17 @@ template <int D, typename T> class RepresentableFunction {
     double *A; ///< Lower bound, NULL if unbounded
     double *B; ///< Upper bound, Null if unbounded
 
-    virtual bool isVisibleAtScale(int scale, int nQuadPts) const { return true; }
-    virtual bool isZeroOnInterval(const double *a, const double *b) const { return false; }
+    virtual bool isVisibleAtScale(int scale, int nQuadPts) const {
+        (void)scale;
+        (void)nQuadPts;
+        return true;
+    }
+
+    virtual bool isZeroOnInterval(const double *a, const double *b) const {
+        (void)a;
+        (void)b;
+        return false;
+    }
 };
 
 /*
diff --git a/src/operators/ConvolutionOperator.cpp b/src/operators/ConvolutionOperator.cpp
index 9d37929aa..11728f82d 100644
--- a/src/operators/ConvolutionOperator.cpp
+++ b/src/operators/ConvolutionOperator.cpp
@@ -82,7 +82,7 @@ template <int D> void ConvolutionOperator<D>::initialize(GaussExp<1> &kernel, do
     TreeBuilder<2> builder;
     OperatorAdaptor adaptor(o_prec, o_mra.getMaxScale());
 
-    for (int i = 0; i < kernel.size(); i++) {
+    for (size_t i = 0; i < kernel.size(); i++) {
         // Rescale Gaussian for D-dim application
         auto *k_func = kernel.getFunc(i).copy();
         k_func->setCoef(std::copysign(std::pow(std::abs(k_func->getCoef()), 1.0 / D), k_func->getCoef()));
diff --git a/src/operators/MWOperator.cpp b/src/operators/MWOperator.cpp
index 225108f48..6a45586f9 100644
--- a/src/operators/MWOperator.cpp
+++ b/src/operators/MWOperator.cpp
@@ -33,7 +33,7 @@ using namespace Eigen;
 namespace mrcpp {
 
 template <int D> void MWOperator<D>::initOperExp(int M) {
-    if (this->raw_exp.size() < M) MSG_ABORT("Incompatible raw expansion");
+    if (this->raw_exp.size() < static_cast<size_t>(M)) MSG_ABORT("Incompatible raw expansion");
     this->oper_exp.clear();
     for (int m = 0; m < M; m++) {
         std::array<OperatorTree *, D> otrees;
@@ -47,14 +47,14 @@ template <int D> void MWOperator<D>::initOperExp(int M) {
 }
 
 template <int D> OperatorTree &MWOperator<D>::getComponent(int i, int d) {
-    if (i < 0 or i >= this->oper_exp.size()) MSG_ERROR("Index out of bounds");
+    if (i < 0 or static_cast<size_t>(i) >= this->oper_exp.size()) MSG_ERROR("Index out of bounds");
     if (d < 0 or d >= D) MSG_ERROR("Dimension out of bounds");
     if (this->oper_exp[i][d] == nullptr) MSG_ERROR("Invalid component");
     return *this->oper_exp[i][d];
 }
 
 template <int D> const OperatorTree &MWOperator<D>::getComponent(int i, int d) const {
-    if (i < 0 or i >= this->oper_exp.size()) MSG_ERROR("Index out of bounds");
+    if (i < 0 or static_cast<size_t>(i) >= this->oper_exp.size()) MSG_ERROR("Index out of bounds");
     if (d < 0 or d >= D) MSG_ERROR("Dimension out of bounds");
     if (this->oper_exp[i][d] == nullptr) MSG_ERROR("Invalid component");
     return *this->oper_exp[i][d];
@@ -64,7 +64,7 @@ template <int D> int MWOperator<D>::getMaxBandWidth(int depth) const {
     int maxWidth = -1;
     if (depth < 0) {
         maxWidth = *std::max_element(this->band_max.begin(), this->band_max.end());
-    } else if (depth < this->band_max.size()) {
+    } else if (static_cast<size_t>(depth) < this->band_max.size()) {
         maxWidth = this->band_max[depth];
     }
     return maxWidth;
diff --git a/src/operators/MWOperator.h b/src/operators/MWOperator.h
index 2dcad2b32..833ccd341 100644
--- a/src/operators/MWOperator.h
+++ b/src/operators/MWOperator.h
@@ -49,7 +49,7 @@ template <int D> class MWOperator {
     MWOperator &operator=(const MWOperator &oper) = delete;
     virtual ~MWOperator() = default;
 
-    int size() const { return this->oper_exp.size(); }
+    size_t size() const { return this->oper_exp.size(); }
     int getMaxBandWidth(int depth = -1) const;
     const std::vector<int> &getMaxBandWidths() const { return this->band_max; }
 
diff --git a/src/operators/OperatorStatistics.cpp b/src/operators/OperatorStatistics.cpp
index f58ae2b0d..8b7be8a8b 100644
--- a/src/operators/OperatorStatistics.cpp
+++ b/src/operators/OperatorStatistics.cpp
@@ -82,6 +82,7 @@ void OperatorStatistics::flushNodeCounters() {
 
 /** Increment g-node usage counter. Needed for load balancing. */
 template <int D, typename T> void OperatorStatistics::incrementGNodeCounters(const MWNode<D, T> &gNode) {
+    (void)gNode;
     int thread = mrcpp_get_thread_num();
     this->gCount[thread]++;
 }
diff --git a/src/treebuilders/AdditionCalculator.h b/src/treebuilders/AdditionCalculator.h
index 9223f1ae6..cf207cced 100644
--- a/src/treebuilders/AdditionCalculator.h
+++ b/src/treebuilders/AdditionCalculator.h
@@ -44,7 +44,7 @@ template <int D, typename T> class AdditionCalculator final : public TreeCalcula
         node_o.zeroCoefs();
         const NodeIndex<D> &idx = node_o.getNodeIndex();
         T *coefs_o = node_o.getCoefs();
-        for (int i = 0; i < this->sum_vec.size(); i++) {
+        for (size_t i = 0; i < this->sum_vec.size(); i++) {
             T c_i = get_coef(this->sum_vec, i);
             FunctionTree<D, T> &func_i = get_func(this->sum_vec, i);
             // This generates missing nodes
diff --git a/src/treebuilders/ConvolutionCalculator.cpp b/src/treebuilders/ConvolutionCalculator.cpp
index 8cfcac1e6..753f09059 100644
--- a/src/treebuilders/ConvolutionCalculator.cpp
+++ b/src/treebuilders/ConvolutionCalculator.cpp
@@ -61,7 +61,7 @@ template <int D, typename T> ConvolutionCalculator<D, T>::~ConvolutionCalculator
     clearTimers();
     this->operStat.flushNodeCounters();
     println(10, this->operStat);
-    for (int i = 0; i < this->bandSizes.size(); i++) { delete this->bandSizes[i]; }
+    for (size_t i = 0; i < this->bandSizes.size(); i++) { delete this->bandSizes[i]; }
 }
 
 template <int D, typename T> void ConvolutionCalculator<D, T>::initTimers() {
@@ -103,7 +103,7 @@ template <int D, typename T> void ConvolutionCalculator<D, T>::printTimers() con
 /** Initialize the number of nodes formally within the bandwidth of an
  operator. The band size is used for thresholding. */
 template <int D, typename T> void ConvolutionCalculator<D, T>::initBandSizes() {
-    for (int i = 0; i < this->oper->size(); i++) {
+    for (size_t i = 0; i < this->oper->size(); i++) {
         // IMPORTANT: only 0-th dimension!
         const OperatorTree &oTree = this->oper->getComponent(i, 0);
         const BandWidth &bw = oTree.getBandWidth();
@@ -248,7 +248,7 @@ template <int D, typename T> void ConvolutionCalculator<D, T>::calcNode(MWNode<D
     os.gThreshold = gThrs;
 
     this->calc_t[mrcpp_get_thread_num()]->resume();
-    for (int n = 0; n < fBand->size(); n++) {
+    for (size_t n = 0; n < fBand->size(); n++) {
         MWNode<D, T> &fNode = *(*fBand)[n];
         NodeIndex<D> &fIdx = idx_band[n];
         os.setFNode(fNode);
@@ -277,7 +277,7 @@ template <int D, typename T> void ConvolutionCalculator<D, T>::calcNode(MWNode<D
 template <int D, typename T> void ConvolutionCalculator<D, T>::applyOperComp(OperatorState<D, T> &os) {
     double fNorm = os.fNode->getComponentNorm(os.ft);
     int o_depth = os.fNode->getScale() - this->oper->getOperatorRoot();
-    for (int i = 0; i < this->oper->size(); i++) {
+    for (size_t i = 0; i < this->oper->size(); i++) {
         // IMPORTANT: only 0-th dimension
         const OperatorTree &ot = this->oper->getComponent(i, 0);
         const BandWidth &bw = ot.getBandWidth();
diff --git a/src/treebuilders/ConvolutionCalculator.h b/src/treebuilders/ConvolutionCalculator.h
index 8ac4b5d34..2b183328b 100644
--- a/src/treebuilders/ConvolutionCalculator.h
+++ b/src/treebuilders/ConvolutionCalculator.h
@@ -59,7 +59,7 @@ template <int D, typename T> class ConvolutionCalculator final : public TreeCalc
 
     OperatorStatistics operStat;
     std::vector<Eigen::MatrixXi *> bandSizes;
-    std::function<double(const NodeIndex<D> &idx)> precFunc = [](const NodeIndex<D> &idx) { return 1.0; };
+    std::function<double(const NodeIndex<D> &idx)> precFunc = [](const NodeIndex<D> &idx) { (void)idx; return 1.0; };
 
     static const int nComp = (1 << D);
     static const int nComp2 = (1 << D) * (1 << D);
diff --git a/src/treebuilders/CopyAdaptor.cpp b/src/treebuilders/CopyAdaptor.cpp
index 8312ebb0f..ddec6227a 100644
--- a/src/treebuilders/CopyAdaptor.cpp
+++ b/src/treebuilders/CopyAdaptor.cpp
@@ -60,7 +60,7 @@ template <int D, typename T> bool CopyAdaptor<D, T>::splitNode(const MWNode<D, T
             for (int bw = -this->bandWidth[d]; bw <= this->bandWidth[d]; bw++) {
                 NodeIndex<D> bwIdx = idx.child(c);
                 bwIdx[d] += bw;
-                for (int i = 0; i < this->tree_vec.size(); i++) {
+                for (size_t i = 0; i < this->tree_vec.size(); i++) {
                     const FunctionTree<D, T> &func_i = get_func(tree_vec, i);
                     const MWNode<D, T> *node_i = func_i.findNode(bwIdx);
                     if (node_i != nullptr) return true;
diff --git a/src/treebuilders/DerivativeCalculator.cpp b/src/treebuilders/DerivativeCalculator.cpp
index f4936c01d..88e69ff93 100644
--- a/src/treebuilders/DerivativeCalculator.cpp
+++ b/src/treebuilders/DerivativeCalculator.cpp
@@ -128,7 +128,7 @@ template <int D, typename T> void DerivativeCalculator<D, T>::calcNode(MWNode<D,
 
     this->calc_t[mrcpp_get_thread_num()].resume();
 
-    for (int n = 0; n < fBand.size(); n++) {
+    for (size_t n = 0; n < fBand.size(); n++) {
         MWNode<D, T> &fNode = *fBand[n];
         NodeIndex<D> &fIdx = idx_band[n];
         os.setFNode(fNode);
diff --git a/src/treebuilders/MultiplicationCalculator.h b/src/treebuilders/MultiplicationCalculator.h
index 49fa67948..fa578ae1b 100644
--- a/src/treebuilders/MultiplicationCalculator.h
+++ b/src/treebuilders/MultiplicationCalculator.h
@@ -44,7 +44,7 @@ template <int D, typename T> class MultiplicationCalculator final : public TreeC
         const NodeIndex<D> &idx = node_o.getNodeIndex();
         T *coefs_o = node_o.getCoefs();
         for (int j = 0; j < node_o.getNCoefs(); j++) { coefs_o[j] = 1.0; }
-        for (int i = 0; i < this->prod_vec.size(); i++) {
+        for (size_t i = 0; i < this->prod_vec.size(); i++) {
             T c_i = get_coef(this->prod_vec, i);
             FunctionTree<D, T> &func_i = get_func(this->prod_vec, i);
             // This generates missing nodes
diff --git a/src/treebuilders/SplitAdaptor.h b/src/treebuilders/SplitAdaptor.h
index 7e81bbe8b..979843dd1 100644
--- a/src/treebuilders/SplitAdaptor.h
+++ b/src/treebuilders/SplitAdaptor.h
@@ -38,7 +38,7 @@ template <int D, typename T = double> class SplitAdaptor final : public TreeAdap
 private:
     bool split;
 
-    bool splitNode(const MWNode<D, T> &node) const override { return this->split; }
+    bool splitNode(const MWNode<D, T> &node) const override { (void)node; return this->split; }
 };
 
 } // namespace mrcpp
diff --git a/src/treebuilders/TimeEvolution_CrossCorrelationCalculator.cpp b/src/treebuilders/TimeEvolution_CrossCorrelationCalculator.cpp
index 844f952d9..4cf51086b 100644
--- a/src/treebuilders/TimeEvolution_CrossCorrelationCalculator.cpp
+++ b/src/treebuilders/TimeEvolution_CrossCorrelationCalculator.cpp
@@ -91,7 +91,7 @@ void TimeEvolution_CrossCorrelationCalculator::applyCcc(MWNode<2> &node) {
             for (int j = 0; j <= node.getOrder(); j++) {
                 // std::min(M, N)  could be used for breaking the following loop
                 // this->cross_correlation->Matrix.size() should be big enough a priori
-                for (int k = 0; 2 * k + p + j < J_power_inetgarls[l_b].size(); k++) {
+                for (size_t k = 0; 2 * k + p + j < J_power_inetgarls[l_b].size(); k++) {
                     double J;
                     if (this->imaginary)
                         J = J_power_inetgarls[l_b][2 * k + p + j].imag();
diff --git a/src/treebuilders/TreeAdaptor.h b/src/treebuilders/TreeAdaptor.h
index 80cecb09e..9e81e6c13 100644
--- a/src/treebuilders/TreeAdaptor.h
+++ b/src/treebuilders/TreeAdaptor.h
@@ -39,7 +39,7 @@ template <int D, typename T> class TreeAdaptor {
     void setMaxScale(int ms) { this->maxScale = ms; }
 
     void splitNodeVector(MWNodeVector<D, T> &out, MWNodeVector<D, T> &inp) const {
-        for (int n = 0; n < inp.size(); n++) {
+        for (size_t n = 0; n < inp.size(); n++) {
             MWNode<D, T> &node = *inp[n];
             // Can be BranchNode in operator application
             if (node.isBranchNode()) continue;
diff --git a/src/treebuilders/TreeBuilder.cpp b/src/treebuilders/TreeBuilder.cpp
index ba0e5d973..f034a886e 100644
--- a/src/treebuilders/TreeBuilder.cpp
+++ b/src/treebuilders/TreeBuilder.cpp
@@ -111,7 +111,7 @@ template <int D, typename T> int TreeBuilder<D, T>::split(MWTree<D, T> &tree, Tr
     MWNodeVector<D, T> *workVec = tree.copyEndNodeTable();
     adaptor.splitNodeVector(newVec, *workVec);
     if (passCoefs) {
-        for (int i = 0; i < workVec->size(); i++) {
+        for (size_t i = 0; i < workVec->size(); i++) {
             MWNode<D, T> &node = *(*workVec)[i];
             if (node.isBranchNode()) { node.giveChildrenCoefs(true); }
         }
@@ -149,7 +149,7 @@ template <int D, typename T> void TreeBuilder<D, T>::calc(MWTree<D, T> &tree, Tr
 
 template <int D, typename T> double TreeBuilder<D, T>::calcScalingNorm(const MWNodeVector<D, T> &vec) const {
     double sNorm = 0.0;
-    for (int i = 0; i < vec.size(); i++) {
+    for (size_t i = 0; i < vec.size(); i++) {
         const MWNode<D, T> &node = *vec[i];
         if (node.getDepth() >= 0) sNorm += node.getScalingNorm();
     }
@@ -158,7 +158,7 @@ template <int D, typename T> double TreeBuilder<D, T>::calcScalingNorm(const MWN
 
 template <int D, typename T> double TreeBuilder<D, T>::calcWaveletNorm(const MWNodeVector<D, T> &vec) const {
     double wNorm = 0.0;
-    for (int i = 0; i < vec.size(); i++) {
+    for (size_t i = 0; i < vec.size(); i++) {
         const MWNode<D, T> &node = *vec[i];
         if (node.getDepth() >= 0) wNorm += node.getWaveletNorm();
     }
diff --git a/src/treebuilders/WaveletAdaptor.h b/src/treebuilders/WaveletAdaptor.h
index 829039bf4..78542e7ee 100644
--- a/src/treebuilders/WaveletAdaptor.h
+++ b/src/treebuilders/WaveletAdaptor.h
@@ -46,7 +46,7 @@ template <int D, typename T> class WaveletAdaptor : public TreeAdaptor<D, T> {
     bool absPrec;
     double prec;
     double splitFac;
-    std::function<double(const NodeIndex<D> &idx)> precFunc = [](const NodeIndex<D> &idx) { return 1.0; };
+    std::function<double(const NodeIndex<D> &idx)> precFunc = [](const NodeIndex<D> &idx) { (void)idx; return 1.0; };
 
     bool splitNode(const MWNode<D, T> &node) const override {
         auto precFac = this->precFunc(node.getNodeIndex()); // returns 1.0 by default
diff --git a/src/treebuilders/add.cpp b/src/treebuilders/add.cpp
index 4ee28cff6..3d2709633 100644
--- a/src/treebuilders/add.cpp
+++ b/src/treebuilders/add.cpp
@@ -91,7 +91,7 @@ template <int D, typename T> void add(double prec, FunctionTree<D, T> &out, T a,
  *
  */
 template <int D, typename T> void add(double prec, FunctionTree<D, T> &out, FunctionTreeVector<D, T> &inp, int maxIter, bool absPrec, bool conjugate) {
-    for (auto i = 0; i < inp.size(); i++)
+    for (size_t i = 0; i < inp.size(); i++)
         if (out.getMRA() != get_func(inp, i).getMRA()) MSG_ABORT("Incompatible MRA");
 
     int maxScale = out.getMRA().getMaxScale();
@@ -108,7 +108,7 @@ template <int D, typename T> void add(double prec, FunctionTree<D, T> &out, Func
     trans_t.stop();
 
     Timer clean_t;
-    for (int i = 0; i < inp.size(); i++) {
+    for (size_t i = 0; i < inp.size(); i++) {
         FunctionTree<D, T> &tree = get_func(inp, i);
         tree.deleteGenerated();
     }
diff --git a/src/treebuilders/apply.cpp b/src/treebuilders/apply.cpp
index edc8db17e..e83065a47 100644
--- a/src/treebuilders/apply.cpp
+++ b/src/treebuilders/apply.cpp
@@ -218,10 +218,10 @@ template <int D, typename T> void apply(double prec, FunctionTree<D, T> &out, Co
 
     // The local precision will be scaled by the maxNorm of the
     // corresponding node(s) in the precTrees vector.
-    for (int i = 0; i < precTrees.size(); i++) get_func(precTrees, i).makeMaxSquareNorms();
+    for (size_t i = 0; i < precTrees.size(); i++) get_func(precTrees, i).makeMaxSquareNorms();
     auto precFunc = [&precTrees](const NodeIndex<D> &idx) -> double {
         auto maxNorm = (precTrees.size()) ? 0.0 : 1.0;
-        for (int i = 0; i < precTrees.size(); i++) {
+        for (size_t i = 0; i < precTrees.size(); i++) {
             auto &pNode = get_func(precTrees, i).getNode(idx);
             maxNorm = std::max(maxNorm, std::sqrt(pNode.getMaxSquareNorm()));
         }
@@ -513,7 +513,7 @@ std::vector<CompFunction<3> *> gradient(DerivativeOperator<3> &oper, CompFunctio
  */
 template <int D, typename T> void divergence(FunctionTree<D, T> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> &inp) {
     if (inp.size() != D) MSG_ABORT("Dimension mismatch");
-    for (auto i = 0; i < inp.size(); i++)
+    for (size_t i = 0; i < inp.size(); i++)
         if (out.getMRA() != get_func(inp, i).getMRA()) MSG_ABORT("Incompatible MRA");
 
     FunctionTreeVector<D, T> tmp_vec;
@@ -530,6 +530,10 @@ template <int D, typename T> void divergence(FunctionTree<D, T> &out, Derivative
 }
 
 template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, FunctionTreeVector<D, T> *inp, const ComplexDouble (*metric)[4]) {
+    (void)out;
+    (void)oper;
+    (void)inp;
+    (void)metric;
     MSG_ABORT("not implemented");
 }
 
@@ -539,6 +543,10 @@ template <int D, typename T> void divergence(FunctionTree<D, T> &out, Derivative
     divergence(out, oper, inp_vec);
 }
 template <int D, typename T> void divergence(CompFunction<D> &out, DerivativeOperator<D> &oper, std::vector<FunctionTree<D, T> *> *inp, const ComplexDouble (*metric)[4]) {
+    (void)out;
+    (void)oper;
+    (void)inp;
+    (void)metric;
     MSG_ABORT("not implemented");
 }
 
diff --git a/src/treebuilders/grid.cpp b/src/treebuilders/grid.cpp
index eaaf3696b..5341bd3c1 100644
--- a/src/treebuilders/grid.cpp
+++ b/src/treebuilders/grid.cpp
@@ -108,12 +108,12 @@ template <int D> void build_grid(FunctionTree<D> &out, const GaussExp<D> &inp, i
         auto maxScale = out.getMRA().getMaxScale();
         TreeBuilder<D> builder;
         DefaultCalculator<D> calculator;
-        for (auto i = 0; i < inp.size(); i++) {
+        for (size_t i = 0; i < inp.size(); i++) {
             AnalyticAdaptor<D> adaptor(inp.getFunc(i), maxScale);
             builder.build(out, calculator, adaptor, maxIter);
         }
     } else {
-        for (auto i = 0; i < inp.size(); i++) {
+        for (size_t i = 0; i < inp.size(); i++) {
             auto *gauss = inp.getFunc(i).copy();
             build_grid(out, *gauss, maxIter);
             delete gauss;
@@ -171,7 +171,7 @@ template <int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTr
  *
  */
 template <int D, typename T> void build_grid(FunctionTree<D, T> &out, FunctionTreeVector<D, T> &inp, int maxIter) {
-    for (auto i = 0; i < inp.size(); i++)
+    for (size_t i = 0; i < inp.size(); i++)
         if (out.getMRA() != get_func(inp, i).getMRA()) MSG_ABORT("Incompatible MRA");
 
     auto maxScale = out.getMRA().getMaxScale();
diff --git a/src/treebuilders/multiply.cpp b/src/treebuilders/multiply.cpp
index 46587ffaa..0a8ef63f7 100644
--- a/src/treebuilders/multiply.cpp
+++ b/src/treebuilders/multiply.cpp
@@ -102,7 +102,7 @@ void multiply(double prec, FunctionTree<D, T> &out, T c, FunctionTree<D, T> &inp
  *
  */
 template <int D, typename T> void multiply(double prec, FunctionTree<D, T> &out, FunctionTreeVector<D, T> &inp, int maxIter, bool absPrec, bool useMaxNorms, bool conjugate) {
-    for (auto i = 0; i < inp.size(); i++)
+    for (size_t i = 0; i < inp.size(); i++)
         if (out.getMRA() != get_func(inp, i).getMRA()) MSG_ABORT("Incompatible MRA");
 
     int maxScale = out.getMRA().getMaxScale();
@@ -110,7 +110,7 @@ template <int D, typename T> void multiply(double prec, FunctionTree<D, T> &out,
     MultiplicationCalculator<D, T> calculator(inp, conjugate);
 
     if (useMaxNorms) {
-        for (int i = 0; i < inp.size(); i++) get_func(inp, i).makeMaxSquareNorms();
+        for (size_t i = 0; i < inp.size(); i++) get_func(inp, i).makeMaxSquareNorms();
         MultiplicationAdaptor<D, T> adaptor(prec, maxScale, inp);
         builder.build(out, calculator, adaptor, maxIter);
     } else {
@@ -124,7 +124,7 @@ template <int D, typename T> void multiply(double prec, FunctionTree<D, T> &out,
     trans_t.stop();
 
     Timer clean_t;
-    for (int i = 0; i < inp.size(); i++) {
+    for (size_t i = 0; i < inp.size(); i++) {
         FunctionTree<D, T> &tree = get_func(inp, i);
         tree.deleteGenerated();
     }
@@ -254,7 +254,7 @@ template <int D, typename T> void dot(double prec, FunctionTree<D, T> &out, Func
     if (inp_a.size() != inp_b.size()) MSG_ABORT("Input length mismatch");
 
     FunctionTreeVector<D, T> tmp_vec;
-    for (int d = 0; d < inp_a.size(); d++) {
+    for (size_t d = 0; d < inp_a.size(); d++) {
         T coef_a = get_coef(inp_a, d);
         T coef_b = get_coef(inp_b, d);
         FunctionTree<D, T> &tree_a = get_func(inp_a, d);
diff --git a/src/trees/BoundingBox.cpp b/src/trees/BoundingBox.cpp
index ff86abf2e..65cc98883 100644
--- a/src/trees/BoundingBox.cpp
+++ b/src/trees/BoundingBox.cpp
@@ -455,7 +455,7 @@ template <> int BoundingBox<1>::getBoxIndex(NodeIndex<1> nIdx) const {
     if (relScale < 0) return -1;
 
     int bIdx = (l >> relScale) - cl;
-    if (bIdx < 0 or bIdx >= this->size()) {
+    if (bIdx < 0 or bIdx >= static_cast<int>(this->size())) {
         return -1;
     } else {
         return bIdx;
diff --git a/src/trees/BoundingBox.h b/src/trees/BoundingBox.h
index 28ad9c052..6043b711f 100644
--- a/src/trees/BoundingBox.h
+++ b/src/trees/BoundingBox.h
@@ -67,7 +67,7 @@ template <int D> class BoundingBox {
     int getBoxIndex(Coord<D> r) const;
     int getBoxIndex(NodeIndex<D> nIdx) const;
 
-    int size() const { return this->totBoxes; }
+    size_t size() const { return this->totBoxes; }
     int size(int d) const { return this->nBoxes[d]; }
     int getScale() const { return this->cornerIndex.getScale(); }
     double getScalingFactor(int d) const { return this->scalingFactor[d]; }
diff --git a/src/trees/FunctionTree.cpp b/src/trees/FunctionTree.cpp
index 7d5b764a9..0b2c2a1d0 100644
--- a/src/trees/FunctionTree.cpp
+++ b/src/trees/FunctionTree.cpp
@@ -239,7 +239,7 @@ template <int D, typename T> void FunctionTree<D, T>::loadTreeTXT(const std::str
     in.close();
     // transform all nodes from quadrature point values to scaling coefficients
     for (int n = nmax; n > -1; n--) {
-        for (int i = 0; i < NodeTable[n].size(); i++) {
+        for (size_t i = 0; i < NodeTable[n].size(); i++) {
             MWNode<D, T> *node = NodeTable[n][i];
             node->cvTransform(Backward);
             node->calcNorms();
@@ -249,7 +249,7 @@ template <int D, typename T> void FunctionTree<D, T>::loadTreeTXT(const std::str
 
     // Transform into scaling and wavelets, starting by leaf nodes and copying scaling into parents
     for (int n = nmax; n > -1; n--) {
-        for (int i = 0; i < NodeTable[n].size(); i++) {
+        for (size_t i = 0; i < NodeTable[n].size(); i++) {
             MWNode<D, T> *node = NodeTable[n][i];
             if (mp[node->getSerialIx()] == nChildren) {
                 // node complete: transform into scaling and wavelets
@@ -438,7 +438,7 @@ template <int D, typename T> void FunctionTree<D, T>::loadTree(const std::string
 template <int D, typename T> T FunctionTree<D, T>::integrate() const {
 
     T result = 0.0;
-    for (int i = 0; i < this->rootBox.size(); i++) {
+    for (size_t i = 0; i < this->rootBox.size(); i++) {
         const FunctionNode<D, T> &fNode = getRootFuncNode(i);
         result += fNode.integrate();
     }
@@ -455,7 +455,7 @@ template <int D, typename T> T FunctionTree<D, T>::integrate() const {
 
 template <int D, typename T> T FunctionTree<D, T>::integrateSide(int dim, bool positiveSide) const {
     T result = 0.0;
-    for (int i = 0; i < this->rootBox.size(); i++) {
+    for (size_t i = 0; i < this->rootBox.size(); i++) {
         const FunctionNode<D, T> &fNode = getRootFuncNode(i);
         if (fNode.getLowerBounds()[dim] == this->rootBox.getLowerBound(dim) && not positiveSide)
             result += fNode.integrate();
@@ -477,7 +477,7 @@ template <int D, typename T> T FunctionTree<D, T>::integrateSide(int dim, bool p
 template <> double FunctionTree<3, double>::integrateEndNodes(RepresentableFunction_M &f) {
     // traverse tree, and treat end nodes only
     std::vector<FunctionNode<3> *> stack; // node from this
-    for (int i = 0; i < this->getRootBox().size(); i++) stack.push_back(&(this->getRootFuncNode(i)));
+    for (size_t i = 0; i < this->getRootBox().size(); i++) stack.push_back(&(this->getRootFuncNode(i)));
     double result = 0.0;
     while (stack.size() > 0) {
         FunctionNode<3> *Node = stack.back();
@@ -888,7 +888,7 @@ template <int D, typename T> std::ostream &FunctionTree<D, T>::print(std::ostrea
  * to the dimension; in practice, it is set to `s=1`.
  */
 template <int D, typename T> int FunctionTree<D, T>::crop(double prec, double splitFac, bool absPrec) {
-    for (int i = 0; i < this->rootBox.size(); i++) {
+    for (size_t i = 0; i < this->rootBox.size(); i++) {
         MWNode<D, T> &root = this->getRootMWNode(i);
         root.crop(prec, splitFac, absPrec);
     }
@@ -916,11 +916,11 @@ void FunctionTree<D, T>::makeCoeffVector(std::vector<T *> &coefs,
     max_index = 0;
     std::vector<MWNode<D, double> *> refstack; // nodes from refTree
     std::vector<MWNode<D, T> *> thisstack;     // nodes from this Tree
-    for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
         refstack.push_back(refTree.getRootBox().getNodes()[rIdx]);
         thisstack.push_back(this->getRootBox().getNodes()[rIdx]);
     }
-    int stack_p = 0;
+    size_t stack_p = 0;
     while (thisstack.size() > stack_p) {
         // refNode and thisNode are the same node in space, but on different trees
         MWNode<D, T> *thisNode = thisstack[stack_p];
@@ -961,7 +961,7 @@ template <int D, typename T> void FunctionTree<D, T>::makeTreefromCoeff(MWTree<D
     int sizecoefW = ((1 << this->getDim()) - 1) * this->getKp1_d();
     this->squareNorm = 0.0;
     this->clearEndNodeTable();
-    for (int rIdx = 0; rIdx < refTree.getRootBox().size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < refTree.getRootBox().size(); rIdx++) {
         MWNode<D, double> *refNode = refTree.getRootBox().getNodes()[rIdx];
         stack.push_back(refNode);
         int ix = ix2coef[refNode->getSerialIx()];
@@ -986,7 +986,7 @@ template <int D, typename T> void FunctionTree<D, T>::makeTreefromCoeff(MWTree<D
             }
         } else {
             // only wavelets are defined in coefVec. Scaling part set below, when creating children
-            if (ix < coefpVec.size() and ix >= 0) {
+            if (ix < static_cast<int>(coefpVec.size()) and ix >= 0) {
                 for (int k = 0; k < size; k++) node->getCoefs()[k + this->getKp1_d()] = coefpVec[ix][k];
             } else {
                 // we do not have W coefficients. Set them to zero
@@ -1037,7 +1037,7 @@ template <int D, typename T> void FunctionTree<D, T>::appendTreeNoCoeff(MWTree<D
     std::vector<MWNode<D, double> *> instack; // node from inTree
     std::vector<MWNode<D, T> *> thisstack;    // node from this Tree
     this->clearEndNodeTable();
-    for (int rIdx = 0; rIdx < inTree.getRootBox().size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < inTree.getRootBox().size(); rIdx++) {
         instack.push_back(inTree.getRootBox().getNodes()[rIdx]);
         thisstack.push_back(this->getRootBox().getNodes()[rIdx]);
     }
@@ -1076,7 +1076,7 @@ template <int D, typename T> void FunctionTree<D, T>::appendTreeNoCoeff(MWTree<D
     std::vector<MWNode<D, ComplexDouble> *> instack; // node from inTree
     std::vector<MWNode<D, T> *> thisstack;           // node from this Tree
     this->clearEndNodeTable();
-    for (int rIdx = 0; rIdx < inTree.getRootBox().size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < inTree.getRootBox().size(); rIdx++) {
         instack.push_back(inTree.getRootBox().getNodes()[rIdx]);
         thisstack.push_back(this->getRootBox().getNodes()[rIdx]);
     }
@@ -1115,17 +1115,17 @@ template <int D, typename T> void FunctionTree<D, T>::deleteGenerated() {
 }
 
 template <int D, typename T> void FunctionTree<D, T>::deleteGeneratedParents() {
-    for (int n = 0; n < this->getRootBox().size(); n++) this->getRootMWNode(n).deleteParent();
+    for (size_t n = 0; n < this->getRootBox().size(); n++) this->getRootMWNode(n).deleteParent();
 }
 
 template <> int FunctionTree<3, double>::saveNodesAndRmCoeff() {
     if (this->isLocal) MSG_INFO("Tree is already in local representation");
     NodesCoeff = new BankAccount; // NB: must be a collective call!
-    int stack_p = 0;
+    size_t stack_p = 0;
     if (mpi::wrk_rank == 0) {
         int sizecoeff = (1 << 3) * this->getKp1_d();
         std::vector<MWNode<3, double> *> stack; // nodes from this Tree
-        for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
+        for (size_t rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
         while (stack.size() > stack_p) {
             MWNode<3, double> *Node = stack[stack_p++];
             NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs());
@@ -1142,12 +1142,12 @@ template <> int FunctionTree<3, double>::saveNodesAndRmCoeff() {
 template <> int FunctionTree<3, ComplexDouble>::saveNodesAndRmCoeff() {
     if (this->isLocal) MSG_INFO("Tree is already in local representation");
     NodesCoeff = new BankAccount; // NB: must be a collective call!
-    int stack_p = 0;
+    size_t stack_p = 0;
     if (mpi::wrk_rank == 0) {
         int sizecoeff = (1 << 3) * this->getKp1_d();
         sizecoeff *= 2;                                // double->ComplexDouble. Saved as twice as many doubles
         std::vector<MWNode<3, ComplexDouble> *> stack; // nodes from this Tree
-        for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
+        for (size_t rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) { stack.push_back(this->getRootBox().getNodes()[rIdx]); }
         while (stack.size() > stack_p) {
             MWNode<3, ComplexDouble> *Node = stack[stack_p++];
             NodesCoeff->put_data(Node->getNodeIndex(), sizecoeff, Node->getCoefs());
@@ -1229,7 +1229,7 @@ FunctionTree<D, ComplexDouble>* FunctionTree<D, T>::CopyTreeToComplex() {
     std::vector<MWNode<D, double> *> instack;         // node from this
     std::vector<MWNode<D, ComplexDouble> *> outstack; // node from outTree
     outTree->clearEndNodeTable();
-    for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
         instack.push_back(this->getRootBox().getNodes()[rIdx]);
         outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
     }
@@ -1272,7 +1272,7 @@ FunctionTree<D, double>* FunctionTree<D, T>::CopyTreeToReal() {
     std::vector<MWNode<D, double> *> instack;  // node from this
     std::vector<MWNode<D, double> *> outstack; // node from outTree
     outTree->clearEndNodeTable();
-    for (int rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < this->getRootBox().size(); rIdx++) {
         instack.push_back(this->getRootBox().getNodes()[rIdx]);
         outstack.push_back(outTree->getRootBox().getNodes()[rIdx]);
     }
diff --git a/src/trees/MWNode.cpp b/src/trees/MWNode.cpp
index 9f40bed23..e85d0d3ce 100644
--- a/src/trees/MWNode.cpp
+++ b/src/trees/MWNode.cpp
@@ -368,6 +368,7 @@ template <int D, typename T> void MWNode<D, T>::giveChildCoefs(int cIdx, bool ov
  * \warning This routine is only used in connection with Periodic Boundary Conditions
  */
 template <int D, typename T> void MWNode<D, T>::giveParentCoefs(bool overwrite) {
+    (void)overwrite;
     MWNode<D, T> node = *this;
     MWNode<D, T> &parent = getMWParent();
     int kp1_d = this->getKp1_d();
@@ -692,6 +693,7 @@ template <int D, typename T> bool MWNode<D, T>::crop(double prec, double splitFa
 }
 
 template <int D, typename T> void MWNode<D, T>::createChildren(bool coefs) {
+    (void)coefs;
     NOT_REACHED_ABORT;
 }
 
@@ -1187,6 +1189,7 @@ template <int D, typename T> bool MWNode<D, T>::hasCoord(const Coord<D> &r) cons
 /** Testing if nodes are compatible wrt NodeIndex and Tree (order, rootScale,
  * relPrec, etc). */
 template <int D, typename T> bool MWNode<D, T>::isCompatible(const MWNode<D, T> &node) {
+    (void)node;
     NOT_IMPLEMENTED_ABORT;
     //    if (nodeIndex != node.nodeIndex) {
     //        println(0, "nodeIndex mismatch" << std::endl);
@@ -1216,6 +1219,7 @@ template <int D, typename T> bool MWNode<D, T>::isAncestor(const NodeIndex<D> &i
 }
 
 template <int D, typename T> bool MWNode<D, T>::isDecendant(const NodeIndex<D> &idx) const {
+    (void)idx;
     NOT_IMPLEMENTED_ABORT;
 }
 
diff --git a/src/trees/MWTree.cpp b/src/trees/MWTree.cpp
index 6a646d33f..da1cc9e82 100644
--- a/src/trees/MWTree.cpp
+++ b/src/trees/MWTree.cpp
@@ -74,7 +74,7 @@ template <int D, typename T> MWTree<D, T>::~MWTree() {
  * when the object is deleted.
  */
 template <int D, typename T> void MWTree<D, T>::deleteRootNodes() {
-    for (int i = 0; i < this->rootBox.size(); i++) {
+    for (size_t i = 0; i < this->rootBox.size(); i++) {
         MWNode<D, T> &root = this->getRootMWNode(i);
         root.deleteChildren();
         root.dealloc();
@@ -91,7 +91,7 @@ template <int D, typename T> void MWTree<D, T>::deleteRootNodes() {
  * immediately available to the new function.
  */
 template <int D, typename T> void MWTree<D, T>::clear() {
-    for (int i = 0; i < this->rootBox.size(); i++) {
+    for (size_t i = 0; i < this->rootBox.size(); i++) {
         MWNode<D, T> &root = this->getRootMWNode(i);
         root.deleteChildren();
         root.clearHasCoefs();
@@ -196,7 +196,7 @@ template <int D, typename T> void MWTree<D, T>::mwTransformDown(bool overwrite)
     tree_utils::make_node_table(*this, nodeTable);
 #pragma omp parallel shared(nodeTable) num_threads(mrcpp_get_num_threads())
     {
-        for (int n = 0; n < nodeTable.size(); n++) {
+        for (size_t n = 0; n < nodeTable.size(); n++) {
             int n_nodes = nodeTable[n].size();
 #pragma omp for schedule(guided)
             for (int i = 0; i < n_nodes; i++) {
@@ -284,9 +284,9 @@ template <int D, typename T> void MWTree<D, T>::decrementNodeCount(int scale) {
 template <int D, typename T> int MWTree<D, T>::getNNodesAtDepth(int depth) const {
     int N = 0;
     if (depth < 0) {
-        if (this->nodesAtNegativeDepth.size() >= -depth) N = this->nodesAtNegativeDepth[-depth];
+        if (this->nodesAtNegativeDepth.size() >= static_cast<size_t>(-depth)) N = this->nodesAtNegativeDepth[-depth];
     } else {
-        if (this->nodesAtDepth.size() > depth) N = this->nodesAtDepth[depth];
+        if (this->nodesAtDepth.size() > static_cast<size_t>(depth)) N = this->nodesAtDepth[depth];
     }
     return N;
 }
@@ -462,10 +462,12 @@ template <int D, typename T> void MWTree<D, T>::resetEndNodeTable() {
 }
 
 template <int D, typename T> int MWTree<D, T>::countBranchNodes(int depth) {
+    (void)depth;
     NOT_IMPLEMENTED_ABORT;
 }
 
 template <int D, typename T> int MWTree<D, T>::countLeafNodes(int depth) {
+    (void)depth;
     NOT_IMPLEMENTED_ABORT;
     //    int nNodes = 0;
     //    TreeIterator<D, T> it(*this);
@@ -482,6 +484,7 @@ template <int D, typename T> int MWTree<D, T>::countLeafNodes(int depth) {
 
 /* Traverse tree and count nodes belonging to this rank. */
 template <int D, typename T> int MWTree<D, T>::countNodes(int depth) {
+    (void)depth;
     NOT_IMPLEMENTED_ABORT;
     //    TreeIterator<D, T> it(*this);
     //    int count = 0;
@@ -499,6 +502,7 @@ template <int D, typename T> int MWTree<D, T>::countNodes(int depth) {
 
 /* Traverse tree and count nodes with allocated coefficients. */
 template <int D, typename T> int MWTree<D, T>::countAllocNodes(int depth) {
+    (void)depth;
     NOT_IMPLEMENTED_ABORT;
     //    TreeIterator<D, T> it(*this);
     //    int count = 0;
@@ -523,8 +527,8 @@ template <int D, typename T> std::ostream &MWTree<D, T>::print(std::ostream &o)
     o << "  nodes: " << this->getNNodes() << std::endl;
     o << "  endNodes: " << this->endNodeTable.size() << std::endl;
     o << "  nodes per scale: " << std::endl;
-    for (int i = this->nodesAtNegativeDepth.size() - 1; i >= 0; i--) { o << "    scale=" << -(i + this->getRootScale() + 1) << "  nodes=" << this->nodesAtNegativeDepth[i] << std::endl; }
-    for (int i = 0; i < this->nodesAtDepth.size(); i++) { o << "    scale=" << i + this->getRootScale() << "  nodes=" << this->nodesAtDepth[i] << std::endl; }
+    for (int i = static_cast<int>(this->nodesAtNegativeDepth.size()) - 1; i >= 0; i--) { o << "    scale=" << -(i + this->getRootScale() + 1) << "  nodes=" << this->nodesAtNegativeDepth[i] << std::endl; }
+    for (size_t i = 0; i < this->nodesAtDepth.size(); i++) { o << "    scale=" << i + this->getRootScale() << "  nodes=" << this->nodesAtDepth[i] << std::endl; }
     return o;
 }
 
@@ -536,7 +540,7 @@ template <int D, typename T> std::ostream &MWTree<D, T>::print(std::ostream &o)
 template <int D, typename T> void MWTree<D, T>::makeMaxSquareNorms() {
     NodeBox<D, T> &rBox = this->getRootBox();
     MWNode<D, T> **roots = rBox.getNodes();
-    for (int rIdx = 0; rIdx < rBox.size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < rBox.size(); rIdx++) {
         // recursively set value of children and descendants
         roots[rIdx]->setMaxSquareNorm();
     }
diff --git a/src/trees/NodeAllocator.cpp b/src/trees/NodeAllocator.cpp
index cd3cf9fcc..9ae9e1962 100644
--- a/src/trees/NodeAllocator.cpp
+++ b/src/trees/NodeAllocator.cpp
@@ -73,6 +73,10 @@ NodeAllocator<2>::NodeAllocator(OperatorTree *tree, SharedMemory<double> *mem, i
 }
 
 template <int D, typename T> NodeAllocator<D, T>::NodeAllocator(OperatorTree *tree, SharedMemory<T> *mem, int coefsPerNode, int nodesPerChunk) {
+    (void)tree;
+    (void)mem;
+    (void)coefsPerNode;
+    (void)nodesPerChunk;
     NOT_REACHED_ABORT;
 }
 
@@ -99,14 +103,14 @@ template <int D, typename T> T *NodeAllocator<D, T>::getCoef_p(int sIdx) {
 }
 
 template <int D, typename T> MWNode<D, T> *NodeAllocator<D, T>::getNodeNoLock(int sIdx) {
-    if (sIdx < 0 or sIdx >= this->stackStatus.size()) return nullptr;
+    if (sIdx < 0 or static_cast<size_t>(sIdx) >= this->stackStatus.size()) return nullptr;
     int chunk = sIdx / this->maxNodesPerChunk; // which chunk
     int cIdx = sIdx % this->maxNodesPerChunk;  // position in chunk
     return this->nodeChunks[chunk] + cIdx;
 }
 
 template <int D, typename T> T *NodeAllocator<D, T>::getCoefNoLock(int sIdx) {
-    if (sIdx < 0 or sIdx >= this->stackStatus.size()) return nullptr;
+    if (sIdx < 0 or static_cast<size_t>(sIdx) >= this->stackStatus.size()) return nullptr;
     int chunk = sIdx / this->maxNodesPerChunk; // which chunk
     int idx = sIdx % this->maxNodesPerChunk;   // position in chunk
     return this->coefChunks[chunk] + idx * this->coefsPerNode;
@@ -121,7 +125,7 @@ template <int D, typename T> int NodeAllocator<D, T>::alloc(int nNodes, bool coe
     if (chunkOverflow) this->topStack = this->maxNodesPerChunk * ((this->topStack + nNodes - 1) / this->maxNodesPerChunk);
 
     // append chunk if necessary
-    int chunk = this->topStack / this->maxNodesPerChunk;
+    size_t chunk = this->topStack / this->maxNodesPerChunk;
     bool needNewChunk = (chunk >= this->nodeChunks.size());
     if (needNewChunk) appendChunk(coefs);
 
@@ -146,7 +150,7 @@ template <int D, typename T> int NodeAllocator<D, T>::alloc(int nNodes, bool coe
 
 template <int D, typename T> void NodeAllocator<D, T>::dealloc(int sIdx) {
     MRCPP_SET_OMP_LOCK();
-    if (sIdx < 0 or sIdx >= this->stackStatus.size()) MSG_ABORT("Invalid serial index: " << sIdx);
+    if (sIdx < 0 or static_cast<size_t>(sIdx) >= this->stackStatus.size()) MSG_ABORT("Invalid serial index: " << sIdx);
     auto *node_p = getNodeNoLock(sIdx);
     node_p->~MWNode();
     this->stackStatus[sIdx] = 0;      // mark as available
@@ -219,7 +223,7 @@ template <int D, typename T> void NodeAllocator<D, T>::appendChunk(bool coefs) {
 template <int D, typename T> int NodeAllocator<D, T>::compress() {
     MRCPP_SET_OMP_LOCK();
     int nNodes = (1 << D);
-    if (this->maxNodesPerChunk * this->nodeChunks.size() <= getTree().getNNodes() + this->maxNodesPerChunk + nNodes - 1) {
+    if (this->maxNodesPerChunk * this->nodeChunks.size() <= static_cast<size_t>(getTree().getNNodes() + this->maxNodesPerChunk + nNodes - 1)) {
         MRCPP_UNSET_OMP_LOCK();
         return 0; // nothing to compress
     }
@@ -370,7 +374,7 @@ template <int D, typename T> void NodeAllocator<D, T>::reassemble() {
     MWNode<D, T> **roots = rootbox.getNodes();
 
     std::stack<MWNode<D, T> *> stack;
-    for (int rIdx = 0; rIdx < rootbox.size(); rIdx++) {
+    for (size_t rIdx = 0; rIdx < rootbox.size(); rIdx++) {
         auto *root_p = getNodeNoLock(rIdx);
         assert(root_p != nullptr);
         stack.push(root_p);
diff --git a/src/trees/NodeBox.cpp b/src/trees/NodeBox.cpp
index bf747c4fc..0e3acfa68 100644
--- a/src/trees/NodeBox.cpp
+++ b/src/trees/NodeBox.cpp
@@ -74,7 +74,7 @@ template <int D, typename T> NodeBox<D, T>::~NodeBox() {
 
 template <int D, typename T> void NodeBox<D, T>::deleteNodes() {
     if (this->nodes == nullptr) { return; }
-    for (int n = 0; n < this->size(); n++) { clearNode(n); }
+    for (size_t n = 0; n < this->size(); n++) { clearNode(n); }
     delete[] this->nodes;
     this->nodes = nullptr;
 }
diff --git a/src/utils/Bank.cpp b/src/utils/Bank.cpp
index 94152c06e..ef26d6094 100644
--- a/src/utils/Bank.cpp
+++ b/src/utils/Bank.cpp
@@ -470,6 +470,9 @@ int Bank::clearAccount(int account, int iclient, MPI_Comm comm) {
     closeAccount(account);
     return openAccount(iclient, comm);
 #else
+    (void)account;
+    (void)iclient;
+    (void)comm;
     return 1;
 #endif
 }
@@ -518,6 +521,8 @@ void Bank::remove_account(int account) {
     for (double *c_p : mem[account]->chunk_p) delete[] c_p;
     mem.erase(account);
     currentsize.erase(account);
+#else
+    (void)account;
 #endif
 }
 
@@ -544,6 +549,9 @@ int Bank::openAccount(int iclient, MPI_Comm comm) {
     } else {
         MPI_Bcast(account_id, 1, MPI_INT, 0, comm);
     }
+#else
+    (void)iclient;
+    (void)comm;
 #endif
     return account_id[0];
 }
@@ -580,6 +588,10 @@ int Bank::openTaskManager(int ntasks, int iclient, MPI_Comm comm) {
     } else {
         MPI_Bcast(&account_id, 1, MPI_INT, 0, comm);
     }
+#else
+    (void)ntasks;
+    (void)iclient;
+    (void)comm;
 #endif
     return account_id;
 }
@@ -592,6 +604,8 @@ void Bank::closeAccount(int account_id) {
     messages[0] = CLOSE_ACCOUNT;
     messages[1] = account_id;
     for (int i = 0; i < bank_size; i++) MPI_Send(messages, 2, MPI_INT, bankmaster[i], 0, comm_bank);
+#else
+    (void)account_id;
 #endif
 }
 
@@ -603,6 +617,8 @@ void Bank::closeTaskManager(int account_id) {
     messages[0] = CLOSE_ACCOUNT;
     messages[1] = account_id;
     MPI_Send(messages, 2, MPI_INT, task_bank, 0, comm_bank);
+#else
+    (void)account_id;
 #endif
 }
 
@@ -667,6 +683,10 @@ int BankAccount::get_func(int id, CompFunction<3> &func, int wait) {
         MPI_Send(messages, 3, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
         recv_function(func, bankmaster[id % bank_size], 1, comm_bank);
     }
+#else
+    (void)id;
+    (void)func;
+    (void)wait;
 #endif
     return 1;
 }
@@ -689,6 +709,9 @@ int BankAccount::get_func_del(int id, CompFunction<3> &orb) {
     } else {
         return 0;
     }
+#else
+    (void)id;
+    (void)orb;
 #endif
     return 1;
 }
@@ -703,6 +726,9 @@ int BankAccount::put_func(int id, CompFunction<3> &func) {
     messages[2] = id;
     MPI_Send(messages, 3, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     send_function(func, bankmaster[id % bank_size], 1, comm_bank);
+#else
+    (void)id;
+    (void)func;
 #endif
     return 1;
 }
@@ -720,6 +746,10 @@ int BankAccount::put_data(int id, int size, double *data) {
     messages[4] = MIN_SCALE; // to indicate that it is defined by id
     MPI_Send(messages, 5, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     MPI_Send(data, size, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank);
+#else
+    (void)id;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -737,6 +767,10 @@ int BankAccount::put_data(int id, int size, ComplexDouble *data) {
     messages[4] = MIN_SCALE; // to indicate that it is defined by id
     MPI_Send(messages, 5, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     MPI_Send(data, size, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank);
+#else
+    (void)id;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -756,6 +790,10 @@ int BankAccount::put_data(NodeIndex<3> nIdx, int size, double *data) {
     int id = std::abs(nIdx.getTranslation(0) + nIdx.getTranslation(1) + nIdx.getTranslation(2));
     MPI_Send(messages, 7, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     MPI_Send(data, size, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank);
+#else
+    (void)nIdx;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -775,6 +813,10 @@ int BankAccount::put_data(NodeIndex<3> nIdx, int size, ComplexDouble *data) {
     int id = std::abs(nIdx.getTranslation(0) + nIdx.getTranslation(1) + nIdx.getTranslation(2));
     MPI_Send(messages, 7, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     MPI_Send(data, size, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank);
+#else
+    (void)nIdx;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -790,6 +832,10 @@ int BankAccount::get_data(int id, int size, double *data) {
     messages[3] = MIN_SCALE;
     MPI_Send(messages, 4, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     MPI_Recv(data, size, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank, &status);
+#else
+    (void)id;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -806,6 +852,10 @@ int BankAccount::get_data(int id, int size, ComplexDouble *data) {
     MPI_Send(messages, 4, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     // fetch as twice as many doubles
     MPI_Recv(data, size * 2, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank, &status);
+#else
+    (void)id;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -825,6 +875,10 @@ int BankAccount::get_data(NodeIndex<3> nIdx, int size, double *data) {
     messages[6] = nIdx.getTranslation(2);
     MPI_Send(messages, 7, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     MPI_Recv(data, size, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank, &status);
+#else
+    (void)nIdx;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -845,6 +899,10 @@ int BankAccount::get_data(NodeIndex<3> nIdx, int size, ComplexDouble *data) {
     MPI_Send(messages, 7, MPI_INT, bankmaster[id % bank_size], 0, comm_bank);
     // fetch as twice as many doubles
     MPI_Recv(data, size * 2, MPI_DOUBLE, bankmaster[id % bank_size], 1, comm_bank, &status);
+#else
+    (void)nIdx;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -861,6 +919,11 @@ int BankAccount::put_nodedata(int id, int nodeid, int size, double *data) {
     messages[4] = size;   // size of this data
     MPI_Send(messages, 5, MPI_INT, bankmaster[nodeid % bank_size], 0, comm_bank);
     MPI_Send(data, size, MPI_DOUBLE, bankmaster[nodeid % bank_size], 1, comm_bank);
+#else
+    (void)id;
+    (void)nodeid;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -878,6 +941,11 @@ int BankAccount::put_nodedata(int id, int nodeid, int size, ComplexDouble *data)
     messages[4] = 2 * size; // size of this data
     MPI_Send(messages, 5, MPI_INT, bankmaster[nodeid % bank_size], 0, comm_bank);
     MPI_Send(data, 2 * size, MPI_DOUBLE, bankmaster[nodeid % bank_size], 1, comm_bank);
+#else
+    (void)id;
+    (void)nodeid;
+    (void)size;
+    (void)data;
 #endif
     return 1;
 }
@@ -895,6 +963,12 @@ int BankAccount::get_nodedata(int id, int nodeid, int size, double *data, std::v
     messages[4] = size;   // expected size of data
     MPI_Send(messages, 5, MPI_INT, bankmaster[nodeid % bank_size], 0, comm_bank);
     MPI_Recv(data, size, MPI_DOUBLE, bankmaster[nodeid % bank_size], 3, comm_bank, &status);
+#else
+    (void)id;
+    (void)nodeid;
+    (void)size;
+    (void)data;
+    (void)idVec;
 #endif
     return 1;
 }
@@ -912,6 +986,12 @@ int BankAccount::get_nodedata(int id, int nodeid, int size, ComplexDouble *data,
     messages[4] = size;   // expected size of data
     MPI_Send(messages, 5, MPI_INT, bankmaster[nodeid % bank_size], 0, comm_bank);
     MPI_Recv(data, size, MPI_DOUBLE, bankmaster[nodeid % bank_size], 3, comm_bank, &status);
+#else
+    (void)id;
+    (void)nodeid;
+    (void)size;
+    (void)data;
+    (void)idVec;
 #endif
     return 1;
 }
@@ -932,6 +1012,10 @@ int BankAccount::get_nodeblock(int nodeid, double *data, std::vector<int> &idVec
     int size = metadata_block[2];
     if (size > 0) MPI_Recv(idVec.data(), metadata_block[1], MPI_INT, bankmaster[nodeid % bank_size], 2, comm_bank, &status);
     if (size > 0) MPI_Recv(data, size, MPI_DOUBLE, bankmaster[nodeid % bank_size], 3, comm_bank, &status);
+#else
+    (void)nodeid;
+    (void)data;
+    (void)idVec;
 #endif
     return 1;
 }
@@ -952,6 +1036,10 @@ int BankAccount::get_nodeblock(int nodeid, ComplexDouble *data, std::vector<int>
     int size = metadata_block[2];
     if (size > 0) MPI_Recv(idVec.data(), metadata_block[1], MPI_INT, bankmaster[nodeid % bank_size], 2, comm_bank, &status);
     if (size > 0) MPI_Recv(data, size, MPI_DOUBLE, bankmaster[nodeid % bank_size], 3, comm_bank, &status);
+#else
+    (void)nodeid;
+    (void)data;
+    (void)idVec;
 #endif
     return 1;
 }
@@ -973,6 +1061,11 @@ int BankAccount::get_orbblock(int orbid, double *&data, std::vector<int> &nodeid
     if (totsize > 0) MPI_Recv(nodeidVec.data(), metadata_block[1], MPI_INT, bankmaster[nodeid % bank_size], 2, comm_bank, &status);
     data = new double[totsize];
     if (totsize > 0) MPI_Recv(data, totsize, MPI_DOUBLE, bankmaster[nodeid % bank_size], 3, comm_bank, &status);
+#else
+    (void)orbid;
+    (void)data;
+    (void)nodeidVec;
+    (void)bankstart;
 #endif
     return 1;
 }
@@ -994,6 +1087,11 @@ int BankAccount::get_orbblock(int orbid, ComplexDouble *&data, std::vector<int>
     if (totsize > 0) MPI_Recv(nodeidVec.data(), metadata_block[1], MPI_INT, bankmaster[nodeid % bank_size], 2, comm_bank, &status);
     data = new ComplexDouble[totsize / 2];
     if (totsize > 0) MPI_Recv(data, totsize, MPI_DOUBLE, bankmaster[nodeid % bank_size], 3, comm_bank, &status);
+#else
+    (void)orbid;
+    (void)data;
+    (void)nodeidVec;
+    (void)bankstart;
 #endif
     return 1;
 }
@@ -1062,6 +1160,9 @@ void TaskManager::put_readytask(int i, int j) {
     messages[2] = i;
     messages[3] = j;
     MPI_Send(messages, message_size, MPI_INT, task_bank, 0, comm_bank);
+#else
+    (void)i;
+    (void)j;
 #endif
 }
 
@@ -1075,6 +1176,9 @@ void TaskManager::del_readytask(int i, int j) {
     messages[2] = i;
     messages[3] = j;
     MPI_Send(messages, message_size, MPI_INT, task_bank, 0, comm_bank);
+#else
+    (void)i;
+    (void)j;
 #endif
 }
 
@@ -1095,6 +1199,9 @@ std::vector<int> TaskManager::get_readytask(int i, int del) {
         readytasks.resize(nready);
         MPI_Recv(readytasks.data(), nready, MPI_INT, task_bank, 845, comm_bank, &status);
     }
+#else
+    (void)i;
+    (void)del;
 #endif
     return readytasks;
 }
diff --git a/src/utils/CompFunction.cpp b/src/utils/CompFunction.cpp
index 3cddb3588..2b681def1 100644
--- a/src/utils/CompFunction.cpp
+++ b/src/utils/CompFunction.cpp
@@ -608,7 +608,7 @@ template <int D> void linear_combination(CompFunction<D> &out, const std::vector
     out.func_ptr->data = inp[0].func_ptr->data;
     out.func_ptr->data.shared = share; // we don' inherit the shareness
     bool iscomplex = false;
-    for (int i = 0; i < inp.size(); i++)
+    for (size_t i = 0; i < inp.size(); i++)
         if (inp[i].iscomplex() or abs(c[i].imag()) > MachineZero) iscomplex = true;
     if (iscomplex) {
         out.func_ptr->data.iscomplex = 1;
@@ -618,7 +618,7 @@ template <int D> void linear_combination(CompFunction<D> &out, const std::vector
     for (int comp = 0; comp < inp[0].Ncomp(); comp++) {
         if (not iscomplex) {
             FunctionTreeVector<D, double> fvec; // one component vector
-            for (int i = 0; i < inp.size(); i++) {
+            for (size_t i = 0; i < inp.size(); i++) {
                 if (std::norm(c[i]) < thrs) continue;
                 if (inp[i].getNNodes() == 0 or inp[i].CompD[comp]->getSquareNorm() < thrs) continue;
                 fvec.push_back(std::make_tuple(c[i].real(), inp[i].CompD[comp]));
@@ -637,7 +637,7 @@ template <int D> void linear_combination(CompFunction<D> &out, const std::vector
             }
         } else {
             FunctionTreeVector<D, ComplexDouble> fvec; // one component vector
-            for (int i = 0; i < inp.size(); i++) {
+            for (size_t i = 0; i < inp.size(); i++) {
                 if (inp[i].isreal()) {
                     inp[i].CompC[comp] = inp[i].CompD[comp]->CopyTreeToComplex();
                     delete inp[i].CompD[comp];
@@ -817,6 +817,9 @@ template <int D> void multiply(CompFunction<D> &out, CompFunction<D> &inp_a, Rep
  */
 template <int D> void multiply(CompFunction<D> &out, CompFunction<D> &inp_a, RepresentableFunction<D, ComplexDouble> &f, double prec, int nrefine, bool conjugate) {
     MSG_ABORT("Not implemented");
+    (void)f;
+    (void)prec;
+    (void)nrefine;
     if (inp_a.Ncomp() > 1) MSG_ABORT("Not implemented");
     if (inp_a.iscomplex() != 1) MSG_ABORT("Not implemented");
     if (conjugate) MSG_ABORT("Not implemented");
@@ -961,7 +964,7 @@ CompFunctionVector::CompFunctionVector(int N)
     vecMRA = defaultCompMRA<3>;
 }
 void CompFunctionVector::distribute() {
-    for (int i = 0; i < this->size(); i++) (*this)[i].func_ptr->rank = i;
+    for (size_t i = 0; i < this->size(); i++) (*this)[i].func_ptr->rank = i;
 }
 
 /** @brief Make a linear combination of functions
@@ -1106,8 +1109,8 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
 
             // 4c) rotate this node
             ComplexMatrix Un(orbjVec.size(), orbiVec.size()); // chunk of U, with reorganized indices
-            for (int i = 0; i < orbiVec.size(); i++) {        // loop over rotated orbitals
-                for (int j = 0; j < orbjVec.size(); j++) { Un(j, i) = U(orbjVec[j], orbiVec[i]); }
+            for (size_t i = 0; i < orbiVec.size(); i++) {        // loop over rotated orbitals
+                for (size_t j = 0; j < orbjVec.size(); j++) { Un(j, i) = U(orbjVec[j], orbiVec[i]); }
             }
             ComplexMatrix rotatedCoeff(csize, orbiVec.size());
             // HERE IT HAPPENS!
@@ -1118,7 +1121,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
             // for now we allocate in buffer, in future could be directly allocated in the final trees
             double thres = prec * prec * scalefac_ref[n] * scalefac_ref[n];
             // make all norms:
-            for (int i = 0; i < orbiVec.size(); i++) {
+            for (size_t i = 0; i < orbiVec.size(); i++) {
                 // check if parent must be split
                 if (parindexVec_ref[n] == -1 or split_serial(orbiVec[i], ix2coef_ref[parindexVec_ref[n]])) {
                     // mark this node for this orbital for later split
@@ -1190,8 +1193,8 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
             ComplexMatrix Un(orbjVec.size(), orbiVec.size());
             ComplexMatrix rotatedCoeff(csize, orbiVec.size());
 
-            for (int i = 0; i < orbiVec.size(); i++) {     // loop over included rotated real and imag part of orbitals
-                for (int j = 0; j < orbjVec.size(); j++) { // loop over input orbital, possibly imaginary parts
+            for (size_t i = 0; i < orbiVec.size(); i++) {     // loop over included rotated real and imag part of orbitals
+                for (size_t j = 0; j < orbjVec.size(); j++) { // loop over input orbital, possibly imaginary parts
                     Un(j, i) = U(orbjVec[j], orbiVec[i]);
                 }
             }
@@ -1202,7 +1205,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
 
             // 3c) find which orbitals need to further refine this node, and store rotated node (after each other while
             // in cache).
-            for (int i = 0; i < orbiVec.size(); i++) { // loop over rotated orbitals
+            for (size_t i = 0; i < orbiVec.size(); i++) { // loop over rotated orbitals
                 needsplit[orbiVec[i]] = -1.0;          // default, do not split
                 // check if this node/orbital needs further refinement
                 double wnorm = 0.0;
@@ -1243,7 +1246,7 @@ void rotate_cplx(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVe
                 nodesRotated.get_orbblock(j, dataVec, nodeidVec, ibank);
                 if (nodeidVec.size() > 0) pointerstodelete.push_back(dataVec);
                 int shift = 0;
-                for (int n = 0; n < nodeidVec.size(); n++) {
+                for (size_t n = 0; n < nodeidVec.size(); n++) {
                     assert(nodeidVec[n] - max_ix >= 0);                // unrotated nodes have been deleted
                     assert(ix2coef.count(nodeidVec[n] - max_ix) == 0); // each nodeid treated once
                     ix2coef[nodeidVec[n] - max_ix] = ix++;
@@ -1395,8 +1398,8 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector
 
             // 4c) rotate this node
             DoubleMatrix Un(orbjVec.size(), orbiVec.size()); // chunk of U, with reorganized indices
-            for (int i = 0; i < orbiVec.size(); i++) {       // loop over rotated orbitals
-                for (int j = 0; j < orbjVec.size(); j++) { Un(j, i) = std::real(U(orbjVec[j], orbiVec[i])); }
+            for (size_t i = 0; i < orbiVec.size(); i++) {       // loop over rotated orbitals
+                for (size_t j = 0; j < orbjVec.size(); j++) { Un(j, i) = std::real(U(orbjVec[j], orbiVec[i])); }
             }
             DoubleMatrix rotatedCoeff(csize, orbiVec.size());
             // HERE IT HAPPENS!
@@ -1406,7 +1409,7 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector
             // for now we allocate in buffer, in future could be directly allocated in the final trees
             double thres = prec * prec * scalefac_ref[n] * scalefac_ref[n];
             // make all norms:
-            for (int i = 0; i < orbiVec.size(); i++) {
+            for (size_t i = 0; i < orbiVec.size(); i++) {
                 // check if parent must be split
                 if (parindexVec_ref[n] == -1 or split_serial(orbiVec[i], ix2coef_ref[parindexVec_ref[n]])) {
                     // mark this node for this orbital for later split
@@ -1478,8 +1481,8 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector
             DoubleMatrix Un(orbjVec.size(), orbiVec.size());
             DoubleMatrix rotatedCoeff(csize, orbiVec.size());
 
-            for (int i = 0; i < orbiVec.size(); i++) {     // loop over included rotated real and imag part of orbitals
-                for (int j = 0; j < orbjVec.size(); j++) { // loop over input orbital, possibly imaginary parts
+            for (size_t i = 0; i < orbiVec.size(); i++) {     // loop over included rotated real and imag part of orbitals
+                for (size_t j = 0; j < orbjVec.size(); j++) { // loop over input orbital, possibly imaginary parts
                     Un(j, i) = std::real(U(orbjVec[j], orbiVec[i]));
                 }
             }
@@ -1489,7 +1492,7 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector
 
             // 3c) find which orbitals need to further refine this node, and store rotated node (after each other while
             // in cache).
-            for (int i = 0; i < orbiVec.size(); i++) { // loop over rotated orbitals
+            for (size_t i = 0; i < orbiVec.size(); i++) { // loop over rotated orbitals
                 needsplit[orbiVec[i]] = -1.0;          // default, do not split
                 // check if this node/orbital needs further refinement
                 double wnorm = 0.0;
@@ -1533,7 +1536,7 @@ void rotate(CompFunctionVector &Phi, const ComplexMatrix &U, CompFunctionVector
                 nodesRotated.get_orbblock(j, dataVec, nodeidVec, ibank);
                 if (nodeidVec.size() > 0) pointerstodelete.push_back(dataVec);
                 int shift = 0;
-                for (int n = 0; n < nodeidVec.size(); n++) {
+                for (size_t n = 0; n < nodeidVec.size(); n++) {
                     assert(nodeidVec[n] - max_ix >= 0);                // unrotated nodes have been deleted
                     assert(ix2coef.count(nodeidVec[n] - max_ix) == 0); // each nodeid treated once
                     ix2coef[nodeidVec[n] - max_ix] = ix++;
@@ -1783,7 +1786,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
             }
 
             // 3f) save multiplied nodes
-            for (int i = 0; i < orbjVec.size(); i++) {
+            for (size_t i = 0; i < orbjVec.size(); i++) {
 #pragma omp critical
                 {
                     ix2coef[orbjVec[i]][node_ix] = coeffpVec[orbjVec[i]].size();
@@ -1839,7 +1842,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
             coeffBlock.conservativeResize(Eigen::NoChange, orbjVec.size()); // keep only used part
             DoubleMatrix MultipliedCoeff(nCoefs, orbjVec.size());
             // 3c) transform to grid
-            for (int j = 0; j < orbjVec.size(); j++) { // TODO: transform all j at once ?
+            for (size_t j = 0; j < orbjVec.size(); j++) { // TODO: transform all j at once ?
                 // TODO: select only nodes that are end nodes?
                 node.attachCoefs(coeffBlock.col(j).data());
                 node.mwTransform(Reconstruction);
@@ -1902,7 +1905,7 @@ CompFunctionVector multiply(CompFunctionVector &Phi, RepresentableFunction<3> &f
                 nodesMultiplied.get_orbblock(j, dataVec, nodeidVec, ibank);
                 if (nodeidVec.size() > 0) pointerstodelete.push_back(dataVec);
                 int shift = 0;
-                for (int n = 0; n < nodeidVec.size(); n++) {
+                for (size_t n = 0; n < nodeidVec.size(); n++) {
                     assert(nodeidVec[n] - max_ix >= 0);                // unmultiplied nodes have been deleted
                     assert(ix2coef.count(nodeidVec[n] - max_ix) == 0); // each nodeid treated once
                     ix2coef[nodeidVec[n] - max_ix] = ix++;
@@ -2064,8 +2067,8 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) {
                 if (orbVec.size() > 0) {
                     ComplexMatrix S_temp(orbVec.size(), orbVec.size());
                     S_temp.noalias() = coeffBlock.transpose().conjugate() * coeffBlock;
-                    for (int i = 0; i < orbVec.size(); i++) {
-                        for (int j = 0; j < orbVec.size(); j++) {
+                    for (size_t i = 0; i < orbVec.size(); i++) {
+                        for (size_t j = 0; j < orbVec.size(); j++) {
                             if (BraKet[orbVec[i]].func_ptr->data.n1[0] != BraKet[orbVec[j]].func_ptr->data.n1[0] and BraKet[orbVec[i]].func_ptr->data.n1[0] != 0 and
                                 BraKet[orbVec[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2081,8 +2084,8 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &BraKet) {
                     ComplexMatrix S_temp(orbVec.size(), orbVec.size());
                     coeffBlock.conservativeResize(Eigen::NoChange, orbVec.size());
                     S_temp.noalias() = coeffBlock.transpose().conjugate() * coeffBlock;
-                    for (int i = 0; i < orbVec.size(); i++) {
-                        for (int j = 0; j < orbVec.size(); j++) {
+                    for (size_t i = 0; i < orbVec.size(); i++) {
+                        for (size_t j = 0; j < orbVec.size(); j++) {
                             if (BraKet[orbVec[i]].func_ptr->data.n1[0] != BraKet[orbVec[j]].func_ptr->data.n1[0] and BraKet[orbVec[i]].func_ptr->data.n1[0] != 0 and
                                 BraKet[orbVec[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2212,8 +2215,8 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &BraKet) {
                 if (orbVec.size() > 0) {
                     ComplexMatrix S_temp(orbVec.size(), orbVec.size());
                     S_temp.noalias() = coeffBlock.transpose() * coeffBlock;
-                    for (int i = 0; i < orbVec.size(); i++) {
-                        for (int j = 0; j < orbVec.size(); j++) {
+                    for (size_t i = 0; i < orbVec.size(); i++) {
+                        for (size_t j = 0; j < orbVec.size(); j++) {
                             if (BraKet[orbVec[i]].func_ptr->data.n1[0] != BraKet[orbVec[j]].func_ptr->data.n1[0] and BraKet[orbVec[i]].func_ptr->data.n1[0] != 0 and
                                 BraKet[orbVec[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2229,8 +2232,8 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &BraKet) {
                     DoubleMatrix S_temp(orbVec.size(), orbVec.size());
                     coeffBlock.conservativeResize(Eigen::NoChange, orbVec.size());
                     S_temp.noalias() = coeffBlock.transpose() * coeffBlock;
-                    for (int i = 0; i < orbVec.size(); i++) {
-                        for (int j = 0; j < orbVec.size(); j++) {
+                    for (size_t i = 0; i < orbVec.size(); i++) {
+                        for (size_t j = 0; j < orbVec.size(); j++) {
                             if (BraKet[orbVec[i]].func_ptr->data.n1[0] != BraKet[orbVec[j]].func_ptr->data.n1[0] and BraKet[orbVec[i]].func_ptr->data.n1[0] != 0 and
                                 BraKet[orbVec[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2284,14 +2287,14 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
         // Here, we keep the real trees of the inputs, as well as creating the complex ones, to avoid copying the trees back after computing the overlap matrix.
         // We restore the original state of the inputs after multiplication.
         if (braisreal) {
-            for (int i = 0; i < Bra.size(); i++) {
+            for (size_t i = 0; i < Bra.size(); i++) {
                 Bra[i].CompC[0] = Bra[i].CompD[0]->CopyTreeToComplex();
                 Bra[i].func_ptr->iscomplex = 1;
                 Bra[i].func_ptr->isreal = 0;
             }
         }
         if (ketisreal) {
-            for (int i = 0; i < Ket.size(); i++) {
+            for (size_t i = 0; i < Ket.size(); i++) {
                 Ket[i].CompC[0] = Ket[i].CompD[0]->CopyTreeToComplex();
                 Ket[i].func_ptr->iscomplex = 1;
                 Ket[i].func_ptr->isreal = 0;
@@ -2438,8 +2441,8 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
                         S_temp.noalias() = coeffBlockBra * coeffBlockKet.transpose();
                     } else
                         MSG_ABORT("Unexpected case");
-                    for (int i = 0; i < orbVecBra.size(); i++) {
-                        for (int j = 0; j < orbVecKet.size(); j++) {
+                    for (size_t i = 0; i < orbVecBra.size(); i++) {
+                        for (size_t j = 0; j < orbVecKet.size(); j++) {
                             if (Bra[orbVecBra[i]].func_ptr->data.n1[0] != Ket[orbVecKet[j]].func_ptr->data.n1[0] and Bra[orbVecBra[i]].func_ptr->data.n1[0] != 0 and
                                 Ket[orbVecKet[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2471,8 +2474,8 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
                     } else
                         MSG_ABORT("Unexpected case");
 
-                    for (int i = 0; i < orbVecBra.size(); i++) {
-                        for (int j = 0; j < orbVecKet.size(); j++) {
+                    for (size_t i = 0; i < orbVecBra.size(); i++) {
+                        for (size_t j = 0; j < orbVecKet.size(); j++) {
                             if (Bra[orbVecBra[i]].func_ptr->data.n1[0] != Ket[orbVecKet[j]].func_ptr->data.n1[0] and Bra[orbVecBra[i]].func_ptr->data.n1[0] != 0 and
                                 Ket[orbVecKet[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2513,7 +2516,7 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
 
     // restore original tree by deleting the temporary complex tree. The real tree still exists, but is not used in the computation.
     if (braisreal) {
-        for (int i = 0; i < Bra.size(); i++) {
+        for (size_t i = 0; i < Bra.size(); i++) {
             delete Bra[i].CompC[0];
             Bra[i].CompC[0] = nullptr;
             Bra[i].func_ptr->iscomplex = 0;
@@ -2521,7 +2524,7 @@ ComplexMatrix calc_overlap_matrix_cplx(CompFunctionVector &Bra, CompFunctionVect
         }
     }
     if (ketisreal) {
-        for (int i = 0; i < Ket.size(); i++) {
+        for (size_t i = 0; i < Ket.size(); i++) {
             delete Ket[i].CompC[0];
             Ket[i].CompC[0] = nullptr;
             Ket[i].func_ptr->iscomplex = 0;
@@ -2657,8 +2660,8 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K
                     DoubleMatrix S_temp(orbVecBra.size(), orbVecKet.size());
                     S_temp.noalias() = coeffBlockBra.transpose() * coeffBlockKet;
 
-                    for (int i = 0; i < orbVecBra.size(); i++) {
-                        for (int j = 0; j < orbVecKet.size(); j++) {
+                    for (size_t i = 0; i < orbVecBra.size(); i++) {
+                        for (size_t j = 0; j < orbVecKet.size(); j++) {
                             if (Bra[orbVecBra[i]].func_ptr->data.n1[0] != Ket[orbVecKet[j]].func_ptr->data.n1[0] and Bra[orbVecBra[i]].func_ptr->data.n1[0] != 0 and
                                 Ket[orbVecKet[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2680,8 +2683,8 @@ ComplexMatrix calc_overlap_matrix(CompFunctionVector &Bra, CompFunctionVector &K
                     coeffBlockBra.conservativeResize(Eigen::NoChange, orbVecBra.size());
                     coeffBlockKet.conservativeResize(Eigen::NoChange, orbVecKet.size());
                     S_temp.noalias() = coeffBlockBra.transpose() * coeffBlockKet;
-                    for (int i = 0; i < orbVecBra.size(); i++) {
-                        for (int j = 0; j < orbVecKet.size(); j++) {
+                    for (size_t i = 0; i < orbVecBra.size(); i++) {
+                        for (size_t j = 0; j < orbVecKet.size(); j++) {
                             if (Bra[orbVecBra[i]].func_ptr->data.n1[0] != Ket[orbVecKet[j]].func_ptr->data.n1[0] and Bra[orbVecBra[i]].func_ptr->data.n1[0] != 0 and
                                 Ket[orbVecKet[j]].func_ptr->data.n1[0] != 0)
                                 continue;
@@ -2751,6 +2754,7 @@ void orthogonalize(double prec, CompFunctionVector &Bra, CompFunctionVector &Ket
  *
  */
 template <int D> void orthogonalize(double prec, CompFunction<D> &Bra, CompFunction<D> &Ket) {
+    (void) prec;
     ComplexDouble overlap = dot(Bra, Ket);
     double sq_norm = Ket.getSquareNorm();
     for (int i = 0; i < Bra.Ncomp(); i++) {
diff --git a/src/utils/CompFunction.h b/src/utils/CompFunction.h
index 7f5cf463d..25b358119 100644
--- a/src/utils/CompFunction.h
+++ b/src/utils/CompFunction.h
@@ -166,7 +166,7 @@ template <int D> class CompFunction {
     void setReal(FunctionTree<D, double> *tree, int i = 0);
     void setCplx(FunctionTree<D, ComplexDouble> *tree, int i = 0);
     void setRank(int i) { func_ptr->rank = i; };
-    const int getRank() const { return func_ptr->rank; };
+    int getRank() const { return func_ptr->rank; };
     void add(ComplexDouble c, CompFunction<D> inp);
 
     int crop(double prec, bool absPrec = true);
@@ -184,7 +184,7 @@ template <int D> class CompFunction {
 
     // NB: All below should be revised. Now only for backwards compatibility to ComplexFunction class
 
-    void free(int type) { free(); }
+    void free(int type) { (void)type; free(); }
     bool hasReal() const { return isreal(); }
     bool hasImag() const { return iscomplex(); }
     bool isShared() const { return share(); }
diff --git a/src/utils/Plotter.cpp b/src/utils/Plotter.cpp
index c29b3ee2e..37bcb3ce9 100644
--- a/src/utils/Plotter.cpp
+++ b/src/utils/Plotter.cpp
@@ -298,16 +298,21 @@ template <int D, typename T> void Plotter<D, T>::writeData(const Eigen::MatrixXd
 
 // Specialized for D=3 below
 template <int D, typename T> void Plotter<D, T>::writeCube(const std::array<int, 3> &npts, const Eigen::Matrix<T, Eigen::Dynamic, 1> &values) {
+    (void)npts;
+    (void)values;
     NOT_IMPLEMENTED_ABORT
 }
 
 // Specialized for D=3 below
 template <int D, typename T> void Plotter<D, T>::writeNodeGrid(const MWNode<D, T> &node, const std::string &color) {
+    (void)node;
+    (void)color;
     NOT_IMPLEMENTED_ABORT
 }
 
 // Specialized for D=3 below
 template <int D, typename T> void Plotter<D, T>::writeGrid(const MWTree<D, T> &tree) {
+    (void)tree;
     NOT_IMPLEMENTED_ABORT
 }
 
@@ -431,11 +436,11 @@ template <> void Plotter<3>::writeGrid(const MWTree<3> &tree) {
     o.precision(6);
     std::string rootColor = " 1 1 1 0 ";
     std::string color = " 0 0 1 1 ";
-    for (auto i = 0; i < tree.getRootBox().size(); i++) {
+    for (size_t i = 0; i < tree.getRootBox().size(); i++) {
         const MWNode<3> &rootNode = tree.getRootMWNode(i);
         writeNodeGrid(rootNode, rootColor);
     }
-    for (auto i = 0; i < tree.getNEndNodes(); i++) {
+    for (int i = 0; i < tree.getNEndNodes(); i++) {
         const MWNode<3> &node = tree.getEndMWNode(i);
         writeNodeGrid(node, color);
     }
diff --git a/src/utils/math_utils.cpp b/src/utils/math_utils.cpp
index 69a13f300..647ae0801 100644
--- a/src/utils/math_utils.cpp
+++ b/src/utils/math_utils.cpp
@@ -95,7 +95,7 @@ double math_utils::binomial_coeff(int n, int j) {
 
 VectorXd math_utils::get_binomial_coefs(unsigned int order) {
     VectorXd coefs = VectorXd::Ones(order + 1);
-    for (int k = 0; k <= order; k++) { coefs[k] = math_utils::binomial_coeff(order, k); }
+    for (unsigned int k = 0; k <= order; k++) { coefs[k] = math_utils::binomial_coeff(order, k); }
     return coefs;
 }
 
diff --git a/src/utils/mpi_utils.cpp b/src/utils/mpi_utils.cpp
index 77526375b..cd8b16939 100644
--- a/src/utils/mpi_utils.cpp
+++ b/src/utils/mpi_utils.cpp
@@ -60,6 +60,9 @@ SharedMemory<T>::SharedMemory(mrcpp::mpi_comm comm, int sh_size)
     MPI_Win_fence(0, this->sh_win);
     this->sh_max_ptr = this->sh_start_ptr + qsize / sizeof(T);
     this->sh_end_ptr = this->sh_start_ptr;
+#else
+    (void)comm;
+    (void)sh_size;
 #endif
 }
 
@@ -105,6 +108,13 @@ template <int D, typename T> void send_tree(FunctionTree<D, T> &tree, int dst, i
         if (coeff) MPI_Send(allocator.getCoefChunk(iChunk), allocator.getCoefChunkSize(), MPI_BYTE, dst, tag + iChunk + 1001, comm);
     }
     println(10, " Time send                   " << std::setw(30) << t1.elapsed());
+#else
+    (void)tree;
+    (void)dst;
+    (void)tag;
+    (void)comm;
+    (void)nChunks;
+    (void)coeff;
 #endif
 }
 
@@ -143,6 +153,13 @@ template <int D, typename T> void recv_tree(FunctionTree<D, T> &tree, int src, i
     Timer t2;
     allocator.reassemble();
     println(10, " Time rewrite pointers       " << std::setw(30) << t2.elapsed());
+#else
+    (void)tree;
+    (void)src;
+    (void)tag;
+    (void)comm;
+    (void)nChunks;
+    (void)coeff;
 #endif
 }
 
@@ -194,6 +211,11 @@ template <int D, typename T> void share_tree(FunctionTree<D, T> &tree, int src,
         }
     }
     println(10, " Time share                  " << std::setw(30) << t1.elapsed());
+#else
+    (void)tree;
+    (void)src;
+    (void)tag;
+    (void)comm;
 #endif
 }
 template class SharedMemory<double>;
diff --git a/src/utils/parallel.cpp b/src/utils/parallel.cpp
index 510e34a1e..f77cd8c42 100644
--- a/src/utils/parallel.cpp
+++ b/src/utils/parallel.cpp
@@ -264,6 +264,8 @@ void finalize() {
 void barrier(MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     MPI_Barrier(comm);
+#else
+    (void)comm;
 #endif
 }
 
@@ -306,6 +308,9 @@ void allreduce_vector(IntVector &vec, MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     int N = vec.size();
     MPI_Allreduce(MPI_IN_PLACE, vec.data(), N, MPI_INT, MPI_SUM, comm);
+#else
+    (void)vec;
+    (void)comm;
 #endif
 }
 
@@ -314,6 +319,9 @@ void allreduce_vector(DoubleVector &vec, MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     int N = vec.size();
     MPI_Allreduce(MPI_IN_PLACE, vec.data(), N, MPI_DOUBLE, MPI_SUM, comm);
+#else
+    (void)vec;
+    (void)comm;
 #endif
 }
 
@@ -322,6 +330,9 @@ void allreduce_vector(ComplexVector &vec, MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     int N = vec.size();
     MPI_Allreduce(MPI_IN_PLACE, vec.data(), N, MPI_C_DOUBLE_COMPLEX, MPI_SUM, comm);
+#else
+    (void)vec;
+    (void)comm;
 #endif
 }
 
@@ -330,6 +341,9 @@ void allreduce_matrix(IntMatrix &mat, MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     int N = mat.size();
     MPI_Allreduce(MPI_IN_PLACE, mat.data(), N, MPI_INT, MPI_SUM, comm);
+#else
+    (void)mat;
+    (void)comm;
 #endif
 }
 
@@ -338,6 +352,9 @@ void allreduce_matrix(DoubleMatrix &mat, MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     int N = mat.size();
     MPI_Allreduce(MPI_IN_PLACE, mat.data(), N, MPI_DOUBLE, MPI_SUM, comm);
+#else
+    (void)mat;
+    (void)comm;
 #endif
 }
 
@@ -346,6 +363,9 @@ void allreduce_matrix(ComplexMatrix &mat, MPI_Comm comm) {
 #ifdef MRCPP_HAS_MPI
     int N = mat.size();
     MPI_Allreduce(MPI_IN_PLACE, mat.data(), N, MPI_C_DOUBLE_COMPLEX, MPI_SUM, comm);
+#else
+    (void)mat;
+    (void)comm;
 #endif
 }
 
@@ -366,6 +386,11 @@ void send_function(const CompFunction<3> &func, int dst, int tag, MPI_Comm comm)
         else
             mrcpp::send_tree(*func.CompC[i], dst, tag, comm, func.Nchunks()[i]);
     }
+#else
+    (void)func;
+    (void)dst;
+    (void)tag;
+    (void)comm;
 #endif
 }
 
@@ -382,6 +407,11 @@ void recv_function(CompFunction<3> &func, int src, int tag, MPI_Comm comm) {
         else
             mrcpp::recv_tree(*func.CompC[i], src, tag, comm, func.Nchunks()[i]);
     }
+#else
+    (void)func;
+    (void)src;
+    (void)tag;
+    (void)comm;
 #endif
 }
 
@@ -395,6 +425,10 @@ void share_function(CompFunction<3> &func, int src, int tag, MPI_Comm comm) {
             else
                 mrcpp::share_tree(*func.CompC[comp], src, tag, comm);
         }
+#else
+        (void)src;
+        (void)tag;
+        (void)comm;
 #endif
     }
 }
@@ -438,6 +472,10 @@ void reduce_function(double prec, CompFunction<3> &func, MPI_Comm comm) {
         fac *= 2;
     }
     MPI_Barrier(comm);
+#else
+    (void)prec;
+    (void)func;
+    (void)comm;
 #endif
 }
 
@@ -479,6 +517,9 @@ template <typename T> void reduce_Tree_noCoeff(mrcpp::FunctionTree<3, T> &tree,
         fac *= 2;
     }
     MPI_Barrier(comm);
+#else
+    (void)tree;
+    (void)comm;
 #endif
 }
 
@@ -498,6 +539,8 @@ template <typename T> void allreduce_Tree_noCoeff(mrcpp::FunctionTree<3, T> &tre
 #ifdef MRCPP_HAS_MPI
     mrcpp::mpi::reduce_Tree_noCoeff(tree, comm_wrk);
     mrcpp::mpi::broadcast_Tree_noCoeff(tree, comm_wrk);
+#else
+    (void)comm;
 #endif
 }
 
@@ -518,6 +561,8 @@ template <typename T> void allreduce_Tree_noCoeff(mrcpp::FunctionTree<3, T> &tre
 #ifdef MRCPP_HAS_MPI
     mrcpp::mpi::reduce_Tree_noCoeff(tree, comm_wrk);
     mrcpp::mpi::broadcast_Tree_noCoeff(tree, comm_wrk);
+#else
+    (void)comm;
 #endif
 }
 
@@ -550,6 +595,9 @@ void broadcast_function(CompFunction<3> &func, MPI_Comm comm) {
         fac /= 2;
     }
     MPI_Barrier(comm);
+#else
+    (void)func;
+    (void)comm;
 #endif
 }
 
@@ -582,6 +630,9 @@ template <typename T> void broadcast_Tree_noCoeff(mrcpp::FunctionTree<3, T> &tre
         fac /= 2;
     }
     MPI_Barrier(comm);
+#else
+    (void)tree;
+    (void)comm;
 #endif
 }
 
diff --git a/src/utils/tree_utils.cpp b/src/utils/tree_utils.cpp
index b01c03347..4eafcf2fd 100644
--- a/src/utils/tree_utils.cpp
+++ b/src/utils/tree_utils.cpp
@@ -89,7 +89,7 @@ template <int D, typename T> void tree_utils::make_node_table(MWTree<D, T> &tree
     while (it.nextParent()) {
         MWNode<D, T> &node = it.getNode();
         if (node.getDepth() == 0) continue;
-        int depth = node.getDepth() + tree.getNNegScales();
+        size_t depth = node.getDepth() + tree.getNNegScales();
         // Add one more element
         if (depth + 1 > table.size()) table.push_back(MWNodeVector<D, T>());
         table[depth].push_back(&node);
@@ -97,7 +97,7 @@ template <int D, typename T> void tree_utils::make_node_table(MWTree<D, T> &tree
     it.init(tree);
     while (it.next()) {
         MWNode<D, T> &node = it.getNode();
-        int depth = node.getDepth() + tree.getNNegScales();
+        size_t depth = node.getDepth() + tree.getNNegScales();
         // Add one more element
         if (depth + 1 > table.size()) table.push_back(MWNodeVector<D, T>());
         table[depth].push_back(&node);
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 257bb1073..365c23f3b 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,4 +1,4 @@
-find_package( Catch2 3.0.1 CONFIG)
+find_package( Catch2 3.8.0 CONFIG)
 if( NOT Catch2_FOUND )
   message( STATUS "Could NOT Find Catch2 (Building v3.8.0 Locally)" )
   include( FetchContent )
diff --git a/tests/operators/helmholtz_operator.cpp b/tests/operators/helmholtz_operator.cpp
index ce669e27c..935af7d52 100644
--- a/tests/operators/helmholtz_operator.cpp
+++ b/tests/operators/helmholtz_operator.cpp
@@ -97,7 +97,7 @@ TEST_CASE("Helmholtz' kernel", "[init_helmholtz], [helmholtz_operator], [mw_oper
                 TreeBuilder<2> builder;
                 OperatorAdaptor adaptor(ccc_prec, oper_mra.getMaxScale());
 
-                for (int i = 0; i < K.size(); i++) {
+                for (size_t i = 0; i < K.size(); i++) {
                     FunctionTree<1> &kern_tree = get_func(K, i);
                     CrossCorrelationCalculator calculator(kern_tree);
 
diff --git a/tests/operators/poisson_operator.cpp b/tests/operators/poisson_operator.cpp
index e6ad04f69..1de777b65 100644
--- a/tests/operators/poisson_operator.cpp
+++ b/tests/operators/poisson_operator.cpp
@@ -95,7 +95,7 @@ TEST_CASE("Initialize Poisson operator", "[init_poisson], [poisson_operator], [m
                 TreeBuilder<2> builder;
                 OperatorAdaptor adaptor(ccc_prec, oper_mra.getMaxScale());
 
-                for (int i = 0; i < kern_vec.size(); i++) {
+                for (size_t i = 0; i < kern_vec.size(); i++) {
                     FunctionTree<1> &kern_tree = get_func(kern_vec, i);
                     CrossCorrelationCalculator calculator(kern_tree);