code: def get_bond_symmetry(site_symmetry,
                        lattice,
                        positions,
                        atom_center,
                        atom_disp,
                        symprec=1e-5):
      
      bond_sym = []
      pos = positions
      for rot in site_symmetry:
          rot_pos = (np.dot(pos[atom_disp] - pos[atom_center], rot.T) +
                     pos[atom_center])
          diff = pos[atom_disp] - rot_pos
          diff -= np.rint(diff)
          dist = np.linalg.norm(np.dot(lattice, diff))
          if dist < symprec:
              bond_sym.append(rot)
  
      return np.array(bond_sym)
docs: Bond symmetry is the symmetry operations that keep the symmetry
      of the cell containing two fixed atoms.