Let us have an integral $I$ over the first Brillouin zone (BZ) of a 2D lattice. $$I = \int\int \Omega(k_x,k_y) dk_x dk_y$$ where $\Omega(k_x,k_y)$ is some function (let's say it's Berry curvature). what is the correct numerical approximation of this integral? I thought Riemann sums would work i.e. $$I \approx \sum_{k_x = k_{x_1}}^{k_{x_2}}\sum_{k_y = k_{y_1}}^{k_{y_2}} \Omega(k_x,k_y)\Delta k_x\Delta k_y$$ where $\Delta k_x = (k_{x_2}-k_{x_1})/N_x$ and $N_x\to\infty$, and $k_{i_j}$ are the limits of BZ. When using this approximation, I calculated a transportation coefficient, the result was not correct. I even increase the $N_x (and\: N_y)$ but the result does not improve. The article that I am trying to replicate says that they used Gaussian meshes with 2500 points.
My question is, how do people generally perform numerical integration over BZ? What are Gaussian meshes? What is a proper way for k points sampling? Is the Riemann approximation that I am using the correct way for BZ integration? Are there other approximations for this kind of integration?
