diff --git a/src/perturbopy/postproc/__init__.py b/src/perturbopy/postproc/__init__.py index 31204d7..fa31f01 100644 --- a/src/perturbopy/postproc/__init__.py +++ b/src/perturbopy/postproc/__init__.py @@ -3,6 +3,7 @@ """ from .calc_modes.calc_mode import CalcMode from .calc_modes.bands import Bands +from .calc_modes.spectral_cumulant import SpectralCumulant from .calc_modes.spins import Spins from .calc_modes.phdisp import Phdisp from .calc_modes.ephmat import Ephmat diff --git a/src/perturbopy/postproc/calc_modes/calc_mode.py b/src/perturbopy/postproc/calc_modes/calc_mode.py index 63df9bd..f4ad917 100644 --- a/src/perturbopy/postproc/calc_modes/calc_mode.py +++ b/src/perturbopy/postproc/calc_modes/calc_mode.py @@ -45,8 +45,9 @@ def __init__(self, pert_dict): """ - # Extract calculation mode name from pert_dict + # Extract calculation mode name and prefix from pert_dict self.calc_mode = pert_dict['input parameters']['after conversion'].pop('calc_mode') + self.prefix = pert_dict['input parameters']['after conversion'].pop('prefix') # Extract basic data from pert_dict self.alat = pert_dict['basic data']['alat'] diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 476ef66..14bb68f 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -34,7 +34,7 @@ def __init__(self, cdyna_file, tet_file, pert_dict): Dictionary containing the inputs and outputs from the dynamics-run calculation. """ - + self.timings = TimingGroup("dynamics-run") super().__init__(pert_dict) @@ -175,7 +175,7 @@ def extract_steady_drift_vel(self, dyna_pp_yaml_path): raise FileNotFoundError(f'File {dyna_pp_yaml_path} not found') dyna_pp_dict = open_yaml(dyna_pp_yaml_path) - + vels = dyna_pp_dict['dynamics-pp']['velocity'] concs = dyna_pp_dict['dynamics-pp']['concentration'] @@ -185,9 +185,9 @@ def extract_steady_drift_vel(self, dyna_pp_yaml_path): steady_conc = [] for irun, dynamics_run in self._data.items(): - + step_number += dynamics_run.num_steps - + if np.allclose(dynamics_run.efield, np.array([0.0, 0.0, 0.0])): steady_drift_vel.append(None) steady_conc.append(concs[step_number]) diff --git a/src/perturbopy/postproc/calc_modes/spectral_cumulant.py b/src/perturbopy/postproc/calc_modes/spectral_cumulant.py new file mode 100644 index 0000000..1d3e7ad --- /dev/null +++ b/src/perturbopy/postproc/calc_modes/spectral_cumulant.py @@ -0,0 +1,132 @@ +import numpy as np +import h5py as h5 +import matplotlib.pyplot as plt +from perturbopy.postproc.calc_modes.calc_mode import CalcMode +from perturbopy.io_utils.io import open_yaml, open_hdf5, close_hdf5 +import os + + +class SpectralCumulant(CalcMode): + """ + Class representation of a Perturbo cumulant calculation + + Parameters + ---------- + prefix : str + Prefix for the HDF5 file containing spectral function data. + + Attributes + ---------- + temp_array : numpy.ndarray + Array of temperatures corresponding to the spectral data. + freq_array : numpy.ndarray + Array of energy value in electron volts (eV). + freq_step : float + Energy step size for the energy grid in eV. + Akw : numpy.ndarray + Spectral function data, indexed by k-point, band, temperature, and ω. + """ + def __init__(self, spectral_file, pert_dict): + """ + Constructor method + + Parameters + ---------- + spectral_file : dict + Dictionary for the prefix_spectral_cumulant.h5 + pert_dict : dict + Dictionary containing the inputs from the spectral-cum calculation. + + """ + super().__init__(pert_dict) + if self.calc_mode != 'spectral-cum': + raise ValueError('Calculation mode for a SpectralCumulantCalcMode object should be "spectral-cum"') + + Akw = spectral_file['spectral_functions'] + w_lower = np.asarray(spectral_file['w_lower_index']) + w_upper = np.asarray(spectral_file['w_upper_index']) + freq_step = np.asarray(spectral_file['wfreq_step_eV']) + self.temp_array = np.asanyarray(spectral_file['temperatures']) + self.freq_array = np.arange(w_lower, w_upper + 1) * freq_step + self.freq_step = freq_step + Akw_np = [] + for key in Akw.keys(): + Akw_np.append(np.asarray(Akw[key])) + close_hdf5(spectral_file) + + self.Akw = np.asarray(Akw_np) + + @classmethod + def from_hdf5_yaml(cls, spectral_path, yaml_path='pert_output.yml'): + """ + Class method to create a SpectralCumulantCalcMode object from the HDF5 file and YAML file + generated by a Perturbo calculation + + Parameters + ---------- + spectral_path : str + Path to the HDF5 file generated by a spectral-cum calculation + yaml_path : str, optional + Path to the YAML file generated by a spectral-cum calculation + + Returns + ------- + + """ + + if not os.path.isfile(spectral_path): + raise FileNotFoundError(f'File {spectral_path} not found') + if not os.path.isfile(yaml_path): + raise FileNotFoundError(f'File {yaml_path} not found') + + spectral_file = open_hdf5(spectral_path) + yaml_dict = open_yaml(yaml_path) + + return cls(spectral_file, yaml_dict) + + def plot_Aw(self, ax, ik=0, it=0, ib=0): + """ + Plots the spectral function A(ω) for a given k-point, temperature, and band. + + Parameters + ---------- + ax : matplotlib.axes.Axes + Axis on which to plot the spectral function. + ik : int, optional + Index of the k-point in the grid. Default is 0. + it : int, optional + Index of the temperature in the temperature array. Default is 0. + ib : int, optional + Index of the band. Default is 0. + + Returns + ------- + matplotlib.axes.Axes + The axis object with the plotted spectral function. + + Notes + ----- + The spectral function is normalized before plotting. The plot includes + labels, legends, and adjusted aesthetics for better visualization. + """ + if it > len(self.temp_array): + raise ValueError('Temperature index is out of range') + if ib > self.Akw[0].shape[0]: + raise ValueError('Band index is out of range') + if ik > len(self.Akw): + raise ValueError('k-point index is out of range') + A0w = self.Akw[ik][ib, it, :] + freq_step = self.freq_step + freq_array = self.freq_array + # normalize + A0w = A0w / (np.sum(A0w) * freq_step) + # plot + ax.plot(freq_array, A0w, lw=2, label=f'T={ int(self.temp_array[it])} K') + ax.legend(fontsize=18) + plt.ylim([0, np.max(A0w) * 1.1]) + plt.xlabel(r'$\omega-\epsilon_{nk}$ (eV)', fontsize=20) + plt.ylabel(r'A($\omega$) (eV$^{-1}$)', fontsize=20) + plt.xticks(fontsize=20) + plt.yticks(fontsize=20) + plt.tight_layout() + return ax diff --git a/tests/refs/sto_spectral-cum.yml b/tests/refs/sto_spectral-cum.yml new file mode 100644 index 0000000..1a7191d --- /dev/null +++ b/tests/refs/sto_spectral-cum.yml @@ -0,0 +1,405 @@ +--- +program: perturbo + +start date and time: + date: 6 December 2024 + time: "14:23:16" + +parallelization: + type: MPI+OpenMP + number of MPI tasks: 1 + number of OpenMP threads: 12 + number of nodes: 1 + +input parameters: + before conversion (before reading from epr file): + band_max: 3 + band_min: 1 + boltz_acc_thr: 0.0000000000000000E+00 + boltz_de: 1.0000000000000000E+00 + boltz_efield: [ 0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, ] + boltz_emax: 9.9990000000000000E+03 + boltz_emin: -9.9990000000000000E+03 + boltz_init_dist: + boltz_init_e0: -9.9990000000000000E+03 + boltz_init_smear: 2.0000000000000000E+01 + boltz_kdim: [ 1, 1, 1, ] + boltz_norm_dist: False + boltz_nstep: 0 + boltz_nstep_min: 1 + boltz_qdim: [ 1, 1, 1, ] + calc_mode: spectral-cum + cauchy_scale: 5.0000000000000003E-02 + scat_impl: std + cum_de: 4.0000000000000002E-01 + cum_inner_emax: 2.0000000000000001E-01 + cum_inner_emin: -2.0000000000000001E-01 + cum_outer_emax: 5.9999999999999998E-01 + cum_outer_emin: -5.9999999999999998E-01 + cum_outer_np: 4 + debug: False + delta_smear: 1.0000000000000000E+01 + find_efermi: False + fklist: band.kpt + fqlist: + ftemper: sto.temper + full_ite: True + hole: False + lmagsym: False + load_scatter_eph: False + nsamples: 100000 + output_nstep: 1 + phfreq_cutoff: 1.0000000000000000E+00 + polar_split: + prefix: sto + sampling: uniform + solver: rk4 + spectral_emax: 5.9999999999999998E-01 + spectral_emin: -5.9999999999999998E-01 + spectral_np: 1 + time_step: 1.0000000000000000E+00 + tmp_dir: ./ + trans_thr: 2.0000000000000000E-03 + scat_impl: std + use_mem: True + yaml_fname: pert_output.yml + after conversion: + band_max: 3 + band_min: 1 + boltz_acc_thr: 0.0000000000000000E+00 + boltz_de: 7.3498617648950546E-05 + boltz_efield: [ 0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, ] + boltz_emax: 7.3491267787185654E+02 + boltz_emin: -7.3491267787185654E+02 + boltz_init_dist: + boltz_init_e0: -7.3491267787185654E+02 + boltz_init_smear: 1.4699723529790110E-03 + boltz_kdim: [ 1, 1, 1, ] + boltz_norm_dist: False + boltz_nstep: 0 + boltz_nstep_min: 1 + boltz_qdim: [ 1, 1, 1, ] + calc_mode: spectral-cum + cauchy_scale: 5.0000000000000003E-02 + scat_impl: std + cum_de: 2.9399447059580220E-05 + cum_inner_emax: 1.4699723529790111E-02 + cum_inner_emin: -1.4699723529790111E-02 + cum_outer_emax: 4.4099170589370330E-02 + cum_outer_emin: -4.4099170589370330E-02 + cum_outer_np: 4 + debug: False + delta_smear: 7.3498617648950549E-04 + find_efermi: False + fklist: band.kpt + fqlist: + ftemper: sto.temper + full_ite: True + hole: False + lmagsym: False + load_scatter_eph: False + nsamples: 100000 + output_nstep: 1 + phfreq_cutoff: 7.3498617648950546E-05 + polar_split: + prefix: sto + sampling: uniform + solver: rk4 + spectral_emax: 4.4099170589370330E-02 + spectral_emin: -4.4099170589370330E-02 + spectral_np: 1 + time_step: 2.0670686467503085E+01 + tmp_dir: ./ + trans_thr: 2.0000000000000000E-03 + scat_impl: std + use_mem: True + yaml_fname: sto_spectral-cum.yml + + +basic data: + alat: 7.3699300000 + alat units: bohr + # a vector is a row + lattice vectors: + - [ 1.00000, 0.00000, 0.00000, ] + - [ 0.00000, 1.00000, 0.00000, ] + - [ 0.00000, 0.00000, 1.00000, ] + lattice vectors units: alat + reciprocal lattice vectors: + # a vector is a row + - [ 1.00000, 0.00000, 0.00000, ] + - [ 0.00000, 1.00000, 0.00000, ] + - [ 0.00000, 0.00000, 1.00000, ] + reciprocal lattice vectors units: 2pi/alat + number of atoms in unit cell: 5 + atomic positions: + # a vector is a row + - [ 0.00000, 0.00000, 0.00000, ] + - [ 0.50000, 0.50000, 0.50000, ] + - [ 0.00000, 0.50000, 0.50000, ] + - [ 0.50000, 0.00000, 0.50000, ] + - [ 0.50000, 0.50000, 0.00000, ] + atomic positions units: alat + volume: 400.3041465593 + volume units: bohr^3 + number of symmetry operations: 48 + symop: + 1: + - [ 1, 0, 0, ] + - [ 0, 1, 0, ] + - [ 0, 0, 1, ] + 2: + - [ -1, 0, 0, ] + - [ 0, -1, 0, ] + - [ 0, 0, 1, ] + 3: + - [ -1, 0, 0, ] + - [ 0, 1, 0, ] + - [ 0, 0, -1, ] + 4: + - [ 1, 0, 0, ] + - [ 0, -1, 0, ] + - [ 0, 0, -1, ] + 5: + - [ 0, 1, 0, ] + - [ 1, 0, 0, ] + - [ 0, 0, -1, ] + 6: + - [ 0, -1, 0, ] + - [ -1, 0, 0, ] + - [ 0, 0, -1, ] + 7: + - [ 0, 1, 0, ] + - [ -1, 0, 0, ] + - [ 0, 0, 1, ] + 8: + - [ 0, -1, 0, ] + - [ 1, 0, 0, ] + - [ 0, 0, 1, ] + 9: + - [ 0, 0, 1, ] + - [ 0, -1, 0, ] + - [ 1, 0, 0, ] + 10: + - [ 0, 0, -1, ] + - [ 0, -1, 0, ] + - [ -1, 0, 0, ] + 11: + - [ 0, 0, 1, ] + - [ 0, 1, 0, ] + - [ -1, 0, 0, ] + 12: + - [ 0, 0, -1, ] + - [ 0, 1, 0, ] + - [ 1, 0, 0, ] + 13: + - [ -1, 0, 0, ] + - [ 0, 0, 1, ] + - [ 0, 1, 0, ] + 14: + - [ -1, 0, 0, ] + - [ 0, 0, -1, ] + - [ 0, -1, 0, ] + 15: + - [ 1, 0, 0, ] + - [ 0, 0, 1, ] + - [ 0, -1, 0, ] + 16: + - [ 1, 0, 0, ] + - [ 0, 0, -1, ] + - [ 0, 1, 0, ] + 17: + - [ 0, 1, 0, ] + - [ 0, 0, 1, ] + - [ 1, 0, 0, ] + 18: + - [ 0, -1, 0, ] + - [ 0, 0, 1, ] + - [ -1, 0, 0, ] + 19: + - [ 0, 1, 0, ] + - [ 0, 0, -1, ] + - [ -1, 0, 0, ] + 20: + - [ 0, -1, 0, ] + - [ 0, 0, -1, ] + - [ 1, 0, 0, ] + 21: + - [ 0, 0, 1, ] + - [ 1, 0, 0, ] + - [ 0, 1, 0, ] + 22: + - [ 0, 0, 1, ] + - [ -1, 0, 0, ] + - [ 0, -1, 0, ] + 23: + - [ 0, 0, -1, ] + - [ -1, 0, 0, ] + - [ 0, 1, 0, ] + 24: + - [ 0, 0, -1, ] + - [ 1, 0, 0, ] + - [ 0, -1, 0, ] + 25: + - [ -1, 0, 0, ] + - [ 0, -1, 0, ] + - [ 0, 0, -1, ] + 26: + - [ 1, 0, 0, ] + - [ 0, 1, 0, ] + - [ 0, 0, -1, ] + 27: + - [ 1, 0, 0, ] + - [ 0, -1, 0, ] + - [ 0, 0, 1, ] + 28: + - [ -1, 0, 0, ] + - [ 0, 1, 0, ] + - [ 0, 0, 1, ] + 29: + - [ 0, -1, 0, ] + - [ -1, 0, 0, ] + - [ 0, 0, 1, ] + 30: + - [ 0, 1, 0, ] + - [ 1, 0, 0, ] + - [ 0, 0, 1, ] + 31: + - [ 0, -1, 0, ] + - [ 1, 0, 0, ] + - [ 0, 0, -1, ] + 32: + - [ 0, 1, 0, ] + - [ -1, 0, 0, ] + - [ 0, 0, -1, ] + 33: + - [ 0, 0, -1, ] + - [ 0, 1, 0, ] + - [ -1, 0, 0, ] + 34: + - [ 0, 0, 1, ] + - [ 0, 1, 0, ] + - [ 1, 0, 0, ] + 35: + - [ 0, 0, -1, ] + - [ 0, -1, 0, ] + - [ 1, 0, 0, ] + 36: + - [ 0, 0, 1, ] + - [ 0, -1, 0, ] + - [ -1, 0, 0, ] + 37: + - [ 1, 0, 0, ] + - [ 0, 0, -1, ] + - [ 0, -1, 0, ] + 38: + - [ 1, 0, 0, ] + - [ 0, 0, 1, ] + - [ 0, 1, 0, ] + 39: + - [ -1, 0, 0, ] + - [ 0, 0, -1, ] + - [ 0, 1, 0, ] + 40: + - [ -1, 0, 0, ] + - [ 0, 0, 1, ] + - [ 0, -1, 0, ] + 41: + - [ 0, -1, 0, ] + - [ 0, 0, -1, ] + - [ -1, 0, 0, ] + 42: + - [ 0, 1, 0, ] + - [ 0, 0, -1, ] + - [ 1, 0, 0, ] + 43: + - [ 0, -1, 0, ] + - [ 0, 0, 1, ] + - [ 1, 0, 0, ] + 44: + - [ 0, 1, 0, ] + - [ 0, 0, 1, ] + - [ -1, 0, 0, ] + 45: + - [ 0, 0, -1, ] + - [ -1, 0, 0, ] + - [ 0, -1, 0, ] + 46: + - [ 0, 0, -1, ] + - [ 1, 0, 0, ] + - [ 0, 1, 0, ] + 47: + - [ 0, 0, 1, ] + - [ 1, 0, 0, ] + - [ 0, -1, 0, ] + 48: + - [ 0, 0, 1, ] + - [ -1, 0, 0, ] + - [ 0, 1, 0, ] + kc dimensions: + - [ 4, 4, 4, ] + spinor: False + polar_alpha: 1.0000000000 + epsil: + # a vector is a row + - [ 6.25856984, 0.00000000, 0.00000000, ] + - [ 0.00000000, 6.25856984, 0.00000000, ] + - [ 0.00000000, 0.00000000, 6.25856984, ] + qc dimensions: + - [ 2, 2, 2, ] + mass: + - 79860.74458118 + - 43628.10158488 + - 14582.19644550 + - 14582.19644550 + - 14582.19644550 + mass units: atomic + zstar: + 1: + - [ 2.55389, 0.00000, 0.00000, ] + - [ 0.00000, 2.55389, 0.00000, ] + - [ 0.00000, 0.00000, 2.55389, ] + 2: + - [ 7.28719, 0.00000, 0.00000, ] + - [ 0.00000, 7.28719, 0.00000, ] + - [ 0.00000, 0.00000, 7.28719, ] + 3: + - [ -5.77445, 0.00000, 0.00000, ] + - [ 0.00000, -2.03331, 0.00000, ] + - [ 0.00000, 0.00000, -2.03331, ] + 4: + - [ -2.03331, 0.00000, 0.00000, ] + - [ 0.00000, -5.77445, 0.00000, ] + - [ 0.00000, 0.00000, -2.03331, ] + 5: + - [ -2.03331, 0.00000, 0.00000, ] + - [ 0.00000, -2.03331, 0.00000, ] + - [ 0.00000, 0.00000, -5.77445, ] + system_2d: False + number of Wannier functions: 3 + wannier_center: + # a vector is a row + - [ 0.50000, 0.50000, 0.50000, ] + - [ 0.50000, 0.50000, 0.50000, ] + - [ 0.50000, 0.50000, 0.50000, ] + wannier_center units: cartesian + wannier_center_cryst: + # a vector is a row + - [ 0.50000, 0.50000, 0.50000, ] + - [ 0.50000, 0.50000, 0.50000, ] + - [ 0.50000, 0.50000, 0.50000, ] + wannier_center_cryst units: crystal + +memory: + physical: 65737300 kB + max rss: 54812 kB + +timings: + + time units: s + + PERTURBO: + cpu: 3.10515600 + wall: 0.35491920 + calls: 1 +... diff --git a/tests/refs/sto_spectral_cumulant.h5 b/tests/refs/sto_spectral_cumulant.h5 new file mode 100644 index 0000000..850ef71 Binary files /dev/null and b/tests/refs/sto_spectral_cumulant.h5 differ diff --git a/tests/test_spectral_cum.py b/tests/test_spectral_cum.py new file mode 100644 index 0000000..8d1af09 --- /dev/null +++ b/tests/test_spectral_cum.py @@ -0,0 +1,55 @@ +import numpy as np +import pytest +import os + +import perturbopy.postproc as ppy + + +@pytest.fixture() +def sto_spectral_cum(): + """ + Method to generate the spectral_cum object corresponding to spectralhdf5. + + Returns + ------- + phdisp : ppy.PhdispCalcMode + + """ + yml_path = os.path.join("refs", "sto_spectral-cum.yml") + spectral_path = os.path.join("refs", "sto_spectral_cumulant.h5") + return ppy.SpectralCumulant.from_hdf5_yaml(spectral_path, yml_path) + +def test_plot_spectral_cum(sto_spectral_cum, plt, with_plt): + """ + Method to test SpectralCumulant.plot_Aw function + + Method to plot the spectral function A(ω) + + Parameters + ---------- + show_qpoint_labels : bool, optional + If true, the q-point labels stored in the labels attribute will be shown on the plot. Default true. + + """ + + if not with_plt: + pytest.skip("Test requires pytest-plt") + + fig, ax = plt.subplots() + ppy.SpectralCumulant.plot_Aw(sto_spectral_cum, ax) + +def test_spectral_cum(): + """ + Method to test SpectralCumulant.plot_Aw function array shape + + Parameters + ---------- + + """ + + yml_path = os.path.join("refs", "sto_spectral-cum.yml") + spectral_path = os.path.join("refs", "sto_spectral_cumulant.h5") + model = ppy.SpectralCumulant.from_hdf5_yaml(spectral_path, yml_path) + np.testing.assert_equal(model.Akw.shape, (1, 3, 2, 3001)) + np.testing.assert_equal(model.freq_array.shape, (3001,)) +