Chapter 3: Thermal Transport Calculations

3.1 Introduction: Beyond the Relaxation Time Approximation

In intermediate-level treatments, phonon thermal conductivity is often calculated using the relaxation time approximation (RTA), which assumes that each phonon mode relaxes independently to equilibrium with a single relaxation time \(\tau_{\mathbf{k}s}\). While RTA provides reasonable estimates for many materials, it has fundamental limitations:

• Momentum conservation in normal (N) processes is ignored
• Collective drift of phonon distributions cannot be captured
• Quantitative accuracy is often insufficient (errors 20-50%)
• Low-temperature behavior in high-purity crystals is incorrect

This chapter presents the full phonon Boltzmann transport equation (PBTE) and its iterative solution methods, which have become the gold standard for thermal transport calculations. We will also explore modern computational tools (ShengBTE, almaBTE, Phono3py) and emerging machine learning approaches.

3.2 The Phonon Boltzmann Transport Equation

3.2.1 General Form of the PBTE

The phonon distribution function \(n_{\mathbf{k}s}(\mathbf{r}, t)\) in non-equilibrium conditions evolves according to the Boltzmann equation:

For steady-state thermal transport with no external forces:

\[ \mathbf{v}_{\mathbf{k}s} \cdot \nabla_{\mathbf{r}} n_{\mathbf{k}s} = \left(\frac{\partial n_{\mathbf{k}s}}{\partial t}\right)_{\text{scatt}} \]

3.2.2 Linearization and Deviation Function

For small temperature gradients, we expand the distribution around equilibrium:

\[ n_{\mathbf{k}s} = n_{\mathbf{k}s}^0(T) + \delta n_{\mathbf{k}s} \]

where \(n_{\mathbf{k}s}^0 = 1/(e^{\hbar\omega_{\mathbf{k}s}/k_BT} - 1)\) is the Bose-Einstein distribution. The temperature gradient creates a deviation:

\[ \mathbf{v}_{\mathbf{k}s} \cdot \nabla T \frac{\partial n^0}{\partial T} = \left(\frac{\partial \delta n_{\mathbf{k}s}}{\partial t}\right)_{\text{scatt}} \]

It is convenient to define a deviation function \(\Phi_{\mathbf{k}s}\):

3.2.3 The Scattering Collision Integral

The collision integral includes all scattering processes. For three-phonon interactions, where phonon \(\mathbf{k}s\) can either absorb a phonon (\(\mathbf{k}'s' \to \mathbf{k}s + \mathbf{k}''s''\)) or emit a phonon (\(\mathbf{k}s \to \mathbf{k}'s' + \mathbf{k}''s''\)):

3.2.4 Relaxation Time Approximation (RTA)

In RTA, the collision integral is approximated as:

\[ \left(\frac{\partial n_{\mathbf{k}s}}{\partial t}\right)_{\text{scatt}} \approx -\frac{\delta n_{\mathbf{k}s}}{\tau_{\mathbf{k}s}} \]

This leads to the simple solution:

\[ \Phi_{\mathbf{k}s}^{\text{RTA}} = \tau_{\mathbf{k}s} \mathbf{v}_{\mathbf{k}s} \cdot \nabla T \]

The thermal conductivity in RTA is:

3.2.5 Iterative Solution of the Full BTE

The full BTE without approximations is a complex linear system. The deviation function satisfies:

\[ \Phi_{\mathbf{k}s} = \Phi_{\mathbf{k}s}^{\text{RTA}} + \Delta\Phi_{\mathbf{k}s} \]

where the correction \(\Delta\Phi\) accounts for collective effects. This is solved iteratively:

3.3 Three-Phonon Scattering from First Principles

3.3.1 Anharmonic Force Constants

The strength of three-phonon interactions is determined by the third-order anharmonic force constants:

\[ \Phi_{\alpha\beta\gamma}(lκ, l'κ', l''κ'') = \frac{\partial^3 E}{\partial u_{\alpha}(lκ) \partial u_{\beta}(l'κ') \partial u_{\gamma}(l''κ'')} \]

where \(u_\alpha(lκ)\) is the \(\alpha\)-component displacement of atom \(κ\) in unit cell \(l\). These can be calculated using:

• Finite differences: Compute forces for displaced atomic configurations
• DFPT (Density Functional Perturbation Theory): Direct calculation
• Machine learning potentials: Fit from MD trajectories

3.3.2 Three-Phonon Interaction Matrix Elements

The interaction strength \(V^{(3)}\) in Fourier space is:

3.3.3 Scattering Rates and Phase Space

The scattering rate for a phonon \(\mathbf{k}s\) is obtained by integrating the collision integral. For emission processes (\(\mathbf{k}s \to \mathbf{k}'s' + \mathbf{k}''s''\)):

\[ \Gamma^{\text{em}}_{\mathbf{k}s} = \frac{\pi}{\hbar} \sum_{\mathbf{k}'s',\mathbf{k}''s''} |V^{(3)}|^2 (n_{\mathbf{k}'s'} + 1)(n_{\mathbf{k}''s''} + 1) \delta(\omega - \omega' - \omega'') \delta_{\mathbf{k},\mathbf{k}'+\mathbf{k}''+\mathbf{G}} \]

The delta functions impose strict conservation laws, restricting the phase space for scattering. The total scattering rate is:

\[ \Gamma_{\mathbf{k}s} = \Gamma^{\text{em}}_{\mathbf{k}s} + \Gamma^{\text{abs}}_{\mathbf{k}s} \]

3.4 Four-Phonon Scattering

3.4.1 When Are Four-Phonon Processes Important?

Traditionally, three-phonon processes were thought to dominate thermal resistance. However, recent work has shown that four-phonon scattering can be crucial in certain materials:

• High temperatures: When \(k_BT \sim \hbar\omega\), four-phonon rates scale as \(T^2\)
• Weak three-phonon scattering: Materials with high symmetry or weak anharmonicity
• Optical phonons: Four-phonon decay channels for high-frequency modes
• 2D materials: Graphene, h-BN where four-phonon effects are enhanced

3.4.2 Four-Phonon Interaction Matrix Elements

Four-phonon processes involve fourth-order force constants:

\[ \Phi_{\alpha\beta\gamma\delta}(lκ, l'κ', l''κ'', l'''κ''') = \frac{\partial^4 E}{\partial u_{\alpha}(lκ) \partial u_{\beta}(l'κ') \partial u_{\gamma}(l''κ'') \partial u_{\delta}(l'''κ''')} \]

The scattering rate for a process \(\mathbf{k}_1s_1 + \mathbf{k}_2s_2 \to \mathbf{k}_3s_3 + \mathbf{k}_4s_4\) is:

\[ \Gamma^{(4)} \propto |V^{(4)}|^2 \delta(\omega_1 + \omega_2 - \omega_3 - \omega_4) \delta(\mathbf{k}_1 + \mathbf{k}_2 - \mathbf{k}_3 - \mathbf{k}_4 - \mathbf{G}) \]

3.4.3 Computational Cost Scaling

The computational cost scales dramatically with scattering order:

Efficient algorithms and high-performance computing are essential for four-phonon calculations.

3.5 Computational Tools for Thermal Transport

3.5.1 ShengBTE: Theory and Workflow

ShengBTE is a widely-used code for solving the phonon BTE iteratively. Developed by Wu Li, it implements the variational approach to the BTE.

&allocations
 nelements=1,
 natoms=2,
 ngrid(:)=25 25 25,
 norientations=0,
/

&crystal
 lfactor=1.0,
 lattvec(:,1)=0.0 2.715 2.715,
 lattvec(:,2)=2.715 0.0 2.715,
 lattvec(:,3)=2.715 2.715 0.0,
 types(1)='Si',
 elements(1)='Si',
 gfactors(1)=1,
 epsilon(1)=0.0,
 born(:,:,1)=0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0,
 scell(:)=5 5 5,
 orientations(:,1)=1 0 0,
 orientations(:,2)=0 1 0,
 orientations(:,3)=0 0 1,
 positions(:,1)=0.0 0.0 0.0,
 positions(:,2)=0.25 0.25 0.25,
 masses(1)=28.0855,
/

¶meters
 T=300.0,
 scalebroad=0.5,
 convergence=1e-5,
 isotopes=.TRUE.,
 nanowires=.FALSE.,
 onlyharmonic=.FALSE.,
/

&flags
 nonanalytic=.FALSE.,
 isotopes=.TRUE.,
 autoisotopes=.TRUE.,
 espresso=.FALSE.,
/

3.5.2 almaBTE: Advanced Features

almaBTE is a modern C++ code offering advanced features:

• Four-phonon scattering: Built-in support for fourth-order processes
• Nanostructures: Suppression functions for thin films, nanowires
• Cumulative analysis: Thermal conductivity vs MFP, frequency cutoffs
• Direct BTE solution: Matrix inversion for small systems
• HDF5 output: Efficient data storage and post-processing

3.5.3 Phono3py: Integration with Phonopy

Phono3py seamlessly integrates with Phonopy for phonon calculations:

• Automated displacement generation for \(\Phi^{(3)}\)
• Compatibility with VASP, QE, LAMMPS, and more
• Built-in RTA thermal conductivity solver
• Spectral thermal conductivity analysis

# 1. Generate displacements for 3rd-order FCs
phono3py --dim="2 2 2" -c POSCAR-unitcell --sym-fc3

# 2. Calculate forces with DFT
# (run VASP or QE for each POSCAR-xxxxx)

# 3. Create FORCES_FC3 file
phono3py --cf3 vasprun.xml-{00001..00110}

# 4. Calculate thermal conductivity
phono3py --mesh="19 19 19" --br --isotope --ts="200 300 400 500 600"

# Output: kappa-mxxx.hdf5 (thermal conductivity data)

3.6 Convergence and Numerical Accuracy

3.6.1 Critical Convergence Parameters

Achieving converged thermal conductivity requires careful attention to:

3.6.2 Common Convergence Issues

3.6.3 Benchmark: Silicon Thermal Conductivity

import numpy as np
import matplotlib.pyplot as plt

# Experimental value
kappa_exp = 145  # W/m·K

# q-grid convergence
q_grids = [15, 20, 25, 30, 35, 40]
kappa_q = [189, 168, 158, 152, 149, 147]  # RTA values

# Smearing convergence (30×30×30 grid)
smearing = [0.1, 0.3, 0.5, 0.7, 1.0, 1.5]  # THz
kappa_smear = [142, 149, 151, 152, 152, 151]

# Iterative BTE convergence (30×30×30, smearing=0.5)
iterations = [0, 5, 10, 20, 50, 100]
kappa_iter = [151, 163, 172, 180, 185, 186]  # 0 = RTA

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# q-grid
axes[0].plot(q_grids, kappa_q, 'o-', label='Calculated')
axes[0].axhline(kappa_exp, color='r', linestyle='--', label='Experiment')
axes[0].set_xlabel('q-grid size')
axes[0].set_ylabel('κ (W/m·K)')
axes[0].set_title('q-grid Convergence')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Smearing
axes[1].plot(smearing, kappa_smear, 's-', label='Calculated')
axes[1].axhline(kappa_exp, color='r', linestyle='--', label='Experiment')
axes[1].set_xlabel('Gaussian smearing (THz)')
axes[1].set_ylabel('κ (W/m·K)')
axes[1].set_title('Smearing Parameter')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Iterations
axes[2].plot(iterations, kappa_iter, '^-', label='Iterative BTE')
axes[2].axhline(kappa_exp, color='r', linestyle='--', label='Experiment')
axes[2].set_xlabel('Iterations')
axes[2].set_ylabel('κ (W/m·K)')
axes[2].set_title('Iterative Convergence')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('si_convergence.png', dpi=300, bbox_inches='tight')
plt.show()

# Final result with 4-phonon
print("Silicon at 300 K:")
print(f"  RTA (3-phonon):            {kappa_iter[0]:.1f} W/m·K")
print(f"  Iterative BTE (3-phonon):  {kappa_iter[-1]:.1f} W/m·K")
print(f"  With 4-phonon:             150 W/m·K")
print(f"  Experiment:                {kappa_exp} W/m·K")

3.7 Mean Free Path Analysis

3.7.1 Phonon Mean Free Path Distribution

The phonon mean free path (MFP) is:

\[ \Lambda_{\mathbf{k}s} = v_{\mathbf{k}s} \tau_{\mathbf{k}s} \]

Different phonon modes contribute differently to thermal conductivity based on their MFP. The cumulative thermal conductivity is:

import numpy as np
import matplotlib.pyplot as plt

def calculate_mfp_distribution(omega, v, tau, kappa_mode):
    """
    Calculate cumulative thermal conductivity vs MFP

    Parameters:
    -----------
    omega : array, phonon frequencies (THz)
    v : array, group velocities (m/s)
    tau : array, relaxation times (ps)
    kappa_mode : array, mode contributions to κ (W/m·K)

    Returns:
    --------
    mfp_bins : array, MFP values (nm)
    kappa_cumulative : array, cumulative κ (W/m·K)
    """
    # Calculate MFP for each mode
    mfp = v * tau * 1e-3  # Convert ps·m/s to nm

    # Sort by MFP
    sort_idx = np.argsort(mfp)
    mfp_sorted = mfp[sort_idx]
    kappa_sorted = kappa_mode[sort_idx]

    # Cumulative sum
    kappa_cumulative = np.cumsum(kappa_sorted)

    return mfp_sorted, kappa_cumulative

# Example: Silicon at 300 K
n_modes = 10000
np.random.seed(42)

# Phonon properties (simplified model)
omega = np.random.lognormal(np.log(10), 0.5, n_modes)  # THz
v = 3000 + 2000*np.random.randn(n_modes)  # m/s
tau = 10 / omega**2  # ps, ω² scattering
tau = np.clip(tau, 0.1, 100)

# Mode contribution (proportional to C·v²·τ)
C_mode = omega**2 / (np.exp(omega/6.25) - 1)**2  # Simplified heat capacity at 300K
kappa_mode = C_mode * v**2 * tau
kappa_mode = kappa_mode / np.sum(kappa_mode) * 150  # Normalize to total κ = 150 W/m·K

# Calculate distribution
mfp, kappa_cum = calculate_mfp_distribution(omega, v, tau, kappa_mode)

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# Cumulative κ vs MFP
ax1.semilogx(mfp, kappa_cum, linewidth=2)
ax1.axhline(150, color='r', linestyle='--', alpha=0.5, label='Total κ')
ax1.axvline(100, color='g', linestyle='--', alpha=0.5, label='100 nm')
ax1.set_xlabel('Mean Free Path (nm)')
ax1.set_ylabel('Cumulative κ (W/m·K)')
ax1.set_title('Cumulative Thermal Conductivity')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Differential contribution
mfp_bins = np.logspace(-1, 4, 50)
kappa_diff, _ = np.histogram(mfp, bins=mfp_bins, weights=kappa_mode)
mfp_centers = np.sqrt(mfp_bins[:-1] * mfp_bins[1:])

ax2.loglog(mfp_centers, kappa_diff, 'o-', linewidth=2)
ax2.set_xlabel('Mean Free Path (nm)')
ax2.set_ylabel('dκ/d(log Λ) (W/m·K)')
ax2.set_title('Differential MFP Contribution')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('mfp_distribution.png', dpi=300, bbox_inches='tight')
plt.show()

# Analysis
mfp_50 = mfp[np.searchsorted(kappa_cum, 75)]  # MFP for 50% of κ
mfp_90 = mfp[np.searchsorted(kappa_cum, 135)]  # MFP for 90% of κ

print(f"MFP Analysis for Si at 300 K:")
print(f"  50% of κ from phonons with MFP < {mfp_50:.1f} nm")
print(f"  90% of κ from phonons with MFP < {mfp_90:.1f} nm")
print(f"  Max MFP: {mfp[-1]:.0f} nm")

3.7.2 Coherent vs Incoherent Transport

The Boltzmann transport equation assumes incoherent (particle-like) phonon transport, where phonons scatter randomly and phase information is lost. However, at small length scales or low temperatures, coherent (wave-like) effects can emerge:

3.8 Machine Learning for Thermal Conductivity

3.8.1 Why Machine Learning?

First-principles BTE calculations are computationally expensive, limiting large-scale screening. Machine learning offers a path to rapid prediction:

• Speed: Predictions in seconds vs days of compute time
• Scalability: Screen millions of materials
• Insights: Identify key descriptors controlling \(\kappa\)
• Design: Inverse design for target thermal properties

3.8.2 Descriptors for Thermal Conductivity

Effective ML models require physically meaningful descriptors:

3.8.3 ML Workflow for Thermal Conductivity

import numpy as np
import pandas as pd
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.metrics import mean_absolute_error, r2_score
import matplotlib.pyplot as plt

# Example dataset (simplified)
# In reality, use databases like Materials Project, JARVIS
data = {
    'material': ['Si', 'Ge', 'Diamond', 'c-BN', 'GaN', 'AlN', 'GaAs', 'PbTe'],
    'mass_avg': [28.1, 72.6, 12.0, 12.4, 41.9, 20.5, 72.3, 167.0],  # amu
    'density': [2.33, 5.32, 3.51, 3.45, 6.15, 3.26, 5.32, 8.24],  # g/cm³
    'debye_temp': [645, 374, 2230, 1900, 600, 950, 344, 136],  # K
    'sound_vel': [6400, 3900, 17500, 14800, 6500, 10500, 4730, 1900],  # m/s
    'gruneisen': [0.99, 1.0, 1.05, 1.2, 1.3, 0.8, 0.87, 1.5],
    'kappa_300K': [145, 60, 2000, 1300, 130, 285, 55, 2.5]  # W/m·K (target)
}

df = pd.DataFrame(data)

# Features and target
features = ['mass_avg', 'density', 'debye_temp', 'sound_vel', 'gruneisen']
X = df[features].values
y = np.log10(df['kappa_300K'].values)  # Log scale for large κ range

# Train-test split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# Train model
model = GradientBoostingRegressor(
    n_estimators=100,
    learning_rate=0.1,
    max_depth=3,
    random_state=42
)

model.fit(X_train, y_train)

# Predictions
y_pred_train = model.predict(X_train)
y_pred_test = model.predict(X_test)

# Convert back from log scale
kappa_pred_train = 10**y_pred_train
kappa_pred_test = 10**y_pred_test
kappa_true_train = 10**y_train
kappa_true_test = 10**y_test

# Metrics
mae_test = mean_absolute_error(kappa_true_test, kappa_pred_test)
r2_test = r2_score(kappa_true_test, kappa_pred_test)

print("Model Performance:")
print(f"  Test MAE: {mae_test:.1f} W/m·K")
print(f"  Test R²:  {r2_test:.3f}")
print()

# Feature importance
importance = model.feature_importances_
for feat, imp in zip(features, importance):
    print(f"  {feat:15s}: {imp:.3f}")

# Visualize
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(kappa_true_train, kappa_pred_train, s=100, alpha=0.6, label='Train')
ax.scatter(kappa_true_test, kappa_pred_test, s=100, alpha=0.6, label='Test')
ax.plot([1, 2000], [1, 2000], 'k--', alpha=0.3)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('True κ (W/m·K)')
ax.set_ylabel('Predicted κ (W/m·K)')
ax.set_title('ML Thermal Conductivity Prediction')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ml_kappa_prediction.png', dpi=300, bbox_inches='tight')
plt.show()

# Predict new material
new_material = np.array([[50, 4.0, 800, 7000, 1.1]])  # Example
kappa_new = 10**model.predict(new_material)[0]
print(f"\nPrediction for new material: κ ≈ {kappa_new:.1f} W/m·K")

3.8.4 Advanced ML Approaches

3.9 High-Throughput Thermal Conductivity Screening

3.9.1 Workflow for Large-Scale Screening

3.9.2 Computational Cost Reduction Strategies

• Pre-screening with ML: Narrow candidates before DFT
• Coarse grids first: Use 15×15×15, refine only for promising materials
• RTA as filter: Full BTE only for \(\kappa^{\text{RTA}} > \text{threshold}\)
• Symmetry exploitation: Reduce irreducible \(\mathbf{q}\)-points
• Parallel computing: Distribute \(\mathbf{q}\)-points across nodes

3.9.3 Case Study: Screening for Low-κ Thermoelectrics

import numpy as np
import matplotlib.pyplot as plt

# Simplified screening example
materials = [
    'PbTe', 'SnTe', 'GeTe', 'PbSe', 'SnSe',
    'Bi2Te3', 'Sb2Te3', 'CoSb3', 'Mg2Si', 'Mg2Sn'
]

# Properties from literature / databases
kappa_300K = np.array([2.5, 3.8, 2.0, 2.0, 0.7, 1.5, 1.8, 4.0, 7.5, 8.0])
seebeck = np.array([170, 150, 180, 160, 230, 200, 180, 100, 300, 280])  # μV/K
sigma = np.array([2e4, 1.5e4, 3e4, 2.5e4, 0.5e4, 1e4, 1.2e4, 5e4, 0.8e4, 1e4])  # S/m

# Thermoelectric figure of merit: ZT = S²σT/κ
T = 300  # K
ZT = (seebeck * 1e-6)**2 * sigma * T / kappa_300K

# Power factor
PF = (seebeck**2) * sigma / 1e3  # mW/m·K²

# Screening criteria
target_kappa_max = 3.0  # W/m·K
target_ZT_min = 0.3

# Filter
mask = (kappa_300K < target_kappa_max) & (ZT > target_ZT_min)
candidates = [m for m, keep in zip(materials, mask) if keep]

print("Thermoelectric Screening Results:")
print("=" * 60)
for mat, k, s, sig, zt in zip(materials, kappa_300K, seebeck, sigma, ZT):
    status = "✓ PASS" if (k < target_kappa_max and zt > target_ZT_min) else "✗ Fail"
    print(f"{mat:10s} κ={k:.2f} W/m·K  S={s:3.0f} μV/K  ZT={zt:.3f}  {status}")

print(f"\nCandidates for detailed study: {', '.join(candidates)}")

# Visualize
fig, ax = plt.subplots(figsize=(8, 6))
colors = ['green' if m else 'gray' for m in mask]
scatter = ax.scatter(kappa_300K, ZT, s=200, c=colors, alpha=0.7, edgecolors='black')

for i, mat in enumerate(materials):
    ax.annotate(mat, (kappa_300K[i], ZT[i]),
                xytext=(5, 5), textcoords='offset points', fontsize=9)

ax.axvline(target_kappa_max, color='r', linestyle='--', label=f'κ < {target_kappa_max} W/m·K')
ax.axhline(target_ZT_min, color='b', linestyle='--', label=f'ZT > {target_ZT_min}')
ax.set_xlabel('Thermal Conductivity κ (W/m·K)', fontsize=12)
ax.set_ylabel('Figure of Merit ZT', fontsize=12)
ax.set_title('Thermoelectric Materials Screening', fontsize=14)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('thermoelectric_screening.png', dpi=300, bbox_inches='tight')
plt.show()

3.10 Practical Examples

3.10.1 Silicon: Benchmark Calculation

Silicon is the standard benchmark for thermal transport codes. Let's walk through a complete calculation workflow.

"""
Silicon thermal conductivity calculation workflow
Requires: VASP, phonopy, phono3py (or ShengBTE)
"""

# Step 1: DFT structure optimization
# INCAR: PREC=Accurate, ENCUT=500, EDIFF=1E-8
# Optimize until forces < 0.001 eV/Å

# Step 2: Harmonic phonons with phonopy
# phonopy --dim="4 4 4" -c POSCAR-unitcell
# (compute forces for displaced structures)
# phonopy --fc vasprun.xml-{001..002}
# phonopy --mesh="30 30 30" -p  # Check no imaginary frequencies

# Step 3: Third-order force constants
# phono3py --dim="2 2 2" --dim-fc2="4 4 4" -c POSCAR-unitcell
# (compute forces for ~100 displaced structures)
# phono3py --cf3 vasprun.xml-{00001..00110}

# Step 4: Calculate thermal conductivity
# phono3py --mesh="19 19 19" --br --isotope --ts="100 200 300 400 500 600"

# Output analysis
import h5py
import numpy as np
import matplotlib.pyplot as plt

# Read results
with h5py.File('kappa-m191919.hdf5', 'r') as f:
    temperatures = f['temperature'][:]
    kappa = f['kappa'][:, :, :]  # (T, 3x3 tensor)
    kappa_iso = np.mean(kappa[:, [0,1,2], [0,1,2]], axis=1)  # Average diagonal

# Experimental data (from literature)
T_exp = np.array([100, 200, 300, 400, 500, 600])
kappa_exp = np.array([900, 280, 145, 90, 65, 50])

# Plot
plt.figure(figsize=(8, 6))
plt.plot(temperatures, kappa_iso, 'o-', linewidth=2, markersize=8, label='Phono3py (RTA)')
plt.plot(T_exp, kappa_exp, 's--', linewidth=2, markersize=8, label='Experiment')
plt.xlabel('Temperature (K)', fontsize=12)
plt.ylabel('Thermal Conductivity (W/m·K)', fontsize=12)
plt.title('Silicon Thermal Conductivity', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('si_kappa_vs_T.png', dpi=300, bbox_inches='tight')
plt.show()

# Analyze MFP distribution at 300 K
with h5py.File('kappa-m191919.hdf5', 'r') as f:
    mfp = f['mean_free_path'][2, :, :]  # T=300K, all modes
    gamma = f['gamma'][2, :, :]  # Scattering rates

print("Silicon at 300 K:")
print(f"  Calculated κ: {kappa_iso[2]:.1f} W/m·K")
print(f"  Experimental: {kappa_exp[2]} W/m·K")
print(f"  Mean MFP: {np.mean(mfp):.1f} nm")
print(f"  Median MFP: {np.median(mfp):.1f} nm")

3.10.2 PbTe: Thermoelectric Material

Lead telluride (PbTe) is a classic thermoelectric material with intrinsically low thermal conductivity due to heavy atoms and strong anharmonicity.

import numpy as np
import matplotlib.pyplot as plt

# Simulated data for PbTe (based on literature)
T = np.array([300, 400, 500, 600, 700, 800])

# Thermal conductivity components
kappa_lattice = np.array([2.5, 2.3, 2.1, 2.0, 1.9, 1.8])  # W/m·K, phonon contribution
kappa_electron = 0.3 * np.ones_like(T)  # Electronic contribution (doped)
kappa_total = kappa_lattice + kappa_electron

# Key features of PbTe
print("PbTe Thermal Transport Characteristics:")
print("=" * 60)
print("  Crystal structure: Rock salt (Fm-3m)")
print("  Lattice parameter: a = 6.46 Å")
print("  Average atomic mass: 167 amu (very heavy)")
print("  Grüneisen parameter: γ ≈ 1.5 (strong anharmonicity)")
print("  Debye temperature: θ_D ≈ 136 K (soft phonons)")
print()
print("Why is κ_lattice so low?")
print("  1. Heavy atoms → low phonon velocities")
print("  2. Strong anharmonicity → short phonon lifetimes")
print("  3. Avoided crossing in acoustic-optical branches")
print("  4. Soft TO phonon mode → strong scattering")
print()

# Temperature dependence
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# κ vs T
ax1.plot(T, kappa_lattice, 's-', linewidth=2, label='Lattice')
ax1.plot(T, kappa_electron, 'o-', linewidth=2, label='Electronic')
ax1.plot(T, kappa_total, '^-', linewidth=2, label='Total')
ax1.set_xlabel('Temperature (K)')
ax1.set_ylabel('Thermal Conductivity (W/m·K)')
ax1.set_title('PbTe Thermal Conductivity')
ax1.legend()
ax1.grid(True, alpha=0.3)

# κ^-1 vs T (approximately linear at high T)
ax2.plot(T, 1/kappa_lattice, 's-', linewidth=2)
ax2.set_xlabel('Temperature (K)')
ax2.set_ylabel('1/κ_lattice (m·K/W)')
ax2.set_title('Inverse Thermal Conductivity\n(High-T Umklapp behavior)')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('pbte_thermal_conductivity.png', dpi=300, bbox_inches='tight')
plt.show()

# Compare to Si
print("\nComparison: PbTe vs Silicon at 300 K")
print("-" * 60)
print(f"{'Property':<30s} {'PbTe':<15s} {'Si':<15s}")
print("-" * 60)
print(f"{'κ_lattice (W/m·K)':<30s} {kappa_lattice[0]:<15.1f} {145:<15.1f}")
print(f"{'Average mass (amu)':<30s} {167:<15.1f} {28.1:<15.1f}")
print(f"{'Debye temperature (K)':<30s} {136:<15.1f} {645:<15.1f}")
print(f"{'Grüneisen parameter':<30s} {1.5:<15.1f} {0.99:<15.1f}")
print(f"{'κ_PbTe / κ_Si':<30s} {kappa_lattice[0]/145:<15.3f} {'—':<15s}")
print("-" * 60)

3.11 Summary and Best Practices

Exercises

Further Reading