diff --git a/doc/sphinx/electrostatics.rst b/doc/sphinx/electrostatics.rst
index e31caff4bd7..df447c6d5d8 100644
--- a/doc/sphinx/electrostatics.rst
+++ b/doc/sphinx/electrostatics.rst
@@ -386,6 +386,4 @@ for an accuracy of :math:`10^{-3}`::
 
 For details of the various methods and their parameters please refer to
 the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a
-shared library. ScaFaCoS can be used only once, either for Coulomb or for
-dipolar interactions.
-
+shared library.
diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst
index f8c7393c8c3..20c197557c9 100644
--- a/doc/sphinx/installation.rst
+++ b/doc/sphinx/installation.rst
@@ -303,7 +303,7 @@ General features
        :ref:`Magnetostatics / Dipolar interactions`
        :ref:`Electrostatics`
 
--  ``SCAFACOS_DIPOLES``
+-  ``SCAFACOS_DIPOLES`` This activates magnetostatics methods of ScaFaCoS.
 
 -  ``ROTATION`` Switch on rotational degrees of freedom for the particles, as well as
    the corresponding quaternion integrator.
diff --git a/doc/sphinx/magnetostatics.rst b/doc/sphinx/magnetostatics.rst
index 52fc41df91d..21fb78f8194 100644
--- a/doc/sphinx/magnetostatics.rst
+++ b/doc/sphinx/magnetostatics.rst
@@ -211,6 +211,4 @@ root mean square accuracy for the magnetic field.
 
 For details of the various methods and their parameters please refer to
 the ScaFaCoS manual. To use this feature, ScaFaCoS has to be built as a
-shared library. ScaFaCoS can be used only once, either for Coulomb or for
-dipolar interactions.
-
+shared library.
diff --git a/src/core/electrostatics_magnetostatics/CMakeLists.txt b/src/core/electrostatics_magnetostatics/CMakeLists.txt
index 498dcf79e30..fcd3037dfbf 100644
--- a/src/core/electrostatics_magnetostatics/CMakeLists.txt
+++ b/src/core/electrostatics_magnetostatics/CMakeLists.txt
@@ -14,6 +14,7 @@ target_sources(
           ${CMAKE_CURRENT_SOURCE_DIR}/p3m-dipolar.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/p3m_gpu.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/scafacos.cpp
+          ${CMAKE_CURRENT_SOURCE_DIR}/ScafacosContext.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/fft.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/coulomb.cpp
           ${CMAKE_CURRENT_SOURCE_DIR}/dipole.cpp
diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.cpp b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp
new file mode 100644
index 00000000000..eca126e5d94
--- /dev/null
+++ b/src/core/electrostatics_magnetostatics/ScafacosContext.cpp
@@ -0,0 +1,247 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+/** @file
+ *  Provide a C-like interface for ScaFaCoS.
+ */
+
+#include "config.hpp"
+
+#if defined(SCAFACOS)
+
+#include "electrostatics_magnetostatics/ScafacosContext.hpp"
+
+#include "cells.hpp"
+#include "communication.hpp"
+#include "electrostatics_magnetostatics/coulomb.hpp"
+#include "electrostatics_magnetostatics/dipole.hpp"
+#include "electrostatics_magnetostatics/scafacos.hpp"
+#include "grid.hpp"
+#include "integrate.hpp"
+#include "particle_data.hpp"
+#include "tuning.hpp"
+
+#include <utils/Vector.hpp>
+
+#include <boost/mpi/collectives.hpp>
+#include <boost/range/algorithm/min_element.hpp>
+#include <boost/range/numeric.hpp>
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <functional>
+#include <limits>
+#include <list>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+namespace Scafacos {
+
+std::string ScafacosContext::get_method_and_parameters() {
+  std::string representation = get_method() + " " + get_parameters();
+  std::replace(representation.begin(), representation.end(), ',', ' ');
+  return representation;
+}
+
+void ScafacosContext::update_system_params() {
+  int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0,
+                box_geo.periodic(2) != 0};
+
+  auto const n_part = boost::mpi::all_reduce(
+      comm_cart, cell_structure.local_particles().size(), std::plus<>());
+  set_common_parameters(box_geo.length().data(), per, n_part);
+}
+
+void ScafacosContextCoulomb::update_particle_data() {
+  positions.clear();
+  charges.clear();
+
+  for (auto const &p : cell_structure.local_particles()) {
+    auto const pos = folded_position(p.r.p, box_geo);
+    positions.push_back(pos[0]);
+    positions.push_back(pos[1]);
+    positions.push_back(pos[2]);
+    charges.push_back(p.p.q);
+  }
+}
+
+#ifdef SCAFACOS_DIPOLES
+void ScafacosContextDipoles::update_particle_data() {
+  positions.clear();
+  dipoles.clear();
+
+  for (auto const &p : cell_structure.local_particles()) {
+    auto const pos = folded_position(p.r.p, box_geo);
+    positions.push_back(pos[0]);
+    positions.push_back(pos[1]);
+    positions.push_back(pos[2]);
+    auto const dip = p.calc_dip();
+    dipoles.push_back(dip[0]);
+    dipoles.push_back(dip[1]);
+    dipoles.push_back(dip[2]);
+  }
+}
+#endif
+
+void ScafacosContextCoulomb::update_particle_forces() const {
+  if (positions.empty())
+    return;
+
+  int it = 0;
+  for (auto &p : cell_structure.local_particles()) {
+    p.f.f += coulomb.prefactor * p.p.q *
+             Utils::Vector3d(Utils::Span<const double>(&(fields[it]), 3));
+    it += 3;
+  }
+
+  /* Check that the particle number did not change */
+  assert(it == fields.size());
+}
+
+#ifdef SCAFACOS_DIPOLES
+void ScafacosContextDipoles::update_particle_forces() const {
+  if (positions.empty())
+    return;
+
+  int it = 0;
+  for (auto &p : cell_structure.local_particles()) {
+    // Indices
+    // 3 "potential" values per particles (see below)
+    int const it_p = 3 * it;
+    // 6 "field" values per particles (see below)
+    int const it_f = 6 * it;
+
+    // The scafacos term "potential" here in fact refers to the magnetic
+    // field. So, the torques are given by m \times B
+    auto const dip = p.calc_dip();
+    auto const t = vector_product(
+        dip,
+        Utils::Vector3d(Utils::Span<const double>(&(potentials[it_p]), 3)));
+    // The force is given by G m, where G is a matrix
+    // which comes from the "fields" output of scafacos like this
+    // 0 1 2
+    // 1 3 4
+    // 2 4 5
+    // where the numbers refer to indices in the "field" output from scafacos
+    auto const G = Utils::Vector<Utils::Vector3d, 3>{
+        {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]},
+        {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]},
+        {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}};
+    auto const f = G * dip;
+
+    // Add to particles
+    p.f.f += dipole.prefactor * f;
+    p.f.torque += dipole.prefactor * t;
+    it++;
+  }
+
+  /* Check that the particle number did not change */
+  assert(it == positions.size() / 3);
+}
+#endif
+
+double ScafacosContextCoulomb::long_range_energy() {
+  update_particle_data();
+  run(charges, positions, fields, potentials);
+  return 0.5 * coulomb.prefactor *
+         boost::inner_product(charges, potentials, 0.0);
+}
+
+#ifdef SCAFACOS_DIPOLES
+double ScafacosContextDipoles::long_range_energy() {
+  update_particle_data();
+  run_dipolar(dipoles, positions, fields, potentials);
+  return -0.5 * dipole.prefactor *
+         boost::inner_product(dipoles, potentials, 0.0);
+}
+#endif
+
+static void set_r_cut_and_tune_local_worker(double r_cut) {
+  set_r_cut_and_tune_local(r_cut);
+}
+
+REGISTER_CALLBACK(set_r_cut_and_tune_local_worker)
+
+/** Determine runtime for a specific cutoff */
+static double time_r_cut(double r_cut) {
+  /* Set cutoff to time */
+  mpi_call_all(set_r_cut_and_tune_local_worker, r_cut);
+
+  return time_force_calc(10);
+}
+
+/** Determine the optimal cutoff by bisection */
+static void tune_r_cut() {
+  const double tune_limit = 1e-3;
+
+  auto const min_box_l = *boost::min_element(box_geo.length());
+  auto const min_local_box_l = *boost::min_element(local_geo.length());
+
+  /* scafacos p3m and Ewald do not accept r_cut 0 for no good reason */
+  double r_min = 1.0;
+  double r_max = std::min(min_local_box_l, min_box_l / 2.0) - skin;
+  double t_min = 0;
+  double t_max = std::numeric_limits<double>::max();
+
+  /* Run bisection */
+  while (std::fabs(r_min - r_max) > tune_limit) {
+    const double dr = 0.5 * (r_max - r_min);
+    const double t_mid = time_r_cut(r_min + dr);
+    t_min = time_r_cut(r_min);
+    t_max = time_r_cut(r_max);
+
+    if (t_min <= 0.0 or t_max <= 0.0) {
+      break;
+    }
+
+    if (t_mid > t_min) {
+      r_max = r_min += dr;
+    } else {
+      r_min += dr;
+    }
+  }
+}
+
+void ScafacosContextCoulomb::tune() {
+  update_particle_data();
+
+  // Check whether we have to do a bisection for the short range cutoff
+  // Check if there is a user-supplied cutoff
+  if (has_near && r_cut() <= 0.0) {
+    // run tuning on the head node
+    if (this_node == 0) {
+      tune_r_cut();
+    }
+  } else {
+    // ESPResSo is not affected by a short range cutoff -> tune in parallel
+    Scafacos::tune(charges, positions);
+  }
+}
+
+void ScafacosContextCoulomb::set_r_cut_and_tune_local(double r_cut) {
+  update_particle_data();
+  set_r_cut(r_cut);
+  Scafacos::tune(charges, positions);
+}
+
+} // namespace Scafacos
+#endif /* SCAFACOS */
diff --git a/src/core/electrostatics_magnetostatics/ScafacosContext.hpp b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp
new file mode 100644
index 00000000000..e22d1039b4a
--- /dev/null
+++ b/src/core/electrostatics_magnetostatics/ScafacosContext.hpp
@@ -0,0 +1,107 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+#ifndef SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP
+#define SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP
+/**
+ * @file
+ * @ref Scafacos::ScafacosContextBase implements the interface of the
+ * ScaFaCoS bridge. It is further derived for the coulombic and dipolar
+ * versions of ScaFaCoS.
+ */
+
+#include "config.hpp"
+
+#if defined(SCAFACOS)
+
+#include "Scafacos.hpp"
+#include "electrostatics_magnetostatics/ScafacosContextBase.hpp"
+
+#include <utils/Vector.hpp>
+
+#include <string>
+#include <vector>
+
+#if defined(SCAFACOS_DIPOLES) && !defined(FCS_ENABLE_DIPOLES)
+#error                                                                         \
+    "SCAFACOS_DIPOLES requires dipoles support in scafacos library (FCS_ENABLE_DIPOLES)."
+#endif
+
+namespace Scafacos {
+
+/** Encapsulation for the particle data needed by ScaFaCoS */
+struct ScafacosContext : ScafacosContextBase, Scafacos {
+  using Scafacos::Scafacos;
+  using ScafacosContextBase::ScafacosContextBase;
+  ~ScafacosContext() override = default;
+  void update_system_params() override;
+  double get_pair_force(double dist) const override {
+    return Scafacos::pair_force(dist);
+  }
+  double get_pair_energy(double dist) const override {
+    return Scafacos::pair_energy(dist);
+  }
+  std::string get_method_and_parameters() override;
+  double get_r_cut() const override { return Scafacos::r_cut(); }
+
+protected:
+  /** Outputs */
+  std::vector<double> fields, potentials;
+};
+
+struct ScafacosContextCoulomb : ScafacosContext {
+  using ScafacosContext::ScafacosContext;
+  void update_particle_data() override;
+  void update_particle_forces() const override;
+  double long_range_energy() override;
+  void add_long_range_force() override {
+    update_particle_data();
+    run(charges, positions, fields, potentials);
+    update_particle_forces();
+  }
+  void tune();
+  void set_r_cut_and_tune_local(double r_cut);
+
+private:
+  /** Inputs */
+  std::vector<double> positions, charges;
+};
+
+#ifdef SCAFACOS_DIPOLES
+struct ScafacosContextDipoles : ScafacosContext {
+  using ScafacosContext::ScafacosContext;
+  void update_particle_data() override;
+  void update_particle_forces() const override;
+  double long_range_energy() override;
+  void add_long_range_force() override {
+    update_particle_data();
+    run_dipolar(dipoles, positions, fields, potentials);
+    update_particle_forces();
+  }
+
+private:
+  /** Inputs */
+  std::vector<double> positions, dipoles;
+};
+#endif
+
+} // namespace Scafacos
+#endif // SCAFACOS
+#endif // SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXT_HPP
diff --git a/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp b/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp
new file mode 100644
index 00000000000..b3f6b32a219
--- /dev/null
+++ b/src/core/electrostatics_magnetostatics/ScafacosContextBase.hpp
@@ -0,0 +1,84 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * 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/>.
+ */
+#ifndef SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXTBASE_HPP
+#define SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXTBASE_HPP
+/**
+ * @file
+ * @ref Scafacos::ScafacosContextBase provides the public interface
+ * of the ScaFaCoS bridge. It relies on type erasure to hide the
+ * ScaFaCoS implementation details from the ESPResSo core. It is
+ * implemented by @ref Scafacos::ScafacosContext.
+ */
+
+#include "config.hpp"
+
+#if defined(SCAFACOS)
+
+#include <utils/Vector.hpp>
+
+#include <string>
+
+namespace Scafacos {
+
+/**
+ * @brief Public interface to the ScaFaCoS context.
+ * Implementations of this class will store particle positions and
+ * charges/dipoles in flat arrays, as required by ScaFaCoS, as well
+ * as the output arrays.
+ */
+struct ScafacosContextBase {
+  virtual ~ScafacosContextBase() = default;
+  /** @brief Collect particle data in continuous arrays. */
+  virtual void update_particle_data() = 0;
+  /** @brief Write forces back to particles. */
+  virtual void update_particle_forces() const = 0;
+  /** @brief Calculate long-range part of the energy. */
+  virtual double long_range_energy() = 0;
+  /** @brief Add long-range part of the forces to particles. */
+  virtual void add_long_range_force() = 0;
+  /** @brief Add near-field pair force. */
+  inline void add_pair_force(double q1q2, Utils::Vector3d const &d, double dist,
+                             Utils::Vector3d &force) {
+    if (dist > get_r_cut())
+      return;
+
+    auto const field = get_pair_force(dist);
+    auto const fak = q1q2 * field / dist;
+    force -= fak * d;
+  }
+  /** @brief Calculate near-field pair energy. */
+  inline double pair_energy(double q1q2, double dist) {
+    if (dist > get_r_cut())
+      return 0.;
+
+    return q1q2 * get_pair_energy(dist);
+  }
+  /** @brief Reinitialize number of particles, box shape and periodicity. */
+  virtual void update_system_params() = 0;
+  virtual double get_r_cut() const = 0;
+  virtual double get_pair_force(double dist) const = 0;
+  virtual double get_pair_energy(double dist) const = 0;
+  virtual std::string get_method_and_parameters() = 0;
+};
+
+} // namespace Scafacos
+#endif // SCAFACOS
+#endif // SRC_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOSCONTEXTBASE_HPP
diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp
index af3a2720a37..1fdf502f41d 100644
--- a/src/core/electrostatics_magnetostatics/coulomb.cpp
+++ b/src/core/electrostatics_magnetostatics/coulomb.cpp
@@ -122,7 +122,7 @@ double cutoff(const Utils::Vector3d &box_l) {
     return rf_params.r_cut;
 #ifdef SCAFACOS
   case COULOMB_SCAFACOS:
-    return Scafacos::get_r_cut();
+    return Scafacos::fcs_coulomb()->get_r_cut();
 #endif
   default:
     return -1.0;
@@ -217,7 +217,7 @@ void on_boxl_change() {
     break;
 #ifdef SCAFACOS
   case COULOMB_SCAFACOS:
-    Scafacos::update_system_params();
+    Scafacos::fcs_coulomb()->update_system_params();
     break;
 #endif
   default:
@@ -291,8 +291,7 @@ void calc_long_range_force(const ParticleRange &particles) {
 #endif
 #ifdef SCAFACOS
   case COULOMB_SCAFACOS:
-    assert(!Scafacos::dipolar());
-    Scafacos::add_long_range_force();
+    Scafacos::fcs_coulomb()->add_long_range_force();
     break;
 #endif
   default:
@@ -350,8 +349,7 @@ double calc_energy_long_range(const ParticleRange &particles) {
 #endif
 #ifdef SCAFACOS
   case COULOMB_SCAFACOS:
-    assert(!Scafacos::dipolar());
-    energy += Scafacos::long_range_energy();
+    energy += Scafacos::fcs_coulomb()->long_range_energy();
     break;
 #endif
   default:
diff --git a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp
index a9b58edad4b..d291eb7f095 100644
--- a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp
+++ b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp
@@ -58,7 +58,7 @@ inline Utils::Vector3d central_force(double const q1q2,
     break;
 #ifdef SCAFACOS
   case COULOMB_SCAFACOS:
-    Scafacos::add_pair_force(q1q2, d.data(), dist, f.data());
+    Scafacos::fcs_coulomb()->add_pair_force(q1q2, d, dist, f);
     break;
 #endif
   default:
@@ -157,7 +157,7 @@ inline double pair_energy(Particle const &p1, Particle const &p2,
 #endif
 #ifdef SCAFACOS
     case COULOMB_SCAFACOS:
-      return Scafacos::pair_energy(q1q2, dist);
+      return Scafacos::fcs_coulomb()->pair_energy(q1q2, dist);
 #endif
     case COULOMB_DH:
       return dh_coulomb_pair_energy(q1q2, dist);
diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp
index f188320be3a..15716439d6a 100644
--- a/src/core/electrostatics_magnetostatics/dipole.cpp
+++ b/src/core/electrostatics_magnetostatics/dipole.cpp
@@ -138,9 +138,9 @@ void on_boxl_change() {
     dp3m_scaleby_box_l();
     break;
 #endif
-#ifdef SCAFACOS
+#ifdef SCAFACOS_DIPOLES
   case DIPOLAR_SCAFACOS:
-    Scafacos::update_system_params();
+    Scafacos::fcs_dipoles()->update_system_params();
     break;
 #endif
   default:
@@ -202,8 +202,7 @@ void calc_long_range_force(const ParticleRange &particles) {
 #endif
 #ifdef SCAFACOS_DIPOLES
   case DIPOLAR_SCAFACOS:
-    assert(Scafacos::dipolar());
-    Scafacos::add_long_range_force();
+    Scafacos::fcs_dipoles()->add_long_range_force();
 #endif
   case DIPOLAR_NONE:
     break;
@@ -249,8 +248,8 @@ double calc_energy_long_range(const ParticleRange &particles) {
 #endif
 #ifdef SCAFACOS_DIPOLES
   case DIPOLAR_SCAFACOS:
-    assert(Scafacos::dipolar());
-    energy = Scafacos::long_range_energy();
+    energy = Scafacos::fcs_dipoles()->long_range_energy();
+    break;
 #endif
   case DIPOLAR_NONE:
     break;
diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp
index 2c8fc72da1d..81c09e77823 100644
--- a/src/core/electrostatics_magnetostatics/scafacos.cpp
+++ b/src/core/electrostatics_magnetostatics/scafacos.cpp
@@ -18,9 +18,6 @@
  * You should have received a copy of the GNU General Public License
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  */
-/** @file
- *  Provide a C-like interface for ScaFaCoS.
- */
 
 #include "config.hpp"
 
@@ -29,373 +26,137 @@
 #include "electrostatics_magnetostatics/scafacos.hpp"
 
 #include "Scafacos.hpp"
-#include "cells.hpp"
 #include "communication.hpp"
+#include "electrostatics_magnetostatics/ScafacosContext.hpp"
 #include "electrostatics_magnetostatics/coulomb.hpp"
 #include "electrostatics_magnetostatics/dipole.hpp"
 #include "errorhandling.hpp"
 #include "event.hpp"
 #include "grid.hpp"
-#include "integrate.hpp"
-#include "particle_data.hpp"
-#include "tuning.hpp"
 
 #include <utils/Vector.hpp>
 
-#include <boost/mpi/collectives.hpp>
-#include <boost/range/algorithm/min_element.hpp>
-
-#include <algorithm>
-#include <cassert>
-#include <cmath>
-#include <functional>
-#include <limits>
 #include <list>
-#include <numeric>
 #include <stdexcept>
 #include <string>
-#include <vector>
-
-#if defined(SCAFACOS_DIPOLES) && !defined(FCS_ENABLE_DIPOLES)
-#error                                                                         \
-    "SCAFACOS_DIPOLES requires dipoles support in scafacos library (FCS_ENABLE_DIPOLES)."
-#endif
 
 namespace Scafacos {
 
-/** Get available scafacos methods */
+/** Get available ScaFaCoS methods */
 std::list<std::string> available_methods() {
   return Scafacos::available_methods();
 }
 
-/** Encapsulation for the particle data needed by scafacos */
-struct ScafacosData {
-  int update_particle_data();
-  void update_particle_forces() const;
-  /** Inputs */
-  std::vector<double> charges, positions, dipoles;
-  /** Outputs */
-  std::vector<double> fields, potentials;
-};
-
-static Scafacos *scafacos = nullptr;
-static ScafacosData particles;
-
-/** \brief Collect particle data in continuous arrays as required by fcs */
-int ScafacosData::update_particle_data() {
-  positions.clear();
-
-  if (!dipolar()) {
-    charges.clear();
-  } else {
+namespace {
 #ifdef SCAFACOS_DIPOLES
-    dipoles.clear();
+ScafacosContextDipoles *dipoles_instance = nullptr;
 #endif
-  }
+ScafacosContextCoulomb *coulomb_instance = nullptr;
+} // namespace
 
-  for (auto const &p : cell_structure.local_particles()) {
-    auto const pos = folded_position(p.r.p, box_geo);
-    positions.push_back(pos[0]);
-    positions.push_back(pos[1]);
-    positions.push_back(pos[2]);
-    if (!dipolar()) {
-      charges.push_back(p.p.q);
-    } else {
 #ifdef SCAFACOS_DIPOLES
-      auto const dip = p.calc_dip();
-      dipoles.push_back(dip[0]);
-      dipoles.push_back(dip[1]);
-      dipoles.push_back(dip[2]);
-#endif
-    }
+ScafacosContextBase *fcs_dipoles() {
+  if (!dipoles_instance) {
+    throw std::runtime_error(
+        "Attempted access to uninitialized Scafacos Dipoles instance.");
   }
-
-  return static_cast<int>(positions.size() / 3);
+  return dipoles_instance;
 }
-
-/** \brief Write forces back to particles */
-void ScafacosData::update_particle_forces() const {
-  int it = 0;
-  if (positions.empty())
-    return;
-
-  for (auto &p : cell_structure.local_particles()) {
-    if (!dipolar()) {
-      p.f.f += coulomb.prefactor * p.p.q *
-               Utils::Vector3d(Utils::Span<const double>(&(fields[it]), 3));
-      it += 3;
-    } else {
-#ifdef SCAFACOS_DIPOLES
-      // Indices
-      // 3 "potential" values per particles (see below)
-      int const it_p = 3 * it;
-      // 6 "field" values per particles (see below)
-      int const it_f = 6 * it;
-
-      // The scafacos term "potential" here in fact refers to the magnetic
-      // field
-      // So, the torques are given by m \times B
-      auto const dip = p.calc_dip();
-      auto const t = vector_product(
-          dip,
-          Utils::Vector3d(Utils::Span<const double>(&(potentials[it_p]), 3)));
-      // The force is given by G m, where G is a matrix
-      // which comes from the "fields" output of scafacos like this
-      // 0 1 2
-      // 1 3 4
-      // 2 4 5
-      // where the numbers refer to indices in the "field" output from scafacos
-      auto const G = Utils::Vector<Utils::Vector3d, 3>{
-          {fields[it_f + 0], fields[it_f + 1], fields[it_f + 2]},
-          {fields[it_f + 1], fields[it_f + 3], fields[it_f + 4]},
-          {fields[it_f + 2], fields[it_f + 4], fields[it_f + 5]}};
-      auto const f = G * dip;
-
-      // Add to particles
-      p.f.f += dipole.prefactor * f;
-      p.f.torque += dipole.prefactor * t;
-      it++;
 #endif
-    }
-  }
 
-  /* Check that the particle number did not change */
-  if (!dipolar()) {
-    assert(it == fields.size());
-  } else {
-    assert(it == positions.size() / 3);
+ScafacosContextBase *fcs_coulomb() {
+  if (!coulomb_instance) {
+    throw std::runtime_error(
+        "Attempted access to uninitialized Scafacos Coulomb instance.");
   }
+  return coulomb_instance;
 }
 
-void add_pair_force(double q1q2, const double *d, double dist, double *force) {
-  if (dist > get_r_cut())
+#ifdef SCAFACOS_DIPOLES
+static void set_parameters_dipoles_worker(const std::string &method,
+                                          const std::string &params) {
+  delete dipoles_instance;
+  dipoles_instance = nullptr;
+
+  auto *instance = new ScafacosContextDipoles(method, comm_cart, params);
+  if (!instance) {
+    runtimeErrorMsg() << "Scafacos Dipoles failed to initialize";
     return;
-
-  assert(scafacos);
-  const double field = scafacos->pair_force(dist);
-  const double fak = q1q2 * field / dist;
-
-  for (int i = 0; i < 3; i++) {
-    force[i] -= fak * d[i];
   }
-}
-
-double pair_energy(double q1q2, double dist) {
-  if (dist <= get_r_cut())
-    return q1q2 * scafacos->pair_energy(dist);
-  return 0.;
-}
-
-void add_long_range_force() {
-  particles.update_particle_data();
+  dipoles_instance = instance;
 
-  if (scafacos) {
-    if (!dipolar()) {
-      scafacos->run(particles.charges, particles.positions, particles.fields,
-                    particles.potentials);
-    } else {
-#ifdef SCAFACOS_DIPOLES
-      scafacos->run_dipolar(particles.dipoles, particles.positions,
-                            particles.fields, particles.potentials);
-#endif
-    }
-  } else
-    throw std::runtime_error(
-        "Scafacos internal error. Instance pointer is not valid.");
+  instance->set_dipolar(true);
+  instance->update_system_params();
 
-  particles.update_particle_forces();
+  dipole.method = DIPOLAR_SCAFACOS;
+  on_coulomb_change();
 }
 
-double long_range_energy() {
-  if (scafacos) {
-    particles.update_particle_data();
-    if (!dipolar()) {
-      scafacos->run(particles.charges, particles.positions, particles.fields,
-                    particles.potentials);
-      return 0.5 * coulomb.prefactor *
-             std::inner_product(particles.charges.begin(),
-                                particles.charges.end(),
-                                particles.potentials.begin(), 0.0);
-    }
-#ifdef SCAFACOS_DIPOLES
-    scafacos->run_dipolar(particles.dipoles, particles.positions,
-                          particles.fields, particles.potentials);
-    return -0.5 * dipole.prefactor *
-           std::inner_product(particles.dipoles.begin(),
-                              particles.dipoles.end(),
-                              particles.potentials.begin(), 0.0);
+REGISTER_CALLBACK(set_parameters_dipoles_worker)
 #endif
-  }
 
-  return 0.0;
-}
-void set_r_cut_and_tune_local(double r_cut) {
-  particles.update_particle_data();
-
-  scafacos->set_r_cut(r_cut);
-  scafacos->tune(particles.charges, particles.positions);
-}
+static void set_parameters_coulomb_worker(const std::string &method,
+                                          const std::string &params) {
+  delete coulomb_instance;
+  coulomb_instance = nullptr;
 
-REGISTER_CALLBACK(set_r_cut_and_tune_local)
+  auto *instance = new ScafacosContextCoulomb(method, comm_cart, params);
+  if (!instance) {
+    runtimeErrorMsg() << "Scafacos Coulomb failed to initialize";
+    return;
+  }
+  coulomb_instance = instance;
 
-/** Determine runtime for a specific cutoff */
-double time_r_cut(double r_cut) {
-  /* Set cutoff to time */
-  mpi_call(set_r_cut_and_tune_local, r_cut);
-  set_r_cut_and_tune_local(r_cut);
+  instance->set_dipolar(false);
+  instance->update_system_params();
 
-  return time_force_calc(10);
+  coulomb.method = COULOMB_SCAFACOS;
+  on_coulomb_change();
+  instance->tune();
 }
 
-/** Determine the optimal cutoff by bisection */
-void tune_r_cut() {
-  const double tune_limit = 1e-3;
+REGISTER_CALLBACK(set_parameters_coulomb_worker)
 
-  auto const min_box_l = *boost::min_element(box_geo.length());
-  auto const min_local_box_l = *boost::min_element(local_geo.length());
-
-  /* scafacos p3m and Ewald do not accept r_cut 0 for no good reason */
-  double r_min = 1.0;
-  double r_max = std::min(min_local_box_l, min_box_l / 2.0) - skin;
-  double t_min = 0;
-  double t_max = std::numeric_limits<double>::max();
-
-  /* Run bisection */
-  while (std::fabs(r_min - r_max) > tune_limit) {
-    const double dr = 0.5 * (r_max - r_min);
-    const double t_mid = time_r_cut(r_min + dr);
-    t_min = time_r_cut(r_min);
-    t_max = time_r_cut(r_max);
-
-    if (t_min <= 0.0 or t_max <= 0.0) {
-      break;
-    }
-
-    if (t_mid > t_min) {
-      r_max = r_min += dr;
-    } else {
-      r_min += dr;
-    }
-  }
+void set_r_cut_and_tune_local(double r_cut) {
+  coulomb_instance->set_r_cut_and_tune_local(r_cut);
 }
 
-void tune() {
-  particles.update_particle_data();
-
-  /* Check whether we have to do a bisection for the short range cutoff */
-  /* Check if there is a user supplied cutoff */
-  if ((scafacos->has_near) && (scafacos->r_cut() <= 0.0)) {
-    // Tuning of r_cut needs to run on the master node because it relies on
-    // master-slave mode communication
-    if (this_node == 0) {
-      tune_r_cut();
-    } else {
-      return; // Tune on the master node will issue mpi calls
-    }
+void free_handle(bool dipolar) {
+  if (this_node == 0)
+    mpi_call(free_handle, dipolar);
+  if (dipolar) {
+#ifdef SCAFACOS_DIPOLES
+    delete dipoles_instance;
+    dipoles_instance = nullptr;
+#endif
   } else {
-    // ESPResSo is not affected by a short range cutoff. Tune in parallel
-    scafacos->tune(particles.charges, particles.positions);
+    delete coulomb_instance;
+    coulomb_instance = nullptr;
   }
 }
 
-static void set_params_safe(const std::string &method,
-                            const std::string &params, bool dipolar_ia,
-                            int n_part) {
-  if (scafacos) {
-    delete scafacos;
-    scafacos = nullptr;
-  }
-
-  scafacos = new Scafacos(method, comm_cart, params);
-
-  int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0,
-                box_geo.periodic(2) != 0};
+REGISTER_CALLBACK(free_handle)
 
-  scafacos->set_dipolar(dipolar_ia);
-#ifdef DIPOLES
-  if (dipolar_ia) {
-    dipole.method = DIPOLAR_SCAFACOS;
-  }
-#endif
-#ifdef ELECTROSTATICS
-  if (!dipolar_ia) {
-    coulomb.method = COULOMB_SCAFACOS;
-  }
+void set_parameters(const std::string &method, const std::string &params,
+                    bool dipolar) {
+  if (dipolar) {
+#ifdef SCAFACOS_DIPOLES
+    mpi_call_all(set_parameters_dipoles_worker, method, params);
 #endif
-  scafacos->set_common_parameters(box_geo.length().data(), per, n_part);
-
-  on_coulomb_change();
-
-  if (!dipolar_ia) {
-    tune();
+  } else {
+    mpi_call_all(set_parameters_coulomb_worker, method, params);
   }
 }
 
-REGISTER_CALLBACK(set_params_safe)
-
-/** Bend result from scafacos back to original format */
-std::string get_method_and_parameters() {
-  if (!scafacos) {
+std::string get_method_and_parameters(bool dipolar) {
+  if (dipolar) {
+#ifdef SCAFACOS_DIPOLES
+    return fcs_dipoles()->get_method_and_parameters();
+#else
     return std::string();
+#endif
   }
-
-  std::string p = scafacos->get_method() + " " + scafacos->get_parameters();
-
-  std::replace(p.begin(), p.end(), ',', ' ');
-
-  return p;
-}
-
-double get_r_cut() {
-  if (scafacos) {
-    if (!scafacos->has_near)
-      return 0;
-    return scafacos->r_cut();
-  }
-  return 0.0;
-}
-
-void set_parameters(const std::string &method, const std::string &params,
-                    bool dipolar_ia) {
-  mpi_call_all(set_params_safe, method, params, dipolar_ia, get_n_part());
-}
-
-bool dipolar() {
-  if (scafacos)
-    return scafacos->dipolar();
-  throw std::runtime_error("Scafacos not initialized");
-}
-
-void set_dipolar(bool d) {
-  if (scafacos) {
-    scafacos->set_dipolar(d);
-  } else
-    runtimeErrorMsg() << "Scafacos not initialized.";
-}
-
-void free_handle() {
-  if (this_node == 0)
-    mpi_call(free_handle);
-  if (scafacos) {
-    delete scafacos;
-    scafacos = nullptr;
-  }
-}
-
-REGISTER_CALLBACK(free_handle)
-
-void update_system_params() {
-  // If scafacos is not active, do nothing
-  if (!scafacos) {
-    throw std::runtime_error("Scafacos object not there");
-  }
-
-  int per[3] = {box_geo.periodic(0) != 0, box_geo.periodic(1) != 0,
-                box_geo.periodic(2) != 0};
-
-  auto const n_part = boost::mpi::all_reduce(
-      comm_cart, cell_structure.local_particles().size(), std::plus<>());
-  scafacos->set_common_parameters(box_geo.length().data(), per, n_part);
+  return fcs_coulomb()->get_method_and_parameters();
 }
 
 } // namespace Scafacos
diff --git a/src/core/electrostatics_magnetostatics/scafacos.hpp b/src/core/electrostatics_magnetostatics/scafacos.hpp
index a1272eeb64e..4587623e96c 100644
--- a/src/core/electrostatics_magnetostatics/scafacos.hpp
+++ b/src/core/electrostatics_magnetostatics/scafacos.hpp
@@ -21,46 +21,41 @@
 
 /** \file
  * This file contains the c-type wrapper interface to the (oop-) scafacos
- * interface. */
+ * interface.
+ */
 
 #ifndef ES_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOS_HPP
 #define ES_CORE_ELECTROSTATICS_MAGNETOSTATICS_SCAFACOS_HPP
 
 #include "config.hpp"
 
+#if defined(SCAFACOS)
+
+#include "electrostatics_magnetostatics/ScafacosContextBase.hpp"
+
+#include <utils/Vector.hpp>
+
 #include <list>
 #include <string>
 
 namespace Scafacos {
-#if defined(SCAFACOS)
-/** Near-field pair force */
-void add_pair_force(double q1q2, const double *d, double dist, double *force);
-/** Near-field pair energy */
-double pair_energy(double q1q2, double dist);
-/** Long range part */
-void add_long_range_force();
-/** Calculate long range energy contribution */
-double long_range_energy();
-/** Get parameters */
-std::string get_method_and_parameters();
-/** Set parameters */
-void set_parameters(const std::string &method, const std::string &params,
-                    bool dipolar);
-double get_r_cut();
-
-/** Is scafacos used for dipolar interactions */
-bool dipolar();
 
-/** Choose whether scafacos is used for dipolar interactions */
-void set_dipolar(bool d);
+/** @brief Access the per-MPI-node ScaFaCoS Coulomb instance */
+ScafacosContextBase *fcs_coulomb();
+#ifdef SCAFACOS_DIPOLES
+/** @brief Access the per-MPI-node ScaFaCoS dipoles instance */
+ScafacosContextBase *fcs_dipoles();
+#endif
 
-/** Reinit scafacos number of particles, box shape and periodicity */
-void update_system_params();
+std::string get_method_and_parameters(bool dipolar);
+void set_parameters(const std::string &method, const std::string &params,
+                    bool dipolar);
+void free_handle(bool dipolar);
 
-#endif /* SCAFACOS */
+void set_r_cut_and_tune_local(double r_cut);
 
 std::list<std::string> available_methods();
 
-void free_handle();
 } // namespace Scafacos
+#endif // SCAFACOS
 #endif
diff --git a/src/core/event.cpp b/src/core/event.cpp
index 9050f774e9c..08f5af1372d 100644
--- a/src/core/event.cpp
+++ b/src/core/event.cpp
@@ -280,12 +280,12 @@ void on_parameter_change(int field) {
 #ifdef SCAFACOS
 #ifdef ELECTROSTATICS
     if (coulomb.method == COULOMB_SCAFACOS) {
-      Scafacos::update_system_params();
+      Scafacos::fcs_coulomb()->update_system_params();
     }
 #endif
-#ifdef DIPOLES
+#ifdef SCAFACOS_DIPOLES
     if (dipole.method == DIPOLAR_SCAFACOS) {
-      Scafacos::update_system_params();
+      Scafacos::fcs_dipoles()->update_system_params();
     }
 #endif
 #endif
diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx
index 86f157782a7..7f64ebad395 100644
--- a/src/python/espressomd/electrostatics.pyx
+++ b/src/python/espressomd/electrostatics.pyx
@@ -691,5 +691,5 @@ IF ELECTROSTATICS:
 
             def _deactivate_method(self):
                 super()._deactivate_method()
-                scafacos.free_handle()
+                scafacos.free_handle(self.dipolar)
                 mpi_bcast_coulomb_params()
diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx
index 85d2d33c4f7..7e3f29f138f 100644
--- a/src/python/espressomd/magnetostatics.pyx
+++ b/src/python/espressomd/magnetostatics.pyx
@@ -345,7 +345,7 @@ IF DIPOLES == 1:
 
             def _deactivate_method(self):
                 dipole.method = DIPOLAR_NONE
-                scafacos.free_handle()
+                scafacos.free_handle(self.dipolar)
                 mpi_bcast_coulomb_params()
 
             def default_params(self):
diff --git a/src/python/espressomd/scafacos.pxd b/src/python/espressomd/scafacos.pxd
index d8c4d8465bc..d215ddea7fa 100644
--- a/src/python/espressomd/scafacos.pxd
+++ b/src/python/espressomd/scafacos.pxd
@@ -25,6 +25,6 @@ from libcpp.list cimport list
 IF SCAFACOS:
     cdef extern from "electrostatics_magnetostatics/scafacos.hpp" namespace "Scafacos":
         void set_parameters(string & method_name, string & params, bool dipolar)
-        string get_method_and_parameters()
+        string get_method_and_parameters(bool dipolar) except+
+        void free_handle(bool dipolar)
         list[string] available_methods_core "Scafacos::available_methods" ()
-        void free_handle()
diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx
index da6133a81e0..d258c267e81 100644
--- a/src/python/espressomd/scafacos.pyx
+++ b/src/python/espressomd/scafacos.pyx
@@ -60,27 +60,24 @@ IF SCAFACOS == 1:
             # Parameters are returned as strings
             # First word is method name, rest are key value pairs
             # which we convert to a dict
-            p = to_str(get_method_and_parameters()).split(" ")
+            p = to_str(get_method_and_parameters(self.dipolar)).split(" ")
             res = {}
-            res["method_name"] = p[0]
+            res["method_name"] = p.pop(0)
 
             method_params = {}
-            i = 1
-            while i < len(p):
-                pname = p[i]
-                i += 1
+            while p:
+                pname = p.pop(0)
 
                 # The first item after the parameter name is always a value
                 # But in the case of array-like properties, there can also
                 # follow several values. Therefore, we treat the next
                 # words as part of the value, if they begin with a digit
-                pvalues = [p[i]]
-                i += 1
-                if i >= len(p):
-                    break
-                while p[i][:1] in "-0123456789":
-                    pvalues.append(p[i])
-                    i += 1
+                pvalues = [p.pop(0)]
+                while p:
+                    if p[0][0] in "-0123456789":
+                        pvalues.append(p.pop(0))
+                    else:
+                        break
 
                 # If there is one value, cast away the list
                 if len(pvalues) == 1:
@@ -89,8 +86,6 @@ IF SCAFACOS == 1:
                     # Cast array elements to strings and join them by commas
                     # to achieve consistency with setting array-likes
                     # such as "pnfft_n":"128,128,128"
-                    for j in range(len(pvalues)):
-                        pvalues[j] = str(pvalues[j])
                     pvalues = ",".join(pvalues)
                 method_params[pname] = pvalues
 
@@ -108,18 +103,6 @@ IF SCAFACOS == 1:
             return res
 
         def _set_params_in_es_core(self):
-            # Verify that scafacos is not used for electrostatics and dipoles
-            # at the same time
-            IF ELECTROSTATICS == 1:
-                if self.dipolar and < int > electrostatics.coulomb.method == <int > electrostatics.COULOMB_SCAFACOS:
-                    raise Exception(
-                        "Scafacos cannot be used for dipoles and charges at the same time")
-
-            IF DIPOLES == 1:
-                if not self.dipolar and < int > magnetostatics.dipole.method == <int > magnetostatics.DIPOLAR_SCAFACOS:
-                    raise Exception(
-                        "Scafacos cannot be used for dipoles and charges at the same time")
-
             # Convert the parameter dictionary to a list of strings
             method_params = self._params["method_params"]
             param_string = ""
diff --git a/src/scafacos/Scafacos.cpp b/src/scafacos/Scafacos.cpp
index a9b2711c8d5..881a7518be4 100644
--- a/src/scafacos/Scafacos.cpp
+++ b/src/scafacos/Scafacos.cpp
@@ -124,10 +124,11 @@ void Scafacos::run(std::vector<double> &charges, std::vector<double> &positions,
   fields.resize(3 * local_n_part);
   potentials.resize(local_n_part);
 
-  handle_error(fcs_tune(handle, local_n_part, &(positions[0]), &(charges[0])));
+  handle_error(
+      fcs_tune(handle, local_n_part, positions.data(), charges.data()));
 
-  handle_error(fcs_run(handle, local_n_part, &(positions[0]), &(charges[0]),
-                       &(fields[0]), &(potentials[0])));
+  handle_error(fcs_run(handle, local_n_part, positions.data(), charges.data(),
+                       fields.data(), potentials.data()));
 }
 
 #ifdef FCS_ENABLE_DIPOLES
@@ -142,8 +143,8 @@ void Scafacos::run_dipolar(std::vector<double> &dipoles,
   potentials.resize(3 * local_n_part);
 
   handle_error(fcs_set_dipole_particles(handle, static_cast<int>(local_n_part),
-                                        &(positions[0]), &(dipoles[0]),
-                                        &(fields[0]), &(potentials[0])));
+                                        positions.data(), dipoles.data(),
+                                        fields.data(), potentials.data()));
   handle_error(fcs_run(handle, 0, nullptr, nullptr, nullptr, nullptr));
 }
 #endif
@@ -151,7 +152,7 @@ void Scafacos::run_dipolar(std::vector<double> &dipoles,
 void Scafacos::tune(std::vector<double> &charges,
                     std::vector<double> &positions) {
   handle_error(
-      fcs_tune(handle, charges.size(), &(positions[0]), &(charges[0])));
+      fcs_tune(handle, charges.size(), positions.data(), charges.data()));
 }
 
 void Scafacos::set_common_parameters(const double *box_l,
@@ -162,15 +163,10 @@ void Scafacos::set_common_parameters(const double *box_l,
   boxa[0] = box_l[0];
   boxb[1] = box_l[1];
   boxc[2] = box_l[2];
-  // Does scafacos calculate the near field part
-  // For charges, if the method supports it, Es calculates near field
-  int sr = 0;
-  if (!dipolar() && has_near) {
-    sr = 0;
-  } else {
-    // Scafacos does near field calc
-    sr = 1;
-  }
+  // Does scafacos calculate the near field part?
+  // For charges, ESPResSo calculates near field if the method supports it.
+  // For dipoles, ScaFaCoS calculates near field.
+  int const sr = dipolar() || !has_near;
   handle_error(fcs_set_common(handle, sr, boxa, boxb, boxc, off, periodicity,
                               total_particles));
 #ifdef FCS_ENABLE_DIPOLES
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index a40807f3964..2f9aba81870 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -109,6 +109,7 @@ python_test(FILE p3m_gpu.py MAX_NUM_PROC 4 LABELS gpu)
 python_test(FILE particle.py MAX_NUM_PROC 4)
 python_test(FILE pressure.py MAX_NUM_PROC 4)
 python_test(FILE scafacos_dipoles_1d_2d.py MAX_NUM_PROC 4)
+python_test(FILE scafacos_interface.py MAX_NUM_PROC 2)
 python_test(FILE tabulated.py MAX_NUM_PROC 2)
 python_test(FILE particle_slice.py MAX_NUM_PROC 4)
 python_test(FILE rigid_bond.py MAX_NUM_PROC 4)
diff --git a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py
index 4d160e2b09f..a2f4ac4ecd0 100644
--- a/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py
+++ b/testsuite/python/dipolar_mdlc_p3m_scafacos_p2nfft.py
@@ -1,5 +1,4 @@
 # Copyright (C) 2010-2019 The ESPResSo project
-
 #
 # This file is part of ESPResSo.
 #
@@ -42,12 +41,14 @@ class Dipolar_p3m_mdlc_p2nfft(ut.TestCase):
     system.time_step = 0.01
     system.cell_system.skin = .4
     system.periodicity = [1, 1, 1]
-    system.thermostat.turn_off()
 
     def tearDown(self):
         self.system.part.clear()
         self.system.actors.clear()
 
+    def vector_error(self, a, b):
+        return np.sum(np.linalg.norm(a - b, axis=1)) / np.sqrt(a.shape[0])
+
     def test_mdlc(self):
         s = self.system
         rho = 0.3
@@ -57,34 +58,25 @@ def test_mdlc(self):
         n_particle = 100
 
         particle_radius = 0.5
-        box_l = pow(((4 * n_particle * np.pi) / (3 * rho)),
-                    1.0 / 3.0) * particle_radius
+        box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius
         s.box_l = 3 * [box_l]
-        f = open(abspath("data/mdlc_reference_data_energy.dat"))
-        ref_E = float(f.readline())
-        f.close()
+        ref_E_path = abspath("data/mdlc_reference_data_energy.dat")
+        ref_E = float(np.genfromtxt(ref_E_path))
 
         # Particles
         data = np.genfromtxt(
             abspath("data/mdlc_reference_data_forces_torques.dat"))
-        for p in data[:, :]:
-            s.part.add(id=int(p[0]), pos=p[1:4], dip=p[4:7])
+        s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
         s.part[:].rotation = (1, 1, 1)
 
         p3m = magnetostatics.DipolarP3M(prefactor=1, mesh=32, accuracy=1E-4)
         dlc = magnetostatic_extensions.DLC(maxPWerror=1E-5, gap_size=2.)
         s.actors.add(p3m)
         s.actors.add(dlc)
-        s.thermostat.turn_off()
         s.integrator.run(0)
-        err_f = np.sum(np.linalg.norm(
-            s.part[:].f - data[:, 7:10], axis=1)) / np.sqrt(data.shape[0])
-        err_t = np.sum(np.linalg.norm(
-            s.part[:].torque_lab - data[:, 10:13], axis=1)) / np.sqrt(data.shape[0])
+        err_f = self.vector_error(s.part[:].f, data[:, 7:10])
+        err_t = self.vector_error(s.part[:].torque_lab, data[:, 10:13])
         err_e = s.analysis.energy()["dipolar"] - ref_E
-        print("Energy difference", err_e)
-        print("Force difference", err_f)
-        print("Torque difference", err_t)
 
         tol_f = 2E-3
         tol_t = 2E-3
@@ -94,9 +86,6 @@ def test_mdlc(self):
         self.assertLessEqual(abs(err_t), tol_t, "Torque difference too large")
         self.assertLessEqual(abs(err_f), tol_f, "Force difference too large")
 
-        del s.actors[0]
-        del s.actors[0]
-
     def test_p3m(self):
         s = self.system
         rho = 0.09
@@ -106,14 +95,12 @@ def test_p3m(self):
         n_particle = 1000
 
         particle_radius = 1
-        box_l = pow(((4 * n_particle * np.pi) / (3 * rho)),
-                    1.0 / 3.0) * particle_radius
+        box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius
         s.box_l = 3 * [box_l]
 
         # Particles
         data = np.genfromtxt(abspath("data/p3m_magnetostatics_system.data"))
-        for p in data[:, :]:
-            s.part.add(id=int(p[0]), pos=p[1:4], dip=p[4:7])
+        s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
         s.part[:].rotation = (1, 1, 1)
 
         p3m = magnetostatics.DipolarP3M(
@@ -122,15 +109,10 @@ def test_p3m(self):
         s.integrator.run(0)
         expected = np.genfromtxt(
             abspath("data/p3m_magnetostatics_expected.data"))[:, 1:]
-        err_f = np.sum(np.linalg.norm(
-            s.part[:].f - expected[:, 0:3], axis=1)) / np.sqrt(data.shape[0])
-        err_t = np.sum(np.linalg.norm(
-            s.part[:].torque_lab - expected[:, 3:6], axis=1)) / np.sqrt(data.shape[0])
+        err_f = self.vector_error(s.part[:].f, expected[:, 0:3])
+        err_t = self.vector_error(s.part[:].torque_lab, expected[:, 3:6])
         ref_E = 5.570
         err_e = s.analysis.energy()["dipolar"] - ref_E
-        print("Energy difference", err_e)
-        print("Force difference", err_f)
-        print("Torque difference", err_t)
 
         tol_f = 2E-3
         tol_t = 2E-3
@@ -150,15 +132,13 @@ def test_scafacos_dipoles(self):
         n_particle = 1000
 
         particle_radius = 1
-        box_l = pow(((4 * n_particle * np.pi) / (3 * rho)),
-                    1.0 / 3.0) * particle_radius
+        box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius
         s.box_l = 3 * [box_l]
 
         # Particles
         data = np.genfromtxt(abspath("data/p3m_magnetostatics_system.data"))
-        for p in data[:, :]:
-            s.part.add(id=int(p[0]), pos=p[1:4],
-                       dip=p[4:7], rotation=(1, 1, 1))
+        s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
+        s.part[:].rotation = (1, 1, 1)
 
         scafacos = magnetostatics.Scafacos(
             prefactor=1,
@@ -177,15 +157,10 @@ def test_scafacos_dipoles(self):
         s.integrator.run(0)
         expected = np.genfromtxt(
             abspath("data/p3m_magnetostatics_expected.data"))[:, 1:]
-        err_f = np.sum(np.linalg.norm(
-            s.part[:].f - expected[:, 0:3], axis=1)) / np.sqrt(data.shape[0])
-        err_t = np.sum(np.linalg.norm(
-            s.part[:].torque_lab - expected[:, 3:6], axis=1)) / np.sqrt(data.shape[0])
+        err_f = self.vector_error(s.part[:].f, expected[:, 0:3])
+        err_t = self.vector_error(s.part[:].torque_lab, expected[:, 3:6])
         ref_E = 5.570
         err_e = s.analysis.energy()["dipolar"] - ref_E
-        print("Energy difference", err_e)
-        print("Force difference", err_f)
-        print("Torque difference", err_t)
 
         tol_f = 2E-3
         tol_t = 2E-3
diff --git a/testsuite/python/scafacos_dipoles_1d_2d.py b/testsuite/python/scafacos_dipoles_1d_2d.py
index 4442236e607..fb9c2e9f5c5 100644
--- a/testsuite/python/scafacos_dipoles_1d_2d.py
+++ b/testsuite/python/scafacos_dipoles_1d_2d.py
@@ -31,7 +31,21 @@
 @utx.skipIfMissingFeatures(["SCAFACOS_DIPOLES"])
 class Scafacos1d2d(ut.TestCase):
 
+    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system.time_step = 0.01
+    system.cell_system.skin = 0.5
+    system.periodicity = [1, 1, 1]
+
+    def tearDown(self):
+        self.system.part.clear()
+        self.system.actors.clear()
+        self.system.periodicity = [1, 1, 1]
+
+    def vector_error(self, a, b):
+        return np.sum(np.linalg.norm(a - b, axis=1)) / np.sqrt(a.shape[0])
+
     def test_scafacos(self):
+        s = self.system
         rho = 0.3
 
         # This is only for box size calculation. The actual particle number is
@@ -40,18 +54,10 @@ def test_scafacos(self):
 
         particle_radius = 0.5
 
-        #################################################
-
-        box_l = pow(((4 * n_particle * np.pi) / (3 * rho)),
-                    1.0 / 3.0) * particle_radius
-        skin = 0.5
-
-        s = espressomd.System(box_l=[1.0, 1.0, 1.0])
-        # give Espresso some parameters
-        s.time_step = 0.01
-        s.cell_system.skin = skin
+        box_l = np.cbrt(4 * n_particle * np.pi / (3 * rho)) * particle_radius
         s.box_l = 3 * [box_l]
-        for dim in 2, 1:
+
+        for dim in (2, 1):
             print("Dimension", dim)
 
             # Read reference data
@@ -62,15 +68,14 @@ def test_scafacos(self):
                 s.periodicity = [1, 0, 0]
                 file_prefix = "data/scafacos_dipoles_1d"
 
-            with open(abspath(file_prefix + "_reference_data_energy.dat")) as f:
-                ref_E = float(f.readline())
+            ref_E_path = abspath(file_prefix + "_reference_data_energy.dat")
+            ref_E = float(np.genfromtxt(ref_E_path))
 
             # Particles
             data = np.genfromtxt(abspath(
                 file_prefix + "_reference_data_forces_torques.dat"))
-            for p in data[:, :]:
-                s.part.add(
-                    id=int(p[0]), pos=p[1:4], dip=p[4:7], rotation=(1, 1, 1))
+            s.part.add(pos=data[:, 1:4], dip=data[:, 4:7])
+            s.part[:].rotation = (1, 1, 1)
 
             if dim == 2:
                 scafacos = magnetostatics.Scafacos(
@@ -92,42 +97,33 @@ def test_scafacos(self):
                 s.box_l = np.array((1, 1, 1.3)) * box_l
 
             else:
-                if dim == 1:
-                    # 1d periodic in x
-                    scafacos = magnetostatics.Scafacos(
-                        prefactor=1,
-                        method_name="p2nfft",
-                        method_params={
-                            "p2nfft_verbose_tuning": 1,
-                            "pnfft_N": "32,128,128",
-                            "pnfft_direct": 0,
-                            "p2nfft_r_cut": 2.855,
-                            "p2nfft_alpha": "1.5",
-                            "p2nfft_intpol_order": "-1",
-                            "p2nfft_reg_kernel_name": "ewald",
-                            "p2nfft_p": "16",
-                            "p2nfft_ignore_tolerance": "1",
-                            "pnfft_window_name": "bspline",
-                            "pnfft_m": "8",
-                            "pnfft_diff_ik": "1",
-                            "p2nfft_epsB": "0.125"})
-                    s.box_l = np.array((1, 1, 1)) * box_l
-                    s.actors.add(scafacos)
-                else:
-                    raise Exception("This shouldn't happen.")
-            s.thermostat.turn_off()
+                # 1d periodic in x
+                scafacos = magnetostatics.Scafacos(
+                    prefactor=1,
+                    method_name="p2nfft",
+                    method_params={
+                        "p2nfft_verbose_tuning": 1,
+                        "pnfft_N": "32,128,128",
+                        "pnfft_direct": 0,
+                        "p2nfft_r_cut": 2.855,
+                        "p2nfft_alpha": "1.5",
+                        "p2nfft_intpol_order": "-1",
+                        "p2nfft_reg_kernel_name": "ewald",
+                        "p2nfft_p": "16",
+                        "p2nfft_ignore_tolerance": "1",
+                        "pnfft_window_name": "bspline",
+                        "pnfft_m": "8",
+                        "pnfft_diff_ik": "1",
+                        "p2nfft_epsB": "0.125"})
+                s.box_l = np.array((1, 1, 1)) * box_l
+                s.actors.add(scafacos)
             s.integrator.run(0)
 
             # Calculate errors
 
-            err_f = np.sum(np.linalg.norm(
-                s.part[:].f - data[:, 7:10], axis=1)) / np.sqrt(data.shape[0])
-            err_t = np.sum(np.linalg.norm(
-                s.part[:].torque_lab - data[:, 10:13], axis=1)) / np.sqrt(data.shape[0])
+            err_f = self.vector_error(s.part[:].f, data[:, 7:10])
+            err_t = self.vector_error(s.part[:].torque_lab, data[:, 10:13])
             err_e = s.analysis.energy()["dipolar"] - ref_E
-            print("Energy difference", err_e)
-            print("Force difference", err_f)
-            print("Torque difference", err_t)
 
             tol_f = 2E-3
             tol_t = 2E-3
@@ -141,7 +137,7 @@ def test_scafacos(self):
                 abs(err_f), tol_f, "Force difference too large")
 
             s.part.clear()
-            del s.actors[0]
+            s.actors.clear()
 
 
 if __name__ == "__main__":
diff --git a/testsuite/python/scafacos_interface.py b/testsuite/python/scafacos_interface.py
new file mode 100644
index 00000000000..6a1653c6523
--- /dev/null
+++ b/testsuite/python/scafacos_interface.py
@@ -0,0 +1,193 @@
+#
+# Copyright (C) 2020 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/>.
+#
+
+import numpy as np
+import unittest as ut
+import unittest_decorators as utx
+import espressomd
+import espressomd.electrostatics
+import espressomd.magnetostatics
+import espressomd.scafacos
+
+
+@utx.skipIfMissingFeatures(["SCAFACOS"])
+class ScafacosInterface(ut.TestCase):
+
+    system = espressomd.System(box_l=3 * [5])
+    system.time_step = 0.01
+    system.cell_system.skin = 0.5
+    system.periodicity = [1, 1, 1]
+
+    def tearDown(self):
+        self.system.part.clear()
+        self.system.actors.clear()
+
+    def test_available_methods(self):
+        # all methods that are accessible from the ScaFaCoS bridge in ESPResSo
+        scafacos_methods = {
+            "direct", "ewald", "fmm", "memd", "mmm1d", "mmm2d",
+            "p2nfft", "p3m", "pepc", "pp3mg", "vmg", "wolf"}
+
+        # all methods that were compiled in when building ScaFaCoS
+        available_methods = espressomd.scafacos.available_methods()
+        self.assertGreaterEqual(len(available_methods), 1)
+        for method in available_methods:
+            self.assertIn(method, scafacos_methods)
+
+    def test_actor_coulomb(self):
+        system = self.system
+
+        system.actors.add(espressomd.electrostatics.Scafacos(
+            prefactor=0.5,
+            method_name="p3m",
+            method_params={
+                "p3m_r_cut": 1.0,
+                "p3m_grid": 32,
+                "p3m_cao": 7}))
+        actor = system.actors[0]
+        params = actor.get_params()
+        self.assertEqual(params["prefactor"], 0.5)
+        self.assertEqual(params["method_name"], "p3m")
+        self.assertEqual(params["method_params"],
+                         {'p3m_cao': '7', 'p3m_r_cut': '1.0', 'p3m_grid': '32'})
+
+    @utx.skipIfMissingFeatures(["SCAFACOS_DIPOLES"])
+    def test_actor_dipoles(self):
+        system = self.system
+
+        method_params = {
+            "p2nfft_verbose_tuning": "0",
+            "pnfft_N": "32,32,32",
+            "pnfft_n": "32,32,32",
+            "pnfft_window_name": "bspline",
+            "pnfft_m": "4",
+            "p2nfft_ignore_tolerance": "1",
+            "pnfft_diff_ik": "0",
+            "p2nfft_r_cut": "11",
+            "p2nfft_alpha": "0.37"}
+        system.actors.add(espressomd.magnetostatics.Scafacos(
+            prefactor=1.2,
+            method_name="p2nfft",
+            method_params=method_params))
+        actor = system.actors[0]
+        params = actor.get_params()
+        self.assertEqual(params["prefactor"], 1.2)
+        self.assertEqual(params["method_name"], "p2nfft")
+        self.assertEqual(params["method_params"], method_params)
+
+    def p3m_data(self):
+        system = self.system
+
+        p3m = espressomd.electrostatics.P3M(
+            prefactor=0.5,
+            accuracy=5e-4,
+            mesh=32,
+            cao=7,
+            r_cut=1.0)
+        system.actors.add(p3m)
+
+        dp3m = espressomd.magnetostatics.DipolarP3M(
+            prefactor=1.0,
+            accuracy=1e-5,
+            cao=7,
+            mesh=48,
+            epsilon="metallic")
+        system.actors.add(dp3m)
+
+        system.integrator.run(0, recalc_forces=True)
+        ref_E_coulomb = system.analysis.energy()["coulomb"]
+        ref_E_dipoles = system.analysis.energy()["dipolar"]
+        ref_forces = np.copy(system.part[:].f)
+        ref_torques = np.copy(system.part[:].torque_lab)
+
+        system.actors.clear()
+
+        return (ref_E_coulomb, ref_E_dipoles, ref_forces, ref_torques)
+
+    def fcs_data(self):
+        system = self.system
+
+        scafacos_coulomb = espressomd.electrostatics.Scafacos(
+            prefactor=0.5,
+            method_name="p3m",
+            method_params={
+                "p3m_r_cut": 1.0,
+                "p3m_grid": 32,
+                "p3m_cao": 7})
+        system.actors.add(scafacos_coulomb)
+
+        scafacos_dipoles = espressomd.magnetostatics.Scafacos(
+            prefactor=1.0,
+            method_name="p2nfft",
+            method_params={
+                "p2nfft_verbose_tuning": 0,
+                "pnfft_N": "32,32,32",
+                "pnfft_n": "32,32,32",
+                "pnfft_window_name": "bspline",
+                "pnfft_m": "4",
+                "p2nfft_ignore_tolerance": "1",
+                "pnfft_diff_ik": "0",
+                "p2nfft_r_cut": "11",
+                "p2nfft_alpha": "0.37"})
+        system.actors.add(scafacos_dipoles)
+
+        system.integrator.run(0, recalc_forces=True)
+        ref_E_coulomb = system.analysis.energy()["coulomb"]
+        ref_E_dipoles = system.analysis.energy()["dipolar"]
+        ref_forces = np.copy(system.part[:].f)
+        ref_torques = np.copy(system.part[:].torque_lab)
+
+        system.actors.clear()
+
+        return (ref_E_coulomb, ref_E_dipoles, ref_forces, ref_torques)
+
+    @utx.skipIfMissingFeatures(["SCAFACOS_DIPOLES", "LENNARD_JONES"])
+    def test_electrostatics_plus_magnetostatics(self):
+        # check that two instances of ScaFaCoS can be used
+        system = self.system
+
+        # add particles
+        N = 100
+        np.random.seed(42)
+        system.part.add(pos=np.random.uniform(0, system.box_l[0], (N, 3)),
+                        dip=np.random.uniform(0, 1, (N, 3)),
+                        q=np.sign((np.arange(N) % 2) * 2 - 1),
+                        rotation=N * [(1, 1, 1)])
+
+        # minimize system
+        system.non_bonded_inter[0, 0].lennard_jones.set_params(
+            epsilon=1.0, sigma=1.0, cutoff=2**(1.0 / 6.0), shift="auto")
+        system.integrator.set_steepest_descent(
+            f_max=10,
+            gamma=0.001,
+            max_displacement=0.01)
+        system.integrator.run(100)
+
+        # compute forces and energies
+        p3m_E_coulomb, p3m_E_dipoles, p3m_forces, p3m_torques = self.p3m_data()
+        fcs_E_coulomb, fcs_E_dipoles, fcs_forces, fcs_torques = self.fcs_data()
+
+        self.assertAlmostEqual(fcs_E_coulomb, p3m_E_coulomb, delta=1e-4)
+        self.assertAlmostEqual(fcs_E_dipoles, p3m_E_dipoles, delta=1e-4)
+        np.testing.assert_allclose(fcs_forces, p3m_forces, rtol=1e-3)
+        np.testing.assert_allclose(fcs_torques, p3m_torques, rtol=1e-3)
+
+
+if __name__ == "__main__":
+    ut.main()