Calculating DOSs can be tricky business. It sounds like for your purposes, you don't have to be super accurate. So, I suggest using a quick and dirty way of calculating the DOS:
- Choose a bunch (at least a million, but feel free to go wild) $\mathbf{q}$ points at random.
- For each one, find $\omega\left(\mathbf{q}\right)$
- Make a histogram of the $\omega$ values you calculated.
That histogram is (proportional to) your DOS.
How does this work? The DOS basically tells you "how many" states have a given energy (or really, how many states are in a small energy window). So, if there are more states at a given energy, then a randomly chosen state is also more likely to have that energy. In short, the DOS is proportional to the probability distribution of the energy for a randomly chosen $\mathbf{q}$. Making a histogram in the way described above is a way to calculate the probability distribution. You'll have to figure out the proportionality constant (or fudge it, I think the Debye model should produce the same results at small $\mathbf{q}$). All this follows from an alternate (but equivilent) definition of the DOS: $g\left(\omega\right) = \int \frac{d\mathbf{q}}{\left(2 \pi \right)^d} \delta\left(\omega - \omega\left(\mathbf{q}\right)\right)$.
Another option is to convert that surface integral into a volume integral using the divergence theorem. This can be nice because finding the surface itself is tricky, but finding if a given $\mathbf{q}$ is inside the volume is easy. Just check that $\omega\left(\mathbf{q}\right) < const$. Computers are much better at volume integrals than surface integrals. This method is tricky, altho it's not totally crazy in your case since you have an analytic form for your dispersion relation. You could take the divergence analytically then do the integral numerically. However, I'd go with the first method unless you need a really accurate answer. There are many more advanced methods out there too, but they're probably overkill for your situation.
(In answer to your question, the surface you'd want to integrate over is an equpotential, so you're right that it's defined by $\omega\left(\mathbf{q}\right)= const$. If you want to do the surface integral, you'd first need to find the surface and then integrate over it. I'd avoid that at all costs if I were you.)
Sorry for the late response but hopefully this can be useful for someone else!
You can reduce the noise using an elliptic integral.
$$D(\varepsilon)=\frac{1}{4 \pi ^3 t}\int_{-\pi }^{\pi }d\phi K\left( \sqrt{1-\left(\frac{\varepsilon +2 t \cos\phi }{4 t}\right)^2}\right)$$
Where K is the Complete Elliptic Integral of the First Kind: http://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html.
It's not trivial to get here. And even from this expression the integral needs to be done numerically with care (it has singularities for many values of $\varepsilon$) but it should give better results. Running for six seconds in Mathematica gives me (with $t=\frac{1}{2\sqrt{3}}$):
Best Answer
I will present the result for the one-dimensional case only. It is readily generalized to your three-dimensional situation. Knowing the phonon dispersion $\omega(q)$, the DOS can be defined as
$$D(\omega) = \frac{1}{2\pi}\int\mathrm{d}q\,\delta(\omega-\omega(q)).$$
Using the formal identity
$$\delta(g(x)) = \sum_{i}\frac{\delta(x-x_i)}{|g'(x_i)|}$$
where $x_i$ are the zeros of $g$ and $g'$ is its derivative, one can write:
$$D(\omega) = \frac{1}{2\pi}\int\mathrm dq \frac{1}{\left|\frac{\mathrm d\omega}{\mathrm dq}\right|}\sum_{\substack{i\\\omega(q_i)=\omega}}\delta(q-q_i) = \frac{1}{2\pi}\sum_{\substack{i\\\omega(q_i)=\omega}}\left|\frac{\mathrm d\omega}{\mathrm dq}(q_i)\right|^{-1}$$
The sum is over all the points where the dispersion reaches the given frequency $\omega$. The flater the dispersion at those points, the more states contribute!