In classical molecular simulations chemical bonds and bond angles have been modeled either as rigid constraints, or as nearly harmonic oscillators. However, neither model is a good description of a chemical bond, which is a quantum oscillator that in most cases occupies the ground state only. A quantum oscillator in the ground state can be represented more faithfully by a flexible constraint. This means that the constraint length adapts itself, in time, to the environment, such that the rotational and potential forces on the constraint cancel out. An accurate algorithm for flexible constraints is presented in this work and applied to study liquid water with the flexible and the polarizable "mobile charge densities in harmonic oscillators" model. The iterations for the flexible constraints are done simultaneously with those for the electronic polarization, resulting in negligible additional computational costs. A comparison with fully flexible and rigidly constrained simulations shows little effect on structure and energetics of the liquid, while the dynamics is somewhat faster with flexible constraints.