diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2bf9b83..4c7ff9c 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,228 +1,228 @@
 stages:
   - docker
   - lint
   - build
   - test
   - deploy
 
 variables:
   IMAGE_TAG: $CI_REGISTRY_IMAGE:$CI_COMMIT_REF_SLUG
   GIT_SUBMODULE_STRATEGY: recursive
   BUILD_DIR: build-release
 
 cache:
   key: "$CI_COMMIT_REF_SLUG"
 
 # ------------------------------------------------------------------------------
 
 .docker_build:
   stage: docker
   image: docker:19.03.12
   services:
     - docker:19.03.12-dind
   variables:
     # Use TLS https://docs.gitlab.com/ee/ci/docker/using_docker_build.html#tls-enabled
     DOCKER_HOST: tcp://docker:2376
     DOCKER_TLS_CERTDIR: "/certs"
     DEFAULT_IMAGE: $CI_REGISTRY_IMAGE:$CI_DEFAULT_BRANCH
   before_script:
     - docker login -u $CI_REGISTRY_USER -p $CI_REGISTRY_PASSWORD $CI_REGISTRY
   script:
     - docker pull $DEFAULT_IMAGE-$IMAGE_NAME || true
     - docker build --cache-from $DEFAULT_IMAGE-$IMAGE_NAME
                    -t $IMAGE_TAG-$IMAGE_NAME -f $DOCKERFILE .
     - docker push $IMAGE_TAG-$IMAGE_NAME
 
 docker build:debian:
   variables:
     IMAGE_NAME: debian-stable
     DOCKERFILE: tests/ci/docker/debian.mpi
   extends: .docker_build
 
 docker build:manylinux:
   variables:
     IMAGE_NAME: manylinux
     DOCKERFILE: tests/ci/docker/manylinux
   extends: .docker_build
 
 docker build:cuda:
   variables:
     IMAGE_NAME: cuda
     DOCKERFILE: tests/ci/docker/ubuntu_lts.cuda
   extends: .docker_build
 
 # ------------------------------------------------------------------------------
 
 .debian_stable:
   variables:
     output: ${CI_COMMIT_REF_SLUG}-debian-stable
   image: ${IMAGE_TAG}-debian-stable
 
 .manylinux:
   variables:
     output: ${CI_COMMIT_REF_SLUG}-manylinux
   image: ${IMAGE_TAG}-manylinux
 
 .cuda:
   variables:
     output: ${CI_COMMIT_REF_SLUG}-cuda
   image: ${IMAGE_TAG}-cuda
 
 # ------------------------------------------------------------------------------
 
 lint:
   stage: lint
   extends:
     - .debian_stable
   allow_failure: true
   script:
     - git remote remove upstream || true
     - git remote add upstream "$CI_PROJECT_URL"
     - git fetch upstream
     - '[ -z "$CI_MERGE_REQUEST_TARGET_BRANCH_NAME" ] && HEAD="$CI_COMMIT_BEFORE_SHA" || HEAD="upstream/$CI_MERGE_REQUEST_TARGET_BRANCH_NAME"'
     - echo "$HEAD"
-    - git-clang-format --diff "$HEAD" > clang_format.patch 2&>1
-    - grep "no modified files to format" clang_format.patch
+    - git-clang-format --diff "$HEAD" > clang_format.patch 2>&1
+    - grep "no modified files to format" clang_format.patch || grep "clang-format did not modify any files" clang_format.patch
   artifacts:
     when: on_failure
     paths:
       - clang_format.patch
 
 
 .build:
   stage: build
   variables:
     COMPILE_LOG: compilation.log
   script:
     - scons 2>&1 | tee $COMPILE_LOG
   artifacts:
     when: always
     paths:
       - build-setup.conf
       - $COMPILE_LOG
       - $BUILD_DIR
       - config.log
 
 .ccache:
   variables:
     CCACHE_BASEDIR: $CI_PROJECT_DIR/$BUILD_DIR
     CCACHE_DIR: $CI_PROJECT_DIR/.ccache
     CCACHE_NOHASDIR: 1
     CCACHE_COMPILERCHECK: content
     CXX: /usr/lib/ccache/g++
     OMPI_CXX: ${CXX}
   cache:
     key: ${output}
     policy: pull-push
     paths:
       - .ccache
   after_script:
     - ccache --show-stats || true
 
 build:mpi:
   extends:
     - .build
     - .debian_stable
     - .ccache
   before_script:
     - ccache --zero-stats || true
     - scons build_tests=True
             use_googletest=True
             build_python=True
             py_exec=python3
             use_mpi=True
             backend=omp
             fftw_threads=omp
             verbose=True -h
 
 build:cuda:
   extends:
     - .build
     - .cuda
   before_script:
     - scons build_tests=True
             use_googletest=True
             build_python=True
             py_exec=python3
             use_mpi=False
             backend=cuda
             fftw_threads=none
             verbose=True -h
 
 # ------------------------------------------------------------------------------
 
 test:mpi:
   stage: test
   dependencies:
     - build:mpi
   extends: .debian_stable
   variables:
     PYTHONPATH: $CI_PROJECT_DIR/$BUILD_DIR/python
     TESTS: $BUILD_DIR/tests
     JUNITXML: results.xml
     TESTS_LOG: tests.log
   script:
     - ls $PYTHONPATH
     - python3 -c 'import sys; print(sys.path)'
     - python3 -m pytest -vvv --last-failed --durations=0 --junitxml=$JUNITXML
       $TESTS 2>&1 | tee $TESTS_LOG
   after_script:
     - python3 -m pytest --cache-show || true
   artifacts:
     when: always
     paths:
       - $JUNITXML
       - $TESTS_LOG
     reports:
       junit:
         - $JUNITXML
   cache:
     key: ${output}-pytest
     policy: pull-push
     paths:
       - .pytest_cache
 
 # ------------------------------------------------------------------------------
 
 .protected_refs:
   extends: .manylinux
   rules:
     - if: '$CI_COMMIT_BRANCH == $CI_DEFAULT_BRANCH || $CI_COMMIT_TAG =~ /^v.*/'
 
 wheels:
   stage: build
   extends:
     - .protected_refs
     - .ccache
   before_script:
     - ccache --zero-stats || true
   script:
     - ./tests/ci/build_wheels.sh
   artifacts:
     paths:
       - dist/wheelhouse
 
 # ------------------------------------------------------------------------------
 
 .deploy_wheels:
   stage: deploy
   extends: .protected_refs
   dependencies:
     - wheels
   script:
     python -m twine upload --verbose dist/wheelhouse/*
 
 package:gitlab:
   extends: .deploy_wheels
   variables:
     TWINE_USERNAME: gitlab-ci-token
     TWINE_PASSWORD: ${CI_JOB_TOKEN}
     TWINE_REPOSITORY_URL: ${CI_API_V4_URL}/projects/${CI_PROJECT_ID}/packages/pypi
 
 package:pypi:
   extends: .deploy_wheels
   variables:
     TWINE_USERNAME: __token__
     TWINE_PASSWORD: ${PYPI_TOKEN}
   rules:
     - if: '$CI_COMMIT_TAG =~ /^v.*/'
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index 6e409f3..86c0018 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1,210 +1,211 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2021 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "statistics.hh"
 #include "fft_engine.hh"
 #include "loop.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <UInt dim>
-Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) { return std::sqrt(surface.var());
+Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) {
+   return std::sqrt(surface.var());
 }
 
 template <UInt dim>
 Real Statistics<dim>::computeSpectralRMSSlope(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(h_size);
   wavevectors *= 2 * M_PI;  // need q for slopes
   const auto psd = computePowerSpectrum(surface);
 
   const Real rms_slope_mean = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(VectorProxy<const Real, dim> q, const Complex& psd_val) {
         // Checking if we're in the zone that does not have hermitian symmetry
         if (std::abs(q.back()) < 1e-15)
           return q.l2squared() * psd_val.real();
         return 2 * q.l2squared() * psd_val.real();
       },
       range<VectorProxy<const Real, dim>>(wavevectors), psd);
 
   return std::sqrt(rms_slope_mean);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 GridHermitian<Real, dim>
 Statistics<dim>::computePowerSpectrum(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   GridHermitian<Real, dim> psd(h_size, surface.getNbComponents());
 
   FFTEngine::makeEngine()->forward(surface, psd);
 
   Real factor = 1. / surface.getGlobalNbPoints();
 
   // Squaring the fourier transform of surface and normalizing
   Loop::loop(
       [factor] CUDA_LAMBDA(Complex & c) {
         c *= factor;
         c *= conj(c);
       },
       psd);
 
   return psd;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Grid<Real, dim>
 Statistics<dim>::computeAutocorrelation(Grid<Real, dim>& surface) {
   Grid<Real, dim> acf(surface.sizes(), surface.getNbComponents());
   auto psd = computePowerSpectrum(surface);
 
   FFTEngine::makeEngine()->backward(acf, psd);
   acf *= acf.getGlobalNbPoints();
   return acf;
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 Real Statistics<dim>::contact(const Grid<Real, dim>& tractions,
                               UInt perimeter) {
   Real points = 0;
   UInt nc = tractions.getNbComponents();
 
   switch (nc) {
   case 1:
     points = Loop::reduce(
         [] CUDA_LAMBDA(const Real& t) -> Real { return t > 0; }, tractions);
     break;
   case 2:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 2> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 2>>(tractions));
     break;
   case 3:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 3> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 3>>(tractions));
     break;
   default:
     TAMAAS_EXCEPTION("Invalid number of components in traction");
   }
 
   auto area = points / tractions.getGlobalNbPoints();
 
   if (dim == 1)
     perimeter = 0;
 
   // Correction from Yastrebov et al. (Trib. Intl., 2017)
   // 10.1016/j.triboint.2017.04.023
   return area - (M_PI - 1 + std::log(2)) /
                     (24. * tractions.getGlobalNbPoints()) * perimeter;
 }
 
 /* -------------------------------------------------------------------------- */
 namespace {
 template <UInt dim>
 class moment_helper {
 public:
   moment_helper(const std::array<UInt, dim>& exp) : exponent(exp) {}
 
   CUDA_LAMBDA Complex operator()(VectorProxy<Real, dim> q,
                                  const Complex& phi) const {
     Real mul = 1;
     for (UInt i = 0; i < dim; ++i)
       mul *= std::pow(q(i), exponent[i]);
 
     // Do not duplicate everything from hermitian symmetry
     if (std::abs(q.back()) < 1e-15)
       return mul * phi;
     return 2 * mul * phi;
   }
 
 private:
   std::array<UInt, dim> exponent;
 };
 }  // namespace
 
 template <>
 std::vector<Real> Statistics<1>::computeMoments(Grid<Real, 1>& surface) {
   constexpr UInt dim = 1;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[1] = Loop::reduce<operation::plus>(moment_helper<dim>{{{2}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[2] = Loop::reduce<operation::plus>(moment_helper<dim>{{{4}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   return moments;
 }
 
 template <>
 std::vector<Real> Statistics<2>::computeMoments(Grid<Real, 2>& surface) {
   constexpr UInt dim = 2;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   auto m02 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m20 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   moments[1] = 0.5 * (m02 + m20);
   auto m22 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m40 = Loop::reduce<operation::plus>(moment_helper<dim>{{{4, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m04 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 4}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   moments[2] = (3 * m22 + m40 + m04) / 3.;
   return moments;
 }
 
 template struct Statistics<1>;
 template struct Statistics<2>;
 
 }  // namespace tamaas