A modification of the constrained equations of motion of Kalibaeva et al. [Mol. Phys. 101, 765 (2003)] in the NPH and NPT ensembles is presented. The modified equations of motion are discretized using central-difference techniques, and the derived integrators are time reversible and conserve the invariant phase space measure. The constraint algorithm builds on the work of Toxvaerd et al. [J. Chem. Phys. 131, 064102 (2009)] in the NVE and NVT ensembles: it thus conserves the holonomic bond constraints at the finite machine precision level in the NPH and NPT ensembles. The algorithms were tested on a system of n=320 ortho-terphenyl molecules, arriving at the target temperature and pressure in a low and high pressure state. Isobaric heat capacities in the NPH and NPT ensembles were calculated for comparison using the fluctuation formulas as well as the thermodynamic definition. The heat capacities agree within the estimated uncertainties.