diff --git a/README.md b/README.md
index cb9621a..c361544 100644
--- a/README.md
+++ b/README.md
@@ -25,8 +25,8 @@ The easiest way is to use the environment file, compatible with `conda`,
 `mamba`, `micromamba` etc.
 
 ```bash
-mamba env create -f environment.yml
-mamba activate rgpotpy
+micromamba env create -f environment.yml
+micromamba activate rgpotpy
 pdm install
 ```
 
diff --git a/meson.build b/meson.build
index 5fbe894..c6530f1 100644
--- a/meson.build
+++ b/meson.build
@@ -61,7 +61,7 @@ py.install_sources([
 # Adapters
 py.install_sources([
     'pypotlib/ase_adapters.py',
-    'pypotlib/aux.py',
+    'pypotlib/_aux.py',
   ],
   pure: false,
   subdir: 'pypotlib'
diff --git a/pypotlib/aux.py b/pypotlib/_aux.py
similarity index 100%
rename from pypotlib/aux.py
rename to pypotlib/_aux.py
diff --git a/pypotlib/systems/cu_slab_h2.py b/pypotlib/systems/cu_slab_h2.py
index 792d006..352ddb3 100644
--- a/pypotlib/systems/cu_slab_h2.py
+++ b/pypotlib/systems/cu_slab_h2.py
@@ -1,6 +1,13 @@
-__all__ = ["CuH2PotSlab", "FuncVal", "PltRange",
-           "PlotPtPos", "contour_plot", "plt_data",
-           "set_slab_dist"]
+__all__ = [
+    "CuH2PotSlab",
+    "FuncVal",
+    "PltRange",
+    "PlotPtPos",
+    "contour_plot",
+    "plt_data",
+    "set_slab_dist",
+    "calculate_hh_and_hcu_distances",
+]
 
 import warnings
 from collections import namedtuple
@@ -8,23 +15,22 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from ase.calculators.calculator import Calculator, all_changes
+import cmcrameri.cm as cmc
 
 from pypotlib import cpot
-from pypotlib.aux import reshape_data
+from pypotlib._aux import reshape_data
 
 FuncVal = namedtuple("FuncVal", ["x", "y", "energy"])
 PltRange = namedtuple("PltRange", ["low", "high"])
-PlotPtPosData = namedtuple('PlotPtPosData', ['h2pos', 'pltpts'])
-PlotPoints = namedtuple('PlotPoints', ['x_npt', 'y_npt'])
+PlotPtPosData = namedtuple("PlotPtPosData", ["h2pos", "pltpts"])
+PlotPoints = namedtuple("PlotPoints", ["x_npt", "y_npt"])
+
 
 class CuH2PotSlab(Calculator):
     implemented_properties = ["energy", "forces"]
     discard_results_on_any_change = True
 
-    def __init__(self,
-                 e_zero=FuncVal(x=0.752542,
-                                y=5.0,
-                                energy=-697.311695)):
+    def __init__(self, e_zero=FuncVal(x=0.752542, y=5.0, energy=-697.311695)):
         Calculator.__init__(self)
         self.cpot = cpot.CuH2Pot()
         self._ezero = e_zero
@@ -43,6 +49,7 @@ def calculate(
         self.results["energy"] = energy - self._ezero.energy
         self.results["forces"] = forces
 
+
 def set_slab_dist(atoms, dist):
     h_ind = np.where(np.asarray(atoms.get_chemical_symbols()) == "H")[0]
     n_hatoms = len(h_ind)
@@ -51,15 +58,50 @@ def set_slab_dist(atoms, dist):
     cpos[-n_hatoms:, 2] = np.repeat(max_slab_z + dist, n_hatoms)
     atoms.set_positions(cpos)
 
-def plt_data(_atms, hh_range=PltRange(low=0.4, high=3),
-                  h2slab_range = PltRange(low=-0.05, high=5),
-                  n_points=PlotPoints(x_npt=60, y_npt=60)):
-    h_dists = np.linspace(
-        hh_range.low, hh_range.high, n_points.x_npt
-    )
+
+def calculate_hh_and_hcu_distances(h_positions, ref_atoms):
+    """
+    Calculate the H-H distance and the specific H(midpoint)-slab distance based on the reference configuration.
+
+    Parameters:
+    - h_positions: A 6-element list or array with the coordinates of the H atoms [x1, y1, z1, x2, y2, z2].
+    - ref_atoms: An ASE Atoms object containing the reference configuration, including H and Cu atoms.
+
+    Returns:
+    - hh_distance: Distance between the two H atoms.
+    - hcu_distance: Distance from the H-H midpoint to the slab, considering the reference configuration.
+    """
+    # Reshape h_positions to a 2x3 matrix for easier manipulation
+    h_positions_matrix = np.array(h_positions).reshape(2, 3)
+
+    # Calculate H-H midpoint
+    midpoint = np.mean(h_positions_matrix, axis=0)
+
+    # Calculate H-H distance
+    hh_distance = np.linalg.norm(h_positions_matrix[0] - h_positions_matrix[1])
+
+    # Determine the maximum Z position of Cu atoms in the reference configuration as the reference slab surface
+    cu_positions = ref_atoms.positions[ref_atoms.numbers == 29]
+    max_cu_z = np.max(cu_positions[:, 2])
+
+    # Calculate the distance from the H-H midpoint to the reference slab surface
+    hcu_distance = midpoint[2] - max_cu_z
+
+    return hh_distance, hcu_distance
+
+
+def plt_data(
+    atms,
+    hh_range=PltRange(low=0.4, high=3.8),
+    h2slab_range=PltRange(low=-0.05, high=5),
+    n_points=PlotPoints(x_npt=60, y_npt=60),
+):
+    h_dists = np.linspace(hh_range.low, hh_range.high, n_points.x_npt)
     slab_dists = np.linspace(
         h2slab_range.low, h2slab_range.high, n_points.y_npt
     )
+    _atms = atms.copy()
+    _atms.calc = atms.calc
     h_ind = np.where(np.asarray(_atms.get_chemical_symbols()) == "H")[0]
     plt_data = []  # x, y, energy
     plt_h2pos = []  # h2 positions
@@ -77,7 +119,8 @@ def plt_data(_atms, hh_range=PltRange(low=0.4, high=3),
             )
     return PlotPtPosData(h2pos=plt_h2pos, pltpts=plt_data)
 
-def contour_plot(data, _max_val = 5, _nlvls=500):
+
+def contour_plot(data, _max_val=5, _nlvls=500, scatter_points=None, title=None):
     x, y, energy = reshape_data(data)
     # Create a contour plot of the energy surface
     fig, ax = plt.subplots()
@@ -87,12 +130,28 @@ def contour_plot(data, _max_val = 5, _nlvls=500):
         y,
         energy,
         levels=np.linspace(0, _max_val, _nlvls),
-        extend="max",
+        extend="both",
+        cmap=cmc.batlow,
     )
+    # Scatter
+    if scatter_points is not None:
+        scatter_points = np.array(scatter_points)
+        ax.scatter(
+            scatter_points[:, 0],
+            scatter_points[:, 1],
+            marker="+",
+            color="lightgray",
+        )
+
     # Add labels and title to the plot
     ax.set_xlabel("H-H distance")
     ax.set_ylabel("Slab distance")
-    ax.set_title('Energy Surface Contour Plot\nShifted by the zero of the energy')
+    if title:
+        ax.set_title(title)
+    else:
+        ax.set_title(
+            "Energy Surface Contour Plot\nShifted by the zero of the energy"
+        )
     # Add a colorbar to the plot
     cbar = fig.colorbar(cs)
     # Show the plot
diff --git a/pyproject.toml b/pyproject.toml
index 61727a4..d1b39a6 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -9,10 +9,12 @@ authors = [
 ]
 # These are only if ase integration is requested, consider a group
 dependencies = [
-    "numpy>=1.24.3",
+    "numpy",
     "ase>=3.22.1",
+    "cmcrameri>=1.8",
+    "matplotlib>=3.9.0.dev0",
 ]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 readme = "README.md"
 license = {text = "MIT"}
 
@@ -49,3 +51,10 @@ all = {composite = ["lint pypotlib/"]}
 [tool.black]
 line-length = 80
 target-version = ['py310']
+
+[[tool.pdm.source]]
+name = "nightly-scientific-python"
+url = "https://pypi.anaconda.org/scientific-python-nightly-wheels/simple/"
+verify_ssl = true
+include_packages = ["scipy", "scipy-*", "numpy", "numpy-*", "matplotlib", "matplotlib-*"]
+exclude_packages = ["*"]