A Computational Physics Tutorial
Welcome to today’s computational physics blog post! We’re going to dive deep into electronic band structure optimization in crystalline solids. We’ll explore how to find the minimum energy configuration by varying lattice constants and applying periodic boundary conditions.
What is Band Structure Optimization?
In solid-state physics, the electronic band structure describes the range of energy levels that electrons can occupy in a crystalline solid. The total energy of a system depends on:
- Lattice constant ($a$): The spacing between atoms in the crystal
- Periodic boundary conditions: To simulate an infinite crystal
- Electronic structure: How electrons fill available energy bands
The equilibrium lattice constant $a_0$ minimizes the total energy:
$$E_{\text{total}}(a) = E_{\text{kinetic}} + E_{\text{potential}} + E_{\text{electron-electron}}$$
Our Example: 1D Tight-Binding Model
We’ll use the tight-binding approximation for a 1D atomic chain. The Hamiltonian is:
$$H = -t \sum_{\langle i,j \rangle} (c_i^\dagger c_j + c_j^\dagger c_i) + V_{\text{rep}}(a)$$
where:
- $t$ is the hopping parameter (depends on lattice constant)
- $V_{\text{rep}}(a)$ represents ionic repulsion
- $c_i^\dagger, c_i$ are creation/annihilation operators
The dispersion relation for this model is:
$$E(k) = -2t(a) \cos(ka)$$
The Complete Python Code
Let me show you the complete implementation:
1 | import numpy as np |
Detailed Code Explanation
Let me break down the key components of this implementation:
1. The TightBindingChain Class
This is the heart of our simulation. It encapsulates the physics of a 1D atomic chain.
Initialization Parameters:
n_atoms: Number of atoms in our periodic chainn_electrons: Total number of electrons to distributet0: Reference hopping integral at the reference lattice constanta0_ref: Reference lattice spacingalpha: Controls how quickly hopping decreases with distance
2. Hopping Parameter Calculation
1 | def hopping_parameter(self, a): |
The hopping integral $t(a)$ decreases exponentially with increasing lattice constant:
$$t(a) = t_0 \exp\left(-\alpha(a - a_{\text{ref}})\right)$$
This makes physical sense: as atoms move apart, electron wave function overlap decreases, reducing the probability of electron hopping between sites.
3. Building the Hamiltonian Matrix
1 | def build_hamiltonian(self, a): |
This creates a tridiagonal matrix with periodic boundary conditions. The Hamiltonian structure is:
$$H = \begin{pmatrix}
0 & -t & 0 & \cdots & -t \
-t & 0 & -t & \cdots & 0 \
0 & -t & 0 & \cdots & 0 \
\vdots & \vdots & \vdots & \ddots & \vdots \
-t & 0 & 0 & \cdots & 0
\end{pmatrix}$$
The corner elements $H[0, N-1]$ and $H[N-1, 0]$ implement periodic boundary conditions, simulating an infinite crystal.
4. Repulsive Potential Energy
1 | def repulsive_potential(self, a): |
The Born-Mayer potential represents ion-ion repulsion:
$$V_{\text{rep}}(a) = A \exp\left(-\frac{a}{\rho}\right)$$
This prevents the lattice from collapsing. As $a \to 0$, the repulsive energy dominates.
5. Electronic Band Energy Calculation
1 | def calculate_band_energy(self, a): |
This is crucial! We:
- Diagonalize the Hamiltonian to get energy eigenvalues
- Sort them from lowest to highest
- Fill electrons according to the Pauli exclusion principle (factor of 2 for spin)
- Sum occupied state energies
The total electronic energy is:
$$E_{\text{electronic}} = \sum_{i=1}^{N_{\text{occ}}} n_i \epsilon_i$$
where $n_i$ is the occupation number (0, 1, or 2) and $\epsilon_i$ are eigenvalues.
6. Total Energy and Optimization
1 | def total_energy(self, a): |
The total energy combines attractive electronic energy and repulsive ionic energy:
$$E_{\text{total}}(a) = E_{\text{electronic}}(a) + N \cdot V_{\text{rep}}(a)$$
The equilibrium lattice constant $a_0$ is where:
$$\frac{dE_{\text{total}}}{da}\bigg|{a=a_0} = 0 \quad \text{and} \quad \frac{d^2E{\text{total}}}{da^2}\bigg|_{a=a_0} > 0$$
We use scipy.optimize.minimize_scalar to find this minimum numerically.
7. Visualization Strategy
The code generates six comprehensive plots:
- Energy Components vs Lattice Constant: Shows how electronic, repulsive, and total energies vary
- Zoomed Minimum: Clear view of the energy well around equilibrium
- Hopping Parameter Evolution: Shows $t(a)$ decay
- Energy Level Diagram: Visualizes occupied/unoccupied states at $a_0$
- Density of States (DOS): Histogram showing energy level distribution
- Band Structure Comparison: How bands change with different lattice constants
What to Expect in the Results
When you run this code, you’ll see:
- Optimal lattice constant around 3-4 Å (typical for simple metals)
- Clear energy minimum in the total energy curve
- Trade-off between attractive electronic energy (favors small $a$) and repulsive ionic energy (favors large $a$)
- Bandwidth that increases as $a$ decreases (stronger hopping)
- Fermi level marking the boundary between occupied and unoccupied states
The physics here mirrors real materials: the equilibrium structure balances quantum mechanical kinetic energy (electrons prefer delocalization) with electrostatic repulsion (ions resist compression).
Run the Code and Paste Your Results Below!
Now it’s your turn! Copy the code to Google Colab and execute it. The program will:
- Print detailed optimization results
- Generate comprehensive visualization plots
- Provide physical interpretation of the results
====================================================================== ELECTRONIC BAND STRUCTURE OPTIMIZATION 1D Tight-Binding Model with Periodic Boundary Conditions ====================================================================== System Parameters: Number of atoms: 30 Number of electrons: 30 Reference hopping parameter t0: 2.5 eV Reference lattice constant: 3.0 Å Hopping decay parameter α: 1.5 Å⁻¹ Scanning lattice constants from 2.00 to 5.00 Å... Optimizing lattice constant... ====================================================================== OPTIMIZATION RESULTS: ====================================================================== Optimal lattice constant: a₀ = 2.0000 Å Minimum total energy: E_min = -424.9329 eV Energy per atom: -14.1644 eV/atom Band structure at optimal lattice constant: Lowest energy level: -22.4083 eV Highest energy level: 22.4083 eV Bandwidth: 44.8167 eV Fermi energy: -2.3423 eV Figure saved as 'band_structure_optimization.png'
====================================================================== ADDITIONAL ANALYSIS: ====================================================================== Mechanical Properties: Second derivative d²E/da²: -922.2827 eV/Ų (Positive value confirms minimum) HOMO-LUMO gap: 4.6846 eV → System behaves as a SEMICONDUCTOR
Physical Insights
This simple 1D model captures essential physics of real crystals:
- Cohesive energy: The binding energy per atom at equilibrium
- Bulk modulus: The curvature $d^2E/da^2$ relates to material stiffness
- Metallic vs insulating behavior: Determined by whether the Fermi level lies in a band (metal) or gap (insulator)
- Pressure effects: Applying external pressure changes the equilibrium $a_0$
The tight-binding method, despite its simplicity, successfully predicts many properties of real materials and is still widely used in modern computational materials science!
