For generating an isobaric–isothermal ensemble in molecular dynamics simulations of atomic systems a correct ensemble distribution can be obtained by the approach of Martyna, Tobias and Klein [J. Chem. Phys. 101, 4177–4189 (1994)]. The constituting equations of the latter approach have been also generalised to molecular systems [M. E. Tuckerman, Statistical Mechanics: Theory and Molecular Simulation, Oxford Graduate Texts (2010)] using the molecular virial instead of the atomic one. An isothermal–isobaric method for systems with holonomic molecular constraints has been also introduced in the past [G. Kalibaeva, M. Ferrario and G. Ciccotti, Mol. Phys. 101 765–778 (2003)] but the factorisation was not worked out completely, resulting in a non-reversible first-order algorithm. Here, we propose the correct factorisation and implement a reversible integrator based on Martyna, Tobias and Klein equations for a system of molecules subjected to holonomic constraints. We add to the algorithm the Suzuki–Yoshida treatment that improves the energy conservation or permits to go to larger time steps. Finally, we test our algorithm by applying it to the dynamics of a Ortotherphenyl model.