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 1s and 0s 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.

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.

I have some function F(kx,ky,kz) that I wish to numerically integrate over a polygon domain - physically, I am integrating over the first Brillouin zone 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 Brillouin zone, 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.

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 1s and 0s 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.