def get_boundary_polyhedra(polyhedra, boundary_x=0, boundary_width=0.5, verbose=True, z_lim=[0, 100]):
    """
    get indices of polyhedra at boundary (assumed to be parallel to x-axis)

    Parameter
    ---------
    polyhedra: dict
        dictionary of all polyhedra
    boundary_x: float
        position of boundary in Angstrom
    boundary_width: float
        width of boundary where center of polyhedra are considered in Angstrom
    verbose: boolean
        optional
    z_lim: list
        upper and lower limit of polyhedra to plot

    Returns
    -------
    boundary_polyhedra: list
        list of polyhedra at boundary
    """
    boundary_polyhedra = []
    for key, polyhedron in polyhedra.items():
        center = polyhedron['vertices'].mean(axis=0)
        if abs(center[0] - boundary_x) < 0.5 and (z_lim[0] < center[2] < z_lim[1]):
            boundary_polyhedra.append(key)
            if verbose:
                print(key, polyhedron['length'], center)

    return boundary_polyhedra