(Fortran) Integrating/summing over complicated 3D domain

Question

I have some function $F(k_x,k_y,k_z)$ that I wish to numerically integrate over a polygon domain - physically, I am integrating over the first Brillouin Zone (BZ) of the FCC lattice (a truncated octahedron).

My goal is to tell Fortran what region of $k$-space and which points to integrate over, given this shape.

I can write plane equations and inequalities for this region and perhaps stick them into a logical structure, then integrate over the whole BZ, but that seems inefficient to me. It seems that I should be able to use any available symmetries to pick out a unique bit of the BZ, then reflect/rotate/etc my answer using those symmetries.

I was told that there already exist routines which can do something like this, but I am unsure if what I googled are appropriate (VASP, etc).

Would anyone be kind enough to suggest an appropriate package(s) if one exists? If not, are there perhaps multidimensional integration methods I should investigate? If so, is there an efficient way to restrict my integration domain?

Another avenue is to simply sum $F(k)$ over a uniform grid within the zone - would it then be more efficient to pre-calculate a rank-3 array with 1's and 0's in my desired shape, call, calculate and sum $F(k)$ only for nonzero elements? Just trying to get a feel for which direction would be most efficient and accurate.

Answer 1

It would seem to me that you can subdivide your domain into a relatively small set of tetrahedra. Then, it is trivial to do the integration because there are many good (and pre-tabulated) quadrature rules on tetrahedra that will yield reasonably high accuracy -- this is what one does in the finite element method all the time, so there is a lot of information out there. The integral over your Brillouin zone is then simply the sum of the integrals over the tetrahedra you have decomposed your zone into.

If the function you want to integrate, $F(k_x,k_y,k_z)$, happens to be highly oscillatory or otherwise not very smooth, then the solution is to split each of your tetrahedra into four smaller tetrahedra, and if necessary repeat the process, and integrate over each of the now smaller cells.

Answer 2

Many plane-wave codes use a very simple k-point weighting scheme:

1. Generate a uniform mesh in k-space, $K$
2. Assign each k-point in the irredicuble wedge a weight $w(k) = \#\{q \in K : \text{k and q are symmetric}\}$
3. Write integrals as $\int_{BZ} F(k) dk = \sum_{k \in K} w(k) F(k)$

This is quick and dirty, but it usually isn't the leading error (for PWDFT, at least).

You can find an example of such an integral being used to sum the electron number given a chemical potential (to find the Fermi energy) in Quantum Espresso's sumkg. I believe the code that finds the irreducible wedge and generates the k-point weights, $w(k)$, is here. Quantum Espresso can also use tetrahedra; the equivalent integral of sumkg is called sumkt. Routines to generate the tetrahedra and uniform meshes in k-space are here.