A new method for obtaining molecular vibrational eigenstates using an efficient basis set made up of semiclassical eigenstates is presented. Basis functions are constructed from a ‘‘primitive’’ basis of Gaussian wave packets distributed uniformly on the phase space manifold defined by a single quasiperiodic classical trajectory (an invariant N‐torus). A uniform distribution is constructed by mapping a grid of points in the Hamilton–Jacobi angle variables, which parametrize the surface of the N‐torus, onto phase space by means of a careful Fourier analysis of the classical dynamics. These primitive Gaussians are contracted to form the semiclassical eigenstates via Fourier transform in a manner similar to that introduced by De Leon and Heller [J. Chem. Phys. 8 1, 5957 (1984)]. Since the semiclassical eigenstates represent an extremely good approximation of the quantum eigenstates, small matrix diagonalizations are sufficient to obtain eigenvalues ‘‘converged’’ to 4–5 significant figures. Such small diagonalizations need not include the ground vibrational state and thus can be used to find accurate eigenstates in select regions of the eigenvalue spectrum. Results for several multidimensional model Hamiltonians are presented.