Comprehensive Python library for phonon calculations that enables lattice dynamics simulations and vibrational property analysis of crystalline materials.
—
Calculate mode Grüneisen parameters that quantify anharmonicity by measuring how phonon frequencies change with volume variations. Grüneisen parameters are crucial for understanding thermal expansion, pressure effects, and anharmonic properties of materials.
Initialize and perform Grüneisen parameter calculations using three Phonopy instances at different volumes.
class PhonopyGruneisen:
def __init__(
self,
phonon: Phonopy,
phonon_plus: Phonopy,
phonon_minus: Phonopy,
delta_strain: float = None
):
"""
Initialize Grüneisen parameter calculation.
Parameters:
- phonon: Reference Phonopy instance at equilibrium volume
- phonon_plus: Phonopy instance at expanded volume (V + ΔV)
- phonon_minus: Phonopy instance at compressed volume (V - ΔV)
- delta_strain: Volume strain used for finite difference (optional)
"""
def get_phonon(self) -> Phonopy:
"""Get reference phonon object."""Example:
from phonopy import Phonopy, PhonopyGruneisen
from phonopy.structure.atoms import PhonopyAtoms
import numpy as np
# Create three Phonopy instances at different volumes
# Reference volume (equilibrium)
ph0 = Phonopy(unitcell_eq, supercell_matrix)
ph0.force_constants = fc_eq
# Expanded volume (+2% strain)
unitcell_plus = unitcell_eq.copy()
unitcell_plus.cell = unitcell_eq.cell * 1.02**(1/3)
ph_plus = Phonopy(unitcell_plus, supercell_matrix)
ph_plus.force_constants = fc_plus
# Compressed volume (-2% strain)
unitcell_minus = unitcell_eq.copy()
unitcell_minus.cell = unitcell_eq.cell * 0.98**(1/3)
ph_minus = Phonopy(unitcell_minus, supercell_matrix)
ph_minus.force_constants = fc_minus
# Calculate Grüneisen parameters
gruneisen = PhonopyGruneisen(ph0, ph_plus, ph_minus)Calculate Grüneisen parameters on regular q-point meshes for Brillouin zone integration.
def set_mesh(
self,
mesh: ArrayLike,
is_gamma_center: bool = False,
shift: ArrayLike = None,
is_time_reversal: bool = True,
is_mesh_symmetry: bool = True
):
"""
Set q-point mesh for Grüneisen parameter calculation.
Parameters:
- mesh: Mesh dimensions [nx, ny, nz]
- is_gamma_center: Center mesh on Gamma point
- shift: Mesh shift from origin [sx, sy, sz]
- is_time_reversal: Apply time-reversal symmetry
- is_mesh_symmetry: Use crystal symmetry to reduce q-points
"""
def get_mesh(self):
"""
Get GruneisenMesh object with calculated parameters.
Returns:
GruneisenMesh object containing frequencies, Grüneisen parameters,
and integration weights for each q-point
"""
def write_yaml_mesh(self, filename: str = "gruneisen_mesh.yaml"):
"""
Write mesh Grüneisen parameters to YAML file.
Parameters:
- filename: Output YAML file name
"""
def write_hdf5_mesh(self, filename: str = "gruneisen_mesh.hdf5"):
"""
Write mesh Grüneisen parameters to HDF5 file.
Parameters:
- filename: Output HDF5 file name
"""Example:
# Calculate Grüneisen parameters on 20x20x20 mesh
gruneisen.set_mesh([20, 20, 20], is_gamma_center=True)
mesh = gruneisen.get_mesh()
# Access mesh data
frequencies = mesh.frequencies # Shape: (nq, nbands)
gruneisen_params = mesh.gruneisen # Shape: (nq, nbands)
weights = mesh.weights # Integration weights
# Save results
gruneisen.write_yaml_mesh("gruneisen_mesh.yaml")
gruneisen.write_hdf5_mesh("gruneisen_mesh.hdf5")Calculate Grüneisen parameters along high-symmetry paths in the Brillouin zone.
def set_band_structure(self, bands: ArrayLike):
"""
Set band structure path for Grüneisen parameter calculation.
Parameters:
- bands: k-point paths as list of segments or 2D array
"""
def get_band_structure(self):
"""
Get GruneisenBandStructure object with calculated parameters.
Returns:
GruneisenBandStructure object containing frequencies, Grüneisen
parameters, and q-point paths
"""
def write_yaml_band_structure(
self,
filename: str = "gruneisen_band.yaml"
):
"""
Write band structure Grüneisen parameters to YAML file.
Parameters:
- filename: Output YAML file name
"""Example:
# Define high-symmetry path
bands = [[[0.0, 0.0, 0.0], # Gamma
[0.5, 0.0, 0.5]], # X
[[0.5, 0.0, 0.5], # X
[0.5, 0.25, 0.75]], # W
[[0.5, 0.25, 0.75], # W
[0.0, 0.0, 0.0]]] # Gamma
# Calculate Grüneisen parameters along path
gruneisen.set_band_structure(bands)
band_structure = gruneisen.get_band_structure()
# Access band structure data
distances = band_structure.distances
frequencies = band_structure.frequencies # Shape: (nq, nbands)
gruneisen_params = band_structure.gruneisen # Shape: (nq, nbands)
# Save results
gruneisen.write_yaml_band_structure("gruneisen_band.yaml")Generate plots of Grüneisen parameters versus frequency and along band structure paths.
def plot_mesh(
self,
cutoff_frequency: float = None,
color_scheme: str = None,
marker: str = 'o',
markersize: float = None,
xlabel: str = None,
ylabel: str = None,
drawing_area: ArrayLike = None,
aspect: str = 'auto'
):
"""
Plot Grüneisen parameters vs frequency from mesh calculation.
Parameters:
- cutoff_frequency: Exclude frequencies below cutoff (THz)
- color_scheme: Color scheme for plotting
- marker: Matplotlib marker style
- markersize: Size of plot markers
- xlabel: X-axis label
- ylabel: Y-axis label
- drawing_area: Plot area bounds [xmin, xmax, ymin, ymax]
- aspect: Axes aspect ratio
"""
def plot_band_structure(
self,
epsilon: float = 1e-4,
color_scheme: str = None,
marker: str = 'o',
markersize: float = None,
xlabel: str = None,
ylabel: str = None,
drawing_area: ArrayLike = None,
aspect: str = 'auto'
):
"""
Plot Grüneisen parameters along band structure path.
Parameters:
- epsilon: Threshold for avoiding small frequencies
- color_scheme: Color scheme for plotting
- marker: Matplotlib marker style
- markersize: Size of plot markers
- xlabel: X-axis label
- ylabel: Y-axis label
- drawing_area: Plot area bounds [xmin, xmax, ymin, ymax]
- aspect: Axes aspect ratio
"""Example:
import matplotlib.pyplot as plt
# Plot Grüneisen parameters vs frequency from mesh
gruneisen.plot_mesh(cutoff_frequency=1.0)
plt.title("Grüneisen Parameters vs Frequency")
plt.show()
# Plot Grüneisen parameters along band structure
gruneisen.plot_band_structure(epsilon=1e-3)
plt.title("Grüneisen Parameters Band Structure")
plt.show()Mesh-based Grüneisen parameter calculation results.
class GruneisenMesh:
@property
def frequencies(self) -> ndarray:
"""Phonon frequencies at mesh q-points (THz)."""
@property
def gruneisen(self) -> ndarray:
"""Grüneisen parameters at mesh q-points."""
@property
def weights(self) -> ndarray:
"""Integration weights for mesh q-points."""
@property
def qpoints(self) -> ndarray:
"""q-point coordinates in reciprocal space."""
def get_average_gruneisen_parameter(
self,
t: float = None,
cutoff_frequency: float = None
) -> float:
"""
Get average Grüneisen parameter.
Parameters:
- t: Temperature for thermal weighting (K)
- cutoff_frequency: Exclude frequencies below cutoff (THz)
Returns:
Average Grüneisen parameter
"""Band structure Grüneisen parameter calculation results.
class GruneisenBandStructure:
@property
def distances(self) -> ndarray:
"""Distances along band structure path."""
@property
def frequencies(self) -> ndarray:
"""Phonon frequencies along path (THz)."""
@property
def gruneisen(self) -> ndarray:
"""Grüneisen parameters along path."""
@property
def qpoints(self) -> ndarray:
"""q-point coordinates along path."""
@property
def labels(self) -> list:
"""High-symmetry point labels."""from phonopy import Phonopy, PhonopyGruneisen
from phonopy.structure.atoms import PhonopyAtoms
import numpy as np
import matplotlib.pyplot as plt
# 1. Create structures at three volumes
strain = 0.02 # 2% volume strain
# Equilibrium structure
unitcell_eq = PhonopyAtoms(...) # Your equilibrium structure
supercell_matrix = [[2, 0, 0], [0, 2, 0], [0, 0, 2]]
# Expanded structure (+strain)
unitcell_plus = unitcell_eq.copy()
scale_plus = (1 + strain)**(1/3)
unitcell_plus.cell = unitcell_eq.cell * scale_plus
# Compressed structure (-strain)
unitcell_minus = unitcell_eq.copy()
scale_minus = (1 - strain)**(1/3)
unitcell_minus.cell = unitcell_eq.cell * scale_minus
# 2. Calculate force constants for each volume
# (This requires separate DFT calculations at each volume)
ph0 = Phonopy(unitcell_eq, supercell_matrix)
ph_plus = Phonopy(unitcell_plus, supercell_matrix)
ph_minus = Phonopy(unitcell_minus, supercell_matrix)
# Set force constants from DFT calculations
ph0.force_constants = fc_equilibrium
ph_plus.force_constants = fc_expanded
ph_minus.force_constants = fc_compressed
# 3. Calculate Grüneisen parameters
gruneisen = PhonopyGruneisen(ph0, ph_plus, ph_minus)
# 4. Mesh calculation
gruneisen.set_mesh([20, 20, 20])
mesh = gruneisen.get_mesh()
# Get average Grüneisen parameter
avg_gruneisen = mesh.get_average_gruneisen_parameter(t=300)
print(f"Average Grüneisen parameter at 300K: {avg_gruneisen:.3f}")
# 5. Band structure calculation
bands = [[[0.0, 0.0, 0.0], [0.5, 0.0, 0.5]], # Gamma-X
[[0.5, 0.0, 0.5], [0.5, 0.5, 0.0]], # X-M
[[0.5, 0.5, 0.0], [0.0, 0.0, 0.0]]] # M-Gamma
gruneisen.set_band_structure(bands)
band_structure = gruneisen.get_band_structure()
# 6. Plotting and analysis
# Plot frequency vs Grüneisen parameter
gruneisen.plot_mesh(cutoff_frequency=1.0)
plt.xlabel("Frequency (THz)")
plt.ylabel("Grüneisen Parameter")
plt.title("Grüneisen Parameters vs Frequency")
plt.grid(True)
plt.show()
# Plot band structure
gruneisen.plot_band_structure()
plt.xlabel("Wave Vector")
plt.ylabel("Grüneisen Parameter")
plt.title("Grüneisen Parameters Band Structure")
plt.show()
# 7. Save results
gruneisen.write_yaml_mesh("gruneisen_mesh.yaml")
gruneisen.write_yaml_band_structure("gruneisen_band.yaml")Grüneisen parameters γ are defined as:
γ = -(V/ω)(∂ω/∂V)
Where:
Interpretation:
Grüneisen parameters are essential for understanding:
Install with Tessl CLI
npx tessl i tessl/pypi-phonopy