  By employing the Dirac-Frenkel time-dependent variational principle, we study
the dynamical properties of the Holstein molecular crystal model with diagonal
and off-diagonal exciton-phonon coupling. A linear combination of the Davydov
D$_1$ (D$_2$) Anstaz, referred to as the multi-D$_1$ Ansatz (multi-D$_2$
Ansatz), is used as the trial state with enhanced accuracy but without
sacrificing efficiency. The time evolution of the exciton probability is found
to be in perfect agreement with that of the hierarchy equations of motion,
demonstrating the promise the multiple Davydov trial states hold as an
efficient, robust description of dynamics of complex quantum systems. In
addition to the linear absorption spectra computed for both diagonal and
off-diagonal cases, for the first time, $2$D spectra have been calculated for
systems with off-diagonal exciton-phonon coupling by employing the multiple
$D_2$ Ansatz to compute the nonlinear response function, testifying to the
great potential of the multiple $D_2$ Ansatz for fast, accurate implementation
of multidimensional spectroscopy. It is found that the signal exhibits a single
peak for weak off-diagonal coupling, while a vibronic multi-peak structure
appears for strong off-diagonal coupling.
