def determine_symmetry_positions(st, element):
    """
    determine non-equivalent positions for atoms of type *element*

    element (str) - name of element, for example Li

    return list of lists -  atom numbers for each non-equivalent position
    """

    from pymatgen.symmetry.analyzer import SpacegroupAnalyzer


    stp = st.convert2pymatgen()

    spg = SpacegroupAnalyzer(stp)

    info = spg.get_symmetry_dataset()

    positions = {}
    for i, (el, pos) in enumerate(zip(st.get_elements(), info['equivalent_atoms'])):
        
        if el == element and pos not in positions:
            positions[pos] = []

        if el == element:
            positions[pos].append(i)


    printlog('I have found ', len(positions), 'non-equivalent positions for', element, ':',positions.keys(), imp = 'y', end = '\n')
    printlog('Atom numbers: ', positions, imp = 'y')
    
    sorted_keys = sorted(list(positions.keys()))
    pos_lists = [positions[key] for key in sorted_keys ]

    return pos_lists

    def test_partial_disorder(self):
        s = Structure.from_file(filename=os.path.join(test_dir, "garnet.cif"))
        a = SpacegroupAnalyzer(s, 0.1)
        prim = a.find_primitive()
        s = prim.copy()
        s["Al3+"] = {"Al3+": 0.5, "Ga3+": 0.5}
        adaptor = EnumlibAdaptor(s, 1, 1, enum_precision_parameter=0.01)
        adaptor.run()
        structures = adaptor.structures
        self.assertEqual(len(structures), 7)
        for s in structures:
            self.assertEqual(s.formula, 'Ca12 Al4 Ga4 Si12 O48')
        s = prim.copy()
        s["Ca2+"] = {"Ca2+": 1/3, "Mg2+": 2/3}
        adaptor = EnumlibAdaptor(s, 1, 1, enum_precision_parameter=0.01)
        adaptor.run()
        structures = adaptor.structures
        self.assertEqual(len(structures), 20)
        for s in structures:
            self.assertEqual(s.formula, 'Ca4 Mg8 Al8 Si12 O48')

        s = prim.copy()
        s["Si4+"] = {"Si4+": 1/3, "Ge4+": 2/3}
        adaptor = EnumlibAdaptor(s, 1, 1, enum_precision_parameter=0.01)
        adaptor.run()
        structures = adaptor.structures
        self.assertEqual(len(structures), 18)
        for s in structures:
            self.assertEqual(s.formula, 'Ca12 Al8 Si4 Ge8 O48')

    def symm_reduce(self, coords_set, threshold=1e-6):
        """
        Reduces the set of adsorbate sites by finding removing
        symmetrically equivalent duplicates

        Args:
            coords_set: coordinate set in cartesian coordinates
            threshold: tolerance for distance equivalence, used
                as input to in_coord_list_pbc for dupl. checking
        """
        surf_sg = SpacegroupAnalyzer(self.slab, 0.1)
        symm_ops = surf_sg.get_symmetry_operations()
        unique_coords = []
        # Convert to fractional
        coords_set = [self.slab.lattice.get_fractional_coords(coords)
                      for coords in coords_set]
        for coords in coords_set:
            incoord = False
            for op in symm_ops:
                if in_coord_list_pbc(unique_coords, op.operate(coords),
                                     atol=threshold):
                    incoord = True
                    break
            if not incoord:
                unique_coords += [coords]
        # convert back to cartesian
        return [self.slab.lattice.get_cartesian_coords(coords)
                for coords in unique_coords]

def unique_symmetry_operations_as_vectors_from_structure( structure, subset = None ):
    """
    Uses pymatgen symmetry analysis to find the minimum complete set of symmetry operations for the space group of a structure.

    Args:
        structure (pymatgen Structure): structure to be analysed.
        subset    (Optional [list]):    list of atom indices to be used for generating the symmetry operations

    Returns:
        a list of lists, containing the symmetry operations as vector mappings.
    """
    symmetry_analyzer = SpacegroupAnalyzer( structure )

    print( "The spacegroup for structure is {}".format(symmetry_analyzer.get_spacegroup_symbol()) )

    symmetry_operations = symmetry_analyzer.get_symmetry_operations()

    mappings = []
    if subset:
        species_subset = [ spec for i,spec in enumerate(structure.species) if i in subset]
        frac_coords_subset = [ coord for i, coord in enumerate( structure.frac_coords ) if i in subset ]
        mapping_structure = Structure( structure.lattice, species_subset, frac_coords_subset ) 
    else:
        mapping_structure = structure

    for symmop in symmetry_operations:
        new_structure = Structure( mapping_structure.lattice, mapping_structure.species, symmop.operate_multi( mapping_structure.frac_coords ) )
        new_mapping = [ x+1 for x in list( coord_list_mapping_pbc( new_structure.frac_coords, mapping_structure.frac_coords ) ) ]
        if new_mapping not in mappings:
            mappings.append( new_mapping )

    return mappings

    def symmetrize(self, structure):
        tensor = self._reduced_tensor

        if self._is_real_space:
            real_lattice = self._lattice
        else:
            real_lattice = self._lattice.reciprocal_lattice

        # I guess this is the reason why tensor.symmetrize (omega) is so slow!
        real_finder = SpacegroupAnalyzer(structure)

        real_symmops = real_finder.get_point_group_operations(cartesian=True)

        cartesian_tensor = self.cartesian_tensor

        sym_tensor = np.zeros((3,3))

        my_tensor = cartesian_tensor

        for real_sym in real_symmops:
            mat = real_sym.rotation_matrix
            prod_sym = np.dot(np.transpose(mat),np.dot(cartesian_tensor,mat))
            sym_tensor = sym_tensor + prod_sym

        sym_tensor = sym_tensor/len(real_symmops)

        self._reduced_tensor = from_cart_to_red(sym_tensor,self._lattice)

def symmetry_reduce(tensors, structure, tol=1e-8, **kwargs):
    """
    Function that converts a list of tensors corresponding to a structure
    and returns a dictionary consisting of unique tensor keys with symmop
    values corresponding to transformations that will result in derivative
    tensors from the original list

    Args:
        tensors (list of tensors): list of Tensor objects to test for
            symmetrically-equivalent duplicates
        structure (Structure): structure from which to get symmetry
        tol (float): tolerance for tensor equivalence
        kwargs: keyword arguments for the SpacegroupAnalyzer

    returns:
        dictionary consisting of unique tensors with symmetry operations
        corresponding to those which will reconstruct the remaining
        tensors as values
    """
    sga = SpacegroupAnalyzer(structure, **kwargs)
    symmops = sga.get_symmetry_operations(cartesian=True)
    unique_tdict = {}
    for tensor in tensors:
        is_unique = True
        for unique_tensor, symmop in itertools.product(unique_tdict, symmops):
            if (np.abs(unique_tensor.transform(symmop) - tensor) < tol).all():
                unique_tdict[unique_tensor].append(symmop)
                is_unique = False
                break
        if is_unique:
            unique_tdict[tensor] = []
    return unique_tdict

    def multiplicity(self):
        """
        Returns the multiplicity of a defect site within the structure (needed for concentration analysis)
        """
        if self._multiplicity is None:
            # generate multiplicity based on space group symmetry operations performed on defect coordinates
            try:
                d_structure = create_saturated_interstitial_structure(self)
            except ValueError:
                logger.debug('WARNING! Multiplicity was not able to be calculated adequately '
                             'for interstitials...setting this to 1 and skipping for now...')
                return 1

            sga = SpacegroupAnalyzer(d_structure)
            periodic_struc = sga.get_symmetrized_structure()
            poss_deflist = sorted(
                periodic_struc.get_sites_in_sphere(self.site.coords, 2, include_index=True),
                key=lambda x: x[1])
            defindex = poss_deflist[0][2]

            equivalent_sites = periodic_struc.find_equivalent_sites(periodic_struc[defindex])
            return len(equivalent_sites)

        else:
            return self._multiplicity

    def __init__(self, structure, min_cell_size=1, max_cell_size=1,
                 symm_prec=0.1, enum_precision_parameter=0.001,
                 refine_structure=False):
        """
        Initializes the adapter with a structure and some parameters.

        Args:
            structure: An input structure.
            min_cell_size (int): The minimum cell size wanted. Defaults to 1.
            max_cell_size (int): The maximum cell size wanted. Defaults to 1.
            symm_prec (float): Symmetry precision. Defaults to 0.1.
            enum_precision_parameter (float): Finite precision parameter for
                enumlib. Default of 0.001 is usually ok, but you might need to
                tweak it for certain cells.
            refine_structure (bool): If you are starting from a structure that
                has been relaxed via some electronic structure code,
                it is usually much better to start with symmetry determination
                and then obtain a refined structure. The refined structure have
                cell parameters and atomic positions shifted to the expected
                symmetry positions, which makes it much less sensitive precision
                issues in enumlib. If you are already starting from an
                experimental cif, refinement should have already been done and
                it is not necessary. Defaults to False.
        """
        if refine_structure:
            finder = SpacegroupAnalyzer(structure, symm_prec)
            self.structure = finder.get_refined_structure()
        else:
            self.structure = structure
        self.min_cell_size = min_cell_size
        self.max_cell_size = max_cell_size
        self.symm_prec = symm_prec
        self.enum_precision_parameter = enum_precision_parameter

def get_unique_site_indices(structure):
    sa = SpacegroupAnalyzer(structure)
    symm_data = sa.get_symmetry_dataset()
    # equivalency mapping for the structure
    # i'th site in the struct equivalent to eq_struct[i]'th site
    eq_atoms = symm_data["equivalent_atoms"]
    return np.unique(eq_atoms).tolist()

    def symm_check(self, ucell, wulff_vertices):
        """
        # Checks if the point group of the Wulff shape matches
        # the point group of its conventional unit cell

        Args:
            ucell (string): Unit cell that the Wulff shape is based on.
            wulff_vertices (list): List of all vertices on the Wulff
                shape. Use wulff.wulff_pt_list to obtain the list
                (see wulff_generator.py).

        return (bool)
        """

        space_group_analyzer = SpacegroupAnalyzer(ucell)
        symm_ops = space_group_analyzer.get_point_group_operations(
            cartesian=True)
        for point in wulff_vertices:
            for op in symm_ops:
                symm_point = op.operate(point)
                if in_coord_list(wulff_vertices, symm_point):
                    continue
                else:
                    return False
        return True

 def set_space_group_from_structure(self, structure):
     spga = SpacegroupAnalyzer(structure=structure)
     self.crystal_system = spga.get_crystal_system()
     self.hall = spga.get_hall()
     self.number = spga.get_space_group_number()
     self.source = "spglib"
     self.symbol = spga.get_space_group_symbol()

    def __init__(self, structure, element):
        """
        Initializes a Substitution Generator
        note: an Antisite is considered a type of substitution
        Args:
            structure(Structure): pymatgen structure object
            element (str or Element or Specie): element for the substitution
        """
        self.structure = structure
        self.element = element

        # Find equivalent site list
        sga = SpacegroupAnalyzer(self.structure)
        self.symm_structure = sga.get_symmetrized_structure()

        self.equiv_sub = []
        for equiv_site_set in list(self.symm_structure.equivalent_sites):
            vac_site = equiv_site_set[0]
            if isinstance(element, str):  # make sure you compare with specie symbol or Element type
                vac_specie = vac_site.specie.symbol
            else:
                vac_specie = vac_site.specie
            if element != vac_specie:
                defect_site = PeriodicSite(element, vac_site.coords, structure.lattice, coords_are_cartesian=True)
                sub = Substitution(structure, defect_site)
                self.equiv_sub.append(sub)

def print_spg(src="POSCAR"):
    srcpos = Poscar.from_file(src)
    finder = SpacegroupAnalyzer(srcpos.structure, symprec=5e-2, angle_tolerance=8)
    spg = finder.get_spacegroup_symbol()
    spg_num = finder.get_spacegroup_number()
    print(spg)
    print(spg_num)

    def test_tensor(self):
        """Initialize Tensor"""

        lattice = Lattice.hexagonal(4,6)
        #rprimd = np.array([[0,0.5,0.5],[0.5,0,0.5],[0.5,0.5,0]])
        #rprimd = rprimd*10
        #lattice = Lattice(rprimd)
        structure = Structure(lattice, ["Ga", "As"],
                                      [[0, 0, 0], [0.5, 0.5, 0.5]])

        #finder = SymmetryFinder(structure)
        finder = SpacegroupAnalyzer(structure)

        spacegroup = finder.get_spacegroup()
        pointgroup = finder.get_point_group()

        cartesian_tensor = [[2,3,1.2],[3,4,1.0],[1.2,1.0,6]]

        tensor = Tensor.from_cartesian_tensor(cartesian_tensor,lattice.reciprocal_lattice,space="g")
        red_tensor = tensor.reduced_tensor
        tensor2 = Tensor(red_tensor,lattice.reciprocal_lattice,space="g")
        assert(((np.abs(tensor2.cartesian_tensor)-np.abs(cartesian_tensor)) < 1E-8).all())

        self.assertTrue(tensor==tensor2)
        print(tensor)

        #print("non-symmetrized cartesian_tensor = ",tensor2.cartesian_tensor)
        tensor2.symmetrize(structure)

        #print("symmetrized_cartesian_tensor = ",tensor2.cartesian_tensor)

        self.serialize_with_pickle(tensor)

    def test_adsorb_both_surfaces(self):

        # Test out for monatomic adsorption
        o = Molecule("O", [[0, 0, 0]])
        adslabs = self.asf_100.adsorb_both_surfaces(o)
        adslabs_one = self.asf_100.generate_adsorption_structures(o)
        self.assertEqual(len(adslabs), len(adslabs_one))
        for adslab in adslabs:
            sg = SpacegroupAnalyzer(adslab)
            sites = sorted(adslab, key=lambda site: site.frac_coords[2])
            self.assertTrue(sites[0].species_string == "O")
            self.assertTrue(sites[-1].species_string == "O")
            self.assertTrue(sg.is_laue())

        # Test out for molecular adsorption
        oh = Molecule(["O", "H"], [[0, 0, 0], [0, 0, 1]])
        adslabs = self.asf_100.adsorb_both_surfaces(oh)
        adslabs_one = self.asf_100.generate_adsorption_structures(oh)
        self.assertEqual(len(adslabs), len(adslabs_one))
        for adslab in adslabs:
            sg = SpacegroupAnalyzer(adslab)
            sites = sorted(adslab, key=lambda site: site.frac_coords[2])
            self.assertTrue(sites[0].species_string in ["O", "H"])
            self.assertTrue(sites[-1].species_string in ["O", "H"])
            self.assertTrue(sg.is_laue())

    def get_sym_eq_kpoints(self, kpoint, cartesian=False, tol=1e-2):
        """
        Returns a list of unique symmetrically equivalent k-points.

        Args:
            kpoint (1x3 array): coordinate of the k-point
            cartesian (bool): kpoint is in cartesian or fractional coordinates
            tol (float): tolerance below which coordinates are considered equal

        Returns:
            ([1x3 array] or None): if structure is not available returns None
        """
        if not self.structure:
            return None
        sg = SpacegroupAnalyzer(self.structure)
        symmops = sg.get_point_group_operations(cartesian=cartesian)
        points = np.dot(kpoint, [m.rotation_matrix for m in symmops])
        rm_list = []
        # identify and remove duplicates from the list of equivalent k-points:
        for i in range(len(points) - 1):
            for j in range(i + 1, len(points)):
                if np.allclose(pbc_diff(points[i], points[j]), [0, 0, 0], tol):
                    rm_list.append(i)
                    break
        return np.delete(points, rm_list, axis=0)

    def from_structure(cls, structure, has_timerev=True, symprec=1e-5, angle_tolerance=5):
        """
        Takes a :class:`Structure` object. Uses spglib to perform various symmetry finding operations.

        Args:
            structure: :class:`Structure` object
            has_timerev: True is time-reversal symmetry is included.
            symprec: Tolerance for symmetry finding
            angle_tolerance: Angle tolerance for symmetry finding.

        .. warning::

            AFM symmetries are not supported.
        """
        # Call spglib to get the list of symmetry operations.
        spga = SpacegroupAnalyzer(structure, symprec=symprec, angle_tolerance=angle_tolerance)
        data = spga.get_symmetry_dataset()

        return cls(
            spgid=data["number"],
            symrel=data["rotations"],
            tnons=data["translations"],
            symafm=len(symrel) * [1],
            has_timerev=has_timerev,
            inord="C",
        )

    def test_structure_transform(self):
        # Test trivial case
        trivial = self.fit_r4.structure_transform(self.structure,
                                                  self.structure.copy())
        self.assertArrayAlmostEqual(trivial, self.fit_r4)

        # Test simple rotation
        rot_symm_op = SymmOp.from_axis_angle_and_translation([1, 1, 1], 55.5)
        rot_struct = self.structure.copy()
        rot_struct.apply_operation(rot_symm_op)
        rot_tensor = self.fit_r4.rotate(rot_symm_op.rotation_matrix)
        trans_tensor = self.fit_r4.structure_transform(self.structure, rot_struct)
        self.assertArrayAlmostEqual(rot_tensor, trans_tensor)

        # Test supercell
        bigcell = self.structure.copy()
        bigcell.make_supercell([2, 2, 3])
        trans_tensor = self.fit_r4.structure_transform(self.structure, bigcell)
        self.assertArrayAlmostEqual(self.fit_r4, trans_tensor)

        # Test rotated primitive to conventional for fcc structure
        sn = self.get_structure("Sn")
        sn_prim = SpacegroupAnalyzer(sn).get_primitive_standard_structure()
        sn_prim.apply_operation(rot_symm_op)
        rotated = self.fit_r4.rotate(rot_symm_op.rotation_matrix)
        transformed = self.fit_r4.structure_transform(sn, sn_prim)
        self.assertArrayAlmostEqual(rotated, transformed)

    def set_miller_family(self):
        """
        get all miller indices for the given maximum index
        get the list of indices that correspond to the given family
        of indices
        """
        recp_structure = Structure(self.recp_lattice, ["H"], [[0, 0, 0]])
        analyzer = SpacegroupAnalyzer(recp_structure, symprec=0.001)
        symm_ops = analyzer.get_symmetry_operations()

        max_index = max(max(m) for m in self.hkl_family)
        r = list(range(-max_index, max_index + 1))
        r.reverse()
        miller_indices = []
        self.all_equiv_millers = []
        self.all_surface_energies = []
        for miller in itertools.product(r, r, r):
            if any([i != 0 for i in miller]):
                d = abs(reduce(gcd, miller))
                miller_index = tuple([int(i / d) for i in miller])
                for op in symm_ops:
                    for i, u_miller in enumerate(self.hkl_family):
                        if in_coord_list(u_miller, op.operate(miller_index)):
                            self.all_equiv_millers.append(miller_index)
                            self.all_surface_energies.append(self.surface_energies[i])

    def set_material_data_from_structure(self, structure, space_group=True, symprec=1e-3, angle_tolerance=5):
        """
        Sets the fields of the Document using a Structure and Spglib to determine the space group properties

        Args:
            structure: A |Structure|
            space_group: if True sets the spacegroup fields using spglib_.
            symprec (float): Tolerance for symmetry finding.
            angle_tolerance (float): Angle tolerance for symmetry finding.
        """

        comp = structure.composition
        el_amt = structure.composition.get_el_amt_dict()
        self.unit_cell_formula = comp.as_dict()
        self.reduced_cell_formula = comp.to_reduced_dict
        self.elements = list(el_amt.keys())
        self.nelements = len(el_amt)
        self.pretty_formula = comp.reduced_formula
        self.anonymous_formula = comp.anonymized_formula
        self.nsites = comp.num_atoms
        self.chemsys = "-".join(sorted(el_amt.keys()))
        if space_group:
            sym = SpacegroupAnalyzer(structure, symprec=symprec, angle_tolerance=angle_tolerance)
            self.spacegroup = SpaceGroupDocument(crystal_system=sym.get_crystal_system(), hall=sym.get_hall(),
                                                 number=sym.get_space_group_number(), point_group=sym.get_point_group_symbol(),
                                                 symbol=sym.get_space_group_symbol(), source="spglib")