Friction generally occurs in the relative motion or the contact interface with the trend of relative motion, which impedes the relative motion and produces energy loss. Micro-scale friction is different from the macro-scale friction due to surface effects and other factors. It is necessary to study the friction behavior on a nano-scale. First-principles method is an important way to study and understand friction on a nano-scale. Nevertheless, the constructing of nearly a thousand models and the processing of a large number of data are very time consuming. In this paper, we establish a high-throughput computational program based on the first-principles method to study the interfacial friction of two-dimensional materials. The program realizes modeling, submitting computation tasks, multi-task concurrent calculation, data collection and processing, and image rendering of calculation results. All of these are done in batch automatically, which greatly saves researchers’ time. In this work, this program is used to simulate the normal load by changing the distance between layers and calculate the potential energy surface of BN/BN and graphene/graphene bilayer sliding systems at a series of interlayer distances, as well as the interlayer friction forces and friction coefficients. The study finds that with the decrease of the interlayer distance, the averaged friction force at BN/BN interface increases approximately linearly, and the friction coefficient is in a range of 0.11–0.17. The friction force at graphene/graphene interface first increases, then decreases, and increases again. The friction coefficient reaches a minimum value (0.014) under a load of 12 nN, and these results are consistent with the previous results, verifying the reliability of the calculation program. In addition, we investigate the effect of surface hydrogenation and fluorination on the tribological property of the BN bilayer and find that the friction at the fluorinated BN/BN interface decreases, which is attributed to the smaller charge transfer at interface. Although the high-throughput calculation method realizes the automation and high-throughput calculation of tribological property at solid interface, there are still some limitations. Firstly, the effect of interlaminar bending is not considered in the process of interlaminar relative sliding. Secondly, the essence of the calculation result is static friction, rather than dynamic friction. In addition, the method does not consider the influence of temperature.