Mobility¶

class
Mobility
(configuration, electron_phonon_coupling=None, kpoints=None, method=None, temperature=None, phonon_modes=None, electron_bands=None, fermi_shift=None, energies=None, inverse_relaxation_time=None, mobility_object=None, integration_method=None, refinement=None, relaxation_time_approximation=None, calculate_hall_coefficients=None, energy_broadening=None)¶ Analysis class for calculating the mobility for a bulk configuration.
Parameters:  configuration (
BulkConfiguration
) – The bulk configuration for which to calculate the mobility.  electron_phonon_coupling (
ElectronPhononCoupling
) – An electronphonon coupling object used for calculating mobilities. Whenmethod=Full
,electron_phonon_coupling
must be provided.  method (
Full
Isotropic
.) – Keyword specifying if the results should be calculated with the full electronphonon coupling including angular dependence (Full
) or with an isotropic scattering rate (Isotropic
). Whenmethod=Full
,electron_phonon_coupling
must be provided. Whenmethod=Isotropic
,kpoints
must be provided together withenergies
andinverse_relaxation_time
.  kpoints (
MonkhorstPackGrid
RegularKpointGrid
 list of list of float) – List of fractional kpoints. Whenmethod=Isotropic
,kpoints
must be provided.
Default:None
.  temperature (PhysicalQuantity of type temperature) – The temperature for the thermal smearing of the Bose and Fermi distributions.
Default:300 * Kelvin
 phonon_modes (list of ints 
All
) – Phonon modes to include as a list of indices orAll
.
Default:All
(Include all phonon modes available in theelectron_phonon_coupling
)  electron_bands (list of int 
All
) – The band indices of the Bloch states to include.
Default:All
(Include all bands available in theelectron_phonon_coupling
)  fermi_shift (PhysicalQuantity of type energy) – The Fermi shift.
Default:0.0 * eV
 energies (PhysicalQuantity of type energy) – List of energies for which to define the scattering rate. This parameter is only used,
if
method=Isotropic
.
Default:numpy.linspace(0.5,0.5,100)*eV
 inverse_relaxation_time (PhysicalQuantity of type energy) – List of scattering rates corresponding to the provided energies, or a constant
scattering rate. If a list is provided the length must be the same as the
energies
. If amobility_object
is also provided, theinverse_relaxation_time
will be added to the scattering rate obtained from themobility_object
. This parameter is only used ifmethod=Isotropic
.
Default:0.0/Second
 mobility_object (
Mobility
) – A mobility object calculated with theFull
method. This can be used to calculate the scattering rate. If ainverse_relaxation_time
is also provided, theinverse_relaxation_time
will be added to the scattering rate obtained from themobility_object
.  integration_method (
GaussianBroadening()
TetrahedronMethod
) – The method to use for calculating the qintegral.
Default:GaussianBroadening(3e2 * eV)
 refinement (Positive float) – The number of times the qgrid is refined in each direction.
Default:1
(no interpolation)  calculate_hall_coefficients (bool.) – Boolean to determine if the Hall coefficients should be calcuated or not.
Default:True
 energy_broadening (PhysicalQuantity of type energy) –
Deprecated since version 2018.0.
The width of the Gaussian function approximating the deltafunction.
Default:3e2 * eV
 relaxation_time_approximation (not used) –
Deprecated since version 2017.0: Should no longer be used.

carrierDensities
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the carrier density. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated carrier density. If ‘contribution’ is Hole
orElectron
, the hole or electron carrier density is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron carrier densities as the first and second element. In case ofSum
, the total carrier density is returned.Return type: PhysicalQuantity with the unit Meter**3

conductivity
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the conductivity. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated conductivity. If ‘contribution’ is Hole
orElectron
, the hole or electron conductivity is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron conductivities as the first and second element. In case ofSum
, the total conductivity is returned.Return type: PhysicalQuantity with the unit Siemens / Meter

conductivityTensor
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the conductivity tensor. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated conductivity tensor. If ‘contribution’ is Hole
orElectron
, the hole or electron conductivity tensor is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron conductivity tensors as the first and second element. In case ofSum
, the total conductivity tensor is returned.Return type: PhysicalQuantity with the unit Siemens / Meter

eigenvaluesK
()¶ Returns: The electron Bloch state energies for each kpoint. The shape is (spin, kpoints, electron bands). Return type: PhysicalQuantity of type energy

electronBands
()¶ Returns: The band indices of the Bloch states included in the electronphonon coupling calculation. Return type: list of nonnegative ints

energies
()¶ Returns: The user inputted energies list with shape (n_energies). Return type: PhysicalQuantity of type energy.

energyBroadening
()¶ Returns: The width of the Gaussian function approximating the deltafunction, if the integration method is GaussianBroadening. If the TetrahedronMethod is used, it returns 0.0*eV. Return type: PhysicalQuantity of type energy

energyZero
()¶ The energy zero is equal to the Fermi level. For fixed spin moment calculations it is the average of the Fermi level for spin up and spin down.
Returns: The energy zero. Return type: PhysicalQuantity of type energy

evaluate
()¶ Returns: The calculated hole and electron mobilities as the first and second element, respectively. Return type: PhysicalQuantity with the unit Meter**2 / (Volt * Second)

fermiLevel
(spin=None)¶

fermiShift
()¶ Returns: The Fermi shift. Return type: PhysicalQuantity of type energy

firstMomentTensor
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the first moment tensor. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated first moment tensor. If ‘contribution’ is Hole
orElectron
, the hole or electron first moment tensor is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron first moment tensors as the first and second element. In case ofSum
, the total first moment tensor is returned.Return type: PhysicalQuantity with the unit Volt * Siemens / (Meter * Kelvin)

generateEnergyDependentInverseRelaxationTime
(energies, energy_broadening=None)¶ Generator method for the inverse relaxation time as a function of energy.
Parameters:  energies (PhysicalQuantity of type energy) – The energies to calculate the scattering rate for, given as an energy list.
 energy_broadening (PhysicalQuantity of type energy) – The energy broadening of deltafuctions.
Default:3.0e2*eV
Returns: The inverse relaxation time. The length is len(
energies
).Return type: PhysicalQuantity of type inverse time.

hallCoefficientTensor
()¶ Returns: The calculated hole and electron Hallcoefficient tensors as the first and second element of the first entry, respectively. The shape is (carrier type, induced electric field direction, applied charge current density direction, applied magnetic field direction). Return type: PhysicalQuantity with the unit Meter**3 / Coulomb

hallConductivityTensor
()¶ Returns: The calculated hole and electron Hallconductivity tensors as the first and second element of the first entry, respectively. The shape is (carrier type, charge current density direction, electric field direction, magnetic field direction). Return type: PhysicalQuantity with the unit Meter / Coulomb * Siemens**2

integrationMethod
()¶ :returns:The method to use for calculating the qintegral. :rtype
GaussianBroadening()
TetrahedronMethod

inverseLifeTime
()¶ Returns: The calculated inverse life time with the shape (spin, phonon modes, kpoints, electron bands). Return type: PhysicalQuantity of type inverse time

inverseRelaxationTime
()¶ Returns: The calculated inverse relaxation time with the shape (spin, phonon modes, kpoints, electron bands). Return type: PhysicalQuantity of type inverse time

inverseRelaxationTimeEnergyDependent
()¶ Returns: The calculated or user inputted inverse relaxation time with the shape (n_energies). Return type: PhysicalQuantity of type inverse time

kpoints
()¶ Returns: The unreduced fractional kpoints. The shape is (kpoints, 3). Return type: list of lists of floats

metatext
()¶ Returns: The metatext of the object or None if no metatext is present. Return type: str  unicode  None

method
()¶ Returns: The method used for the calculation. Either Full
orIsotropic
.Return type: NLFlag

mobility
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the mobility. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated mobility. If ‘contribution’ is Hole
orElectron
, the hole or electron mobility is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron mobilities as the first and second element. In case ofSum
, the total mobility is returned. Notice that the summed mobility is not necessarily equal to the sum of Hole and Electron contributions.Return type: PhysicalQuantity with the unit Meter**2 / (Volt * Second)

mobilityTensor
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the mobility tensor. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated mobility tensors. If ‘contribution’ is Hole
orElectron
, the hole or electron mobility tensor is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron mobility tensors as the first and second element. In case ofSum
, the total mobility tensor is returned. Notice that the summed mobility tensor is not necessarily equal to the sum of Hole and Electron contributions.Return type: PhysicalQuantity with the unit Meter**2 / (Volt * Second)

nlprint
(stream=None)¶ Print a string containing an ASCII table useful for plotting the AnalysisSpin object.
Parameters: stream (python stream) – The stream the table should be written to.
Default:NLPrintLogger()

phononEnergies
()¶ Returns: The phonon energies. The shape is (phonon modes in electronphonon coupling, qpoints). Return type: PhysicalQuantity of type energy

phononModes
()¶ Returns: The phonon modes included in the mobility calculation. Return type: list of nonnegative ints

qpoints
()¶ Returns: The unreduced fractional qpoints. The shape is (qpoints, 3). Return type: list of lists of floats

refinement
()¶ Returns: The number of times the qgrid is refined in each direction. Return type: Positive float

resistivityTensor
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the resistivity tensor. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated resistivity tensor. If ‘contribution’ is Hole
orElectron
, the hole or electron resistivity tensor is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron resistivity tensors as the first and second element. In case ofSum
, the total resistivity tensor is returned.Return type: PhysicalQuantity with the unit Meter / Siemens

seebeckCoefficient
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the Seebeck coefficient. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated Seebeck coefficient. If ‘contribution’ is Hole
orElectron
, the hole or electron Seebeck coefficient is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron Seebeck coefficient as the first and second element. In case ofSum
, the total Seebeck coefficient is returned.Return type: PhysicalQuantity with the unit Volt / Kelvin

seebeckCoefficientTensor
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the Seebeck coefficient tensor. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated Seebeck coefficient tensor. If ‘contribution’ is Hole
orElectron
, the hole or electron Seebeck coefficient tensor is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron Seebeck coefficient tensor as the first and second element. In case ofSum
, the total Seebeck coefficient tensor is returned.Return type: PhysicalQuantity with the unit Volt / Kelvin

setMetatext
(metatext)¶ Set a given metatext string on the object.
Parameters: metatext (str  unicode  None) – The metatext string that should be set. A value of “None” can be given to remove the current metatext.

temperature
()¶ Returns: The temperature for the thermal smearing of the Bose and Fermi distributions. Return type: PhysicalQuantity of type temperature

thermalConductivity
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the electronic thermal conductivity. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated electronic thermal conductivity. If ‘contribution’ is Hole
orElectron
, the hole or electron thermal conductivity is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron thermal conductivities as the first and second element. In case ofSum
, the total thermal conductivity is returned.Return type: PhysicalQuantity with the unit J / (Second * Meter * Kelvin)

thermalConductivityTensor
(contribution=None)¶ Parameters: contribution ( All
Sum
Hole
Electron
.) – The contribution to the electronic thermal conductivity tensor. This can be eitherAll
,Sum
,Hole
, orElectron
.
Default:All
Returns: The calculated electronic thermal conductivity tensor. If ‘contribution’ is Hole
orElectron
, the hole or electron thermal conductivity tensor is returned. If ‘contribution’ isAll
, a list is returned with the calculated hole and electron thermal conductivity tensors as the first and second element. In case ofSum
, the total thermal conductivity tensor is returned.Return type: PhysicalQuantity with the unit J / (Second * Meter * Kelvin)

updateTransportCoefficients
(energies=None, inverse_relaxation_time=None, mobility_object=None, temperature=None, fermi_shift=None, energy_broadening=None)¶ Function for updating the transport coefficients (mobility, conductivity, resistivity, Seebeck coefficien, Hall conductivity). This function can only be used if method=Isotropic.
Note that by calling this function you will change variables on the object itself.
Parameters:  energies (PhysicalQuantity of type energy) – List of energies for which to define the scattering rate.
 inverse_relaxation_time (PhysicalQuantity of type energy) – List of scattering rates corrresponding to the provided energies, or a constant
scattering rate. If a list is provided the length must be the same as the
energies
. If amobility_object
is also provided, theinverse_relaxation_time
will be added to the scattering rate obtained from themobility_object
.  mobility_object (
Mobility
) – A mobility object calculated with theFull
method. This can be used to calculate the scattering rate. If ainverse_relaxation_time
is also provided, theinverse_relaxation_time
will be added to the scattering rate obtained from themobility_object
.  temperature (PhysicalQuantity of type temperature) – The temperature to use.
 fermi_shift (PhysicalQuantity of type energy) – The Fermi shift.
 energy_broadening (PhysicalQuantity of type energy.) – The energy broadening of deltafuctions to use
Returns: Dictionary with all the updated transport coefficients.
Return type: dict.

velocitiesK
()¶ Returns: The electron Bloch velocities for each kpoint. The shape is (spin, kpoints, electron bands, velocity components). Return type: PhysicalQuantity with the unit Meter / Second
.
 configuration (
Usage Examples¶
Calculate the mobility of graphene including the full angular dependence of the scattering rates and with the electronphonon coupling and relaxation time calculated from firstprinciples:
# 
# Bulk configuration
# 
lattice = Hexagonal(2.4612*Angstrom, 6.709*Angstrom)
elements = [Carbon, Carbon]
fractional_coordinates = [[0.0 , 0.0 , 0.5],
[0.333333333333, 0.666666666667, 0.5]]
bulk_configuration = BulkConfiguration(
bravais_lattice=lattice,
elements=elements,
fractional_coordinates=fractional_coordinates
)
# 
# Calculator
# 
numerical_accuracy_parameters = NumericalAccuracyParameters(
k_point_sampling=(5, 5, 1),
)
calculator = LCAOCalculator(
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
bulk_configuration.setCalculator(calculator)
bulk_configuration.update()
# 
# Hamiltonian derivatives
# 
hamiltonian_derivatives = HamiltonianDerivatives(
bulk_configuration,
repetitions=(5, 5, 1),
)
# 
# Dynamical matrix
# 
dynamical_matrix = DynamicalMatrix(
bulk_configuration,
repetitions=(5, 5, 1),
max_interaction_range=10*Angstrom,
)
# 
# Electron Phonon Coupling
# 
k_a = numpy.linspace(0.323333, 0.343333, 11)
k_b = numpy.linspace(0.323333, 0.343333, 11)
k_c = numpy.linspace(0, 0, 1)
kpoints = [[a, b, c] for a in k_a for b in k_b for c in k_c]
q_a = numpy.linspace(0.02, 0.02, 11)
q_b = numpy.linspace(0.02, 0.02, 11)
q_c = numpy.linspace(0, 0, 1)
qpoints = [[a, b, c] for a in q_a for b in q_b for c in q_c]
electron_phonon_coupling = ElectronPhononCoupling(
configuration=bulk_configuration,
dynamical_matrix=dynamical_matrix,
hamiltonian_derivatives=hamiltonian_derivatives,
kpoints_fractional=kpoints,
qpoints_fractional=qpoints,
electron_bands=[3, 4],
phonon_modes=All,
energy_tolerance=0.03*eV,
initial_state_energy_range=[100,100]*eV,
store_dense_coupling_matrices=False,
)
# 
# Mobility
# 
mobility = Mobility(
configuration=bulk_configuration,
method=Full,
electron_phonon_coupling=electron_phonon_coupling,
)
nlsave('graphene_method_full.nc', mobility)
Note that the above script only shows the principal work flow for calculating the mobility. For converged calculations, one will often need to use more kpoints and qpoints in the ElectronPhononCoupling object  see the paper [GMSB16] for details. In order to perform a mobility calculation using the Full
method, one needs to input an ElectronPhononCoupling object. While this will result in accurate calculations including the full angular dependece of the scattering rates and relaxation time, it is also a numerically challenging task because the ElectronPhononCoupling object needs to be calculated with many k and qpoints and in addition one needs the DynamicalMatrix and HamiltonianDerivatives objects.
An alternative approach to mobility calculations is available in QuantumATK. By specifying the method keyword to Isotropic
and providing an inverse relaxation time like:
mobility = Mobility(
bulk_configuration,
method=Isotropic,
kpoints=MonkhorstPackGrid(10, 10, 10),
energies=numpy.linspace(0.5, 0.5, 100)*eV,
inverse_relaxation_time=1e10*Second**1)
it is possible to perform much faster calculations, since it is not required to have an ElectronPhononCoupling. Using the Isotropic
method the calculation of both the DynamicalMatrix and HamiltonianDerivatives objects can also be avoided.
Using the Isotropic
method will require knowledge of the inverse relaxation time (or scattering rate). This can either be specified as a constant inverse_relaxation_time=1e10 * Second**1
or as an array with the same length as the provided energies inverse_relaxation_time=numpy.linspace(0, 1e12, 100) * Second**1
. Knowledge of the inverse relaxation time can be obtained from the literature or it can be approximated from a previous Mobility calculation using the Full
method, i.e. calculated from an ElectronPhononCoupling object. This can be obtained as:
# Load a mobility object calculated with the Full method.
mobility = nlread('graphene_method_full.nc', Mobility)[0]
# Define an energy range.
energies=numpy.linspace(0.5, 0.5, 100)*eV
# Calculate the inverse relaxation time.
inverse_relaxation_time = mobility.generateEnergyDependentInverseRelaxationTime(energies=energies)
The obtained inverse relaxation time can be input directly to a new Mobility calculation with the Isotropic
method. Alternatively, the old Mobility object can be provided as input to a new Mobility calculation:
# Load a mobility object calculated with the Full method.
mobility_full = nlread('graphene_method_full.nc', Mobility)[0]
# Perform and new mobility calculation with the Isotropic method.
mobility = Mobility(
bulk_configuration,
method=Isotropic,
energies=numpy.linspace(0.5,0.5,100)*eV,
mobility_object=mobility_full)
If both the inverse_relaxation_time
and mobility_object
are provided, the resulting inverse relaxation time will be the sum of the provided and the one generated from the mobility_object
.
A full script for a Mobility calculation using the Isotropic
method is provided below using the constant inverse relaxation time approximation:
# 
# Bulk configuration
# 
lattice = Hexagonal(2.4612*Angstrom, 6.709*Angstrom)
elements = [Carbon, Carbon]
fractional_coordinates = [[0.0 , 0.0 , 0.5],
[0.333333333333, 0.666666666667, 0.5]]
bulk_configuration = BulkConfiguration(
bravais_lattice=lattice,
elements=elements,
fractional_coordinates=fractional_coordinates
)
# 
# Calculator
# 
numerical_accuracy_parameters = NumericalAccuracyParameters(
k_point_sampling=(5, 5, 1),
)
calculator = LCAOCalculator(
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
bulk_configuration.setCalculator(calculator)
bulk_configuration.update()
# 
# Mobility
# 
kpoint_grid = RegularKpointGrid(
ka_range=[0.32, 0.34],
kb_range=[0.32, 0.34],
kc_range=[0.0, 0.0],
na=15,
nb=15,
nc=1,
)
mobility_constant_relaxation_time = Mobility(
configuration=bulk_configuration,
method=Isotropic,
energies=numpy.linspace(0.05, 0.05, 100)*eV,
inverse_relaxation_time=1e12*Second**1,
electron_bands=[3,4],
kpoints=kpoint_grid,
temperature=300*Kelvin,
)
nlsave('graphene_constant_relaxation_time.nc', mobility_constant_relaxation_time)
mobility_isotropic_scattering_rate.py
The script below shows an example of how to calculate the mobility with an energy dependent inverse relaxation time:
# 
# Bulk configuration
# 
lattice = Hexagonal(2.4612*Angstrom, 6.709*Angstrom)
elements = [Carbon, Carbon]
fractional_coordinates = [[0.0 , 0.0 , 0.5],
[0.333333333333, 0.666666666667, 0.5]]
bulk_configuration = BulkConfiguration(
bravais_lattice=lattice,
elements=elements,
fractional_coordinates=fractional_coordinates
)
# 
# Calculator
# 
numerical_accuracy_parameters = NumericalAccuracyParameters(
k_point_sampling=(5, 5, 1),
)
calculator = LCAOCalculator(
numerical_accuracy_parameters=numerical_accuracy_parameters,
)
bulk_configuration.setCalculator(calculator)
bulk_configuration.update()
# 
# Mobility
# 
kpoint_grid = RegularKpointGrid(
ka_range=[0.32, 0.34],
kb_range=[0.32, 0.34],
kc_range=[0.0, 0.0],
na=15,
nb=15,
nc=1,
)
# Setup energy list.
energies = numpy.linspace(0.05, 0.05, 100)*eV,
# Setup inverse relaxation time as an array with the same length as the energies.
inverse_relaxation_time = abs(numpy.linspace(1e12, 1e+12, 100))*Second**1
mobility_linear_relaxation_time = Mobility(
configuration=bulk_configuration,
method=Isotropic,
energies=energies,
inverse_relaxation_time=inverse_relaxation_time,
electron_bands=[3,4],
kpoints=kpoint_grid,
temperature=300*Kelvin,
)
nlsave('graphene_linear_relaxation_time.nc', mobility_linear_relaxation_time)
Notes¶
The mobility module in QuantumATK enables you to calculate the phononlimited mobility using the semiclassical Boltzmann transport equation (BTE). There are two principle methods to use:
In order to calculate the mobility from firstprinciples including the full angular dependence of the electronphonon scattering, you must use the
Full
method. In that case you need a BulkConfiguration and a ElectronPhononCoupling object. See Ref. [GMSB16] to get a description of the background theory and results obtained with the Mobility object.Assuming an isotropic inverse relaxation time, one can avoid integrations over k,qspace and instead integrate over energy. This leads to much faster mobility calculations since the ElectronPhononCoupling object is no longer needed. Instead you need to provide an inverse relaxation time to the Mobility calculation. The implementation in QuantumATK of the isotropic scattering rate method essentially follows BoltzTrap, Ref. [MS06]. Major differences to the description in Ref. [MS06] are:
 QuantumATK allows for an energy dependent inverse relaxation time.
 The band velocities in QuantumATK are calculated using first order pertubation theory and are thus exact and without any problems at band crossings.
 The energy integrations are performed using the tetrahedron method [BJA94].
[BJA94]  P. E. Blöchl, O. Jepsen, and O. K. Andersen. Improved tetrahedron method for brillouinzone integrations. Phys. Rev. B, 49:16223–16233, Jun 1994. doi:10.1103/PhysRevB.49.16223. 
[GMSB16]  (1, 2) T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge. Firstprinciples method for electronphonon coupling and electron mobility: Applications to twodimensional materials. Phys. Rev. B, 93:035414, Jan 2016. doi:10.1103/PhysRevB.93.035414. 
[MS06]  (1, 2) Georg K.H. Madsen and David J. Singh. Boltztrap. a code for calculating bandstructure dependent quantities. Computer Physics Communications, 175(1):67 – 71, 2006. doi:http://dx.doi.org/10.1016/j.cpc.2006.03.007. 