Returns all symmetrically distinct indices below a certain max-index for
    a given structure. Analysis is based on the symmetry of the reciprocal
    lattice of the structure.
    Args:
        structure (Structure): input structure.
        max_index (int): The maximum index. For example, a max_index of 1
            means that (100), (110), and (111) are returned for the cubic
            structure. All other indices are equivalent to one of these.

def get_symmetrically_distinct_miller_indices(structure, max_index):
    """
    Returns all symmetrically distinct indices below a certain max-index for
    a given structure. Analysis is based on the symmetry of the reciprocal
    lattice of the structure.
    Args:
        structure (Structure): input structure.
        max_index (int): The maximum index. For example, a max_index of 1
            means that (100), (110), and (111) are returned for the cubic
            structure. All other indices are equivalent to one of these.
    """

    r = list(range(-max_index, max_index + 1))
    r.reverse()

    # First we get a list of all hkls for conventional (including equivalent)
    conv_hkl_list = [miller for miller in itertools.product(r, r, r) if any([i != 0 for i in miller])]

    sg = SpacegroupAnalyzer(structure)
    # Get distinct hkl planes from the rhombohedral setting if trigonal
    if sg.get_crystal_system() == "trigonal":
        transf = sg.get_conventional_to_primitive_transformation_matrix()
        miller_list = [hkl_transformation(transf, hkl) for hkl in conv_hkl_list]
        prim_structure = SpacegroupAnalyzer(structure).get_primitive_standard_structure()
        symm_ops = get_recp_symmetry_operation(prim_structure)
    else:
        miller_list = conv_hkl_list
        symm_ops = get_recp_symmetry_operation(structure)

    unique_millers, unique_millers_conv = [], []

    def is_already_analyzed(miller_index):
        for op in symm_ops:
            if in_coord_list(unique_millers, op.operate(miller_index)):
                return True
        return False

    for i, miller in enumerate(miller_list):
        d = abs(reduce(gcd, miller))
        miller = tuple([int(i / d) for i in miller])
        if not is_already_analyzed(miller):
            if sg.get_crystal_system() == "trigonal":
                # Now we find the distinct primitive hkls using
                # the primitive symmetry operations and their
                # corresponding hkls in the conventional setting
                unique_millers.append(miller)
                d = abs(reduce(gcd, conv_hkl_list[i]))
                cmiller = tuple([int(i / d) for i in conv_hkl_list[i]])
                unique_millers_conv.append(cmiller)
            else:
                unique_millers.append(miller)
                unique_millers_conv.append(miller)

    return unique_millers_conv