[Wannier] BoltzWann module for thermoelectric study

[Giovanni Pizzi]

Hi,
I didn’t check the code but I’m ready to bet it’s doing the right thing.
However, before trying to use it, I encourage you think: what would be the density of dopants (number of atoms per unit cell) for a 10^24 cm^-3 doping?

I suggest you compare with the density, e.g., of Si atoms in the unit cell (or the density of atoms in your system).

Then, try to think if this is or not a physically meaningful density of dopants.

Best,
Giovanni


On 27 Jun 2025, at 15:46, Abhijeet Jaysingrao kale ic39253 <ic39253 at imail.iitm.ac.in> wrote:

Dear Prof. Giovanni Pizzi and Wannier90 community,

Kindly let me explain what I am attempting to do. Using the BoltzWann module, I want to reproduce the result computed by BoltzTraP code by varying carrier concentration n (cm-3) at constant temperature T (example<https://pubs.acs.org/cms/10.1021/cm504244b/asset/images/large/cm-2014-04244b_0004.jpeg>). However, in the Wannier90 input file, we can only input the chemical potential μ instead of n (cm-3). So, I need to convert the desired n (cm-3) into μ, for which I have come up with the below python code. Please excuse my coding inaccuracies as I am a beginner at it and kindly fix any mathematical mistakes (possibly in bisection algorithm) and coding errors. If it's all correct, it fails to give converged μ at higher n (e.g. 10^24 cm-3).  I have also attached DOS output computed using quantum espresso code for reference.

Code:

import numpy as np
from scipy import constants

# Constants
T = 300  # Temperature in Kelvin
k_B = constants.k / constants.e  # Boltzmann constant in eV/K (≈8.617333e-5)

# --- Convert unit cell volume to cm³ ---
V_cell_au3 = 1805.5096  # Unit cell volume in (a.u.)³
bohr_to_angstrom = 0.52917721
V_cell_A3 = V_cell_au3 * (bohr_to_angstrom ** 3)
V_cell_cm3 = V_cell_A3 * 1e-24

print(f"Unit cell volume = {V_cell_cm3:.4e} cm³")

# Fermi-Dirac function
def fermi_dirac(E, mu, T):
    kT = k_B * T
    x = (E - mu) / kT

    # Handle large values to prevent overflow
    return np.piecewise(x,
                       [x < -100, x > 100],
                       [1.0, 0.0, lambda x: 1/(1 + np.exp(x))])

# Load DOS data
data = np.loadtxt('case_dos.dat')
energy = data[:, 0]  # Energy in eV
dos = data[:, 1]     # DOS in states/eV/unit cell

# Compute carrier concentration using trapezoid rule
def carrier_concentration(mu, energy, dos, T, V_cm3):
    product = dos * fermi_dirac(energy, mu, T)
    n_unit_cell = np.trapezoid(product, energy)
    n_cm3 = n_unit_cell / V_cm3
    return n_cm3

# Bisection algorithm with auto-bound detection
def find_mu(target_n, energy, dos, T, V_cm3, tol=1e-3, max_iter=100):
    # Auto-detect proper bounds
    mu_low, mu_high = -100, 100  # Start with wide range
    n_low = carrier_concentration(mu_low, energy, dos, T, V_cm3)
    n_high = carrier_concentration(mu_high, energy, dos, T, V_cm3)

    # adjust lower bound
    while n_low > target_n:
        mu_low -= 50
        n_low = carrier_concentration(mu_low, energy, dos, T, V_cm3)
        print(f"Expanding lower bound to μ = {mu_low} eV, n = {n_low:.4e} cm⁻³")

    # adjustupper bound
    while n_high < target_n:
        mu_high += 50
        n_high = carrier_concentration(mu_high, energy, dos, T, V_cm3)
        print(f"Expanding upper bound to μ = {mu_high} eV, n = {n_high:.4e} cm⁻³")

    print(f"\nStarting bisection with bounds: μ = [{mu_low}, {mu_high}] eV")
    print(f"Initial concentrations: n_low = {n_low:.4e} cm⁻³, n_high = {n_high:.4e} cm⁻³")

    # Bisection loop
    for i in range(max_iter):
        mu_mid = (mu_low + mu_high) / 2
        n_mid = carrier_concentration(mu_mid, energy, dos, T, V_cm3)

        rel_error = abs(n_mid - target_n) / target_n
        print(f"Iter {i+1}: μ = {mu_mid:.6f} eV, n = {n_mid:.4e} cm⁻³, error = {rel_error:.2%}")

        if rel_error < tol:
            print("\nConverged!")
            return mu_mid

        if n_mid < target_n:
            mu_low = mu_mid
        else:
            mu_high = mu_mid

    print(f"Warning: Max iterations reached ({max_iter})")
    return mu_mid

if __name__ == "__main__":
    # Target carrier concentration (cm⁻³)
    target_n = 1e18

    # Find chemical potential
    print(f"\nSearching for μ for n = {target_n:.1e} cm⁻³")
    mu_target = find_mu(target_n, energy, dos, T, V_cell_cm3)

    # Final verification
    n_final = carrier_concentration(mu_target, energy, dos, T, V_cell_cm3)
    print(f"\nResult: μ = {mu_target:.6f} eV gives n = {n_final:.4e} cm⁻³")

On Mon, Jun 23, 2025 at 4:02 PM Abhijeet Jaysingrao kale ic39253 <ic39253 at imail.iitm.ac.in<mailto:ic39253 at imail.iitm.ac.in>> wrote:
Dear Prof. Giovanni Pizzi,

Thank you very much for the information. I tried obtaining a Fermi-Dirac distribution but did not succeed. Can I get further help regarding how to obtain Fermi-Dirac distribution function from Wannier output and then finally compute carrier concentration?

Regards,
Abhijeet.

On Thu, Jun 19, 2025 at 8:26 PM Giovanni Pizzi <giovanni.pizzi at epfl.ch<mailto:giovanni.pizzi at epfl.ch>> wrote:
Hi,

1. You should convert the carrier concentration into a chemical potential by integrating the product of the Fermi-Dirac distribution function and of the DOS, and looking for the chemical potential that gives you the expected number of electrons (e.g. with a bisection algorithm).

2. Internally it should work, but indeed the code only prints the upper triangle for sigma, Sigma*S and kappa (electronic contribution) to thermal conductivity. If you want to check if there is a non-symmetric part, you need to adapt the code to print all coordinates.
(Note that for the Seebeck coefficient, all 9 components are reported instead).

Best,
Giovanni Pizzi

Group leader, Materials Software and Data Group, PSI
https://www.psi.ch/en/lms/people/giovanni-pizzi




On 18 Jun 2025, at 12:23, Abhijeet Jaysingrao kale ic39253 <ic39253 at imail.iitm.ac.in<mailto:ic39253 at imail.iitm.ac.in>> wrote:

Dear Wannier90 community,

I am Abhijeet from IIT Madras (India). I am using the BoltzWann module to compute thermoelectric properties of monolayer systems. I'd be grateful if anyone could suggest me on following questions:

1. How to input carrier concentration range like chemical potential and temperature?
2. Can we use BoltzWann for thermoelectric properties of magnetic materials even though it outputs 6 components of electrical conductivity when it requires all 9?

Thanks in advance.
Regards,
Abhijeet.

_______________________________________________
Wannier mailing list
Wannier at lists.quantum-espresso.org<mailto:Wannier at lists.quantum-espresso.org>
https://lists.quantum-espresso.org/mailman/listinfo/wannier

<case_dos.dat>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.quantum-espresso.org/pipermail/wannier/attachments/20250627/8eb8a268/attachment.html>