We derive a c-number generalized Langevin equation (GLE) describing the evolution of the expectation values ⟨xi⟩t of the atomic position operators xi of an open system. The latter is coupled linearly to a harmonic bath kept at a fixed temperature. The equations of motion contain a non-Markovian friction term with the classical kernel [L. Kantorovich, Phys. Rev. B 78, 094304 (2008)] and a zero mean non-Gaussian random force with correlation functions that depend on the initial preparation of the open system. We used a density operator formalism without assuming that initially the combined system was decoupled. The only approximation made in deriving quantum GLE consists of assuming that the Hamiltonian of the open system at time t can be expanded up to the second order with respect to operators of atomic displacements ui=xi−⟨xi⟩t in the open system around their exact atomic positions ⟨xi⟩t (the “harmonization” approximation). The noise is introduced to ensure that sampling many quantum GLE trajectories yields exactly the average one. An explicit expression for the pair correlation function of the noise, consistent with the classical limit, is also proposed. Unlike the usually considered quantum operator GLE, the proposed c-number quantum GLE can be used in direct molecular dynamic simulations of open systems under general equilibrium or nonequilibrium conditions.