Thermalized Hamiltonian system (using Hamiltonian Monte Carlo)
We now simulate the same systems with thermalization.
In principle, we could simulate thermalization in a completely
Hamiltonian system by introducing one thousand or so additional
small particles that act as a heat bath. But for speed of simulation,
we replace this pure Hamiltonian approach by what we know its
outcome would be: collisions in a Hamiltonian bath lead to
every particle's momentum being repeatedly randomized from
a Gaussian distribution of variance proportional to (m kT).
Now that we have introduced a value for kT it is meaningful
to mention how big the external force on the piston is.
We define kT = 1, and define the external force to be 0.1, so
that a horizontal piston displacement from x=0 to x=20
causes work of 2kT to be done.
In the animation below, the piston is started well over to the left,
and there is just one particle.
7 
 Piston mass = 10, thermalized

The simulation uses the Hamiltonian Monte Carlo method.
For a review of Hamiltonian Monte Carlo, see
chapter 30 of Information Theory,
Inference, and Learning Algorithms.
Hamiltonian Monte Carlo uses the leapfrog method to
reversibly simulate the Hamiltonian dynamics between thermalization events,
and corrects for inaccuracies of the leapfrog method using
the Metropolis method to reject trajectories in which the
Hamiltonian (which should be conserved) increases by too much.
The result is an algorithm that asymptotically samples from the
equilibrium distribution of the system.