In this article, we investigate the convergence of quantum dynamics calculations with coupled variationally optimized gaussian product basis functions, describing wavepacket motion on regions of molecular potential energy surfaces calculated on the fly. As a benchmark system, we model the radiationless decay of fulvene from its first electronic excited state through an extended S(1)∕S(0) conical intersection seam and monitor two associated properties: the spatial extent to which the conical intersection seam is sampled and the timescale and stepwise nature of the population transfer. We suggest that the fully variational description reviewed here (direct dynamics-variational multi-configuration gaussian) provides a way to balance accuracy against computational cost for molecules of comparable sizes by choosing the number of coupled gaussian product basis functions and a middle way forward between grid based and trajectory surface hopping approaches to non-adiabatic molecular quantum dynamics calculations.