Implemented in Stan language.

Transform unormalized posterior into a Hamiltonian/Energy by $H=U(\theta)+K(p)$ $=-\log{[p(X|\theta) p(\theta)]}+\frac{p^T p}{2}$.

For each iteration:

Give a random momentum to each of an ensemble of particles, then making it evolve according to Hamiltonian dynamics for some period of time, and then accept new position according to an acceptance rule, like that of the Metropolis-Hastings algorithm (with a symmetric Proposal distribution), which is that we accept if $H_{\text{final}} > H_{\text{initial}}$, and otherwise we use the the ratio of Boltzmann distribution as the probability of accepting.

As explained here, the symmetric proposal distribution is ensured by using a time-reversible (approximate) Hamiltonian dynamics.

As said here, if the Hamiltonian dynamics simulation is perfect, the acceptance rule always accept because $H$ is conserved. Therefore, we use approximate Hamiltonian dynamics (like Leapfrog algorithm) to explore more.

Note that the symmetry of the proposal distribution means that it is not the Boltzmann distribution (as it would be if we simulated Langevin dynamics). In fact, if using a deterministic algo like the leapfrog algo, our proposal is deterministic. However, perhaps it still has, over many iterations, properties like the Boltzmann dist (like a derandomized Langevin dynamics), or at least it "tends" to go downhill, and that's why it's useful?

Hamiltonian Monte Carlo explained

McKay's book chapter

another book chapter

Leapfrog algorithm is a common way of solving the Hamiltonian dynamics (e.g. http://diffsharp.github.io/DiffSharp/examples-hamiltonianmontecarlo.html)

Scalable Bayesian Inference with Hamiltonian Monte Carlo

Efficient Bayesian inference with Hamiltonian Monte Carlo -- Michael Betancourt (Part 1)