This repository implements the (1,1) rational integration method discussed in the following paper:

Kaprzyk, S. "An algebraic formulation of the simplex linear method for multiple integrals over the d-dimensional domain." Computer Physics Communications 183, no. 1 (2012): 20-25.

What this paper discusses is an integration routine for integrals of the form

I(z) = \int_V [d^D k] \frac{f(k)}{z - g(k)} over a volume V, as well as computing $K(z) = \int [dz] I(z)$ efficiently. f(k), g(k), and z are allowed to be complex. The user inputs these functions evaluated on a discrete grid of points that are space-filling for the volume V. The user can also input weights for each k-point, for example if they have symmetry-reduced their necessary integration region.

Applications of these integrals include:

• Computing spectral functions (density of states) of condensed matter systems
• Computing the Fermi energy of solids
• Computing response functions (transport, susceptibility, etc) of systems
• Hilbert transforms

Motivation:

We note due to the structure of the integrand, if there are k such that g(k) approaches z, these will be singular contributions. Regular integration methods such as the trapezoid method, even iteratively refined, requires many points centered around these singularities (or near singularities). Additionally, these functions f(k), g(k) may be the outputs of expensive numerical methods that cannot reasonably produce values at arbitrarily many points.

Although versions of this method exist for 3D regions exist, called the "tetrahedron method", most of them assume the numerator $f(k) = 1$. The ones that allow nonzero f(k) don't allow for complex-valued f(k), g(k), z. Our goal is to ameliorate that, as well as to extend the method to arbitrary dimensional integrals. Our first goal is to do it up to 6D integrals. Realistically I don't see it extending past 8D integrals.

Details:

The linear simplex method takes the integration region (specified by the user's points), constructs simplices filling the space between the points, and takes the values f(k), g(k) and constructs an N-linear approximation for them inside the simplices. This is not a linear approximation to the integrand (which is what the trapezoid method is). We are separately linearly interpolating the numerator and denominator. The integral is approximated by summing a series of weights for each simplex, and inside each tetrahedron the integrand F_n(k)/(z-G_n(k)) is analytically integrated, where F_n(k), G_n(k) are now the linear approximations to f(k), g(k) inside each simplex n.

Necessary Prerequisites:

• A modern Fortran compiler (I know gfortran version 9.3.0+ works; I know gfortran version 4.8.5 doesn't work.)
• scipy (for scipy.Delaunay)
• numpy version 1.21.4+ (to use f2py.get_include()) That's it!

Optional Prerequisites: To use some of the helper scripts and run some of the advanced tests, you may need:

• spglib
• Atomic Simulation Environment (ASE)
• Quantum Espresso
• (TODO) CUDA is required to compile the GPU-accelerated version.

However, there are some tests that can be run that do not require these packages.

Install:

• Go into the /src directory. Modify the compile_scripts/compile_applications.sh file. Things to modify are the compiler/compiler options and the python-related variables.
• From the /src directory, run the compile_scripts/compile_applications.sh, perhaps as 'source compile_scripts/compile_applications.sh'.
• This should compile the code and also create the python-importable libraries, which end up in the ROOT/shared_libs directory.

Python Usage:

• Set LD_LIBRARY_PATH and PYTHONPATH to include the ROOT/shared_libs folder.
• Open python, and simply 'import [module]'. For example, 'import integrate_mod' will import the integration routine. An example of this is in Examples/Pb_Quantum_Espresso/python_calling_ex.py.

Fortran Usage:

• The user has a couple choices. The first is to treat this suite of modules as their own code, and simply link the .o files in /src/core_modules/obj/.o to their own program.
• The second choice is to link the shared libraries in src/shared_libs.

Acceleration: (TODO) There are also parallel versions of the executables one can make.

• CPU parallelization is accomplished with Fortran coarrays. In theory, this means you do not have to run it with mpirun.
• GPU parallelization is accomplished with CUDA.

Application Scripts:

• Finding Fermi energy for a system (translation: finding the real z such that K(z) = N, for a given N)
• Preprocessing interpolation script using radial basis functions. In case the input data is expensive to generate, preprocessing it with an interpolation scheme that's higher-order than linear may lead to more accurate or smoother results.

Examples:

• Finding the density of states of a k^2 energy dispersion example over a cubic cell. Three cases are: -- With unit numerator -- With a Fermi function numerator -- With a slightly randomized mesh (illustrating the power of this code to handle non-uniform meshes)
• Finding the density of states and Fermi energy of a nearest-neighbor (hyper)cubic tight-binding model. Three cases: -- 2D case -- 3D case -- 4D case
• Computing the charge susceptibility of a model
• Computing the electron-phonon perturbative renormalization of electron energies: -- Of a model system -- Of a real system whose bandstructure was computed with QE, alongside a dispersionless phonon and uniform electron-phonon coupling -- Of a real system with a momentum-dependent electron-phonon coupling