# coding: utf-8
# Copyright (c) Pymatgen Development Team.
# Distributed under the terms of the MIT License.
"""
This module implements the Quasi-harmonic Debye approximation that can
be used to compute thermal properties.
See the following papers for more info:
http://doi.org/10.1016/j.comphy.2003.12.001 (2004)
http://doi.org/10.1103/PhysRevB.90.174107 (2014)
"""
from collections import defaultdict
import numpy as np
from scipy.constants import physical_constants
from scipy.integrate import quadrature
from scipy.misc import derivative
from scipy.optimize import minimize
from pymatgen.core.units import FloatWithUnit
from pymatgen.analysis.eos import EOS, PolynomialEOS
__author__ = "Kiran Mathew, Brandon Bocklund"
__credits__ = "Cormac Toher"
[docs]class QuasiharmonicDebyeApprox:
"""
Args:
energies (list): list of DFT energies in eV
volumes (list): list of volumes in Ang^3
structure (Structure):
t_min (float): min temperature
t_step (float): temperature step
t_max (float): max temperature
eos (str): equation of state used for fitting the energies and the
volumes.
options supported by pymatgen: "quadratic", "murnaghan", "birch",
"birch_murnaghan", "pourier_tarantola", "vinet",
"deltafactor", "numerical_eos"
pressure (float): in GPa, optional.
poisson (float): poisson ratio.
use_mie_gruneisen (bool): whether or not to use the mie-gruneisen
formulation to compute the gruneisen parameter.
The default is the slater-gamma formulation.
anharmonic_contribution (bool): whether or not to consider the anharmonic
contribution to the Debye temperature. Cannot be used with
use_mie_gruneisen. Defaults to False.
"""
def __init__(self, energies, volumes, structure, t_min=300.0, t_step=100,
t_max=300.0, eos="vinet", pressure=0.0, poisson=0.25,
use_mie_gruneisen=False, anharmonic_contribution=False):
self.energies = energies
self.volumes = volumes
self.structure = structure
self.temperature_min = t_min
self.temperature_max = t_max
self.temperature_step = t_step
self.eos_name = eos
self.pressure = pressure
self.poisson = poisson
self.use_mie_gruneisen = use_mie_gruneisen
self.anharmonic_contribution = anharmonic_contribution
if self.use_mie_gruneisen and self.anharmonic_contribution:
raise ValueError('The Mie-Gruneisen formulation and anharmonic contribution are circular referenced and '
'cannot be used together.')
self.mass = sum([e.atomic_mass for e in self.structure.species])
self.natoms = self.structure.composition.num_atoms
self.avg_mass = physical_constants["atomic mass constant"][0] * self.mass / self.natoms # kg
self.kb = physical_constants["Boltzmann constant in eV/K"][0]
self.hbar = physical_constants["Planck constant over 2 pi in eV s"][0]
self.gpa_to_ev_ang = 1./160.21766208 # 1 GPa in ev/Ang^3
self.gibbs_free_energy = [] # optimized values, eV
# list of temperatures for which the optimized values are available, K
self.temperatures = []
self.optimum_volumes = [] # in Ang^3
# fit E and V and get the bulk modulus(used to compute the Debye
# temperature)
print("Fitting E and V")
self.eos = EOS(eos)
self.ev_eos_fit = self.eos.fit(volumes, energies)
self.bulk_modulus = self.ev_eos_fit.b0_GPa # in GPa
self.optimize_gibbs_free_energy()
[docs] def optimize_gibbs_free_energy(self):
"""
Evaluate the gibbs free energy as a function of V, T and P i.e
G(V, T, P), minimize G(V, T, P) wrt V for each T and store the
optimum values.
Note: The data points for which the equation of state fitting fails
are skipped.
"""
temperatures = np.linspace(
self.temperature_min, self.temperature_max,
int(np.ceil((self.temperature_max - self.temperature_min) / self.temperature_step) + 1))
for t in temperatures:
try:
G_opt, V_opt = self.optimizer(t)
except Exception:
if len(temperatures) > 1:
print("EOS fitting failed, so skipping this data point, {}".format(t))
continue
else:
raise
self.gibbs_free_energy.append(G_opt)
self.temperatures.append(t)
self.optimum_volumes.append(V_opt)
[docs] def optimizer(self, temperature):
"""
Evaluate G(V, T, P) at the given temperature(and pressure) and
minimize it wrt V.
1. Compute the vibrational helmholtz free energy, A_vib.
2. Compute the gibbs free energy as a function of volume, temperature
and pressure, G(V,T,P).
3. Preform an equation of state fit to get the functional form of
gibbs free energy:G(V, T, P).
4. Finally G(V, P, T) is minimized with respect to V.
Args:
temperature (float): temperature in K
Returns:
float, float: G_opt(V_opt, T, P) in eV and V_opt in Ang^3.
"""
G_V = [] # G for each volume
# G = E(V) + PV + A_vib(V, T)
for i, v in enumerate(self.volumes):
G_V.append(self.energies[i] +
self.pressure * v * self.gpa_to_ev_ang +
self.vibrational_free_energy(temperature, v))
# fit equation of state, G(V, T, P)
eos_fit = self.eos.fit(self.volumes, G_V)
# minimize the fit eos wrt volume
# Note: the ref energy and the ref volume(E0 and V0) not necessarily
# the same as minimum energy and min volume.
volume_guess = eos_fit.volumes[np.argmin(eos_fit.energies)]
min_wrt_vol = minimize(eos_fit.func, volume_guess)
# G_opt=G(V_opt, T, P), V_opt
return min_wrt_vol.fun, min_wrt_vol.x[0]
[docs] def vibrational_free_energy(self, temperature, volume):
"""
Vibrational Helmholtz free energy, A_vib(V, T).
Eq(4) in doi.org/10.1016/j.comphy.2003.12.001
Args:
temperature (float): temperature in K
volume (float)
Returns:
float: vibrational free energy in eV
"""
y = self.debye_temperature(volume) / temperature
return self.kb * self.natoms * temperature * (
9./8. * y + 3 * np.log(1 - np.exp(-y)) - self.debye_integral(y))
[docs] def vibrational_internal_energy(self, temperature, volume):
"""
Vibrational internal energy, U_vib(V, T).
Eq(4) in doi.org/10.1016/j.comphy.2003.12.001
Args:
temperature (float): temperature in K
volume (float): in Ang^3
Returns:
float: vibrational internal energy in eV
"""
y = self.debye_temperature(volume) / temperature
return self.kb * self.natoms * temperature * (9./8. * y +
3*self.debye_integral(y))
[docs] def debye_temperature(self, volume):
"""
Calculates the debye temperature.
Eq(6) in doi.org/10.1016/j.comphy.2003.12.001. Thanks to Joey.
Eq(6) above is equivalent to Eq(3) in doi.org/10.1103/PhysRevB.37.790
which does not consider anharmonic effects. Eq(20) in the same paper
and Eq(18) in doi.org/10.1016/j.commatsci.2009.12.006 both consider
anharmonic contributions to the Debye temperature through the Gruneisen
parameter at 0K (Gruneisen constant).
The anharmonic contribution is toggled by setting the anharmonic_contribution
to True or False in the QuasiharmonicDebyeApprox constructor.
Args:
volume (float): in Ang^3
Returns:
float: debye temperature in K
"""
term1 = (2./3. * (1. + self.poisson) / (1. - 2. * self.poisson))**1.5
term2 = (1./3. * (1. + self.poisson) / (1. - self.poisson))**1.5
f = (3. / (2. * term1 + term2))**(1. / 3.)
debye = 2.9772e-11 * (volume / self.natoms) ** (-1. / 6.) * f * np.sqrt(self.bulk_modulus/self.avg_mass)
if self.anharmonic_contribution:
gamma = self.gruneisen_parameter(0, self.ev_eos_fit.v0) # 0K equilibrium Gruneisen parameter
return debye * (self.ev_eos_fit.v0 / volume) ** (gamma)
else:
return debye
[docs] @staticmethod
def debye_integral(y):
"""
Debye integral. Eq(5) in doi.org/10.1016/j.comphy.2003.12.001
Args:
y (float): debye temperature/T, upper limit
Returns:
float: unitless
"""
# floating point limit is reached around y=155, so values beyond that
# are set to the limiting value(T-->0, y --> \infty) of
# 6.4939394 (from wolfram alpha).
factor = 3. / y ** 3
if y < 155:
integral = quadrature(lambda x: x ** 3 / (np.exp(x) - 1.), 0, y)
return list(integral)[0] * factor
else:
return 6.493939 * factor
[docs] def gruneisen_parameter(self, temperature, volume):
"""
Slater-gamma formulation(the default):
gruneisen paramter = - d log(theta)/ d log(V)
= - ( 1/6 + 0.5 d log(B)/ d log(V) )
= - (1/6 + 0.5 V/B dB/dV),
where dB/dV = d^2E/dV^2 + V * d^3E/dV^3
Mie-gruneisen formulation:
Eq(31) in doi.org/10.1016/j.comphy.2003.12.001
Eq(7) in Blanco et. al. Joumal of Molecular Structure (Theochem)
368 (1996) 245-255
Also se J.P. Poirier, Introduction to the Physics of the Earth’s
Interior, 2nd ed. (Cambridge University Press, Cambridge,
2000) Eq(3.53)
Args:
temperature (float): temperature in K
volume (float): in Ang^3
Returns:
float: unitless
"""
if isinstance(self.eos, PolynomialEOS):
p = np.poly1d(self.eos.eos_params)
# first derivative of energy at 0K wrt volume evaluated at the
# given volume, in eV/Ang^3
dEdV = np.polyder(p, 1)(volume)
# second derivative of energy at 0K wrt volume evaluated at the
# given volume, in eV/Ang^6
d2EdV2 = np.polyder(p, 2)(volume)
# third derivative of energy at 0K wrt volume evaluated at the
# given volume, in eV/Ang^9
d3EdV3 = np.polyder(p, 3)(volume)
else:
func = self.ev_eos_fit.func
dEdV = derivative(func, volume, dx=1e-3)
d2EdV2 = derivative(func, volume, dx=1e-3, n=2, order=5)
d3EdV3 = derivative(func, volume, dx=1e-3, n=3, order=7)
# Mie-gruneisen formulation
if self.use_mie_gruneisen:
p0 = dEdV
return (self.gpa_to_ev_ang * volume *
(self.pressure + p0 / self.gpa_to_ev_ang) /
self.vibrational_internal_energy(temperature, volume))
# Slater-gamma formulation
# first derivative of bulk modulus wrt volume, eV/Ang^6
dBdV = d2EdV2 + d3EdV3 * volume
return -(1./6. + 0.5 * volume * dBdV /
FloatWithUnit(self.ev_eos_fit.b0_GPa, "GPa").to("eV ang^-3"))
[docs] def thermal_conductivity(self, temperature, volume):
"""
Eq(17) in 10.1103/PhysRevB.90.174107
Args:
temperature (float): temperature in K
volume (float): in Ang^3
Returns:
float: thermal conductivity in W/K/m
"""
gamma = self.gruneisen_parameter(temperature, volume)
theta_d = self.debye_temperature(volume) # K
theta_a = theta_d * self.natoms**(-1./3.) # K
prefactor = (0.849 * 3 * 4**(1./3.)) / (20. * np.pi**3)
# kg/K^3/s^3
prefactor = prefactor * (self.kb/self.hbar)**3 * self.avg_mass
kappa = prefactor / (gamma**2 - 0.514 * gamma + 0.228)
# kg/K/s^3 * Ang = (kg m/s^2)/(Ks)*1e-10
# = N/(Ks)*1e-10 = Nm/(Kms)*1e-10 = W/K/m*1e-10
kappa = kappa * theta_a**2 * volume**(1./3.) * 1e-10
return kappa
[docs] def get_summary_dict(self):
"""
Returns a dict with a summary of the computed properties.
"""
d = defaultdict(list)
d["pressure"] = self.pressure
d["poisson"] = self.poisson
d["mass"] = self.mass
d["natoms"] = int(self.natoms)
d["bulk_modulus"] = self.bulk_modulus
d["gibbs_free_energy"] = self.gibbs_free_energy
d["temperatures"] = self.temperatures
d["optimum_volumes"] = self.optimum_volumes
for v, t in zip(self.optimum_volumes, self.temperatures):
d["debye_temperature"].append(self.debye_temperature(v))
d["gruneisen_parameter"].append(self.gruneisen_parameter(t, v))
d["thermal_conductivity"].append(self.thermal_conductivity(t, v))
return d