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?