"""Module for processing elastic tensor."""
# Reference: https://doi.org/10.1103/PhysRevB.98.014107
import numpy as np
from collections import OrderedDict
[docs]class ElasticTensor(object):
"""Module for processing elastic tensor."""
def __init__(self, et_tensor=[]):
"""Initialize class."""
self.et_tensor = et_tensor
@property
def voigt_modulus(self):
"""Get Voigt modulus."""
c = self.et_tensor
Kv = float(
(c[0][0] + c[1][1] + c[2][2]) + 2 * (c[0][1] + c[1][2] + c[2][0])
) / float(9)
Gv = float(
(c[0][0] + c[1][1] + c[2][2])
- (c[0][1] + c[1][2] + c[2][0])
+ 3 * (c[3][3] + c[4][4] + c[5][5])
) / float(15)
return [Kv, Gv]
@property
def compliance_tensor(self):
"""Get compliance."""
return np.linalg.inv(self.et_tensor)
@property
def reuss_modulus(self):
"""Get Reuss modulus."""
c = self.compliance_tensor
Kr = 1 / float(
(c[0][0] + c[1][1] + c[2][2]) + 2 * (c[0][1] + c[1][2] + c[2][0])
)
Gr = 15 / (
4 * (c[0][0] + c[1][1] + c[2][2])
- 4 * (c[0][1] + c[1][2] + c[2][0])
+ 3 * (c[3][3] + c[4][4] + c[5][5])
)
return [Kr, Gr]
[docs] def long_velocity(self, atoms=None):
"""Longitudinal velocity using Navier equation."""
# density = atoms.density
weight = float(atoms.composition.weight)
volume = atoms.volume
mass_density = 1.6605e3 * weight / volume
avg_mod = self.average_modulus
k_vrh = avg_mod[0]
g_vrh = avg_mod[1]
# 1e9:= GPa to Pascal (kg/ms^2)
vel = np.sqrt(1e9 * (k_vrh + 4.0 / 3.0 * g_vrh) / mass_density)
return vel
[docs] def trans_velocity(self, atoms=None):
"""Transverse velocity."""
avg_mod = self.average_modulus
g_vrh = avg_mod[1]
volume = atoms.volume
weight = float(atoms.composition.weight)
mass_density = 1.6605e3 * weight / volume
vel = np.sqrt(1e9 * g_vrh / mass_density)
return vel
[docs] def velocity_average(self, atoms=None):
"""Average velocity."""
vt = self.trans_velocity(atoms=atoms)
vl = self.long_velocity(atoms=atoms)
return 1.0 / (
np.cbrt(
(1.0 / 3.0) * (2.0 / (vt * vt * vt) + 1.0 / (vl * vl * vl))
)
)
[docs] def debye_temperature(self, atoms=None):
"""Debye temperature."""
const = 1.05457e-34 / 1.38065e-23 # (h/kb)
v0 = atoms.volume * 1e-30 / atoms.num_atoms
vl = self.long_velocity(atoms=atoms)
vt = self.trans_velocity(atoms=atoms)
vm = 3 ** (1.0 / 3.0) * (1 / vl ** 3 + 2 / vt ** 3) ** (-1.0 / 3.0)
theta = const * vm * (6 * np.pi ** 2 / v0) ** (1.0 / 3.0)
return theta
@property
def average_modulus(self):
"""Get average modulus."""
return (
np.array(self.voigt_modulus) + np.array(self.reuss_modulus)
) / 2
@property
def pugh_ratio_voigt(self):
"""Get Voigt Pugh ratio."""
Kv, Gv = self.voigt_modulus
return Gv / Kv
@property
def pettifor_criteria(self):
"""Get Pettifor criteria."""
c = self.et_tensor
return c[0][1] - c[3][3]
@property
def is_brittle(self):
"""Check if brittle material."""
return self.pugh_ratio_voigt > 0.571 and self.pettifor_criteria < 0
@property
def is_ductile(self):
"""Check if ductile material."""
return self.pugh_ratio_voigt < 0.571 and self.pettifor_criteria > 0
@property
def melting_temperature_metals(self):
"""Get crude Melting temp. estimate."""
# https://doi.org/10.1016/0036-9748(84)90267-9
avg_mod = self.average_modulus
k_vrh = avg_mod[0]
return 607 + 9.3 * k_vrh
@property
def cauchy_pressure(self):
"""Get Cauchy pressure."""
# >0 ionic bonding
# <0 covalent bonding
c = self.et_tensor
return c[0][1] - c[3][3]
@property
def poisson_ratio(self):
"""Get poisson's ratio."""
k, g = self.average_modulus
return (3 * k - 2 * g) / (6 * k + 2 * g)
@property
def universal_ansiotropy_ratio(self):
"""Get universal ansiotropy ratio."""
Kv, Gv = self.voigt_modulus
Kr, Gr = self.reuss_modulus
return 5 * (Gv / Gr) + (Kv / Kr) - 6
@property
def youngs_modulus(self):
"""Get Youngs modulus."""
k, g = self.average_modulus
return 9e9 * k * g / (3 * k + g)
[docs] def to_dict(self):
"""Get dictionary representation."""
d = OrderedDict()
d["voigt_bulk_modulus"] = self.voigt_modulus[0]
d["voigt_shear_modulus"] = self.voigt_modulus[1]
d["reuss_bulk_modulus"] = self.reuss_modulus[0]
d["reuss_shear_modulus"] = self.reuss_modulus[1]
d["poisson_ratio"] = self.poisson_ratio
d["youngs_modulus"] = self.youngs_modulus
d["universal_ansiotropy_ratio"] = self.universal_ansiotropy_ratio
if not isinstance(self.et_tensor, list):
et_tensor = self.et_tensor.tolist()
else:
et_tensor = self.et_tensor
d["raw_et_tensor"] = et_tensor
return d
"""
from jarvis.io.vasp.outputs import Vasprun,Outcar
o=Outcar('../../examples/vasp/SiOptb88/SiOptb88/MAIN-ELASTIC-bulk@mp_149/OUTCAR')
print (o.elastic_props())
p=Atoms.from_poscar('../../examples/vasp/SiOptb88/SiOptb88/MAIN-ELASTIC-bulk@mp_149/POSCAR')
et=ElasticTensor(o.elastic_props()['cij'])
#print (et.voigt_modulus) #[87.26666666666667, 63.28]
#print (et.reuss_modulus) #[87.26666666666665, 60.24546397096941]
#print (et.average_modulus) #[87.26666667 61.76273199]
#print (et.poisson_ratio) #0.21367500388646996
#print (et.universal_ansiotropy_ratio) #0.21367500388646996
#print ('j_vel',(1e9*(et.average_modulus[0])/p.density)**.5)
#print ('j_vel',(1e9*
(et.average_modulus[0]+4/3*et.average_modulus[1])/p.density)**.5)
k,g=et.average_modulus
mass_density=p.density
#print ((1e9 * (k + 4. / 3. * g) / mass_density/1.6605e3 ) ** 0.5)
print (et.to_dict())
from pymatgen.analysis.elasticity.elastic import ElasticTensor
et=ElasticTensor.from_voigt(o.elastic_props()['cij'])
from pymatgen.core.structure import Structure
pmg=Structure.from_file('../../examples/vasp/
SiOptb88/SiOptb88/MAIN-ELASTIC-bulk@mp_149/POSCAR')
#print (et.k_vrh,et.g_vrh)
#print (et.long_v(pmg))
#print ('density_j,density_p',p.density,pmg.density)
"""