diff --git a/pyproject.toml b/pyproject.toml
index 63ff2ad93..87880e9a1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -11,7 +11,7 @@
 [build-system]
 requires = [
     "setuptools_scm[toml]>=6.2",
-    "scikit-build-core[pyproject]",
+    "scikit-build-core[pyproject]==0.5.1",
     "Cython>=0.29.28",
     # see https://github.com/scipy/oldest-supported-numpy/
     # should fix #310
diff --git a/src/sisl/geom/flat.py b/src/sisl/geom/flat.py
index 3977f899f..3e55fdddd 100644
--- a/src/sisl/geom/flat.py
+++ b/src/sisl/geom/flat.py
@@ -74,13 +74,39 @@ def graphene(bond=1.42, atoms=None, orthogonal=False):
         atoms = Atom(Z=6, R=bond * 1.01)
     return honeycomb(bond, atoms, orthogonal)
 
-def honeycomb_flake(shells: int, bond: float, atoms) -> Geometry:
+def honeycomb_flake(shells: int, bond: float, atoms, vaccuum: float = 20.) -> Geometry:
+    """Hexagonal flake of a honeycomb lattice, with zig-zag edges.
 
+    Parameters
+    ----------
+    shells:
+        Number of shells in the flake. 0 means a single hexagon, and subsequent
+        shells add hexagon units surrounding the previous shell.
+    bond:
+        bond length between atoms (*not* lattice constant)
+    atoms:
+        the atom (or atoms) that the honeycomb lattice consists of
+    vaccuum:
+        Amount of vacuum to add to the cell on all directions
+    """
+
+    # Function that generates one of the six triangular portions of the
+    # hexagonal flake. The rest of the portions are obtained by rotating
+    # this one by 60, 120, 180, 240 and 300 degrees.
     def _minimal_op(shells):
 
+        # The function is based on the horizontal lines of the hexagon,
+        # which are made of a pair of atoms.
+        # For each shell, we first need to complete the incomplete horizontal
+        # lines of the previous shell, and then branch them up and down to create
+        # the next horizontal lines.
+
+        # Displacement from the end of one horizontal pair to the beggining of the next
         branch_displ_x = bond * np.cos(np.radians(60))
         branch_displ_y = bond * np.sin(np.radians(60))
 
+        # Iterate over shells. We also keep track of the atom types, in case
+        # we have two different atoms in the honeycomb lattice.
         op = np.array([[bond, 0, 0]])
         types = np.array([0])
         for shell in range(shells):
@@ -102,24 +128,52 @@ def _minimal_op(shells):
 
         return op, types
     
+    # Get the coordinates of 1/6 of the hexagon for the requested number of shells.
     op, types = _minimal_op(shells)
 
-    # Get the first two of the six portions of the hexagon. If there are two 
-    # atom types, then we need to reverse the order of the atoms of the second portion.
-    single_atom_type = isinstance(str, atoms) or len(atoms) == 1
+    single_atom_type = isinstance(atoms, (str, Atom)) or len(atoms) == 1
 
+    # Create a geometry from the coordinates.
     ats = atoms  if single_atom_type else np.asarray(atoms)[types]
     geom = Geometry(op, atoms=ats)
 
+    # The second portion of the hexagon is obtained by rotating the first one by 60 degrees.
+    # However, if there are two different atoms in the honeycomb lattice, we need to reverse the types.
     next_triangle = geom if single_atom_type else Geometry(op, atoms=np.asarray(atoms)[types - 1])
     geom += next_triangle.rotate(60, [0,0,1])
 
-    # Then just rotate the two triangles by 120 and 240 degrees to get the full hexagon
+    # Then just rotate the two triangles by 120 and 240 degrees to get the full hexagon.
     geom += geom.rotate(120, [0,0,1]) + geom.rotate(240, [0,0,1])
 
-    return geom
+    # Set the cell according to the requested vacuum
+    max_x = np.max(geom.xyz[:,0])
+    geom.cell[0,0] = max_x * 2 + vaccuum
+    geom.cell[1,1] = max_x * 2 + vaccuum
+    geom.cell[2,2] = 20.
+
+    # Center the flake and return
+    return geom.translate(geom.center(what="cell"))
+
+def graphene_flake(shells: int, bond: float = 1.42, atoms=None, vaccuum: float = 20.) -> Geometry:
+    """Hexagonal flake of graphene, with zig-zag edges.
 
-def graphene_flake(shells: int, bond: float = 1.42, atoms=None) -> Geometry:
+    Parameters
+    ----------
+    shells:
+        Number of shells in the flake. 0 means a single hexagon, and subsequent
+        shells add hexagon units surrounding the previous shell.
+    bond:
+        bond length between atoms (*not* lattice constant)
+    atoms:
+        the atom (or atoms) that the honeycomb lattice consists of.
+        Default to Carbon atom.
+    vaccuum:
+        Amount of vacuum to add to the cell on all directions
+
+    See Also
+    --------
+    honeycomb_flake: the equivalent of this, but with non-default atoms.
+    """
     if atoms is None:
         atoms = Atom(Z=6, R=bond * 1.01)
-    return honeycomb_flake(shells, bond, atoms)
+    return honeycomb_flake(shells, bond, atoms, vaccuum)