tbplas.Lindhard

class tbplas.Lindhard(model: PrimitiveCell, overlap: PrimitiveCell = None, echo_details: bool = True)

Lindhard function calculator.

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/hbar in 2d case and e**2/(hbar*nm) in 3d case.

2. 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.

3. 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

4. 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.

5. 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/(hbar*nm) / [e**2/(eV*nm) * eV/hbar] = e**2/(hbar*nm) / [e**2/(hbar*nm)] = 1

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

__init__(model: PrimitiveCell, overlap: PrimitiveCell = None, echo_details: bool = True) → None

Parameters:

• model – model for which properties will be calculated

• overlap – model holding the overlap of basis functions

• echo_details – whether to output parallelization details

Methods

__init__(model[, overlap, echo_details])
calc_ac_cond()	Calculate AC conductivity using Kubo-Greenwood formula.
calc_bands()	Calculate band structure along given k_path. :return: (k_len, bands) k_len: (num_kpt,) float64 array, length of k-path bands: (num_kpt, num_bands) float64 array, band structure.
calc_bands_scipy()	Calculate band structure along given k_path.
calc_dos()	Calculate density of states for given energy range and step. :return: (energies, dos) energies: (num_eng,) float64 array energies determined by e_min, e_max and e_step dos: (num_eng,) float64 array density of states in states/eV.
calc_dyn_pol()	Calculate dynamic polarizability for arbitrary q-points. :return: (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(dyn_pol)	Calculate dielectric function from dynamic polarization.
calc_epsilon_q0(omegas, ac_cond)	Calculate dielectric function from AC conductivity for q=0. :param omegas: (num_omega,) float64 array energies in eV :param ac_cond: (num_omega,) complex128 array AC conductivity in e**2/(hbar*nm) in 3d case :return: (num_omega,) complex128 array relative dielectric function.
calc_states()	Calculate bands as well as eigenstates.
k_cart2frac(cart_coord[, unit])	Convert CARTESIAN coordinates of k-points in 1/unit to FRACTIONAL coordinates. :param cart_coord: (num_kpt, 3) float64 array CARTESIAN coordinates in 1/unit :param unit: scaling factor from unit to NANOMETER, e.g. unit=0.1 for ANGSTROM :return: (num_kpt, 3) float64 array FRACTIONAL coordinates of k-points.
k_frac2cart(frac_coord[, unit])	Convert FRACTIONAL coordinates of k-points to CARTESIAN coordinates in 1/unit. :param frac_coord: (num_kpt, 3) float64 array FRACTIONAL coordinates of k-points :param unit: scaling factor from unit to NANOMETER, e.g. unit=0.1 for ANGSTROM :return: (num_kpt, 3) float64 array CARTESIAN coordinates of k-points in 1/unit.

Attributes

is_master

Wrapper for determine the master process for compatibility.

__init__(model: PrimitiveCell, overlap: PrimitiveCell = None, echo_details: bool = True) → None

Parameters:

• model – model for which properties will be calculated

• overlap – model holding the overlap of basis functions

• echo_details – whether to output parallelization details

calc_ac_cond() → 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.

Returns:

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

calc_bands() → Tuple[ndarray, ndarray]

Calculate band structure along given k_path. :return: (k_len, bands)

k_len: (num_kpt,) float64 array, length of k-path bands: (num_kpt, num_bands) float64 array, band structure

calc_bands_scipy() → Tuple[ndarray, ndarray]

Calculate band structure along given k_path.

This function calls C++ interface to set up the Hamiltonian. However, the diagonalization is done using scipy. For now this function is intended to be called by the ‘ParamFit’ class. It also works as an example for obtaining the Hamiltonian for post-process.

Returns:

(k_len, bands) k_len: (num_kpt,) float64 array, length of k-path bands: (num_kpt, num_bands) float64 array, band structure

calc_dos() → Tuple[ndarray, ndarray]

Calculate density of states for given energy range and step. :return: (energies, dos)

energies: (num_eng,) float64 array energies determined by e_min, e_max and e_step dos: (num_eng,) float64 array density of states in states/eV

calc_dyn_pol() → Tuple[ndarray, ndarray]

Calculate dynamic polarizability for arbitrary q-points. :return: (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(dyn_pol: ndarray) → ndarray

Calculate dielectric function from dynamic polarization.

Parameters:

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. :param omegas: (num_omega,) float64 array

energies in eV

Parameters:

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

Returns:

(num_omega,) complex128 array relative dielectric function

calc_states() → Tuple[ndarray, ndarray]

Calculate bands as well as eigenstates.

Returns:

(bands, states) bands: (num_kpt, num_bands) float64 array

band structure

states: (num_kpt, num_bands, num_orb) complex128 array

eigenstates at each k-point, with each ROW being a state

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

Convert CARTESIAN coordinates of k-points in 1/unit to FRACTIONAL coordinates. :param cart_coord: (num_kpt, 3) float64 array

CARTESIAN coordinates in 1/unit

Parameters:

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

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

Convert FRACTIONAL coordinates of k-points to CARTESIAN coordinates in 1/unit. :param frac_coord: (num_kpt, 3) float64 array

FRACTIONAL coordinates of k-points

Parameters:

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

__weakref__

list of weak references to the object

property is_master: bool

Wrapper for determine the master process for compatibility.