changing ext_force_density on lb_gpu kills the velocity field
diff --git a/src/python/espressomd/lb.pxd b/src/python/espressomd/lb.pxd
index f50b9b970..161dd2f94 100644
--- a/src/python/espressomd/lb.pxd
+++ b/src/python/espressomd/lb.pxd
@@ -195,7 +195,6 @@ cdef inline python_lbfluid_set_gamma_even(gamma_even):
     lb_lbfluid_set_gamma_even(c_gamma_even)
 
 ###############################################
-
 cdef inline python_lbfluid_set_ext_force_density(p_ext_force_density, p_agrid, p_tau):
 
     cdef Vector3d c_ext_force_density
@@ -254,21 +253,14 @@ cdef inline python_lbfluid_get_gamma(p_gamma):
         p_gamma = c_gamma
 
 
-cdef inline python_lbfluid_get_ext_force_density(p_ext_force_density, p_agrid, p_tau):
-
+cdef inline Vector3d python_lbfluid_get_ext_force_density(p_agrid, p_tau):
     cdef Vector3d c_ext_force_density
     # call c-function
     c_ext_force_density = lb_lbfluid_get_ext_force_density()
     # unit conversion LB -> MD
-    p_ext_force_density[
-        0] = c_ext_force_density[
-            0] / p_agrid / p_agrid / p_tau / p_tau
-    p_ext_force_density[
-        1] = c_ext_force_density[
-            1] / p_agrid / p_agrid / p_tau / p_tau
-    p_ext_force_density[
-        2] = c_ext_force_density[
-            2] / p_agrid / p_agrid / p_tau / p_tau
+    for i in range(3):
+        c_ext_force_density[i] /= p_agrid * p_agrid * p_tau * p_tau
+    return c_ext_force_density
 
 cdef inline void python_lbnode_set_velocity(Vector3i node, Vector3d velocity):
     cdef double inv_lattice_speed = lb_lbfluid_get_tau() / lb_lbfluid_get_agrid()
diff --git a/src/python/espressomd/lb.pyx b/src/python/espressomd/lb.pyx
index 606a0d47b..a6d4918d4 100644
--- a/src/python/espressomd/lb.pyx
+++ b/src/python/espressomd/lb.pyx
@@ -169,8 +169,7 @@ cdef class HydrodynamicInteraction(Actor):
             raise Exception("lb_lbfluid_set_agrid error")
 
         if not self._params["ext_force_density"] == default_params["ext_force_density"]:
-            if python_lbfluid_get_ext_force_density(self._params["ext_force_density"], self._params["agrid"], self._params["tau"]):
-                raise Exception("lb_lbfluid_set_ext_force_density error")
+            self._params["ext_force_density"] = self.ext_force_density
 
         return self._params
 
@@ -261,6 +260,19 @@ cdef class HydrodynamicInteraction(Actor):
         def __set__(self, value):
             raise NotImplementedError
 
+    property ext_force_density:
+        def __get__(self):
+            cdef Vector3d res
+            res = python_lbfluid_get_ext_force_density(
+                self._params["agrid"], self._params["tau"])
+            return make_array_locked(res)
+
+        def __set__(self, ext_force_density):
+            python_lbfluid_set_ext_force_density(
+    ext_force_density,
+     self._params["agrid"],
+     self._params["tau"])
+
     def nodes(self):
         """Provides a generator for iterating over all lb nodes"""
         
@@ -327,9 +339,9 @@ IF CUDA:
             length = positions.shape[0]
             velocities = np.empty_like(positions)
             if three_point:
-                quadratic_velocity_interpolation( < double * >np.PyArray_GETPTR2(positions, 0, 0), < double * >np.PyArray_GETPTR2(velocities, 0, 0), length)
+                quadratic_velocity_interpolation(< double * >np.PyArray_GETPTR2(positions, 0, 0), < double * >np.PyArray_GETPTR2(velocities, 0, 0), length)
             else:
-                linear_velocity_interpolation( < double * >np.PyArray_GETPTR2(positions, 0, 0), < double * >np.PyArray_GETPTR2(velocities, 0, 0), length)
+                linear_velocity_interpolation(< double * >np.PyArray_GETPTR2(positions, 0, 0), < double * >np.PyArray_GETPTR2(velocities, 0, 0), length)
             return velocities * lb_lbfluid_get_lattice_speed()
 
 cdef class LBFluidRoutines(object):
diff --git a/testsuite/python/lb.py b/testsuite/python/lb.py
index 90c3dbdad..8917412fd 100644
--- a/testsuite/python/lb.py
+++ b/testsuite/python/lb.py
@@ -221,6 +221,17 @@ class TestLB(object):
         self.assertAlmostEqual(self.lbf[0, 0, 0].density, density, delta=1e-4)
 
         self.assertEqual(self.lbf[3, 2, 1].index, (3, 2, 1))
+        ext_force_density = [0.1, 0.2, 1.2]
+        self.lbf[1, 2, 3].velocity = v_fluid
+        self.lbf.ext_force_density = ext_force_density
+        np.testing.assert_allclose(
+            np.copy(self.lbf[1, 2, 3].velocity),
+            v_fluid,
+            atol=1e-4)
+        np.testing.assert_allclose(
+            np.copy(self.lbf.ext_force_density),
+            ext_force_density,
+            atol=1e-4)
 
     def test_parameter_change_without_seed(self):
         self.system.actors.clear()
['src/python/espressomd/lb.pxd', 'src/python/espressomd/lb.pyx', 'testsuite/python/lb.py']