diff --git a/CMakeLists.txt b/CMakeLists.txt
index f7631db2788..60d8384b5d2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -72,6 +72,7 @@ option(WITH_GSL "Build with GSL support" ON)
 option(WITH_CUDA   "Build with GPU support"  ON)
 option(WITH_HDF5   "Build with HDF5 support" ON)
 option(WITH_TESTS  "Enable tests"            ON)
+option(WITH_BENCHMARKS "Enable benchmarks"   OFF)
 option(WITH_SCAFACOS "Build with Scafacos support" ON)
 option(WITH_VALGRIND_INSTRUMENTATION "Build with valgrind instrumentation markers" OFF)
 if( CMAKE_VERSION VERSION_GREATER 3.5.2 AND CMAKE_CXX_COMPILER_ID STREQUAL "Clang" )
@@ -477,6 +478,11 @@ if(WITH_TESTS)
   add_subdirectory(testsuite)
 endif(WITH_TESTS)
 
+if(WITH_BENCHMARKS)
+  add_custom_target(benchmark)
+  add_subdirectory(maintainer/benchmarks)
+endif(WITH_BENCHMARKS)
+
 #######################################################################
 # Subdirectories
 #######################################################################
diff --git a/maintainer/benchmarks/CMakeLists.txt b/maintainer/benchmarks/CMakeLists.txt
new file mode 100644
index 00000000000..c438d02aaae
--- /dev/null
+++ b/maintainer/benchmarks/CMakeLists.txt
@@ -0,0 +1,94 @@
+if(NOT DEFINED TEST_NP)
+  include(ProcessorCount)
+  ProcessorCount(NP)
+  math(EXPR TEST_NP "${NP}/2 + 1")
+endif()
+
+if(EXISTS ${MPIEXEC})
+  # OpenMPI 3.0 and higher checks the number of processes against the number of CPUs
+  execute_process(COMMAND ${MPIEXEC} --version RESULT_VARIABLE mpi_version_result OUTPUT_VARIABLE mpi_version_output ERROR_VARIABLE mpi_version_output)
+  if (mpi_version_result EQUAL 0 AND mpi_version_output MATCHES "\\(Open(RTE| MPI)\\) ([3-9]\\.|1[0-9])")
+    set(MPIEXEC_OVERSUBSCRIBE "-oversubscribe")
+  else()
+    set(MPIEXEC_OVERSUBSCRIBE "")
+  endif()
+endif()
+
+function(PYTHON_BENCHMARK)
+  cmake_parse_arguments(BENCHMARK "" "FILE;RUN_WITH_MPI;MIN_NUM_PROC;MAX_NUM_PROC" "ARGUMENTS;DEPENDENCIES" ${ARGN})
+  get_filename_component(BENCHMARK_NAME ${BENCHMARK_FILE} NAME_WE)
+  foreach(argument IN LISTS BENCHMARK_ARGUMENTS)
+    string(REGEX REPLACE "[^-a-zA-Z0-9_\\.]+" "_" argument ${argument})
+    string(REGEX REPLACE "^[-_]+" "" argument ${argument})
+    set(BENCHMARK_NAME "${BENCHMARK_NAME}__${argument}")
+  endforeach(argument)
+  configure_file(${BENCHMARK_FILE} ${CMAKE_CURRENT_BINARY_DIR}/${BENCHMARK_FILE})
+  foreach(dependency IN LISTS BENCHMARK_DEPENDENCIES)
+    configure_file(${dependency} ${CMAKE_CURRENT_BINARY_DIR}/${dependency})
+  endforeach(dependency)
+  set(BENCHMARK_FILE "${CMAKE_CURRENT_BINARY_DIR}/${BENCHMARK_FILE}")
+  list(APPEND BENCHMARK_ARGUMENTS "--output=${CMAKE_BINARY_DIR}/benchmarks.csv.part")
+
+  # default values
+  if (NOT DEFINED BENCHMARK_RUN_WITH_MPI)
+    set(BENCHMARK_RUN_WITH_MPI TRUE)
+  endif()
+  if (NOT DEFINED BENCHMARK_MIN_NUM_PROC)
+    set(BENCHMARK_MIN_NUM_PROC 1)
+  endif()
+  if (NOT DEFINED BENCHMARK_MAX_NUM_PROC)
+    set(BENCHMARK_MAX_NUM_PROC ${NP})
+  endif()
+  # parallel schemes
+  if(EXISTS ${MPIEXEC} AND ${BENCHMARK_RUN_WITH_MPI})
+    set(BENCHMARK_CONFIGURATIONS "0")
+    if(${NP} GREATER 0 AND ${BENCHMARK_MAX_NUM_PROC} GREATER 0 AND ${BENCHMARK_MIN_NUM_PROC} LESS 2)
+      list(APPEND BENCHMARK_CONFIGURATIONS 1)
+    endif()
+    if(${NP} GREATER 1 AND ${BENCHMARK_MAX_NUM_PROC} GREATER 1 AND ${BENCHMARK_MIN_NUM_PROC} LESS 3)
+      list(APPEND BENCHMARK_CONFIGURATIONS 2)
+    endif()
+    if(${NP} GREATER 3 AND ${BENCHMARK_MAX_NUM_PROC} GREATER 3 AND ${BENCHMARK_MIN_NUM_PROC} LESS 5)
+      list(APPEND BENCHMARK_CONFIGURATIONS 4)
+    endif()
+    if(${NP} GREATER 7 AND ${BENCHMARK_MAX_NUM_PROC} GREATER 7 AND ${BENCHMARK_MIN_NUM_PROC} LESS 9)
+      list(APPEND BENCHMARK_CONFIGURATIONS 8)
+    endif()
+    if(${NP} GREATER 15 AND ${BENCHMARK_MAX_NUM_PROC} GREATER 15 AND ${BENCHMARK_MIN_NUM_PROC} LESS 17)
+      list(APPEND BENCHMARK_CONFIGURATIONS 16)
+    endif()
+    list(REMOVE_AT BENCHMARK_CONFIGURATIONS 0)
+    foreach(nproc IN LISTS BENCHMARK_CONFIGURATIONS)
+      add_test(NAME benchmark__${BENCHMARK_NAME}__parallel_${nproc}
+               COMMAND ${MPIEXEC} ${MPIEXEC_OVERSUBSCRIBE} ${MPIEXEC_NUMPROC_FLAG} ${nproc}
+                       ${CMAKE_BINARY_DIR}/pypresso ${BENCHMARK_FILE} ${BENCHMARK_ARGUMENTS}
+               CONFIGURATIONS "parallel")
+    endforeach(nproc)
+  else()
+      add_test(NAME benchmark__${BENCHMARK_NAME}__serial
+               COMMAND ${CMAKE_BINARY_DIR}/pypresso ${BENCHMARK_FILE} ${BENCHMARK_ARGUMENTS}
+               CONFIGURATIONS "serial")
+  endif()
+endfunction(PYTHON_BENCHMARK)
+
+python_benchmark(FILE lj.py ARGUMENTS "--particles_per_core=1000;--volume_fraction=0.50")
+python_benchmark(FILE lj.py ARGUMENTS "--particles_per_core=1000;--volume_fraction=0.02")
+python_benchmark(FILE lj.py ARGUMENTS "--particles_per_core=10000;--volume_fraction=0.50")
+python_benchmark(FILE lj.py ARGUMENTS "--particles_per_core=10000;--volume_fraction=0.02")
+python_benchmark(FILE p3m.py ARGUMENTS "--particles_per_core=1000;--volume_fraction=0.25;--bjerrum_length=4")
+python_benchmark(FILE p3m.py ARGUMENTS "--particles_per_core=10000;--volume_fraction=0.25;--bjerrum_length=4")
+
+add_custom_target(benchmark_python_serial COMMAND ${CMAKE_CTEST_COMMAND} $(ARGS) -C serial --output-on-failure)
+add_dependencies(benchmark_python_serial pypresso)
+
+add_custom_target(benchmark_python_parallel COMMAND ${CMAKE_CTEST_COMMAND} $(ARGS) -C parallel --output-on-failure)
+add_dependencies(benchmark_python_parallel pypresso)
+
+add_custom_target(benchmark_python)
+if(EXISTS ${MPIEXEC})
+  add_dependencies(benchmark_python pypresso benchmark_python_parallel)
+else()
+  add_dependencies(benchmark_python pypresso benchmark_python_serial)
+endif()
+
+add_dependencies(benchmark benchmark_python)
diff --git a/maintainer/benchmarks/lj.py b/maintainer/benchmarks/lj.py
new file mode 100644
index 00000000000..1b77fd0dea9
--- /dev/null
+++ b/maintainer/benchmarks/lj.py
@@ -0,0 +1,198 @@
+#
+# Copyright (C) 2013-2018 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo 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 General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+from __future__ import print_function
+import os
+import sys
+import numpy as np
+from time import time, sleep
+import argparse
+
+parser = argparse.ArgumentParser(description="Benchmark LJ simulations. "
+                                 "Save the results to a CSV file.")
+parser.add_argument("--particles_per_core", metavar="N", action="store",
+                    type=int, default=1000, required=False,
+                    help="Number of particles in the simulation box")
+parser.add_argument("--volume_fraction", metavar="FRAC", action="store",
+                    type=float, default=0.50, required=False,
+                    help="Fraction of the simulation box volume occupied by "
+                    "particles (range: [0.01-0.74], default: 0.50)")
+group = parser.add_mutually_exclusive_group()
+group.add_argument("--output", metavar="FILEPATH", action="store",
+                   type=str, required=False, default="benchmarks.csv",
+                   help="Output file (default: benchmarks.csv)")
+group.add_argument("--visualizer", action="store_true",
+                   help="Starts the visualizer (for debugging purposes)")
+
+args = parser.parse_args()
+
+# process and check arguments
+n_proc = int(os.environ.get("OMPI_COMM_WORLD_SIZE", 1))
+n_part = n_proc * args.particles_per_core
+measurement_steps = int(np.round(5e6 / args.particles_per_core, -2))
+assert args.volume_fraction > 0, "volume_fraction must be a positive number"
+assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \
+    "volume_fraction exceeds the physical limit of sphere packing (~0.74)"
+if not args.visualizer:
+    assert(measurement_steps >= 100), \
+        "{} steps per tick are too short".format(measurement_steps)
+
+
+import espressomd
+from espressomd import thermostat
+if args.visualizer:
+    from espressomd import visualization
+    from threading import Thread
+
+required_features = ["LENNARD_JONES"]
+espressomd.assert_features(required_features)
+
+print(espressomd.features())
+
+# Interaction parameters (Lennard-Jones)
+#############################################################
+
+lj_eps = 1.0  # LJ epsilon
+lj_sig = 1.0  # particle diameter
+lj_cut = lj_sig * 2**(1. / 6.)  # cutoff distance
+
+# System parameters
+#############################################################
+
+# volume of N spheres with radius r: N * (4/3*pi*r^3)
+box_l = (n_part * 4. / 3. * np.pi * (lj_sig / 2.)**3
+         / args.volume_fraction)**(1. / 3.)
+
+# System
+#############################################################
+system = espressomd.System(box_l=3 * (box_l,))
+# PRNG seeds
+#############################################################
+system.random_number_generator_state = list(range(
+    n_proc * (system._get_PRNG_state_size() + 1)))
+#np.random.seed(1)
+# Integration parameters
+#############################################################
+system.time_step = 0.01
+system.cell_system.skin = 0.5
+system.thermostat.turn_off()
+
+
+#############################################################
+#  Setup System                                             #
+#############################################################
+
+# Interaction setup
+#############################################################
+system.non_bonded_inter[0, 0].lennard_jones.set_params(
+    epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
+
+print("LJ-parameters:")
+print(system.non_bonded_inter[0, 0].lennard_jones.get_params())
+
+# Particle setup
+#############################################################
+
+for i in range(n_part):
+    system.part.add(id=i, pos=np.random.random(3) * system.box_l)
+
+#############################################################
+#  Warmup Integration                                       #
+#############################################################
+
+system.integrator.set_steepest_descent(
+    f_max=0,
+    gamma=0.001,
+    max_displacement=0.01)
+
+# warmup
+while system.analysis.energy()["total"] > 3 * n_part:
+    print("minimization: {:.1f}".format(system.analysis.energy()["total"]))
+    system.integrator.run(10)
+print()
+system.integrator.set_vv()
+
+system.thermostat.set_langevin(kT=1.0, gamma=1.0)
+
+# tune skin
+print("Tune skin: {}".format(system.cell_system.tune_skin(
+    min_skin=0.2, max_skin=1, tol=0.05, int_steps=100)))
+system.integrator.run(min(5 * measurement_steps, 60000))
+print("Tune skin: {}".format(system.cell_system.tune_skin(
+    min_skin=0.2, max_skin=1, tol=0.05, int_steps=100)))
+system.integrator.run(min(10 * measurement_steps, 60000))
+
+print(system.non_bonded_inter[0, 0].lennard_jones)
+
+if not args.visualizer:
+    # print initial energies
+    energies = system.analysis.energy()
+    print(energies)
+
+    # time integration loop
+    print("Timing every {} steps".format(measurement_steps))
+    main_tick = time()
+    all_t = []
+    for i in range(30):
+        tick = time()
+        system.integrator.run(measurement_steps)
+        tock = time()
+        t = (tock - tick) / measurement_steps
+        print("step {}, time = {:.2e}, verlet: {:.2f}"
+              .format(i, t, system.cell_system.get_state()["verlet_reuse"]))
+        all_t.append(t)
+    main_tock = time()
+    # average time
+    all_t = np.array(all_t)
+    avg = np.average(all_t)
+    ci = 1.96 * np.std(all_t) / np.sqrt(len(all_t) - 1)
+    print("average: {:.3e} +/- {:.3e} (95% C.I.)".format(avg, ci))
+
+    # print final energies
+    energies = system.analysis.energy()
+    print(energies)
+
+    # write report
+    cmd = " ".join(x for x in sys.argv[1:] if not x.startswith("--output"))
+    report = ('"{script}","{arguments}",{cores},"{mpi}",{mean:.3e},'
+              '{ci:.3e},{n},{dur:.1f},{E1:.5e},{E2:.5e},{E3:.5e}\n'.format(
+                  script=os.path.basename(sys.argv[0]), arguments=cmd,
+                  cores=n_proc, dur=main_tock - main_tick, n=measurement_steps,
+                  mpi="OMPI_COMM_WORLD_SIZE" in os.environ, mean=avg, ci=ci,
+                  E1=system.analysis.energy()["total"],
+                  E2=system.analysis.energy()["kinetic"],
+                  E3=system.analysis.energy()["non_bonded"]))
+    if not os.path.isfile(args.output):
+        report = ('"script","arguments","cores","MPI","mean","ci",'
+                  '"steps_per_tick","duration","E1","E2","E3"\n' + report)
+    with open(args.output, "a") as f:
+        f.write(report)
+else:
+    # use visualizer
+    visualizer = visualization.openGLLive(system)
+
+    def main_thread():
+        while True:
+            system.integrator.run(1)
+            visualizer.update()
+            sleep(1 / 60.)  # limit framerate to at most 60 FPS
+
+    t = Thread(target=main_thread)
+    t.daemon = True
+    t.start()
+    visualizer.start()
diff --git a/maintainer/benchmarks/p3m.py b/maintainer/benchmarks/p3m.py
new file mode 100644
index 00000000000..819648f901e
--- /dev/null
+++ b/maintainer/benchmarks/p3m.py
@@ -0,0 +1,209 @@
+#
+# Copyright (C) 2013-2018 The ESPResSo project
+#
+# This file is part of ESPResSo.
+#
+# ESPResSo is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ESPResSo 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 General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+from __future__ import print_function
+import os
+import sys
+import numpy as np
+from time import time
+import argparse
+
+parser = argparse.ArgumentParser(description="Benchmark P3M simulations. "
+                                 "Save the results to a CSV file.")
+parser.add_argument("--particles_per_core", metavar="N", action="store",
+                    type=int, default=1000, required=False,
+                    help="Number of particles in the simulation box")
+parser.add_argument("--volume_fraction", metavar="FRAC", action="store",
+                    type=float, default=0.25, required=False,
+                    help="Fraction of the simulation box volume occupied by "
+                    "particles (range: [0.01-0.74], default: 0.25)")
+parser.add_argument("--bjerrum_length", metavar="LENGTH", action="store",
+                    type=float, default=4., required=False,
+                    help="Bjerrum length (default: 4)")
+group = parser.add_mutually_exclusive_group()
+group.add_argument("--output", metavar="FILEPATH", action="store",
+                   type=str, required=False, default="benchmarks.csv",
+                   help="Output file (default: benchmarks.csv)")
+group.add_argument("--visualizer", action="store_true",
+                   help="Starts the visualizer (for debugging purposes)")
+
+args = parser.parse_args()
+
+# process and check arguments
+n_proc = int(os.environ.get("OMPI_COMM_WORLD_SIZE", 1))
+n_part = n_proc * args.particles_per_core
+measurement_steps = int(np.round(5e5 / args.particles_per_core, -1))
+assert args.bjerrum_length > 0, "bjerrum_length must be a positive number"
+assert args.volume_fraction > 0, "volume_fraction must be a positive number"
+assert args.volume_fraction < np.pi / (3 * np.sqrt(2)), \
+    "volume_fraction exceeds the physical limit of sphere packing (~0.74)"
+if not args.visualizer:
+    assert(measurement_steps >= 50), \
+        "{} steps per tick are too short".format(measurement_steps)
+
+
+import espressomd
+from espressomd import thermostat
+from espressomd import electrostatics
+if args.visualizer:
+    from espressomd import visualization
+    from threading import Thread
+
+required_features = ["ELECTROSTATICS", "LENNARD_JONES", "MASS"]
+espressomd.assert_features(required_features)
+
+print(espressomd.features())
+
+# Interaction parameters (Lennard-Jones, Coulomb)
+#############################################################
+
+species = ["anion", "cation"]
+types = {"anion": 0, "cation": 0}
+charges = {"anion": -1.0, "cation": 1.0}
+lj_sigmas = {"anion": 1.0, "cation": 1.0}
+lj_epsilons = {"anion": 1.0, "cation": 1.0}
+WCA_cut = 2.**(1. / 6.)
+lj_cuts = {"anion": WCA_cut * lj_sigmas["anion"],
+           "cation": WCA_cut * lj_sigmas["cation"]}
+masses = {"anion": 1.0, "cation": 1.0}
+
+# System parameters
+#############################################################
+
+# volume of N spheres with radius r: N * (4/3*pi*r^3)
+lj_sig = (lj_sigmas["cation"] + lj_sigmas["anion"]) / 2
+box_l = (n_part * 4. / 3. * np.pi * (lj_sig / 2.)**3
+         / args.volume_fraction)**(1. / 3.)
+
+# System
+#############################################################
+system = espressomd.System(box_l=3 * (box_l,))
+system.cell_system.set_domain_decomposition(use_verlet_lists=True)
+# PRNG seeds
+#############################################################
+system.random_number_generator_state = list(range(
+    n_proc * (system._get_PRNG_state_size() + 1)))
+# Integration parameters
+#############################################################
+system.time_step = 0.01
+system.cell_system.skin = .4 
+system.thermostat.turn_off()
+
+
+#############################################################
+#  Setup System                                             #
+#############################################################
+
+# Interaction setup
+#############################################################
+
+for i in range(len(species)):
+    ion1 = species[i]
+    for j in range(i, len(species)):
+        ion2 = species[j]
+        lj_sig = (lj_sigmas[ion1] + lj_sigmas[ion2]) / 2
+        lj_cut = (lj_cuts[ion1] + lj_cuts[ion2]) / 2
+        lj_eps = (lj_epsilons[ion1] * lj_epsilons[ion2])**(1. / 2.)
+        system.non_bonded_inter[types[ion1],
+                                types[ion2]].lennard_jones.set_params(
+            epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
+
+# Particle setup
+#############################################################
+
+for i in range(0, n_part, len(species)):
+    for t in species:
+        system.part.add(pos=np.random.random(3) * system.box_l,
+                        q=charges[t], type=types[t], mass=masses[t])
+
+#############################################################
+#  Warmup Integration                                       #
+#############################################################
+
+energy = system.analysis.energy()
+print("Before Minimization: E_total = {}".format(energy["total"]))
+system.minimize_energy.init(f_max=1000, gamma=30.0,
+                            max_steps=1000, max_displacement=0.05)
+system.minimize_energy.minimize()
+system.minimize_energy.minimize()
+energy = system.analysis.energy()
+print("After Minimization: E_total = {}".format(energy["total"]))
+
+
+system.integrator.set_vv()
+system.thermostat.set_langevin(kT=1.0, gamma=1.0)
+
+system.integrator.run(min(3 * measurement_steps, 1000))
+print("Tune skin: {}".format(system.cell_system.tune_skin(
+    min_skin=0.4, max_skin=1.6, tol=0.05, int_steps=100)))
+system.integrator.run(min(3 * measurement_steps, 3000))
+print("Tune p3m")
+p3m = electrostatics.P3M(prefactor=args.bjerrum_length, accuracy=1e-4)
+system.actors.add(p3m)
+system.integrator.run(min(3 * measurement_steps, 3000))
+print("Tune skin: {}".format(system.cell_system.tune_skin(
+    min_skin=1.0, max_skin=1.6, tol=0.05, int_steps=100)))
+
+
+if not args.visualizer:
+    # print initial energies
+    energies = system.analysis.energy()
+    print(energies)
+
+    # time integration loop
+    print("Timing every {} steps".format(measurement_steps))
+    main_tick = time()
+    all_t = []
+    for i in range(30):
+        tick = time()
+        system.integrator.run(measurement_steps)
+        tock = time()
+        t = (tock - tick) / measurement_steps
+        print("step {}, time = {:.2e}, verlet: {:.2f}"
+              .format(i, t, system.cell_system.get_state()["verlet_reuse"]))
+        all_t.append(t)
+    main_tock = time()
+    # average time
+    all_t = np.array(all_t)
+    avg = np.average(all_t)
+    ci = 1.96 * np.std(all_t) / np.sqrt(len(all_t) - 1)
+    print("average: {:.3e} +/- {:.3e} (95% C.I.)".format(avg, ci))
+
+    # print final energies
+    energies = system.analysis.energy()
+    print(energies)
+
+    # write report
+    cmd = " ".join(x for x in sys.argv[1:] if not x.startswith("--output"))
+    report = ('"{script}","{arguments}",{cores},"{mpi}",{mean:.3e},'
+              '{ci:.3e},{n},{dur:.1f},{E1:.5e},{E2:.5e},{E3:.5e}\n'.format(
+                  script=os.path.basename(sys.argv[0]), arguments=cmd,
+                  cores=n_proc, dur=main_tock - main_tick, n=measurement_steps,
+                  mpi="OMPI_COMM_WORLD_SIZE" in os.environ, mean=avg, ci=ci,
+                  E1=system.analysis.energy()["total"],
+                  E2=system.analysis.energy()["coulomb"],
+                  E3=system.analysis.energy()["non_bonded"]))
+    if not os.path.isfile(args.output):
+        report = ('"script","arguments","cores","MPI","mean","ci",'
+                  '"steps_per_tick","duration","E1","E2","E3"\n' + report)
+    with open(args.output, "a") as f:
+        f.write(report)
+else:
+    # use visualizer
+    visualizer = visualization.openGLLive(system)
+    visualizer.run(1)
diff --git a/maintainer/benchmarks/runner.sh b/maintainer/benchmarks/runner.sh
new file mode 100644
index 00000000000..245b10ca9b5
--- /dev/null
+++ b/maintainer/benchmarks/runner.sh
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+cd "$(git rev-parse --show-toplevel)"
+mkdir -p build
+cd build
+
+# manage headers files with different features
+configs="myconfig-minimal.hpp myconfig-default.hpp myconfig-maxset.hpp"
+cat > myconfig-minimal.hpp << EOF
+#define ELECTROSTATICS
+#define LENNARD_JONES
+#define MASS
+EOF
+cp ../src/core/myconfig-default.hpp myconfig-default.hpp
+sed 's/#define ADDITIONAL_CHECKS//' ../maintainer/configs/maxset.hpp > myconfig-maxset.hpp
+
+# prepare build area
+rm -rf src/ maintainer/
+cmake -DWITH_BENCHMARKS=ON ..
+cat > benchmarks.csv << EOF
+"config","script","arguments","cores","MPI","mean","ci","steps_per_tick","duration","E1","E2","E3"
+EOF
+
+# run benchmarks
+for config in ${configs}
+do
+  echo "### ${config}" >> benchmarks.log
+  cp ${config} myconfig.hpp
+  make -j$(nproc)
+  rm -f benchmarks.csv.part
+  touch benchmarks.csv.part
+  make benchmark 2>&1 | tee -a benchmarks.log
+  sed -ri "s/^/\"$(basename ${config})\",/" benchmarks.csv.part
+  cat benchmarks.csv.part >> benchmarks.csv
+done
+
+rm benchmarks.csv.part
+
diff --git a/maintainer/benchmarks/suite.sh b/maintainer/benchmarks/suite.sh
new file mode 100644
index 00000000000..0fd59db15f5
--- /dev/null
+++ b/maintainer/benchmarks/suite.sh
@@ -0,0 +1,30 @@
+#!/bin/bash
+
+# list of commits to benchmark
+commits="HEAD"
+
+cd "$(git rev-parse --show-toplevel)"
+mkdir -p build
+cd build
+
+# prepare output files
+rm -f benchmarks.log
+cat > benchmarks_suite.csv << EOF
+"commit","config","script","arguments","cores","MPI","mean","ci","steps_per_tick","duration","E1","E2","E3"
+EOF
+
+# run benchmarks
+for commit in ${commits}
+do
+  echo "### commit ${commit}" >> benchmarks.log
+  git checkout ${commit} ../src ../libs
+  bash ../maintainer/benchmarks/runner.sh
+  sed -ri "s/^/\"${commit}\",/" benchmarks.csv
+  tail -n +2 benchmarks.csv >> benchmarks_suite.csv
+done
+
+rm benchmarks.csv
+
+# restore files
+git checkout HEAD ../src ../libs
+