By linear operation of csl lattice vectors to get the best corresponding
        slab lattice. That is the area of a,b vectors (within the surface plane)
        is the smallest, the c vector first, has shortest length perpendicular
        to surface [h,k,l], second, has shortest length itself.

        Args:
            csl (3 by 3 integer array):
                    input csl lattice.
            surface (list of three integers, e.g. h, k, l):
                    the miller index of the surface, with the format of [h,k,l]
            normal (logic):
                    determine if the c vector needs to perpendicular to surface
            trans_cry (3 by 3 array):
                    transform matrix from crystal system to orthogonal system
            max_search (int): max search for the GB lattice vectors that give the smallest GB
                lattice. If normal is true, also max search the GB c vector that perpendicular
                to the plane.
            quick_gen (bool): whether to quickly generate a supercell, no need to find the smallest
                cell if set to true.

        Returns:
            t_matrix: a slab lattice ( 3 by 3 integer array):

def slab_from_csl(csl, surface, normal, trans_cry, max_search=20, quick_gen=False):
        """
        By linear operation of csl lattice vectors to get the best corresponding
        slab lattice. That is the area of a,b vectors (within the surface plane)
        is the smallest, the c vector first, has shortest length perpendicular
        to surface [h,k,l], second, has shortest length itself.

        Args:
            csl (3 by 3 integer array):
                    input csl lattice.
            surface (list of three integers, e.g. h, k, l):
                    the miller index of the surface, with the format of [h,k,l]
            normal (logic):
                    determine if the c vector needs to perpendicular to surface
            trans_cry (3 by 3 array):
                    transform matrix from crystal system to orthogonal system
            max_search (int): max search for the GB lattice vectors that give the smallest GB
                lattice. If normal is true, also max search the GB c vector that perpendicular
                to the plane.
            quick_gen (bool): whether to quickly generate a supercell, no need to find the smallest
                cell if set to true.

        Returns:
            t_matrix: a slab lattice ( 3 by 3 integer array):
        """

        # set the transform matrix in real space
        trans = trans_cry
        # transform matrix in reciprocal space
        ctrans = np.linalg.inv(trans.T)

        t_matrix = csl.copy()
        # vectors constructed from csl that perpendicular to surface
        ab_vector = []
        # obtain the miller index of surface in terms of csl.
        miller = np.matmul(surface, csl.T)
        if reduce(gcd, miller) != 1:
            miller = [int(round(x / reduce(gcd, miller))) for x in miller]
        miller_nonzero = []
        # quickly generate a supercell, normal is not work in this way
        if quick_gen:
            scale_factor = []
            eye = np.eye(3, dtype=np.int)
            for i, j in enumerate(miller):
                if j == 0:
                    scale_factor.append(eye[i])
                else:
                    miller_nonzero.append(i)
            if len(scale_factor) < 2:
                index_len = len(miller_nonzero)
                for i in range(index_len):
                    for j in range(i + 1, index_len):
                        lcm_miller = lcm(miller[miller_nonzero[i]], miller[miller_nonzero[j]])
                        l = [0, 0, 0]
                        l[miller_nonzero[i]] = -int(round(lcm_miller / miller[miller_nonzero[i]]))
                        l[miller_nonzero[j]] = int(round(lcm_miller / miller[miller_nonzero[j]]))
                        scale_factor.append(l)
                        if len(scale_factor) == 2:
                            break
            t_matrix[0] = np.array(np.dot(scale_factor[0], csl))
            t_matrix[1] = np.array(np.dot(scale_factor[1], csl))
            t_matrix[2] = csl[miller_nonzero[0]]
            if abs(np.linalg.det(t_matrix)) > 1000:
                warnings.warn('Too large matrix. Suggest to use quick_gen=False')
            return t_matrix

        for i, j in enumerate(miller):
            if j == 0:
                ab_vector.append(csl[i])
            else:
                c_index = i
                miller_nonzero.append(j)

        if len(miller_nonzero) > 1:
            t_matrix[2] = csl[c_index]
            index_len = len(miller_nonzero)
            lcm_miller = []
            for i in range(index_len):
                for j in range(i + 1, index_len):
                    com_gcd = gcd(miller_nonzero[i], miller_nonzero[j])
                    mil1 = int(round(miller_nonzero[i] / com_gcd))
                    mil2 = int(round(miller_nonzero[j] / com_gcd))
                    lcm_miller.append(max(abs(mil1), abs(mil2)))
            lcm_sorted = sorted(lcm_miller)
            if index_len == 2:
                max_j = lcm_sorted[0]
            else:
                max_j = lcm_sorted[1]
        else:
            if not normal:
                t_matrix[0] = ab_vector[0]
                t_matrix[1] = ab_vector[1]
                t_matrix[2] = csl[c_index]
                return t_matrix
            else:
                max_j = abs(miller_nonzero[0])
        if max_j > max_search:
            max_j = max_search
        # area of a, b vectors
        area = None
        # length of c vector
        c_norm = np.linalg.norm(np.matmul(t_matrix[2], trans))
        # c vector length along the direction perpendicular to surface
        c_length = np.abs(np.dot(t_matrix[2], surface))
        # check if the init c vector perpendicular to the surface
        if normal:
            c_cross = np.cross(np.matmul(t_matrix[2], trans), np.matmul(surface, ctrans))
            if np.linalg.norm(c_cross) < 1.e-8:
                normal_init = True
            else:
                normal_init = False

        j = np.arange(0, max_j + 1)
        combination = []
        for i in itertools.product(j, repeat=3):
            if sum(abs(np.array(i))) != 0:
                combination.append(list(i))
            if len(np.nonzero(i)[0]) == 3:
                for i1 in range(3):
                    new_i = list(i).copy()
                    new_i[i1] = -1 * new_i[i1]
                    combination.append(new_i)
            elif len(np.nonzero(i)[0]) == 2:
                new_i = list(i).copy()
                new_i[np.nonzero(i)[0][0]] = -1 * new_i[np.nonzero(i)[0][0]]
                combination.append(new_i)
        for i in combination:
            if reduce(gcd, i) == 1:
                temp = np.dot(np.array(i), csl)
                if abs(np.dot(temp, surface) - 0) < 1.e-8:
                    ab_vector.append(temp)
                else:
                    # c vector length along the direction perpendicular to surface
                    c_len_temp = np.abs(np.dot(temp, surface))
                    # c vector length itself
                    c_norm_temp = np.linalg.norm(np.matmul(temp, trans))
                    if normal:
                        c_cross = np.cross(np.matmul(temp, trans), np.matmul(surface, ctrans))
                        if np.linalg.norm(c_cross) < 1.e-8:
                            if normal_init:
                                if c_norm_temp < c_norm:
                                    t_matrix[2] = temp
                                    c_norm = c_norm_temp
                            else:
                                c_norm = c_norm_temp
                                normal_init = True
                                t_matrix[2] = temp
                    else:
                        if c_len_temp < c_length or \
                                (abs(c_len_temp - c_length) < 1.e-8 and c_norm_temp < c_norm):
                            t_matrix[2] = temp
                            c_norm = c_norm_temp
                            c_length = c_len_temp

        if normal and (not normal_init):
            logger.info('Did not find the perpendicular c vector, increase max_j')
            while (not normal_init):
                if max_j == max_search:
                    warnings.warn('Cannot find the perpendicular c vector, please increase max_search')
                    break
                max_j = 3 * max_j
                if max_j > max_search:
                    max_j = max_search
                j = np.arange(0, max_j + 1)
                combination = []
                for i in itertools.product(j, repeat=3):
                    if sum(abs(np.array(i))) != 0:
                        combination.append(list(i))
                    if len(np.nonzero(i)[0]) == 3:
                        for i1 in range(3):
                            new_i = list(i).copy()
                            new_i[i1] = -1 * new_i[i1]
                            combination.append(new_i)
                    elif len(np.nonzero(i)[0]) == 2:
                        new_i = list(i).copy()
                        new_i[np.nonzero(i)[0][0]] = -1 * new_i[np.nonzero(i)[0][0]]
                        combination.append(new_i)
                for i in combination:
                    if reduce(gcd, i) == 1:
                        temp = np.dot(np.array(i), csl)
                        if abs(np.dot(temp, surface) - 0) > 1.e-8:
                            c_cross = np.cross(np.matmul(temp, trans), np.matmul(surface, ctrans))
                            if np.linalg.norm(c_cross) < 1.e-8:
                                # c vetor length itself
                                c_norm_temp = np.linalg.norm(np.matmul(temp, trans))
                                if normal_init:
                                    if c_norm_temp < c_norm:
                                        t_matrix[2] = temp
                                        c_norm = c_norm_temp
                                else:
                                    c_norm = c_norm_temp
                                    normal_init = True
                                    t_matrix[2] = temp
                if normal_init:
                    logger.info('Found perpendicular c vector')

        # find the best a, b vectors with their formed area smallest and average norm of a,b smallest.
        for i in itertools.combinations(ab_vector, 2):
            area_temp = np.linalg.norm(np.cross(np.matmul(i[0], trans),
                                                np.matmul(i[1], trans)))
            if abs(area_temp - 0) > 1.e-8:
                ab_norm_temp = np.linalg.norm(np.matmul(i[0], trans)) + \
                               np.linalg.norm(np.matmul(i[1], trans))
                if area is None:
                    area = area_temp
                    ab_norm = ab_norm_temp
                    t_matrix[0] = i[0]
                    t_matrix[1] = i[1]
                elif area_temp < area:
                    t_matrix[0] = i[0]
                    t_matrix[1] = i[1]
                    area = area_temp
                    ab_norm = ab_norm_temp
                elif abs(area - area_temp) < 1.e-8 and ab_norm_temp < ab_norm:
                    t_matrix[0] = i[0]
                    t_matrix[1] = i[1]
                    area = area_temp
                    ab_norm = ab_norm_temp

        # make sure we have a left-handed crystallographic system
        if np.linalg.det(np.matmul(t_matrix, trans)) < 0:
            t_matrix *= -1

        if normal and abs(np.linalg.det(t_matrix)) > 1000:
            warnings.warn('Too large matrix. Suggest to use Normal=False')
        return t_matrix