A method to compute solubilities for molecular systems using atomistic simulations, based on an extension of the Einstein crystal method, has recently been presented [Li et al., J. Chem. Phys. 146, 214110 (2017)]. This methodology is particularly appealing to compute solubilities in cases of practical importance including, but not limited to, solutions where the solute is sparingly soluble and molecules of importance for the pharmaceutical industry, which are often characterized by strong polar interactions and slow relaxation time scales. The mathematical derivation of this methodology hinges on a factorization of the partition function which is not necessarily applicable in the case of a system subject to holonomic molecular constraints. We show here that, although the mathematical procedure to derive it is slightly different, essentially the same mathematical relation for calculating the solubility can be safely applied for computing the solubility of systems subject to constraints, which are the majority of the systems used for practical molecular simulations.