1. A computer-implemented method of simulating behaviour of a thermodynamic system over time, the thermodynamic system having potential energy that can be split into more quickly varying parts and more slowly varying parts and having a state described by collective vectors of position and momentum at any given time, the method comprising: a momentum refreshment process and a conservative dynamics process performed by a processor, wherein the momentum refreshment process comprises carrying out an operation of mixing the collective momentum vector with a noise vector and carrying out an acceptance and rejection operation; the conservative dynamics process comprises applying a multiple time stepping conservative dynamics operation to a current state, in which operation calculations for forces corresponding to more slowly varying energy parts in the thermodynamic system undergo an averaging procedure and are carried out at a larger time step than calculations for forces corresponding to more quickly varying energy parts; and carrying out an acceptance and rejection operation; wherein the acceptance and rejection operations are based on an approximation of the system energy expressed using shadow Hamiltonians and comprise accepting a current state or returning a replacement state, wherein the states in the method are denoted by Ω i =(Y i T ,t i ) T , i=0,. .. , I, where I is a given integer, Y i =(X i T ,1,P i T ,b i ) T , X i is a collective vector of atomic positions, P i is a collective vector of atomic momenta, b i is a scalar, and t i is time and wherein the mixing operation comprises: given a current state, mixing its collective atomic momentum vector P i with an independent and identically distributed normal noise vector Ξ i of dimension 3N, so that 
 {tilde over (P)} i =cos(φ) P i +sin(φ)Ξ i , 
 {tilde over (Ξ)} i =cos(φ)Ξ i −sin(φ) P i , where i is a given integer, N is the number of particles in the system, 0<φ≦π/2 is a given angle, Ξ i ˜N[0,βM −1 ], N[0,βM −1 ] denotes the (3N)-dimensional normal distribution with zero mean and covariance matrix βM −1 , M is the diagonal mass matrix of the molecular system, and β=1/k B T is the inverse temperature, proposed vector being denoted by {tilde over (Ω)} i =({tilde over (Y)} i T ,t i ) T , {tilde over (Y)} i =(X i T ,1,{tilde over (P)} i T ,b i ) T ; and preferably wherein the subsequent acceptance/rejection operation comprises obtaining the resulting state Ω i ; through a Metropolis accept/reject criterion: ⁢ Ω _ i = { Ω ~ i with ⁢ ⁢ probability ⁢ ⁢ min ( 1 , exp ⁡ ( - βΔ ⁢ ⁢ H Δ ⁢ ⁢ t e ) Ω i otherwise , ⁢ ⁢ ⁢ with Δ ⁢ ⁢ H Δ ⁢ ⁢ t e := H Δ ⁢ ⁢ t ⁡ ( X i , P ~ i ) + 1 2 ⁢ Ξ ~ i T ⁢ M - 1 ⁢ Ξ ~ i - H Δ ⁢ ⁢ t ⁡ ( X i , P i ) - 1 2 ⁢ Ξ i T ⁢ M - 1 ⁢ Ξ i , in which H Δt is the shadow Hamiltonian and ΔH Δt e is the change in shadow Hamiltonian due to the mixing operation.