Some comments to the setup in your script file:
1). The k-mesh is not set properly. The size of k-mesh should be properly scaled according to the length of the lattice vectors. Since the lattice constant c is roughly 4 times larger than the lattice constant a, the size of k-mesh along the kz direction should be smaller than the one along the kx and ky direction. So, "k_point_sampling=(11, 11, 3)," may be more reasonable than "k_point_sampling=(11, 11, 11)".
2). To reproduce the results in literature properly, one had better choose the same exchange-correlation functional with the one used in literature, use the more accurate basis set such as DZP. In addition, the same structural parameters should be used.