We studied the accuracy, robustness, and self-consistency of pixel-domain simulations of the gravitational lensing effect on the primordial cosmic microwave background (CMB) anisotropies due to the large-scale structure of the Universe. In particular, we investigated the dependence of the precision of the results precision on some crucial parameters of these techniques and propose a semi-analytic framework to determine their values so that the required precision is a priori assured and the numerical workload simultaneously optimized. Our focus was on the B-mode signal, but we also discuss other CMB observables, such as the total intensity, T, and E-mode polarization, emphasizing differences and similarities between all these cases. Our semi-analytic considerations are backed up by extensive numerical results. Those are obtained using a code, nicknamed lenS2HAT – for lensing using scalable spherical harmonic transforms (S2HAT) – which we have developed in the course of this work. The code implements a version of the previously described pixel-domain approach and permits performing the simulations at very high resolutions and data volumes, thanks to its efficient parallelization provided by the S2HAT library – a parallel library for calculating of the spherical harmonic transforms. The code is made publicly available.