We present an algorithm for the simulation of nonlinear 2D spectra of molecular systems in the UV-vis spectral region from atomistic molecular dynamics trajectories subject to nonadiabatic relaxation. We combine the nonlinear exciton propagation (NEP) protocol, that relies on a quasiparticle approach with the surface hopping methodology to account for quantum-classical feedback during the dynamics. Phenomena, such as dynamic Stokes shift due to nuclear relaxation, spectral diffusion, and population transfer among electronic states, are thus naturally included and benchmarked on a model of two electronic states coupled to a harmonic coordinate and a classical heatbath. The capabilities of the algorithm are further demonstrated for the bichromophore diphenylmethane that is described in a fully microscopic fashion including all 69 classical nuclear degrees of freedom. We demonstrate that simulated 2D signals are especially sensitive to the applied theoretical approximations (i.e., choice of active space in the CASSCF method) where population dynamics appears comparable.