The exploration of electronically excited states and the study of diverse photochemical and photophysical processes are the main goals of molecular electronic spectroscopy. Exact quantum-mechanical simulation of such experiments is, however, beyond current computational capabilities for the majority of systems. Semiclassical methods, which only require local information about the potential energy surface, provide a practical alternative, because they allow on-the-fly calculations using high-level ab initio electronic structure methods. In this thesis, we explore the on-the-fly ab initio thawed Gaussian approximation one of the simplest yet most practical semiclassical approximation. Despite its simplicity, the thawed Gaussian approximation requires a rather expensive evaluation of the Hessian along the trajectory. To accelerate this method, we symmetrically compose the standard second-order geometric integrator, a semiclassical generalization of the Störmer-Verlet algorithm, and obtain arbitrary even-order geometric integrators. Using a 20-dimensional coupled Morse potential, we show that the high-order integrators are more efficient than the standard second-order method while preserving the norm and time-reversibility. To illustrate that the efficiency gains and preservation of geometric properties are not limited to model systems with analytical potentials, we perform on-the-fly ab initio dynamics of ammonia on its first excited-state surface. Although the thawed Gaussian approximation is generally regarded as only a crude approximation to the exact solution of the Schrödinger equation, we show that it captures the isotope effects in the absorption spectra of ammonia isotopologues. In contrast, the popular global harmonic models fail to properly account for anharmonicity, which is essential for reproducing experimental results. Further analysis reveals that the single progression seen in the experiment can be explained by the commensurate evolution of normal modes and the missing mode effect. The choice of normal-mode coordinates, needed for propagating the thawed Gaussian wavepacket in ab initio applications, is not unique. This choice influences the region of the phase space where the rovibrational coupling is minimized. We show that whereas the choice of normal modes has a negligible effect on the electronic spectra, it has important implications for the dynamical interpretation of these spectra. The normal-mode coordinates can be associated with the global harmonic model that is constructed about the reference configuration. We show that while such harmonic models can give entirely wrong spectra, the same normal-mode coordinates used in the thawed Gaussian approximation can yield satisfactory results.