def get_smallest_vectors(supercell_bases,
                         supercell_pos,
                         primitive_pos,
                         symprec=1e-5):
    """Find shortest atomic pair vectors

    Note
    ----
    Shortest vectors from an atom in primitive cell to an atom in
    supercell in the fractional coordinates of primitive cell. If an
    atom in supercell is on the border centered at an atom in
    primitive and there are multiple vectors that have the same
    distance (up to tolerance) and different directions, several
    shortest vectors are stored.
    In fact, this method is not limited to search shortest vectors between
    sueprcell atoms and primitive cell atoms, but can be used to measure
    shortest vectors between atoms in periodic supercell lattice frame.

    Parameters
    ----------
    supercell_bases : ndarray
        Supercell basis vectors as row vectors, (a, b, c)^T. Must be order='C'.
        dtype='double'
        shape=(3, 3)
    supercell_pos : array_like
        Atomic positions in fractional coordinates of supercell.
        dtype='double'
        shape=(size_super, 3)
    primitive_pos : array_like
        Atomic positions in fractional coordinates of supercell. Note that not
        in fractional coodinates of primitive cell.
        dtype='double'
        shape=(size_prim, 3)
    symprec : float, optional, default=1e-5
        Tolerance to find equal distances of vectors

    Returns
    -------
    shortest_vectors : ndarray
        Shortest vectors in supercell coordinates. The 27 in shape is the
        possible maximum number of elements.
        dtype='double'
        shape=(size_super, size_prim, 27, 3)
    multiplicities : ndarray
        Number of equidistance shortest vectors
        dtype='intc'
        shape=(size_super, size_prim)

    """

    reduced_bases = get_reduced_bases(supercell_bases,
                                      method='delaunay',
                                      tolerance=symprec)
    trans_mat_float = np.dot(supercell_bases, np.linalg.inv(reduced_bases))
    trans_mat = np.rint(trans_mat_float).astype(int)
    assert (np.abs(trans_mat_float - trans_mat) < 1e-8).all()
    trans_mat_inv_float = np.linalg.inv(trans_mat)
    trans_mat_inv = np.rint(trans_mat_inv_float).astype(int)
    assert (np.abs(trans_mat_inv_float - trans_mat_inv) < 1e-8).all()

    # Reduce all positions into the cell formed by the reduced bases.
    supercell_fracs = np.dot(supercell_pos, trans_mat)
    supercell_fracs -= np.rint(supercell_fracs)
    supercell_fracs = np.array(supercell_fracs, dtype='double', order='C')
    primitive_fracs = np.dot(primitive_pos, trans_mat)
    primitive_fracs -= np.rint(primitive_fracs)
    primitive_fracs = np.array(primitive_fracs, dtype='double', order='C')

    # For each vector, we will need to consider all nearby images in the
    # reduced bases.
    lattice_points = np.array([[i, j, k]
                               for i in (-1, 0, 1)
                               for j in (-1, 0, 1)
                               for k in (-1, 0, 1)],
                              dtype='intc', order='C')

    # Here's where things get interesting.
    # We want to avoid manually iterating over all possible pairings of
    # supercell atoms and primitive atoms, because doing so creates a
    # tight loop in larger structures that is difficult to optimize.
    #
    # Furthermore, it seems wise to call numpy.dot on as large of an array
    # as possible, since numpy can shell out to BLAS to handle the
    # real heavy lifting.

    shortest_vectors = np.zeros(
        (len(supercell_fracs), len(primitive_fracs), 27, 3),
        dtype='double', order='C')
    multiplicity = np.zeros((len(supercell_fracs), len(primitive_fracs)),
                            dtype='intc', order='C')
    import phonopy._phonopy as phonoc
    phonoc.gsv_set_smallest_vectors(
        shortest_vectors,
        multiplicity,
        supercell_fracs,
        primitive_fracs,
        lattice_points,
        np.array(reduced_bases.T, dtype='double', order='C'),
        np.array(trans_mat_inv.T, dtype='intc', order='C'),
        symprec)

    # # For every atom in the supercell and every atom in the primitive cell,
    # # we want 27 images of the vector between them.
    # #
    # # 'None' is used to insert trivial axes to make these arrays broadcast.
    # #
    # # shape: (size_super, size_prim, 27, 3)
    # candidate_fracs = (
    #     supercell_fracs[:, None, None, :]    # shape: (size_super, 1, 1, 3)
    #     - primitive_fracs[None, :, None, :]  # shape: (1, size_prim, 1, 3)
    #     + lattice_points[None, None, :, :]   # shape: (1, 1, 27, 3)
    # )

    # # To compute the lengths, we want cartesian positions.
    # #
    # # Conveniently, calling 'numpy.dot' between a 4D array and a 2D array
    # # does vector-matrix multiplication on each row vector in the last axis
    # # of the 4D array.
    # #
    # # shape: (size_super, size_prim, 27)
    # lengths = np.array(np.sqrt(
    #     np.sum(np.dot(candidate_fracs, reduced_bases)**2, axis=-1)),
    #                    dtype='double', order='C')

    # # Create the output, initially consisting of all candidate vectors scaled
    # # by the primitive cell.
    # #
    # # shape: (size_super, size_prim, 27, 3)
    # candidate_vectors = np.array(np.dot(candidate_fracs, trans_mat_inv),
    #                              dtype='double', order='C')

    # # The last final bits are done in C.
    # #
    # # We will gather the shortest ones from each list of 27 vectors.
    # shortest_vectors = np.zeros_like(candidate_vectors,
    #                                  dtype='double', order='C')
    # multiplicity = np.zeros(shortest_vectors.shape[:2], dtype='intc',
    #                         order='C')

    # import phonopy._phonopy as phonoc
    # phonoc.gsv_copy_smallest_vectors(shortest_vectors,
    #                                  multiplicity,
    #                                  candidate_vectors,
    #                                  lengths,
    #                                  symprec)

    return shortest_vectors, multiplicity
Find shortest atomic pair vectors

    Note
    ----
    Shortest vectors from an atom in primitive cell to an atom in
    supercell in the fractional coordinates of primitive cell. If an
    atom in supercell is on the border centered at an atom in
    primitive and there are multiple vectors that have the same
    distance (up to tolerance) and different directions, several
    shortest vectors are stored.
    In fact, this method is not limited to search shortest vectors between
    sueprcell atoms and primitive cell atoms, but can be used to measure
    shortest vectors between atoms in periodic supercell lattice frame.

    Parameters
    ----------
    supercell_bases : ndarray
        Supercell basis vectors as row vectors, (a, b, c)^T. Must be order='C'.
        dtype='double'
        shape=(3, 3)
    supercell_pos : array_like
        Atomic positions in fractional coordinates of supercell.
        dtype='double'
        shape=(size_super, 3)
    primitive_pos : array_like
        Atomic positions in fractional coordinates of supercell. Note that not
        in fractional coodinates of primitive cell.
        dtype='double'
        shape=(size_prim, 3)
    symprec : float, optional, default=1e-5
        Tolerance to find equal distances of vectors

    Returns
    -------
    shortest_vectors : ndarray
        Shortest vectors in supercell coordinates. The 27 in shape is the
        possible maximum number of elements.
        dtype='double'
        shape=(size_super, size_prim, 27, 3)
    multiplicities : ndarray
        Number of equidistance shortest vectors
        dtype='intc'
        shape=(size_super, size_prim)