class PEXSISolver(number_of_poles, fermi_level_bracket=None, number_of_electrons_tolerance=None, number_of_inertia_steps=None, maximum_number_of_pexsi_iterations=None, temperature=None, options=None, processes_per_pole=None)

Calculate the density matrix using the PEXSI algorithm.

  • number_of_poles (int) – The number of poles to use. Must be even.
  • fermi_level_bracket (PhysicalQuantity of type energy) – The bracket in which to search for the Fermi level.
    Default: (10, 10)*Hartree
  • number_of_electrons_tolerance (tuple of floats) – The bracket for the tolerance for the number of electrons calculated by the PEXSI algorithm, given as (minimum, maximum).
    Default: (1.e-4, 1.0)
  • number_of_inertia_steps (int) – The number of SCF steps that should use inertia counting.
    Default: 5
  • maximum_number_of_pexsi_iterations (int) – The maximum number of PEXSI iterations per SCF step.
    Default: 10
  • temperature (PhysicalQuantity of type temperature) – The smearing of the Fermi distribution.
    Default: 300*Kelvin
  • processes_per_pole (int > 0) – The number of processes to use per pole. In case the number of processes exceeds the number of poles a larger number of processes per pole will be used.
    Default: 1
Returns:The Fermi level bracket.
Return type:PhysicalQuantity of type energy

Returns the maximum number of PEXSI iterations per SCF cycle.

Returns:The maximum number of PEXSI iterations per SCF cycle.
Return type:int
Returns:The tolerance for the number of electrons as (final, initial) tolerance.
Return type:tuple of floats
Returns:The maximum number of PEXSI iterations per SCF cycle.
Return type:int
Returns:The number of poles.
Return type:int
Returns:The number of processes to use per pole.
Return type:int
Returns:The temperature of the Fermi distribution.
Return type:PhysicalQuantity of type temperature

Usage Example

Define a density matrix method that uses the PEXSI solver, supplying the number of poles to use in the pole expansion.

density_matrix_method = PEXSISolver(number_of_poles=64)

For more control, one can specify various other parameters.

density_matrix_method = PEXSISolver(
    fermi_level_bracket=(-10.0, 10.0) * Hartree,
    number_of_electrons_tolerance=(1.e-4, 1.0),

To use the PEXSISolver when using equivalent bulk to generate the initial density for a device calculation, set it on the AlgorithmParameters passed to the EquivalentBulk object.

algorithm_parameters = AlgorithmParameters(

initial_density_type = EquivalentBulk(algorithm_parameters=algorithm_parameters)

device_algorithm_parameters = DeviceAlgorithmParameters(

device_calculator = DeviceLCAOCalculator(


The implemented pole expansion and selective inversion (PEXSI) method is based on the work by Lin et al., [LLY+09], [LCYH13].

The method calculates the density matrix \(D\) by means of a pole expansion

\[D = \Im \sum_{l=1}^{P} \omega_l(\beta) \left(H - (z_l + \mu)S\right)^{-1}\]

with expansion weights \(\omega_l\) as a function of the inverse temperature \(\beta\) and the shifts \(z_l\). The method calculates only the required elements \(\left(H - (z_l + \mu)S\right)_{ij}^{-1}\) using selective inversion.

Note that for the moment, the implemented method is only usable for the \(\Gamma\)-point and unpolarized calculations.

The computational scaling of this method depends on the dimensionality of the system. For quasi-1D systems like nanotubes, it scales linearly with regards to the number of electrons, i.e. \(\mathcal{O}(N_e)\). For quasi-2D systems it scales with \(\mathcal{O}(N_e^{1.5})\), and for general 3D systems it scales with \(\mathcal{O}(N_e^2)\).

The implemented method is multi-level parallel: On the first level, it parallelizes over all poles in the pole expansion. On the second level, it can also distribute the calculation of each pole over multiple processors. This distribution is done automatically, as soon as the calculation runs on enough processors.

The PEXSI algorithm can be roughly divided into two parts. First, an inertia counting scheme is followed to identify a trial Fermi level within the supplied fermi_level_bracket. Afterwards, the actual pole expansion is done to calculate the density matrix for this trial Fermi level. From this, the number of electrons corresponding to this density matrix can be calculated, allowing the trial Fermi level to be updated once more by doing a Newton step. This PEXSI cycle is repeated at most maximum_number_of_pexsi_iterations times, until the calculated number of electrons corresponds to the exact number of electrons within the current tolerance as explained below. After number_of_inertia_steps SCF steps have elapsed, subsequent steps will not use the inertia counting scheme anymore, as by this point the calculated Fermi level should be close to the true Fermi level.

In this implementation of the PEXSI method, the tolerance on \(\Delta N_e(\mu)\) - which is the deviation from the required number of electrons for a given trial Fermi level - is adapted dynamically based on the current state of the SCF cycle. In the first iteration, the tolerance starts out at the maximum, which is the upper bound of the given number_of_electrons_tolerance bracket. After each iteration, the largest absolute change in the density matrix is computed, and the tolerance for the following iteration is set to that value, with the lower bound of the tolerance as a minimum. By following this scheme, a SCF iteration typically only requires very few PEXSI iterations.

[LCYH13]Lin Lin, Mohan Chen, Chao Yang, and Lixin He. Accelerating atomic orbital-based electronic structure calculation via pole expansion and selected inversion. Journal of Physics: Condensed Matter, 25(29):295501, 2013. doi:10.1088/0953-8984/25/29/295501.
[LLY+09]Lin Lin, Jianfeng Lu, Lexing Ying, Roberto Car, and Weinan E. Fast algorithm for extracting the diagonal of the inverse matrix with application to the electronic structure analysis of metallic systems. Commun. Math. Sci., 7(3):755–777, 09 2009.