tessl/pypi-phonopy

Comprehensive Python library for phonon calculations that enables lattice dynamics simulations and vibrational property analysis of crystalline materials.

—

Pending

Overview

Eval results

Files

Quasi-Harmonic Approximation (QHA)

Name: tessl/pypi-phonopy
Author: tessl

Perform quasi-harmonic approximation calculations to study temperature and pressure effects on material properties through phonon calculations at multiple volumes. QHA enables prediction of thermal expansion, bulk modulus, heat capacity, and other thermodynamic properties by combining phonon calculations with equation of state fitting.

Capabilities

QHA Initialization

Initialize QHA calculations with volume-dependent phonon data and optional electronic contributions.

class PhonopyQHA:
    def __init__(
        self,
        volumes: ArrayLike = None,
        electronic_energies: ArrayLike = None,
        temperatures: ArrayLike = None,
        free_energy: ArrayLike = None,
        cv: ArrayLike = None,
        entropy: ArrayLike = None,
        t_max: float = None,
        energy_plot_factor: float = None,
        pressure: float = None,
        eos: str = "vinet",
        verbose: bool = False
    ):
        """
        Initialize quasi-harmonic approximation calculation.
        
        Parameters:
        - volumes: Crystal volumes for each calculation (Å³)
        - electronic_energies: Electronic energies at each volume (eV)
        - temperatures: Temperature range for calculations (K)
        - free_energy: Phonon free energies [nvols, ntemps] (kJ/mol)
        - cv: Phonon heat capacities [nvols, ntemps] (J/K/mol)
        - entropy: Phonon entropies [nvols, ntemps] (J/K/mol)
        - t_max: Maximum temperature for analysis (K)
        - energy_plot_factor: Scaling factor for energy plots
        - pressure: Applied pressure (GPa)
        - eos: Equation of state ('vinet', 'murnaghan', 'birch_murnaghan', 'p3')
        - verbose: Print detailed calculation information
        """

Example:

from phonopy import Phonopy, PhonopyQHA
import numpy as np

# Volume range (typically ±10% around equilibrium)
volumes = np.array([95, 100, 105, 110, 115])  # Å³

# Electronic energies from DFT at each volume
electronic_energies = np.array([-10.5, -10.8, -10.6, -10.2, -9.7])  # eV

# Temperature range
temperatures = np.arange(0, 1000, 10)  # K

# Collect phonon data from Phonopy calculations at each volume
free_energies = []
heat_capacities = []
entropies = []

for i, vol in enumerate(volumes):
    # Create Phonopy instance for this volume
    ph = Phonopy(unitcell_scaled[i], supercell_matrix)
    ph.force_constants = force_constants[i]
    
    # Calculate thermal properties
    ph.run_mesh([20, 20, 20])
    ph.run_thermal_properties(temperatures)
    thermal = ph.get_thermal_properties_dict()
    
    free_energies.append(thermal['free_energy'])
    heat_capacities.append(thermal['heat_capacity'])
    entropies.append(thermal['entropy'])

# Create QHA instance
qha = PhonopyQHA(
    volumes=volumes,
    electronic_energies=electronic_energies,
    temperatures=temperatures,
    free_energy=np.array(free_energies).T,  # [ntemps, nvols]
    cv=np.array(heat_capacities).T,
    entropy=np.array(entropies).T,
    eos='vinet'
)

Bulk Modulus Analysis

Calculate bulk modulus as a function of volume and temperature.

def get_bulk_modulus(self) -> tuple:
    """
    Get bulk modulus data from equation of state fitting.
    
    Returns:
    Tuple containing:
    - bulk_modulus: Bulk modulus values (GPa)
    - volumes_for_fit: Volumes used in EOS fitting (Å³)  
    - parameters: EOS fitting parameters
    """

def get_bulk_modulus_temperature(self) -> tuple:
    """
    Get bulk modulus as function of temperature.
    
    Returns:
    Tuple of (temperatures, bulk_modulus_values)
    """

def write_bulk_modulus_temperature(
    self, 
    filename: str = "bulk_modulus-temperature.dat"
):
    """
    Write bulk modulus vs temperature to file.
    
    Parameters:
    - filename: Output file name
    """

def plot_bulk_modulus(self, thin_number: int = 10):
    """
    Plot bulk modulus vs volume with EOS fitting.
    
    Parameters:
    - thin_number: Thinning factor for temperature curves
    """

Example:

# Get bulk modulus analysis
bulk_modulus, volumes_fit, eos_params = qha.get_bulk_modulus()
print(f"Bulk modulus at 0K: {bulk_modulus[0]:.1f} GPa")

# Temperature-dependent bulk modulus
temps, bulk_mod_T = qha.get_bulk_modulus_temperature()

# Plot bulk modulus
qha.plot_bulk_modulus(thin_number=5)
plt.title("Bulk Modulus vs Volume")
plt.show()

# Save data
qha.write_bulk_modulus_temperature("bulk_modulus_T.dat")

Helmholtz Free Energy

Calculate and analyze Helmholtz free energy as function of volume and temperature.

def get_helmholtz_volume(self) -> tuple:
    """
    Get Helmholtz free energy vs volume data.
    
    Returns:
    Tuple containing:
    - temperatures: Temperature array (K)
    - volumes: Volume array (Å³)
    - helmholtz_energies: Free energies [ntemps, nvols] (eV)
    """

def write_helmholtz_volume(
    self, 
    filename: str = "helmholtz-volume.dat"
):
    """
    Write Helmholtz free energy vs volume to file.
    
    Parameters:
    - filename: Output file name
    """

Volume-Temperature Relations

Calculate equilibrium volume as function of temperature and thermal expansion.

def get_volume_temperature(self) -> tuple:
    """
    Get equilibrium volume vs temperature.
    
    Returns:
    Tuple of (temperatures, equilibrium_volumes)
    """

def get_thermal_expansion(self) -> tuple:
    """
    Get thermal expansion coefficient vs temperature.
    
    Returns:
    Tuple of (temperatures, thermal_expansion_coefficients)
    """

def write_volume_temperature(
    self, 
    filename: str = "volume-temperature.dat"
):
    """
    Write volume vs temperature to file.
    
    Parameters:
    - filename: Output file name
    """

def write_thermal_expansion(
    self, 
    filename: str = "thermal_expansion.dat"
):
    """
    Write thermal expansion vs temperature to file.
    
    Parameters:
    - filename: Output file name
    """

def plot_volume_temperature(self, exp_data: ArrayLike = None):
    """
    Plot volume vs temperature with optional experimental comparison.
    
    Parameters:
    - exp_data: Experimental data as [temperatures, volumes] for comparison
    """

Example:

# Volume-temperature relationship
temps, volumes = qha.get_volume_temperature()
print(f"Volume at 300K: {volumes[30]:.2f} Å³")  # T[30] = 300K

# Thermal expansion coefficient
temps, alpha = qha.get_thermal_expansion()
print(f"Thermal expansion at 300K: {alpha[30]:.2e} K⁻¹")

# Plot with experimental data (if available)
exp_temps = [300, 400, 500, 600]
exp_volumes = [100.2, 100.8, 101.5, 102.3]
qha.plot_volume_temperature(exp_data=[exp_temps, exp_volumes])
plt.show()

# Save results
qha.write_volume_temperature("volume_T.dat")
qha.write_thermal_expansion("alpha_T.dat")

Gibbs Free Energy and Heat Capacity

Calculate temperature-dependent thermodynamic properties at equilibrium volume.

def get_gibbs_temperature(self) -> tuple:
    """
    Get Gibbs free energy vs temperature at equilibrium volume.
    
    Returns:
    Tuple of (temperatures, gibbs_free_energies)
    """

def get_heat_capacity_P_numerical(self) -> tuple:
    """
    Get heat capacity at constant pressure (numerical derivative).
    
    Returns:
    Tuple of (temperatures, heat_capacities_P)
    """

def get_heat_capacity_P_polyfit(self) -> tuple:
    """
    Get heat capacity at constant pressure (polynomial fit derivative).
    
    Returns:
    Tuple of (temperatures, heat_capacities_P)
    """

def write_gibbs_temperature(
    self, 
    filename: str = "gibbs-temperature.dat"
):
    """
    Write Gibbs free energy vs temperature to file.
    
    Parameters:
    - filename: Output file name
    """

def write_heat_capacity_P_numerical(
    self, 
    filename: str = "Cp-temperature.dat"
):
    """
    Write heat capacity at constant pressure to file.
    
    Parameters:
    - filename: Output file name
    """

Grüneisen Temperature

Calculate average Grüneisen parameter as function of temperature.

def get_gruneisen_temperature(self) -> tuple:
    """
    Get average Grüneisen parameter vs temperature.
    
    Returns:
    Tuple of (temperatures, gruneisen_parameters)
    """

def write_gruneisen_temperature(
    self, 
    filename: str = "gruneisen-temperature.dat"
):
    """
    Write Grüneisen parameter vs temperature to file.
    
    Parameters:
    - filename: Output file name
    """

Comprehensive Plotting

Generate comprehensive QHA analysis plots.

def plot_qha(
    self,
    thin_number: int = 10,
    volume_temp_experimental: ArrayLike = None,
    thermal_expansion_experimental: ArrayLike = None,
    cp_experimental: ArrayLike = None,
    cv_experimental: ArrayLike = None
):
    """
    Generate comprehensive QHA plots including volume-temperature, 
    thermal expansion, heat capacity, and Grüneisen parameter.
    
    Parameters:
    - thin_number: Thinning factor for curves
    - volume_temp_experimental: Experimental volume-temperature data
    - thermal_expansion_experimental: Experimental thermal expansion data
    - cp_experimental: Experimental Cp data
    - cv_experimental: Experimental Cv data
    """

def plot_pdf_qha(
    self,
    filename: str = "qha.pdf",
    thin_number: int = 10
):
    """
    Generate PDF file with all QHA plots.
    
    Parameters:
    - filename: Output PDF file name
    - thin_number: Thinning factor for curves
    """

Properties Access

Access calculated QHA properties directly.

@property
def bulk_modulus(self) -> ndarray:
    """Bulk modulus values (GPa)."""

@property  
def thermal_expansion(self) -> ndarray:
    """Thermal expansion coefficients (K⁻¹)."""

@property
def helmholtz_volume(self) -> ndarray:
    """Helmholtz free energies at each volume and temperature (eV)."""

@property
def volume_temperature(self) -> ndarray:
    """Equilibrium volumes at each temperature (Å³)."""

@property
def gibbs_temperature(self) -> ndarray:
    """Gibbs free energies at equilibrium volume vs temperature (eV)."""

@property
def bulk_modulus_temperature(self) -> ndarray:
    """Bulk modulus vs temperature (GPa)."""

@property
def heat_capacity_P_numerical(self) -> ndarray:
    """Heat capacity at constant pressure, numerical (J/K/mol)."""

@property
def heat_capacity_P_polyfit(self) -> ndarray:
    """Heat capacity at constant pressure, polynomial fit (J/K/mol)."""

@property
def gruneisen_temperature(self) -> ndarray:
    """Average Grüneisen parameter vs temperature."""

Complete QHA Workflow

from phonopy import Phonopy, PhonopyQHA
from phonopy.structure.atoms import PhonopyAtoms
import numpy as np
import matplotlib.pyplot as plt

# 1. Define volume range (typically ±10% around equilibrium)
V0 = 100.0  # Equilibrium volume (Å³)
volume_strains = np.array([-0.08, -0.04, 0.0, 0.04, 0.08])
volumes = V0 * (1 + volume_strains)

# Electronic energies from DFT calculations at each volume
electronic_energies = np.array([-10.85, -10.92, -10.89, -10.81, -10.68])

# 2. Temperature range for thermal properties
temperatures = np.arange(0, 1000, 10)  # 0 to 1000 K in 10 K steps

# 3. Calculate phonon thermal properties at each volume
phonon_free_energies = []
phonon_heat_capacities = []  
phonon_entropies = []

for i, vol in enumerate(volumes):
    # Create scaled unit cell for this volume
    scale_factor = (vol / V0)**(1/3)
    unitcell_scaled = unitcell.copy()
    unitcell_scaled.cell *= scale_factor
    
    # Phonopy calculation
    ph = Phonopy(unitcell_scaled, [[2,0,0],[0,2,0],[0,0,2]])
    ph.force_constants = force_constants_at_volume[i]  # From DFT
    
    # Calculate thermal properties
    ph.run_mesh([20, 20, 20])
    ph.run_thermal_properties(temperatures)
    thermal = ph.get_thermal_properties_dict()
    
    phonon_free_energies.append(thermal['free_energy'])
    phonon_heat_capacities.append(thermal['heat_capacity'])
    phonon_entropies.append(thermal['entropy'])

# Convert to proper array shape [ntemps, nvols]
free_energy = np.array(phonon_free_energies).T
cv = np.array(phonon_heat_capacities).T
entropy = np.array(phonon_entropies).T

# 4. Initialize QHA calculation
qha = PhonopyQHA(
    volumes=volumes,
    electronic_energies=electronic_energies,
    temperatures=temperatures,
    free_energy=free_energy,
    cv=cv,
    entropy=entropy,
    eos='vinet',  # Vinet equation of state
    verbose=True
)

# 5. Extract thermodynamic properties

# Bulk modulus analysis
bulk_modulus_0K = qha.bulk_modulus[0]
print(f"Bulk modulus at 0 K: {bulk_modulus_0K:.1f} GPa")

# Volume-temperature relationship  
temps, eq_volumes = qha.get_volume_temperature()
V_300K = eq_volumes[30]  # Volume at 300 K
print(f"Volume at 300 K: {V_300K:.2f} Å³")

# Thermal expansion
temps, alpha = qha.get_thermal_expansion()  
alpha_300K = alpha[30]  # Thermal expansion at 300 K
print(f"Thermal expansion at 300 K: {alpha_300K:.2e} K⁻¹")

# Heat capacity at constant pressure
temps, Cp = qha.get_heat_capacity_P_numerical()
Cp_300K = Cp[30]  # Heat capacity at 300 K
print(f"Heat capacity (Cp) at 300 K: {Cp_300K:.1f} J/K/mol")

# Average Grüneisen parameter
temps, gamma_avg = qha.get_gruneisen_temperature()
gamma_300K = gamma_avg[30]
print(f"Average Grüneisen parameter at 300 K: {gamma_300K:.2f}")

# 6. Generate comprehensive plots
qha.plot_qha(thin_number=5)
plt.suptitle("Quasi-Harmonic Approximation Results")
plt.tight_layout()
plt.show()

# Plot specific properties
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Volume vs temperature
axes[0,0].plot(temps, eq_volumes)
axes[0,0].set_xlabel("Temperature (K)")
axes[0,0].set_ylabel("Volume (Å³)")
axes[0,0].set_title("Equilibrium Volume vs Temperature")
axes[0,0].grid(True)

# Thermal expansion
axes[0,1].plot(temps[1:], alpha[1:])  # Skip T=0
axes[0,1].set_xlabel("Temperature (K)")  
axes[0,1].set_ylabel("Thermal Expansion (K⁻¹)")
axes[0,1].set_title("Thermal Expansion Coefficient")
axes[0,1].grid(True)

# Heat capacity
axes[1,0].plot(temps, Cp, label='Cp (constant P)')
axes[1,0].plot(temps, cv.mean(axis=1), label='Cv (constant V)')
axes[1,0].set_xlabel("Temperature (K)")
axes[1,0].set_ylabel("Heat Capacity (J/K/mol)")
axes[1,0].set_title("Heat Capacity")
axes[1,0].legend()
axes[1,0].grid(True)

# Bulk modulus vs temperature
temps_BM, BM = qha.get_bulk_modulus_temperature()
axes[1,1].plot(temps_BM, BM)
axes[1,1].set_xlabel("Temperature (K)")
axes[1,1].set_ylabel("Bulk Modulus (GPa)")
axes[1,1].set_title("Bulk Modulus vs Temperature")
axes[1,1].grid(True)

plt.tight_layout()
plt.show()

# 7. Save results to files
qha.write_volume_temperature("volume_temperature.dat")
qha.write_thermal_expansion("thermal_expansion.dat")  
qha.write_heat_capacity_P_numerical("Cp_temperature.dat")
qha.write_bulk_modulus_temperature("bulk_modulus_temperature.dat")
qha.write_gruneisen_temperature("gruneisen_temperature.dat")

# Save comprehensive plot to PDF
qha.plot_pdf_qha("qha_analysis.pdf")
print("QHA analysis complete. Results saved to files.")

Physical Interpretation

The quasi-harmonic approximation assumes:

Harmonic phonons at each volume
Volume-dependent frequencies captured through multiple calculations
No phonon-phonon interactions (pure harmonic limit)

Key outputs and their significance:

Thermal expansion α: How much material expands per unit temperature increase
Bulk modulus B: Material's resistance to uniform compression
Heat capacity Cp: Heat required to raise temperature at constant pressure
Grüneisen parameter γ: Measure of anharmonicity and mode coupling

Limitations: