An improved constraint approach reduces the energy drift rate to acceptable levels. In an embodiment of this approach, massively parallel constrained velocity Verlet NVE (constant particle number, constant volume, constant energy) MD simulations can be run using single precision arithmetic with very low energy drift (e.g., ˜1 Kelvin per microsecond simulated time) using large timesteps (e.g., 2.5 fs) for typical systems and MD force fields.