Molecular quantum dynamics simulations are essential for understanding many fundamental phenomena in physics and chemistry. They often require solving the time-dependent Schrödinger equation for molecular nuclei, which is challenging even for medium-sized systems. Therefore, for the majority of molecules, it is necessary to use approximate methods to evaluate the nuclear quantum dynamics. In this thesis, we comprehensively investigate the semiclassical variational thawed Gaussian wavepacket dynamics. This method is the most accurate of the single-trajectory Gaussian-based techniques. Moreover, it is symplectic, time-reversible, and norm- and energy-conserving. In one-dimensional double-well potentials, we numerically demonstrate that the variational method may include quantum tunneling approximately. To model anharmonic multidimensional systems, we develop a non-separable coupled Morse potential for which the expectation value in a Gaussian wavepacket can be evaluated analytically. To show that the variational method, unlike the exact quantum grid-based techniques, can simulate high-dimensional systems, we do most of our analysis on a twenty-dimensional coupled Morse system. The equations of motion in the variational method require expectation values of the potential energy, its gradient, and its Hessian, which can be prohibitively expensive to be evaluated for real molecular systems. One way to overcome this problem is to Taylor expand the potential energy around the center of the Gaussian wavepacket. Using the local harmonic approximation of the potential in the variational method yields Heller's thawed Gaussian approximation, which conserves neither the symplectic structure nor the effective energy. Therefore, we employ the local cubic approximation of the potential. The local cubic variational thawed Gaussian wavepacket dynamics is symplectic, conserves the effective energy, and improves the accuracy over the local harmonic method. Although they are more accurate than Heller's thawed Gaussian approximation, both the fully variational and local cubic variational methods are also much more expensive. In order to reduce the computational cost, the high-order geometric integrators are extensively investigated in this thesis. We implement the second-order geometric integrator developed by Faou and Lubich and apply various composition methods to obtain high-order geometric integrators. The high-order integrators can be significantly more efficient than the second-order algorithm, while at the same time being much more accurate. In addition, the geometric integrators preserve the geometric properties of the variational methods, such as norm conservation, symplecticity, time reversibility, and, for sufficiently small time steps, energy or effective energy conservation. To include non-zero temperature effects in the simulation of vibronic spectra, we use the concept of thermo-field dynamics. By computing the vibronic spectra of a harmonic potential at different temperatures, we demonstrate the validity of the thermo-field thawed Gaussian wavepacket dynamics. Furthermore, for the calculation of vibronic spectra of systems with low anharmonicity at finite temperatures, we show the superior accuracy of the fully variational and local cubic variational thermo-field thawed Gaussian wavepacket dynamics over the local harmonic thermo-field method.