We present a novel approach to numerical modelling of the crystalline solid as a heat bath. The approach allows bringing together a small and a large crystalline domain, and model accurately the resultant interface, using harmonic assumptions for the larger domain, which is excluded from the explicit model and viewed only as a hypothetic heat bath. Such an interface is non-reflective for the elastic waves, as well as providing thermostatting conditions for the small domain. The small domain can be modelled with a standard molecular dynamics approach, and its interior may accommodate arbitrary non-linearities. The formulation involves a normal decomposition for the random thermal motion term R(t) in the generalized Langevin equation for solid–solid interfaces. Heat bath temperature serves as a parameter for the distribution of the normal mode amplitudes found from the Gibbs canonical distribution for the phonon gas. Spectral properties of the normal modes (polarization vectors and normal frequencies) are derived from the interatomic potential. Approach results in a physically motivated random force term R(t) derived consistently to represent the correlated thermal motion of lattice atoms. We describe the method in detail, and demonstrate applications to one- and two-dimensional lattice structures. Copyright © 2006 John Wiley & Sons, Ltd.