We compute a light front wave function for heavy vector mesons based on long-distance matrix elements constrained by decay width analyses in the nonrelativistic QCD framework. Our approach provides a systematic expansion of the wave function in quark velocity. The first relativistic correction included in our calculation is found to be significant and crucial for a good description of the HERA exclusive J / ψ production data. When looking at cross section ratios between nuclear and proton targets, the wave function dependence does not cancel out exactly. In particular the fully nonrelativistic limit is found not to be a reliable approximation even in this ratio. The important role of the Melosh rotation to express the rest frame wave function on the light front is illustrated.