# coding: utf-8
# Copyright (c) Max-Planck-Institut für Eisenforschung GmbH - Computational Materials Design (CM) Department
# Distributed under the terms of "New BSD License", see the LICENSE file.
from __future__ import print_function
import numpy as np
import scipy.integrate
import scipy.optimize as spy
import scipy.constants
import warnings
from pyiron.atomistics.structure.atoms import Atoms, ase_to_pyiron
from pyiron.atomistics.master.parallel import AtomisticParallelMaster
from pyiron.base.master.parallel import JobGenerator
__author__ = "Joerg Neugebauer, Jan Janssen"
__copyright__ = (
"Copyright 2020, Max-Planck-Institut für Eisenforschung GmbH - "
"Computational Materials Design (CM) Department"
)
__version__ = "1.0"
__maintainer__ = "Jan Janssen"
__email__ = "janssen@mpie.de"
__status__ = "production"
__date__ = "Sep 1, 2017"
eV_div_A3_to_GPa = (
1e21 / scipy.constants.physical_constants["joule-electron volt relationship"][0]
)
def _debye_kernel(xi):
return xi ** 3 / (np.exp(xi) - 1)
[docs]def debye_integral(x):
return scipy.integrate.quad(_debye_kernel, 0, x)[0]
[docs]def debye_function(x):
if hasattr(x, "__len__"):
return np.array([3 / xx ** 3 * debye_integral(xx) for xx in x])
return 3 / x ** 3 * debye_integral(x)
# https://gitlab.com/ase/ase/blob/master/ase/eos.py
[docs]def birchmurnaghan_energy(V, E0, B0, BP, V0):
"BirchMurnaghan equation from PRB 70, 224107"
eta = (V0 / V) ** (1 / 3)
return E0 + 9 * B0 * V0 / 16 * (eta ** 2 - 1) ** 2 * (
6 + BP * (eta ** 2 - 1) - 4 * eta ** 2
)
[docs]def vinet_energy(V, E0, B0, BP, V0):
"Vinet equation from PRB 70, 224107"
eta = (V / V0) ** (1 / 3)
return E0 + 2 * B0 * V0 / (BP - 1) ** 2 * (
2 - (5 + 3 * BP * (eta - 1) - 3 * eta) * np.exp(-3 * (BP - 1) * (eta - 1) / 2)
)
[docs]def murnaghan(V, E0, B0, BP, V0):
"From PRB 28,5480 (1983"
E = E0 + B0 * V / BP * (((V0 / V) ** BP) / (BP - 1) + 1) - V0 * B0 / (BP - 1)
return E
[docs]def birch(V, E0, B0, BP, V0):
"""
From Intermetallic compounds: Principles and Practice, Vol. I: Principles
Chapter 9 pages 195-210 by M. Mehl. B. Klein, D. Papaconstantopoulos
paper downloaded from Web
case where n=0
"""
E = (
E0
+ 9 / 8 * B0 * V0 * ((V0 / V) ** (2 / 3) - 1) ** 2
+ 9 / 16 * B0 * V0 * (BP - 4) * ((V0 / V) ** (2 / 3) - 1) ** 3
)
return E
[docs]def pouriertarantola(V, E0, B0, BP, V0):
"Pourier-Tarantola equation from PRB 70, 224107"
eta = (V / V0) ** (1 / 3)
squiggle = -3 * np.log(eta)
E = E0 + B0 * V0 * squiggle ** 2 / 6 * (3 + squiggle * (BP - 2))
return E
[docs]def fitfunction(parameters, vol, fittype="vinet"):
"""
Fit the energy volume curve
Args:
parameters (list): [E0, B0, BP, V0] list of fit parameters
vol (float/numpy.dnarray): single volume or a vector of volumes as numpy array
fittype (str): on of the following ['birch', 'birchmurnaghan', 'murnaghan', 'pouriertarantola', 'vinet']
Returns:
(float/numpy.dnarray): single energy as float or a vector of energies as numpy array
"""
[E0, b0, bp, V0] = parameters
# Unit correction
B0 = b0 / eV_div_A3_to_GPa
BP = bp
V = vol
if fittype.lower() == "birchmurnaghan":
return birchmurnaghan_energy(V, E0, B0, BP, V0)
elif fittype.lower() == "vinet":
return vinet_energy(V, E0, B0, BP, V0)
elif fittype.lower() == "murnaghan":
return murnaghan(V, E0, B0, BP, V0)
elif fittype.lower() == "pouriertarantola":
return pouriertarantola(V, E0, B0, BP, V0)
elif fittype.lower() == "birch":
return birch(V, E0, B0, BP, V0)
else:
raise ValueError
[docs]def fit_leastsq(p0, datax, datay, fittype="vinet"):
"""
Least square fit
Args:
p0 (list): [E0, B0, BP, V0] list of fit parameters
datax (float/numpy.dnarray): volumes to fit
datay (float/numpy.dnarray): energies corresponding to the volumes
fittype (str): on of the following ['birch', 'birchmurnaghan', 'murnaghan', 'pouriertarantola', 'vinet']
Returns:
list: [E0, B0, BP, V0], [E0_err, B0_err, BP_err, V0_err]
"""
# http://stackoverflow.com/questions/14581358/getting-standard-errors-on-fitted-parameters-using-the-optimize-leastsq-method-i
errfunc = lambda p, x, y, fittype: fitfunction(p, x, fittype) - y
pfit, pcov, infodict, errmsg, success = spy.leastsq(
errfunc, p0, args=(datax, datay, fittype), full_output=1, epsfcn=0.0001
)
if (len(datay) > len(p0)) and pcov is not None:
s_sq = (errfunc(pfit, datax, datay, fittype) ** 2).sum() / (
len(datay) - len(p0)
)
pcov = pcov * s_sq
else:
pcov = np.inf
error = []
for i in range(len(pfit)):
try:
error.append(np.absolute(pcov[i][i]) ** 0.5)
except:
error.append(0.00)
pfit_leastsq = pfit
perr_leastsq = np.array(error)
return pfit_leastsq, perr_leastsq
[docs]class DebyeModel(object):
"""
Calculate Thermodynamic Properties based on the Murnaghan output
"""
def __init__(self, murnaghan, num_steps=50):
self._murnaghan = murnaghan
# self._atoms_per_cell = len(murnaghan.structure)
self._v_min = None
self._v_max = None
self._num_steps = None
self._volume = None
self._init_volume()
self.num_steps = num_steps
self._fit_volume = None
self._debye_T = None
def _init_volume(self):
vol = self._murnaghan["output/volume"]
self._v_min, self._v_max = np.min(vol), np.max(vol)
def _set_volume(self):
if self._v_min and self._v_max and self._num_steps:
self._volume = np.linspace(self._v_min, self._v_max, self._num_steps)
self._reset()
# print ('set_volume: ', self._num_steps)
@property
def num_steps(self):
return self._num_steps
@num_steps.setter
def num_steps(self, val):
self._num_steps = val
self._set_volume()
@property
def volume(self):
if self._volume is None:
self._init_volume()
self._set_volume()
return self._volume
@volume.setter
def volume(self, volume_lst):
self._volume = volume_lst
self._v_min = np.min(volume_lst)
self._v_max = np.max(volume_lst)
self._reset()
def _reset(self):
self._debye_T = None
[docs] def polynomial(self, poly_fit=None, volumes=None):
if poly_fit is None:
self._murnaghan.fit_polynomial() # TODO: include polyfit in output
poly_fit = self._murnaghan.fit_dict["poly_fit"]
p_fit = np.poly1d(poly_fit)
if volumes is None:
return p_fit(self.volume)
return p_fit(volumes)
@property
def debye_temperature(self):
if self._debye_T is not None:
return self._debye_T
GPaTokBar = 10
Ang3_to_Bohr3 = (
scipy.constants.angstrom ** 3
/ scipy.constants.physical_constants["Bohr radius"][0] ** 3
)
convert = 67.48 # conversion factor, Moruzzi Eq. (4)
empirical = 0.617 # empirical factor, Moruzzi Eq. (6)
gamma_low, gamma_high = 1, 2 / 3 # low/high T gamma
out = self._murnaghan["output"]
V0 = out["equilibrium_volume"]
B0 = out["equilibrium_bulk_modulus"]
Bp = out["equilibrium_b_prime"]
vol = self.volume
mass = set(self._murnaghan.structure.get_masses())
if len(mass) > 1:
raise NotImplementedError(
"Debye temperature only for single species systems!"
)
mass = list(mass)[0]
r0 = (3 * V0 * Ang3_to_Bohr3 / (4 * np.pi)) ** (1.0 / 3.0)
debye_zero = empirical * convert * np.sqrt(r0 * B0 * GPaTokBar / mass)
# print('r0, B0, Bp, mass, V0', r0, B0, Bp, mass, V0)
# print('gamma_low, gamma_high: ', gamma_low, gamma_high)
# print('debye_zero, V0: ', debye_zero, V0)
if vol is None:
print("WARNING: vol: ", vol)
debye_low = debye_zero * (V0 / vol) ** (-gamma_low + 0.5 * (1 + Bp))
debye_high = debye_zero * (V0 / vol) ** (-gamma_high + 0.5 * (1 + Bp))
self._debye_T = (debye_low, debye_high)
return self._debye_T
[docs] def energy_vib(self, T, debye_T=None, low_T_limit=True):
kB = 0.086173422 / 1000 # eV/K
if debye_T is None:
if low_T_limit:
debye_T = self.debye_temperature[0] # low
else:
debye_T = self.debye_temperature[1] # high
if hasattr(debye_T, "__len__"):
val = [
9.0 / 8.0 * kB * d_T
+ T * kB * (3 * np.log(1 - np.exp(-d_T / T)) - debye_function(d_T / T))
for d_T in debye_T
]
val = np.array(val)
else:
val = 9.0 / 8.0 * kB * debye_T + T * kB * (
3 * np.log(1 - np.exp(-debye_T / T)) - debye_function(debye_T / T)
)
atoms_per_cell = len(self._murnaghan.structure)
return atoms_per_cell * val
[docs]class MurnaghanJobGenerator(JobGenerator):
@property
def parameter_list(self):
"""
Returns:
(list)
"""
parameter_lst = []
for strain in np.linspace(
1 - self._job.input["vol_range"],
1 + self._job.input["vol_range"],
self._job.input["num_points"],
):
basis = self._job.ref_job.structure.copy()
basis.set_cell(basis.cell * strain ** (1.0 / 3.0), scale_atoms=True)
parameter_lst.append([np.round(strain, 7), basis])
return parameter_lst
[docs] @staticmethod
def job_name(parameter):
return "strain_" + str(parameter[0]).replace(".", "_")
[docs] def modify_job(self, job, parameter):
job.structure = parameter[1]
return job
[docs]class EnergyVolumeFit(object):
"""
Fit energy volume curves
Args:
volume_lst (list/numpy.dnarray): vector of volumes
energy_lst (list/numpy.dnarray): vector of energies
Attributes:
.. attribute:: volume_lst
vector of volumes
.. attribute:: energy_lst
vector of energies
.. attribute:: fit_dict
dictionary of fit parameters
"""
def __init__(self, volume_lst=None, energy_lst=None):
self._volume_lst = volume_lst
self._energy_lst = energy_lst
self._fit_dict = None
@property
def volume_lst(self):
return self._volume_lst
@volume_lst.setter
def volume_lst(self, vol_lst):
self._volume_lst = vol_lst
@property
def energy_lst(self):
return self._energy_lst
@energy_lst.setter
def energy_lst(self, eng_lst):
self._energy_lst = eng_lst
@property
def fit_dict(self):
return self._fit_dict
def _get_volume_and_energy_lst(self, volume_lst=None, energy_lst=None):
"""
Internal function to get the vector of volumes and the vector of energies
Args:
volume_lst (list/numpy.dnarray/None): vector of volumes
energy_lst (list/numpy.dnarray/None): vector of energies
Returns:
list: vector of volumes and vector of energies
"""
if volume_lst is None:
if self._volume_lst is None:
raise ValueError("Volume list not set.")
volume_lst = self._volume_lst
if energy_lst is None:
if self._energy_lst is None:
raise ValueError("Volume list not set.")
energy_lst = self._energy_lst
return volume_lst, energy_lst
[docs] def fit_eos_general_intern(self, fittype="birchmurnaghan"):
self._fit_dict = self.fit_eos_general(
volume_lst=self._volume_lst, energy_lst=self._energy_lst, fittype=fittype
)
[docs] def fit_eos_general(
self, volume_lst=None, energy_lst=None, fittype="birchmurnaghan"
):
"""
Fit on of the equations of state
Args:
volume_lst (list/numpy.dnarray/None): vector of volumes
energy_lst (list/numpy.dnarray/None): vector of energies
fittype (str): on of the following ['birch', 'birchmurnaghan', 'murnaghan', 'pouriertarantola', 'vinet']
Returns:
dict: dictionary with fit results
"""
volume_lst, energy_lst = self._get_volume_and_energy_lst(
volume_lst=volume_lst, energy_lst=energy_lst
)
fit_dict = {}
pfit_leastsq, perr_leastsq = self._fit_leastsq(
volume_lst=volume_lst, energy_lst=energy_lst, fittype=fittype
)
fit_dict["fit_type"] = fittype
fit_dict["volume_eq"] = pfit_leastsq[3]
fit_dict["energy_eq"] = pfit_leastsq[0]
fit_dict["bulkmodul_eq"] = pfit_leastsq[1]
fit_dict["b_prime_eq"] = pfit_leastsq[2]
fit_dict["least_square_error"] = perr_leastsq # [e0, b0, bP, v0]
return fit_dict
[docs] def fit_polynomial(self, volume_lst=None, energy_lst=None, fit_order=3):
"""
Fit a polynomial
Args:
volume_lst (list/numpy.dnarray/None): vector of volumes
energy_lst (list/numpy.dnarray/None): vector of energies
fit_order (int): Degree of the polynomial
Returns:
dict: dictionary with fit results
"""
volume_lst, energy_lst = self._get_volume_and_energy_lst(
volume_lst=volume_lst, energy_lst=energy_lst
)
fit_dict = {}
# compute a polynomial fit
z = np.polyfit(volume_lst, energy_lst, fit_order)
p_fit = np.poly1d(z)
fit_dict["poly_fit"] = z
# get equilibrium lattice constant
# search for the local minimum with the lowest energy
p_deriv_1 = np.polyder(p_fit, 1)
roots = np.roots(p_deriv_1)
# volume_eq_lst = np.array([np.real(r) for r in roots if np.abs(np.imag(r)) < 1e-10])
volume_eq_lst = np.array(
[
np.real(r)
for r in roots
if (
abs(np.imag(r)) < 1e-10
and r >= min(volume_lst)
and r <= max(volume_lst)
)
]
)
e_eq_lst = p_fit(volume_eq_lst)
arg = np.argsort(e_eq_lst)
# print ("v_eq:", arg, e_eq_lst)
if len(e_eq_lst) == 0:
return None
e_eq = e_eq_lst[arg][0]
volume_eq = volume_eq_lst[arg][0]
# get bulk modulus at equ. lattice const.
p_2deriv = np.polyder(p_fit, 2)
p_3deriv = np.polyder(p_fit, 3)
a2 = p_2deriv(volume_eq)
a3 = p_3deriv(volume_eq)
b_prime = -(volume_eq * a3 / a2 + 1)
fit_dict["fit_type"] = "polynomial"
fit_dict["fit_order"] = fit_order
fit_dict["volume_eq"] = volume_eq
fit_dict["energy_eq"] = e_eq
fit_dict["bulkmodul_eq"] = eV_div_A3_to_GPa * volume_eq * a2
fit_dict["b_prime_eq"] = b_prime
fit_dict["least_square_error"] = self.get_error(volume_lst, energy_lst, p_fit)
return fit_dict
def _fit_leastsq(self, volume_lst, energy_lst, fittype="birchmurnaghan"):
"""
Internal helper function for the least square fit
Args:
volume_lst (list/numpy.dnarray/None): vector of volumes
energy_lst (list/numpy.dnarray/None): vector of energies
fittype (str): on of the following ['birch', 'birchmurnaghan', 'murnaghan', 'pouriertarantola', 'vinet']
Returns:
list: [E0, B0, BP, V0], [E0_err, B0_err, BP_err, V0_err]
"""
vol_lst = np.array(volume_lst).flatten()
eng_lst = np.array(energy_lst).flatten()
a, b, c = np.polyfit(vol_lst, eng_lst, 2)
v0 = -b / (2 * a)
pfit_leastsq, perr_leastsq = fit_leastsq(
[a * v0 ** 2 + b * v0 + c, 2 * a * v0 * eV_div_A3_to_GPa, 4, v0],
vol_lst,
eng_lst,
fittype,
)
return pfit_leastsq, perr_leastsq # [e0, b0, bP, v0]
[docs] @staticmethod
def get_error(x_lst, y_lst, p_fit):
"""
Args:
x_lst:
y_lst:
p_fit:
Returns:
numpy.dnarray
"""
y_fit_lst = np.array(p_fit(x_lst))
error_lst = (y_lst - y_fit_lst) ** 2
return np.mean(error_lst)
[docs] def fit_energy(self, volume_lst):
"""
Gives the energy value for the corresponding energy volume fit defined in the fit dictionary.
Args:
volume_lst: list of volumes
Returns:
list of energies
"""
if not self._fit_dict:
return ValueError("parameter 'fit_dict' has to be defined!")
v = volume_lst
e0 = self._fit_dict["energy_eq"]
b0 = self._fit_dict["bulkmodul_eq"] / eV_div_A3_to_GPa
b_p = self._fit_dict["b_prime_eq"]
v0 = self._fit_dict["volume_eq"]
if self._fit_dict["fit_type"] == "birch":
return self.birch(v, e0, b0, b_p, v0)
elif self._fit_dict["fit_type"] == "birchmurnaghan":
return self.birchmurnaghan_energy(v, e0, b0, b_p, v0)
elif self._fit_dict["fit_type"] == "murnaghan":
return self.murnaghan(v, e0, b0, b_p, v0)
elif self._fit_dict["fit_type"] == "pouriertarantola":
return self.pouriertarantola(v, e0, b0, b_p, v0)
else:
return self.vinet_energy(v, e0, b0, b_p, v0)
[docs] @staticmethod
def birchmurnaghan_energy(V, E0, B0, BP, V0):
"BirchMurnaghan equation from PRB 70, 224107"
return birchmurnaghan_energy(V, E0, B0, BP, V0)
[docs] @staticmethod
def vinet_energy(V, E0, B0, BP, V0):
"Vinet equation from PRB 70, 224107"
return vinet_energy(V, E0, B0, BP, V0)
[docs] @staticmethod
def murnaghan(V, E0, B0, BP, V0):
"From PRB 28,5480 (1983"
return murnaghan(V, E0, B0, BP, V0)
[docs] @staticmethod
def birch(V, E0, B0, BP, V0):
"""
From Intermetallic compounds: Principles and Practice, Vol. I: Principles
Chapter 9 pages 195-210 by M. Mehl. B. Klein, D. Papaconstantopoulos
paper downloaded from Web
case where n=0
"""
return birch(V, E0, B0, BP, V0)
[docs] @staticmethod
def pouriertarantola(V, E0, B0, BP, V0):
return pouriertarantola(V, E0, B0, BP, V0)
# ToDo: not all abstract methods implemented
[docs]class Murnaghan(AtomisticParallelMaster):
def __init__(self, project, job_name="murnaghan"):
"""
Args:
project:
job_name:
"""
super(Murnaghan, self).__init__(project, job_name)
self.__name__ = "Murnaghan"
self.__version__ = "0.3.0"
# print ("h5_path: ", self.project_hdf5._h5_path)
# define default input
self.input["num_points"] = (11, "number of sample points")
self.input["fit_type"] = (
"polynomial",
"['polynomial', 'birch', 'birchmurnaghan', 'murnaghan', 'pouriertarantola', 'vinet']",
)
self.input["fit_order"] = (3, "order of the fit polynom")
self.input["vol_range"] = (
0.1,
"relative volume variation around volume defined by ref_ham",
)
self.debye_model = DebyeModel(self)
self.fit_module = EnergyVolumeFit()
self.fit_dict = None
self._debye_T = None
self._job_generator = MurnaghanJobGenerator(self)
@property
def fit(self):
return self.debye_model
@property
def equilibrium_volume(self):
return self.fit_dict["volume_eq"]
@property
def equilibrium_energy(self):
return self.fit_dict["energy_eq"]
[docs] def fit_polynomial(self, fit_order=3, vol_erg_dic=None):
return self.poly_fit(fit_order=fit_order, vol_erg_dic=vol_erg_dic)
[docs] def fit_murnaghan(self, vol_erg_dic=None):
return self._fit_eos_general(vol_erg_dic=vol_erg_dic, fittype="murnaghan")
[docs] def fit_birch_murnaghan(self, vol_erg_dic=None):
return self._fit_eos_general(vol_erg_dic=vol_erg_dic, fittype="birchmurnaghan")
[docs] def fit_vinet(self, vol_erg_dic=None):
return self._fit_eos_general(vol_erg_dic=vol_erg_dic, fittype="vinet")
def _fit_eos_general(self, vol_erg_dic=None, fittype="birchmurnaghan"):
self._set_fit_module(vol_erg_dic=vol_erg_dic)
fit_dict = self.fit_module.fit_eos_general(fittype=fittype)
self.input["fit_type"] = fit_dict["fit_type"]
self.input["fit_order"] = 0
with self.project_hdf5.open("input") as hdf5_input:
self.input.to_hdf(hdf5_input)
with self.project_hdf5.open("output") as hdf5:
hdf5["equilibrium_energy"] = fit_dict["energy_eq"]
hdf5["equilibrium_volume"] = fit_dict["volume_eq"]
hdf5["equilibrium_bulk_modulus"] = fit_dict["bulkmodul_eq"]
hdf5["equilibrium_b_prime"] = fit_dict["b_prime_eq"]
self.fit_dict = fit_dict
return fit_dict
def _fit_leastsq(self, volume_lst, energy_lst, fittype="birchmurnaghan"):
return self.fit_module._fit_leastsq(
volume_lst=volume_lst, energy_lst=energy_lst, fittype=fittype
)
def _set_fit_module(self, vol_erg_dic=None):
if vol_erg_dic is not None:
if "volume" in vol_erg_dic.keys() and "energy" in vol_erg_dic.keys():
self.fit_module = EnergyVolumeFit(
volume_lst=vol_erg_dic["volume"], energy_lst=vol_erg_dic["energy"]
)
else:
raise KeyError
else:
df = self.output_to_pandas()
self.fit_module = EnergyVolumeFit(
volume_lst=df["volume"].values, energy_lst=df["energy"].values
)
[docs] def poly_fit(self, fit_order=3, vol_erg_dic=None):
self._set_fit_module(vol_erg_dic=vol_erg_dic)
fit_dict = self.fit_module.fit_polynomial(fit_order=fit_order)
if fit_dict is None:
self._logger.warning("Minimum could not be found!")
else:
self.input["fit_type"] = fit_dict["fit_type"]
self.input["fit_order"] = fit_dict["fit_order"]
with self.project_hdf5.open("input") as hdf5_input:
self.input.to_hdf(hdf5_input)
with self.project_hdf5.open("output") as hdf5:
hdf5["equilibrium_energy"] = fit_dict["energy_eq"]
hdf5["equilibrium_volume"] = fit_dict["volume_eq"]
hdf5["equilibrium_bulk_modulus"] = fit_dict["bulkmodul_eq"]
hdf5["equilibrium_b_prime"] = fit_dict["b_prime_eq"]
with self._hdf5.open("output") as hdf5:
structure = self.get_structure(iteration_step=-1)
if not isinstance(structure, Atoms):
structure = ase_to_pyiron(structure)
structure.to_hdf(hdf5)
self.fit_dict = fit_dict
return fit_dict
[docs] def list_structures(self):
if self.ref_job.structure is not None:
return [parameter[1] for parameter in self._job_generator.parameter_list]
else:
return []
[docs] def collect_output(self):
if self.server.run_mode.interactive:
ham = self.project_hdf5.inspect(self.child_ids[0])
erg_lst = ham["output/generic/energy_tot"]
vol_lst = ham["output/generic/volume"]
arg_lst = np.argsort(vol_lst)
self._output["volume"] = vol_lst[arg_lst]
self._output["energy"] = erg_lst[arg_lst]
else:
erg_lst, vol_lst, err_lst, id_lst = [], [], [], []
for job_id in self.child_ids:
ham = self.project_hdf5.inspect(job_id)
print("job_id: ", job_id, ham.status)
if "energy_tot" in ham["output/generic"].list_nodes():
energy = ham["output/generic/energy_tot"][-1]
elif "energy_pot" in ham["output/generic"].list_nodes():
energy = ham["output/generic/energy_pot"][-1]
else:
raise ValueError('Neither energy_pot or energy_tot was found.')
volume = ham["output/generic/volume"][-1]
erg_lst.append(np.mean(energy))
err_lst.append(np.var(energy))
vol_lst.append(volume)
id_lst.append(job_id)
vol_lst = np.array(vol_lst)
erg_lst = np.array(erg_lst)
err_lst = np.array(err_lst)
id_lst = np.array(id_lst)
arg_lst = np.argsort(vol_lst)
self._output["volume"] = vol_lst[arg_lst]
self._output["energy"] = erg_lst[arg_lst]
self._output["error"] = err_lst[arg_lst]
self._output["id"] = id_lst[arg_lst]
with self.project_hdf5.open("output") as hdf5_out:
for key, val in self._output.items():
hdf5_out[key] = val
if self.input["fit_type"] == "polynomial":
self.fit_polynomial(fit_order=self.input["fit_order"])
else:
self._fit_eos_general(fittype=self.input["fit_type"])
[docs] def plot(self, num_steps=100, plt_show=True):
try:
import matplotlib.pylab as plt
except ImportError:
import matplotlib.pyplot as plt
if not self.fit_dict:
if self.input["fit_type"] == "polynomial":
self.fit_polynomial(fit_order=self.input["fit_order"])
else:
self._fit_eos_general(fittype=self.input["fit_type"])
df = self.output_to_pandas()
vol_lst, erg_lst = df["volume"].values, df["energy"].values
x_i = np.linspace(np.min(vol_lst), np.max(vol_lst), num_steps)
color = "blue"
if self.fit_dict is not None:
if self.input["fit_type"] == "polynomial":
p_fit = np.poly1d(self.fit_dict["poly_fit"])
least_square_error = self.fit_module.get_error(vol_lst, erg_lst, p_fit)
plt.title("Murnaghan: error: " + str(least_square_error))
plt.plot(
x_i,
p_fit(x_i),
"-",
label=self.input["fit_type"],
color=color,
linewidth=3,
)
else:
V0 = self.fit_dict["volume_eq"]
E0 = self.fit_dict["energy_eq"]
B0 = self.fit_dict["bulkmodul_eq"]
BP = self.fit_dict["b_prime_eq"]
eng_fit_lst = fitfunction(parameters=[E0, B0, BP, V0],
vol=x_i,
fittype=self.input["fit_type"])
plt.plot(
x_i,
eng_fit_lst,
"-",
label=self.input["fit_type"],
color=color,
linewidth=3,
)
plt.plot(vol_lst, erg_lst, "x", color=color, markersize=20)
plt.legend()
plt.xlabel("Volume ($\AA^3$)")
plt.ylabel("energy (eV)")
if plt_show:
plt.show()
[docs] def get_structure(self, iteration_step=-1):
"""
Returns: Structure with equilibrium volume
"""
if not (self.structure is not None):
raise AssertionError()
if iteration_step == -1:
snapshot = self.structure.copy()
old_vol = snapshot.get_volume()
new_vol = self["output/equilibrium_volume"]
k = (new_vol / old_vol) ** (1.0 / 3.0)
new_cell = snapshot.cell * k
snapshot.set_cell(new_cell, scale_atoms=True)
return snapshot
elif iteration_step == 0:
return self.structure
else:
raise ValueError("iteration_step should be either 0 or -1.")