diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index 479f0ca8..d78f64b5 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -6,6 +6,8 @@ on:
       - '**'
     tags:
       - '**'
+  pull_request:
+    branches: [ master ]
 
 jobs:
   tests:
@@ -24,6 +26,14 @@ jobs:
           # Upgrade to latest meson and ninja
           sudo pip install --upgrade meson ninja
 
+      # needs fenics only to run test_poisson_nernst_planck_solver_fenics
+      - name: install_fenics
+        run: |
+          sudo apt-get install -y --no-install-recommends software-properties-common
+          sudo add-apt-repository -y ppa:fenics-packages/fenics
+          sudo apt-get update -qy
+          sudo apt-get install -y fenics python3-dolfin python3-mshr
+
       - name: build_c
         run: |
           sudo pip install .[test]
diff --git a/.github/workflows/tests_fenics.yml b/.github/workflows/tests_fenics.yml
deleted file mode 100644
index 645257ea..00000000
--- a/.github/workflows/tests_fenics.yml
+++ /dev/null
@@ -1,66 +0,0 @@
-name: Test electrochemistry examples including fenics functionality
-
-on:
-  push:
-    branches:
-      - '**'
-    tags:
-      - '**'
-
-jobs:
-  tests:
-
-    runs-on: ubuntu-latest
-
-    steps:
-      - uses: actions/checkout@v3
-        with:
-          fetch-depth: 0
-
-      - name: install_python
-        run: |
-          sudo apt-get update -qy
-          sudo apt-get install -y python3-dev python3-pip python3-venv libxml2-dev libxslt-dev zlib1g-dev
-          # Upgrade to latest meson and ninja
-          sudo pip install --upgrade meson ninja
-          sudo pip install --upgrade jupyter jupytext pandas
-
-      - name: install_fenics
-        run: |
-          sudo apt-get install -y --no-install-recommends software-properties-common
-          sudo add-apt-repository -y ppa:fenics-packages/fenics
-          sudo apt-get update -qy
-          sudo apt-get install -y fenics python3-dolfin python3-mshr
-
-      - name: install matscipy
-        run: |
-          sudo pip install .
-
-      - name: build sample ipynb and html
-        run: |
-          cd examples/electrochemistry
-                    
-          # check that all sample html files are produced successfully
-          html_files=(./*.html)
-          rm *.html
-          
-          ls -1 *.py | grep -v to_html.py | bash to_html.sh
-          
-          for f in "${html_files[@]}"; do
-            if ! test -f "$f"; then
-              echo "Sample $f creation failed."
-              exit 1
-            fi
-          done
-
-    # use tar here, see https://github.com/actions/upload-artifact/issues/38
-      - name: Pack build assets
-        run: |
-          bash -c "tar -cvf electrochemistry-sample-ipynb.tar examples/electrochemistry/*.ipynb examples/electrochemistry/*.html"
-
-      - name: Upload sample ipynb artifact
-        uses: actions/upload-artifact@v3
-        with:
-          name: electrochemistry-sample-ipynb
-          path: electrochemistry-sample-ipynb.tar
-          if-no-files-found: error
\ No newline at end of file
diff --git a/docs/applications/electrochemistry_1.py b/docs/applications/electrochemistry_1.py
index 6dad6907..edb6d9ef 100644
--- a/docs/applications/electrochemistry_1.py
+++ b/docs/applications/electrochemistry_1.py
@@ -55,10 +55,10 @@
 pnp = PoissonNernstPlanckSystem(c=c, z=z, L=L, delta_u=delta_u)
 
 # %% [markdown]
-# The method `useStandardInterfaceBC` will apply boundary conditions as shown in <a hef="#figure1">Figure 1</a>, with Dirichlet boundary conditions on the potential and zero total flux boundary conditions on ion diffusion and electromigration.
+# The method `use_standard_interface_bc` will apply boundary conditions as shown in <a hef="#figure1">Figure 1</a>, with Dirichlet boundary conditions on the potential and zero total flux boundary conditions on ion diffusion and electromigration.
 
 # %%
-pnp.useStandardInterfaceBC()
+pnp.use_standard_interface_bc()
 ui, nij, _ = pnp.solve()
 
 # %% [markdown]
@@ -188,7 +188,7 @@ def make_patch_spines_invisible(ax):
 # The only difference is the application of another set of predefined boundary conditions.
 
 # %%
-pnp.useStandardCellBC()
+pnp.use_standard_cell_bc()
 
 # %% [markdown]
 # These boundary conditions assume no Stern layer, i.e. $\lambda_S = 0$
@@ -309,7 +309,7 @@ def make_patch_spines_invisible(ax):
 # The following applies zero flux and, particularly, Robin boundary conditions as shown in <a href="#figure2">Figure 2</a>.
 
 # %%
-pnp.useSternLayerCellBC()
+pnp.use_stern_layer_cell_bc()
 
 pnp.output = True
 xij = pnp.solve()
diff --git a/docs/applications/electrochemistry_2.py b/docs/applications/electrochemistry_2.py
index b5bc6af3..d29af178 100644
--- a/docs/applications/electrochemistry_2.py
+++ b/docs/applications/electrochemistry_2.py
@@ -63,12 +63,12 @@ def make_patch_spines_invisible(ax):
     c, z, L, delta_u=delta_u, N=N,
     solver="hybr", options={'xtol':1e-12})
 
-pnp['std_interface'].useStandardInterfaceBC()
+pnp['std_interface'].use_standard_interface_bc()
 uij, nij, lamj = pnp['std_interface'].solve()
 
 # define desired system
 pnp['fenics_interface'] = PoissonNernstPlanckSystemFEniCS(c, z, L, delta_u=delta_u, N=N)
-pnp['fenics_interface'].useStandardInterfaceBC()
+pnp['fenics_interface'].use_standard_interface_bc()
 uij, nij, _ = pnp['fenics_interface'].solve()
 
 # %% [markdown]
@@ -184,7 +184,7 @@ def make_patch_spines_invisible(ax):
 pnp['std_interface_high_potential'] = PoissonNernstPlanckSystem(
     c, z, L, delta_u=delta_u,N=N,
     solver="hybr", options={'xtol':1e-14})
-pnp['std_interface_high_potential'].useStandardInterfaceBC()
+pnp['std_interface_high_potential'].use_standard_interface_bc()
 uij, nij, lamj = pnp['std_interface_high_potential'].solve()
 
 # %% [markdown]
@@ -192,7 +192,7 @@ def make_patch_spines_invisible(ax):
 
 # %%
 pnp['fenics_interface_high_potential'] = PoissonNernstPlanckSystemFEniCS(c, z, L, delta_u=delta_u, N=200)
-pnp['fenics_interface_high_potential'].useStandardInterfaceBC()
+pnp['fenics_interface_high_potential'].use_standard_interface_bc()
 uij, nij, _ = pnp['fenics_interface_high_potential'].solve()
 
 # %% [markdown]
diff --git a/matscipy/cli/electrochemistry/poisson_nernst_planck_solver_cli.py b/matscipy/cli/electrochemistry/poisson_nernst_planck_solver_cli.py
index 83c234fc..add7614b 100644
--- a/matscipy/cli/electrochemistry/poisson_nernst_planck_solver_cli.py
+++ b/matscipy/cli/electrochemistry/poisson_nernst_planck_solver_cli.py
@@ -164,6 +164,37 @@ class ArgumentDefaultsAndRawDescriptionHelpFormatter(
     logger = logging.getLogger()
     logger.setLevel(loglevel)
 
+    # use FEniCS finite element solver if available,
+    # otherwise own controlled volume scheme
+    try:
+        import fenics
+        from matscipy.electrochemistry.poisson_nernst_planck_solver_fenics \
+            import PoissonNernstPlanckSystemFEniCS as PoissonNernstPlanckSystem
+        import dolfin
+        import ffc
+        dolfin.cpp.log.set_log_level(loglevel)
+
+        if not args.outfile:
+            # fenics writes some log messages to stdout. If we pipe the output,
+            # we don't want that, hence here we have to suppress fenics logging
+            # if no output file has been specified
+            ffc.log.set_level(logging.ERROR)
+            fenics.set_log_level(logging.ERROR)
+            fenics.set_log_active(False)
+            dolfin.cpp.log.set_log_level(logging.ERROR)
+            dolfin.cpp.log.set_log_active(False)
+
+            logging.getLogger('UFL').setLevel(logging.ERROR)
+            logging.getLogger('FFC').setLevel(logging.ERROR)
+
+        logger.info("Will use FEniCS finite element solver.")
+    except ModuleNotFoundError:
+        logger.warning(
+            "No FEniCS finite element solver found,"
+            " falling back to internal controlled-volume implementation."
+            " ATTENTION: Number conservation not exact.")
+        from matscipy.electrochemistry import PoissonNernstPlanckSystem
+
     # remove all handlers
     for h in logger.handlers:
         logger.removeHandler(h)
@@ -181,22 +212,6 @@ class ArgumentDefaultsAndRawDescriptionHelpFormatter(
         fh.setLevel(loglevel)
         logger.addHandler(fh)
 
-    # use FEniCS finite element solver if available,
-    # otherwise own controlled volume scheme
-    try:
-        import fenics
-        from matscipy.electrochemistry.poisson_nernst_planck_solver_fenics \
-            import PoissonNernstPlanckSystemFEniCS as PoissonNernstPlanckSystem
-        import dolfin
-        dolfin.cpp.log.set_log_level(loglevel)
-        logger.info("Will use FEniCS finite element solver.")
-    except ModuleNotFoundError:
-        logger.warning(
-            "No FEniCS finite element solver found,"
-            " falling back to internal controlled-volume implementation."
-            " ATTENTION: Number conservation not exact.")
-        from matscipy.electrochemistry import PoissonNernstPlanckSystem
-
     # set up system
     pnp = PoissonNernstPlanckSystem(
         c=np.array(args.concentrations, dtype=float),
@@ -211,22 +226,16 @@ class ArgumentDefaultsAndRawDescriptionHelpFormatter(
         relative_permittivity=float(args.relative_permittivity))
 
     if args.boundary_conditions in ('cell-robin') and float(args.lambda_S) > 0:
-        pnp.useSternLayerCellBC()
+        pnp.use_stern_layer_cell_bc()
     elif ((args.boundary_conditions in ('cell-robin') and float(args.lambda_S) == 0)
             or args.boundary_conditions == 'cell'):
-        pnp.useStandardCellBC()
+        pnp.use_standard_cell_bc()
     elif args.boundary_conditions == 'interface':
-        pnp.useStandardInterfaceBC()
+        pnp.use_standard_interface_bc()
     else:
         raise ValueError("Boundary conditions '{}' not implemented!".format(
                          args.boundary_conditions))
 
-    if not args.outfile:
-        # fenics writes some log messages to stdout. If we pipe the output,
-        # we don't want that, hence here we have to suppress fenics logging
-        # if no output file has been specified
-        dolfin.cpp.log.set_log_active(False)
-
     pnp.solve()
 
     extra_kwargs = {}
diff --git a/matscipy/electrochemistry/continuous2discrete.py b/matscipy/electrochemistry/continuous2discrete.py
index ac174e6d..09cb27f2 100644
--- a/matscipy/electrochemistry/continuous2discrete.py
+++ b/matscipy/electrochemistry/continuous2discrete.py
@@ -289,7 +289,7 @@ def generate_structure(distribution, box=np.array([50, 50, 100]), count=100, n_g
             Z, _ = integrate.quad(d, support[-1][0], support[-1][1])
         else:  # discrete support
             support.append(np.linspace(0, box[k], n_gridpoints[k]))
-            Z = np.sum(d(support[-1])) # Normalization constant
+            Z = np.sum(d(support[-1]))  # Normalization constant
 
         logger.info("Normalizing 'distribution' {} by {}.".format(d, Z))
         normalized_distribution.append(
diff --git a/matscipy/electrochemistry/poisson_boltzmann_distribution.py b/matscipy/electrochemistry/poisson_boltzmann_distribution.py
index 415b4829..6ffa9770 100644
--- a/matscipy/electrochemistry/poisson_boltzmann_distribution.py
+++ b/matscipy/electrochemistry/poisson_boltzmann_distribution.py
@@ -86,8 +86,7 @@ def debye(c, z,
     lambda_D : float
         Debye length, sqrt( epsR*eps*R*T/(2*F^2*I) ) [m]
     """
-    I = ionic_strength(c, z)
-    return np.sqrt(relative_permittivity*vacuum_permittivity*R*T/(2.0*F**2*I))
+    return np.sqrt(relative_permittivity*vacuum_permittivity*R*T/(2.0*F**2*ionic_strength(c, z)))
 
 
 def gamma(u, T=298.15):
diff --git a/matscipy/electrochemistry/poisson_nernst_planck_solver.py b/matscipy/electrochemistry/poisson_nernst_planck_solver.py
index 0e3887d9..78f88d40 100644
--- a/matscipy/electrochemistry/poisson_nernst_planck_solver.py
+++ b/matscipy/electrochemistry/poisson_nernst_planck_solver.py
@@ -225,7 +225,6 @@ def solver_callback(self, xij, *_):
             self.convergenceStepRelative = []
             self.convergenceResidualAbsolute = []
 
-
             dxij = np.zeros(self.N)
 
         self.xij.append(xij)
@@ -257,7 +256,7 @@ def discretize(self):
         """Sets up discretization scheme and initial value"""
         # indices
         self.Ni = self.N+1
-        I = np.arange(self.Ni)
+        index_i = np.arange(self.Ni)
 
         self.logger.info('{:<{lwidth}s} {:> 8.4g}'.format(
           'discretization segments N', self.N, lwidth=self.label_width))
@@ -271,7 +270,7 @@ def discretize(self):
           'dx', self.dx, lwidth=self.label_width))
 
         # positions (scaled)
-        self.X = self.x0_scaled + I*self.dx
+        self.X = self.x0_scaled + index_i * self.dx
 
         # Boundary & initial values
 
@@ -289,9 +288,6 @@ def discretize(self):
         if self.zi0 is None:
             self.zi0 = np.kron(self.z, np.ones((self.Ni, 1))).T  # does not change
 
-
-        # self.initial_values()
-
     def initial_values(self):
         """
         Solves decoupled linear system to get initial potential distribution.
@@ -379,7 +375,7 @@ def solve(self):
         return self.uij, self.nij, self.lamj
 
     # standard sets of boundary conditions:
-    def useStandardInterfaceBC(self):
+    def use_standard_interface_bc(self):
         """Interface at left hand side and open bulk at right hand side"""
         self.boundary_conditions = []
 
@@ -392,24 +388,24 @@ def useStandardInterfaceBC(self):
             'Right hand side Dirichlet boundary condition:                              u1 = {:> 8.4g}'.format(self.u1))
 
         self.boundary_conditions.extend([
-            lambda x: self.leftPotentialDirichletBC(x, self.u0),
-            lambda x: self.rightPotentialDirichletBC(x, self.u1)])
+            lambda x: self.left_potential_dirichlet_bc(x, self.u0),
+            lambda x: self.right_potential_dirichlet_bc(x, self.u1)])
 
         for k in range(self.M):
             self.logger.info(
                 'Ion species {:02d} left hand side concentration Flux boundary condition:       j0 = {:> 8.4g}'.format(
-                    k,0))
+                    k, 0))
             self.logger.info(
                 'Ion species {:02d} right hand side concentration Dirichlet boundary condition: c1 = {:> 8.4g}'.format(
                     k, self.c_scaled[k]))
             self.boundary_conditions.extend([
-                lambda x, k=k: self.leftControlledVolumeSchemeFluxBC(x, k),
-                lambda x, k=k: self.rightDirichletBC(x, k, self.c_scaled[k])])
+                lambda x, k=k: self.left_controlled_volume_scheme_flux_bc(x, k),
+                lambda x, k=k: self.right_dirichlet_bc(x, k, self.c_scaled[k])])
         # counter-intuitive behavior of lambda in loop:
         # https://stackoverflow.com/questions/2295290/what-do-lambda-function-closures-capture
         # workaround: default parameter k=k
 
-    def useStandardCellBC(self):
+    def use_standard_cell_bc(self):
         """Interfaces at left hand side and right hand side"""
         self.boundary_conditions = []
 
@@ -421,8 +417,8 @@ def useStandardCellBC(self):
         self.logger.info('{:>{lwidth}s} u1 = {:< 8.4g}'.format(
           'Right hand side Dirichlet boundary condition', self.u1, lwidth=self.label_width))
         self.boundary_conditions.extend([
-            lambda x: self.leftPotentialDirichletBC(x, self.u0),
-            lambda x: self.rightPotentialDirichletBC(x, self.u1)])
+            lambda x: self.left_potential_dirichlet_bc(x, self.u0),
+            lambda x: self.right_potential_dirichlet_bc(x, self.u1)])
 
         N0 = self.L_scaled*self.c_scaled  # total amount of species in cell
         for k in range(self.M):
@@ -434,10 +430,10 @@ def useStandardCellBC(self):
                 N0[k], lwidth=self.label_width))
 
             self.boundary_conditions.extend([
-                lambda x, k=k: self.leftControlledVolumeSchemeFluxBC(x, k),
-                lambda x, k=k, N0=N0[k]: self.numberConservationConstraint(x, k, N0)])
+                lambda x, k=k: self.left_controlled_volume_scheme_flux_bc(x, k),
+                lambda x, k=k, N0=N0[k]: self.number_conservation_constraint(x, k, N0)])
 
-    def useSternLayerCellBC(self, implicit=False):
+    def use_stern_layer_cell_bc(self, implicit=False):
         """Interfaces at left hand side and right hand side,
         Stern layer either by prescribing linear potential regime between cell
         boundary and outer Helmholtz plane (OHP), or by applying Robin BC;
@@ -466,8 +462,8 @@ def useSternLayerCellBC(self, implicit=False):
             self.logger.info('{:>{lwidth}s} u1 + lambda_S*dudx = {:< 8.4g}'.format(
                 'Right hand side Robin boundary condition', self.u1, lwidth=self.label_width))
             self.boundary_conditions.extend([
-                lambda x: self.leftPotentialRobinBC(x, self.lambda_S_scaled, self.u0),
-                lambda x: self.rightPotentialRobinBC(x, self.lambda_S_scaled, self.u1)])
+                lambda x: self.left_potential_robin_bc(x, self.lambda_S_scaled, self.u0),
+                lambda x: self.right_potential_robin_bc(x, self.lambda_S_scaled, self.u1)])
 
         else:  # explicitly treat Stern layer via linear regime
             self.logger.info('Explicitly treating Stern layer as uniformly charged regions')
@@ -481,9 +477,8 @@ def useSternLayerCellBC(self, implicit=False):
             self.logger.info('{:>{lwidth}s} u1 = {:< 8.4g}'.format(
                 'Right hand side Dirichlet boundary condition', self.u1, lwidth=self.label_width))
             self.boundary_conditions.extend([
-                lambda x: self.leftPotentialDirichletBC(x, self.u0),
-                lambda x: self.rightPotentialDirichletBC(x, self.u1)])
-
+                lambda x: self.left_potential_dirichlet_bc(x, self.u0),
+                lambda x: self.right_potential_dirichlet_bc(x, self.u1)])
 
         N0 = self.L_scaled*self.c_scaled  # total amount of species in cell
         for k in range(self.M):
@@ -495,11 +490,11 @@ def useSternLayerCellBC(self, implicit=False):
                 N0[k], lwidth=self.label_width))
 
             self.boundary_conditions.extend([
-                lambda x, k=k: self.leftControlledVolumeSchemeFluxBC(x, k),
-                lambda x, k=k, N0=N0[k]: self.numberConservationConstraint(x, k, N0)])
+                lambda x, k=k: self.left_controlled_volume_scheme_flux_bc(x, k),
+                lambda x, k=k, N0=N0[k]: self.number_conservation_constraint(x, k, N0)])
 
     # TODO: meaningful test for Dirichlet BC
-    def useStandardDirichletBC(self):
+    def use_standard_dirichlet_bc(self):
         """Dirichlet BC for all variables at all boundaries"""
         self.boundary_conditions = []
 
@@ -513,22 +508,22 @@ def useStandardDirichletBC(self):
 
         # set up boundary conditions
         self.boundary_conditions.extend([
-          lambda x: self.leftPotentialDirichletBC(x, self.u0),
-          lambda x: self.rightPotentialDirichletBC(x, self.u1)])
+          lambda x: self.left_potential_dirichlet_bc(x, self.u0),
+          lambda x: self.right_potential_dirichlet_bc(x, self.u1)])
 
         for k in range(self.M):
-          self.logger.info(
-              'Ion species {:02d} left hand side concentration Dirichlet boundary condition:  c0 = {:> 8.4g}'.format(
-                  k, self.c_scaled[k]))
-          self.logger.info(
-              'Ion species {:02d} right hand side concentration Dirichlet boundary condition: c1 = {:> 8.4g}'.format(
-                  k, self.c_scaled[k]))
-          self.boundary_conditions.extend([
-              lambda x, k=k: self.leftDirichletBC(x, k, self.c_scaled[k]),
-              lambda x, k=k: self.rightDirichletBC(x, k, self.c_scaled[k])])
+            self.logger.info(
+                'Ion species {:02d} left hand side concentration Dirichlet boundary condition:  c0 = {:> 8.4g}'.format(
+                    k, self.c_scaled[k]))
+            self.logger.info(
+                'Ion species {:02d} right hand side concentration Dirichlet boundary condition: c1 = {:> 8.4g}'.format(
+                    k, self.c_scaled[k]))
+            self.boundary_conditions.extend([
+                lambda x, k=k: self.left_dirichlet_bc(x, k, self.c_scaled[k]),
+                lambda x, k=k: self.right_dirichlet_bc(x, k, self.c_scaled[k])])
 
     # boundary conditions and constraints building blocks:
-    def leftFiniteDifferenceSchemeFluxBC(self, x, k, j0=0):
+    def left_finite_difference_scheme_flux_bc(self, x, k, j0=0):
         """
         Parameters
         ----------
@@ -559,9 +554,9 @@ def leftFiniteDifferenceSchemeFluxBC(self, x, k, j0=0):
                 dndx, self.zi0[k, 0], nijk[0], dudx, self.dx, j0))
         return bcval
 
-    def rightFiniteDifferenceSchemeFluxBC(self, x, k, j0=0):
+    def right_finite_difference_scheme_flux_bc(self, x, k, j0=0):
         """
-        See ```leftFiniteDifferenceSchemeFluxBC```
+        See ```left_finite_difference_scheme_flux_bc```
         """
         uij = x[:self.Ni]
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
@@ -579,7 +574,7 @@ def rightFiniteDifferenceSchemeFluxBC(self, x, k, j0=0):
                 dndx, self.zi0[k, -1], nijk[-1], dudx, self.dx, j0))
         return bcval
 
-    def leftControlledVolumeSchemeFluxBC(self, x, k, j0=0):
+    def left_controlled_volume_scheme_flux_bc(self, x, k, j0=0):
         """
         Compute left hand side flux boundary condition residual in accord with
         controlled volume scheme.
@@ -609,10 +604,10 @@ def leftControlledVolumeSchemeFluxBC(self, x, k, j0=0):
             'CV flux BC F[0]  = n1*B(z(u0-u1)) - n0*B(z(u1-u0)) - j0*dx = {:> 8.4g}'.format(bcval))
         return bcval
 
-    def rightControlledVolumeSchemeFluxBC(self, x, k, j0=0):
+    def right_controlled_volume_scheme_flux_bc(self, x, k, j0=0):
         """
         Compute right hand side flux boundary condition residual in accord with
-        controlled volume scheme. See ``leftControlledVolumeSchemeFluxBC``
+        controlled volume scheme. See ``left_controlled_volume_scheme_flux_bc``
         """
         uij = x[:self.Ni]
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
@@ -626,25 +621,25 @@ def rightControlledVolumeSchemeFluxBC(self, x, k, j0=0):
             'CV flux BC F[-1]  = n[-1]*B(z(u[-2]-u[-1])) - n[-2]*B(z(u[-1]-u[-2])) - j0*dx = {:> 8.4g}'.format(bcval))
         return bcval
 
-    def leftPotentialDirichletBC(self, x, u0=0):
-        return self.leftDirichletBC(x, -1, u0)
+    def left_potential_dirichlet_bc(self, x, u0=0):
+        return self.left_dirichlet_bc(x, -1, u0)
 
-    def leftDirichletBC(self, x, k, x0=0):
+    def left_dirichlet_bc(self, x, k, x0=0):
         """Construct Dirichlet BC at left boundary"""
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
         return nijk[0] - x0
 
-    def rightPotentialDirichletBC(self, x, x0=0):
-        return self.rightDirichletBC(x, -1, x0)
+    def right_potential_dirichlet_bc(self, x, x0=0):
+        return self.right_dirichlet_bc(x, -1, x0)
 
-    def rightDirichletBC(self, x, k, x0=0):
+    def right_dirichlet_bc(self, x, k, x0=0):
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
         return nijk[-1] - x0
 
-    def leftPotentialRobinBC(self, x, lam, u0=0):
-        return self.leftRobinBC(x, -1, lam, u0)
+    def left_potential_robin_bc(self, x, lam, u0=0):
+        return self.left_robin_bc(x, -1, lam, u0)
 
-    def leftRobinBC(self, x, k, lam, x0=0):
+    def left_robin_bc(self, x, k, lam, x0=0):
         """
         Compute left hand side Robin (u + lam*dudx = u0 ) BC at in accord with
         2nd order finite difference scheme.
@@ -671,20 +666,20 @@ def leftRobinBC(self, x, k, lam, x0=0):
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
         return nijk[0] + lam/(2*self.dx) * (3.0*nijk[0] - 4.0*nijk[1] + nijk[2]) - x0
 
-    def rightPotentialRobinBC(self, x, lam, u0=0):
-        return self.rightRobinBC(x, -1, lam, u0)
+    def right_potential_robin_bc(self, x, lam, u0=0):
+        return self.right_robin_bc(x, -1, lam, u0)
 
-    def rightRobinBC(self, x, k, lam, x0=0):
+    def right_robin_bc(self, x, k, lam, x0=0):
         """Construct Robin (u + lam*dudx = u0 ) BC at right boundary."""
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
         return nijk[-1] + lam/(2*self.dx) * (3.0*nijk[-1] - 4.0*nijk[-2] + nijk[-3]) - x0
 
-    def numberConservationConstraint(self, x, k, N0):
+    def number_conservation_constraint(self, x, k, N0):
         """N0: total amount of species, k: ion species"""
         nijk = x[(k+1)*self.Ni:(k+2)*self.Ni]
 
-        ## TODO: this integration scheme assumes constant concentrations within
-        ## an interval. Adapt to controlled volume scheme!
+        # TODO: this integration scheme assumes constant concentrations within
+        # an interval. Adapt to controlled volume scheme!
 
         # rescale to fit interval
         N = np.sum(nijk*self.dx) * self.N / self.Ni
@@ -696,7 +691,7 @@ def numberConservationConstraint(self, x, k, N0):
         return constraint_val
 
     # TODO: remove or standardize
-    # def leftNeumannBC(self,x,j0):
+    # def left_neumann_bc(self,x,j0):
     #   """Construct finite difference Neumann BC (flux BC) at left boundary"""
     #   # right hand side first derivative of second order error
     #   # df0dx = 1 / (2*dx) * (-3 f0 + 4 f1 - f2 ) + O(dx^2) = j0
@@ -705,7 +700,7 @@ def numberConservationConstraint(self, x, k, N0):
     #     'Neumann BC F[0]  = -3*x[0]  + 4*x[1]  - x[2]  = {:> 8.4g}'.format(bcval))
     #   return bcval
     #
-    # def rightNeumannBC(self,x,j0):
+    # def right_neumann_bc(self,x,j0):
     #   """Construct finite difference Neumann BC (flux BC) at right boundary"""
     #   # left hand side first derivative of second order error
     #   # dfndx = 1 / (2*dx) * (+3 fn - 4 fn-1 + fn-2 ) + O(dx^2) = 0
@@ -798,7 +793,7 @@ def nernst_planck_pde(self, x):
         for k in range(self.M):
             self.logger.debug(
                 'ion species {:02d} concentration range [c_min, c_max] = [ {:>.4g}, {:>.4g} ]'.format(
-                k, np.min(nijk1[k, :]), np.max(nijk1[k, :])))
+                    k, np.min(nijk1[k, :]), np.max(nijk1[k, :])))
 
         Fn = np.zeros([self.M, self.Ni])
         # loop over k = 1..M reduced Nernst-Planck equations:
@@ -858,12 +853,12 @@ def G(self, x):
         return F
 
     @property
-    def I(self):  # ionic strength
+    def ionic_strength(self):  # ionic strength
         """Compute the system's ionic strength from charges and concentrations.
 
         Returns
         -------
-        I : float
+        ionic_strength : float
             ionic strength ( 1/2 * sum(z_i^2*c_i) )
             [concentration unit, i.e. mol m^-3]
         """
@@ -880,7 +875,7 @@ def lambda_D(self):
         """
         return np.sqrt(
             self.relative_permittivity*self.vacuum_permittivity*self.R*self.T/(
-                2.0*self.F**2*self.I))
+                2.0*self.F**2*self.ionic_strength))
 
     # default 0.1 mM (i.e. mol/m^3) NaCl aqueous solution
     def init(self,
@@ -938,7 +933,7 @@ def init(self,
             absolute tolerance for Newton solver convergence (default: 1e-10)
         maxit : int, optional
             maximum number of Newton iterations (default: 20)
-        solver: func( funx(x), x0), optional
+        solver: func( func(x), x0), optional
             solver to use (default: None, will use own simple Newton solver)
         potential0: (N+1,) ndarray, optional (default: None)
             potential initial values
@@ -981,7 +976,6 @@ def init(self,
         self.x0 = x0  # reference position
         self.delta_u = delta_u  # potential difference
 
-
         self.relative_permittivity = relative_permittivity
         self.vacuum_permittivity = vacuum_permittivity
         # R = N_A * k_B
@@ -1025,7 +1019,7 @@ def init(self,
         self.l_unit = self.lambda_D
 
         # concentration unit is ionic strength
-        self.c_unit = self.I
+        self.c_unit = self.ionic_strength
 
         # no time unit for now, only steady state
         # self.t_unit = self.l_unit**2 / self.Dn # fixes Dn_scaled = 1
diff --git a/matscipy/electrochemistry/poisson_nernst_planck_solver_fenics.py b/matscipy/electrochemistry/poisson_nernst_planck_solver_fenics.py
index 064e68dc..34a9054b 100644
--- a/matscipy/electrochemistry/poisson_nernst_planck_solver_fenics.py
+++ b/matscipy/electrochemistry/poisson_nernst_planck_solver_fenics.py
@@ -35,7 +35,7 @@
 
 class Boundary(fn.SubDomain):
     """Mark a point to be within the domain boundary for fenics."""
-    # Boundary causes crash kernel if __init__ doe not call super().__init__()
+    # Boundary causes crash kernel if __init__ does not call super().__init__()
     def __init__(self, x0=0, tol=1e-14):
         super().__init__()
         self.tol = tol
@@ -119,44 +119,44 @@ def boundary_C(self, x, on_boundary):
         """Mark domain center."""
         return fn.near(x[0], (self.x0_scaled + self.x1_scaled)/2., self.bctol)
 
-    def applyLeftPotentialDirichletBC(self, u0):
+    def apply_left_potential_dirichlet_bc(self, u0):
         """FEniCS Dirichlet BC u0 for potential at left boundary."""
         self.boundary_conditions.extend([
             fn.DirichletBC(self.W.sub(0), u0,
                            lambda x, on_boundary: self.boundary_L(x, on_boundary))])
 
-    def applyRightPotentialDirichletBC(self, u0):
+    def apply_right_potential_dirichlet_bc(self, u0):
         """FEniCS Dirichlet BC u0 for potential at right boundary."""
         self.boundary_conditions.extend([
             fn.DirichletBC(self.W.sub(0), u0,
                            lambda x, on_boundary: self.boundary_R(x, on_boundary))])
 
-    def applyLeftConcentrationDirichletBC(self, k, c0):
+    def apply_left_concentration_dirichlet_bc(self, k, c0):
         """FEniCS Dirichlet BC c0 for k'th ion species at left boundary."""
         self.boundary_conditions.extend([
             fn.DirichletBC(self.W.sub(k+1), c0,
                            lambda x, on_boundary: self.boundary_L(x, on_boundary))])
 
-    def applyRightConcentrationDirichletBC(self, k, c0):
+    def apply_right_concentration_dirichlet_bc(self, k, c0):
         """FEniCS Dirichlet BC c0 for k'th ion species at right boundary."""
         self.boundary_conditions.extend([
             fn.DirichletBC(self.W.sub(k+1), c0,
                            lambda x, on_boundary: self.boundary_R(x, on_boundary))])
 
-    def applyCentralReferenceConcentrationConstraint(self, k, c0):
+    def apply_central_reference_concentration_constraint(self, k, c0):
         """FEniCS Dirichlet BC c0 for k'th ion species at right boundary."""
         self.boundary_conditions.extend([
             fn.DirichletBC(self.W.sub(k+1), c0,
                            lambda x, on_boundary: self.boundary_C(x, on_boundary))])
 
     # TODO: Robin BC!
-    def applyLeftPotentialRobinBC(self, u0, lam0):
+    def apply_left_potential_robin_bc(self, u0, lam0):
         self.logger.warning("Not implemented!")
 
-    def applyRightPotentialRobinBC(self, u0, lam0):
+    def apply_right_potential_robin_bc(self, u0, lam0):
         self.logger.warning("Not implemented!")
 
-    def applyNumberConservationConstraint(self, k, c0):
+    def apply_number_conservation_constraint(self, k, c0):
         """
         Enforce number conservation constraint via Lagrange multiplier.
         See https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/neumann-poisson/python/documentation.html
@@ -164,17 +164,17 @@ def applyNumberConservationConstraint(self, k, c0):
         self.constraints += self.lam[k]*self.q[k]*fn.dx \
             + (self.p[k]-c0)*self.mu[k]*fn.dx
 
-    def applyPotentialDirichletBC(self, u0, u1):
+    def apply_potential_dirichlet_bc(self, u0, u1):
         """Potential Dirichlet BC u0 and u1 on left and right boundary."""
-        self.applyLeftPotentialDirichletBC(u0)
-        self.applyRightPotentialDirichletBC(u1)
+        self.apply_left_potential_dirichlet_bc(u0)
+        self.apply_right_potential_dirichlet_bc(u1)
 
-    def applyPotentialRobinBC(self, u0, u1, lam0, lam1):
+    def apply_potential_robin_bc(self, u0, u1, lam0, lam1):
         """Potential Robin BC on left and right boundary."""
-        self.applyLeftPotentialRobinBC(u0, lam0)
-        self.applyRightPotentialRobinBC(u1, lam1)
+        self.apply_left_potential_robin_bc(u0, lam0)
+        self.apply_right_potential_robin_bc(u1, lam1)
 
-    def useStandardInterfaceBC(self):
+    def use_standard_interface_bc(self):
         """Interface at left hand side and open bulk at right hand side"""
         self.boundary_conditions = []
 
@@ -185,14 +185,14 @@ def useStandardInterfaceBC(self):
         self.logger.info('Left hand side Dirichlet boundary condition:  u0 = {:> 8.4g}'.format(self.u0))
         self.logger.info('Right hand side Dirichlet boundary condition: u1 = {:> 8.4g}'.format(self.u1))
 
-        self.applyPotentialDirichletBC(self.u0, self.u1)
+        self.apply_potential_dirichlet_bc(self.u0, self.u1)
 
         for k in range(self.M):
             self.logger.info(('Ion species {:02d} right hand side concentration '
                               'Dirichlet boundary condition: c1 = {:> 8.4g}').format(k, self.c_scaled[k]))
-            self.applyRightConcentrationDirichletBC(k, self.c_scaled[k])
+            self.apply_right_concentration_dirichlet_bc(k, self.c_scaled[k])
 
-    def useStandardCellBC(self):
+    def use_standard_cell_bc(self):
         """
         Interfaces at left hand side and right hand side, species-wise
         number conservation within interval."""
@@ -211,7 +211,7 @@ def useStandardCellBC(self):
         self.logger.info('{:>{lwidth}s} u1 = {:< 8.4g}'.format(
           'Right hand side Dirichlet boundary condition', self.u1, lwidth=self.label_width))
 
-        self.applyPotentialDirichletBC(self.u0, self.u1)
+        self.apply_potential_dirichlet_bc(self.u0, self.u1)
 
         # Number conservation constraints
         self.constraints = 0
@@ -220,9 +220,9 @@ def useStandardCellBC(self):
             self.logger.info('{:>{lwidth}s} N0 = {:<8.4g}'.format(
                 'Ion species {:02d} number conservation constraint'.format(k),
                 N0[k], lwidth=self.label_width))
-            self.applyNumberConservationConstraint(k, self.c_scaled[k])
+            self.apply_number_conservation_constraint(k, self.c_scaled[k])
 
-    def useCentralReferenceConcentrationBasedCellBC(self):
+    def use_central_reference_concentration_based_cell_bc(self):
         """
         Interfaces at left hand side and right hand side, species-wise
         concentration fixed at cell center."""
@@ -241,15 +241,15 @@ def useCentralReferenceConcentrationBasedCellBC(self):
         self.logger.info('{:>{lwidth}s} u1 = {:< 8.4g}'.format(
           'Right hand side Dirichlet boundary condition', self.u1, lwidth=self.label_width))
 
-        self.applyPotentialDirichletBC(self.u0, self.u1)
+        self.apply_potential_dirichlet_bc(self.u0, self.u1)
 
         for k in range(self.M):
             self.logger.info(
                 'Ion species {:02d} reference concentration condition: c1 = {:> 8.4g} at cell center'.format(
                     k, self.c_scaled[k]))
-            self.applyCentralReferenceConcentrationConstraint(k, self.c_scaled[k])
+            self.apply_central_reference_concentration_constraint(k, self.c_scaled[k])
 
-    def useSternLayerCellBC(self):
+    def use_stern_layer_cell_bc(self):
         """
         Interfaces at left hand side and right hand side, species-wise
         number conservation within interval."""
@@ -298,7 +298,7 @@ def useSternLayerCellBC(self):
             self.logger.info('{:>{lwidth}s} N0 = {:<8.4g}'.format(
                 'Ion species {:02d} number conservation constraint'.format(k),
                 N0[k], lwidth=self.label_width))
-            self.applyNumberConservationConstraint(k, self.c_scaled[k])
+            self.apply_number_conservation_constraint(k, self.c_scaled[k])
 
     def discretize(self):
         """Builds function space, call again after introducing constraints"""
diff --git a/matscipy/electrochemistry/steric_correction.py b/matscipy/electrochemistry/steric_correction.py
index 1a8d5d64..72bcee6e 100644
--- a/matscipy/electrochemistry/steric_correction.py
+++ b/matscipy/electrochemistry/steric_correction.py
@@ -642,7 +642,7 @@ def neigh_list_based_target_function(x, r=1.0, constraints=None, Dij=None):
     r: float or (N,) ndarray, optional (default=1.0)
         steric radii of particles
     constraints: callable, returns (float, (N,dim) ndarray)
-        constraint function value and gradien
+        constraint function value and gradient
     Dij: (N, N) ndarray, optional (default=None)
         pairwise minimum allowed distance matrix, overrides r
 
@@ -870,6 +870,8 @@ def apply_steric_correction(
         Forwarded to scipy minimzer.
         https://docs.scipy.org/doc/scipy/reference/optimize.minimize-bfgs.html
         (default: {'gtol':1.e-5,'maxiter':10,'disp':True,'eps':1.e-8})
+    method: str
+        scipy.optimize.minimize method
     target_function: func, optional
         One of the target functions within this submodule, or function
         of same signature. (default: neigh_list_based_target_function)
diff --git a/tests/test_electrochemistry_cli.py b/tests/test_electrochemistry_cli.py
index dd028da6..b7453e04 100644
--- a/tests/test_electrochemistry_cli.py
+++ b/tests/test_electrochemistry_cli.py
@@ -182,15 +182,29 @@ def test_pnp_c2d_pipeline_mode(self):
                 stdout=subprocess.PIPE, stderr=subprocess.PIPE,
                 cwd=tmpdir, encoding='utf-8', env=self.myenv)
 
-            c2d = subprocess.Popen([ 'c2d' ],
-                stdin=pnp.stdout, stdout=subprocess.PIPE, stderr=subprocess.PIPE,
+            # detour via python string
+            pnp_output, pnp_error = pnp.communicate()
+
+            c2d = subprocess.Popen(['c2d'],
+                stdin=subprocess.PIPE, stdout=subprocess.PIPE, stderr=subprocess.PIPE,
                 cwd=tmpdir, encoding='utf-8', env=self.myenv)
-            pnp.stdout.close()  # Allow pnp to receive a SIGPIPE if p2 exits.
+
+            # pnp.stdout.close()  # Allow pnp to receive a SIGPIPE if p2 exits.
+
+            c2d_output, c2d_error = c2d.communicate(input=pnp_output)
+
+            print("  poisson-nernst-planck stdout")
+            print(pnp_output)
+            print("  poisson-nernst-planck stderr")
+            print(pnp_error)
+            print("  continuous2discrete output")
+            print(c2d_output)
+            print("  continuous2discrete stderr")
+            print(c2d_error)
 
             self.assertEqual(c2d.wait(),0)
             self.assertEqual(pnp.wait(),0)
 
-            c2d_output = c2d.communicate()[0]
             with io.StringIO(c2d_output) as xyz_instream:
                 xyz = ase.io.read(xyz_instream, format='extxyz')
             self.assertEqual(len(xyz),len(self.ref_xyz))
diff --git a/tests/test_poisson_nernst_planck_solver.py b/tests/test_poisson_nernst_planck_solver.py
index 1f898e60..d622f523 100644
--- a/tests/test_poisson_nernst_planck_solver.py
+++ b/tests/test_poisson_nernst_planck_solver.py
@@ -39,7 +39,7 @@ def test_poisson_nernst_planck_solver_std_interface_bc(self):
         pnp = PoissonNernstPlanckSystem(
             c=[0.1,0.1], z=[1,-1], L=1e-7, delta_u=0.05,
             N=200, e=1e-12, maxit=20)
-        pnp.useStandardInterfaceBC()
+        pnp.use_standard_interface_bc()
         pnp.solve()
 
         self.assertArrayAlmostEqual(pnp.grid, self.ref_data ['x'])
diff --git a/tests/test_poisson_nernst_planck_solver_fenics.py b/tests/test_poisson_nernst_planck_solver_fenics.py
index 0279d08d..d9be8b19 100644
--- a/tests/test_poisson_nernst_planck_solver_fenics.py
+++ b/tests/test_poisson_nernst_planck_solver_fenics.py
@@ -26,12 +26,11 @@
 
 try:
     import fenics
+    from matscipy.electrochemistry.poisson_nernst_planck_solver_fenics \
+        import PoissonNernstPlanckSystemFEniCS as PoissonNernstPlanckSystem
 except ImportError:
     print("fenics not found: skipping fenics-dependent tests")
 
-from matscipy.electrochemistry.poisson_nernst_planck_solver_fenics \
-    import PoissonNernstPlanckSystemFEniCS as PoissonNernstPlanckSystem
-
 
 class PoissonNernstPlanckSolverTest(matscipytest.MatSciPyTestCase):
 
@@ -49,7 +48,7 @@ def test_poisson_nernst_planck_solver_fenics_std_interface_bc(self):
         pnp = PoissonNernstPlanckSystem(
             c=[0.1,0.1], z=[1,-1], L=1e-7, delta_u=0.05,
             N=200, e=1e-12, maxit=20)
-        pnp.useStandardInterfaceBC()
+        pnp.use_standard_interface_bc()
         pnp.solve()
 
         # Reference data has been generated with controlled-volume solver and is slightly off the FEM results,