0
# Quasi-Harmonic Approximation (QHA)
1

2
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.
3

4
## Capabilities
5

6
### QHA Initialization
7

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

10
```python { .api }
11
class PhonopyQHA:
12
    def __init__(
13
        self,
14
        volumes: ArrayLike = None,
15
        electronic_energies: ArrayLike = None,
16
        temperatures: ArrayLike = None,
17
        free_energy: ArrayLike = None,
18
        cv: ArrayLike = None,
19
        entropy: ArrayLike = None,
20
        t_max: float = None,
21
        energy_plot_factor: float = None,
22
        pressure: float = None,
23
        eos: str = "vinet",
24
        verbose: bool = False
25
    ):
26
        """
27
        Initialize quasi-harmonic approximation calculation.
28
        
29
        Parameters:
30
        - volumes: Crystal volumes for each calculation (ų)
31
        - electronic_energies: Electronic energies at each volume (eV)
32
        - temperatures: Temperature range for calculations (K)
33
        - free_energy: Phonon free energies [nvols, ntemps] (kJ/mol)
34
        - cv: Phonon heat capacities [nvols, ntemps] (J/K/mol)
35
        - entropy: Phonon entropies [nvols, ntemps] (J/K/mol)
36
        - t_max: Maximum temperature for analysis (K)
37
        - energy_plot_factor: Scaling factor for energy plots
38
        - pressure: Applied pressure (GPa)
39
        - eos: Equation of state ('vinet', 'murnaghan', 'birch_murnaghan', 'p3')
40
        - verbose: Print detailed calculation information
41
        """
42
```
43

44
**Example:**
45

46
```python
47
from phonopy import Phonopy, PhonopyQHA
48
import numpy as np
49

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

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

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

59
# Collect phonon data from Phonopy calculations at each volume
60
free_energies = []
61
heat_capacities = []
62
entropies = []
63

64
for i, vol in enumerate(volumes):
65
    # Create Phonopy instance for this volume
66
    ph = Phonopy(unitcell_scaled[i], supercell_matrix)
67
    ph.force_constants = force_constants[i]
68
    
69
    # Calculate thermal properties
70
    ph.run_mesh([20, 20, 20])
71
    ph.run_thermal_properties(temperatures)
72
    thermal = ph.get_thermal_properties_dict()
73
    
74
    free_energies.append(thermal['free_energy'])
75
    heat_capacities.append(thermal['heat_capacity'])
76
    entropies.append(thermal['entropy'])
77

78
# Create QHA instance
79
qha = PhonopyQHA(
80
    volumes=volumes,
81
    electronic_energies=electronic_energies,
82
    temperatures=temperatures,
83
    free_energy=np.array(free_energies).T,  # [ntemps, nvols]
84
    cv=np.array(heat_capacities).T,
85
    entropy=np.array(entropies).T,
86
    eos='vinet'
87
)
88
```
89

90
### Bulk Modulus Analysis
91

92
Calculate bulk modulus as a function of volume and temperature.
93

94
```python { .api }
95
def get_bulk_modulus(self) -> tuple:
96
    """
97
    Get bulk modulus data from equation of state fitting.
98
    
99
    Returns:
100
    Tuple containing:
101
    - bulk_modulus: Bulk modulus values (GPa)
102
    - volumes_for_fit: Volumes used in EOS fitting (ų)  
103
    - parameters: EOS fitting parameters
104
    """
105

106
def get_bulk_modulus_temperature(self) -> tuple:
107
    """
108
    Get bulk modulus as function of temperature.
109
    
110
    Returns:
111
    Tuple of (temperatures, bulk_modulus_values)
112
    """
113

114
def write_bulk_modulus_temperature(
115
    self, 
116
    filename: str = "bulk_modulus-temperature.dat"
117
):
118
    """
119
    Write bulk modulus vs temperature to file.
120
    
121
    Parameters:
122
    - filename: Output file name
123
    """
124

125
def plot_bulk_modulus(self, thin_number: int = 10):
126
    """
127
    Plot bulk modulus vs volume with EOS fitting.
128
    
129
    Parameters:
130
    - thin_number: Thinning factor for temperature curves
131
    """
132
```
133

134
**Example:**
135

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

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

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

149
# Save data
150
qha.write_bulk_modulus_temperature("bulk_modulus_T.dat")
151
```
152

153
### Helmholtz Free Energy
154

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

157
```python { .api }
158
def get_helmholtz_volume(self) -> tuple:
159
    """
160
    Get Helmholtz free energy vs volume data.
161
    
162
    Returns:
163
    Tuple containing:
164
    - temperatures: Temperature array (K)
165
    - volumes: Volume array (ų)
166
    - helmholtz_energies: Free energies [ntemps, nvols] (eV)
167
    """
168

169
def write_helmholtz_volume(
170
    self, 
171
    filename: str = "helmholtz-volume.dat"
172
):
173
    """
174
    Write Helmholtz free energy vs volume to file.
175
    
176
    Parameters:
177
    - filename: Output file name
178
    """
179
```
180

181
### Volume-Temperature Relations
182

183
Calculate equilibrium volume as function of temperature and thermal expansion.
184

185
```python { .api }
186
def get_volume_temperature(self) -> tuple:
187
    """
188
    Get equilibrium volume vs temperature.
189
    
190
    Returns:
191
    Tuple of (temperatures, equilibrium_volumes)
192
    """
193

194
def get_thermal_expansion(self) -> tuple:
195
    """
196
    Get thermal expansion coefficient vs temperature.
197
    
198
    Returns:
199
    Tuple of (temperatures, thermal_expansion_coefficients)
200
    """
201

202
def write_volume_temperature(
203
    self, 
204
    filename: str = "volume-temperature.dat"
205
):
206
    """
207
    Write volume vs temperature to file.
208
    
209
    Parameters:
210
    - filename: Output file name
211
    """
212

213
def write_thermal_expansion(
214
    self, 
215
    filename: str = "thermal_expansion.dat"
216
):
217
    """
218
    Write thermal expansion vs temperature to file.
219
    
220
    Parameters:
221
    - filename: Output file name
222
    """
223

224
def plot_volume_temperature(self, exp_data: ArrayLike = None):
225
    """
226
    Plot volume vs temperature with optional experimental comparison.
227
    
228
    Parameters:
229
    - exp_data: Experimental data as [temperatures, volumes] for comparison
230
    """
231
```
232

233
**Example:**
234

235
```python
236
# Volume-temperature relationship
237
temps, volumes = qha.get_volume_temperature()
238
print(f"Volume at 300K: {volumes[30]:.2f} ų")  # T[30] = 300K
239

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

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

250
# Save results
251
qha.write_volume_temperature("volume_T.dat")
252
qha.write_thermal_expansion("alpha_T.dat")
253
```
254

255
### Gibbs Free Energy and Heat Capacity
256

257
Calculate temperature-dependent thermodynamic properties at equilibrium volume.
258

259
```python { .api }
260
def get_gibbs_temperature(self) -> tuple:
261
    """
262
    Get Gibbs free energy vs temperature at equilibrium volume.
263
    
264
    Returns:
265
    Tuple of (temperatures, gibbs_free_energies)
266
    """
267

268
def get_heat_capacity_P_numerical(self) -> tuple:
269
    """
270
    Get heat capacity at constant pressure (numerical derivative).
271
    
272
    Returns:
273
    Tuple of (temperatures, heat_capacities_P)
274
    """
275

276
def get_heat_capacity_P_polyfit(self) -> tuple:
277
    """
278
    Get heat capacity at constant pressure (polynomial fit derivative).
279
    
280
    Returns:
281
    Tuple of (temperatures, heat_capacities_P)
282
    """
283

284
def write_gibbs_temperature(
285
    self, 
286
    filename: str = "gibbs-temperature.dat"
287
):
288
    """
289
    Write Gibbs free energy vs temperature to file.
290
    
291
    Parameters:
292
    - filename: Output file name
293
    """
294

295
def write_heat_capacity_P_numerical(
296
    self, 
297
    filename: str = "Cp-temperature.dat"
298
):
299
    """
300
    Write heat capacity at constant pressure to file.
301
    
302
    Parameters:
303
    - filename: Output file name
304
    """
305
```
306

307
### Grüneisen Temperature
308

309
Calculate average Grüneisen parameter as function of temperature.
310

311
```python { .api }
312
def get_gruneisen_temperature(self) -> tuple:
313
    """
314
    Get average Grüneisen parameter vs temperature.
315
    
316
    Returns:
317
    Tuple of (temperatures, gruneisen_parameters)
318
    """
319

320
def write_gruneisen_temperature(
321
    self, 
322
    filename: str = "gruneisen-temperature.dat"
323
):
324
    """
325
    Write Grüneisen parameter vs temperature to file.
326
    
327
    Parameters:
328
    - filename: Output file name
329
    """
330
```
331

332
### Comprehensive Plotting
333

334
Generate comprehensive QHA analysis plots.
335

336
```python { .api }
337
def plot_qha(
338
    self,
339
    thin_number: int = 10,
340
    volume_temp_experimental: ArrayLike = None,
341
    thermal_expansion_experimental: ArrayLike = None,
342
    cp_experimental: ArrayLike = None,
343
    cv_experimental: ArrayLike = None
344
):
345
    """
346
    Generate comprehensive QHA plots including volume-temperature, 
347
    thermal expansion, heat capacity, and Grüneisen parameter.
348
    
349
    Parameters:
350
    - thin_number: Thinning factor for curves
351
    - volume_temp_experimental: Experimental volume-temperature data
352
    - thermal_expansion_experimental: Experimental thermal expansion data
353
    - cp_experimental: Experimental Cp data
354
    - cv_experimental: Experimental Cv data
355
    """
356

357
def plot_pdf_qha(
358
    self,
359
    filename: str = "qha.pdf",
360
    thin_number: int = 10
361
):
362
    """
363
    Generate PDF file with all QHA plots.
364
    
365
    Parameters:
366
    - filename: Output PDF file name
367
    - thin_number: Thinning factor for curves
368
    """
369
```
370

371
### Properties Access
372

373
Access calculated QHA properties directly.
374

375
```python { .api }
376
@property
377
def bulk_modulus(self) -> ndarray:
378
    """Bulk modulus values (GPa)."""
379

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

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

388
@property
389
def volume_temperature(self) -> ndarray:
390
    """Equilibrium volumes at each temperature (ų)."""
391

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

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

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

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

408
@property
409
def gruneisen_temperature(self) -> ndarray:
410
    """Average Grüneisen parameter vs temperature."""
411
```
412

413
## Complete QHA Workflow
414

415
```python
416
from phonopy import Phonopy, PhonopyQHA
417
from phonopy.structure.atoms import PhonopyAtoms
418
import numpy as np
419
import matplotlib.pyplot as plt
420

421
# 1. Define volume range (typically ±10% around equilibrium)
422
V0 = 100.0  # Equilibrium volume (ų)
423
volume_strains = np.array([-0.08, -0.04, 0.0, 0.04, 0.08])
424
volumes = V0 * (1 + volume_strains)
425

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

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

432
# 3. Calculate phonon thermal properties at each volume
433
phonon_free_energies = []
434
phonon_heat_capacities = []  
435
phonon_entropies = []
436

437
for i, vol in enumerate(volumes):
438
    # Create scaled unit cell for this volume
439
    scale_factor = (vol / V0)**(1/3)
440
    unitcell_scaled = unitcell.copy()
441
    unitcell_scaled.cell *= scale_factor
442
    
443
    # Phonopy calculation
444
    ph = Phonopy(unitcell_scaled, [[2,0,0],[0,2,0],[0,0,2]])
445
    ph.force_constants = force_constants_at_volume[i]  # From DFT
446
    
447
    # Calculate thermal properties
448
    ph.run_mesh([20, 20, 20])
449
    ph.run_thermal_properties(temperatures)
450
    thermal = ph.get_thermal_properties_dict()
451
    
452
    phonon_free_energies.append(thermal['free_energy'])
453
    phonon_heat_capacities.append(thermal['heat_capacity'])
454
    phonon_entropies.append(thermal['entropy'])
455

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

461
# 4. Initialize QHA calculation
462
qha = PhonopyQHA(
463
    volumes=volumes,
464
    electronic_energies=electronic_energies,
465
    temperatures=temperatures,
466
    free_energy=free_energy,
467
    cv=cv,
468
    entropy=entropy,
469
    eos='vinet',  # Vinet equation of state
470
    verbose=True
471
)
472

473
# 5. Extract thermodynamic properties
474

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

479
# Volume-temperature relationship  
480
temps, eq_volumes = qha.get_volume_temperature()
481
V_300K = eq_volumes[30]  # Volume at 300 K
482
print(f"Volume at 300 K: {V_300K:.2f} ų")
483

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

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

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

499
# 6. Generate comprehensive plots
500
qha.plot_qha(thin_number=5)
501
plt.suptitle("Quasi-Harmonic Approximation Results")
502
plt.tight_layout()
503
plt.show()
504

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

508
# Volume vs temperature
509
axes[0,0].plot(temps, eq_volumes)
510
axes[0,0].set_xlabel("Temperature (K)")
511
axes[0,0].set_ylabel("Volume (ų)")
512
axes[0,0].set_title("Equilibrium Volume vs Temperature")
513
axes[0,0].grid(True)
514

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

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

531
# Bulk modulus vs temperature
532
temps_BM, BM = qha.get_bulk_modulus_temperature()
533
axes[1,1].plot(temps_BM, BM)
534
axes[1,1].set_xlabel("Temperature (K)")
535
axes[1,1].set_ylabel("Bulk Modulus (GPa)")
536
axes[1,1].set_title("Bulk Modulus vs Temperature")
537
axes[1,1].grid(True)
538

539
plt.tight_layout()
540
plt.show()
541

542
# 7. Save results to files
543
qha.write_volume_temperature("volume_temperature.dat")
544
qha.write_thermal_expansion("thermal_expansion.dat")  
545
qha.write_heat_capacity_P_numerical("Cp_temperature.dat")
546
qha.write_bulk_modulus_temperature("bulk_modulus_temperature.dat")
547
qha.write_gruneisen_temperature("gruneisen_temperature.dat")
548

549
# Save comprehensive plot to PDF
550
qha.plot_pdf_qha("qha_analysis.pdf")
551
print("QHA analysis complete. Results saved to files.")
552
```
553

554
## Physical Interpretation
555

556
The quasi-harmonic approximation assumes:
557
1. **Harmonic phonons** at each volume
558
2. **Volume-dependent frequencies** captured through multiple calculations  
559
3. **No phonon-phonon interactions** (pure harmonic limit)
560

561
**Key outputs and their significance:**
562

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

568
**Limitations:**
569
- Fails at high temperatures where anharmonic effects dominate
570
- Cannot capture temperature-dependent phonon linewidths
571
- May be inaccurate near phase transitions
572
- Electronic contributions must be included separately
573

574
**Applications:**
575
- Predicting thermal expansion of materials
576
- Understanding pressure effects on lattice dynamics
577
- Calculating thermodynamic properties for materials design
578
- Comparing with experimental thermal expansion measurements