We present a method based on centroid molecular dynamics (CMD) to calculate nonlinear quantum force correlation functions important in the golden rule approach for studying vibrational energy relaxation (VER) in condensed phases. We consider a model of a diatomic molecule in a two-dimensional neon liquid and also a diatomic coupled to a small Helium cluster. The predictions of the theory for the neon bath are compared and found in close agreement with available theories for VER based on the Egelstaff correction factor and Feynman-Kleinert variational theory. For the Helium cluster, the force spectrum obtained from CMD is found to be in slightly better agreement with the exact result than a method based on a cumulant approach. The results support the use of CMD in condensed phase studies of VER when quantum effects are important.