Liouville operators as follows: where Defining where α = A, B, C , one possible reversible measure-preserving integration algorithm for the BK thermostat is then Using the so-called direct translation technique [36] we can expand the above symmetric break-up of the Liouville operator into a pseudo-code form, ready to be implemented on the computer: L BK A ≡ ˜ L BK 1 + ˜ L BK 4 = F ( q ) ∂ ∂p + F$_{p}$$_{ξ}$ ∂ ∂p$_{ξ}$ (24a) L BK B ≡ ˜ L BK 2 + ˜ L BK 3 = p m ∂ ∂q + F$_{p}$$_{ζ}$ ∂ ∂p$_{ζ}$ (24b) L BK C ≡ ˜ L BK 5 + ˜ L BK 6 = − p$_{ζ}$ m$_{ζ}$ p ∂ ∂p − p$_{ξ}$ m$_{ξ}$ q ∂ ∂q + p$_{ζ}$ m$_{ζ}$ ∂ ∂ζ + p$_{ξ}$ m$_{ξ}$ ∂ ∂ξ (24c) F ( q ) = − ∂V /∂q (25a) U BK α ( τ ) = exp [ τ ˜ L BK α ] , (26) F$_{p}$$_{ζ}$ = p 2 m − k$_{B}$T . (25c) F$_{p}$$_{ξ}$ = q ∂V ∂q − k$_{B}$T (25b) U ( τ ) BK = U BK B ( τ 4 ) U BK C ( τ 2 ) U BK B ( τ 4 ) × U BK A ( τ ) × U BK B ( τ 4 ) U BK C ( τ 2 ) U BK B ( τ 4 ) . (27) • q → q + τ 4 p m p$_{ζ}$ → p$_{ζ}$ + τ $_{4}$F$_{p}$$_{ζ}$    : U BK B ( τ 4 ) 10