Gibbs free energy#
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
from phonopy import PhonopyQHA
from phonopy.units import CastepToTHz
import profess
import ase_tools
Introduction#
Lithium adopts the BCC structure in ambient conditions, but it favors the FCC structure at low temperatures. See the following paper for details:
Ackland, et al., “Quantum and isotope effects in lithium metal,” Science 356, 1254 (2017).
In this example, we compute the Gibbs free energy for the two structures across a range of temperatures, eventually predicting the phase transition.
Static lattice approximation#
First, we define a helper function that creates a profess system with prescribed cell vectors and ion coordinates.
def build_system(box_vecs, xyz_coords):
system = (
profess.System.create(box_vecs, 1200, ['a','ev'])
.add_ions('potentials/li.gga.recpot', xyz_coords, 'a')
.add_electrons()
.add_wang_teter_functional() # kinetic energy
.add_hartree_functional() # electrostatic energies
.add_ion_electron_functional()
.add_ion_ion_interaction()
.add_perdew_burke_ernzerhof_functional() # exchange-correlation energy
)
return system.minimize_energy()
Then, we find the equilibrium geometries for BCC Li and FCC Li, observing that BCC is roughly 1 meV/atom higher in energy than FCC. These calculations assume T=0 and a static lattice.
# initial guess for volume [A^3]
vol_per_atom = 20
# determine bcc energy
box_vec = (2*vol_per_atom)**(1.0/3.0) * np.array([[-0.5,0.5,0.5],[0.5,-0.5,0.5],[0.5,0.5,-0.5]])
system = build_system(box_vec, [(0,0,0)])
print('----------------------\nRelaxing BCC structure\n----------------------')
ase_tools.minimize_forces_stress(system, 'BFGSLineSearch', 1e-4)
print('\n Final volume: {:4.2f} A3'.format(system.volume('a3')))
print(' Final energy: {:4.3f} eV\n'.format(system.energy('ev')))
ene_bcc = system.energy('ev')
# determine fcc energy
box_vec = (4*vol_per_atom)**(1.0/3.0) * np.array([[0.5,0.5,0.0],[0.0,0.5,0.5],[0.5,0.0,0.5]])
system = build_system(box_vec, [(0,0,0)])
print('----------------------\nRelaxing FCC structure\n----------------------')
ase_tools.minimize_forces_stress(system, 'BFGSLineSearch', 1e-4)
print('\n Final volume: {:4.2f} A3'.format(system.volume('a3')))
print(' Final energy: {:4.3f} eV\n'.format(system.energy('ev')))
ene_fcc = system.energy('ev')
----------------------
Relaxing BCC structure
----------------------
Step[ FC] Time Energy fmax
BFGSLineSearch: 0[ 0] 09:52:11 -0.270481 0.0007
BFGSLineSearch: 1[ 4] 09:52:13 -0.270484 0.0000
Final volume: 20.16 A3
Final energy: -7.360 eV
----------------------
Relaxing FCC structure
----------------------
Step[ FC] Time Energy fmax
BFGSLineSearch: 0[ 0] 09:52:13 -0.270516 0.0009
BFGSLineSearch: 1[ 3] 09:52:15 -0.270520 0.0004
BFGSLineSearch: 2[ 4] 09:52:15 -0.270521 0.0000
Final volume: 20.21 A3
Final energy: -7.361 eV
Gibbs free energy with Phonopy#
Finally, we use Phonopy to compute Gibbs free energies with the quasi-harmonic approximation.
See the Phonopy documentation for more details.
The thermal analysis explains why the bcc structure is observed at room temperature, despite fcc having slighly lower energy as $T\to 0$. The difference between the energies obtained in the previous section and the $T\to 0$ energies in the plot are the zero-point energies.
vols_per_atom = np.linspace(20,23,6)
temperatures = np.linspace(0,300,100)
for structure in ['bcc', 'fcc']:
# create empty arrays
volumes = np.zeros(vols_per_atom.size)
energies = np.zeros(vols_per_atom.size)
free_energies = np.zeros([temperatures.size, vols_per_atom.size])
entropies = np.zeros([temperatures.size, vols_per_atom.size])
heat_capacities = np.zeros([temperatures.size, vols_per_atom.size])
for i, vol_per_atom in enumerate(vols_per_atom):
# set box vectors
if structure=='bcc':
box_vec = (2*vol_per_atom)**(1.0/3.0) * np.array([[-0.5,0.5,0.5],[0.5,-0.5,0.5],[0.5,0.5,-0.5]])
elif structure=='fcc':
box_vec = (4*vol_per_atom)**(1.0/3.0) * np.array([[0.5,0.5,0.0],[0.0,0.5,0.5],[0.5,0.0,0.5]])
# compute volume and energy
system = build_system(box_vec, [(0,0,0)])
volumes[i] = system.volume('a3')
energies[i] = system.energy('ev')
# generate supercells with phonopy
unitcell = PhonopyAtoms(symbols=(['Li']*1), cell=box_vec, scaled_positions=[(0,0,0)])
supercell_matrix = [[2,0,0], [0,2,0], [0,0,2]]
phonon = Phonopy(unitcell, supercell_matrix, primitive_matrix=np.identity(3), factor=CastepToTHz)
phonon.generate_displacements(distance=0.01)
disps = phonon.displacements
supercells = phonon.supercells_with_displacements
# compute forces for all supercells
set_of_forces = []
for supercell in supercells:
system = build_system(supercell.cell, supercell.positions)
forces = np.array(system.forces())*27.2114/0.529177
drift_force = forces.sum(axis=0)
for j, _ in enumerate(forces):
forces[j] -= drift_force / forces.shape[0]
set_of_forces.append(forces)
# phonopy post-process
phonon.produce_force_constants(forces=set_of_forces)
# get thermal properties
phonon.run_mesh([48,48,48])
phonon.run_thermal_properties(temperatures=temperatures)
free_energies[:,i] = phonon.get_thermal_properties_dict()['free_energy']
entropies[:,i] = phonon.get_thermal_properties_dict()['entropy']
heat_capacities[:,i] = phonon.get_thermal_properties_dict()['heat_capacity']
# conduct quasi-harmonic analysis
qha = PhonopyQHA(volumes=volumes,
electronic_energies=energies,
temperatures=temperatures,
free_energy=free_energies,
cv=heat_capacities,
entropy=entropies,
eos='vinet',
verbose=False)
# extract gibbs energies
if structure=='bcc':
gibbs_bcc = qha.gibbs_temperature
elif structure=='fcc':
gibbs_fcc = qha.gibbs_temperature
plt.plot(temperatures[:-1], gibbs_bcc, label='bcc')
plt.plot(temperatures[:-1], gibbs_fcc, label='fcc')
plt.xlabel('Temperature [K]')
plt.ylabel('Gibbs Free Energy [eV/atom]')
plt.legend()
plt.show()
