How to calculate the quantum density of states?
If you really do want to calculate the density of states then you can try direct counting to get the total number of states N and then calculate the density of states $\rho =[N(E+dE)-N(E)]/dE$. This, of course, will only work if you have a few vibrational normal modes as all the combinations of quanta have to be included, e.g. in water with 3 modes, quanta (in the order of increasing energy) are: $(0,0,0), (0,0,1),(0,0,2),(0,1,0), (1,0,0),(0,0,3),(0,1,1)$, etc.  But for a molecule with just a few more atoms, it is quite probable you will get $10^6$ vibrational states /wavenumber at an excess energy of, say,  $5000 \pu{~cm^{-1}}$.

In the field of unimolecular reactions (RRKM models etc.) there has been much effort in calculating total and density of states. There is the Beyer-Swinehart algorithm; Beyer, T., &amp; Swinehart, D. F. 1973, Commun. Assoc. Comput. Machinery, 16, 379 and Stein S &amp; Rabinovitch B. J. Chem. Phys. v58, 2438, 1973. Whitten and Rabinovitch have also modified the classical sum of states to give good agreement with the direct counting (quantum) method. (J. Chem. Phys. 41, 1883, 1964, and 48, 1427 1968). Their modification $^{[1]}$ gives the total sum of states as:

$$ N(E)=\frac{(E+\alpha E_{zpe})^s}{s!\prod_{i=1}^s h\nu_i}$$

where s is the number of (harmonic) oscillators, $E_{zpe}$ is the zero point energy for the i th oscillator, 

$$\alpha = 1-\beta w(E')$$ 

$$\beta = \frac{s-1}{s}\frac{&lt;\nu^2&gt;}{&lt;\nu&gt;^2}$$ 

and w depends on the range of $E'=E/E_{zpe}$. When $0.1 \lt E' \lt 1; w(E')=(5E'+2.73\sqrt(E')+3.51)^{-1}$, and when $1\lt E' \lt 8; w(E')=\exp(-2.4191(E')^{1/4})$.
The equation in N can be differentiated to find the density of states.

If, however, you want to calculate a rate constant using the Fermi Golden rule it is better not to start with the density of states separated out, but to calculate the Franck-Condon weighted density of states. 

The rate constant has the Golden Rule form 
$$ k= \frac{2\pi}{\hbar}|&lt;av''|H|bv'&gt;|^2\delta(E_{av'}-E_{bv''})$$

and then this is modified depending on what actual situation is to be analysed, e.g. thermally averaged over the Boltzmann distribution if the molecules are in solution.  

In this case the best place to start is probably with the papers 

(a)  Heller, Freed &amp; Gelbart, J. Chem. Phys. 56, 2309 (1972),

(b)  Bixon &amp; Jortner, J. Chem. Phys. 48, 715-726, 1968 and 

(c) S. H. Lin; J. Chem. Phys. 58, 5760 (1973). 

Each deal with non-adiabatic interactions causing radiationless transitions in molecules, i.e. from one vibronic level into a continuum of states, such as singlet to triplet intersystem crossing. There is a summary of different methods in Eyring, Lin &amp; Lin 'Basic Chemical Kinetics' ch 7.

$^{[1]}$ equations given in Steinfeld, Francisco &amp; Hase, 'Chemical Kinetics &amp; Dynamics'. 2nd ed. p333.