A novel and efficient pseudospectral method for performing fully coupled six-dimensional bound state dynamics calculations is presented, including overall rotational effects. A Lanczos based iterative diagonalization scheme produces the energy levels in increasing energies. This scheme, which requires repetitively acting the Hamiltonian operator on a vector, circumvents the problem of constructing the full matrix. This permits the use of ultralarge molecular basis sets (up to over one million states for a given symmetry) in order to fully converge the calculations. The Lanczos scheme was conducted in a symmetry adapted spectral representation, containing Wigner functions attached to each monomer. The Hamiltonian operator has been split into different terms, each corresponding to an associated diagonal or nearly diagonal representation. The potential term is evaluated by a pseudospectral scheme of Gaussian accuracy, which guarantees the variational principle. Spectroscopic properties are computed with this method for four of the most widely used water dimer potentials, and compared against recent terahertz laser spectroscopy results. Comparisons are also made with results from other dynamics methods, including quantum Monte Carlo (QMC) and reversed adiabatic approximation calculations. None of the potential surfaces produces an acceptable agreement with experiments. While QMC methods yield good results for ground (nodeless) states, they are highly inaccurate for excited states.