#!/usr/bin/env python
import numpy as np
import argparse
import dynaphopy.interface.iofile as reading
from phonopy import PhonopyQHA
from phonopy.file_IO import read_v_e, write_FORCE_CONSTANTS
from dynaphopy.interface.phonopy_link import get_force_constants_from_file
import yaml
parser = argparse.ArgumentParser(description='qha_quasiparticles options')
parser.add_argument('input_file', metavar='data_file', type=str,
help='input file containing structure related data')
parser.add_argument('-cv_data', metavar='data_files', type=str, nargs='*', required=True,
help='quasiparticle data at constant volume')
parser.add_argument('-ct_data', metavar='data_files', type=str, nargs='*', required=True,
help='quasiparticle data at constant temperature')
parser.add_argument('-temperatures', metavar='temperatures', type=float, nargs='*', required=True,
help='list of temperatures')
parser.add_argument('--ref_temperature', metavar='temperature', type=float, default=None,
help='Temperature at volume expansion')
parser.add_argument('--order', metavar='order', type=int, default=2,
help='fitting polynomial order')
parser.add_argument('-ev', metavar='data', type=str, required=True,
help='Unit cell volume vs energy file in ang.^3 and eV')
parser.add_argument('--pressure', metavar='p', type=float, required=False, default=0.0,
help='external pressure in GPa (default: 0 GPa)')
args = parser.parse_args()
class ForceConstantsFitting():
def __init__(self,
structure,
files_volume,
files_temperature,
temperatures=None,
volumes=None,
mesh=(40, 40, 40),
ref_temperature=None,
fitting_order=2,
tmin=None,
tmax=None,
use_NAC=False):
self.structure = structure
self.files_volume = files_volume
self.files_temperature = files_temperature
self.supercell = structure.get_supercell_phonon()
self._mesh = mesh
self.ref_temperature = ref_temperature
self.fitting_order = fitting_order
self._temperatures = temperatures
self._volumes = volumes
self._tmin = tmin
self._tmax = tmax
self._shift_temperature = None
self._shift_volume = None
self._shift_matrix_temperature = None
self._shift_matrix_volume = None
self._eigenvectors = None
self._nac = use_NAC # Use_Nac requires a BORN named file (phonopy style) in the work directory
# Properties
@property
def temperatures(self):
if self._temperatures is None:
return range(len(self.files_temperature))
else:
return self._temperatures
@property
def volumes(self):
if self._volumes is None:
return range(len(self.files_volume))
else:
return self._volumes
def get_temperature_range(self, step=10):
if self._tmin is None:
self._tmin = self.temperatures[0]
if self._tmax is None:
self._tmax = self.temperatures[-1]
return np.arange(self._tmin, self._tmax, step)
def get_shift_matrix_temperature(self):
if self._shift_matrix_temperature is None:
h_frequencies, ev = self.get_harmonic_frequencies_and_eigenvectors()
shift_matrix = []
list_t = self.temperatures
for i, t in enumerate(list_t):
with open(self.files_temperature[i], 'r') as stream:
data = yaml.load(stream)
renormalized_frequencies = []
for j, qpoint in enumerate(data):
renormalized_frequencies.append(qpoint['frequencies'])
shift_matrix.append(np.array(renormalized_frequencies) - np.array(h_frequencies))
self._shift_matrix_temperature = np.array(shift_matrix).swapaxes(0, 2)
return self._shift_matrix_temperature
def get_shift_matrix_volume(self):
if self._shift_matrix_volume is None:
#ref_data = self.get_reference_temperature_data()
h_frequencies, ev = self.get_harmonic_frequencies_and_eigenvectors()
ref_frequencies = self.get_fitted_shifts_temperature(temperature=self.ref_temperature) + h_frequencies
shift_matrix = []
list_v = self.volumes
for i, t in enumerate(list_v):
with open(self.files_volume[i], 'r') as stream:
data = yaml.load(stream)
renormalized_frequencies = []
for qpoint in data:
renormalized_frequencies.append(qpoint['frequencies'])
shift_matrix.append(np.array(renormalized_frequencies) - ref_frequencies)
self._shift_matrix_volume = np.array(shift_matrix).swapaxes(0, 2)
return self._shift_matrix_volume
def get_interpolated_shifts_temperature(self, temperature, kind='quadratic'):
from scipy.interpolate import interp1d
shift_matrix = self.get_shift_matrix_temperature()
shift_temperature = interp1d(self.temperatures, shift_matrix, kind=kind)
interpolated_shifts = shift_temperature(temperature).T
return interpolated_shifts
def get_fitted_shifts_temperature(self, temperature):
shift_matrix = self.get_shift_matrix_temperature()
if self._shift_temperature is None:
p = []
for i, row in enumerate(shift_matrix):
p2 = []
for j, r in enumerate(row):
p2.append(np.polyfit(self.temperatures, r, self.fitting_order))
p.append(p2)
self._shift_temperature = p
interpolated_shifts = []
for p in self._shift_temperature:
row = []
for p2 in p:
row.append(np.poly1d(p2)(temperature))
interpolated_shifts.append(row)
interpolated_shifts = np.array(interpolated_shifts).T
return interpolated_shifts
def get_interpolated_shifts_volume(self, volume, kind='quadratic'):
from scipy.interpolate import interp1d
shift_matrix = self.get_shift_matrix_volume()
shift_volume = interp1d(self.volumes, shift_matrix, kind=kind)
interpolated_shifts = shift_volume(volume).T
return interpolated_shifts
def get_fitted_shifts_volume(self, volume):
shift_matrix = self.get_shift_matrix_volume()
if self._shift_volume is None:
p = []
for i, row in enumerate(shift_matrix):
p2 = []
for j, r in enumerate(row):
p2.append(np.polyfit(self.volumes, r, self.fitting_order))
p.append(p2)
self._shift_volume = p
interpolated_shifts = []
for p in self._shift_volume:
row = []
for p2 in p:
row.append(np.poly1d(p2)(volume))
interpolated_shifts.append(row)
interpolated_shifts = np.array(interpolated_shifts).T
return interpolated_shifts
def get_harmonic_frequencies_and_eigenvectors(self):
from dynaphopy.interface.phonopy_link import obtain_eigenvectors_and_frequencies
if self._eigenvectors is None:
data = self.get_data_temperature(0)
com_ev = []
com_freq = []
for qpoint in data:
arranged_ev, frequencies = obtain_eigenvectors_and_frequencies(self.structure,
qpoint['reduced_wave_vector'],
print_data=False)
com_ev.append(arranged_ev)
com_freq.append(frequencies)
self._eigenvectors = com_ev
self._h_frequencies = com_freq
return self._h_frequencies, self._eigenvectors
def get_data_temperature(self, index):
with open(self.files_temperature[index], 'r') as stream:
data = yaml.load(stream)
return data
def get_data_volume(self, index):
with open(self.files_volume[index], 'r') as stream:
data = yaml.load(stream)
return data
def _get_renormalized_force_constants(self, renormalized_frequencies):
h_frequencies, eigenvectors = self.get_harmonic_frequencies_and_eigenvectors()
from dynaphopy.interface.phonopy_link import get_renormalized_force_constants
fc = get_renormalized_force_constants(renormalized_frequencies,
eigenvectors, self.structure,
self.supercell,
symmetrize=False)
return fc
def get_total_shifts(self, volume, temperature, interpolate_temperature=False, interpolate_volume=False):
if volume is None:
volume = self.volumes[0]
if interpolate_temperature:
temperature_shifts = self.get_interpolated_shifts_temperature(temperature)
else:
temperature_shifts = self.get_fitted_shifts_temperature(temperature)
if interpolate_volume:
volume_shifts = self.get_fitted_shifts_volume(volume)
else:
volume_shifts = self.get_interpolated_shifts_volume(volume)
return temperature_shifts + volume_shifts
def get_total_force_constants(self, temperature=300, volume=None):
total_shifts = self.get_total_shifts(volume, temperature)
h_frequencies, ev = self.get_harmonic_frequencies_and_eigenvectors()
total_frequency = h_frequencies + total_shifts
return self._get_renormalized_force_constants(total_frequency)
def get_thermal_properties(self, volume=None):
from dynaphopy.interface.phonopy_link import obtain_phonopy_thermal_properties
if volume is None:
volume = self.volumes[0]
print ('use volume: {}'.format(volume))
free_energy_list = []
entropy_list = []
cv_list = []
print ('temperature free energy(KJ/K/mol) entropy(KJ/K/mol) cv (J/K/mol)')
for t in self.get_temperature_range():
fc = self.get_total_force_constants(temperature=t, volume=volume)
free_energy, entropy, cv = obtain_phonopy_thermal_properties(self.structure,
temperature=t,
mesh=self._mesh,
force_constants=fc, NAC=self._nac)
print (' {:.1f} {:14.8f} {:14.8f} {:14.8f}'.format(t, free_energy, entropy, cv))
free_energy_list.append(free_energy)
entropy_list.append(entropy)
cv_list.append(cv)
return free_energy_list, entropy_list, cv_list
def plot_density_of_states(self, volume_index=0, temperature=300):
from dynaphopy.interface.phonopy_link import obtain_phonopy_dos
import matplotlib.pyplot as plt
fc = self.get_total_force_constants(temperature=temperature, volume=self.volumes[volume_index])
dos = obtain_phonopy_dos(self.structure, mesh=self._mesh, force_constants=fc)
plt.plot(dos[0], dos[1])
plt.show()
def plot_thermal_properties(self, volume_index=0):
import matplotlib.pyplot as plt
free_energy, entropy, cv = self.get_thermal_properties(volume=self.volumes[volume_index])
temperature = self.get_temperature_range()
plt.plot(temperature, free_energy, label='free energy')
plt.plot(temperature, entropy, label='entropy')
plt.plot(temperature, cv, label='Cv')
plt.grid()
plt.legend()
plt.show()
def qpoint_to_index(self, qpoint, silent=False):
for qindex, qp in enumerate(self.get_data_temperature(0)):
if np.allclose(qp['reduced_wave_vector'], qpoint):
return qindex
if not silent:
print ('{} is not a commensurate point'.format(qpoint))
exit()
return None
def plot_fitted_shifts_temperature(self, qpoint=(0, 0, 0), branch=None, interpolated=False):
import matplotlib.pyplot as plt
qindex = self.qpoint_to_index(qpoint)
shift_matrix = self.get_shift_matrix_temperature()
chk_list = np.arange(self.temperatures[0], self.temperatures[-1], 10)
if interpolated:
chk_shift_matrix = np.array([self.get_interpolated_shifts_temperature(t) for t in chk_list]).T
else:
chk_shift_matrix = np.array([self.get_fitted_shifts_temperature(t) for t in chk_list]).T
plt.suptitle('Frequency shift at wave vector={} (relative to harmonic)'.format(qpoint))
plt.xlabel('Temperature [K]')
plt.ylabel('Frequency shift [THz]')
if branch is None:
plt.plot(chk_list, chk_shift_matrix[:, qindex].T, '-')
plt.plot(self.temperatures, shift_matrix[:, qindex].T, 'o')
else:
plt.title('Branch {}'.format(branch))
plt.plot(chk_list, chk_shift_matrix[branch, qindex].T, '-')
plt.plot(self.temperatures, shift_matrix[branch, qindex].T, 'o')
plt.show()
def plot_fitted_shifts_volumes(self, qpoint=(0, 0, 0), branch=None, interpolate=False):
import matplotlib.pyplot as plt
qindex = self.qpoint_to_index(qpoint)
shift_matrix = self.get_shift_matrix_volume()
chk_list = np.arange(self.volumes[0], self.volumes[-1], 0.1)
if interpolate:
chk_shift_matrix = np.array([self.get_interpolated_shifts_volume(v) for v in chk_list]).T
else:
chk_shift_matrix = np.array([self.get_fitted_shifts_volume(v) for v in chk_list]).T
plt.suptitle('Shifts at wave vector={} (relative to {} K)'.format(qpoint, self.ref_temperature))
plt.xlabel('Volumes [Ang.^3]')
plt.ylabel('Frequency [THz]')
if branch is None:
plt.plot(chk_list, chk_shift_matrix[:, qindex].T, '-')
plt.plot(self.volumes, shift_matrix[:, qindex].T, 'o')
else:
plt.title('Branch {}'.format(branch))
plt.plot(chk_list, chk_shift_matrix[branch, qindex].T, '-')
plt.plot(self.volumes, shift_matrix[branch, qindex].T, 'o')
plt.show()
def plot_linewidth_volume(self, qpoint=(0, 0, 0)):
import matplotlib.pyplot as plt
qindex = self.qpoint_to_index(qpoint)
if qindex is not None:
linewidth = []
for v in range(len(self.volumes)):
data = self.get_data_volume(v)
linewidth.append(data[qindex]['linewidths'])
plt.title('Linewidths at wave vector={}'.format(qpoint))
plt.xlabel('Volume [Angs.^3]')
plt.ylabel('Linewidth [THz]')
plt.plot(self.volumes, linewidth, linestyle='--', marker='o')
plt.show()
else:
print ('qpoint not found!')
return None
def plot_linewidth_temperature(self, qpoint=(0, 0, 0)):
import matplotlib.pyplot as plt
qindex = self.qpoint_to_index(qpoint)
if qindex is not None:
linewidth = []
for t in range(len(self.temperatures)):
data = self.get_data_temperature(t)
linewidth.append(data[qindex]['linewidths'])
plt.title('Linewidths at wave vector={}'.format(qpoint))
plt.xlabel('Temperature [K]')
plt.ylabel('Linewidth [THz]')
plt.plot(self.temperatures, linewidth, linestyle='--', marker='o')
plt.show()
else:
print ('qpoint not found!')
return None
def plot_shift_volume(self, qpoint=(0, 0, 0), branch=None):
import matplotlib.pyplot as plt
qindex = self.qpoint_to_index(qpoint)
shift = []
for v in range(len(self.volumes)):
data = self.get_data_volume(v)
shift.append(data[qindex]['frequency_shifts'])
plt.suptitle('Frequency shift at wave vector={} (relative to {} K)'.format(qpoint, self.ref_temperature))
plt.xlabel('Volume [Angs.^3]')
plt.ylabel('Frequency shift [THz]')
if branch is None:
plt.plot(self.volumes, shift)
else:
plt.title('Branch {}'.format(branch))
plt.plot(self.volumes, np.array(shift).T[branch].T)
plt.show()
def get_shift_temperature(self, qpoint=(0, 0, 0)):
qindex = self.qpoint_to_index(qpoint)
shift = []
for t in range(len(self.temperatures)):
data = self.get_data_temperature(t)
shift.append(data[qindex]['frequency_shifts'])
return np.array(shift).T
def plot_shift_temperature(self, qpoint=(0, 0, 0), branch=None):
import matplotlib.pyplot as plt
qindex = self.qpoint_to_index(qpoint)
shift = []
for t in range(len(self.temperatures)):
data = self.get_data_temperature(t)
shift.append(data[qindex]['frequency_shifts'])
plt.suptitle('Frequency shift at wave vector={} (relative to harmonic)'.format(qpoint))
plt.xlabel('Temperature [K]')
plt.ylabel('Frequency shift [THz]')
if branch is None:
plt.plot(self.temperatures, shift)
else:
plt.title('Branch {}'.format(branch))
plt.plot(self.temperatures, np.array(shift).T[branch].T)
plt.show()
def get_band_structure(self, temperature, volume, band_ranges=None):
from dynaphopy.interface.phonopy_link import obtain_phonon_dispersion_bands
if band_ranges is None:
band_ranges = self.structure.get_path_using_seek_path()['ranges']
fc = self.get_total_force_constants(temperature=temperature, volume=volume)
return obtain_phonon_dispersion_bands(self.structure, bands_ranges=band_ranges,
force_constants=fc, NAC=self._nac,
band_resolution=30, band_connection=False)
def get_dos(self, temperature, volume):
from dynaphopy.interface.phonopy_link import obtain_phonopy_dos
fc = self.get_total_force_constants(temperature=temperature, volume=volume)
return obtain_phonopy_dos(self.structure, mesh=self._mesh, force_constants=fc, NAC=self._nac)
def plot_dos(self, temperature, volume):
import matplotlib.pyplot as plt
dos = self.get_dos(temperature, volume)
plt.plot(dos[0], dos[1])
plt.show()
class QuasiparticlesQHA():
def __init__(self, args, load_data=False, verbose=False, tmin=None):
input_parameters = reading.read_parameters_from_input_file(args.input_file)
if 'structure_file_name_outcar' in input_parameters:
structure = reading.read_from_file_structure_outcar(input_parameters['structure_file_name_outcar'])
else:
structure = reading.read_from_file_structure_poscar(input_parameters['structure_file_name_poscar'])
structure.get_data_from_dict(input_parameters)
if 'supercell_phonon' in input_parameters:
supercell_phonon = input_parameters['supercell_phonon']
else:
supercell_phonon = np.identity(3)
structure.set_force_constants(get_force_constants_from_file(file_name=input_parameters['force_constants_file_name'],
fc_supercell=supercell_phonon))
if '_mesh_phonopy' in input_parameters:
mesh = input_parameters['_mesh_phonopy']
else:
mesh = [20, 20, 20] # Default
print ('mesh set to: {}'.format(mesh))
if 'bands' in input_parameters is None:
self._bands = structure.get_path_using_seek_path()
else:
self._bands = input_parameters['_band_ranges']
volumes, energies = read_v_e(args.ev, factor=1.0, volume_factor=1.0, pressure=args.pressure)
self._fc_fit = ForceConstantsFitting(structure,
files_temperature=args.ct_data,
files_volume=args.cv_data,
temperatures=args.temperatures,
volumes=volumes,
mesh=mesh,
ref_temperature=args.ref_temperature,
fitting_order=args.order,
tmin=tmin,
use_NAC=True)
if not load_data:
temperatures = self._fc_fit.get_temperature_range()
free_energy = []
entropy = []
cv = []
for v, e in zip(volumes, energies):
print ('Volume: {} Ang. Energy(U): {} eV'.format(v, e))
tp_data = self._fc_fit.get_thermal_properties(volume=v)
free_energy.append(tp_data[0])
entropy.append(tp_data[1])
cv.append(tp_data[2])
free_energy = np.array(free_energy).T
entropy = np.array(entropy).T
cv = np.array(cv).T
np.save('free_energy.npy', free_energy)
np.save('temperatures.npy', temperatures)
np.save('cv.npy', cv)
np.save('entropy.npy', entropy)
else:
free_energy = np.load('free_energy.npy')
temperatures = np.load('temperatures.npy')
cv = np.load('cv.npy')
entropy = np.load('entropy.npy')
self.phonopy_qha = PhonopyQHA(volumes,
energies,
eos="vinet", # options: 'vinet', 'murnaghan' or 'birch_murnaghan'
temperatures=temperatures,
free_energy=free_energy,
cv=cv,
entropy=entropy,
t_max=self.fc_fit.get_temperature_range()[-1],
verbose=False)
# Write data files to disk
self.phonopy_qha.write_bulk_modulus_temperature()
self.phonopy_qha.write_gibbs_temperature()
self.phonopy_qha.write_heat_capacity_P_numerical()
self.phonopy_qha.write_gruneisen_temperature()
self.phonopy_qha.write_thermal_expansion()
self.phonopy_qha.write_helmholtz_volume()
self.phonopy_qha.write_volume_expansion()
self.phonopy_qha.write_volume_temperature()
if verbose:
self.phonopy_qha.plot_qha().show()
# Designed for test only
def volume_shift(self, volume_range=np.arange(-2.0, 2.0, 0.1)):
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 1)
volumes = self.phonopy_qha._qha._volumes
energies = self.phonopy_qha._qha._electronic_energies
free_energy = np.load('free_energy.npy')
temperatures = np.load('temperatures.npy')
cv = np.load('cv.npy')
entropy = np.load('entropy.npy')
for i in volume_range:
volumesi = np.array(volumes) + i
print volumesi
phonopy_qha = PhonopyQHA(volumesi,
energies,
eos="vinet", # options: 'vinet', 'murnaghan' or 'birch_murnaghan'
temperatures=temperatures,
free_energy=free_energy,
cv=cv,
entropy=entropy,
t_max=self.fc_fit.get_temperature_range()[-1],
verbose=False)
cp = phonopy_qha.get_heat_capacity_P_numerical()
import matplotlib.pyplot as plt
import matplotlib.colors as colors
cNorm = colors.Normalize(vmin=volume_range[0], vmax=volume_range[-1])
scalarMap = plt.cm.ScalarMappable(norm=cNorm, cmap=plt.cm.get_cmap('plasma'))
ax.plot(phonopy_qha._qha._temperatures[:-3], cp, label='{}'.format(i), color=scalarMap.to_rgba(i))
import matplotlib as mpl
ax2 = fig.add_axes([0.93, 0.1, 0.02, 0.8])
mpl.colorbar.ColorbarBase(ax2, cmap=plt.cm.get_cmap('plasma'), norm=cNorm,
spacing='proportional', ticks=volume_range,
boundaries=None, format='%1i')
plt.show()
def plot_dos_gradient(self):
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.colorbar as colorbar
volumes = self.phonopy_qha.get_volume_temperature()
temperatures = self.fc_fit.get_temperature_range()
fig, ax = plt.subplots(1,1)
for t, v in zip(temperatures[::40], volumes[::20]):
print ('temperature: {} K'.format(t))
dos = self.fc_fit.get_dos(t, v)
cNorm = colors.Normalize(vmin=temperatures[0], vmax=temperatures[-1])
scalarMap = plt.cm.ScalarMappable(norm=cNorm, cmap=plt.cm.get_cmap('plasma'))
ax.plot(dos[0], dos[1], color=scalarMap.to_rgba(t))
plt.suptitle('Phonon density of states')
plt.xlabel('Frequency [THz]')
ax2 = fig.add_axes([0.93, 0.1, 0.02, 0.8])
colorbar.ColorbarBase(ax2, cmap=plt.cm.get_cmap('plasma'), norm=cNorm,
spacing='proportional', ticks=temperatures[::40],
boundaries=None, format='%1i')
plt.show()
def plot_band_structure_gradient(self, tmin=300, tmax=1600, tstep=100):
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.colorbar as colorbar
def replace_list(text_string):
substitutions = {'GAMMA': u'$\Gamma$',
}
for item in substitutions.iteritems():
text_string = text_string.replace(item[0], item[1])
return text_string
volumes = self.phonopy_qha.get_volume_temperature()
cNorm = colors.Normalize(vmin=tmin, vmax=tmax)
fig, ax = plt.subplots(1,1)
for t in np.arange(tmin, tmax, tstep):
print ('temperature: {} K'.format(t))
v = self.get_volume_at_temperature(t)
scalarMap = plt.cm.ScalarMappable(norm=cNorm, cmap=plt.cm.get_cmap('plasma'))
band_data = self.fc_fit.get_band_structure(t, v, band_ranges=self._bands['ranges'])
for i, freq in enumerate(band_data[1]):
ax.plot(band_data[1][i], band_data[2][i], color=scalarMap.to_rgba(t))
# plt.axes().get_xaxis().set_visible(False)
#plt.axes().get_xaxis().set_ticks([])
plt.ylabel('Frequency [THz]')
plt.xlabel('Wave vector')
plt.xlim([0, band_data[1][-1][-1]])
plt.axhline(y=0, color='k', ls='dashed')
plt.suptitle('Phonon dispersion')
if 'labels' in self._bands:
plt.rcParams.update({'mathtext.default': 'regular'})
labels = self._bands['labels']
labels_e = []
x_labels = []
for i in range(len(band_data[1])):
if labels[i][0] == labels[i - 1][1]:
labels_e.append(replace_list(labels[i][0]))
else:
labels_e.append(
replace_list(labels[i - 1][1]) + '/' + replace_list(labels[i][0]))
x_labels.append(band_data[1][i][0])
x_labels.append(band_data[1][-1][-1])
labels_e.append(replace_list(labels[-1][1]))
labels_e[0] = replace_list(labels[0][0])
plt.xticks(x_labels, labels_e, rotation='horizontal')
ax2 = fig.add_axes([0.93, 0.1, 0.02, 0.8])
colorbar.ColorbarBase(ax2, cmap=plt.cm.get_cmap('plasma'), norm=cNorm,
spacing='proportional', ticks=np.arange(tmin, tmax, tstep),
boundaries=None, format='%1i')
plt.show()
def get_volume_at_temperature(self, temperature):
temperatures = self.get_qha_temperatures()
volumes = self.phonopy_qha.get_volume_temperature()
volume = np.interp(temperature, temperatures, volumes)
return volume
def plot_band_structure_constant_pressure(self, temperature=300, external_data=None):
import matplotlib.pyplot as plt
def replace_list(text_string):
substitutions = {'GAMMA': u'$\Gamma$',
}
for item in substitutions.iteritems():
text_string = text_string.replace(item[0], item[1])
return text_string
volume = self.get_volume_at_temperature(temperature)
fig, ax = plt.subplots(1,1)
band_data = self.fc_fit.get_band_structure(temperature, volume, band_ranges=self._bands['ranges'])
for i, freq in enumerate(band_data[1]):
ax.plot(band_data[1][i], band_data[2][i], color='r')
#plt.axes().get_xaxis().set_ticks([])
plt.ylabel('Frequency [THz]')
plt.xlabel('Wave vector')
plt.xlim([0, band_data[1][-1][-1]])
plt.axhline(y=0, color='k', ls='dashed')
plt.suptitle('Phonon dispersion')
if 'labels' in self._bands:
plt.rcParams.update({'mathtext.default': 'regular'})
labels = self._bands['labels']
labels_e = []
x_labels = []
for i in range(len(band_data[1])):
if labels[i][0] == labels[i - 1][1]:
labels_e.append(replace_list(labels[i][0]))
else:
labels_e.append(
replace_list(labels[i - 1][1]) + '/' + replace_list(labels[i][0]))
x_labels.append(band_data[1][i][0])
x_labels.append(band_data[1][-1][-1])
labels_e.append(replace_list(labels[-1][1]))
labels_e[0] = replace_list(labels[0][0])
plt.xticks(x_labels, labels_e, rotation='horizontal')
ax.plot(external_data[:, 0], external_data[:, 1], 'o', color='b')
plt.show()
def get_qha_temperatures(self):
max_t_index = self.phonopy_qha._qha._get_max_t_index(self.phonopy_qha._qha._temperatures)
temperatures = self.phonopy_qha._qha._temperatures[:max_t_index]
return temperatures
def get_FC_constant_pressure(self):
temperatures = self.get_qha_temperatures()
volumes = self.phonopy_qha.get_volume_temperature()
for t, v in zip(temperatures[::20], volumes[::20]):
fc = self.fc_fit.get_total_force_constants(temperature=t, volume=v)
write_FORCE_CONSTANTS(fc.get_array(), filename='FC_{}'.format(t))
def get_total_shift_constant_pressure(self, qpoint=(0, 0, 0)):
qindex = self.fc_fit.qpoint_to_index(qpoint)
volumes = self.phonopy_qha.get_volume_temperature()
temperatures = self.get_qha_temperatures()
h_frequencies, ev = self.fc_fit.get_harmonic_frequencies_and_eigenvectors()
chk_shift_matrix = []
for v, t in zip(volumes, temperatures):
total_shifts = self.fc_fit.get_total_shifts(volume=v, temperature=t)
chk_shift_matrix.append(total_shifts)
chk_shift_matrix = np.array(chk_shift_matrix).T
return chk_shift_matrix[:, qindex]
def plot_total_shift_constant_pressure(self, qpoint=(0, 0, 0), branch=None):
import matplotlib.pyplot as plt
temperatures = self.get_qha_temperatures()
chk_shift_matrix = self.get_total_shift_constant_pressure(qpoint=qpoint)
plt.suptitle('Total frequency shift at wave vector={} (relative to {} K)'.format(qpoint, self.fc_fit.ref_temperature))
plt.xlabel('Temperature [K]')
plt.ylabel('Frequency shift [THz]')
if branch is None:
plt.plot(temperatures, chk_shift_matrix.T, '-')
else:
plt.title('Branch {}'.format(branch))
plt.plot(temperatures, chk_shift_matrix[branch].T, '-')
plt.show()
@property
def fc_fit(self):
return self._fc_fit
# FROM HERE DIRTY DEVELOPMENT
qp = QuasiparticlesQHA(args, load_data=False, verbose=True, tmin=200)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5])
qp.fc_fit.plot_fitted_shifts_volumes(qpoint=[0.5, 0.5, 0.5])
experiment_data = np.loadtxt('FC/exp_bands')
qp.plot_band_structure_constant_pressure(external_data=experiment_data)
#qp.plot_band_structure_gradient()
qp.fc_fit.plot_shift_volume(qpoint=[0.5, 0.5, 0.5], branch=0)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5], branch=0)
qp.plot_total_shift_constant_pressure(qpoint=[0.5, 0.5, 0.5], branch=0)
qp.plot_total_shift_constant_pressure(qpoint=[0.5, 0.5, 0.5], branch=2)
np.savetxt('shift_r', qp.get_total_shift_constant_pressure(qpoint=[0.5, 0.5, 0.5]).T)
np.savetxt('tshift_r', qp.fc_fit.get_shift_temperature(qpoint=[0.5, 0.5, 0.5]).T)
print qp.fc_fit.temperatures
print qp.get_qha_temperatures()
exit()
qp.plot_total_shift_constant_pressure(qpoint=[0.0, 0.0, 0.0], branch=4)
qp.plot_total_shift_constant_pressure(qpoint=[0.0, 0.0, 0.0])
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.0, 0.0, 0.0])
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5], branch=0)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5], branch=1)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5], branch=2)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5], branch=3)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5], branch=4)
qp.fc_fit.plot_fitted_shifts_temperature(qpoint=[0.5, 0.5, 0.5])
qp.fc_fit.plot_fitted_shifts_volumes(qpoint=[0.5, 0.5, 0.5])
qp.fc_fit.plot_linewidth_temperature(qpoint=[0.5, 0.5, 0.5])
qp.fc_fit.plot_linewidth_volume(qpoint=[0.5, 0.5, 0.5])
qp.plot_dos_gradient()
qp.fc_fit.plot_linewidth_volume([0, 0, 0])
