The surrogate Hamiltonian is a general scheme to simulate the many body quantum dynamics composed of a primary system coupled to a bath. The method has been based on a representative bath Hamiltonian composed of two-level systems that is able to mimic the true system-bath dynamics up to a prespecified time. The original surrogate Hamiltonian method is limited to short time dynamics since the size of the Hilbert space required to obtain convergence grows exponentially with time. By randomly swapping bath modes with a secondary thermal reservoir, the method can simulate quantum dynamics of the primary system from short times to thermal equilibrium. By averaging a small number of realizations converged values of the system observables are obtained avoiding the exponential increase in resources. The method is demonstrated for the equilibration of a molecular oscillator with a thermal bath.