Coherent states in imaginary time are used to represent the Boltzmann operator in terms of a classical chain of N coupled beads. The resulting representation is similar to the path integral molecular dynamics (PIMD) approach, except for additional nearest-neighbor coupling terms in the effective Hamiltonian. Each time step of the molecular dynamics evolution of the effective classical system provides appropriately thermalized Gaussians that can be used as distribution functions to sample initial conditions for real-time semiclassical initial value representation (IVR) methods. The approach is illustrated with numerical calculations of thermal expectation values and canonical time correlation functions for several one-dimensional quantum systems.