diff --git a/creationfunctions.py b/creationfunctions.py
index f94cbfa70d5849b2c9ce997dfaba4acfe38f3c82..17526adde5c07e0ce2da539e0fb0916a0fb4c6cb 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -11,7 +11,7 @@ General:
 
 - ``stencil='D2Q9'``: stencil name e.g. 'D2Q9', 'D3Q19'. See :func:`pystencils.stencils.get_stencil` for details
 - ``method='srt'``: name of lattice Boltzmann method. This determines the selection and relaxation pattern of
-   moments/cumulants, i.e. which moment/cumulant basis is chosen, and which of the basis vectors are relaxed together
+  moments/cumulants, i.e. which moment/cumulant basis is chosen, and which of the basis vectors are relaxed together
     - ``srt``: single relaxation time (:func:`lbmpy.methods.create_srt`)
     - ``trt``: two relaxation time, first relaxation rate is for even moments and determines the viscosity (as in SRT),
       the second relaxation rate is used for relaxing odd moments, and controls the bulk viscosity.
@@ -51,7 +51,7 @@ General:
   fields. In each timestep the corresponding quantities are written to the given fields.
 - ``velocity_input``: symbolic field where the velocities are read from (for advection diffusion LBM)
 - ``density_input``: symbolic field or field access where to read density from. When passing this parameter,
-   ``velocity_input`` has to be passed as well
+  ``velocity_input`` has to be passed as well
 - ``kernel_type``: supported values: 'stream_pull_collide' (default), 'collide_only'
 
 
diff --git a/methods/creationfunctions.py b/methods/creationfunctions.py
index e9e0ae5b70c2b5ebae1fec21fc1a53963817f101..859b03a4780b13e0dbfd6bcf28acbf9828222317 100644
--- a/methods/creationfunctions.py
+++ b/methods/creationfunctions.py
@@ -315,8 +315,8 @@ def create_trt_kbc(dim, shear_relaxation_rate, higher_order_relaxation_rate, met
     energy_transport_tensor = list(poly_repr([a for a in moments_of_order(3, dim, True)
                                               if 3 not in a]))
 
-    explicitly_defined = set(rho + velocity + shear_tensor_off_diagonal +
-                             shear_tensor_diagonal + energy_transport_tensor)
+    explicitly_defined = set(rho + velocity + shear_tensor_off_diagonal
+                             + shear_tensor_diagonal + energy_transport_tensor)
     rest = list(set(poly_repr(moments_up_to_component_order(2, dim))) - explicitly_defined)
     assert len(rest) + len(explicitly_defined) == 3**dim
 
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 20b7171e14b6e2bc9749aea4176396cb2f293f50..ffbb4336fbae5f17eea4aca4c3505b83f8c2c2ec 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -167,8 +167,8 @@ def free_energy_functional_n_phases(num_phases=None, surface_tensions=symmetric_
         n = num_phases - 1
         if k == l:
             assert surface_tensions(l, l) == 0
-        return 3 / sp.sqrt(2) * interface_width * (surface_tensions(k, n) +
-                                                   surface_tensions(l, n) - surface_tensions(k, l))
+        return 3 / sp.sqrt(2) * interface_width * (surface_tensions(k, n)
+                                                   + surface_tensions(l, n) - surface_tensions(k, l))
 
     def bulk_term(i, j):
         return surface_tensions(i, j) / 2 * (f1(phi[i]) + f1(phi[j]) - f2(phi[i] + phi[j]))
diff --git a/phasefield/contact_angle_circle_fitting.py b/phasefield/contact_angle_circle_fitting.py
index 52681f6a96423f09b0ddb530a44923cfb9950351..070193efe7f3bab92954cff3a56810d31e68b83d 100644
--- a/phasefield/contact_angle_circle_fitting.py
+++ b/phasefield/contact_angle_circle_fitting.py
@@ -23,12 +23,12 @@ def circle_intersections(midpoint0, midpoint1, radius0, radius1):
         yy = -(-b2 * k1 + b1 * k2) / (k1 - k2)
         return [(xx, yy)]
     elif np.abs(radius1 - radius0) < midpoint_distance < radius0 + radius1:
-        xx1 = (-b2 * k2 + x1 + k2 * y1 - np.sqrt(-b2 ** 2 + radius1 ** 2 + k2 ** 2 * radius1 ** 2 - 2 * b2 *
-               k2 * x1 - k2 ** 2 * x1 ** 2 + 2 * b2 * y1 + 2 * k2 * x1 * y1 - y1 ** 2)) / (1 + k2 ** 2)
+        xx1 = (-b2 * k2 + x1 + k2 * y1 - np.sqrt(-b2 ** 2 + radius1 ** 2 + k2 ** 2 * radius1 ** 2 - 2 * b2
+               * k2 * x1 - k2 ** 2 * x1 ** 2 + 2 * b2 * y1 + 2 * k2 * x1 * y1 - y1 ** 2)) / (1 + k2 ** 2)
         yy1 = k2 * xx1 + b2
 
-        xx2 = (-b2 * k2 + x1 + k2 * y1 + np.sqrt(-b2 ** 2 + radius1 ** 2 + k2 ** 2 * radius1 ** 2 - 2 * b2 *
-               k2 * x1 - k2 ** 2 * x1 ** 2 + 2 * b2 * y1 + 2 * k2 * x1 * y1 - y1 ** 2)) / (1 + k2 ** 2)
+        xx2 = (-b2 * k2 + x1 + k2 * y1 + np.sqrt(-b2 ** 2 + radius1 ** 2 + k2 ** 2 * radius1 ** 2 - 2 * b2
+               * k2 * x1 - k2 ** 2 * x1 ** 2 + 2 * b2 * y1 + 2 * k2 * x1 * y1 - y1 ** 2)) / (1 + k2 ** 2)
         yy2 = k2 * xx2 + b2
 
         return [(xx1, yy1), (xx2, yy2)]
diff --git a/phasefield/high_density_ratio_model.py b/phasefield/high_density_ratio_model.py
index 5d11f4dcc08009bbc19b06afdef3516611b8bb87..819808f1689892752de4b0a4a0642a20f02d7b71 100644
--- a/phasefield/high_density_ratio_model.py
+++ b/phasefield/high_density_ratio_model.py
@@ -42,10 +42,10 @@ def free_energy_high_density_ratio(eos, density, density_gas, density_liquid, c_
     def f(c):
         return c ** 2 * (1 - c) ** 2
 
-    f_bulk = (lambdas[0] / 2 * (psi_eos - psi_0) +
-              lambdas[1] / 2 * f(c_liquid_1) +
-              lambdas[2] / 2 * f(c_liquid_2))
-    f_interface = (kappas[0] / 2 * Diff(density) ** 2 +
-                   kappas[1] / 2 * Diff(c_liquid_1) ** 2 +
-                   kappas[2] / 2 * Diff(c_liquid_2) ** 2)
+    f_bulk = (lambdas[0] / 2 * (psi_eos - psi_0)
+              + lambdas[1] / 2 * f(c_liquid_1)
+              + lambdas[2] / 2 * f(c_liquid_2))
+    f_interface = (kappas[0] / 2 * Diff(density) ** 2
+                   + kappas[1] / 2 * Diff(c_liquid_1) ** 2
+                   + kappas[2] / 2 * Diff(c_liquid_2) ** 2)
     return f_bulk + f_interface
diff --git a/phasefield/n_phase_boyer.py b/phasefield/n_phase_boyer.py
index a57b98aa40d2a49cb16320ad4c5ea2a9f297a335..bdcac888c9c2f90a2b47b5e20bed0ec3b161ece9 100644
--- a/phasefield/n_phase_boyer.py
+++ b/phasefield/n_phase_boyer.py
@@ -123,9 +123,9 @@ def capital_f0(c, surface_tension, f=lambda c: c ** 2 * (1 - c) ** 2):
 def free_energy(c, epsilon, surface_tensions, stabilization_factor):
     alpha, _ = diffusion_coefficients(surface_tensions)
 
-    capital_f = (capital_f0(c, surface_tensions) +
-                 correction_g(c, surface_tensions) +
-                 stabilization_factor * stabilization_term(c, alpha))
+    capital_f = (capital_f0(c, surface_tensions)
+                 + correction_g(c, surface_tensions)
+                 + stabilization_factor * stabilization_term(c, alpha))
 
     def f(x):
         return x ** 2 * (1 - x) ** 2
@@ -200,9 +200,9 @@ def capital_gamma(sigma, i, index_tuple):
     assert tuple(sorted(index_tuple)) == index_tuple
     j, k, m = index_tuple
     alpha, gamma = diffusion_coefficients(sigma)
-    return -6 * (alpha[i, j] * (sigma[j, k] + sigma[j, m]) +
-                 alpha[i, k] * (sigma[j, k] + sigma[k, m]) +
-                 alpha[i, m] * (sigma[j, m] + sigma[k, m]) - gamma[i])
+    return -6 * (alpha[i, j] * (sigma[j, k] + sigma[j, m])
+                 + alpha[i, k] * (sigma[j, k] + sigma[k, m])
+                 + alpha[i, m] * (sigma[j, m] + sigma[k, m]) - gamma[i])
 
 
 def capital_lambda(surface_tensions, index_tuple):
diff --git a/phasefield/phasefieldstep_direct.py b/phasefield/phasefieldstep_direct.py
index 29f3a19600f2082e6648a88524e22b6fcae1f9d8..3937e37f8873cfb616178aca60844194201442aa 100644
--- a/phasefield/phasefieldstep_direct.py
+++ b/phasefield/phasefieldstep_direct.py
@@ -119,8 +119,8 @@ class PhaseFieldStepDirect:
         # Sync functions
         self.phi_sync = dh.synchronization_function([self.phi_field.name])
         self.mu_sync = dh.synchronization_function([self.mu_field.name])
-        self.pdf_sync = dh.synchronization_function([self.hydro_pdfs[0].name] +
-                                                    [src.name for src, _ in self.ch_pdfs])
+        self.pdf_sync = dh.synchronization_function([self.hydro_pdfs[0].name]
+                                                    + [src.name for src, _ in self.ch_pdfs])
 
         self.reset()