Andersen thermostat in MD with rigid bond constraints I have been writing a custom MD code for some time. I have been using the Andersen thermostat, mainly because its ease of use, and the fact that I am not concerned with gathering dynamical data, but rather free energies. I have implemented the SHAKE algorithm found in LAMMPS. However, the thermostat seems to be undercooling now, when constraints are present. For instance, if I make the set point temperature 350 K, my system will equilibrate at around 330 K.

I am calculating the number of degrees of freedom for the temperature calculation as $3N-3-$number of bond constraints. Has anyone else had a similar experience with the Andersen thermostat? Would a Langevin thermostat work better here? Why is the Andersen thermostat not playing well with the SHAKE algorithm?
Hünenberger’s Thermostat algorithms for molecular dynamics simulations (PDF available here) is a generally concise and accurate reference on the topic. On stochastical thermostats such as Langevin, it does not use the number of degrees of freedom (i.e. reduced by the number of constraints) but actually uses the number of atoms.

One would have to follow the proper derivation of the thermostat's target ensemble (from the original Andersen paper) to make sure that is correct, but it seems likely to fix your problem of having too low kinetic energy.