diff --git a/chapman_enskog/chapman_enskog.py b/chapman_enskog/chapman_enskog.py
index 75ed83341b5b9dccd904c692c3bedd9a26bda6bd..9ca4509f881977eb3ba13960cee0a50357e37360 100644
--- a/chapman_enskog/chapman_enskog.py
+++ b/chapman_enskog/chapman_enskog.py
@@ -331,8 +331,8 @@ def substitute_collision_operator_moments(expr, lb_moment_computation, collision
     return expr.subs(subs_dict)
 
 
-def take_moments(eqn, pdf_to_moment_name=(('f', '\Pi'), ('\Omega f', '\\Upsilon')), velocity_name='c', max_expansion=5,
-                 use_one_neighborhood_aliasing=False):
+def take_moments(eqn, pdf_to_moment_name=(('f', '\\Pi'), ('\\Omega f', '\\Upsilon')), velocity_name='c',
+                 max_expansion=5, use_one_neighborhood_aliasing=False):
 
     pdf_symbols = [tuple(expanded_symbol(name, superscript=i) for i in range(max_expansion))
                    for name, _ in pdf_to_moment_name]
@@ -469,7 +469,7 @@ def remove_error_terms(expr):
 # ----------------------------------------------------------------------------------------------------------------------
 
 
-def get_taylor_expanded_lb_equation(pdf_symbol_name="f", pdfs_after_collision_operator="\Omega f", velocity_name="c",
+def get_taylor_expanded_lb_equation(pdf_symbol_name="f", pdfs_after_collision_operator="\\Omega f", velocity_name="c",
                                     dim=3, taylor_order=2, shift=True):
     dim_labels = [sp.Rational(i, 1) for i in range(dim)]
 
@@ -502,8 +502,8 @@ def get_taylor_expanded_lb_equation(pdf_symbol_name="f", pdfs_after_collision_op
 
 
 def chapman_enskog_ansatz(equation, time_derivative_orders=(1, 3), spatial_derivative_orders=(1, 2),
-                          pdfs=(['f', 0, 3], ['\Omega f', 1, 3]), commutative=True):
-    """Uses a Chapman Enskog Ansatz to expand given equation.
+                          pdfs=(['f', 0, 3], ['\\Omega f', 1, 3]), commutative=True):
+    r"""Uses a Chapman Enskog Ansatz to expand given equation.
 
     Args:
         equation: equation to expand
diff --git a/continuous_distribution_measures.py b/continuous_distribution_measures.py
index 39b173e5a3c724d64c375dff9e0e6f8ecbbf0459..3edee60ebeefa5e71afd7335d4d5efda3b07081e 100644
--- a/continuous_distribution_measures.py
+++ b/continuous_distribution_measures.py
@@ -11,19 +11,22 @@ from pystencils.sympyextensions import complete_the_squares_in_exp
 
 @memorycache()
 def moment_generating_function(generating_function, symbols, symbols_in_result):
-    """
+    r"""
     Computes the moment generating function of a probability distribution. It is defined as:
 
     .. math ::
         F[f(\mathbf{x})](\mathbf{t}) = \int e^{<\mathbf{x}, \mathbf{t}>} f(x)\; dx
 
-    :param generating_function: sympy expression
-    :param symbols: a sequence of symbols forming the vector x
-    :param symbols_in_result: a sequence forming the vector t
-    :return: transformation result F: an expression that depends now on symbols_in_result
-             (symbols have been integrated out)
+    Args:
+        generating_function: sympy expression
+        symbols: a sequence of symbols forming the vector x
+        symbols_in_result: a sequence forming the vector t
+
+    Returns:
+        transformation result F: an expression that depends now on symbols_in_result
+        (symbols have been integrated out)
 
-    .. note::
+    Note:
          This function uses sympys symbolic integration mechanism, which may not work or take a large
          amount of time for some functions.
          Therefore this routine does some transformations/simplifications on the function first, which are
diff --git a/methods/conservedquantitycomputation.py b/methods/conservedquantitycomputation.py
index 84e2df2c42d9cbbea65cc4f41a2caf1c887a2e13..53958b269b7f8c1ab86c47f6862e7e64adafd708 100644
--- a/methods/conservedquantitycomputation.py
+++ b/methods/conservedquantitycomputation.py
@@ -5,7 +5,7 @@ from pystencils import Field, Assignment, AssignmentCollection
 
 
 class AbstractConservedQuantityComputation(abc.ABC):
-    """
+    r"""
 
     This class defines how conserved quantities are computed as functions of the pdfs.
     Conserved quantities are used for output and as input to the equilibrium in the collision step
@@ -221,7 +221,7 @@ class DensityVelocityComputation(AbstractConservedQuantityComputation):
 
 def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs,
                                                     symbolic_zeroth_moment, symbolic_first_moments):
-    """
+    r"""
     Returns an equation system that computes the zeroth and first order moments with the least amount of operations
 
     The first equation of the system is equivalent to
@@ -274,7 +274,7 @@ def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs,
 
 
 def divide_first_order_moments_by_rho(assignment_collection, dim):
-    """
+    r"""
     Assumes that the equations of the passed equation collection are the following
         - rho = f_0  + f_1 + ...
         - u_0 = ...
@@ -290,7 +290,7 @@ def divide_first_order_moments_by_rho(assignment_collection, dim):
 
 
 def add_density_offset(assignment_collection, offset=sp.Rational(1, 1)):
-    """
+    r"""
     Assumes that first equation is the density (zeroth moment). Changes the density equations by adding offset to it.
     """
     old_eqs = assignment_collection.main_assignments
diff --git a/moments.py b/moments.py
index 188fd23bdd6e4022e55e1887b511203c8ab4b4ad..1460ad164e260b226c6c3d781cc7645583c1ee4f 100644
--- a/moments.py
+++ b/moments.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-"""
+r"""
 Module Overview
 ~~~~~~~~~~~~~~~
 
@@ -285,7 +285,7 @@ is_shear_moment.shear_moments = set([c[0] * c[1] for c in itertools.combinations
 
 @memorycache(maxsize=512)
 def discrete_moment(func, moment, stencil):
-    """
+    r"""
     Computes discrete moment of given distribution function
 
     .. math ::
@@ -340,15 +340,17 @@ def moment_matrix(moments, stencil):
 
 
 def gram_schmidt(moments, stencil, weights=None):
-    """
+    r"""
     Computes orthogonal set of moments using the method by Gram-Schmidt
 
-    :param moments: sequence of moments, either in tuple or polynomial form
-    :param stencil: stencil as sequence of directions
-    :param weights: optional weights, that define the scalar product which is used for normalization.
-                    Scalar product :math:`< a,b > = \sum a_i b_i w_i` with weights :math:`w_i`.
-                    Passing no weights sets all weights to 1.
-    :return: set of orthogonal moments in polynomial form
+    Args:
+        moments: sequence of moments, either in tuple or polynomial form
+        stencil: stencil as sequence of directions
+        weights: optional weights, that define the scalar product which is used for normalization.
+                 Scalar product :math:`< a,b > = \sum a_i b_i w_i` with weights :math:`w_i`.
+                 Passing no weights sets all weights to 1.
+    Returns:
+        set of orthogonal moments in polynomial form
     """
     if weights is None:
         weights = sp.eye(len(stencil))
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 6c17ba5a429d077eb57d18445b15ba6fbf02a1ea..48b561e2b551d4a83bd8d9a7774215489f7306d0 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -186,7 +186,7 @@ def separate_into_bulk_and_interface(free_energy):
 
 
 def analytic_interface_profile(x, interface_width=interface_width_symbol):
-    """Analytic expression for a 1D interface normal to x with given interface width.
+    r"""Analytic expression for a 1D interface normal to x with given interface width.
 
     The following doctest shows that the returned analytical solution is indeed a solution of the ODE that we
     get from the condition :math:`\mu_0 = 0` (thermodynamic equilibrium) for a situation with only a single order
diff --git a/phasefield/cahn_hilliard_lbm.py b/phasefield/cahn_hilliard_lbm.py
index 8481e4ab08580d0b0835422f587a55463b378c4f..02dd3ebc5a88d783b4c70517967e55ca4aaf70c1 100644
--- a/phasefield/cahn_hilliard_lbm.py
+++ b/phasefield/cahn_hilliard_lbm.py
@@ -6,7 +6,7 @@ from lbmpy.maxwellian_equilibrium import get_weights
 
 
 def cahn_hilliard_lb_method(stencil, mu, relaxation_rate=sp.Symbol("omega"), gamma=1):
-    """Returns LB equilibrium that solves the Cahn Hilliard equation.
+    r"""Returns LB equilibrium that solves the Cahn Hilliard equation.
 
     ..math ::
 
diff --git a/relaxationrates.py b/relaxationrates.py
index d298a2b163d88f51bc0f5ac4d50b10ab0c5b7f55..15c6b40a588612d68fd3a3823e60a224afb412db 100644
--- a/relaxationrates.py
+++ b/relaxationrates.py
@@ -3,12 +3,12 @@ from lbmpy.moments import is_shear_moment, get_order
 
 
 def relaxation_rate_from_lattice_viscosity(nu):
-    """Computes relaxation rate from lattice viscosity: :math:`\omega = \frac{1}{3\nu_L + \frac{1}{2}}`"""
+    r"""Computes relaxation rate from lattice viscosity: :math:`\omega = \frac{1}{3\nu_L + \frac{1}{2}}`"""
     return 2 / (6 * nu + 1)
 
 
 def lattice_viscosity_from_relaxation_rate(omega):
-    """Computes lattice viscosity from relaxation rate:
+    r"""Computes lattice viscosity from relaxation rate:
     :math:`\nu_L=\frac{1}{3}\left(\frac{1}{\omega}-\frac{1}{2}\right)`"""
     return (2 - omega) / (6 * omega)