GPU Electrostatics methods do not obey temperature change
diff --git a/doc/sphinx/installation.rst b/doc/sphinx/installation.rst
index f74b0fc23..bc73c07be 100644
--- a/doc/sphinx/installation.rst
+++ b/doc/sphinx/installation.rst
@@ -560,7 +560,7 @@ system. To activate the P3M electrostatics solver, execute
 ::
 
     from espressomd.electrostatics import P3M p3m=P3M(accuracy=1E-3,
-    bjerrum_length=1) system.actors.add(p3m)
+    prefactor=1) system.actors.add(p3m)
 
 The integrator uses by default the velocity verlet algorithm and is
 created by the system class. To perform an integration step, execute
diff --git a/doc/sphinx/inter.rst b/doc/sphinx/inter.rst
index 0ac391907..216383e6e 100644
--- a/doc/sphinx/inter.rst
+++ b/doc/sphinx/inter.rst
@@ -1315,14 +1315,17 @@ The Coulomb (or electrostatic) interaction is defined as
 follows. For a pair of particles at distance :math:`r` with charges
 :math:`q_1` and :math:`q_2`, the interaction is given by
 
-.. math:: U_C(r)=l_B k_B T\frac{q_1 q_2}{r}.
+.. math:: U_C(r)=C \cdot \frac{q_1 q_2}{r}.
 
-where :math:`l_B = e_o^2 / (4 \pi \epsilon k_B T)` denotes the Bjerrum
-length, which measures the strength of the electrostatic interaction. As
-a special case, when the thermostat is switched off, the value of
-Bjerrum length you enter is treated as :math:`l_B k_B T` rather than
-:math:`l_B`. This is used to perform an NVE integration (see also
-section :ref:`\`\`thermostat\`\`\: Setting up the thermostat`).
+where 
+
+.. math::  
+   C=\frac{1}{4\pi \epsilon_0 \epsilon_r} 
+   :label: coulomb_prefactor
+    
+is a prefactor which can be set by the user. 
+The commonly used Bjerrum length :math:`l_B = e_o^2 / (4 \pi \e\epsilon k_B T)` is the length at which the Coulomb energy between two unit charges is equal to the thermal energy :math:`k_B T`.
+Based on the this length, the prefactor is given by :math:`C=l_B *k_B T`.
 
 Computing electrostatic interactions is computationally very expensive.
 |es| features some state-of-the-art algorithms to deal with these
@@ -1330,10 +1333,9 @@ interactions as efficiently as possible, but almost all of them require
 some knowledge to use them properly. Uneducated use can result in
 completely unphysical simulations.
 
-Coulomb interactions have to be added to the actors of the system to become
-active. This prevents the simultaneous use of multiple electrostatic solvers.
-
-.. todo:: Document missing implementation for actor.remove()
+Coulomb interactions have to be added to the list of active actors of the system object to become
+active. This is done by calling the add-method of :attr:`espressomd.system.System.actors`.
+Only one electrostatics method can be active at any time.
 
 Note that using the electrostatic interaction also requires assigning charges to
 the particles via the particle property
@@ -1349,16 +1351,19 @@ activation to achieve the given accuracy::
     from espressomd import electrostatics
     
     system = espressomd.System()
-    solver = electrostatics.<SOLVER>(bjerrum_length = 1.0, <ADDITIONAL REQUIRED PARAMETERS>)
+    solver = electrostatics.<SOLVER>(prefactor = C, <ADDITIONAL REQUIRED PARAMETERS>)
     system.actors.add(solver)
 
+where the prefactor :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor`
+
+
 Coulomb P3M
 ~~~~~~~~~~~
 
 :class:`espressomd.electrostatics.P3M`
 
 Required parameters:
-    * bjerrum_length
+    * prefactor
     * accuracy
 
 For this feature to work, you need to have the ``fftw3`` library
@@ -1383,7 +1388,7 @@ To manually tune or retune P3M, call :meth:`espresso.electrostatics.P3M.Tune`.
 Note, however, that this is a method the P3M object inherited from
 :attr:`espressomd.electrostatics.ElectrostaticInteraction`. 
 All parameters passed to the method are fixed in the tuning routine. If not
-specified in the ``Tune()`` method, the parameters ``bjerrum_length`` and
+specified in the ``Tune()`` method, the parameters ``prefactor`` and
 ``accuracy`` are reused.
 
 It is not easy to calculate the various parameters of the P3M method
@@ -1415,7 +1420,7 @@ Coulomb P3M on GPU
 :class:`espressomd.electrostatics.P3MGPU`
 
 Required parameters:
-    * bjerrum_length
+    * prefactor
     * accuracy
 
 The GPU implementation of P3M calculates the far field portion on the GPU. 
@@ -1426,14 +1431,12 @@ with only single precision which limits the maximum precision. The algorithm
 does not work in combination with the electrostatic extensions :ref:`ICC` and
 :ref:`ELC`.
 
-.. todo:: Check P3MGPU for non-cubic boxes, and also for cubic.
-
 Coulomb Ewald GPU
 ~~~~~~~~~~~~~~~~~
 
 
 Required parameters:
-    * bjerrum_length
+    * prefactor
     * accuracy
     * precision
     * K_max
@@ -1481,8 +1484,9 @@ For a list of all parameters see :attr:`espressomd.electrostatics.DH` or :attr:`
 
 Uses the Debye-Hückel electrostatic potential defined by
 
-  .. math:: U^{C-DH} = l_B k_B T \frac{q_1 q_2 exp(-\kappa r)}{r}\quad \mathrm{for}\quad r<r_{\mathrm{cut}}
+  .. math:: U^{C-DH} = C \cdot \frac{q_1 q_2 exp(-\kappa r)}{r}\quad \mathrm{for}\quad r<r_{\mathrm{cut}}
 
+where :math:`C` is defined as in Eqn. :eq:`coulomb_prefactor`.
 The Debye-Hückel potential is an approximate method for calculating
 electrostatic interactions, but technically it is treated as other
 short-ranged non-bonding potentials. For :math:`r>r_{\mathrm cut}` it is
@@ -1501,9 +1505,9 @@ introduces three new parameters :math:`\varepsilon_\mathrm{int}`,
 
    U(r)^{C-DHC} = 
      \begin{cases} 
-       \frac{l_B k_B T q_1 q_2}{\varepsilon_{\text{int}} r} & \text{if } r < r_0, \\ 
-       \frac{l_B k_B T q_1 q_2 e^{-\alpha (r - r_0)}}{\varepsilon_{\text{int}} r} & \text{if } r_0 < r < r_1,  \\
-       \frac{l_B k_B T q_1 q_2 e^{-\kappa r}}{\varepsilon_{\text{ext}} r} & \text{if } r_{\text{cut}} > r > r_1,  \\
+       \frac{C q_1 q_2}{\varepsilon_{\text{int}} r} & \text{if } r < r_0, \\ 
+       \frac{C q_1 q_2 e^{-\alpha (r - r_0)}}{\varepsilon_{\text{int}} r} & \text{if } r_0 < r < r_1,  \\
+       \frac{C q_1 q_2 e^{-\kappa r}}{\varepsilon_{\text{ext}} r} & \text{if } r_{\text{cut}} > r > r_1,  \\
        0 & \text{if } r > r_{\text{cut}}.
      \end{cases}
 
@@ -1538,9 +1542,10 @@ parameters.
 
 ::
 
-    mmm1d = MMM1D(bjerrum_length=lb, far_switch_radius = fr, maxPWerror=err, tune=False, bessel_cutoff=bc)
-    mmm1d = MMM1D(bjerrum_length=lb, maxPWerror=err)
+    mmm1d = MMM1D(prefactor=C, far_switch_radius = fr, maxPWerror=err, tune=False, bessel_cutoff=bc)
+    mmm1d = MMM1D(prefactor=C, maxPWerror=err)
 
+where the prefactor :math:`C` is defined in Eqn. :eq:`coulomb_prefactor`.
 MMM1D coulomb method for systems with periodicity 0 0 1. Needs the
 nsquared cell system (see section :ref:`cellsystem`). The first form sets parameters
 manually. The switch radius determines at which xy-distance the force
@@ -1554,8 +1559,8 @@ test force calculations.
 
 ::
 
-    mmm1d_gpu = MMM1DGPU(bjerrum_length=lb, far_switch_radius = fr, maxPWerror=err, tune=False, bessel_cutoff=bc)
-    mmm1d_gpu = MMM1DGPU(bjerrum_length=lb, maxPWerror=err)
+    mmm1d_gpu = MMM1DGPU(prefactor=C, far_switch_radius = fr, maxPWerror=err, tune=False, bessel_cutoff=bc)
+    mmm1d_gpu = MMM1DGPU(prefactor=C, maxPWerror=err)
 
 MMM1D is also available in a GPU implementation. Unlike its CPU
 counterpart, it does not need the nsquared cell system. The first form
@@ -1587,10 +1592,10 @@ the given bound. Note thate the user has to take care that the particles don't
 leave the box in the nonperiodic z-direction e.g. with constraints. By default,
 no dielectric contrast is set and it is used as::
 
-	mmm2d = electrostatics.MMM2D(bjerrum_length = 1.0, maxPWerror = 1e-3)
+	mmm2d = electrostatics.MMM2D(prefactor=C, maxPWerror = 1e-3)
 	system.actors.add(mmm2d)
 
-
+where the prefactor :math:`C` is defined in Eqn. :eq:`coulomb_prefactor`.
 For a detailed list of parameters see :attr:`espressomd.electrostatics.MMM2D`. 
 The last two, mutually exclusive parameters `dielectric` and
 `dielectric_constants_on` allow to specify dielectric contrasts at the
@@ -1600,7 +1605,7 @@ however is only used to calculate the contrasts. That is, specifying
 :math:`\epsilon_t=\epsilon_m=\epsilon_b=\text{const}` is always
 identical to :math:`\epsilon_t=\epsilon_m=\epsilon_b=1`::
 
-	mmm2d = electrostatics.MMM2D(bjerrum_length = 1.0, maxPWerror = 1e-3, dielectric = 1, top = 1, mid = 1, bot = 1)
+	mmm2d = electrostatics.MMM2D(prefactor = C, maxPWerror = 1e-3, dielectric = 1, top = 1, mid = 1, bot = 1)
 
 The second form specifies only the dielectric contrasts at the boundaries,
 that is :math:`\Delta_t=\frac{\epsilon_m-\epsilon_t}{\epsilon_m+\epsilon_t}`
@@ -1608,15 +1613,15 @@ and :math:`\Delta_b=\frac{\epsilon_m-\epsilon_b}{\epsilon_m+\epsilon_b}`.
 Using this form allows to choose :math:`\Delta_{t/b}=-1`, corresponding
 to metallic boundary conditions::
 
-	mmm2d = electrostatics.MMM2D(bjerrum_length = 1.0, maxPWerror = 1e-3, dielectric_contrast_on = 1, delta_mid_top = -1, delta_mid_bot = -1)
+	mmm2d = electrostatics.MMM2D(prefactor = C, maxPWerror = 1e-3, dielectric_contrast_on = 1, delta_mid_top = -1, delta_mid_bot = -1)
 
 Using `const_pot` allows to maintain a constant electric potential difference `pot_diff`
 between the xy-planes at :math:`z=0` and :math:`z=L`, where :math:`L`
 denotes the box length in :math:`z`-direction::
 	
-	mmm2d = electrostatics.MMM2D(bjerrum_length = 100.0, maxPWerror = 1e-3, const_pot = 1, pot_diff = 100.0)
+	mmm2d = electrostatics.MMM2D(prefactor = 100.0, maxPWerror = 1e-3, const_pot = 1, pot_diff = 100.0)
 
-This is done by countering the total dipol moment of the system with the
+This is done by countering the total dipole moment of the system with the
 electric field :math:`E_{induced}` and superposing a homogeneous electric field
 :math:`E_{applied} = \frac{U}{L}` to retain :math:`U`. This mimics the
 induction of surface charges :math:`\pm\sigma = E_{induced} \cdot \epsilon_0`
@@ -2013,7 +2018,7 @@ fast Coulomb solvers* library. The specific methods available depend on the comp
 To use SCAFACOS, create an instance of :attr:`espressomd.electrostatics.Scafacos` and add it to the list of active actors. Three parameters have to be specified:
 * method_name: name of the SCAFACOS method being used.
 * method_params: dictionary containing the method-specific parameters
-* bjerrum_length
+* prefactor: Coulomb prefactor as defined in :eq:`coulomb_prefactor`.
 The method-specific parameters are described in the SCAFACOS manual.
 Additionally, methods supporting tuning have the parameter ``tolerance_field`` which sets the desired root mean square accuracy for the electric field 
 
@@ -2022,7 +2027,7 @@ cutoff to :math:`1.5` and tune the other parameters for an accuracy of
 :math:`10^{-3}`, use::
 
   from espressomd.electrostatics import Scafacos
-  scafacos=Scafacos(bjerrum_length=1,method_name="ewald", 
+  scafacos=Scafacos(prefactor=1,method_name="ewald", 
     method_params={"ewald_r_cut":1.5, "tolerance_field":1e-3})
   system.actors.add(scafacos)
   
@@ -2040,14 +2045,15 @@ is defined as follows:
 
 .. math::
 
-   U^{D-P3M}(\vec{r}) = l_{B} k_B T \left( \frac{(\vec{\mu}_i \cdot \vec{\mu}_j)}{r^3} 
+   U^{D-P3M}(\vec{r}) = D \cdot \left( \frac{(\vec{\mu}_i \cdot \vec{\mu}_j)}{r^3} 
      - \frac{3  (\vec{\mu}_i \cdot \vec{r})  (\vec{\mu}_j \cdot \vec{r}) }{r^5} \right)
 
 where :math:`r=|\vec{r}|`.
+The prefactor :math:`D` is can be set by the user and is given by
+
+.. math::
 
-:math:`l_{B}` is a dimensionless parameter similar to the Bjerrum length
-in electrostatics which helps to tune the effect of the medium on the
-magnetic interaction between two magnetic dipoles.
+  D =\frac{\mu_0 \mu}{4\pi}
 
 Computing magnetostatic interactions is computationally very expensive.
 features some state-of-the-art algorithms to deal with these
diff --git a/doc/sphinx/io.rst b/doc/sphinx/io.rst
index d8ef8fb06..d8f7ee94b 100644
--- a/doc/sphinx/io.rst
+++ b/doc/sphinx/io.rst
@@ -109,7 +109,7 @@ To give an example::
     checkpoint.register("system.part")
 
     # ... set charges of particles here ... from espressomd import
-    electrostatics p3m = electrostatics.P3M(bjerrum_length=1.0, accuracy=1e-2)
+    electrostatics p3m = electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
     system.actors.add(p3m)
     checkpoint.register("p3m")
 
diff --git a/samples/billard.py b/samples/billard.py
index a6b9d13ba..18950fced 100644
--- a/samples/billard.py
+++ b/samples/billard.py
@@ -159,7 +159,7 @@ def main():
     ball.fix = [1,1,1]
     system.thermostat.set_langevin(kT=0, gamma=0.8)
     #ELECTROSTATICS
-    #	p3m = electrostatics.P3M(bjerrum_length=50, accuracy=1e-2)
+    #	p3m = electrostatics.P3M(prefactor=50, accuracy=1e-2)
     #	system.actors.add(p3m)
 
     cleared_balls= [0,0]
diff --git a/samples/coulomb_debye_hueckel.py b/samples/coulomb_debye_hueckel.py
index cbf052951..59bbead78 100644
--- a/samples/coulomb_debye_hueckel.py
+++ b/samples/coulomb_debye_hueckel.py
@@ -113,14 +113,14 @@ for i in range(n_part / 2):
 #    print("Particle {} has charge {} and is of type {}.".format(i,system.part[i].q,system.part[i].type))
 
 # Activating the Debye-Hueckel interaction
-# The Bjerrum length is set to one. Assuming the solvent is water, this
+# The Coulomb prefactor is set to one. Assuming the solvent is water, this
 # means that lj_sig is 0.714 nm in SI units.
-l_B = 1
+coulomb_prefactor = 1
 # inverse Debye length for 1:1 electrolyte in water at room temperature (nm)
 dh_kappa = numpy.sqrt(mol_dens) / 0.304
 # convert to MD units
 dh_kappa = dh_kappa / 0.714
-dh = electrostatics.CDH(bjerrum_length=l_B, kappa=dh_kappa,
+dh = electrostatics.CDH(prefactor=coulomb_prefactor, kappa=dh_kappa,
                         r_cut=5 / dh_kappa, eps_int=1, eps_ext=1, r0=0.5, r1=1, alpha=1)
 system.actors.add(dh)
 print(system.actors)
diff --git a/samples/debye_hueckel.py b/samples/debye_hueckel.py
index 9755db68f..8fd5ac7f2 100644
--- a/samples/debye_hueckel.py
+++ b/samples/debye_hueckel.py
@@ -115,15 +115,15 @@ for i in range(n_part / 2):
 #    print("Particle {} has charge {} and is of type {}.".format(i,system.part[i].q,system.part[i].type))
 
 # Activating the Debye-Hueckel interaction
-# The Bjerrum length is set to one. Assuming the solvent is water, this
+# The Coulomb prefactor is set to one. Assuming the solvent is water, this
 # means that lj_sig is 0.714 nm in SI units.
-l_B = 1
+coulomb_prefactor =1
 # inverse Debye length for 1:1 electrolyte in water at room temperature (nm)
 dh_kappa = numpy.sqrt(mol_dens) / 0.304
 # convert to MD units
 dh_kappa = dh_kappa / 0.714
 dh = electrostatics.DH(
-    bjerrum_length=l_B, kappa=dh_kappa, r_cut=int(5 / dh_kappa))
+    prefactor=coulomb_prefactor, kappa=dh_kappa, r_cut=int(5 / dh_kappa))
 system.actors.add(dh)
 print(system.actors)
 
diff --git a/samples/ekboundaries.py b/samples/ekboundaries.py
index e926b2ce3..fde3684fd 100644
--- a/samples/ekboundaries.py
+++ b/samples/ekboundaries.py
@@ -7,7 +7,7 @@ system.time_step = 0.1
 
 
 ek = electrokinetics.Electrokinetics(
-    lb_density=1, friction=1, agrid=1, viscosity=1, T=1, bjerrum_length=1)
+    lb_density=1, friction=1, agrid=1, viscosity=1, T=1, prefactor=1)
 
 pos = electrokinetics.Species(
     density=0.05, D=0.1, valency=1, ext_force=[0, 0, 1.])
diff --git a/samples/electrophoresis.py b/samples/electrophoresis.py
index 2d328cbfa..0f72359de 100644
--- a/samples/electrophoresis.py
+++ b/samples/electrophoresis.py
@@ -156,7 +156,7 @@ if os.path.isfile("p3m_checkpoint") and read_checkpoint == True:
     print("reading p3m from file")
     p3m = pickle.load(open("p3m_checkpoint", "r"))
 else:
-    p3m = electrostatics.P3M(bjerrum_length=1.0, accuracy=1e-2)
+    p3m = electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
     print("Tuning P3M")
 
 system.actors.add(p3m)
diff --git a/samples/minimal-charged-particles.py b/samples/minimal-charged-particles.py
index 856d18318..4abde4b7f 100644
--- a/samples/minimal-charged-particles.py
+++ b/samples/minimal-charged-particles.py
@@ -100,7 +100,7 @@ system.force_cap = lj_cap
 
 # P3M setup after charge assigned
 #############################################################
-p3m = electrostatics.P3M(bjerrum_length=1.0, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
 system.actors.add(p3m)
 
 
diff --git a/samples/p3m.py b/samples/p3m.py
index cf62db2c5..29fc838bc 100644
--- a/samples/p3m.py
+++ b/samples/p3m.py
@@ -114,8 +114,8 @@ for i in range(n_part / 2 - 1):
 #############################################################
 
 print("\nSCRIPT--->Create p3m\n")
-#p3m = electrostatics.P3M_GPU(bjerrum_length=2.0, accuracy=1e-2)
-p3m = electrostatics.P3M(bjerrum_length=2.0, accuracy=1e-2)
+#p3m = electrostatics.P3M_GPU(prefactor=2.0, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=2.0, accuracy=1e-2)
 
 print("\nSCRIPT--->Add actor\n")
 system.actors.add(p3m)
diff --git a/samples/save_checkpoint.py b/samples/save_checkpoint.py
index 999944248..65d3a8f62 100644
--- a/samples/save_checkpoint.py
+++ b/samples/save_checkpoint.py
@@ -59,7 +59,7 @@ checkpoint.register("system.part")
 for i in range(n_part / 2 - 1):
     system.part[2 * i].q = -1.0
     system.part[2 * i + 1].q = 1.0
-p3m = electrostatics.P3M(bjerrum_length=1.0, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
 system.actors.add(p3m)
 
 checkpoint.register("p3m")
diff --git a/samples/store_properties.py b/samples/store_properties.py
index 8ff885c26..37801da2e 100644
--- a/samples/store_properties.py
+++ b/samples/store_properties.py
@@ -111,7 +111,7 @@ for i in range(n_part / 2 - 1):
 
 # P3M setup after charge assigned
 #############################################################
-p3m = electrostatics.P3M(bjerrum_length=1.0, accuracy=1e-2)
+p3m = electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
 system.actors.add(p3m)
 
 
diff --git a/samples/visualization_mmm2d.py b/samples/visualization_mmm2d.py
index 2abbf913a..622c103f9 100755
--- a/samples/visualization_mmm2d.py
+++ b/samples/visualization_mmm2d.py
@@ -40,7 +40,7 @@ print("After Minimization: E_total=", energy['total'])
 
 system.thermostat.set_langevin(kT=0.1, gamma=1.0)
 
-mmm2d = electrostatics.MMM2D(bjerrum_length = 100.0, maxPWerror = 1e-3, const_pot = 1, pot_diff = 50.0)
+mmm2d = electrostatics.MMM2D(prefactor = 10.0, maxPWerror = 1e-3, const_pot = 1, pot_diff = 50.0)
 system.actors.add(mmm2d)
 
 def main():
diff --git a/src/core/elc.hpp b/src/core/elc.hpp
index bfcc59401..1515563a4 100755
--- a/src/core/elc.hpp
+++ b/src/core/elc.hpp
@@ -84,7 +84,7 @@ extern ELC_struct elc_params;
 /** set parameters for ELC.
     @param maxPWerror the required accuracy of the potential and the force. Note
    that this counts for the
-    plain 1/r contribution alone, without the Bjerrum length and the charge
+    plain 1/r contribution alone, without the prefactor and the charge
    prefactor.
     @param min_dist   the gap size.
     @param far_cut    the cutoff of the exponential sum. If -1, the cutoff is
diff --git a/src/core/electrokinetics.hpp b/src/core/electrokinetics.hpp
index bce36b6b2..342feb30a 100755
--- a/src/core/electrokinetics.hpp
+++ b/src/core/electrokinetics.hpp
@@ -60,7 +60,7 @@ typedef struct {
   float gamma_even;
   float friction;
   float T;
-  float bjerrumlength;
+  float prefactor;
   float lb_force[3];
   unsigned int number_of_species;
   int reaction_species[3];
@@ -164,7 +164,7 @@ int ek_set_lb_density(double lb_density);
 int ek_set_viscosity(double viscosity);
 int ek_set_friction(double friction);
 int ek_set_T(double T);
-int ek_set_bjerrumlength(double bjerrumlength);
+int ek_set_prefactor(double prefactor);
 #ifdef EK_ELECTROSTATIC_COUPLING
 int ek_set_electrostatics_coupling( bool electrostatics_coupling );
 void ek_calculate_electrostatic_coupling();
diff --git a/src/core/electrokinetics_cuda.cu b/src/core/electrokinetics_cuda.cu
index 4d8ae594a..33e0298c2 100755
--- a/src/core/electrokinetics_cuda.cu
+++ b/src/core/electrokinetics_cuda.cu
@@ -2440,10 +2440,10 @@ int ek_init() {
   if( ek_parameters.agrid < 0.0 ||
       ek_parameters.viscosity < 0.0 ||
       ek_parameters.T < 0.0 ||
-      ek_parameters.bjerrumlength < 0.0 ) 
+      ek_parameters.prefactor < 0.0 ) 
   {
       
-    fprintf( stderr, "ERROR: invalid agrid, viscosity, T or bjerrum_length\n" );
+    fprintf( stderr, "ERROR: invalid agrid, viscosity, T or prefactor\n" );
     
     return 1;
   }
@@ -2557,7 +2557,7 @@ int ek_init() {
     if(electrostatics != nullptr)
       delete electrostatics;
 
-    FdElectrostatics::InputParameters es_parameters = {ek_parameters.bjerrumlength, ek_parameters.T, int(ek_parameters.dim_x), int(ek_parameters.dim_y), int(ek_parameters.dim_z), ek_parameters.agrid};
+    FdElectrostatics::InputParameters es_parameters = {ek_parameters.prefactor, int(ek_parameters.dim_x), int(ek_parameters.dim_y), int(ek_parameters.dim_z), ek_parameters.agrid};
     try {
       electrostatics = new FdElectrostatics(es_parameters, stream[0]);
     }
@@ -3362,7 +3362,7 @@ void ek_print_parameters() {
   printf( "  float gamma_even = %f;\n",                 ek_parameters.gamma_even );
   printf( "  float friction = %f;\n",                   ek_parameters.friction );
   printf( "  float T = %f;\n",                          ek_parameters.T );
-  printf( "  float bjerrumlength = %f;\n",              ek_parameters.bjerrumlength );
+  printf( "  float prefactor = %f;\n",              ek_parameters.prefactor );
   printf( "  float lb_force[] = {%f, %f, %f};\n",       ek_parameters.lb_force[0], 
                                                         ek_parameters.lb_force[1], 
                                                         ek_parameters.lb_force[2] );
@@ -3494,9 +3494,9 @@ int ek_set_lb_density( double lb_density ) {
 }
 
 
-int ek_set_bjerrumlength( double bjerrumlength ) {
+int ek_set_prefactor( double prefactor ) {
 
-  ek_parameters.bjerrumlength = bjerrumlength;
+  ek_parameters.prefactor = prefactor;
   return 0;
 }
 #ifdef EK_ELECTROSTATIC_COUPLING
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
index 430e7412c..892795222 100755
--- a/src/core/energy.cpp
+++ b/src/core/energy.cpp
@@ -75,7 +75,7 @@ void init_energies(Observable_stat *stat) {
 
   switch (coulomb.Dmethod) {
   case DIPOLAR_NONE:
-    n_dipolar = 1;
+    n_dipolar = 1; // because there may be an external magnetic field
     break;
   case DIPOLAR_MDLC_P3M:
     n_dipolar = 3;
diff --git a/src/core/fd-electrostatics.hpp b/src/core/fd-electrostatics.hpp
index 78e5df050..924fc9616 100644
--- a/src/core/fd-electrostatics.hpp
+++ b/src/core/fd-electrostatics.hpp
@@ -20,7 +20,7 @@ typedef void cufftReal;
 class FdElectrostatics {
   public:
     struct InputParameters {
-      float bjerrum_length, kT;
+      float prefactor;
       int dim_x, dim_y, dim_z;
       float agrid;
     };
diff --git a/src/core/fd-electrostatics_cuda.cu b/src/core/fd-electrostatics_cuda.cu
index 742b83162..80480598f 100644
--- a/src/core/fd-electrostatics_cuda.cu
+++ b/src/core/fd-electrostatics_cuda.cu
@@ -133,8 +133,8 @@ __global__ void createGreensfcn() {
       fde_parameters_gpu.greensfcn[index] = 0.0f;
     } else {
       fde_parameters_gpu.greensfcn[index] =
-          -4.0f * PI_FLOAT * fde_parameters_gpu.bjerrum_length *
-          fde_parameters_gpu.kT * fde_parameters_gpu.agrid *
+          -4.0f * PI_FLOAT * fde_parameters_gpu.prefactor *
+          fde_parameters_gpu.agrid *
           fde_parameters_gpu.agrid * 0.5f /
           (cos(2.0f * PI_FLOAT * coord[0] /
                (cufftReal)fde_parameters_gpu.dim_x) +
diff --git a/src/core/initialize.cpp b/src/core/initialize.cpp
index 788087735..1cc2f1831 100755
--- a/src/core/initialize.cpp
+++ b/src/core/initialize.cpp
@@ -314,8 +314,6 @@ void on_coulomb_change() {
   EVENT_TRACE(fprintf(stderr, "%d: on_coulomb_change\n", this_node));
   invalidate_obs();
 
-  recalc_coulomb_prefactor();
-
 #ifdef ELECTROSTATICS
   switch (coulomb.method) {
   case COULOMB_DH:
@@ -561,9 +559,6 @@ void on_temperature_change() {
   }
 #endif
 
-#ifdef ELECTROSTATICS
-  recalc_coulomb_prefactor();
-#endif
 }
 
 void on_parameter_change(int field) {
diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp
index 186ad512e..de2ead758 100755
--- a/src/core/integrate.cpp
+++ b/src/core/integrate.cpp
@@ -1194,26 +1194,20 @@ int integrate_set_npt_isotropic(double ext_pressure, double piston, int xdir,
 
 /* Sanity Checks */
 #ifdef ELECTROSTATICS
-  if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.bjerrum > 0) {
-    fprintf(stderr, "WARNING: If electrostatics is being used you must use the "
-                    "-cubic_box option!\n");
-    fprintf(stderr,
-            "Automatically reverting to a cubic box for npt integration.\n");
-    fprintf(stderr, "Be aware though that all of the coulombic pressure is "
-                    "added to the x-direction only!\n");
-    nptiso.cubic_box = 1;
-  }
+  if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.prefactor > 0) {
+    runtimeErrorMsg() << "WARNING: If electrostatics is being used you must use the the cubic box npt." ;
+    integ_switch = INTEG_METHOD_NVT;
+    mpi_bcast_parameter(FIELD_INTEG_SWITCH);
+    return ES_ERROR;
+ }
 #endif
 
 #ifdef DIPOLES
-  if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.Dbjerrum > 0) {
-    fprintf(stderr, "WARNING: If magnetostatics is being used you must use the "
-                    "-cubic_box option!\n");
-    fprintf(stderr,
-            "Automatically reverting to a cubic box for npt integration.\n");
-    fprintf(stderr, "Be aware though that all of the magnetostatic pressure is "
-                    "added to the x-direction only!\n");
-    nptiso.cubic_box = 1;
+  if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.Dprefactor > 0) {
+    runtimeErrorMsg() << "WARNING: If magnetostatics is being used you must use the the cubic box npt." ;
+    integ_switch = INTEG_METHOD_NVT;
+    mpi_bcast_parameter(FIELD_INTEG_SWITCH);
+    return ES_ERROR;
   }
 #endif
 
diff --git a/src/core/interaction_data.cpp b/src/core/interaction_data.cpp
index fd70b6781..4467e1539 100755
--- a/src/core/interaction_data.cpp
+++ b/src/core/interaction_data.cpp
@@ -78,10 +78,10 @@ std::vector<IA_parameters> ia_params;
 #if defined(ELECTROSTATICS) || defined(DIPOLES)
 Coulomb_parameters coulomb = {
 #ifdef ELECTROSTATICS
-    0.0, 0.0, COULOMB_NONE,
+  0.0, COULOMB_NONE,
 #endif
 #ifdef DIPOLES
-    0.0, 0.0, DIPOLAR_NONE,
+  0.0, DIPOLAR_NONE,
 #endif
 };
 #endif
@@ -585,19 +585,30 @@ void set_dipolar_method_local(DipolarInteraction method) {
 /*                                 electrostatics */
 /********************************************************************************/
 
-int coulomb_set_bjerrum(double bjerrum) {
-  if (bjerrum < 0.0)
+int coulomb_set_prefactor(double prefactor)
+{
+  if (prefactor < 0.0) {
+    runtimeErrorMsg() << "Coulomb prefactor has to be >=0";
     return ES_ERROR;
+  }
+  
+  coulomb.prefactor=prefactor;
+  mpi_bcast_coulomb_params();
 
-  coulomb.bjerrum = bjerrum;
+ 
+  return ES_OK;
+}
 
-  if (coulomb.bjerrum == 0.0) {
-    switch (coulomb.method) {
+/** @brief Deactivates the current Coulomb mhthod 
+    This was part of coulomb_set_bjerrum()
+*/
+void deactivate_coulomb_method() {
+coulomb.prefactor =0;
+switch (coulomb.method) {
 #ifdef P3M
     case COULOMB_ELC_P3M:
     case COULOMB_P3M_GPU:
     case COULOMB_P3M:
-      p3m_set_bjerrum();
       break;
 #endif
     case COULOMB_DH:
@@ -619,64 +630,34 @@ int coulomb_set_bjerrum(double bjerrum) {
     mpi_bcast_coulomb_params();
     coulomb.method = COULOMB_NONE;
     mpi_bcast_coulomb_params();
-  }
-
-  return ES_OK;
 }
 
+
+
+
+
+
 /* =========================================================
    ========================================================= */
 #endif /*ifdef ELECTROSTATICS */
 
 #ifdef DIPOLES
 
-int dipolar_set_Dbjerrum(double bjerrum) {
-  if (bjerrum < 0.0)
+int dipolar_set_Dprefactor(double prefactor)
+{
+  if (prefactor < 0.0){
+    runtimeErrorMsg() << "Dipolar prefactor has to be >=0";
     return ES_ERROR;
-
-  coulomb.Dbjerrum = bjerrum;
-
-  if (coulomb.Dbjerrum == 0.0) {
-    switch (coulomb.Dmethod) {
-#ifdef DP3M
-    case DIPOLAR_MDLC_P3M:
-    // fall through
-    case DIPOLAR_P3M:
-      coulomb.Dbjerrum = bjerrum;
-      dp3m_set_bjerrum();
-      break;
-#endif
-    case DIPOLAR_SCAFACOS:;
-    // Fall through
-    default:
-      break;
-    }
-
-    mpi_bcast_coulomb_params();
-    set_dipolar_method_local(DIPOLAR_NONE);
-    mpi_bcast_coulomb_params();
   }
+  
+  coulomb.Dprefactor = prefactor;
 
+  mpi_bcast_coulomb_params();
   return ES_OK;
 }
 
 #endif /* ifdef  DIPOLES */
 
-void recalc_coulomb_prefactor() {
-#ifdef ELECTROSTATICS
-  if (temperature > 0.0)
-    coulomb.prefactor = coulomb.bjerrum * temperature;
-  else
-    coulomb.prefactor = coulomb.bjerrum;
-#endif
-
-#ifdef DIPOLES
-  if (temperature > 0.0)
-    coulomb.Dprefactor = coulomb.Dbjerrum * temperature;
-  else
-    coulomb.Dprefactor = coulomb.Dbjerrum;
-#endif
-}
 
 int virtual_set_params(int bond_type) {
   if (bond_type < 0)
diff --git a/src/core/interaction_data.hpp b/src/core/interaction_data.hpp
index f9556d17c..3bf4400eb 100755
--- a/src/core/interaction_data.hpp
+++ b/src/core/interaction_data.hpp
@@ -467,8 +467,6 @@ struct IA_parameters {
 typedef struct {
 
 #ifdef ELECTROSTATICS
-  /** Bjerrum length. */
-  double bjerrum;
   /** bjerrum length times temperature. */
   double prefactor;
 
@@ -477,7 +475,6 @@ typedef struct {
 #endif
 
 #ifdef DIPOLES
-  double Dbjerrum;
   double Dprefactor;
   DipolarInteraction Dmethod;
 #endif
@@ -852,11 +849,19 @@ extern int ia_excl;
  ************************************************/
 
 #ifdef ELECTROSTATICS
-int coulomb_set_bjerrum(double bjerrum);
+/** @brief Set the electrostatics prefactor */
+int coulomb_set_prefactor(double prefactor);
+
+
+/** @brief Deactivates the current Coulomb mhthod 
+    This was part of coulomb_set_bjerrum()
+*/
+void deactivate_coulomb_method();
 #endif
 
 #ifdef DIPOLES
-int dipolar_set_Dbjerrum(double bjerrum);
+/** @brief Set the dipolar prefactor */
+int dipolar_set_Dprefactor(double prefactor);
 #endif
 
 /** get interaction parameters between particle sorts i and j */
@@ -896,9 +901,6 @@ void realloc_ia_params(int nsize);
     value is used in the verlet pair list algorithm. */
 void recalc_maximal_cutoff();
 
-/** call when the temperature changes, for Bjerrum length adjusting. */
-void recalc_coulomb_prefactor();
-
 /** check whether all force calculation routines are properly initialized. */
 int interactions_sanity_checks();
 
diff --git a/src/core/maggs.cpp b/src/core/maggs.cpp
index df51a4878..ed37cc378 100755
--- a/src/core/maggs.cpp
+++ b/src/core/maggs.cpp
@@ -154,7 +154,7 @@ typedef struct {
 /*******************************/
 
 /* create system structure with all zeros. Filled in maggs_set_parameters(); */
-MAGGS_struct maggs = { 1, 0, 1.0, 1.0, 0. , 0. , 0. , 0. , 0. , 0 , 0. , 0. , {{0.},{0.}}};
+MAGGS_struct maggs = { 1, 0, 1.0, 1.0, 0. , 0. , 0. , 0. ,  0 , 0. , 0. , {{0.},{0.}}};
 
 /* local mesh. */
 static lattice_parameters lparams;
@@ -195,7 +195,7 @@ static t_dirs* neighbor;
 // void maggs_update_plaquette(int mue, int nue, int* Neighbor, int index, double delta); /* update fields on plaquette */
 
 /****** setup everything: ******/
-// int maggs_set_parameters(double bjerrum, double f_mass, int mesh); /* set the system parameters */
+// int maggs_set_parameters(double prefactor, double f_mass, int mesh); /* set the system parameters */
 // void maggs_setup_neighbors(); /* setup the site neighbors */
 // void maggs_setup_local_lattice(); /* set lattice parameters */
 // void maggs_calc_surface_patches(t_surf_patch* surface_patch); /* prepare for communication */
@@ -393,8 +393,8 @@ int maggs_sanity_checks()
 
   FOR3D(d) if(node_grid[d] > max_node_grid) max_node_grid = node_grid[d];
 	
-  if (maggs.bjerrum == 0.) {
-      runtimeErrorMsg() <<"MEMD: bjerrum length is zero.";
+  if (maggs.prefactor == 0.) {
+      runtimeErrorMsg() <<"MEMD: prefactor is zero.";
     ret = -1;
   }
   else if ( (box_l[0] != box_l[1]) || (box_l[1] != box_l[2]) ) {
@@ -443,9 +443,8 @@ int maggs_sanity_checks()
 
 void maggs_compute_dipole_correction() {
   /*
-  maggs.prefactor  = sqrt(4. * M_PI * maggs.bjerrum * temperature);
+  maggs.pref1  = sqrt(4. * M_PI * maggs.prefactor);
                    = sqrt(1/epsilon)
-  maggs.pref2      = maggs.bjerrum * temperature;
                    = 1/4*pi*epsilon
   */
 
@@ -476,7 +475,7 @@ void maggs_compute_dipole_correction() {
 /****** setup everything: ******/
 /*******************************/
 
-int maggs_set_parameters(double bjerrum, double f_mass, int mesh, int finite_epsilon_flag, double epsilon_infty)
+int maggs_set_parameters(double prefactor, double f_mass, int mesh, int finite_epsilon_flag, double epsilon_infty)
 {
   if (f_mass <=0.) {
     return -1;
@@ -488,7 +487,7 @@ int maggs_set_parameters(double bjerrum, double f_mass, int mesh, int finite_eps
   maggs.mesh           = mesh;
   maggs.finite_epsilon_flag  = finite_epsilon_flag;
   maggs.epsilon_infty  = epsilon_infty;
-  maggs.bjerrum        = bjerrum;
+  maggs.prefactor        = prefactor;
   maggs.f_mass         = f_mass; 
   maggs.invsqrt_f_mass = 1./sqrt(f_mass); 	
 	
@@ -733,7 +732,7 @@ int maggs_get_mesh_1D()
  @param node_y              index of the node in y direction
  @param node_z              index of the node in z direction
  @param direction           direction in which the link points from the node. 0 is for x, 1 is for y, 2 is for z
- @param relative_epsilon    permittivity to set, relative to the background permittivity set by the bjerrum length
+ @param relative_epsilon    permittivity to set, relative to the background permittivity set by via the prefactor
  */
 double maggs_set_permittivity(int node_x, int node_y, int node_z, int direction, double relative_epsilon)
 {
@@ -1188,7 +1187,7 @@ void maggs_calc_self_energy_coeffs()
 
   factor = M_PI / maggs.mesh;
   prefac = 1. / maggs.mesh;
-  prefac = 0.5 * prefac * prefac * prefac * SQR(maggs.prefactor);
+  prefac = 0.5 * prefac * prefac * prefac * SQR(maggs.pref1);
   
   for(i=0;i<8;i++) 
     {
@@ -1423,13 +1422,13 @@ void maggs_calc_init_e_field()
 
     MPI_Allreduce(&localqz, &qz, 1, MPI_DOUBLE, MPI_SUM, zplane);
     qz = qz/(maggs.mesh*maggs.mesh);
-    qplane = qz*maggs.prefactor*invasq;
+    qplane = qz*maggs.pref1*invasq;
     /*    if(fabs(qplane) >= 0.01*ROUND_ERROR_PREC) { */
     for(iy=lparams.inner_left_down[1];iy<lparams.inner_up_right[1];iy++) {
       for(ix=lparams.inner_left_down[0];ix<lparams.inner_up_right[0];ix++) {  
 	index = maggs_get_linear_index(ix, iy, iz, lparams.dim);
 	Dfield[3*index+ZPLUS]  = Dfield[3*neighbor[index][ZMINUS]+ZPLUS] + qplane;
-	/*	    + qz*maggs.prefactor*invasq; */
+	/*	    + qz*maggs.pref1*invasq; */
       }
     }
     /*    } */
@@ -1466,12 +1465,12 @@ void maggs_calc_init_e_field()
       MPI_Allreduce(&localqy, &qy, 1, MPI_DOUBLE, MPI_SUM, yline);
 			
       qy = qy/maggs.mesh;
-      qline = (qy-qz)*maggs.prefactor*invasq;
+      qline = (qy-qz)*maggs.pref1*invasq;
       /*      if(fabs(qy-qz)>=ROUND_ERROR_PREC) { */
       for(ix=lparams.inner_left_down[0];ix<lparams.inner_up_right[0];ix++) {  
 	index = maggs_get_linear_index(ix, iy, iz, lparams.dim);
 	Dfield[3*index+YPLUS]  = Dfield[3*neighbor[index][YMINUS]+YPLUS] + qline;
-	/*	    (qy-qz)*maggs.prefactor*invasq; */
+	/*	    (qy-qz)*maggs.pref1*invasq; */
       }
       /*      } */
 			
@@ -1497,7 +1496,7 @@ void maggs_calc_init_e_field()
 	index = maggs_get_linear_index(ix, iy, iz, lparams.dim);
 	qx = lattice[index].charge / lattice[index].permittivity[0]; 
 	Dfield[3*index+XPLUS] = Dfield[3*neighbor[index][XMINUS]+XPLUS] + 
-	  (qx-qy)*maggs.prefactor*invasq;
+	  (qx-qy)*maggs.pref1*invasq;
       }
 			
       if(ix>=lparams.inner_up_right[0]-1) {
@@ -1675,7 +1674,7 @@ void maggs_calc_charge_fluxes_1D(double q, double *help, double *flux, int dir)
   int l,m; 
   double q_scaled;
 	
-  q_scaled = q * maggs.prefactor*maggs.inva;
+  q_scaled = q * maggs.pref1*maggs.inva;
   index = 0;
 	
   maggs_calc_directions(dir, &dir1, &dir2);   
@@ -2137,7 +2136,7 @@ void maggs_calc_self_influence(Particle* P)
   FOR3D(k) {
     local_force[k] = pow(((0.5 - relative_position[k])*SELF_FACTOR_1), 3.0) +
       (0.5 - relative_position[k]) * SELF_FACTOR_2;
-    local_force[k] *= maggs.bjerrum;
+    local_force[k] *= maggs.prefactor;
     local_force[k] *= particle_charge * particle_charge;
     local_force[k] *= invasq;
   }
@@ -2198,7 +2197,7 @@ void maggs_calc_part_link_forces(Particle *p, int index, double *grad)
   }  
   /* Attention! Here, the INTERLACING is done! */
   FOR3D(j){
-    p->f.f[j] += maggs.prefactor * local_force[j];
+    p->f.f[j] += maggs.pref1 * local_force[j];
   }
 }
 
@@ -2325,17 +2324,10 @@ void maggs_init()
 {
     maggs.inva  = (double) maggs.mesh/box_l[0]; 
     maggs.a     = 1.0/maggs.inva;
-    if (temperature>0.0) {
-      maggs.prefactor  = sqrt(4. * M_PI * maggs.bjerrum * temperature);
-      maggs.pref2      = maggs.bjerrum * temperature;
-    } else {
-      maggs.prefactor  = sqrt(4. * M_PI * maggs.bjerrum);
-      maggs.pref2      = maggs.bjerrum;
-    }			
+    maggs.pref1  = sqrt(4. * M_PI * maggs.prefactor);
     
     if (maggs_sanity_checks()) {
-      maggs.bjerrum = 0.0;
-      coulomb.bjerrum = 0.0;
+      runtimeErrorMsg() << "MEMD sanity checks failed.";
       return;
     }
     		
diff --git a/src/core/maggs.hpp b/src/core/maggs.hpp
index 26b6d8bac..77982cfc6 100755
--- a/src/core/maggs.hpp
+++ b/src/core/maggs.hpp
@@ -67,12 +67,9 @@ typedef struct {
   double f_mass;
   /** inverse of square root of f_mass. */
   double invsqrt_f_mass;
-  /** prefactor to convert field to force. */
   double prefactor;
-  /** prefactor / (4*pi) */
-  double pref2;
-  /** bjerrum length of the system. */
-  double bjerrum;
+  /** prefactor to convert field to force. */
+  double pref1;
   /** mesh size in one dimension */
   int    mesh;
   /** = 1/a = mesh / box_length */
@@ -97,13 +94,13 @@ extern MAGGS_struct maggs;
 void maggs_init(); /** called from: initialize.cpp */
 
 /** set the main parameters for the algorithm.
-    @param bjerrum   Bjerrum length for the system
+    @param prefactor   Electrostatics prefactor for the system
     @param f_mass    parameter to tune the speed of light (1/c^2)
     @param mesh      Mesh size in one dimension
     @param finite_epsilon_flag whether to do epsilon-at-infinity-correction
     @param epsilon_infty epsilon-at-infinity
  */
-int maggs_set_parameters(double bjerrum, double f_mass, int mesh, int finite_epsilon_flag, double epsilon_infty);
+int maggs_set_parameters(double prefactor, double f_mass, int mesh, int finite_epsilon_flag, double epsilon_infty);
 
 /** get lattice size in one dimension
  @return mesh in 1D
@@ -115,7 +112,7 @@ int maggs_get_mesh_1D();
  @param node_y              index of the node in y direction
  @param node_z              index of the node in z direction
  @param direction           direction in which the link points from the node. 0 is for x, 1 is for y, 2 is for z
- @param relative_epsilon    permittivity to set, relative to the background permittivity set by the bjerrum length
+ @param relative_epsilon    permittivity to set, relative to the background permittivity set by the electrostatics prefactor
  */
 double maggs_set_permittivity(int node_x, int node_y, int node_z, int direction, double relative_epsilon);
 
diff --git a/src/core/p3m-dipolar.cpp b/src/core/p3m-dipolar.cpp
index 3c4a748b1..0e973c929 100755
--- a/src/core/p3m-dipolar.cpp
+++ b/src/core/p3m-dipolar.cpp
@@ -310,7 +310,7 @@ void dp3m_pre_init(void) {
   dfft_pre_init();
 }
 
-void dp3m_set_bjerrum() {
+void dp3m_set_prefactor() {
   dp3m.params.alpha = 0.0;
   dp3m.params.alpha_L = 0.0;
   dp3m.params.r_cut = 0.0;
@@ -324,13 +324,12 @@ void dp3m_set_bjerrum() {
 void dp3m_init() {
   int n;
 
-  if (coulomb.Dbjerrum == 0.0) {
-    if (coulomb.Dbjerrum == 0.0) {
+  if (coulomb.Dprefactor <= 0.0) {
       dp3m.params.r_cut = 0.0;
       dp3m.params.r_cut_iL = 0.0;
-      if (this_node == 0)
+      if (this_node == 0) {
         P3M_TRACE(
-            fprintf(stderr, "0: dp3m_init: dipolar Bjerrum length is zero.\n");
+            fprintf(stderr, "0: dp3m_init: dipolar prefactor is zero.\n");
             fprintf(stderr, "   Magnetostatics of dipoles switched off!\n"));
     }
   } else {
@@ -527,8 +526,6 @@ void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha,
   if (n_interpol != -1)
     dp3m.params.inter = n_interpol;
 
-  coulomb.Dprefactor =
-      (temperature > 0) ? temperature * coulomb.Dbjerrum : coulomb.Dbjerrum;
 }
 
 /*****************************************************************************/
@@ -1576,9 +1573,6 @@ double dp3m_perform_aliasing_sums_energy(int n[3], double nominator[1]) {
 #define P3M_TUNE_MAX_CUTS 50
 /** Tune dipolar P3M parameters to desired accuracy.
 
-    Usage:
-    \verbatim inter dipolar <bjerrum> p3m tune accuracy <value> [r_cut <value>
-   mesh <value> cao <value>] \endverbatim
 
     The parameters are tuned to obtain the desired accuracy in best
     time, by running mpi_integrate(0) for several parameter sets.
@@ -1982,9 +1976,9 @@ int dp3m_adaptive_tune(char **logger) {
   mpi_bcast_event(P3M_COUNT_DIPOLES);
 
   /* Print Status */
-  sprintf(b, "Dipolar P3M tune parameters: Accuracy goal = %.5e Bjerrum Length "
+  sprintf(b, "Dipolar P3M tune parameters: Accuracy goal = %.5e prefactor "
              "= %.5e\n",
-          dp3m.params.accuracy, coulomb.Dbjerrum);
+          dp3m.params.accuracy, coulomb.Dprefactor);
   *logger = strcat_alloc(*logger, b);
   sprintf(b, "System: box_l = %.5e # charged part = %d Sum[q_i^2] = %.5e\n",
           box_l[0], dp3m.sum_dip_part, dp3m.sum_mu2);
@@ -2551,7 +2545,8 @@ void dp3m_calc_send_mesh() {
 /************************************************/
 
 void dp3m_scaleby_box_l() {
-  if (coulomb.Dbjerrum == 0.0) {
+  if (coulomb.Dprefactor < 0.0) {
+    runtimeErrorMsg() << "Dipolar prefactor has to be >=0" ;
     return;
   }
 
diff --git a/src/core/p3m-dipolar.hpp b/src/core/p3m-dipolar.hpp
index 713713d53..f6477d25e 100755
--- a/src/core/p3m-dipolar.hpp
+++ b/src/core/p3m-dipolar.hpp
@@ -121,7 +121,7 @@ int dp3m_set_eps(double eps);
  */
 void dp3m_init(void);
 
-void dp3m_set_bjerrum(void);
+void dp3m_set_prefactor(void);
 
 /** Updates \ref p3m_parameter_struct::alpha and \ref
  * p3m_parameter_struct::r_cut if \ref box_l changed. */
@@ -136,8 +136,8 @@ int dp3m_sanity_checks();
     Dcur_ca_frac. */
 void dp3m_dipole_assign(void);
 
-/** set bjerrum length for dipolar p3m */
-void dp3m_set_bjerrum(void);
+/** set prefactor for dipolar p3m */
+void dp3m_set_prefactor(void);
 
 int dp3m_adaptive_tune(char **log);
 
diff --git a/src/core/p3m.cpp b/src/core/p3m.cpp
index 5d0752552..dbbde3ed5 100755
--- a/src/core/p3m.cpp
+++ b/src/core/p3m.cpp
@@ -299,7 +299,7 @@ void p3m_free() {
     free(p3m.int_caf[i]);
 }
 
-void p3m_set_bjerrum() {
+void p3m_set_prefactor() {
   p3m.params.alpha = 0.0;
   p3m.params.alpha_L = 0.0;
   p3m.params.r_cut = 0.0;
@@ -311,12 +311,12 @@ void p3m_set_bjerrum() {
 }
 
 void p3m_init() {
-  if (coulomb.bjerrum == 0.0) {
+  if (coulomb.prefactor <= 0.0) {
     p3m.params.r_cut = 0.0;
     p3m.params.r_cut_iL = 0.0;
 
     if (this_node == 0)
-      P3M_TRACE(fprintf(stderr, "0: P3M_init: Bjerrum length is "
+      P3M_TRACE(fprintf(stderr, "0: P3M_init: prefactor is "
                                 "zero.\nElectrostatics switched off!\n"););
 
   } else {
@@ -412,8 +412,6 @@ void p3m_set_tune_params(double r_cut, int mesh[3], int cao, double alpha,
   if (n_interpol != -1)
     p3m.params.inter = n_interpol;
 
-  coulomb.prefactor =
-      (temperature > 0) ? temperature * coulomb.bjerrum : coulomb.bjerrum;
 }
 
 /*@}*/
@@ -1748,8 +1746,8 @@ int p3m_adaptive_tune(char **log) {
   mpi_bcast_event(P3M_COUNT_CHARGES);
   /* Print Status */
   sprintf(b,
-          "P3M tune parameters: Accuracy goal = %.5e Bjerrum Length = %.5e \n",
-          p3m.params.accuracy, coulomb.bjerrum);
+          "P3M tune parameters: Accuracy goal = %.5e prefactor = %.5e \n",
+          p3m.params.accuracy, coulomb.prefactor);
   *log = strcat_alloc(*log, b);
   sprintf(b, "System: box_l = %.5e # charged part = %d Sum[q_i^2] = %.5e\n",
           box_l[0], p3m.sum_qpart, p3m.sum_q2);
@@ -2296,7 +2294,8 @@ void p3m_calc_send_mesh() {
 /************************************************/
 
 void p3m_scaleby_box_l() {
-  if (coulomb.bjerrum == 0.0) {
+  if (coulomb.prefactor < 0.0) {
+    runtimeErrorMsg() << "The Coulomb prefactor has to be >=0";
     return;
   }
 
diff --git a/src/core/p3m.hpp b/src/core/p3m.hpp
index 401bac2fe..ba3d12779 100755
--- a/src/core/p3m.hpp
+++ b/src/core/p3m.hpp
@@ -125,8 +125,6 @@ extern p3m_data_struct p3m;
 
 void p3m_pre_init(void);
 
-void p3m_set_bjerrum(void);
-
 int p3m_adaptive_tune(char **log);
 
 /** Initialize all structures, parameters and arrays needed for the 
@@ -162,9 +160,6 @@ enum P3M_TUNE_ERROR { P3M_TUNE_FAIL = 1, P3M_TUNE_NOCUTOFF = 2, P3M_TUNE_CAOTOLA
 
 /** Tune P3M parameters to desired accuracy.
 
-    Usage:
-    \verbatim inter coulomb <bjerrum> p3m tune accuracy <value> [r_cut <value> mesh <value> cao <value>] \endverbatim
-
     The parameters are tuned to obtain the desired accuracy in best
     time, by running mpi_integrate(0) for several parameter sets.
 
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index 9f0b4ac2d..ae16c40d6 100755
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1192,7 +1192,7 @@ void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded,
 
   // Allocate mem for the double list
   stat->data.resize(total);
-
+  
   // Number of doubles per interaction (pressure=1, stress tensor=9,...)
   stat->chunk_size = c_size;
 
diff --git a/src/core/unit_tests/List_test.cpp b/src/core/unit_tests/List_test.cpp
index 50f9e03a2..0dacc38f3 100644
--- a/src/core/unit_tests/List_test.cpp
+++ b/src/core/unit_tests/List_test.cpp
@@ -240,6 +240,12 @@ BOOST_AUTO_TEST_CASE(resize) {
     l.resize(127);
     BOOST_CHECK(l.size() == 127);
     BOOST_CHECK(l.capacity() == 127);
+    l.resize(42);
+    BOOST_CHECK(l.size() == 42);
+    BOOST_CHECK(l.capacity() == 42);
+    l.resize(0);
+    BOOST_CHECK(l.size() == 0);
+    BOOST_CHECK(l.capacity() == 0);
   }
 
   /* List::resize() with size > capacity */
diff --git a/src/python/espressomd/electrokinetics.pxd b/src/python/espressomd/electrokinetics.pxd
index cccb69e28..38719ccd7 100644
--- a/src/python/espressomd/electrokinetics.pxd
+++ b/src/python/espressomd/electrokinetics.pxd
@@ -26,7 +26,7 @@ IF ELECTROKINETICS and CUDA:
             float gamma_even
             float friction
             float T
-            float bjerrumlength
+            float prefactor 
             float lb_force[3]
             unsigned int number_of_species
             int reaction_species[3]
@@ -80,7 +80,7 @@ IF ELECTROKINETICS and CUDA:
         int ek_set_viscosity(double viscosity)
         int ek_set_friction(double friction)
         int ek_set_T(double T)
-        int ek_set_bjerrumlength(double bjerrumlength)
+        int ek_set_prefactor(double prefactor)
         int ek_set_bulk_viscosity(double bulk_viscosity)
         int ek_set_gamma_odd(double gamma_odd)
         int ek_set_gamma_even(double gamma_even)
diff --git a/src/python/espressomd/electrokinetics.pyx b/src/python/espressomd/electrokinetics.pyx
index 3bed51395..09b3672ed 100644
--- a/src/python/espressomd/electrokinetics.pyx
+++ b/src/python/espressomd/electrokinetics.pyx
@@ -29,10 +29,10 @@ IF ELECTROKINETICS:
 
 
         def valid_keys(self):
-            return "agrid", "lb_density", "viscosity", "friction", "bulk_viscosity", "gamma_even", "gamma_odd", "T", "bjerrum_length", "stencil", "advection", "fluid_coupling"
+            return "agrid", "lb_density", "viscosity", "friction", "bulk_viscosity", "gamma_even", "gamma_odd", "T", "prefactor", "stencil", "advection", "fluid_coupling"
 
         def required_keys(self):
-            return ["agrid", "lb_density", "viscosity", "friction", "T", "bjerrum_length"]
+            return ["agrid", "lb_density", "viscosity", "friction", "T", "prefactor"]
 
         def default_params(self):
             return {"agrid": -1,
@@ -43,7 +43,7 @@ IF ELECTROKINETICS:
                     "gamma_even": 0.0,
                     "friction": 0.0,
                     "T": -1,
-                    "bjerrum_length": -1,
+                    "prefactor": -1,
                     "stencil": "linkcentered",
                     "advection": True,
                     "fluid_coupling": "friction"}
@@ -71,7 +71,7 @@ IF ELECTROKINETICS:
                     "gamma_even": ek_parameters.gamma_even,
                     "friction": ek_parameters.friction,
                     "T": ek_parameters.T,
-                    "bjerrum_length":ek_parameters.bjerrumlength,
+                    "prefactor":ek_parameters.prefactor,
                     "stencil": stencil,
                     "advection": ek_parameters.advection,
                     "fluid_coupling": fluid_coupling}
@@ -95,7 +95,7 @@ IF ELECTROKINETICS:
             ek_set_viscosity(self._params["viscosity"])
             ek_set_friction(self._params["friction"])
             ek_set_T(self._params["T"])
-            ek_set_bjerrumlength(self._params["bjerrum_length"])
+            ek_set_prefactor(self._params["prefactor"])
             ek_set_bulk_viscosity(self._params["bulk_viscosity"])
             ek_set_gamma_odd(self._params["gamma_odd"])
             ek_set_gamma_even(self._params["gamma_even"])
diff --git a/src/python/espressomd/electrostatics.pxd b/src/python/espressomd/electrostatics.pxd
index e65b75879..b9c564429 100644
--- a/src/python/espressomd/electrostatics.pxd
+++ b/src/python/espressomd/electrostatics.pxd
@@ -59,10 +59,10 @@ IF ELECTROSTATICS:
                 COULOMB_EK, \
                 COULOMB_SCAFACOS
 
-        int coulomb_set_bjerrum(double bjerrum)
+        int coulomb_set_prefactor(double prefactor)
+        void deactivate_coulomb_method()
 
         ctypedef struct Coulomb_parameters:
-            double bjerrum
             double prefactor
             CoulombMethod method
 
@@ -240,7 +240,7 @@ IF ELECTROSTATICS:
         int mmm1d_tune(char **log);
 
     cdef extern from "interaction_data.hpp":
-        int coulomb_set_bjerrum(double bjerrum)
+        int coulomb_set_prefactor(double prefactor)
 
         ctypedef enum CoulombMethod :
             COULOMB_NONE, 
@@ -258,7 +258,6 @@ IF ELECTROSTATICS:
             COULOMB_EK 
 
         ctypedef struct Coulomb_parameters:
-            double bjerrum
             double prefactor
             CoulombMethod method
 
diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx
index 8ce457176..86b95f8b5 100644
--- a/src/python/espressomd/electrostatics.pyx
+++ b/src/python/espressomd/electrostatics.pyx
@@ -40,8 +40,8 @@ IF ELECTROSTATICS == 1:
                 "Subclasses of ElectrostaticInteraction must define the _set_params_in_es_core() method.")
     
         def _deactivate_method(self):
-            coulomb.method = COULOMB_NONE
-            mpi_bcast_coulomb_params()
+            deactivate_coulomb_method()
+            handle_errors("Coulom method deactivation")
     
         def Tune(self, **subsetTuneParams):
     
@@ -69,8 +69,8 @@ IF COULOMB_DEBYE_HUECKEL:
 
         Parameters
         ----------
-        bjerrum_length      : float
-                              Bjerrum length
+        prefactor      : float
+                             Electrostatics prefactor (see :eq:`coulomb_prefactor`)
         kappa               : float
                               Inverse Debye screening length
         r_cut               : float 
@@ -88,8 +88,8 @@ IF COULOMB_DEBYE_HUECKEL:
 
         """
         def validate_params(self):
-            if (self._params["bjerrum_length"] <= 0):
-                raise ValueError("Bjerrum_length should be a positive double")
+            if (self._params["prefactor"] <= 0):
+                raise ValueError("Prefactor should be a positive float")
             if (self._params["kappa"] < 0):
                 raise ValueError("kappa should be a non-negative double")
             if (self._params["r_cut"] < 0):
@@ -108,13 +108,13 @@ IF COULOMB_DEBYE_HUECKEL:
                 raise ValueError("alpha should be a non-negative double")
 
         def valid_keys(self):
-            return "bjerrum_length", "kappa", "r_cut", "eps_int", "eps_ext", "r0", "r1", "alpha"
+            return "prefactor", "kappa", "r_cut", "eps_int", "eps_ext", "r0", "r1", "alpha"
 
         def required_keys(self):
-            return "bjerrum_length", "kappa", "r_cut", "eps_int", "eps_ext", "r0", "r1", "alpha"
+            return "prefactor", "kappa", "r_cut", "eps_int", "eps_ext", "r0", "r1", "alpha"
 
         def _set_params_in_es_core(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             dh_set_params_cdh(self._params["kappa"], self._params["r_cut"], self._params[
                               "eps_int"], self._params["eps_ext"], self._params["r0"], self._params["r1"], self._params["alpha"])
 
@@ -128,7 +128,7 @@ IF COULOMB_DEBYE_HUECKEL:
             self._set_params_in_es_core()
 
         def default_params(self):
-            return {"bjerrum_length": -1,
+            return {"prefactor": -1,
                     "kappa": -1,
                     "r_cut": -1,
                     "eps_int": -1,
@@ -146,8 +146,8 @@ ELSE:
 
             Parameters
             ----------
-            bjerrum_length  : float
-                              Bjerrum length
+            prefactor  : float
+                             Electrostatics prefactor (see :eq:`coulomb_prefactor`)
             kappa           : float
                               Inverse Debye sreening length
             r_cut           : float
@@ -155,22 +155,22 @@ ELSE:
 
             """
             def validate_params(self):
-                if (self._params["bjerrum_length"] <= 0):
+                if (self._params["prefactor"] <= 0):
                     raise ValueError(
-                        "Bjerrum_length should be a positive double")
+                        "prefactor should be a positive float")
                 if (self._params["kappa"] < 0):
                     raise ValueError("kappa should be a non-negative double")
                 if (self._params["r_cut"] < 0):
                     raise ValueError("r_cut should be a non-negative double")
 
             def valid_keys(self):
-                return "bjerrum_length", "kappa", "r_cut"
+                return "prefactor", "kappa", "r_cut"
 
             def required_keys(self):
-                return "bjerrum_length", "kappa", "r_cut"
+                return "prefactor", "kappa", "r_cut"
 
             def _set_params_in_es_core(self):
-                coulomb_set_bjerrum(self._params["bjerrum_length"])
+                coulomb_set_prefactor(self._params["prefactor"])
                 dh_set_params(self._params["kappa"], self._params["r_cut"])
 
             def _get_params_from_es_core(self):
@@ -183,7 +183,7 @@ ELSE:
                 self._set_params_in_es_core()
 
             def default_params(self):
-                return {"bjerrum_length": -1,
+                return {"prefactor": -1,
                         "kappa": -1,
                         "r_cut": -1}
 
@@ -200,10 +200,8 @@ IF P3M == 1:
 
             Parameters
             ----------
-            bjerrum_length : float
-                             The separation at which the electrostatic interaction
-                             energy and thermal energy are comparable. Defined as
-                             :math:`l_B = e_o^2 / (4 \pi \epsilon k_B T)` 
+            prefactor : float
+                    Electrostatics prefactor (see :eq:`coulomb_prefactor`)
             accuracy : float
                        P3M tunes its parameters to provide this target accuracy.
 
@@ -240,8 +238,8 @@ IF P3M == 1:
 
         def validate_params(self):
             default_params = self.default_params()
-            if not (self._params["bjerrum_length"] > 0.0):
-                raise ValueError("Bjerrum_length should be a positive double")
+            if not (self._params["prefactor"] > 0.0):
+                raise ValueError("prefactor should be a positive float")
 
             if not (self._params["r_cut"] >= 0 or self._params["r_cut"] == default_params["r_cut"]):
                 raise ValueError("P3M r_cut has to be >=0")
@@ -280,10 +278,10 @@ IF P3M == 1:
                     "alpha should be positive")
 
         def valid_keys(self):
-            return "mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut", "bjerrum_length", "tune"
+            return "mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut", "prefactor", "tune"
 
         def required_keys(self):
-            return ["bjerrum_length", "accuracy"]
+            return ["prefactor", "accuracy"]
 
         def default_params(self):
             return {"cao": 0,
@@ -299,13 +297,13 @@ IF P3M == 1:
         def _get_params_from_es_core(self):
             params = {}
             params.update(p3m.params)
-            params["bjerrum_length"] = coulomb.bjerrum
+            params["prefactor"] = coulomb.prefactor
             params["tune"] = self._params["tune"]
             return params
 
         def _set_params_in_es_core(self):
             #Sets lb, bcast, resets vars to zero if lb=0
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             #Sets cdef vars and calls p3m_set_params() in core 
             python_p3m_set_params(self._params["r_cut"],
                         self._params["mesh"], self._params["cao"],
@@ -320,7 +318,7 @@ IF P3M == 1:
             python_p3m_set_mesh_offset(self._params["mesh_off"])
 
         def _tune(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             python_p3m_set_tune_params(self._params["r_cut"], self._params["mesh"], self._params[
                                        "cao"], -1.0, self._params["accuracy"], self._params["inter"])
             resp = python_p3m_adaptive_tune()
@@ -346,10 +344,8 @@ IF P3M == 1:
 
                 Parameters
                 ----------
-                bjerrum_length : float
-                                 The separation at which the electrostatic interaction
-                                 energy and thermal energy are comparable. Defined as
-                                 :math:`l_B = e_o^2 / (4 \pi \epsilon k_B T)` 
+                prefactor : float
+                           Electrostatics prefactor (see :eq:`coulomb_prefactor`)
                 accuracy : float
                            P3M tunes its parameters to provide this target accuracy.
 
@@ -421,10 +417,10 @@ IF P3M == 1:
                         "mesh_off should be a list of length 3 and values between 0.0 and 1.0")
 
             def valid_keys(self):
-                return "mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut", "bjerrum_length", "tune"
+                return "mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut", "prefactor", "tune"
 
             def required_keys(self):
-                return ["bjerrum_length", "accuracy"]
+                return ["prefactor", "accuracy"]
 
             def default_params(self):
                 return {"cao": 0,
@@ -440,12 +436,12 @@ IF P3M == 1:
             def _get_params_from_es_core(self):
                 params = {}
                 params.update(p3m.params)
-                params["bjerrum_length"] = coulomb.bjerrum
+                params["prefactor"] = coulomb.prefactor
                 params["tune"] = self._params["tune"]
                 return params
 
             def _tune(self):
-                coulomb_set_bjerrum(self._params["bjerrum_length"])
+                coulomb_set_prefactor(self._params["prefactor"])
                 python_p3m_set_tune_params(self._params["r_cut"], self._params["mesh"], self._params[
                                            "cao"], -1.0, self._params["accuracy"], self._params["inter"])
                 resp = python_p3m_adaptive_tune()
@@ -463,12 +459,13 @@ IF P3M == 1:
                 self._set_params_in_es_core()
 
             def _set_params_in_es_core(self):
-                coulomb_set_bjerrum(self._params["bjerrum_length"])
+                coulomb_set_prefactor(self._params["prefactor"])
                 python_p3m_set_params(self._params["r_cut"], self._params["mesh"], self._params[
                                       "cao"], self._params["alpha"], self._params["accuracy"])
                 p3m_set_eps(self._params["epsilon"])
                 p3m_set_ninterpol(self._params["inter"])
                 python_p3m_set_mesh_offset(self._params["mesh_off"])
+                handle_errors("p3m gpu init" )
 
 IF ELECTROSTATICS and CUDA and EWALD_GPU:
     cdef class EwaldGpu(ElectrostaticInteraction):
@@ -482,10 +479,8 @@ IF ELECTROSTATICS and CUDA and EWALD_GPU:
 
             Parameters
             ----------
-            bjerrum_length : float
-                             The separation at which the electrostatic interaction
-                             energy and thermal energy are comparable. Defined as
-                             :math:`l_B = e_o^2 / (4 \pi \epsilon k_B T)` 
+            prefactor : float
+                    Electrostatics prefactor (see :eq:`coulomb_prefactor`)
             accuracy : float
                        Maximal allowed root mean square error regarding the forces
 
@@ -524,10 +519,10 @@ IF ELECTROSTATICS and CUDA and EWALD_GPU:
             del self.thisptr
 
         def valid_keys(self):
-            return "bjerrum_length", "rcut", "num_kx", "num_ky", "num_kz",  "K_max", "alpha", "accuracy", "precision", "time_calc_steps"
+            return "prefactor", "rcut", "num_kx", "num_ky", "num_kz",  "K_max", "alpha", "accuracy", "precision", "time_calc_steps"
 
         def default_params(self):
-            return {"bjerrum_length": -1,
+            return {"prefactor": -1,
                     "rcut": -1,
                     "num_kx": -1,
                     "num_ky": -1,
@@ -542,8 +537,8 @@ IF ELECTROSTATICS and CUDA and EWALD_GPU:
 
         def validate_params(self):
             default_params = self.default_params()
-            if self._params["bjerrum_length"] <= 0.0 and self._params["bjerrum_length"] != default_params["bjerrum_length"]:
-                raise ValueError("Bjerrum_length should be a positive double")
+            if self._params["prefactor"] <= 0.0:
+                raise ValueError("prefactor should be a positive float")
             if isinstance(self._params["K_max"], (list, np.ndarray)):
                 if is_valid_type(self._params["K_max"], int) and len(self._params["K_max"]) == 3:
                     self._params["num_kx"] = self._params["K_max"][0]
@@ -569,10 +564,10 @@ IF ELECTROSTATICS and CUDA and EWALD_GPU:
                 raise ValueError("alpha has to be a positive double")
 
         def required_keys(self):
-            return "bjerrum_length", "accuracy", "precision", "K_max"
+            return "prefactor", "accuracy", "precision", "K_max"
 
         def _tune(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             default_params = self.default_params()
             if self._params["time_calc_steps"] == default_params["time_calc_steps"]:
                 self._params[
@@ -588,14 +583,14 @@ IF ELECTROSTATICS and CUDA and EWALD_GPU:
                 print(self.log)
 
         def _set_params_in_es_core(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             self.thisptr.set_params(self._params["rcut"], self._params[
                                     "num_kx"], self._params["num_ky"], self._params["num_kz"], self._params["alpha"])
 
         def _get_params_from_es_core(self):
             params = {}
             params.update(ewaldgpu_params)
-            params["bjerrum_length"] = coulomb.bjerrum
+            params["prefactor"] = coulomb.prefactor
             return params
 
         def _activate_method(self):
@@ -615,8 +610,8 @@ IF ELECTROSTATICS:
 
         Parameters
         ----------
-        bjerrum_length      : float
-                              Bjerrum length
+        prefactor      : float
+                             Electrostatics prefactor (see :eq:`coulomb_prefactor`)
         maxWPerror          : float
                               Maximal pairwise error
         far_switch_radius   : float, optional
@@ -628,8 +623,8 @@ IF ELECTROSTATICS:
 
         def validate_params(self):
             default_params = self.default_params()
-            if self._params["bjerrum_length"] < 0:
-                raise ValueError("Bjerrum_length should be a positive double")
+            if self._params["prefactor"] <= 0.:
+                raise ValueError("Prefactor should be a positive float")
             if self._params["maxPWerror"] < 0 and self._params["maxPWerror"] != default_params["maxPWerror"]:
                 raise ValueError("maxPWerror should be a positive double")
             if self._params["far_switch_radius"] < 0 and self._params["far_switch_radius"] != default_params["far_switch_radius"]:
@@ -638,29 +633,29 @@ IF ELECTROSTATICS:
                 raise ValueError("bessel_cutoff should be a positive integer")
 
         def default_params(self):
-            return {"bjerrum_length": -1,
+            return {"prefactor": -1,
                     "maxPWerror": -1,
                     "far_switch_radius": -1,
                     "bessel_cutoff": -1,
                     "tune": True}
 
         def valid_keys(self):
-            return "bjerrum_length", "maxPWerror", "far_switch_radius", "bessel_cutoff", "tune"
+            return "prefactor", "maxPWerror", "far_switch_radius", "bessel_cutoff", "tune"
 
         def required_keys(self):
-            return ["bjerrum_length", "maxPWerror"]
+            return ["prefactor", "maxPWerror"]
 
         def _get_params_from_es_core(self):
             params = {}
             params.update(mmm1d_params)
             params["far_switch_radius"] = np.sqrt(params["far_switch_radius_2"])
             del params["far_switch_radius_2"]
-            params["bjerrum_length"] = coulomb.bjerrum
+            params["prefactor"] = coulomb.prefactor
             params["tune"] = self._params["tune"]
             return params
 
         def _set_params_in_es_core(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             MMM1D_set_params(
                 self._params["far_switch_radius"], self._params["maxPWerror"])
 
@@ -687,8 +682,8 @@ IF ELECTROSTATICS and MMM1D_GPU:
 
         Parameters
         ----------
-        bjerrum_length      : float
-                              Bjerrum length
+        prefactor      : float
+                             Electrostatics prefactor (see :eq:`coulomb_prefactor`)
         maxWPerror          : float
                               Maximal pairwise error
         far_switch_radius   : float, optional
@@ -714,8 +709,8 @@ IF ELECTROSTATICS and MMM1D_GPU:
 
         def validate_params(self):
             default_params = self.default_params()
-            if self._params["bjerrum_length"] < 0:
-                raise ValueError("Bjerrum_length should be a positive double")
+            if self._params["prefactor"] <= 0:
+                raise ValueError("prefactor should be a positive float")
             if self._params["maxPWerror"] < 0 and self._params["maxPWerror"] != default_params["maxPWerror"]:
                 raise ValueError("maxPWerror should be a positive double")
             if self._params["far_switch_radius"] < 0 and self._params["far_switch_radius"] != default_params["far_switch_radius"]:
@@ -724,29 +719,29 @@ IF ELECTROSTATICS and MMM1D_GPU:
                 raise ValueError("bessel_cutoff should be a positive integer")
 
         def default_params(self):
-            return {"bjerrum_length": -1,
+            return {"prefactor": -1,
                     "maxPWerror": -1.0,
                     "far_switch_radius": -1.0,
                     "bessel_cutoff": -1,
                     "tune": True}
 
         def valid_keys(self):
-            return "bjerrum_length", "maxPWerror", "far_switch_radius", "bessel_cutoff", "tune"
+            return "prefactor", "maxPWerror", "far_switch_radius", "bessel_cutoff", "tune"
 
         def required_keys(self):
-            return ["bjerrum_length", "maxPWerror"]
+            return ["prefactor", "maxPWerror"]
 
         def _get_params_from_es_core(self):
             params = {}
             params.update(mmm1d_params)
-            params["bjerrum_length"] = coulomb.bjerrum
+            params["prefactor"] = coulomb.prefactor
             return params
 
         def _set_params_in_es_core(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             default_params = self.default_params()
 
-            self.thisptr.set_params(globals.box_l[2], globals.temperature * coulomb.bjerrum, self._params[
+            self.thisptr.set_params(globals.box_l[2], coulomb.prefactor, self._params[
                                     "maxPWerror"], self._params["far_switch_radius"], self._params["bessel_cutoff"])
 
         def _tune(self):
@@ -770,8 +765,8 @@ IF ELECTROSTATICS:
 
         Parameters
         ----------
-        bjerrum_length          : float
-                                  Bjerrum length
+        prefactor          : float
+                                  Electrostatics prefactor (see :eq:`coulomb_prefactor`)
         maxWPerror              : float
                                   Maximal pairwise error
         dielectric              : int, optional
@@ -818,8 +813,8 @@ IF ELECTROSTATICS:
         """
         def validate_params(self):
             default_params = self.default_params()
-            if self._params["bjerrum_length"] < 0:
-                raise ValueError("Bjerrum_length should be a positive double")
+            if self._params["prefactor"] <= 0:
+                raise ValueError("prefactor should be a positive float")
             if self._params["maxPWerror"] < 0 and self._params["maxPWerror"] != default_params["maxPWerror"]:
                 raise ValueError("maxPWerror should be a positive double")
             if self._params["dielectric"] == 1 and (self._params["top"] < 0 or self._params["mid"] < 0 or self._params["bot"] < 0):
@@ -831,7 +826,7 @@ IF ELECTROSTATICS:
                     "dielectric and dielectric_contrast are mutually exclusive!")
 
         def default_params(self):
-            return {"bjerrum_length": -1,
+            return {"prefactor": -1,
                     "maxPWerror": -1,
                     "far_cut": -1,
                     "top": 0,
@@ -845,15 +840,15 @@ IF ELECTROSTATICS:
                     "pot_diff": 0}
 
         def required_keys(self):
-            return ["bjerrum_length", "maxPWerror"]
+            return ["prefactor", "maxPWerror"]
 
         def valid_keys(self):
-            return "bjerrum_length", "maxPWerror", "top", "mid", "bot", "delta_mid_top", "delta_mid_bot", "pot_diff", "dielectric", "dielectric_contrast_on", "const_pot", "far_cut"
+            return "prefactor", "maxPWerror", "top", "mid", "bot", "delta_mid_top", "delta_mid_bot", "pot_diff", "dielectric", "dielectric_contrast_on", "const_pot", "far_cut"
 
         def _get_params_from_es_core(self):
             params = {}
             params.update(mmm2d_params)
-            params["bjerrum_length"] = coulomb.bjerrum
+            params["prefactor"] = coulomb.prefactor
             if params["dielectric_contrast_on"] == 1 or params["const_pot"] == 1:
                 params["dielectric"] = 1
             else:
@@ -861,7 +856,7 @@ IF ELECTROSTATICS:
             return params
 
         def _set_params_in_es_core(self):
-            coulomb_set_bjerrum(self._params["bjerrum_length"])
+            coulomb_set_prefactor(self._params["prefactor"])
             if self._params["dielectric"]:
                 self._params["delta_mid_top"] = (self._params[
                                                  "mid"] - self._params["top"]) / (self._params["mid"] + self._params["top"])
@@ -902,7 +897,7 @@ IF ELECTROSTATICS:
                 actors.Actor.__init__(self, *args, **kwargs)
 
             def _activate_method(self):
-                coulomb_set_bjerrum(self._params["bjerrum_length"])
+                coulomb_set_prefactor(self._params["prefactor"])
                 self._set_params_in_es_core()
                 mpi_bcast_coulomb_params()
 
diff --git a/src/python/espressomd/magnetostatics.pxd b/src/python/espressomd/magnetostatics.pxd
index 975cc5583..81c3d6dca 100644
--- a/src/python/espressomd/magnetostatics.pxd
+++ b/src/python/espressomd/magnetostatics.pxd
@@ -16,11 +16,10 @@ IF DIPOLES == 1:
             DIPOLAR_MDLC_DS,
             DIPOLAR_SCAFACOS
 
-        int dipolar_set_Dbjerrum(double bjerrum)
+        int dipolar_set_Dprefactor(double prefactor)
 
         ctypedef struct coulomb_parameters "Coulomb_parameters":
             double Dprefactor
-            double Dbjerrum
             dipolar_interaction Dmethod
         cdef extern coulomb_parameters coulomb
 
diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx
index a50792872..3d92ecaa4 100644
--- a/src/python/espressomd/magnetostatics.pyx
+++ b/src/python/espressomd/magnetostatics.pyx
@@ -34,15 +34,8 @@ IF DIPOLES == 1:
 
         Attributes
         ----------
-        bjerrum_length : :obj:`float`
-            Similar to the Bjerrum length in electrostatics.
-            Gives the separation at which the dipolar interaction energy
-            corrensponds to a thermal energy of :math:`1 k_B T`.
-            If not set, an explicit choice for the `prefactor` has to be made.
-
-        prefactor : :obj:`float`, optional
-            If given, the magnetostatic prefactor is explicitly set to
-            this value. Non-optional if `bjerrum_length` is not set.
+        prefactor : :obj:`float`
+                Magnetostatics prefactor
 
         """
 
@@ -50,44 +43,17 @@ IF DIPOLES == 1:
             """Check validity of given parameters.
 
             """
-            if not (("bjerrum_length" in self._params) ^ ("prefactor" in self._params)):
-                raise ValueError(
-                    "Either the bjerrum length or the explicit prefactor has to be given")
-
-            if "bjerrum_length" in self._params:
-                if not (self._params["bjerrum_length"] > 0.0):
-                    raise ValueError(
-                        "Bjerrum_length should be a positive double")
-            if "prefactor" in self._params:
-                if not self._params["prefactor"] > 0:
-                    raise ValueError("prefactor should be a positive double")
+            if not self._params["prefactor"] >= 0:
+                raise ValueError("prefactor should be a positive float")
 
         def set_magnetostatics_prefactor(self):
             """
-            Set the magnetostatics prefactor using either the `bjerrum_length`
-            or, if given, `prefactor`.
-
+            Set the magnetostatics prefactor 
+            
             """
-            if "bjerrum_length" in self._params:
-                if temperature == 0:
+            if dipolar_set_Dprefactor(self._params["prefactor"]):
                     raise Exception(
-                        "Bjerrum length is not defined, if temperature is zero")
-                if dipolar_set_Dbjerrum(self._params["bjerrum_length"]):
-                    raise Exception(
-                        "Could not set magnetostatic bjerrum length")
-                return True
-            if "prefactor" in self._params:
-                if temperature == 0.:
-                    if dipolar_set_Dbjerrum(self._params["prefactor"]):
-                        raise Exception(
-                            "Could not set magnetostatic prefactor")
-                else:
-                    if dipolar_set_Dbjerrum(self._params["prefactor"] / temperature):
-                        raise Exception(
-                            "Could not set magnetostatic prefactor")
-                    else:
-                        self._params["bjerrum_length"] = self._params[
-                            "prefactor"] / temperature
+                        "Could not set magnetostatic prefactor")
             # also necessary on 1 CPU or GPU, does more than just broadcasting
             mpi_bcast_coulomb_params()
 
@@ -99,7 +65,7 @@ IF DIPOLES == 1:
             return coulomb.Dmethod
 
         def _deactivate_method(self):
-            dipolar_set_Dbjerrum(0.0)
+            dipolar_set_Dprefactor(0.0)
             coulomb.Dmethod = DIPOLAR_NONE
             mpi_bcast_coulomb_params()
 
@@ -167,7 +133,7 @@ IF DP3M == 1:
                     "mesh_off should be a list of length 3 and values between 0.0 and 1.0")
 
         def valid_keys(self):
-            return "prefactor", "alpha_L", "r_cut_iL", "mesh", "mesh_off", "cao", "inter", "accuracy", "epsilon", "cao_cut", "a", "ai", "alpha", "r_cut", "inter2", "cao3", "additional_mesh", "bjerrum_length", "tune"
+            return "prefactor", "alpha_L", "r_cut_iL", "mesh", "mesh_off", "cao", "inter", "accuracy", "epsilon", "cao_cut", "a", "ai", "alpha", "r_cut", "inter2", "cao3", "additional_mesh", "tune"
 
         def required_keys(self):
             return ["accuracy", ]
@@ -278,7 +244,7 @@ IF DIPOLES == 1:
             return ()
 
         def valid_keys(self):
-            return ("bjerrum_length", "prefactor")
+            return ("prefactor",) 
 
         def _get_params_from_es_core(self):
             return {"prefactor": coulomb.Dprefactor}
@@ -313,7 +279,7 @@ IF DIPOLES == 1:
             return ("n_replica",)
 
         def valid_keys(self):
-            return ("bjerrum_length", "prefactor", "n_replica")
+            return ("prefactor", "n_replica")
 
         def _get_params_from_es_core(self):
             return {"prefactor": coulomb.Dprefactor, "n_replica": Ncut_off_magnetic_dipolar_direct_sum}
@@ -341,9 +307,8 @@ IF DIPOLES == 1:
                 Actor.__init__(self, *args, **kwargs)
 
             def _activate_method(self):
-                dipolar_set_Dbjerrum(self._params["bjerrum_length"])
+                dipolar_set_Dprefactor(self._params["prefactor"])
                 self._set_params_in_es_core()
-                mpi_bcast_coulomb_params()
             
             def _deactivate_method(self):
                 coulomb.Dmethod = DIPOLAR_NONE
@@ -371,7 +336,7 @@ IF DIPOLES == 1:
                 return ()
     
             def valid_keys(self):
-                return ("bjerrum_length", "prefactor")
+                return ("prefactor",)
     
             def _get_params_from_es_core(self):
                 return {"prefactor": coulomb.Dprefactor}
diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx
index 28a1885c9..6f3dab6f1 100644
--- a/src/python/espressomd/scafacos.pyx
+++ b/src/python/espressomd/scafacos.pyx
@@ -30,10 +30,10 @@ IF SCAFACOS == 1:
             self._set_params_in_es_core()
 
         def valid_keys(self):
-            return "method_name", "method_params", "bjerrum_length"
+            return "method_name", "method_params", "prefactor"
 
         def required_keys(self):
-            return "method_name", "method_params", "bjerrum_length"
+            return "method_name", "method_params", "prefactor"
 
         def validate_params(self):
             return True
@@ -79,11 +79,11 @@ IF SCAFACOS == 1:
             
             res["method_params"]=method_params
             
-            # Re-add the Bjerrum length to the parameter set
+            # Re-add the prefactor to the parameter set
             if self.dipolar: 
-                res["bjerrum_length"] =magnetostatics.coulomb.Dbjerrum
+                res["prefactor"] =magnetostatics.coulomb.Dprefactor
             else:
-                res["bjerrum_length"] =electrostatics.coulomb.bjerrum
+                res["prefactor"] =electrostatics.coulomb.prefactor
             return res
 
         def _set_params_in_es_core(self):
diff --git a/testsuite/analyze_energy.py b/testsuite/analyze_energy.py
index 890d704c1..fbf7278b0 100644
--- a/testsuite/analyze_energy.py
+++ b/testsuite/analyze_energy.py
@@ -141,7 +141,7 @@ class AnalyzeEnergy(ut.TestCase):
         self.system.part[1].pos = [3, 2, 2]
         self.system.part[0].q = 1
         self.system.part[1].q = -1
-        p3m = electrostatics.P3M(bjerrum_length=1.0,
+        p3m = electrostatics.P3M(prefactor=1.0,
                                  accuracy=9.910945054074526e-08,
                                  mesh=[22,22,22],
                                  cao=7,
diff --git a/testsuite/coulomb_cloud_wall.py b/testsuite/coulomb_cloud_wall.py
index e843ff616..f768b4a95 100644
--- a/testsuite/coulomb_cloud_wall.py
+++ b/testsuite/coulomb_cloud_wall.py
@@ -62,20 +62,20 @@ class CoulombCloudWall(ut.TestCase):
             self.S.part.add(id=int(id), pos=pos, q=q)
             self.forces[id] = f
 
-    def compare(self, method_name, energy=True):
+    def compare(self, method_name, energy=True,prefactor=None):
         # Compare forces and energy now in the system to stored ones
 
         # Force
         force_abs_diff = 0.
         for p in self.S.part:
-            force_abs_diff += abs(np.sqrt(sum((p.f - self.forces[p.id])**2)))
+            force_abs_diff += abs(np.sqrt(sum((p.f/prefactor - self.forces[p.id])**2)))
         force_abs_diff /= len(self.S.part)
 
         print(method_name, "force difference", force_abs_diff)
 
         # Energy
         if energy:
-            energy_abs_diff = abs(self.S.analysis.energy()["total"] - self.reference_energy)
+            energy_abs_diff = abs(self.S.analysis.energy()["total"]/prefactor - self.reference_energy)
             print(method_name, "energy difference", energy_abs_diff)
             self.assertLessEqual(
                 energy_abs_diff,
@@ -96,16 +96,16 @@ class CoulombCloudWall(ut.TestCase):
 
     if espressomd.has_features(["P3M"]):
         def test_p3m(self):
-            self.S.actors.add(P3M(bjerrum_length=1, r_cut=1.001, accuracy=1e-3,
+            self.S.actors.add(P3M(prefactor=3, r_cut=1.001, accuracy=1e-3,
                                   mesh=64, cao=7, alpha=2.70746, tune=False))
             self.S.integrator.run(0)
-            self.compare("p3m", energy=True)
+            self.compare("p3m", energy=True,prefactor=3)
 
     if espressomd.has_features(["ELECTROSTATICS", "CUDA"]):
         def test_p3m_gpu(self):
             self.S.actors.add(
                 P3MGPU(
-                    bjerrum_length=1,
+                    prefactor=2.2,
                     r_cut=1.001,
                     accuracy=1e-3,
                     mesh=64,
@@ -113,14 +113,14 @@ class CoulombCloudWall(ut.TestCase):
                     alpha=2.70746,
                     tune=False))
             self.S.integrator.run(0)
-            self.compare("p3m_gpu", energy=False)
+            self.compare("p3m_gpu", energy=False,prefactor=2.2)
 
     if espressomd.has_features(["SCAFACOS"]):
         if "p3m" in scafacos.available_methods():
             def test_scafacos_p3m(self):
                 self.S.actors.add(
                     Scafacos(
-                        bjerrum_length=1,
+                        prefactor=0.5,
                         method_name="p3m",
                         method_params={
                             "p3m_r_cut": 1.001,
@@ -128,23 +128,27 @@ class CoulombCloudWall(ut.TestCase):
                             "p3m_cao": 7,
                             "p3m_alpha": 2.70746}))
                 self.S.integrator.run(0)
-                self.compare("scafacos_p3m", energy=True)
+                self.compare("scafacos_p3m", energy=True,prefactor=0.5)
 
         if "p2nfft" in scafacos.available_methods():
             def test_scafacos_p2nfft(self):
                 self.S.actors.add(
                     Scafacos(
-                        bjerrum_length=1,
+                        prefactor=2.8,
                         method_name="p2nfft",
                         method_params={
                             "p2nfft_r_cut": 1.001,
                             "tolerance_field": 1E-4}))
                 self.S.integrator.run(0)
-                self.compare("scafacos_p2nfft", energy=True)
+                self.compare("scafacos_p2nfft", energy=True,prefactor=2.8)
 
     def test_zz_deactivation(self):
-        # Is the energy 0, if no methods active
+        # Is the energy and force 0, if no methods active
         self.assertEqual(self.S.analysis.energy()["total"], 0.0)
+        self.S.integrator.run(0,recalc_forces=True)
+        for p in self.S.part:
+            self.assertAlmostEqual(np.linalg.norm(p.f),0,places=11)
+
 
 
 if __name__ == "__main__":
diff --git a/testsuite/coulomb_cloud_wall_duplicated.py b/testsuite/coulomb_cloud_wall_duplicated.py
index bd4dc8135..91017a250 100644
--- a/testsuite/coulomb_cloud_wall_duplicated.py
+++ b/testsuite/coulomb_cloud_wall_duplicated.py
@@ -82,7 +82,7 @@ class CoulombCloudWall(ut.TestCase):
 
         if "P3M" in espressomd.features():
             def test_p3m(self):
-                self.S.actors.add(P3M(bjerrum_length=1, r_cut=1.001, accuracy = 1e-3,
+                self.S.actors.add(P3M(prefactor=1, r_cut=1.001, accuracy = 1e-3,
                                       mesh=[64, 64, 128], cao=7, alpha=2.70746, tune=False))
                 self.S.integrator.run(0)
                 self.compare("p3m", energy=True)
@@ -91,7 +91,7 @@ class CoulombCloudWall(ut.TestCase):
             def test_p3m_gpu(self):
                 self.S.actors.add(
                     P3MGPU(
-                        bjerrum_length=1,
+                        prefactor=1,
                         r_cut=1.001,
                         accuracy=1e-3,
                         mesh=[64, 64, 128],
diff --git a/testsuite/coulomb_mixed_periodicity.py b/testsuite/coulomb_mixed_periodicity.py
index e8d70b9c8..b3aefbc42 100644
--- a/testsuite/coulomb_mixed_periodicity.py
+++ b/testsuite/coulomb_mixed_periodicity.py
@@ -113,7 +113,7 @@ class CoulombMixedPeriodicity(ut.TestCase):
             self.S.periodicity=1,1,1
             self.S.box_l = (10, 10, 10)
 
-            p3m=el.P3M(bjerrum_length=1, accuracy = 1e-6,mesh=(64,64,64))
+            p3m=el.P3M(prefactor=1, accuracy = 1e-6,mesh=(64,64,64))
             
             self.S.actors.add(p3m)
             elc=el_ext.ELC(maxPWerror=1E-6,gap_size=1)
@@ -126,7 +126,7 @@ class CoulombMixedPeriodicity(ut.TestCase):
         self.S.box_l = (10, 10, 10)
         self.S.cell_system.set_layered(n_layers=10,use_verlet_lists=False)
         self.S.periodicity=1,1,0
-        mmm2d=(el.MMM2D(bjerrum_length=1,maxPWerror=1E-7))
+        mmm2d=(el.MMM2D(prefactor=1,maxPWerror=1E-7))
             
         self.S.actors.add(mmm2d)
         self.S.integrator.run(0)
@@ -146,7 +146,7 @@ class CoulombMixedPeriodicity(ut.TestCase):
                 self.S.box_l=10,10,10
                 
                 scafacos=el.Scafacos(
-                        bjerrum_length=1,
+                        prefactor=1,
                         method_name="p2nfft",
                         method_params={
                             "tolerance_field": 1E-5})
diff --git a/testsuite/coulomb_tuning.py b/testsuite/coulomb_tuning.py
index 8da6b0bdc..b410643b4 100644
--- a/testsuite/coulomb_tuning.py
+++ b/testsuite/coulomb_tuning.py
@@ -77,7 +77,7 @@ class CoulombCloudWallTune(ut.TestCase):
         def test_p3m(self):
             # We have to add some tolerance here, because the reference
             # system is not homogeneous
-            self.S.actors.add(P3M(bjerrum_length=1., accuracy=5e-4,
+            self.S.actors.add(P3M(prefactor=1., accuracy=5e-4,
                                   tune=True))
             self.S.integrator.run(0)
             self.compare("p3m")
@@ -86,7 +86,7 @@ class CoulombCloudWallTune(ut.TestCase):
         def test_p3m_gpu(self):
             # We have to add some tolerance here, because the reference
             # system is not homogeneous
-            self.S.actors.add(P3MGPU(bjerrum_length=1., accuracy=5e-4,
+            self.S.actors.add(P3MGPU(prefactor=1., accuracy=5e-4,
                                       tune=True))
             self.S.integrator.run(0)
             self.compare("p3m_gpu")
diff --git a/testsuite/dawaanr-and-dds-gpu.py b/testsuite/dawaanr-and-dds-gpu.py
index a0bc8f7c5..9e6f04065 100644
--- a/testsuite/dawaanr-and-dds-gpu.py
+++ b/testsuite/dawaanr-and-dds-gpu.py
@@ -87,7 +87,7 @@ class DDSGPUTest(ut.TestCase):
             # and torque
             self.es.thermostat.set_langevin(kT=1.297, gamma=0.0)
 
-            dds_cpu = DipolarDirectSumCpu(bjerrum_length=pf_dawaanr)
+            dds_cpu = DipolarDirectSumCpu(prefactor=pf_dawaanr)
             self.es.actors.add(dds_cpu)
             self.es.integrator.run(steps=0, recalc_forces=True)
 
@@ -104,7 +104,7 @@ class DDSGPUTest(ut.TestCase):
                 self.es.actors.remove(self.es.actors.active_actors[i])
 
             self.es.integrator.run(steps=0, recalc_forces=True)
-            dds_gpu = DipolarDirectSumGpu(bjerrum_length=pf_dds_gpu)
+            dds_gpu = DipolarDirectSumGpu(prefactor=pf_dds_gpu)
             self.es.actors.add(dds_gpu)
             self.es.integrator.run(steps=0, recalc_forces=True)
 
diff --git a/testsuite/dipolar_mdlc_p3m_scafacos_p2nfft.py b/testsuite/dipolar_mdlc_p3m_scafacos_p2nfft.py
index c2eb8a4b7..48870f9b5 100644
--- a/testsuite/dipolar_mdlc_p3m_scafacos_p2nfft.py
+++ b/testsuite/dipolar_mdlc_p3m_scafacos_p2nfft.py
@@ -162,7 +162,7 @@ class Dipolar_p3m_mdlc_p2nfft(ut.TestCase):
             s.part.add(id=int(p[0]), pos=p[1:4], dip=p[4:7],rotation=(1,1,1))
 
         scafacos = magnetostatics.Scafacos(
-            bjerrum_length=1,
+            prefactor=1,
             method_name="p2nfft",
             method_params={
                 "p2nfft_verbose_tuning": 0,
diff --git a/testsuite/ek_eof_one_species_x.py b/testsuite/ek_eof_one_species_x.py
index 4f16457b6..4b245a918 100644
--- a/testsuite/ek_eof_one_species_x.py
+++ b/testsuite/ek_eof_one_species_x.py
@@ -152,7 +152,7 @@ class ek_eof_one_species_x(ut.TestCase):
             viscosity=viscosity_kinematic,
             friction=friction,
             T=temperature,
-            bjerrum_length=bjerrum_length,
+            prefactor=bjerrum_length*temperature,
             stencil="linkcentered")
 
         counterions = electrokinetics.Species(
diff --git a/testsuite/ek_eof_one_species_x_nonlinear.py b/testsuite/ek_eof_one_species_x_nonlinear.py
index a158b029f..dc9fefc07 100644
--- a/testsuite/ek_eof_one_species_x_nonlinear.py
+++ b/testsuite/ek_eof_one_species_x_nonlinear.py
@@ -154,7 +154,7 @@ class ek_eof_one_species_x(ut.TestCase):
             viscosity=viscosity_kinematic,
             friction=friction,
             T=temperature,
-            bjerrum_length=bjerrum_length,
+            prefactor=bjerrum_length*temperature,
             stencil="nonlinear")
 
         counterions = electrokinetics.Species(
diff --git a/testsuite/ek_eof_one_species_y.py b/testsuite/ek_eof_one_species_y.py
index dda46159f..8cc68d06a 100644
--- a/testsuite/ek_eof_one_species_y.py
+++ b/testsuite/ek_eof_one_species_y.py
@@ -152,7 +152,7 @@ class ek_eof_one_species_x(ut.TestCase):
             viscosity=viscosity_kinematic,
             friction=friction,
             T=temperature,
-            bjerrum_length=bjerrum_length,
+            prefactor=bjerrum_length*temperature,
             stencil="linkcentered")
 
         counterions = electrokinetics.Species(
diff --git a/testsuite/ek_eof_one_species_y_nonlinear.py b/testsuite/ek_eof_one_species_y_nonlinear.py
index 226cfa221..e3346cccf 100644
--- a/testsuite/ek_eof_one_species_y_nonlinear.py
+++ b/testsuite/ek_eof_one_species_y_nonlinear.py
@@ -153,7 +153,7 @@ class ek_eof_one_species_x(ut.TestCase):
             viscosity=viscosity_kinematic,
             friction=friction,
             T=temperature,
-            bjerrum_length=bjerrum_length,
+            prefactor=bjerrum_length*temperature,
             stencil="nonlinear")
 
         counterions = electrokinetics.Species(
diff --git a/testsuite/ek_eof_one_species_z.py b/testsuite/ek_eof_one_species_z.py
index 43f223d19..9d236d193 100644
--- a/testsuite/ek_eof_one_species_z.py
+++ b/testsuite/ek_eof_one_species_z.py
@@ -151,7 +151,7 @@ class ek_eof_one_species_x(ut.TestCase):
             viscosity=viscosity_kinematic,
             friction=friction,
             T=temperature,
-            bjerrum_length=bjerrum_length,
+            prefactor=bjerrum_length*temperature,
             stencil="linkcentered")
 
         counterions = electrokinetics.Species(
diff --git a/testsuite/ek_eof_one_species_z_nonlinear.py b/testsuite/ek_eof_one_species_z_nonlinear.py
index d9b3df4ef..f5c5cca51 100644
--- a/testsuite/ek_eof_one_species_z_nonlinear.py
+++ b/testsuite/ek_eof_one_species_z_nonlinear.py
@@ -154,7 +154,7 @@ class ek_eof_one_species_x(ut.TestCase):
             viscosity=viscosity_kinematic,
             friction=friction,
             T=temperature,
-            bjerrum_length=bjerrum_length,
+            prefactor=bjerrum_length*temperature,
             stencil="nonlinear")
 
         counterions = electrokinetics.Species(
diff --git a/testsuite/elc_vs_mmm2d_neutral.py b/testsuite/elc_vs_mmm2d_neutral.py
index 766955783..6745c3187 100644
--- a/testsuite/elc_vs_mmm2d_neutral.py
+++ b/testsuite/elc_vs_mmm2d_neutral.py
@@ -30,11 +30,11 @@ class ELC_vs_MMM2D_neutral(ut.TestCase):
         }
 
         mmm2d_param_sets = { 
-        "inert":         { "bjerrum_length": 1.0, "maxPWerror": self.acc },
-        "dielectric":    { "bjerrum_length": 1.0, "maxPWerror": self.acc, "dielectric_contrast_on": 1, "delta_mid_bot": 0.1, "delta_mid_top": 0.9 },
-        "const_pot_0":   { "bjerrum_length": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 0.0}, 
-        "const_pot_1":   { "bjerrum_length": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 1.0}, 
-        "const_pot_m1":  { "bjerrum_length": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": -1.0} 
+        "inert":         { "prefactor": 1.0, "maxPWerror": self.acc },
+        "dielectric":    { "prefactor": 1.0, "maxPWerror": self.acc, "dielectric_contrast_on": 1, "delta_mid_bot": 0.1, "delta_mid_top": 0.9 },
+        "const_pot_0":   { "prefactor": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 0.0}, 
+        "const_pot_1":   { "prefactor": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 1.0}, 
+        "const_pot_m1":  { "prefactor": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": -1.0} 
         }
 
         self.system.box_l = [self.box_l, self.box_l, self.box_l]
@@ -74,7 +74,7 @@ class ELC_vs_MMM2D_neutral(ut.TestCase):
         self.system.cell_system.set_domain_decomposition(use_verlet_lists = True)
         self.system.cell_system.node_grid = buf_node_grid
         self.system.periodicity = [1, 1, 1]
-        p3m = P3M(bjerrum_length=1.0, accuracy=self.acc, mesh = [16,16,24], cao = 6)
+        p3m = P3M(prefactor=1.0, accuracy=self.acc, mesh = [16,16,24], cao = 6)
         self.system.actors.add(p3m)
         
         elc = electrostatic_extensions.ELC(**elc_param_sets["inert"])
diff --git a/testsuite/elc_vs_mmm2d_nonneutral.py b/testsuite/elc_vs_mmm2d_nonneutral.py
index 71fb8c853..31cb891f3 100644
--- a/testsuite/elc_vs_mmm2d_nonneutral.py
+++ b/testsuite/elc_vs_mmm2d_nonneutral.py
@@ -29,10 +29,10 @@ class ELC_vs_MMM2D_neutral(ut.TestCase):
         }
 
         mmm2d_param_sets = { 
-        "inert":         { "bjerrum_length": 1.0, "maxPWerror": self.acc },
-        "const_pot_0":   { "bjerrum_length": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 0.0}, 
-        "const_pot_1":   { "bjerrum_length": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 1.0}, 
-        "const_pot_m1":  { "bjerrum_length": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": -1.0} 
+        "inert":         { "prefactor": 1.0, "maxPWerror": self.acc },
+        "const_pot_0":   { "prefactor": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 0.0}, 
+        "const_pot_1":   { "prefactor": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": 1.0}, 
+        "const_pot_m1":  { "prefactor": 1.0, "maxPWerror": self.acc, "const_pot": 1, "pot_diff": -1.0} 
         }
 
         self.system.box_l = [self.box_l, self.box_l, self.box_l]
@@ -68,7 +68,7 @@ class ELC_vs_MMM2D_neutral(ut.TestCase):
         self.system.cell_system.set_domain_decomposition(use_verlet_lists = True)
         self.system.cell_system.node_grid = buf_node_grid
         self.system.periodicity = [1, 1, 1]
-        p3m = P3M(bjerrum_length=1.0, accuracy=self.acc, mesh = [20,20,32], cao = 7)
+        p3m = P3M(prefactor=1.0, accuracy=self.acc, mesh = [20,20,32], cao = 7)
         self.system.actors.add(p3m)
         
         elc = electrostatic_extensions.ELC(**elc_param_sets["inert"])
diff --git a/testsuite/electrostaticInteractions.py b/testsuite/electrostaticInteractions.py
index c4df853ff..677af6b75 100644
--- a/testsuite/electrostaticInteractions.py
+++ b/testsuite/electrostaticInteractions.py
@@ -49,7 +49,7 @@ class ElectrostaticInteractionsTests(ut.TestCase):
         u = np.zeros_like(r)
         # r<r_cut
         i = np.where(r < df_params['r_cut'])[0]
-        u[i] = df_params['bjerrum_length'] * kT * q1 * \
+        u[i] = df_params['prefactor'] * kT * q1 * \
             q2 * np.exp(-df_params['kappa'] * r[i]) / r[i]
         return u
 
@@ -71,7 +71,7 @@ class ElectrostaticInteractionsTests(ut.TestCase):
                  r_cut=8.906249999999998,
                  alpha=0.387611049779351,
                  tune=False))
-        p3m = espressomd.electrostatics.P3M(bjerrum_length=1.0,
+        p3m = espressomd.electrostatics.P3M(prefactor=1.0,
                                             accuracy=9.910945054074526e-08,
                                             mesh=[22, 22, 22],
                                             cao=7,
@@ -83,16 +83,16 @@ class ElectrostaticInteractionsTests(ut.TestCase):
                                p3m_energy)
         # need to update forces
         self.system.integrator.run(0)
-        self.assertTrue(np.allclose(self.system.part[0].f,
-                                    [p3m_force, 0, 0]))
-        self.assertTrue(np.allclose(self.system.part[1].f,
-                                    [-p3m_force, 0, 0]))
+        np.testing.assert_allclose(np.copy(self.system.part[0].f),
+                                    [p3m_force, 0, 0],atol=1E-5)
+        np.testing.assert_allclose(np.copy(self.system.part[1].f),
+                                    [-p3m_force, 0, 0],atol=1E-10)
         self.system.actors.remove(p3m)
 
     @ut.skipIf( espressomd.has_features(["COULOMB_DEBYE_HUECKEL"]),
            "Features not available, skipping test!")
     def test_dh(self):
-        dh_params = dict(bjerrum_length=1.0,
+        dh_params = dict(prefactor=1.0,
                          kappa=2.0,
                          r_cut=2.0)
         test_DH = generate_test_for_class(
@@ -100,8 +100,8 @@ class ElectrostaticInteractionsTests(ut.TestCase):
             electrostatics.DH,
             dh_params)
         dh = espressomd.electrostatics.DH(
-            bjerrum_length=dh_params[
-                'bjerrum_length'],
+            prefactor=dh_params[
+                'prefactor'],
                                            kappa=dh_params['kappa'],
                                            r_cut=dh_params['r_cut'])
         self.system.actors.add(dh)
@@ -124,12 +124,12 @@ class ElectrostaticInteractionsTests(ut.TestCase):
             u_dh_core[i] = self.system.analysis.energy()['coulomb']
             f_dh_core[i] = self.system.part[0].f[0]
 
-        self.assertTrue(np.allclose(u_dh_core,
+        np.testing.assert_allclose(u_dh_core,
                                     u_dh,
-                                    atol=1e-7))
-        self.assertTrue(np.allclose(f_dh_core,
+                                    atol=1e-7)
+        np.testing.assert_allclose(f_dh_core,
                                     -f_dh,
-                                    atol=1e-2))
+                                    atol=1e-2)
         self.system.actors.remove(dh)
 
 
diff --git a/testsuite/ewald_gpu.py b/testsuite/ewald_gpu.py
index 6b35a5963..41a552a77 100644
--- a/testsuite/ewald_gpu.py
+++ b/testsuite/ewald_gpu.py
@@ -32,7 +32,7 @@ class ewald_GPU_test(ut.TestCase):
 
         es = espressomd.System()
         test_params = {}
-        test_params["bjerrum_length"] = 2
+        test_params["prefactor"] = 2
         test_params["num_kx"] = 2
         test_params["num_ky"] = 2
         test_params["num_kz"] = 2
diff --git a/testsuite/exclusions.py b/testsuite/exclusions.py
index 145767ed8..e0f7daebe 100644
--- a/testsuite/exclusions.py
+++ b/testsuite/exclusions.py
@@ -127,7 +127,7 @@ class Exclusions(ut.TestCase):
         self.s.part.add(id=1, pos=[1, 0, 0], type=0, q=-1.)
 
         # Small alpha means large short-range contribution
-        self.s.actors.add(P3M(bjerrum_length=1, r_cut=3.0, accuracy=1e-3,
+        self.s.actors.add(P3M(prefactor=1, r_cut=3.0, accuracy=1e-3,
                                   mesh=32, cao=7, alpha=0.1, tune=False))
 
         # Only short-range part of the coulomb energy
diff --git a/testsuite/icc.py b/testsuite/icc.py
index 2b31debe4..977ad6194 100644
--- a/testsuite/icc.py
+++ b/testsuite/icc.py
@@ -55,7 +55,7 @@ class test_icc(ut.TestCase):
         S.part.add(pos=[b2, b2, b2 + q_dist / 2], q=-q_test, fix=[1, 1, 1])
 
         # Actors
-        p3m = P3M(bjerrum_length=1, mesh=32, cao=7, accuracy=1e-5)
+        p3m = P3M(prefactor=1, mesh=32, cao=7, accuracy=1e-5)
         icc = ICC(
             n_icc=nicc_tot,
             convergence=1e-6,
diff --git a/testsuite/mmm1d.py b/testsuite/mmm1d.py
index 7b16c55d3..f4349e12d 100644
--- a/testsuite/mmm1d.py
+++ b/testsuite/mmm1d.py
@@ -84,7 +84,7 @@ class ElectrostaticInteractionsTests(ut.TestCase):
         return func
 
     if espressomd.has_features("ELECTROSTATICS", "PARTIAL_PERIODIC"):
-        test_mmm1d = generateTestForElectrostaticInteraction(MMM1D, dict(bjerrum_length=2.0,
+        test_mmm1d = generateTestForElectrostaticInteraction(MMM1D, dict(prefactor=2.0,
                                                                     maxPWerror= 0.001, 
                                                                     far_switch_radius = 3, 
                                                                     tune=False))
diff --git a/testsuite/p3m_gpu.py b/testsuite/p3m_gpu.py
index f06631a6e..c34a6ddb6 100644
--- a/testsuite/p3m_gpu.py
+++ b/testsuite/p3m_gpu.py
@@ -31,8 +31,9 @@ class P3MGPU_test(ut.TestCase):
         from espressomd.electrostatics import P3MGPU
 
         es = espressomd.System()
+        es.box_l=10,10,10
         test_params = {}
-        test_params["bjerrum_length"] = 2
+        test_params["prefactor"] = 2
         test_params["cao"] = 2
         test_params["r_cut"] = 0.9
         test_params["accuracy"] = 1e-1
diff --git a/testsuite/scafacos_dipoles_1d_2d.py b/testsuite/scafacos_dipoles_1d_2d.py
index 8318294ff..1b914e7e9 100644
--- a/testsuite/scafacos_dipoles_1d_2d.py
+++ b/testsuite/scafacos_dipoles_1d_2d.py
@@ -79,7 +79,7 @@ class Scafacos1d2d(ut.TestCase):
 
             if dim == 2:
                 scafacos = magnetostatics.Scafacos(
-                    bjerrum_length=1,
+                    prefactor=1,
                     method_name="p2nfft",
                     method_params={
                         "p2nfft_verbose_tuning": 0,
@@ -100,7 +100,7 @@ class Scafacos1d2d(ut.TestCase):
                 if dim == 1:
                     # 1d periodic in x
                     scafacos = magnetostatics.Scafacos(
-                        bjerrum_length=1,
+                        prefactor=1,
                         method_name="p2nfft",
                         method_params={
                             "p2nfft_verbose_tuning": 1,
diff --git a/testsuite/tests_common.py b/testsuite/tests_common.py
index 57f69e8fb..89a54f564 100644
--- a/testsuite/tests_common.py
+++ b/testsuite/tests_common.py
@@ -63,9 +63,15 @@ def params_match(inParams, outParams):
         if k not in outParams:
             print(k, "missing from returned parameters")
             return False
-        if outParams[k] != inParams[k]:
-            print("Mismatch in parameter ", k, inParams[k], outParams[k])
-            return False
+        if type(inParams[k]) == float:
+          if abs(outParams[k] -inParams[k])>=1E-14:
+              print("Mismatch in parameter ", k, inParams[k], outParams[k],type(inParams[k]),type(outParams[k]),abs(inParams[k]-outParams[k]))
+              return False
+        else:
+          if outParams[k] !=inParams[k]:
+              print("Mismatch in parameter ", k, inParams[k], outParams[k],type(inParams[k]),type(outParams[k]))
+              return False
+
 
     return True
 
['doc/sphinx/io.rst', 'testsuite/ek_eof_one_species_z.py', 'samples/coulomb_debye_hueckel.py', 'src/core/interaction_data.cpp', 'src/core/interaction_data.hpp', 'testsuite/ek_eof_one_species_x_nonlinear.py', 'doc/sphinx/installation.rst', 'testsuite/coulomb_mixed_periodicity.py', 'testsuite/ek_eof_one_species_y.py', 'testsuite/ek_eof_one_species_y_nonlinear.py', 'samples/store_properties.py', 'src/core/p3m.hpp', 'src/python/espressomd/electrostatics.pxd', 'testsuite/analyze_energy.py', 'testsuite/ewald_gpu.py', 'testsuite/icc.py', 'samples/ekboundaries.py', 'testsuite/tests_common.py', 'samples/visualization_mmm2d.py', 'src/core/p3m.cpp', 'testsuite/elc_vs_mmm2d_nonneutral.py', 'src/core/elc.hpp', 'src/core/electrokinetics_cuda.cu', 'src/python/espressomd/electrokinetics.pyx', 'testsuite/coulomb_tuning.py', 'src/core/fd-electrostatics_cuda.cu', 'testsuite/scafacos_dipoles_1d_2d.py', 'doc/sphinx/inter.rst', 'src/core/p3m-dipolar.cpp', 'testsuite/mmm1d.py', 'samples/electrophoresis.py', 'testsuite/ek_eof_one_species_z_nonlinear.py', 'src/core/integrate.cpp', 'testsuite/coulomb_cloud_wall.py', 'samples/debye_hueckel.py', 'src/core/electrokinetics.hpp', 'src/python/espressomd/scafacos.pyx', 'testsuite/ek_eof_one_species_x.py', 'testsuite/p3m_gpu.py', 'src/core/p3m-dipolar.hpp', 'testsuite/electrostaticInteractions.py', 'samples/billard.py', 'src/python/espressomd/magnetostatics.pxd', 'testsuite/dawaanr-and-dds-gpu.py', 'src/python/espressomd/electrokinetics.pxd', 'src/python/espressomd/electrostatics.pyx', 'src/core/initialize.cpp', 'testsuite/exclusions.py', 'testsuite/elc_vs_mmm2d_neutral.py', 'src/core/energy.cpp', 'samples/minimal-charged-particles.py', 'src/core/statistics.cpp', 'src/core/fd-electrostatics.hpp', 'testsuite/coulomb_cloud_wall_duplicated.py', 'src/core/unit_tests/List_test.cpp', 'samples/p3m.py', 'samples/save_checkpoint.py', 'src/core/maggs.hpp', 'src/core/maggs.cpp', 'testsuite/dipolar_mdlc_p3m_scafacos_p2nfft.py', 'src/python/espressomd/magnetostatics.pyx']