import numpy as np
from InterPhon.util import MatrixLike, AtomType, SelectIndex, FilePath, File, KptPath
from InterPhon.util import k_points, Symmetry2D
from InterPhon.core import UnitCell
from InterPhon.core import SuperCell
from InterPhon.core import PostArgument
from InterPhon.core import PreProcess
from InterPhon.inout import vasp, aims, espresso
[docs]class PostProcess(PreProcess):
"""
Post process class to control post-process.
This child class is inherited from the :class:`core.PreProcess` parent class.
The information required in post-process is stored in the instance variables of this class.
The instance variables are set by the :class:`core.PostProcess.set_user_arg`, :class:`core.PostProcess.set_reciprocal_lattice`,
:class:`core.PostProcess.set_force_constant`, and :class:`core.PostProcess.set_k_points` methods.
From the instance variables, eigen-frequency and eigen-mode are calculated using the :class:`core.PostProcess.eval_phonon` method,
which will be further analyzed by employing the modules in analysis sub-package.
:param in_file_unit_cell: Path of unit cell input file
:type in_file_unit_cell: str
:param in_file_super_cell: Path of super cell input file
:type in_file_super_cell: str
:param code_name: Specification of the file-format by a DFT program, defaults to vasp
:type code_name: str
:param user_arg: Instance of PostArgument class
:type user_arg: :class:`core.PostArgument`
:param unit_cell: Instance of UnitCell class
:type unit_cell: :class:`core.UnitCell`
:param super_cell: Instance of SuperCell class
:type super_cell: :class:`core.SuperCell`
"""
num_figure = 0
def __init__(self, in_file_unit_cell: FilePath,
in_file_super_cell: FilePath,
code_name: str = 'vasp',
user_arg=PostArgument(),
unit_cell=UnitCell(),
super_cell=SuperCell()):
"""
Constructor of PostProcess class.
"""
super(PostProcess, self).__init__(user_arg,
unit_cell,
super_cell)
self.unit_cell.read_unit_cell(in_file_unit_cell, code_name=code_name)
self.unit_cell.set_mass_true()
self.super_cell.read_unit_cell(in_file_super_cell, code_name=code_name)
self.reciprocal_matrix = np.empty((3, 3))
self.force_constant = np.empty((len(self.super_cell.atom_type) * 3,
len(self.unit_cell.atom_true) * 3))
self.k_points: KptPath = []
self.auto_k_points = []
self.dyn_matrix = np.empty((len(self.k_points),
len(self.unit_cell.xyz_true),
len(self.unit_cell.xyz_true)), dtype=complex)
self.w_q = np.empty((len(self.k_points),
len(self.unit_cell.xyz_true)), dtype=float)
self.v_q = np.empty((len(self.k_points),
len(self.unit_cell.xyz_true),
len(self.unit_cell.xyz_true)), dtype=complex)
[docs] def set_user_arg(self, dict_args: dict) -> None:
"""
Set a **PostArgument** instance from the information given by user.
:param dict_args: Argument dictionary given by user
:type dict_args: dict
"""
super(PostProcess, self).set_user_arg(dict_args)
self.user_arg.check_match_argument(self.unit_cell, self.super_cell)
self.super_cell.set_super_ind_true(self.unit_cell, self.user_arg)
[docs] def set_reciprocal_lattice(self) -> None:
"""
Set the instance variable (**self.reciprocal_matrix**).
"""
_volume = np.dot(self.unit_cell.lattice_matrix[0, 0:3],
np.cross(self.unit_cell.lattice_matrix[1, 0:3],
self.unit_cell.lattice_matrix[2, 0:3]))
for i in range(3):
self.reciprocal_matrix[i, 0:3] = 2 * np.pi * np.cross(self.unit_cell.lattice_matrix[(i + 1) % 3, 0:3],
self.unit_cell.lattice_matrix[(i + 2) % 3, 0:3]) / _volume
[docs] def set_force_constant(self, force_files: FilePath,
code_name: str = 'vasp',
sym_flag: bool = True) -> None:
"""
Set the instance variable (**self.force_constant**).
:param force_files: Path of DFT output files which contain atomic forces
:type force_files: str
:param code_name: Specification of the file-format by a DFT program, defaults to vasp
:type code_name: str
:param sym_flag: Specify whether to use symmetry operation, defaults to `True`
:type sym_flag: bool
"""
_ind_pbc = self.user_arg.periodicity.nonzero()[0]
if sym_flag:
if _ind_pbc.shape[0] != 2:
print('Caution:')
print('Current version supports symmetry functionality only for 2D periodic systems.')
print('"-sym" is changed from "{0}" to "False".'.format(sym_flag))
sym_flag = False
if code_name == 'vasp':
if sym_flag:
self.sym = Symmetry2D(self.unit_cell, self.super_cell, self.user_arg)
_, _, _ = self.sym.search_point_group()
_, _, _, _ = self.sym.search_image_atom()
self.sym.search_self_image_atom()
self.sym.search_independent_displacement()
self.sym.gen_additional_displacement()
force_ind = 0
for i, require in enumerate(self.sym.require_atom):
num_of_calculation = 2 * (len(self.sym.independent_additional_displacement_cart[i]) + 1)
_independent_displace = [self.sym.independent_by_single_displacement_cart[i][_v] for _v in range(len(self.sym.independent_by_single_displacement_cart[i]))]
_additional_displace = self.sym.independent_additional_displacement_cart[i]
for __v in range(len(_additional_displace)):
_independent_displace.append(_additional_displace[__v])
to_cart_displace = np.transpose(np.array(_independent_displace))
to_random_displace = np.linalg.inv(to_cart_displace)
forward_force = np.empty((len(self.super_cell.atom_type) * 3, 3))
backward_force = np.empty((len(self.super_cell.atom_type) * 3, 3))
for j in range(num_of_calculation):
if force_ind % 2 == 0:
_forward_matrix = vasp.read_output_lines(force_files[force_ind], len(self.super_cell.atom_type))
__forward_matrix = _forward_matrix.copy()
elif force_ind % 2 == 1:
_backward_matrix = vasp.read_output_lines(force_files[force_ind], len(self.super_cell.atom_type))
__backward_matrix = _backward_matrix.copy()
if j // 2 == 0:
for k, (W_index, W_cart) in enumerate(zip(self.sym.independent_by_W_index[i], self.sym.independent_by_W_displacement_cart[i]), start=0):
_forward_rot = W_cart @ np.transpose(_forward_matrix)
_backward_rot = W_cart @ np.transpose(_backward_matrix)
_image_sindex = [self.super_cell.atom_true[v] for i, v in enumerate(self.sym.same_supercell_index_select[W_index][require])]
_original_sindex = [self.super_cell.atom_true[i] for i, v in enumerate( self.sym.same_supercell_index_select[W_index][require])]
__forward_matrix[_image_sindex] = np.transpose(_forward_rot)[_original_sindex]
__backward_matrix[_image_sindex] = np.transpose(_backward_rot)[_original_sindex]
forward_force[:, j // 2 + k] = __forward_matrix.reshape([self.force_constant.shape[0], ])
backward_force[:, j // 2 + k] = __backward_matrix.reshape([self.force_constant.shape[0], ])
kk = k
else:
for _, _ in enumerate(self.sym.independent_additional_displacement_cart[i]):
forward_force[:, j // 2 + kk] = __forward_matrix.reshape([self.force_constant.shape[0], ])
backward_force[:, j // 2 + kk] = __backward_matrix.reshape([self.force_constant.shape[0], ])
force_ind += 1
_dif_force = - (forward_force - backward_force) @ to_random_displace / (2 * self.user_arg.displacement * 10 ** (-10))
self.force_constant[:, 3 * require: 3 * (require + 1)] = _dif_force
_original_basis = np.transpose(self.unit_cell.lattice_matrix.copy())
to_cart_coord = _original_basis / np.linalg.norm(_original_basis, axis=0)
to_direct_coord = np.linalg.inv(to_cart_coord)
for _point_group_ind, _not_require in zip(self.sym.point_group_ind, self.sym.not_require_atom):
W_in_cart = to_cart_coord @ self.sym.W_select[_point_group_ind] @ to_direct_coord
for _super_index, _super_same_index in enumerate(self.sym.same_supercell_index_select[_point_group_ind][_not_require]):
self.force_constant[3 * self.super_cell.atom_true[_super_index]: 3 * (self.super_cell.atom_true[_super_index] + 1),
3 * _not_require: 3 * (_not_require + 1)] \
= np.linalg.inv(W_in_cart) \
@ self.force_constant.copy()[3 * self.super_cell.atom_true[_super_same_index]: 3 * (self.super_cell.atom_true[_super_same_index] + 1),
3 * self.sym.same_index_select[_point_group_ind][0][_not_require]: 3 * (self.sym.same_index_select[_point_group_ind][0][_not_require] + 1)] \
@ W_in_cart
else:
for _ind_file, _force_file in enumerate(force_files):
if _ind_file % 2 == 0:
_forward_matrix = vasp.read_output_lines(_force_file, len(self.super_cell.atom_type))
elif _ind_file % 2 == 1:
_backward_matrix = vasp.read_output_lines(_force_file, len(self.super_cell.atom_type))
_dif_force = - (_forward_matrix - _backward_matrix) / (2 * self.user_arg.displacement * 10 ** (-10))
self.force_constant[:, _ind_file // 2] = _dif_force.reshape([self.force_constant.shape[0], ])
elif code_name == 'espresso':
if sym_flag:
self.sym = Symmetry2D(self.unit_cell, self.super_cell, self.user_arg)
_, _, _ = self.sym.search_point_group()
_, _, _, _ = self.sym.search_image_atom()
self.sym.search_self_image_atom()
self.sym.search_independent_displacement()
self.sym.gen_additional_displacement()
force_ind = 0
for i, require in enumerate(self.sym.require_atom):
num_of_calculation = 2 * (len(self.sym.independent_additional_displacement_cart[i]) + 1)
_independent_displace = [self.sym.independent_by_single_displacement_cart[i][_v] for _v in range(len(self.sym.independent_by_single_displacement_cart[i]))]
_additional_displace = self.sym.independent_additional_displacement_cart[i]
for __v in range(len(_additional_displace)):
_independent_displace.append(_additional_displace[__v])
to_cart_displace = np.transpose(np.array(_independent_displace))
to_random_displace = np.linalg.inv(to_cart_displace)
forward_force = np.empty((len(self.super_cell.atom_type) * 3, 3))
backward_force = np.empty((len(self.super_cell.atom_type) * 3, 3))
for j in range(num_of_calculation):
if force_ind % 2 == 0:
_forward_matrix = espresso.read_output_lines(force_files[force_ind], len(self.super_cell.atom_type))
__forward_matrix = _forward_matrix.copy()
elif force_ind % 2 == 1:
_backward_matrix = espresso.read_output_lines(force_files[force_ind], len(self.super_cell.atom_type))
__backward_matrix = _backward_matrix.copy()
if j // 2 == 0:
for k, (W_index, W_cart) in enumerate(zip(self.sym.independent_by_W_index[i], self.sym.independent_by_W_displacement_cart[i]), start=0):
_forward_rot = W_cart @ np.transpose(_forward_matrix)
_backward_rot = W_cart @ np.transpose(_backward_matrix)
_image_sindex = [self.super_cell.atom_true[v] for i, v in enumerate(self.sym.same_supercell_index_select[W_index][require])]
_original_sindex = [self.super_cell.atom_true[i] for i, v in enumerate(self.sym.same_supercell_index_select[W_index][require])]
__forward_matrix[_image_sindex] = np.transpose(_forward_rot)[_original_sindex]
__backward_matrix[_image_sindex] = np.transpose(_backward_rot)[_original_sindex]
forward_force[:, j // 2 + k] = __forward_matrix.reshape([self.force_constant.shape[0], ])
backward_force[:, j // 2 + k] = __backward_matrix.reshape([self.force_constant.shape[0], ])
kk = k
else:
for _, _ in enumerate(self.sym.independent_additional_displacement_cart[i]):
forward_force[:, j // 2 + kk] = __forward_matrix.reshape([self.force_constant.shape[0], ])
backward_force[:, j // 2 + kk] = __backward_matrix.reshape([self.force_constant.shape[0], ])
force_ind += 1
_dif_force = - (forward_force - backward_force) @ to_random_displace / (2 * self.user_arg.displacement * 10 ** (-10))
self.force_constant[:, 3 * require: 3 * (require + 1)] = _dif_force
_original_basis = np.transpose(self.unit_cell.lattice_matrix.copy())
to_cart_coord = _original_basis / np.linalg.norm(_original_basis, axis=0)
to_direct_coord = np.linalg.inv(to_cart_coord)
for _point_group_ind, _not_require in zip(self.sym.point_group_ind, self.sym.not_require_atom):
W_in_cart = to_cart_coord @ self.sym.W_select[_point_group_ind] @ to_direct_coord
for _super_index, _super_same_index in enumerate(self.sym.same_supercell_index_select[_point_group_ind][_not_require]):
self.force_constant[3 * self.super_cell.atom_true[_super_index]: 3 * (self.super_cell.atom_true[_super_index] + 1),
3 * _not_require: 3 * (_not_require + 1)] \
= np.linalg.inv(W_in_cart) \
@ self.force_constant.copy()[3 * self.super_cell.atom_true[_super_same_index]: 3 * (self.super_cell.atom_true[_super_same_index] + 1),
3 * self.sym.same_index_select[_point_group_ind][0][_not_require]: 3 * (self.sym.same_index_select[_point_group_ind][0][_not_require] + 1)] \
@ W_in_cart
else:
for _ind_file, _force_file in enumerate(force_files):
if _ind_file % 2 == 0:
_forward_matrix = espresso.read_output_lines(_force_file, len(self.super_cell.atom_type))
elif _ind_file % 2 == 1:
_backward_matrix = espresso.read_output_lines(_force_file, len(self.super_cell.atom_type))
_dif_force = - (_forward_matrix - _backward_matrix) / (2 * self.user_arg.displacement * 10 ** (-10))
self.force_constant[:, _ind_file // 2] = _dif_force.reshape([self.force_constant.shape[0], ])
elif code_name == 'aims':
if sym_flag:
self.sym = Symmetry2D(self.unit_cell, self.super_cell, self.user_arg)
_, _, _ = self.sym.search_point_group()
_, _, _, _ = self.sym.search_image_atom()
self.sym.search_self_image_atom()
self.sym.search_independent_displacement()
self.sym.gen_additional_displacement()
force_ind = 0
for i, require in enumerate(self.sym.require_atom):
num_of_calculation = 2 * (len(self.sym.independent_additional_displacement_cart[i]) + 1)
_independent_displace = [self.sym.independent_by_single_displacement_cart[i][_v] for _v in range(len(self.sym.independent_by_single_displacement_cart[i]))]
_additional_displace = self.sym.independent_additional_displacement_cart[i]
for __v in range(len(_additional_displace)):
_independent_displace.append(_additional_displace[__v])
to_cart_displace = np.transpose(np.array(_independent_displace))
to_random_displace = np.linalg.inv(to_cart_displace)
forward_force = np.empty((len(self.super_cell.atom_type) * 3, 3))
backward_force = np.empty((len(self.super_cell.atom_type) * 3, 3))
for j in range(num_of_calculation):
if force_ind % 2 == 0:
_forward_matrix = aims.read_output_lines(force_files[force_ind], len(self.super_cell.atom_type))
__forward_matrix = _forward_matrix.copy()
elif force_ind % 2 == 1:
_backward_matrix = aims.read_output_lines(force_files[force_ind], len(self.super_cell.atom_type))
__backward_matrix = _backward_matrix.copy()
if j // 2 == 0:
for k, (W_index, W_cart) in enumerate(zip(self.sym.independent_by_W_index[i], self.sym.independent_by_W_displacement_cart[i]), start=0):
_forward_rot = W_cart @ np.transpose(_forward_matrix)
_backward_rot = W_cart @ np.transpose(_backward_matrix)
_image_sindex = [self.super_cell.atom_true[v] for i, v in enumerate(self.sym.same_supercell_index_select[W_index][require])]
_original_sindex = [self.super_cell.atom_true[i] for i, v in enumerate(self.sym.same_supercell_index_select[W_index][require])]
__forward_matrix[_image_sindex] = np.transpose(_forward_rot)[_original_sindex]
__backward_matrix[_image_sindex] = np.transpose(_backward_rot)[_original_sindex]
forward_force[:, j // 2 + k] = __forward_matrix.reshape([self.force_constant.shape[0], ])
backward_force[:, j // 2 + k] = __backward_matrix.reshape([self.force_constant.shape[0], ])
kk = k
else:
for _, _ in enumerate(self.sym.independent_additional_displacement_cart[i]):
forward_force[:, j // 2 + kk] = __forward_matrix.reshape([self.force_constant.shape[0], ])
backward_force[:, j // 2 + kk] = __backward_matrix.reshape([self.force_constant.shape[0], ])
force_ind += 1
_dif_force = - (forward_force - backward_force) @ to_random_displace / (2 * self.user_arg.displacement * 10 ** (-10))
self.force_constant[:, 3 * require: 3 * (require + 1)] = _dif_force
_original_basis = np.transpose(self.unit_cell.lattice_matrix.copy())
to_cart_coord = _original_basis / np.linalg.norm(_original_basis, axis=0)
to_direct_coord = np.linalg.inv(to_cart_coord)
for _point_group_ind, _not_require in zip(self.sym.point_group_ind, self.sym.not_require_atom):
W_in_cart = to_cart_coord @ self.sym.W_select[_point_group_ind] @ to_direct_coord
for _super_index, _super_same_index in enumerate(self.sym.same_supercell_index_select[_point_group_ind][_not_require]):
self.force_constant[3 * self.super_cell.atom_true[_super_index]: 3 * (self.super_cell.atom_true[_super_index] + 1),
3 * _not_require: 3 * (_not_require + 1)] \
= np.linalg.inv(W_in_cart) \
@ self.force_constant.copy()[3 * self.super_cell.atom_true[_super_same_index]: 3 * (self.super_cell.atom_true[_super_same_index] + 1),
3 * self.sym.same_index_select[_point_group_ind][0][_not_require]: 3 * (self.sym.same_index_select[_point_group_ind][0][_not_require] + 1)] \
@ W_in_cart
else:
for _ind_file, _force_file in enumerate(force_files):
if _ind_file % 2 == 0:
_forward_matrix = aims.read_output_lines(_force_file, len(self.super_cell.atom_type))
elif _ind_file % 2 == 1:
_backward_matrix = aims.read_output_lines(_force_file, len(self.super_cell.atom_type))
_dif_force = - (_forward_matrix - _backward_matrix) / (2 * self.user_arg.displacement * 10 ** (-10))
self.force_constant[:, _ind_file // 2] = _dif_force.reshape([self.force_constant.shape[0], ])
[docs] def set_k_points(self, k_file: FilePath) -> None:
"""
Set the instance variable (**self.k_points**) by reading **KPOINTS** file given in VASP format.
The following four samplings of the wave vector k are available:
**1) Gamma-centered**
**2) Monkhorst-Pack**
**3) Line-path generation** (for band plot)
**4) Explicit designation**
For more details, see the following reference:
"The Basics of Electronic Structure Theory for Periodic Systems, Frontiers in chemistry 7, 1 (2019)".
:param k_file: Path of KPOINTS file
:type k_file: str
"""
try:
with open(k_file, 'r') as infile:
lines = infile.readlines()
except IOError:
print("\nFail to open '{0}' file".format(k_file))
print("Please check the path of k-points file\n")
raise
self.k_points: KptPath = []
self.auto_k_points = []
_ind_pbc = self.user_arg.periodicity.nonzero()[0]
if lines[2].split()[0][0] in ('G', 'g'):
# Gamma-centered grid of k-points
self.k_points, self.auto_k_points = k_points.gamma_centered(lines, _ind_pbc)
elif lines[2].split()[0][0] in ('M', 'm'):
# Monkhorst-Pack grid of k-points
self.k_points, self.auto_k_points = k_points.monkhorst_pack(lines, _ind_pbc)
elif lines[2].split()[0][0] in ('L', 'l'):
# Line-path of k-points
self.k_points = k_points.line_path(lines, _ind_pbc)
elif lines[2].split()[0][0] in ('R', 'r'):
# Explicit list of k-points
self.k_points = k_points.explicit_reciprocal(lines)
[docs] def eval_phonon(self) -> None:
"""
Set the instance variables, eigen-frequency (**self.w_q**) and corresponding eigen-mode (**self.v_q**).
"""
if len(self.k_points) != 0:
self.dyn_matrix = np.empty((len(self.k_points),
len(self.unit_cell.xyz_true),
len(self.unit_cell.xyz_true)), dtype=complex)
self.w_q = np.empty((len(self.k_points),
len(self.unit_cell.xyz_true)), dtype=float)
self.v_q = np.empty((len(self.k_points),
len(self.unit_cell.xyz_true),
len(self.unit_cell.xyz_true)), dtype=complex)
_ind_k = 0
_enlarge = 1
for ind, value in enumerate(self.user_arg.periodicity):
if value:
_enlarge = _enlarge * self.user_arg.enlargement[ind]
_ind_pbc = self.user_arg.periodicity.nonzero()[0]
for k_point in self.k_points:
q = np.dot(k_point, self.reciprocal_matrix)
_dyn_matrix = np.zeros([len(self.unit_cell.xyz_true), len(self.unit_cell.xyz_true)], dtype=complex)
for s1, satom_ind in enumerate(self.super_cell.xyz_true):
for s2, atom_ind in enumerate(self.unit_cell.xyz_true):
pos_vector = self.super_cell.atom_cart[satom_ind, 0:3] - self.unit_cell.atom_cart[atom_ind, 0:3]
if _ind_pbc.shape[0] == 0:
pass
elif _ind_pbc.shape[0] == 1:
for first in (-1, 0, 1):
_tmp_pos = pos_vector + first * self.super_cell.lattice_matrix[_ind_pbc[0], 0:3]
if np.dot(pos_vector, pos_vector) > np.dot(_tmp_pos, _tmp_pos):
pos_vector = _tmp_pos.copy()
elif _ind_pbc.shape[0] == 2:
for first in (-1, 0, 1):
for second in (-1, 0, 1):
_tmp_pos = pos_vector + first * self.super_cell.lattice_matrix[_ind_pbc[0], 0:3] \
+ second * self.super_cell.lattice_matrix[_ind_pbc[1], 0:3]
if np.dot(pos_vector, pos_vector) > np.dot(_tmp_pos, _tmp_pos):
pos_vector = _tmp_pos.copy()
elif _ind_pbc.shape[0] == 3:
for first in (-1, 0, 1):
for second in (-1, 0, 1):
for third in (-1, 0, 1):
_tmp_pos = pos_vector + first * self.super_cell.lattice_matrix[_ind_pbc[0], 0:3] \
+ second * self.super_cell.lattice_matrix[_ind_pbc[1], 0:3] \
+ third * self.super_cell.lattice_matrix[_ind_pbc[2], 0:3]
if np.dot(pos_vector, pos_vector) > np.dot(_tmp_pos, _tmp_pos):
pos_vector = _tmp_pos.copy()
_dyn_matrix[(s1 // (_enlarge * 3)) * 3 + (s1 % 3), s2] = \
_dyn_matrix[(s1 // (_enlarge * 3)) * 3 + (s1 % 3), s2] \
+ complex(self.force_constant[satom_ind * 3 + (s1 % 3), s2], 0) \
* np.exp(1j * complex(np.dot(q, pos_vector), 0))
_mass_true = np.sqrt(np.dot(self.unit_cell.mass_true.reshape([len(self.unit_cell.xyz_true), 1]),
self.unit_cell.mass_true.reshape([1, len(self.unit_cell.xyz_true)])))
self.dyn_matrix[_ind_k, :, :] = _dyn_matrix / _mass_true
_eig_w, self.v_q[_ind_k, :, :] = np.linalg.eig(self.dyn_matrix[_ind_k, :, :])
self.w_q[_ind_k, :] = (np.sqrt(_eig_w).real - np.abs(np.sqrt(_eig_w).imag)) / (2 * np.pi) / 10 ** 12 # THz
self.w_q[_ind_k, :] = self.w_q[_ind_k, np.argsort(self.w_q[_ind_k, :])]
self.v_q[_ind_k, :, :] = self.v_q[_ind_k, np.argsort(self.w_q[_ind_k, :]), :]
_ind_k = _ind_k + 1
if _ind_pbc.shape[0] == 0:
with open('./freqency_at_gamma_point.dat', 'w') as outfile:
comment = "Discrete frequency of non-periodic system"
outfile.write("%s" % comment + '\n')
outfile.write("%s" % ' Frequency (THz)' + '\n')
for _freq in self.w_q[0]:
line = ' %16.9f ' % _freq
outfile.write("%s" % line + '\n')
# def write_dos(self, out_file: FilePath = 'total_dos.dat', sigma: float = 0.0, num_dos: int = 200,
# partial_dos: bool = False, plot: bool = False) -> File:
# """
# Method of PostProcess class.
# Process to write and plot phonon DOS.
#
# usage:
# " >>> instance_of_PostProcess.write_dos(out_file='total_dos.dat', sigma=0.1, num_dos=1000,
# partial_dos=True, plot=True)"
#
# :param out_file: (str) Name of output file for total dos.
# :param sigma: (float) Value of sigma for gaussian smearing,
# if the sigma value is 0.0, DOS will be presented by the linear tetrahedron method.
# :param num_dos: (int) The number of energy points evaluated for DOS.
# :param partial_dos: (bool) Write (True) output file for partial dos or not (False).
# :param plot: (bool) Plot (True) phonon DOS or not (False).
# :return: (File)
# """
# minimum_freq = np.min(self.w_q)
# maximum_freq = np.max(self.w_q)
# _freq = np.arange(minimum_freq - 2, maximum_freq + 2, (maximum_freq - minimum_freq + 4) / num_dos)
# _pdos = np.zeros((len(self.unit_cell.xyz_true), _freq.shape[0]))
#
# if sigma == 0.0:
# # Linear Tetrahedron Method for Brillouin zone integration
# _ind_pbc = self.user_arg.periodicity.nonzero()[0]
# if _ind_pbc.shape[0] == 0:
# _freq = self.w_q[0, :]
# _tdos = np.ones(_freq.shape[0])
#
# elif _ind_pbc.shape[0] == 1:
# _pdos = tetrahedron_1d(_freq, _pdos, self.k_points, self.w_q, self.v_q, self.auto_k_points, _ind_pbc)
# _tdos = _pdos.sum(axis=0)
#
# elif _ind_pbc.shape[0] == 2:
# _pdos = tetrahedron_2d(_freq, _pdos, self.k_points, self.w_q, self.v_q, self.auto_k_points, _ind_pbc)
# _tdos = _pdos.sum(axis=0)
#
# elif _ind_pbc.shape[0] == 3:
# _pdos = tetrahedron_3d(_freq, _pdos, self.k_points, self.w_q, self.v_q, self.auto_k_points, _ind_pbc)
# _tdos = _pdos.sum(axis=0)
#
# else:
# # Gaussian Smearing Method for Brillouin zone integration
# for eig_freqs, eig_modes in zip(self.w_q, self.v_q):
# for ind_freq, eig_freq in enumerate(eig_freqs):
# for ind_mode, _ in enumerate(self.unit_cell.xyz_true):
# _pdos[ind_mode, :] = _pdos[ind_mode, :] \
# + 1 / (sigma * np.sqrt(2 * np.pi)) \
# * np.exp(- (_freq - eig_freq) ** 2 / (2 * sigma ** 2)) \
# / len(self.k_points) \
# * (abs(eig_modes[ind_freq, ind_mode]) ** 2)
#
# _tdos = _pdos.sum(axis=0)
#
# with open(out_file + '/total_dos.dat', 'w') as outfile:
# if sigma == 0.0:
# comment = "Total phonon DOS by Linear Tetrahedron Method"
# outfile.write("%s" % comment + '\n')
#
# elif sigma != 0.0:
# comment = "Total phonon DOS by Gaussian Smearing with Sigma = %f" % sigma
# outfile.write("%s" % comment + '\n')
#
# outfile.write("%s" %
# ' Frequency (THz)' +
# ' Total_Density_of_State' + '\n')
#
# for x, y in zip(_freq, _tdos):
# line = ' %16.9f ' % x + ' %16.9f ' % y
# outfile.write("%s" % line + '\n')
#
# if partial_dos is True:
# with open(out_file + '/partial_dos.dat', 'w') as outfile:
# if sigma == 0.0:
# comment = "Partial phonon DOS by Linear Tetrahedron Method"
# outfile.write("%s" % comment + '\n')
#
# elif sigma != 0.0:
# comment = "Partial phonon DOS by Gaussian Smearing with Sigma = %f" % sigma
# outfile.write("%s" % comment + '\n')
#
# line = "%s" % ' Frequency (THz)'
# for xyz_true in self.unit_cell.xyz_true:
# line += ' pdos-atom-{0:0>3} '.format(xyz_true)
# outfile.write(line + '\n')
#
# for x_ind, x in enumerate(_freq):
# line = ' %16.9f ' % x
# for ind, _ in enumerate(self.unit_cell.xyz_true):
# line += ' %16.9f ' % _pdos[ind, x_ind]
# outfile.write(line + '\n')
#
# if plot is True:
# from matplotlib import pyplot as plt
# self.num_figure = self.num_figure + 1
# plt.figure(self.num_figure)
# plt.plot(_freq, _tdos, color='black', label='total dos')
# for ind, xyz_true in enumerate(self.unit_cell.xyz_true):
# plt.plot(_freq, _pdos[ind], linestyle='--', label='atom-{0:0>4}'.format(xyz_true))
# plt.xlabel('Frequency (THz)')
# plt.ylabel('DOS (a.u.)')
# plt.legend(fontsize='x-large')
# plt.show()
#
# elif plot is True:
# from matplotlib import pyplot as plt
# self.num_figure = self.num_figure + 1
# plt.figure(self.num_figure)
# plt.plot(_freq, _tdos, color='black', label='total dos')
# plt.xlabel('Frequency (THz)')
# plt.ylabel('DOS (a.u.)')
# plt.legend(fontsize='x-large')
# plt.show()
#
# def write_band(self, out_file: FilePath = 'band.dat', plot: bool = False) -> File:
# """
# Method of PostProcess class.
# Process to write and plot phonon band.
#
# usage:
# " >>> instance_of_PostProcess.write_band(out_file='band.dat', plot=True)"
#
# :param out_file: (str) Name of output file for phonon band.
# :param plot: (bool) Plot (True) phonon band or not (False).
# :return: (File)
# """
#
# with open(out_file, 'w') as outfile:
# comment = "Phonon Band"
# outfile.write("%s" % comment + '\n')
# outfile.write("%s" %
# ' K_Points_Path' +
# ' Frequency (THz)' + '\n')
#
# k_points_length = np.zeros(len(self.k_points))
# for ind, kpt in enumerate(self.k_points[1:], 1):
# k_points_length[ind] = k_points_length[ind - 1] \
# + np.sqrt(np.dot(self.k_points[ind] - self.k_points[ind - 1],
# self.k_points[ind] - self.k_points[ind - 1]))
#
# for x, y_set in zip(k_points_length, self.w_q):
# line = ' %16.9f ' % x
# for y in y_set:
# line = line + ' %16.9f ' % y
# outfile.write("%s" % line + '\n')
#
# if plot is True:
# from matplotlib import pyplot as plt
# self.num_figure = self.num_figure + 1
# plt.figure(self.num_figure)
# plt.plot(k_points_length, self.w_q, color='black')
# plt.xlabel('K points')
# plt.ylabel('Frequency (THz)')
# plt.show()
#
# def write_mode(self, out_file: FilePath = 'mode',
# mode_inds: SelectIndex = [0],
# k_point: MatrixLike = np.array([0.0, 0.0, 0.0]),
# num_images: int = 30,
# plot: bool = False) -> File:
# """
# Method of PostProcess class.
# Process to write and plot phonon mode.
# Plotting phonon mode is supported with the help of Atomic Simulation Environment (ASE) package.
# (see, https://wiki.fysik.dtu.dk/ase/index.html)
# Therefore, in order to see phonon mode, the ASE package must be installed and
# the position of the package has to be automatically found by the Python interpreter (plot=True).
#
# usage:
# " >>> instance_of_PostProcess.write_mode(out_file='mode', mode_inds=[0,1],
# k_point=[0.0, 0.0, 0.0], num_images=10, plot=True)"
#
# :param out_file: (str) Name of output file for phonon mode.
# :param mode_inds: (List[int]) Mode index where '0' indicates the mode of the lowest frequency.
# :param k_point: (np.ndarray[float]) '(3,) size' k-point where corresponding phonon mode will be evaluated.
# :param num_images: (int) The number of images to see the oscillating phonon mode.
# :param plot: (bool) Plot (True) phonon mode or not (False).
# :return: (File)
# """
# k_point = np.array(k_point)
#
# _ind_k = None
# for ind, kpt in enumerate(self.k_points):
# if np.allclose(k_point, kpt):
# _ind_k = ind
# if _ind_k is None:
# raise error.Not_Specified_Kpath_Error
#
# mode = self.v_q[_ind_k, :, :]
#
# _current_position_true = self.unit_cell.atom_cart.copy()[self.unit_cell.atom_true]
# for mode_ind in mode_inds:
# for ind, x in enumerate(np.linspace(0, 2 * np.pi, num_images, endpoint=False), 1):
# self.unit_cell.atom_cart[self.unit_cell.atom_true] = _current_position_true + \
# np.sin(x) * mode[mode_ind, :].reshape((-1, 3)).real
#
# if ind == 1:
# lines = vasp.write_input_lines(self.unit_cell, "unknown system")
# lines[7] = "Cartesian configuration= %4d" % ind + '\n'
# else:
# line = vasp.write_input_lines(self.unit_cell, "unknown system")
# line[7] = "Cartesian configuration= %4d" % ind + '\n'
# lines.extend(line[7:])
# with open("XDATCAR" + "_%s_%d" % (out_file, mode_ind), 'w') as outfile:
# outfile.write("%s" % "".join(lines))
#
# if plot is True:
# try:
# from ase.io.trajectory import Trajectory
# from ase.io import read, iread
# from ase.visualize import view
# except ImportError:
# print("\nThe parent directory of ase package must be included in 'sys.path'")
# print("Example: 'sys.path.append('C:\\Users\\user\\Documents')")
# print("if the parent directory of ase package is 'C:\\Users\\user\\Documents'\n")
# raise
#
# atom = read('POSCAR')
# _current_position_true = atom.positions.copy()[self.unit_cell.atom_true]
# for mode_ind in mode_inds:
# traj = Trajectory("%s.%d.traj" % (out_file, mode_ind), 'w')
# for ind, x in enumerate(np.linspace(0, 2 * np.pi, num_images, endpoint=False), 1):
# atom.positions[self.unit_cell.atom_true] = _current_position_true + \
# np.sin(x) * mode[mode_ind, :].reshape((-1, 3)).real
# traj.write(atom)
# traj.close()
# atoms = iread("%s.%d.traj" % (out_file, mode_ind))
# view(atoms)
#
# def write_thermal_properties(self, out_file: FilePath = 'thermal_properties.dat',
# temp: MatrixLike = np.arange(0, 1001, 10),
# plot: bool = False) -> File:
# """
# Method of PostProcess class.
# Process to write and plot thermal properties.
#
# usage:
# " >>> instance_of_PostProcess.write_thermal_properties(out_file='thermal_properties.dat',
# temp=np.arange(0, 1201, 10), plot=True)"
#
# :param out_file: (str) Name of output file for thermal properties.
# :param temp: (np.ndarray) Temperature array where thermal properties will be evaluated.
# :param plot: (bool) Plot (True) thermal properties or not (False).
# :return: (File)
# """
# with open(out_file, 'w') as outfile:
# comment = "Thermal Properties / atom"
# outfile.write("%s" % comment + '\n')
# outfile.write("%s" %
# ' Temperature (K)' +
# ' Free_energy (eV/atom)' +
# ' Entropy (meV/K/atom)' +
# ' Heat_capacity (meV/K/atom)' + '\n')
#
# kb = 1.38 * 10 ** (-23) / (1.602 * 10 ** (-19))
# h = 6.626 * 10 ** (-34) / (1.602 * 10 ** (-19))
#
# temp = np.array(temp)
# free_energy = np.zeros(temp.shape)
# entropy = np.zeros(temp.shape)
# heat_capacity = np.zeros(temp.shape)
#
# for eig_freqs in self.w_q:
# for eig_freq in eig_freqs:
# for ind, t in enumerate(temp):
# if t < 10**(-6):
# free_energy[ind] = free_energy[ind] + (1.0 / 2.0 * h * eig_freq) / len(self.k_points) / (
# self.w_q.shape[1] / 3)
# else:
# free_energy[ind] = free_energy[ind] + (1.0 / 2.0 * h * eig_freq
# + kb * t * np.log(1 - np.exp(-h * eig_freq / (kb * t)))) / \
# len(self.k_points) / (self.w_q.shape[1] / 3)
# entropy[ind] = entropy[ind] + (-kb * np.log(1 - np.exp(-h * eig_freq / (kb * t)))
# + 1 / t * h * eig_freq / (np.exp(h * eig_freq / (kb * t)) - 1)) / \
# len(self.k_points) / (self.w_q.shape[1] / 3)
# heat_capacity[ind] = heat_capacity[ind] + kb * np.power(h * eig_freq / (kb * t), 2) \
# * np.exp(h * eig_freq / (kb * t)) / \
# np.power(np.exp(h * eig_freq / (kb * t)) - 1, 2) / \
# len(self.k_points) / (self.w_q.shape[1] / 3)
#
# for x, y1, y2, y3 in zip(temp, free_energy, entropy, heat_capacity):
# line = ' %20.9f ' % x + ' %20.9f ' % y1 + ' %20.9f ' % (y2 * 1000) + ' %20.9f ' % (y3 * 1000)
# outfile.write("%s" % line + '\n')
#
# if plot is True:
# from matplotlib import pyplot as plt
# self.num_figure = self.num_figure + 1
# plt.figure(self.num_figure)
# plt.plot(temp, free_energy, label='Free energy (eV/atom)')
# plt.plot(temp, entropy * 1000, label='Entropy (meV/K/atom)')
# plt.plot(temp, heat_capacity * 1000, label='Heat capacity (meV/K/atom)')
# plt.xlabel('Temperature (K)')
# plt.ylabel('Thermal Properties')
# plt.legend(fontsize='x-large')
# plt.show()