tbplas.Lindhard

class tbplas.Lindhard(cell: Any, kmesh_size: Tuple[int, int, int], energy_min: float = 0.0, energy_max: float = 10.0, energy_step: int = 1000, mu: float = 0.0, temperature: float = 300.0, g_s: int = 1, back_epsilon: float = 1.0, dimension: int = 3, delta: float = 0.005, **kwargs)

Lindhard function calculator.

_cell

primitive cell for which properties will be evaluated, kept as an alias to the ‘model’ attribute of ‘DiagSolver’ class for compatibility

Type:

‘PrimitiveCell’ instance

_omegas

energy grid in eV

Type:

(num_omega,) float64 array

_kmesh_size

dimension of mesh grid in the 1st Brillouin zone

Type:

(3,) int64 array

_kmesh_grid

grid coordinates of k-points on the mesh grid

Type:

(num_kpt, 3) int64 array

_mu

chemical potential in eV

Type:

float

_beta

1/kBT in 1/eV

Type:

float

_g_s

spin degeneracy

Type:

int

_back_epsilon

relative background dielectric constant

Type:

float

_dimension

dimension of the system

Type:

int

_delta

broadening parameter in eV that appears in the denominator of Lindhard function

Type:

float

Notes

  1. Units

‘Lindhard’ class uses eV for energies, elementary charge of electron for charges and nm for lengths. So the unit for dynamic polarization is 1 / (eV * nm**2) in 2d case and 1 / (eV * nm**3) in 3d case. Relative dielectric function is dimensionless. The unit for AC conductivity is e**2/h_bar in 2d case and e**2/(h_bar*nm) in 3d case.

  1. Coulomb potential and dielectric function

The Coulomb potential in Lindhard class takes the following form:
V(r) = 1 / epsilon * e**2 / r

= 1 / (epsilon_0 * epsilon_r) * e**2 / r

where epsilon is the ABSOLUTE dielectric electric constant of background.

Now we evaluate the absolute dielectric constant of vacuum in the units of eV/nm/q. Suppose that we have two electrons separated at the distance of 1nm, then their potential energy in SI units is:

V = 1 / (4 * pi * epsilon_0) * q**2 / 1e-9 [Joule]

= 1 / (4 * pi * epsilon_0) * q**2 / 1e-9 * 1 / q [eV] = q * 1e9 / (4 * pi * epsilon_0) [eV]

The counterpart in eV/nm/q units is:

V = 1 / epsilon_0_arb [eV]

So we have:

1 / epsilon_0_arb = q * 1e9 / (4 * pi * epsilon_0)

where the right part are the values in SI. Finally, we have:

epsilon_0_arb = 0.6944615417149689

in eV/nm/q units. That’s EPSILON0 in constants module.

  1. Fourier transform of Coulomb potential

The 3D Fourier transform of Coulomb potential V(r) as defined in section 2 takes the following form:

V(q) = 4 * pi * e**2 / (epsilon_0 * epsilon_r * q**2)

For the derivation, see pp. 19 of Many-Particle Physics by Gerald D. Mahan.

The 2D Fourier transform takes the form:

V(q) = 2 * pi * e**2 / (epsilon_0 * epsilon_r * q)

See the following url for the derivation: https://math.stackexchange.com/questions/3627267/fourier-transform-of-the-2d-coulomb-potential

  1. System dimension

‘Lindhard’ class deals with system dimension in two approaches. The first approach is to treat all systems as 3-dimensional. Supercell technique is required in this approach, with vacuum layers added on non-periodic directions. Also, the component(s) of kmesh_size should be set to 1 accordingly.

The second approach utilizes dimension-specific formula whenever possible. For now, only 2-dimensional case has been implemented. This approach requires that the system should be periodic in xOy plane, i.e. the non-periodic direction should be along ‘c’ axis. Otherwise, the results will be wrong.

The first approach suffers from the problem that dynamic polarization and AC conductivity scale inversely proportional to the product of supercell lengths, i.e. |c| in 2d case and |a|*|b| in 1d case. This is due to the elementary volume in reciprocal space (dnk) in the summation in Lindhard function. On the contrary, the second approach has no such issue. If the supercell lengths of non-periodic directions are set to 1 nm, then the first approach yields the same results as the second approach.

For the dielectric function, the situation is more complicated. From the equation

epsilon(q) = 1 - V(q) * dyn_pol(q)

we can see it is also affected by the Coulomb potential V(q). Since

V(q) = 4 * pi * e**2 / (epsilon_0 * epsilon_r * q**2) (3d)

and

V(q) = 2 * pi * e**2 / (epsilon_0 * epsilon_r * q) (2d)

the influence of system dimension is q-dependent. Setting supercell lengths to 1 nm will NOT produce the same result as the second approach.

  1. Dielectric function for q = 0

The dielectric function for q = 0 cannot be evaluated from dynamic polarizability as v(q) diverges. Instead, we calculate it from ac conductivity. We have two formulae. The first one is eqn. (63)

epsilon = 1 + 1j * 4 * pi / omega * sigma

in https://www.openmx-square.org/tech_notes/Dielectric_Function_YTL.pdf, while the second one is eqn. (1.10)

epsilon = 1 + 1j / (epsilon_0 * omega) * sigma

in Plasmonics: Fundamentals and Applications by Stefan A. Maier. The former is in Gaussian units, while the latter is in SI units. Unfortunately, we cannot apply them directly.

As we take the Coulomb potential to be
V(r) = 1 / epsilon * e**2 / r

= 1 / (epsilon_0 * epsilon_r) * e**2 / r

there should be a 4pi factor as in Gaussian units. However, epsilon_0 is not 1 as we are not using cm/g/s as the fundamental units. The actual formula is

epsilon = 1 + 1j * 4 * pi / (epsilon_0 * omega) * sigma

which can be regarded as a mixture of the formulae in Gaussian and SI units.

The unit of epsilon follows:
[epsilon] = [sigma] / [epsilon_0 * omega]

= e**2/(h_bar*nm) / [e**2/(eV*nm) * eV/h_bar] = e**2/(h_bar*nm) / [e**2/(h_bar*nm)] = 1

However, this equation holds for 3d systems only. For 2d systems there will be a remaining unit of nm.

__init__(cell: Any, kmesh_size: Tuple[int, int, int], energy_min: float = 0.0, energy_max: float = 10.0, energy_step: int = 1000, mu: float = 0.0, temperature: float = 300.0, g_s: int = 1, back_epsilon: float = 1.0, dimension: int = 3, delta: float = 0.005, **kwargs) → None
Parameters:

  • cell – primitive cell for which properties will be evaluated

  • kmesh_size – dimension of mesh grid in 1st Brillouin zone

  • energy_min – low bound of energy grid for evaluating response properties

  • energy_max – upper bound of energy grid for evaluating response properties

  • energy_step – resolution of energy grid

  • mu – chemical potential of the cell in eV

  • temperature – temperature in Kelvin

  • g_s – spin degeneracy

  • back_epsilon – relative background dielectric constant

  • dimension – dimension of the system

  • delta – broadening parameter in eV

  • kwargs – arguments for DiagSolver.__init__

Raises:

ValueError – if kmesh_size and dimension are not properly set

Methods

__init__(cell, kmesh_size[, energy_min, ...])

param cell:

primitive cell for which properties will be evaluated

all_average(data_local)

Average results over random samples broadcast to all process.

all_reduce(data_local)

Reduce local data and broadcast to all processes.

average(data_local)

Average results over random samples and store results to master process.

barrier()

Wrapper for self.comm.Barrier.

bcast(data_local)

Broadcast data from master to other processes.

calc_ac_cond([component, use_fortran])

Calculate AC conductivity using Kubo-Greenwood formula.

calc_bands(k_points[, convention, solver, ...])

Calculate band structure along given k_path.

calc_dos(k_points[, e_min, e_max, e_step, ...])

Calculate density of states for given energy range and step.

calc_dyn_pol_arbitrary(q_points[, use_fortran])

Calculate dynamic polarizability for arbitrary q-points.

calc_dyn_pol_regular(q_points[, use_fortran])

Calculate dynamic polarizability for regular q-points on k-mesh.

calc_epsilon(q_points, dyn_pol)

Calculate dielectric function for given q-points from dynamic polarization.

calc_epsilon_q0(omegas, ac_cond)

Calculate dielectric function from AC conductivity for q=0.

calc_states(k_points[, convention, solver, ...])

Calculate energies and wave functions on k-points.

cart2frac(cart_coord[, unit])

Convert CARTESIAN coordinates of k-points in 1/unit to FRACTIONAL coordinates.

dist_bound(n_max)

Same as dist_range, but returns the lower and upper bounds.

dist_list(raw_list[, algorithm])

Distribute a list over processes.

dist_range(n_max)

Distribute range(n_max) over processes.

frac2cart(frac_coord[, unit])

Convert FRACTIONAL coordinates of k-points to CARTESIAN coordinates in 1/unit.

grid2cart(grid_coord[, unit])

Convert GRID coordinates of k-points to CARTESIAN coordinates in 1/unit.

grid2frac(grid_coord)

Convert GRID coordinates of k-points to FRACTIONAL coordinates.

log([event, fmt])

Log the date and time of event.

print([text])

Print text on master process.

reduce(data_local)

Reduce local data to master process.

wrap_frac(k_points)

Wrap fractional coordinates of k-points back into 1st Brillouin zone.

wrap_grid(k_points)

Wrap grid coordinates of k-points back into 1st Brillouin zone.

Attributes

is_master

Determine whether this is the master process.

lat_vec

Get the lattice vectors of the model.

model_is_pc

Check if the model is a primitive cell.

mpi_enabled

Determine whether MPI is enabled.

num_orb

Get the number of orbitals in the model.

rank

Interface for the '__rank' attribute.

recip_lat_vec

Get the reciprocal lattice vectors of the model in 1/nm.

size

Interface for the '__size' attribute.
__init__(cell: Any, kmesh_size: Tuple[int, int, int], energy_min: float = 0.0, energy_max: float = 10.0, energy_step: int = 1000, mu: float = 0.0, temperature: float = 300.0, g_s: int = 1, back_epsilon: float = 1.0, dimension: int = 3, delta: float = 0.005, **kwargs) → None
Parameters:

  • cell – primitive cell for which properties will be evaluated

  • kmesh_size – dimension of mesh grid in 1st Brillouin zone

  • energy_min – low bound of energy grid for evaluating response properties

  • energy_max – upper bound of energy grid for evaluating response properties

  • energy_step – resolution of energy grid

  • mu – chemical potential of the cell in eV

  • temperature – temperature in Kelvin

  • g_s – spin degeneracy

  • back_epsilon – relative background dielectric constant

  • dimension – dimension of the system

  • delta – broadening parameter in eV

  • kwargs – arguments for DiagSolver.__init__

Raises:

ValueError – if kmesh_size and dimension are not properly set

all_average(data_local: ndarray) → ndarray

Average results over random samples broadcast to all process.

Parameters:

data_local – local results on each process

Returns:

averaged data from data_local

all_reduce(data_local: ndarray) → ndarray

Reduce local data and broadcast to all processes.

Parameters:

data_local – local results on each process

Returns:

summed data from data_local

average(data_local: ndarray) → ndarray

Average results over random samples and store results to master process.

Parameters:

data_local – local results on each process

Returns:

averaged data from data_local

barrier() → None

Wrapper for self.comm.Barrier.

bcast(data_local: ndarray) → None

Broadcast data from master to other processes.

Parameters:

data_local – local results on each process

Returns:

None

calc_ac_cond(component: str = 'xx', use_fortran: bool = True) → Tuple[ndarray, ndarray]

Calculate AC conductivity using Kubo-Greenwood formula.

Reference: section 12.2 of Wannier90 user guide. NOTE: there is not such g_s factor in the reference.

Parameters:

  • component – which component of conductivity to evaluate should be in “xx”, “xy”, “xz”, “yx”, “yy”, “yz”, “zx”, “zy” and “zz”

  • use_fortran – choose between FORTRAN and Cython backends

Returns:

(omegas, ac_cond) omegas: (num_omega,) float64 array energies in eV ac_cond: (num_omega,) complex128 array AC conductivity in e**2/(h_bar*nm) in 3d case and e**2/h_bar in 2d case

Raises:

ValueError – if component is illegal

calc_bands(k_points: ndarray, convention: int = 1, solver: str = 'lapack', orbital_indices: Union[Iterable[int], ndarray] = None, **kwargs) → namedtuple

Calculate band structure along given k_path.

The result is a named tuple with the following attributes: k_len: (num_kpt,) float64 array, length of k-path bands: (num_kpt, num_bands) float64 array, band structure proj: (num_kpt, num_bands) float64 array, projection on selected orbitals

Parameters:

  • k_points – (num_kpt, 3) float64 array FRACTIONAL coordinates of the k-points

  • convention – convention for setting up the Hamiltonian

  • solver – solver type, should be “lapack” or “arpack”

  • orbital_indices – indices of orbitals for evaluating projection, take all orbitals into consideration if no orbitals are specified

  • kwargs – arguments for the arpack solver

Returns:

k_len, band structure and projection packed in named tuple

Raises:

ValueError – if solver is neither lapack nor arpack

calc_dos(k_points: ndarray, e_min: float = None, e_max: float = None, e_step: float = 0.05, sigma: float = 0.05, basis: str = 'Gaussian', g_s: int = 1, **kwargs) → Tuple[ndarray, ndarray]

Calculate density of states for given energy range and step.

Parameters:

  • k_points – (num_kpt, 3) float64 array FRACTIONAL coordinates of k-points

  • e_min – lower bound of the energy range in eV

  • e_max – upper hound of the energy range in eV

  • e_step – energy step in eV

  • sigma – broadening parameter in eV

  • basis – basis function to approximate the Delta function

  • g_s – spin degeneracy

  • kwargs – arguments for ‘calc_bands’

Returns:

(energies, dos) energies: (num_grid,) float64 array energy grid corresponding to e_min, e_max and e_step dos: (num_grid,) float64 array density of states in states/eV

Raises:

ValueError – if basis is neither Gaussian nor Lorentzian, or the solver is neither lapack nor arpack

calc_dyn_pol_arbitrary(q_points: ndarray, use_fortran: bool = True) → Tuple[ndarray, ndarray]

Calculate dynamic polarizability for arbitrary q-points.

Parameters:

  • q_points – (num_qpt, 3) float64 array CARTESIAN coordinates of q-points in 1/NM

  • use_fortran – choose between FORTRAN and Cython backends

Returns:

(omegas, dyn_pol) omegas: (num_omega,) float64 array energies in eV dyn_pol: (num_qpt, num_omega) complex128 array dynamic polarizability in 1/(eV*nm**2) or 1/(eV*nm**3)

calc_dyn_pol_regular(q_points: ndarray, use_fortran: bool = True) → Tuple[ndarray, ndarray]

Calculate dynamic polarizability for regular q-points on k-mesh.

Parameters:

  • q_points – (num_qpt, 3) int64 array GRID coordinates of q-points

  • use_fortran – choose between FORTRAN and Cython backends

Returns:

(omegas, dyn_pol) omegas: (num_omega,) float64 array energies in eV dyn_pol: (num_qpt, num_omega) complex128 array dynamic polarizability in 1/(eV*nm**2) or 1/(eV*nm**3)

calc_epsilon(q_points: ndarray, dyn_pol: ndarray) → ndarray

Calculate dielectric function for given q-points from dynamic polarization.

Parameters:

  • q_points – (num_qpt, 3) float64 array CARTESIAN coordinates of q-points in 1/NM

  • dyn_pol – (num_qpt, num_omega) complex128 array dynamic polarization from calc_dyn_pol_*

Returns:

(num_qpt, num_omega) complex128 array relative dielectric function

calc_epsilon_q0(omegas: ndarray, ac_cond: ndarray) → ndarray

Calculate dielectric function from AC conductivity for q=0.

Parameters:

  • omegas – (num_omega,) float64 array energies in eV

  • ac_cond – (num_omega,) complex128 array AC conductivity in e**2/(h_bar*nm) in 3d case

Returns:

(num_omega,) complex128 array relative dielectric function

Raises:

ValueError – if dimension is not 3

calc_states(k_points: ndarray, convention: int = 1, solver: str = 'lapack', all_reduce: bool = True, **kwargs) → Tuple[ndarray, ndarray]

Calculate energies and wave functions on k-points.

Parameters:

  • k_points – (num_kpt, 3) float64 array FRACTIONAL coordinates of the k-points

  • convention – convention for setting up the Hamiltonian

  • solver – solver type, should be “lapack” or “arpack”

  • all_reduce – whether to call MPI_Allreduce to synchronize the results on each process

  • kwargs – arguments for the arpack solver

Returns:

(bands, states) bands: (num_kpt, num_bands) float64 array energies on the k-points states: (num_kpt, num_bands, num_orb) complex128 array wave functions on the k-points, each ROW of states[i_k] is a wave function

Raises:

ValueError – if solver is neither lapack nor arpack

cart2frac(cart_coord: ndarray, unit: float = 1.0) → ndarray

Convert CARTESIAN coordinates of k-points in 1/unit to FRACTIONAL coordinates.

Parameters:

  • cart_coord – (num_kpt, 3) float64 array CARTESIAN coordinates in 1/unit

  • unit – scaling factor from unit to NANOMETER, e.g. unit=0.1 for ANGSTROM

Returns:

(num_kpt, 3) float64 array FRACTIONAL coordinates of k-points

dist_bound(n_max: int) → Tuple[int, int]

Same as dist_range, but returns the lower and upper bounds. Both of the bounds are close, i.e. [i_min, i_max].

Parameters:

n_max – upper bound of range

Returns:

lower and upper bounds of subrange assigned to this process

dist_list(raw_list: List[Any], algorithm: str = 'range') → List[Any]

Distribute a list over processes.

Parameters:

  • raw_list – raw list to distribute

  • algorithm – distribution algorithm, should be either “remainder” or “range”

Returns:

sublist assigned to this process

dist_range(n_max: int) → range

Distribute range(n_max) over processes.

Parameters:

n_max – upper bound of the range

Returns:

subrange assigned to this process

frac2cart(frac_coord: ndarray, unit: float = 1.0) → ndarray

Convert FRACTIONAL coordinates of k-points to CARTESIAN coordinates in 1/unit.

Parameters:

  • frac_coord – (num_kpt, 3) float64 array FRACTIONAL coordinates of k-points

  • unit – scaling factor from unit to NANOMETER, e.g. unit=0.1 for ANGSTROM

Returns:

(num_kpt, 3) float64 array CARTESIAN coordinates of k-points in 1/unit

grid2cart(grid_coord: ndarray, unit: float = 1.0) → ndarray

Convert GRID coordinates of k-points to CARTESIAN coordinates in 1/unit.

Parameters:

  • grid_coord – (num_kpt, 3) int64 array GRID coordinates of k-points

  • unit – scaling factor from unit to NANOMETER, e.g. unit=0.1 for ANGSTROM

Returns:

(num_kpt, 3) float64 array CARTESIAN coordinates of k-points in 1/unit

grid2frac(grid_coord: ndarray) → ndarray

Convert GRID coordinates of k-points to FRACTIONAL coordinates.

Parameters:

grid_coord – (num_kpt, 3) int64 array GRID coordinates of k-points

Returns:

(num_kpt, 3) float64 array FRACTIONAL coordinates of k-points

log(event: str = '', fmt: str = '%x %X') → None

Log the date and time of event.

Parameters:

  • event – notice of the event

  • fmt – date and time format

Returns:

None.

print(text: str = '') → None

Print text on master process.

NOTE: flush=True is essential for some MPI implementations, e.g. MPICH3.

Parameters:

text – text to print

Returns:

None

reduce(data_local: ndarray) → ndarray

Reduce local data to master process.

Parameters:

data_local – local results on each process

Returns:

summed data from data_local

static wrap_frac(k_points: ndarray) → ndarray

Wrap fractional coordinates of k-points back into 1st Brillouin zone.

Parameters:

k_points – (num_kpt, 3) float64 array FRACTIONAL coordinates of k-points

Returns:

(num_kpt, 3) float64 array wrapped FRACTIONAL coordinates of k-points

wrap_grid(k_points: ndarray) → ndarray

Wrap grid coordinates of k-points back into 1st Brillouin zone.

Parameters:

k_points – (num_kpt, 3) int64 array GRID coordinates of k-points

Returns:

(num_kpt, 3) int64 array wrapped GRID coordinates of k-points

__weakref__

list of weak references to the object (if defined)

property is_master: bool

Determine whether this is the master process.

property lat_vec: ndarray

Get the lattice vectors of the model.

Returns:

(3, 3) float64 array lattice vectors of in nm

property model_is_pc: bool

Check if the model is a primitive cell.

Returns:

whether the model is a primitive cell

property mpi_enabled: bool

Determine whether MPI is enabled.

property num_orb: int

Get the number of orbitals in the model.

Returns:

number of orbitals

property rank: int

Interface for the ‘__rank’ attribute.

Returns:

rank of this MPI process

property recip_lat_vec: ndarray

Get the reciprocal lattice vectors of the model in 1/nm.

Returns:

(3, 3) float64 array reciprocal lattice vectors of in nm

property size: int

Interface for the ‘__size’ attribute.

Returns:

number of MPI processes