How to use the pymatgen.core.operations.SymmOp function in pymatgen

# C1 symmetry breaks assumptions in the algorithm afterwards
        return symmops
    else:
        full = list(generators)

        for g in full:
            for s in generators:
                op = np.dot(g, s)
                d = np.abs(full - op) < tol
                if not np.any(np.all(np.all(d, axis=2), axis=1)):
                    full.append(op)

        d = np.abs(full - UNIT) < tol
        if not np.any(np.all(np.all(d, axis=2), axis=1)):
            full.append(UNIT)
        return [SymmOp(op) for op in full]

# Normalize the inertia tensor so that it does not scale with size
            # of the system.  This mitigates the problem of choosing a proper
            # comparison tolerance for the eigenvalues.
            inertia_tensor /= total_inertia
            eigvals, eigvecs = np.linalg.eig(inertia_tensor)
            self.principal_axes = eigvecs.T
            self.eigvals = eigvals
            v1, v2, v3 = eigvals
            eig_zero = abs(v1 * v2 * v3) < self.eig_tol ** 3
            eig_all_same = abs(v1 - v2) < self.eig_tol and abs(
                v1 - v3) < self.eig_tol
            eig_all_diff = abs(v1 - v2) > self.eig_tol and abs(
                v1 - v2) > self.eig_tol and abs(v2 - v3) > self.eig_tol

            self.rot_sym = []
            self.symmops = [SymmOp(np.eye(4))]

            if eig_zero:
                logger.debug("Linear molecule detected")
                self._proc_linear()
            elif eig_all_same:
                logger.debug("Spherical top molecule detected")
                self._proc_sph_top()
            elif eig_all_diff:
                logger.debug("Asymmetric top molecule detected")
                self._proc_asym_top()
            else:
                logger.debug("Symmetric top molecule detected")
                self._proc_sym_top()

def _check_rot_sym(self, axis):
        """
        Determines the rotational symmetry about supplied axis.  Used only for
        symmetric top molecules which has possible rotational symmetry
        operations > 2.
        """
        min_set = self._get_smallest_set_not_on_axis(axis)
        max_sym = len(min_set)
        for i in xrange(max_sym, 0, -1):
            if max_sym % i != 0:
                continue
            op = SymmOp.from_axis_angle_and_translation(axis, 360 / i)
            rotvalid = self.is_valid_op(op)
            if rotvalid:
                self.symmops.append(op)
                self.rot_sym.append((axis, i))
                return i
        return 1

Args:
		structure (pymatgen.core.structure.Structure)
		symprec (number): symmetry precision for pymatgen SpacegroupAnalyzer
	"""

	sga = SpacegroupAnalyzer(structure, symprec * max(structure.lattice.abc))
	symmops = sga.get_symmetry_operations(cartesian = True)
	lattice = structure.lattice.matrix
	invlattice = structure.lattice.inv_matrix
	newops = []
	for op in symmops:
		newrot = np.dot(lattice, op.rotation_matrix)
		newrot = np.dot(newrot, invlattice)
		newtrans = np.dot(op.translation_vector, invlattice)
		newops.append(SymmOp.from_rotation_and_translation(
			newrot, newtrans))
	return newops

def _test_rotation(self, rot, origin, fixed, to_fit, tol_atoms):
        tol_atoms_plus = 1.1 * tol_atoms
        found_map = False
        mapping_op = None
        biggest_dist = 0
        logger.debug("Trying candidate rotation : \n" + str(rot))
        for site in fixed:
            if time.time() - self._start_time > self._timeout:
                logger.debug("Timeout reached when testing rotations.")
                break
            if site.species_and_occu == origin.species_and_occu:
                shift = site.coords
                op = SymmOp.from_rotation_and_translation(
                                                    rot.rotation_matrix, shift)
                nstruct = apply_operation(to_fit, op)
                correspondance = OrderedDict()
                all_match = True
                biggest_dist = 0
                # check to see if transformed struct matches fixed structure
                for trans in nstruct:
                    cands = fixed.get_sites_in_sphere(trans.coords,
                                                      tol_atoms_plus)
                    if len(cands) == 0:
                        logger.debug("No candidates found.")
                        all_match = False
                        break
                    cands = sorted(cands, key=lambda a: a[1])
                    (closest, closest_dist) = cands[0]
                    if closest_dist > tol_atoms or \

Takes a symmetry operation and transforms it.
        :param symmop: SymmOp or MagSymmOp
        :return: 
        """
        W = symmop.rotation_matrix
        w = symmop.translation_vector
        Q = np.linalg.inv(self.P)
        W_ = np.matmul(np.matmul(Q, W), self.P)
        I = np.identity(3)
        w_ = np.matmul(Q, (w + np.matmul(W - I, self.p)))
        if isinstance(symmop, MagSymmOp):
            return MagSymmOp.from_rotation_and_translation_and_time_reversal(
                rotation_matrix=W_, translation_vec=w_,
                time_reversal=symmop.time_reversal, tol=symmop.tol)
        elif isinstance(symmop, SymmOp):
            return SymmOp.from_rotation_and_translation(
                rotation_matrix=W_, translation_vec=w_, tol=symmop.tol)

yield (shells[0][x], shells[1][y], shells[2][z])
            test_rotations = random_rot()
        det = np.linalg.det
        inv = np.linalg.inv
        fixed_vol = fixed.volume
        for pool in test_rotations:
            if all([nn.species_and_occu == origin.species_and_occu \
                    for nn in pool]):
                # Can a unitary transformation bring the cell vectors together?
                coords = [nn.coords - origin.coords for nn in pool]
                cell_v = np.array(coords).transpose()
                d = det(cell_v)
                if abs(d) < 0.001 or abs(abs(d) - fixed_vol) > 0.001:
                    continue
                rot = np.dot(fixed_basis, inv(cell_v))
                r = SymmOp.from_rotation_and_translation(rot,
                                                         np.array([0, 0, 0]))

                if r not in cand_rot:
                    transf = r.rotation_matrix
                    transf = np.dot(transf.transpose(), transf)

                    # This resolves a very very strange bug in numpy that
                    # causes random eigenvectors to be returned for matrices
                    # very very close to the identity matrix.
                    transf = np.eye(3) if np.allclose(transf, np.eye(3)) \
                                       else transf
                    pbis = sqrt_matrix(transf)
                    if shear_invariant(pbis) < tol_shear:
                        cand_rot[r] = shear_invariant(pbis)
            if time.time() - self._start_time > self._timeout:
                logger.debug("Timeout reached when generating rotations.")

def site_symm(point, gen_pos, tol=1e-3, lattice=Euclidean_lattice):
    '''
    Given gen_pos (a list of SymmOps), return the list of symmetry operations
    leaving a point (coordinate or SymmOp) invariant.
    '''
    #Convert point into a SymmOp
    if type(point) != SymmOp:
        point = SymmOp.from_rotation_and_translation([[0,0,0],[0,0,0],[0,0,0]], point)
    symmetry = []
    for op in gen_pos:
        is_symmetry = True
        #Calculate the effect of applying op to point
        difference = SymmOp((op*point).affine_matrix - point.affine_matrix)
        #Check that the rotation matrix is unaltered by op
        if not np.allclose(difference.rotation_matrix, np.zeros((3,3)), rtol = 1e-3, atol = 1e-3):
            is_symmetry = False
        #Check that the displacement is less than tol
        displacement = difference.translation_vector
        if distance(displacement, lattice) > tol:
            is_symmetry = False
        if is_symmetry:
            '''The actual site symmetry's translation vector may vary from op by
            a factor of +1 or -1 (especially when op contains +-1/2).
            We record this to distinguish between special Wyckoff positions.
            As an example, consider the point (-x+1/2,-x,x+1/2) in position 16c
            of space group Ia-3(206). The site symmetry includes the operations
            (-z+1,x-1/2,-y+1/2) and (y+1/2,-z+1/2,-x+1). These operations are
            not listed in the general position, but correspond to the operations
            (-z,x+1/2,-y+1/2) and (y+1/2,-z+1/2,-x), respectively, just shifted

hkl=hkl_pair['hkl'][1],
                      min_thick=hkl_pair['thickness'][1],
                      min_vac=0,
                      primitive=False, from_ase=True)

    # find top atoms reference of lower part of gb
    substrate_top_z = np.max(np.array([site.coords for site in lower])[:, 2])

    # define twist and tilt vectors
    twist_shift_normal = lower.lattice.matrix[2, :] / \
                         np.linalg.norm(lower.lattice.matrix[2, :])
    tilt_normal = upper.lattice.matrix[1, :] / \
                  np.linalg.norm(upper.lattice.matrix[2, :])

    # define twist operation SymmOp object
    twist_op = SymmOp.from_axis_angle_and_translation(axis=twist_shift_normal, \
                                                      angle=twist,
                                                      angle_in_radians=False,
                                                      translation_vec=(
                                                          0, 0, 0))
    # define tilt operation SymmOp object
    tilt_op = SymmOp.from_axis_angle_and_translation(axis=tilt_normal, \
                                                     angle=tilt,
                                                     angle_in_radians=False,
                                                     translation_vec=(0, 0, 0))
    upper.apply_operation(twist_op)
    upper.to(fmt="poscar", filename="POSCAR_upper.vasp")
    upper.apply_operation(tilt_op)

    # define shift separation along twist vector normal to upper plane
    shift = -1 * twist_shift_normal / np.linalg.norm(
        twist_shift_normal) * separation

val = float(data)
                            parameters.append(val)
                        except ValueError:
                            if data.startswith("-"):
                                parameters.append(-paras[data[1:]])
                            else:
                                parameters.append(paras[data])
                    if len(nn) == 1:
                        coords.append(np.array([0, 0, parameters[0]]))
                    elif len(nn) == 2:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        axis = [0, 1, 0]
                        op = SymmOp.from_origin_axis_angle(coords1, axis,
                                                           angle, False)
                        coord = op.operate(coords2)
                        vec = coord - coords1
                        coord = vec * bl / np.linalg.norm(vec) + coords1
                        coords.append(coord)
                    elif len(nn) == 3:
                        coords1 = coords[nn[0] - 1]
                        coords2 = coords[nn[1] - 1]
                        coords3 = coords[nn[2] - 1]
                        bl = parameters[0]
                        angle = parameters[1]
                        dih = parameters[2]
                        v1 = coords3 - coords2
                        v2 = coords1 - coords2
                        axis = np.cross(v1, v2)
                        op = SymmOp.from_origin_axis_angle(