2H NMR spin-lattice relaxation was used to probe small-amplitude torsional vibrations (librations) of the organic terephthalate linkers in metal-organic framework (MOF) UiO-66(Zr) saturated with benzene molecules. In UiO-66 (Zr) the mobile phenylene fragments exhibit a complex rotational dynamics of the phenylene rings with fast librations and much slower π-flips around the C2 symmetry axis. We show that due to the intrinsic broad distribution of the π-flips rate, the relaxation process for the deuterium in the C-D group of phenylene fragment is multiexponential. Two main modes of T1 relaxation are clearly detected, corresponding to the fast T1 fast and the slow T1 slow relaxation. Based on the experimental observation of two-exponential relaxation, a computational model for this T1 relaxation behavior capable to reproduce the peculiarities of the MOF linkers dynamics was built. Computational analysis allows to establish that the librational motion affects mostly the T1 slow, while T1 fast remains unaffected by this motion. Simulation of the T1 slow dependence on the libration rate klib shows that in the range of the librational frequencies of 106-109 Hz the T1 slow is not sensitive to the klib variation, and therefore a precise correspondence between T1 slow and klib cannot be established. T1 slow exhibits a specific "peak-like-shape" dependence of klib in the range of 109-1012 Hz. In this range of libration frequencies an unambiguous relation between T1 slow and klib exists only in a very narrow frequency window of 0.1 × 1010-5 × 1010 Hz. The best conditions to characterize the librational motion by means of T1 relaxation analysis are met when the flipping motion is almost frozen (kflip < 103 Hz) because T1 slow becomes extremely sensitive to the variation of klib.