diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 04141d74035bcf96ec5e30aa9c14b634e0513423..93a6234a898ed36aa1cc3c1ff4b30bef9f6da6df 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -343,22 +343,9 @@ def pressure_tensor_interface_component_new(free_energy, order_parameters, dim,
     return result
 
 
-def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix=None,
+def pressure_tensor_from_free_energy(free_energy, order_parameters, dim,
                                      include_bulk=True, include_interface=True):
-
-    def transform(f, matrix, from_symbols, to_symbols):
-        transformed_vals = matrix * sp.Matrix(to_symbols)
-        substitutions = {a: b for a, b in zip(from_symbols, transformed_vals)}
-        return f.subs(substitutions)
-
-    c = sp.symbols("C_:{}".format(len(order_parameters)))
-
-    if transformation_matrix:
-        phi_to_c = partial(transform, matrix=transformation_matrix, from_symbols=order_parameters, to_symbols=c)
-        free_energy = expand_diff_full(phi_to_c(free_energy))
-        op = c
-    else:
-        op = order_parameters
+    op = order_parameters
 
     def get_entry(i, j):
         p_if = pressure_tensor_interface_component(free_energy, op, dim, i, j) if include_interface else 0
@@ -369,11 +356,6 @@ def pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transfo
         return sp.expand(p_if + p_b)
 
     result = sp.Matrix(dim, dim, get_entry)
-
-    if transformation_matrix:
-        c_to_phi = partial(transform, matrix=transformation_matrix.inv(), from_symbols=c, to_symbols=order_parameters)
-        result = result.applyfunc(lambda e: expand_diff_full(c_to_phi(e)))
-
     return result
 
 
diff --git a/phasefield/kerneleqs.py b/phasefield/kerneleqs.py
index fce9bbede63e26fc0d6aa0c4123b087b4c10979d..75f421d9238e7a186555f8824ad85709baecb13d 100644
--- a/phasefield/kerneleqs.py
+++ b/phasefield/kerneleqs.py
@@ -29,10 +29,10 @@ def force_kernel_using_mu(force_field, phi_field, mu_field, dx=1, discretization
 
 
 def pressure_tensor_kernel(free_energy, order_parameters, phi_field, pressure_tensor_field,
-                           transformation_matrix=None, dx=1, discretization='standard'):
+                           dx=1, discretization='standard'):
     dim = phi_field.spatial_dimensions
 
-    p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim, transformation_matrix)
+    p = pressure_tensor_from_free_energy(free_energy, order_parameters, dim)
 
     p = p.subs({op: phi_field(i) for i, op in enumerate(order_parameters)})
     index_map = symmetric_tensor_linearization(dim)
diff --git a/phasefield/phasefieldstep.py b/phasefield/phasefieldstep.py
index 706e0c0776f581bcb390b072d85c4cddc3e01ca7..31b070d9c0f0733ecb995b571797cc5dbbb213cc 100644
--- a/phasefield/phasefieldstep.py
+++ b/phasefield/phasefieldstep.py
@@ -6,7 +6,7 @@ import numpy as np
 from lbmpy.lbstep import LatticeBoltzmannStep
 from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
 from lbmpy.phasefield.kerneleqs import mu_kernel, CahnHilliardFDStep, pressure_tensor_kernel, \
-    force_kernel_using_pressure_tensor, pressure_tensor_kernel_pbs
+    force_kernel_using_pressure_tensor
 from pystencils import create_kernel, create_data_handling
 from lbmpy.phasefield.analytical import chemical_potentials_from_free_energy, symmetric_tensor_linearization
 from pystencils.boundaries.boundaryhandling import FlagInterface
@@ -98,20 +98,9 @@ class PhaseFieldStep:
         self.phi_sync = data_handling.synchronization_function([self.phi_field_name], target=target)
         self.mu_eqs = mu_kernel(F, phi, self.phi_field, self.mu_field, dx)
 
-        pbs = False
-        if pbs:
-            self.pbs_field_name = name + "_pbs"
-            self.pbs_field = dh.add_array(self.pbs_field_name, gpu=gpu)
-            self.pressure_tensor_eqs = pressure_tensor_kernel_pbs(self.free_energy, order_parameters, self.phi_field,
-                                                                  self.pressure_tensor_field, self.pbs_field, dx=dx,
-                                                                  density_field=None, discretization=discretization)
-            # TODO get current density! not last one
-            # TODO call post-run on hydro-lbm before computing pbs to store the latest density
-        else:
-            self.pressure_tensor_eqs = pressure_tensor_kernel(self.free_energy, order_parameters,
-                                                              self.phi_field, self.pressure_tensor_field, dx=dx,
-                                                              transformation_matrix=transformation_matrix,
-                                                              discretization=discretization)
+        self.pressure_tensor_eqs = pressure_tensor_kernel(self.free_energy, order_parameters,
+                                                          self.phi_field, self.pressure_tensor_field, dx=dx,
+                                                          discretization=discretization)
         mu_and_pressure_tensor_eqs = self.mu_eqs + self.pressure_tensor_eqs
         mu_and_pressure_tensor_eqs = apply_neumann_boundaries(mu_and_pressure_tensor_eqs)
         self.mu_and_pressure_tensor_kernel = create_kernel(sympy_cse_on_assignment_list(mu_and_pressure_tensor_eqs),