diff --git a/CHANGELOG.md b/CHANGELOG.md
index 136fac2e4e..003fe67383 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -8,6 +8,8 @@ we hit release version 1.0.0.
 ## [0.14.3] - YYYY-MM-DD
 
 ### Added
+- Creation of honeycomb flakes (`sisl.geom.honeycomb_flake`,
+`sisl.geom.graphene_flake`). #636
 - added `Geometry.as_supercell` to create the supercell structure,
   thanks to @pfebrer for the suggestion
 - added `Lattice.to` and `Lattice.new` to function the same
diff --git a/docs/api/default_geom.rst b/docs/api/default_geom.rst
index 466bb36351..c2436278e0 100644
--- a/docs/api/default_geom.rst
+++ b/docs/api/default_geom.rst
@@ -67,6 +67,8 @@ Surfaces (slabs)
    honeycomb
    bilayer
    graphene
+   honeycomb_flake
+   graphene_flake
 
 
 Helpers
diff --git a/src/sisl/geom/flat.py b/src/sisl/geom/flat.py
index 42e39e2351..36f2fe1adc 100644
--- a/src/sisl/geom/flat.py
+++ b/src/sisl/geom/flat.py
@@ -8,7 +8,7 @@
 
 from ._common import geometry_define_nsc
 
-__all__ = ['honeycomb', 'graphene']
+__all__ = ['honeycomb', 'graphene', 'honeycomb_flake', 'graphene_flake']
 
 
 @set_module("sisl.geom")
@@ -73,3 +73,114 @@ def graphene(bond=1.42, atoms=None, orthogonal=False):
     if atoms is None:
         atoms = Atom(Z=6, R=bond * 1.01)
     return honeycomb(bond, atoms, orthogonal)
+
+@set_module("sisl.geom")
+def honeycomb_flake(shells: int, bond: float, atoms, vacuum: 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
+    vacuum:
+        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 * 0.5 # cos(60) = 0.5
+        branch_displ_y = bond * 3 ** 0.5 / 2 # sin(60) = sqrt(3)/2
+
+        # 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):
+            n_new_branches = 2 + shell
+            prev_edge = branch_displ_y * (shell)
+            
+            sat = np.zeros((shell + 1, 3))
+            sat[:, 0] = op[-1, 0] + bond
+            sat[:, 1] = np.linspace(-prev_edge, prev_edge, shell + 1)
+            
+            edge = branch_displ_y * (shell + 1)
+
+            branches = np.zeros((n_new_branches, 3))
+            branches[:, 0] = sat[0, 0] + branch_displ_x
+            branches[:, 1] = np.linspace(-edge, edge, n_new_branches )
+
+            op = np.concatenate([op, sat, branches])
+            types = np.concatenate([types, np.full(len(sat), 1), np.full(len(branches), 0)])
+
+        return op, types
+    
+    # Get the coordinates of 1/6 of the hexagon for the requested number of shells.
+    op, types = _minimal_op(shells)
+
+    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.
+    geom += geom.rotate(120, [0,0,1]) + geom.rotate(240, [0,0,1])
+
+    # Set the cell according to the requested vacuum
+    max_x = np.max(geom.xyz[:,0])
+    geom.cell[0,0] = max_x * 2 + vacuum
+    geom.cell[1,1] = max_x * 2 + vacuum
+    geom.cell[2,2] = 20.
+
+    # Center the flake
+    geom = geom.translate(geom.center(what="cell"))
+
+    # Set boundary conditions
+    geometry_define_nsc(geom, [False, False, False])
+
+    return geom
+
+@set_module("sisl.geom")
+def graphene_flake(shells: int, bond: float = 1.42, atoms=None, vacuum: float = 20.) -> Geometry:
+    """Hexagonal flake of graphene, 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.
+        Default to Carbon atom.
+    vacuum:
+        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, vacuum)
diff --git a/src/sisl/geom/tests/test_geom.py b/src/sisl/geom/tests/test_geom.py
index 0ba7806829..ac97c5eef0 100644
--- a/src/sisl/geom/tests/test_geom.py
+++ b/src/sisl/geom/tests/test_geom.py
@@ -51,6 +51,26 @@ def test_flat():
     a = graphene(orthogonal=True)
     assert is_right_handed(a)
 
+def test_flat_flakes():
+    g = graphene_flake(shells=0, bond=1.42)
+    assert g.na == 6
+    # All atoms are close to the center
+    assert len(g.close(g.center(), 1.44)) == g.na
+    # All atoms have two neighbours
+    assert len(g.axyz(AtomNeighbours(min=2, max=2, R=1.44))) == g.na
+
+    g = graphene_flake(shells=1, bond=1.42)
+    assert g.na == 24
+    assert len(g.close(g.center(), 4)) == g.na
+    assert len(g.axyz(AtomNeighbours(min=2, max=2, R=1.44))) == 12
+    assert len(g.axyz(AtomNeighbours(min=3, max=3, R=1.44))) == 12
+
+    bn = honeycomb_flake(shells=1, atoms=['B', 'N'], bond=1.42)
+    assert bn.na == 24
+    assert np.allclose(bn.xyz, g.xyz)
+    # Check that atoms are alternated.
+    assert len(bn.axyz(AtomZ(5) & AtomNeighbours(min=1, R=1.44, neighbour=AtomZ(5)))) == 0
+
 def test_nanotube():
     a = nanotube(1.42)
     assert is_right_handed(a)