The far infrared spectra of three deuterated isotopologues of ethylene glycol, CH₂OD–CH₂OD, CH₂OD–CH₂OH and CH₂OH–CH₂OD, where the latter two species differ in their intramolecular hydrogen-bonding arrangement, are studied using highly correlated ab initio methods, vibrational second order perturbation theory and a variational procedure of reduced dimensionality. New subroutines suitable for the study of large systems with more than two interacting large amplitude motions were implemented and applied. The molecular symmetry of ethylene glycol decays by the formation of weak intramolecular bonds producing very asymmetrical stable structures. Three internal rotations contribute to the formation of a very anisotropic potential energy surface and to the puzzling distribution of the rovibrational energy levels. The ground vibrational state rotational constants and the centrifugal distortion constants (S-reduction, I^r representation) corresponding to the aGg′ (G1) and gGg′ (G2) forms are provided for the studied isotopologues. The low-lying vibrational levels up to 550 cm⁻¹ are obtained variationally for J = 0. Two series of sublevels of the ground vibrational state are obtained: eight components localized in G1 lying between 0.0 and 0.3 cm⁻¹ and eight sublevels localized in G2 lying between 138.1 and 138.4 cm⁻¹. The gap between both sets is lower in CH₂OD–CH₂OD and more dispersed in the monodeuterated variety.