Is it possible to compute the colour of *any* molecule? Above are some spectral emission lines for hydrogen, helium and neon. 

Using the Schroedinger equation, it's possible to derive the colours that hydrogen will emit when light is fired at it. A little more complicated (although beyond my rudimentary knowledge of quantum mechanics), the same can be done for helium, and so on.

Of course as we progress to molecules, even the relatively small $C_5H_{11}Br$, it's (I'm told) very hard to predict the colours emitted from first principles, and even then incorrect assumptions must be made (like the Born-Oppenheimer approximation), so upon progressing to larger molecules, we find that these simplifications diverge from experiment.

My question is this: with the advent of computers, are we able to (in practice, not theoretically in a few decades) accurately compute the colours given off by a molecule  by numerically solving the laws of quantum mechanics for complicated systems? Much like nonlinear differential equations can now be 'solved' with computers without any simplifying assumptions (like $\sin(x) \approx x$), is it in humanity's capacity today to approximate colours from first principles, so that computing time is the only restraint? Note the difference between an incorrect model of reality (Born-Oppenheimer), and one that converges to the correct one (even if slowly).
It is possible to compute electronic spectra with decent accuracy. Time dependent density functional theory (TDDFT) permits fairly rapid calculation of UV/vis spectra at a quantum mechanical level of theory. Whilst DFT is not truly a first principles method, it can furnish electronic spectra accurate enough to be useful, depending on care taken in selecting appropriate functionals, accounting for solvatochromic effects, etc. Obviously, DFT uses the Born-Oppenheimer approximation, but the beauty of the approximation is that it's actually very close to reality.

See, for instance, this paper by Jacquemin et al., where they calculate $\lambda_{max}$ for a wide range of dyes using TDDFT and get a MAE of 22 nm using CAM-B3LYP (a long-range hybrid functional), reducing that to 10 nm with an empirical calibration curve. The paper only really cares about accurately reproducing the location of $\lambda_{max}$, however this paper by Quartarolo and Russo computes complete spectra with fairly reasonable accuracy. Not perfect, but the general shape of the spectra are reproduced adequately. Part of the divergence between the experimental and computed spectra probably comes down to higher order broadening effects not reflected in the simple Gaussian convolution of the computed data, but I'm not a spectroscopist.

Additional: as for computing the colour of any molecule, well, that does depend heavily on the size of your molecule and any mitigating factors such as symmetry elements. Most QM techniques scale fairly viciously with the number of electrons in your molecule, and whilst you can always reduce the size of your basis set and heap on less and less conservative approximations, accuracy will suffer.