From 6036d5919e771d546bc169794ec05c2e204c61c1 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Tue, 24 Dec 2024 18:56:06 +0100 Subject: [PATCH 01/21] added link to the pump pulse HDF5 file --- src/perturbopy/yml_files/files_description.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/perturbopy/yml_files/files_description.yml b/src/perturbopy/yml_files/files_description.yml index 365d356..09b64e9 100644 --- a/src/perturbopy/yml_files/files_description.yml +++ b/src/perturbopy/yml_files/files_description.yml @@ -158,7 +158,8 @@ prefix_popu.h5: pump_pulse.h5: type: HDF5 - description: "[OPTIONAL] HDF5 file containing the pump pulse information, used if pump_pulse=.true." + description: "[OPTIONAL] HDF5 file containing the pump pulse information, used if pump_pulse=.true.. The name of the pump pulse file can be specified with the pump_pulse_fname input parameter." + format example: mydoc_dynamics.html#pump_pulse_h5_file prefix.trans_coef: type: text From 08f8294cb14131372336b11a849a627e7eff8a75 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sat, 28 Dec 2024 19:18:01 +0100 Subject: [PATCH 02/21] added a script to plot pump pulse occupations on bands --- .../postproc/utils/spectra_generate_pulse.py | 12 +--- .../postproc/utils/spectra_plots.py | 67 +++++++++++++++++++ 2 files changed, 69 insertions(+), 10 deletions(-) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 8bcf6b0..a0c999f 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -20,14 +20,6 @@ from . import spectra_plots -def gaussian(x, mu, sig, hole_nband, elec_nband): - """ - Gaussian function normalized to unity max occupation - """ - - return np.exp(-0.5 * ((x - mu) / sig)**2) / (hole_nband * elec_nband) - - def sigma_from_fwhm(fwhm): """ Comupte Gaussian sigma from the Full Width at Half Maximum (FWHM) parameter: @@ -308,8 +300,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, for jband in range(hole_nband): # diff. between elec and hole for a given elec branch iband and hole branch jband delta = pump_factor * \ - gaussian(elec_energy_array[ekidx, iband] - hole_energy_array[hkidx, jband], - pump_energy, pump_energy_broadening, hole_nband, elec_nband) + spectra_plots.gaussian(elec_energy_array[ekidx, iband] - hole_energy_array[hkidx, jband], + pump_energy, pump_energy_broadening, hole_nband, elec_nband) # Only for the intersected k points, we add the delta occupation elec_occs_amplitude[ekidx, iband] += delta diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index d707ac3..ec2e952 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -22,9 +22,18 @@ # Plotting parameters from .plot_tools import plotparams + plt.rcParams.update(plotparams) +def gaussian(x, mu, sig, hole_nband, elec_nband): + """ + Gaussian function normalized to unity max occupation + """ + + return np.exp(-0.5 * ((x - mu) / sig)**2) / (hole_nband * elec_nband) + + def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, h_occs, hole_kpoint_array, hole_energy_array, pump_energy): """ @@ -264,3 +273,61 @@ def plot_trans_abs_map(ax, time_grid, energy_grid, trans_abs, cbar.set_label(r'$\Delta A$ (arb. units)', fontsize=fsize) return ax + + +def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, + pump_energy, pump_energy_broadening, + scale=0.2): + """ + Plot the occupations created by a pump pulse on the bands. + """ + + + # Total number of bands + nband = len(bands) + + # Total number of k-points + npoint = len(bands[1]) + + # Occupation amplitude for each bands + occs_amplitude_bands = np.zeros((nband, npoint)) + + elec_band_list = list(range(first_el_band_idx, nband)) + hole_band_list = list(range(1, first_el_band_idx)) + + elec_nband = len(elec_band_list) + hole_nband = len(hole_band_list) + + # First, compute the occupations for every valence-conduction pair + for i_elec_band in elec_band_list: + for i_hole_band in hole_band_list: + elec_band = bands[i_elec_band] + hole_band = bands[i_hole_band] + + occ = gaussian(elec_band - hole_band, + pump_energy, pump_energy_broadening, hole_nband, elec_nband) + + occs_amplitude_bands[i_elec_band - 1, :] += occ + occs_amplitude_bands[i_hole_band - 1, :] += occ + + max_occ = np.max(occs_amplitude_bands.ravel()) + factor = scale / max_occ + occs_amplitude_bands *= factor + + # Plot the occupations on bands with fill_between + # use kpath + for iband in range(1, nband): + + band = bands[iband] + + occ_top = band + occs_amplitude_bands[iband - 1, :] + occ_bottom = band - occs_amplitude_bands[iband - 1, :] + + if iband >= first_el_band_idx: + color = 'red' + else: + color = 'blue' + + ax.fill_between(kpath, occ_top, occ_bottom, color=color) + + return ax From 243487381abeb1b901c1aa78b7c05e48e64eaf49 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sun, 29 Dec 2024 12:26:34 +0100 Subject: [PATCH 03/21] Updated pump pulse HDF5 generation: removed unnecessary datasets --- src/perturbopy/postproc/utils/spectra_generate_pulse.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index a0c999f..c396b44 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -350,19 +350,15 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, h5f['optional_params'] = optional_params - h5f.create_dataset('finite_width', data=finite_width) if finite_width: h5f.create_dataset('time_window', data=pump_time_window) h5f.create_dataset('num_steps', data=num_steps) h5f.create_dataset('time_step', data=pump_time_step) - h5f.create_dataset('pulse_type', data=f'gaussian; FWHM: {pump_fwhm:.3f} fs; ' - f'pump_factor: {pump_factor:.3f}') else: h5f.create_dataset('time_window', data=elec_dyna_time_step) h5f.create_dataset('num_steps', data=1) h5f.create_dataset('time_step', data=elec_dyna_time_step) - h5f.create_dataset('pulse_type', data='step') h5f['time_window'].attrs['units'] = 'fs' h5f['time_step'].attrs['units'] = 'fs' @@ -375,10 +371,6 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, elec_pump_pulse_file.create_dataset('hole', data=0) hole_pump_pulse_file.create_dataset('hole', data=1) - # Add carrier attribute to the pump_pulse_snaps group - elec_pump_pulse_file.create_dataset('carrier', data='electrons') - hole_pump_pulse_file.create_dataset('carrier', data='holes') - # We save the delta_occs_array for animation. If it takes too much space, remove it. elec_delta_occs_array = \ gaussian_excitation(elec_pump_pulse_file, elec_occs_amplitude, From 376bc40f57e83a5f65f75525f65a3f24ee8416ea Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sun, 29 Dec 2024 12:39:10 +0100 Subject: [PATCH 04/21] Docstring for occupation plotting --- .../postproc/utils/spectra_plots.py | 32 +++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index ec2e952..3f2a452 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -280,6 +280,38 @@ def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, scale=0.2): """ Plot the occupations created by a pump pulse on the bands. + Takes ax with the bands plotted and fills the bands with the occupations. + Occupations are set up with a Gaussian for each valence-conduction pair. + + Parameters + ---------- + + ax : matplotlib.axes.Axes + Axis for plotting the bands with occupations. + + kpath : numpy.ndarray + K-path for the bands, same as the one used for plotting the bands. + + bands : numpy.ndarray + Bands for the material. Shape: (num_bands, num_kpoints). Same as the one used for plotting the bands. + + first_el_band_idx : int + Index of the first electron (conduction) band. + + pump_energy : float + Pump energy in eV. + + pump_energy_broadening : float + Pump energy broadening in eV. + + scale : float + Scale factor for the occupation amplitude. For plotting purposes. + + Returns + ------- + + ax : matplotlib.axes.Axes + Axis with the bands filled with the occupations. """ From c3f0fadb75301d4ea8c446af9ec9810eb1c44fce Mon Sep 17 00:00:00 2001 From: imaliyov Date: Mon, 30 Dec 2024 13:03:50 +0100 Subject: [PATCH 05/21] Use energy broadening FWHM and not sigma as input parameter --- .../postproc/utils/spectra_generate_pulse.py | 20 ++++++++++++++----- .../postproc/utils/spectra_plots.py | 1 - 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index c396b44..60439e7 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -25,6 +25,8 @@ def sigma_from_fwhm(fwhm): Comupte Gaussian sigma from the Full Width at Half Maximum (FWHM) parameter: .. math:: f(t) = A / ( sigma * sqrt(2pi)) exp(-(t-t0)^2 / (2sigma^2)) + + For fwhm = 1, sigma = 0.42466. """ sigma = 1.0 / (2.0 * np.sqrt(2.0 * np.log(2.0))) * fwhm @@ -35,8 +37,9 @@ def sigma_from_fwhm(fwhm): def delta_occs_pulse_coef(t, dt, tw, sigma): """ Additional occupation due to the pulse excitation. - Assuming the Gaussian pulse shape, the occupation is increasing in - time according to the error function. + To compute precisely the additional occupation for dt, + we use the integral of the Gaussian function from t - dt to t, + which is the difference of two error functions. Parameters ---------- @@ -50,7 +53,7 @@ def delta_occs_pulse_coef(t, dt, tw, sigma): Time window during which the pump pulse will be applied. sigma : float - Gaussian sigma broadening. + Gaussian sigma for pulse duration. Returns ------- @@ -166,7 +169,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_energy, pump_time_step=1.0, pump_fwhm=20.0, - pump_energy_broadening=0.040, + pump_energy_broadening=0.090, pump_time_window=50.0, finite_width=True, animate=True @@ -177,6 +180,10 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, We use raw k-point and energy arrays as read from the HDF5 files for efficiency. All energies in eV, k points in crystal coordinates. + Pulse duration and energy broadening FWHM are related, typically, + pump_energy_broadening = 1.8 / pump_fwhm. + However, here, we leave the possibility to set them independently. + Parameters ---------- elec_pump_pulse_path : str @@ -204,7 +211,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, The full width at half maximum of the pump pulse (fs). pump_energy_broadening : float, optional - Energy broadening of the pump pulse (eV). + Energy broadening FWHM of the pump pulse (eV). pump_time_window : float, optional The time window for the pump pulse (fs). @@ -292,6 +299,9 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, hole_kpoint_array.view('float64,float64,float64').reshape(-1), return_indices=True)[1:] + # Convert the pump energy FWHM to sigma + pump_energy_broadening_sigma = sigma_from_fwhm(pump_energy_broadening) + # Iteratre over electron and hole bands # Even though, the for loops are inefficient in Python, # the number of bands is rarely more than 10, therefore diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index 3f2a452..b8d9916 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -314,7 +314,6 @@ def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, Axis with the bands filled with the occupations. """ - # Total number of bands nband = len(bands) From 16e02ae814c1cc8aaf63f90274d454f84b325371 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Mon, 30 Dec 2024 13:09:40 +0100 Subject: [PATCH 06/21] Changed variable names for pump duration and spectral width FWHM --- .../postproc/utils/spectra_generate_pulse.py | 28 +++++++++---------- .../postproc/utils/spectra_plots.py | 6 ++-- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 60439e7..f58f920 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -168,8 +168,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, elec_dyna_run, hole_dyna_run, pump_energy, pump_time_step=1.0, - pump_fwhm=20.0, - pump_energy_broadening=0.090, + pump_duration_fwhm=20.0, + pump_spectral_width_fwhm=0.090, pump_time_window=50.0, finite_width=True, animate=True @@ -181,7 +181,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, All energies in eV, k points in crystal coordinates. Pulse duration and energy broadening FWHM are related, typically, - pump_energy_broadening = 1.8 / pump_fwhm. + pump_spectral_width_fwhm = 1.8 / pump_duration_fwhm. However, here, we leave the possibility to set them independently. Parameters @@ -207,10 +207,10 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, Time step in fs for the pump pulse generation. Note: the pump_time_step MUST match the one used in the dynamics run in Perturbo! - pump_fwhm : float, optional + pump_duration_fwhm : float, optional The full width at half maximum of the pump pulse (fs). - pump_energy_broadening : float, optional + pump_spectral_width_fwhm : float, optional Energy broadening FWHM of the pump pulse (eV). pump_time_window : float, optional @@ -226,11 +226,11 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, print(f"{'PUMP PULSE PARAMETERS':*^70}") print(f"{'Pump pulse energy (eV):':>40} {pump_energy:.3f}") - print(f"{'Pump pulse energy broadening (eV):':>40} {pump_energy_broadening:.4f}") + print(f"{'Pump pulse energy broadening (eV):':>40} {pump_spectral_width_fwhm:.4f}") print(f"{'Finite width:':>40} {finite_width}") if finite_width: print(f"{'Pump pulse time step (fs):':>40} {pump_time_step:.3f}") - print(f"{'Pump pulse FWHM (fs):':>40} {pump_fwhm:.3f}") + print(f"{'Pump pulse FWHM (fs):':>40} {pump_duration_fwhm:.3f}") print(f"{'Pump pulse time window (fs):':>40} {pump_time_window:.3f}") print("") @@ -300,7 +300,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, return_indices=True)[1:] # Convert the pump energy FWHM to sigma - pump_energy_broadening_sigma = sigma_from_fwhm(pump_energy_broadening) + pump_energy_broadening_sigma = sigma_from_fwhm(pump_spectral_width_fwhm) # Iteratre over electron and hole bands # Even though, the for loops are inefficient in Python, @@ -311,7 +311,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, # diff. between elec and hole for a given elec branch iband and hole branch jband delta = pump_factor * \ spectra_plots.gaussian(elec_energy_array[ekidx, iband] - hole_energy_array[hkidx, jband], - pump_energy, pump_energy_broadening, hole_nband, elec_nband) + pump_energy, pump_spectral_width_fwhm, hole_nband, elec_nband) # Only for the intersected k points, we add the delta occupation elec_occs_amplitude[ekidx, iband] += delta @@ -347,7 +347,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, optional_params[0] = pump_factor if finite_width: - optional_params[1] = pump_fwhm + optional_params[1] = pump_duration_fwhm for h5f in [elec_pump_pulse_file, hole_pump_pulse_file]: h5f.create_group('pump_pulse_snaps') @@ -355,8 +355,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, h5f.create_dataset('pump_energy', data=pump_energy) h5f['pump_energy'].attrs['units'] = 'eV' - h5f.create_dataset('pump_energy_broadening', data=pump_energy_broadening) - h5f['pump_energy_broadening'].attrs['units'] = 'eV' + h5f.create_dataset('pump_spectral_width_fwhm', data=pump_spectral_width_fwhm) + h5f['pump_spectral_width_fwhm'].attrs['units'] = 'eV' h5f['optional_params'] = optional_params @@ -385,13 +385,13 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, elec_delta_occs_array = \ gaussian_excitation(elec_pump_pulse_file, elec_occs_amplitude, time_grid, - pump_time_window, pump_fwhm, pump_time_step, + pump_time_window, pump_duration_fwhm, pump_time_step, hole=False, finite_width=finite_width) hole_delta_occs_array = \ gaussian_excitation(hole_pump_pulse_file, hole_occs_amplitude, time_grid, - pump_time_window, pump_fwhm, pump_time_step, + pump_time_window, pump_duration_fwhm, pump_time_step, hole=True, finite_width=finite_width) # Close diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index b8d9916..3ab2b33 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -276,7 +276,7 @@ def plot_trans_abs_map(ax, time_grid, energy_grid, trans_abs, def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, - pump_energy, pump_energy_broadening, + pump_energy, pump_spectral_width_fwhm, scale=0.2): """ Plot the occupations created by a pump pulse on the bands. @@ -301,7 +301,7 @@ def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, pump_energy : float Pump energy in eV. - pump_energy_broadening : float + pump_spectral_width_fwhm : float Pump energy broadening in eV. scale : float @@ -336,7 +336,7 @@ def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, hole_band = bands[i_hole_band] occ = gaussian(elec_band - hole_band, - pump_energy, pump_energy_broadening, hole_nband, elec_nband) + pump_energy, pump_spectral_width_fwhm, hole_nband, elec_nband) occs_amplitude_bands[i_elec_band - 1, :] += occ occs_amplitude_bands[i_hole_band - 1, :] += occ From df268ab94ee0735f127469ca61f0f50a88bb2b69 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Wed, 1 Jan 2025 13:58:33 +0100 Subject: [PATCH 07/21] PumpPulse advanced functionality --- setup.cfg | 2 +- src/perturbopy/postproc/__init__.py | 2 +- .../postproc/calc_modes/dyna_run.py | 102 ++++++++++++++++-- src/perturbopy/postproc/dbs/units_dict.py | 3 +- .../postproc/utils/spectra_generate_pulse.py | 84 +++++++++++---- .../postproc/utils/spectra_plots.py | 74 ++++++++++++- 6 files changed, 230 insertions(+), 37 deletions(-) diff --git a/setup.cfg b/setup.cfg index 617006d..73a5e3b 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,3 +1,3 @@ [flake8] max-line-length = 160 -ignore = E114, E127, E124, E221, W293, E721, W503, F841, F401, E202, W504 +ignore = E114, E127, E124, E221, W293, E721, W503, F841, F401, E202, W504, E731 diff --git a/src/perturbopy/postproc/__init__.py b/src/perturbopy/postproc/__init__.py index fb521c0..56eb6d2 100644 --- a/src/perturbopy/postproc/__init__.py +++ b/src/perturbopy/postproc/__init__.py @@ -12,7 +12,7 @@ from .calc_modes.trans import Trans from .calc_modes.imsigma import Imsigma from .calc_modes.imsigma_spin import ImsigmaSpin -from .calc_modes.dyna_run import DynaRun +from .calc_modes.dyna_run import DynaRun, PumpPulse from .calc_modes.dyna_pp import DynaPP from .dbs.units_dict import UnitsDict diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 36e1519..99dd6ba 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -1,5 +1,6 @@ -import numpy as np import os +import numpy as np +import warnings from perturbopy.postproc.calc_modes.calc_mode import CalcMode from perturbopy.postproc.dbs.recip_pt_db import RecipPtDB from perturbopy.postproc.calc_modes.dyna_indiv_run import DynaIndivRun @@ -7,6 +8,12 @@ from perturbopy.postproc.utils.timing import Timing, TimingGroup from perturbopy.postproc.dbs.units_dict import UnitsDict +# for plotting +import matplotlib.pyplot as plt +from perturbopy.postproc.utils.spectra_plots import find_fwhm +from perturbopy.postproc.utils.plot_tools import plotparams +plt.rcParams.update(plotparams) + class DynaRun(CalcMode): """ @@ -256,15 +263,24 @@ class PumpPulse(): pump_energy : float Energy of the pump pulse excitation in eV. - energy_broadening : float - Energy broadening of the pump pulse excitation in eV. + spectral_width_fwhm : float + Energy broadening FWHM of the pump pulse excitation in eV. + + pump_duration_fwhm : float + Duration FWHM of the pump pulse excitation in fs. num_steps : int Number of steps in the pump pulse excitation. - time_step : float + pump_factor : float + Factor for the pump pulse excitation. + + pump_time_step : float Time step of the pump pulse excitation in fs. + num_bands : int + Number of bands in the pump pulse excitation. Tailored for the Perturbo dynamics-run calculation. + num_kpoints : int Number of k-points in the pump pulse excitation. Tailored for the Perturbo dynamics-run calculation. @@ -278,6 +294,12 @@ class PumpPulse(): hole : bool Flag to indicate if the pump pulse excitation is for the hole. Must be the same as in the ultrafast simulation. + + time_profile : np.ndarray + Array of time profile for the pump pulse excitation. First column is time in fs, second column is the time profile. + + energy_profile : np.ndarray + Array of energy profile for the pump pulse excitation. First column is energy in eV, second column is the energy profile """ def __init__(self, pump_dict): @@ -298,14 +320,20 @@ def __init__(self, pump_dict): self.pump_energy = pump_dict['pump_energy'] self.pump_energy_units = pump_dict['pump_energy units'] - self.energy_broadening = pump_dict['energy_broadening'] - self.energy_broadening_units = pump_dict['energy_broadening units'] + self.spectral_width_fwhm = pump_dict['pump_spectral_width_fwhm'] + self.spectral_width_fwhm_units = pump_dict['pump_spectral_width_fwhm units'] + + self.pump_duration_fwhm = pump_dict['pump_duration_fwhm'] + self.pump_duration_fwhm_units = pump_dict['pump_duration_fwhm units'] self.num_steps = pump_dict['num_steps'] - self.time_step = pump_dict['time_step'] - self.time_step_units = pump_dict['time_step units'] + self.pump_factor = pump_dict['pump_factor'] + self.pump_time_step = pump_dict['pump_time_step'] + self.pump_time_step_units = pump_dict['pump_time_step units'] + + self.num_bands = pump_dict['num_bands'] self.num_kpoints = pump_dict['num_kpoints'] self.carrier_number_array = \ @@ -318,6 +346,10 @@ def __init__(self, pump_dict): self.hole = pump_dict['hole'] + # arrays of time and energy profile + self.time_profile = np.array(pump_dict['time_profile']) + self.energy_profile = np.array(pump_dict['energy_profile']) + def __str__(self): """ Method to print the pump pulse excitation parameters. @@ -325,9 +357,59 @@ def __str__(self): text = 'Pump pulse excitation parameters:\n' text += f"{'Pump energy':>30}: {self.pump_energy} {self.pump_energy_units}\n" - text += f"{'Energy broadening':>30}: {self.energy_broadening} {self.energy_broadening_units}\n" + text += f"{'Pump duration FWHM':>30}: {self.pump_duration_fwhm} {self.pump_duration_fwhm_units}\n" + text += f"{'Energy broadening FWHM':>30}: {self.spectral_width_fwhm} {self.spectral_width_fwhm_units}\n" + text += f"{'Pump factor':>30}: {self.pump_factor}\n" text += f"{'Number of steps':>30}: {self.num_steps}\n" - text += f"{'Time step':>30}: {self.time_step} {self.time_step_units}\n" + text += f"{'Time step':>30}: {self.pump_time_step} {self.pump_time_step_units}\n" text += f"{'Hole':>30}: {self.hole}\n" return text + + def plot_time_profile(self, ax=None): + """ + Plot the pump pulse time profile. + """ + + if self.time_profile is None: + warnings.warn('No time profile found for the pump pulse') + return None + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=(10, 8)) + + # Find the FWHM and half-maximum of the time profile + time_left_FWHM, time_right_FWHM, time_half_max = find_fwhm(self.time_profile[:, 0], self.time_profile[:, 1]) + + ax.plot(self.time_profile[:, 0], self.time_profile[:, 1]) + ax.plot([time_left_FWHM, time_right_FWHM], [time_half_max, time_half_max], marker='o', color='tab:red', lw=3, label='FWHM') + ax.set_xlabel('Time (fs)', fontsize=24) + ax.set_ylabel('Time Gaussian', fontsize=24) + ax.legend() + + return ax + + def plot_energy_profile(self, ax=None): + """ + Plot the pump pulse energy profile. + """ + + if self.energy_profile is None: + warnings.warn('No energy profile found for the pump pulse') + return None + + if ax is None: + fig, ax = plt.subplots(1, 1, figsize=(10, 8)) + + # Find the FWHM and half-maximum of the energy profile + energy_left_FWHM, energy_right_FWHM, energy_half_max = find_fwhm(self.energy_profile[:, 0], self.energy_profile[:, 1]) + + ax.plot(self.energy_profile[:, 0], self.energy_profile[:, 1]) + ax.plot([energy_left_FWHM, energy_right_FWHM], [energy_half_max, energy_half_max], marker='o', color='tab:red', lw=3, label='FWHM') + ax.axvline(self.pump_energy, color='gray', lw=3, label='Pump energy') + ax.set_xlabel('Energy (eV)', fontsize=24) + ax.set_ylabel('Energy Gaussian', fontsize=24) + ax.legend(loc='upper right') + + return ax + diff --git a/src/perturbopy/postproc/dbs/units_dict.py b/src/perturbopy/postproc/dbs/units_dict.py index e2231f9..7b97d75 100644 --- a/src/perturbopy/postproc/dbs/units_dict.py +++ b/src/perturbopy/postproc/dbs/units_dict.py @@ -4,7 +4,8 @@ class UnitsDict(dict): """ - This is a class representation of a set of physical quantities with units organized in a dictionary. + Modified dictionary class to store physical quantities with units. + Class representation of a set of physical quantities with units organized in a dictionary. The physical quantities may either be arrays or floats. Attributes diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index f58f920..d1c0e39 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -13,6 +13,8 @@ from matplotlib.animation import FuncAnimation +from perturbopy.postproc import PumpPulse + from .memory import get_size from .timing import TimingGroup from .constants import energy_conversion_factor @@ -69,7 +71,7 @@ def delta_occs_pulse_coef(t, dt, tw, sigma): return 0.5 * numer / denom -def gaussian_excitation(pump_file, occs_amplitude, time_grid, time_window, pulse_fwhm, time_step, +def gaussian_excitation(pump_file, occs_amplitude, time_grid, time_window, pump_duration, time_step, hole=False, finite_width=True): """ The Gaussian pulse excitation for the Perturbo dynamics. @@ -91,7 +93,7 @@ def gaussian_excitation(pump_file, occs_amplitude, time_grid, time_window, pulse time_window : float The time window for the pump pulse (fs). During this time window, the pump pulse will be active. - pulse_fwhm : float + pump_duration : float The full width at half maximum of the pump pulse (fs). time_step : float @@ -122,7 +124,7 @@ def gaussian_excitation(pump_file, occs_amplitude, time_grid, time_window, pulse num_steps = time_grid.shape[0] # Compute the Gaussian sigma from the full width at half maximum (FWHM) - sigma = sigma_from_fwhm(pulse_fwhm) + sigma = sigma_from_fwhm(pump_duration) # Number of k-points and bands num_k, num_bands = occs_amplitude.shape @@ -138,10 +140,10 @@ def gaussian_excitation(pump_file, occs_amplitude, time_grid, time_window, pulse delta_occs_array = np.zeros((num_bands, num_steps, num_k)) # A pre-factor for the pump pulse, num_steps steps - pulse_coef = delta_occs_pulse_coef(time_grid[:], time_step, time_window, sigma) + pump_time_profile = delta_occs_pulse_coef(time_grid[:], time_step, time_window, sigma) # t - time, j - band, k - k-point - delta_occs_array[:, :, :] = np.einsum('t,jk->ktj', pulse_coef[:], occs_amplitude[:, :]) + delta_occs_array[:, :, :] = np.einsum('t,jk->ktj', pump_time_profile[:], occs_amplitude[:, :]) # Print the size of the array since it can be large: num_bands * num_steps * num_k get_size(delta_occs_array, f'delta(fnk(t)) {carrier}', dump=True) @@ -171,9 +173,9 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_duration_fwhm=20.0, pump_spectral_width_fwhm=0.090, pump_time_window=50.0, + pump_factor=0.3, finite_width=True, - animate=True - ): + animate=True): """ Setup the Gaussian pump pulse excitation for electrons and holes. Write into the pump_pulse.h5 HDF5 file. @@ -219,6 +221,9 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, The number of the snapshots in the pump_pulse.h5 file will be equal to pump_time_window / pump_time_step. + pump_factor : float, optional + The amplitude of the pump pulse. The maximum occupation will be pump_factor. + finite_width : bool, optional If True, the pulse is finite in time. If False, the pulse is a step function and occs_amplitude will be set as initial occupation. @@ -282,9 +287,6 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, elec_occs_amplitude = np.zeros_like(elec_energy_array) hole_occs_amplitude = np.zeros_like(hole_energy_array) - # A pre-factor for the pump pulse - pump_factor = 1.0 - print('Computing the optically excited occupations...') t_occs = TimingGroup('Setup occupations') t_occs.add('total', level=3).start() @@ -311,7 +313,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, # diff. between elec and hole for a given elec branch iband and hole branch jband delta = pump_factor * \ spectra_plots.gaussian(elec_energy_array[ekidx, iband] - hole_energy_array[hkidx, jband], - pump_energy, pump_spectral_width_fwhm, hole_nband, elec_nband) + pump_energy, pump_energy_broadening_sigma) # Only for the intersected k points, we add the delta occupation elec_occs_amplitude[ekidx, iband] += delta @@ -329,6 +331,11 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, print(f"{'Difference:':>40} {sum(elec_occs_amplitude.ravel()) - sum(hole_occs_amplitude.ravel()):.4e}") print("") + # Create the energy profile, only for reference + en_grid = np.linspace(pump_energy - 3 * pump_energy_broadening_sigma, + pump_energy + 3 * pump_energy_broadening_sigma, 200) + en_profile = spectra_plots.gaussian(en_grid, pump_energy, pump_energy_broadening_sigma) + # Write to HDF5 file, replacing snap_t_1 # Electron if os.path.exists(elec_pump_pulse_path): @@ -343,11 +350,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, hole_pump_pulse_file = open_hdf5(hole_pump_pulse_path, mode='w') # Optional pump pulse parameters specific to a given shape of pulse - optional_params = np.zeros(10, dtype=float) - optional_params[0] = pump_factor - - if finite_width: - optional_params[1] = pump_duration_fwhm + optional_params = np.ones(10, dtype=float) * (-9999) for h5f in [elec_pump_pulse_file, hole_pump_pulse_file]: h5f.create_group('pump_pulse_snaps') @@ -358,20 +361,34 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, h5f.create_dataset('pump_spectral_width_fwhm', data=pump_spectral_width_fwhm) h5f['pump_spectral_width_fwhm'].attrs['units'] = 'eV' - h5f['optional_params'] = optional_params + h5f.create_dataset('pump_factor', data=pump_factor) + h5f.create_dataset('optional_params', data=optional_params) if finite_width: h5f.create_dataset('time_window', data=pump_time_window) h5f.create_dataset('num_steps', data=num_steps) - h5f.create_dataset('time_step', data=pump_time_step) + h5f.create_dataset('pump_time_step', data=pump_time_step) + h5f.create_dataset('pump_duration_fwhm', data=pump_duration_fwhm) + + pump_time_profile = delta_occs_pulse_coef(time_grid[:], pump_time_step, + pump_time_window, sigma_from_fwhm(pump_duration_fwhm)) + h5f.create_dataset('time_profile', data=np.c_[time_grid, pump_time_profile]) else: + pump_time_profile = None h5f.create_dataset('time_window', data=elec_dyna_time_step) h5f.create_dataset('num_steps', data=1) - h5f.create_dataset('time_step', data=elec_dyna_time_step) + h5f.create_dataset('pump_time_step', data=elec_dyna_time_step) + h5f.create_dataset('pump_duration_fwhm', data=elec_dyna_time_step) + h5f.create_dataset('time_profile', data=np.c_[0, 1]) + h5f['pump_duration_fwhm'].attrs['units'] = 'fs' h5f['time_window'].attrs['units'] = 'fs' - h5f['time_step'].attrs['units'] = 'fs' + h5f['pump_time_step'].attrs['units'] = 'fs' + + h5f.create_dataset('energy_profile', data=np.c_[en_grid, en_profile]) + h5f['energy_profile'].attrs['units'] = 'eV' + h5f['time_profile'].attrs['units'] = 'fs' elec_pump_pulse_file.create_dataset('num_bands', data=elec_nband) hole_pump_pulse_file.create_dataset('num_bands', data=hole_nband) @@ -409,3 +426,30 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, elec_delta_occs_array, elec_kpoint_array, elec_energy_array, hole_delta_occs_array, hole_kpoint_array, hole_energy_array, pump_energy) + + # Setup the pump pulse object + pump_dict = { + 'pump_energy': pump_energy, + 'pump_energy units': 'eV', + 'pump_spectral_width_fwhm': pump_spectral_width_fwhm, + 'pump_spectral_width_fwhm units': 'eV', + 'pump_duration_fwhm': pump_duration_fwhm, + 'pump_duration_fwhm units': 'fs', + 'num_steps': num_steps, + 'pump_factor': pump_factor, + 'pump_time_step': pump_time_step, + 'pump_time_step units': 'fs', + 'num_bands': elec_nband, + 'num_kpoints': elec_kpoint_array.shape[0], + 'optional_params': optional_params, + 'hole': None, + 'time_profile': np.c_[time_grid, pump_time_profile], + 'energy_profile': np.c_[en_grid, en_profile], + # These quantities will be computed at the perturbo.x stage + 'pump pulse carrier number': None, + 'carrier_number units': None, + } + + pump_pulse = PumpPulse(pump_dict) + + return pump_pulse diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index 3ab2b33..7da1f39 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -7,7 +7,8 @@ import numpy as np import matplotlib.pyplot as plt import warnings -from scipy import special +from scipy.interpolate import Akima1DInterpolator +from scipy.optimize import brentq from perturbopy.io_utils.io import open_hdf5, close_hdf5 @@ -26,12 +27,77 @@ plt.rcParams.update(plotparams) -def gaussian(x, mu, sig, hole_nband, elec_nband): +def gaussian(x, mu, sigma): """ - Gaussian function normalized to unity max occupation + Gaussian function. Not normalized. """ - return np.exp(-0.5 * ((x - mu) / sig)**2) / (hole_nband * elec_nband) + return np.exp(-0.5 * ((x - mu) / sigma)**2) + + +def find_fwhm(x, y, num_interp_points=2000): + """ + Find the Full Width at Half Maximum (FWHM) for a single prominent peak y(x), + using an Akima1DInterpolator to create a smooth spline, and root-finding + (brentq) for a precise half-max crossing. + + Parameters + ---------- + x : array_like + 1D array of x-values (assumed sorted in ascending order). + y : array_like + 1D array of y-values corresponding to x. + num_interp_points : int, optional + Number of points for evaluating the Akima spline; default is 2000. + + Returns + ------- + x_left : float + x-value at the left FWHM crossing. + + x_right : float + x-value at the right FWHM crossing. + + half_max : float + Half-maximum value of the peak. + """ + + # Create an Akima spline over the original data + akima = Akima1DInterpolator(x, y) + + # Evaluate on a fine grid to locate the approximate peak + x_fine = np.linspace(x[0], x[-1], num_interp_points) + y_fine = akima(x_fine) + + # 1. Find approximate peak index on the fine grid + i_max = np.argmax(y_fine) + x_max = x_fine[i_max] + y_max = y_fine[i_max] + half_max = y_max / 2.0 + + # f(x) = akima(x) - half_max + func = lambda xx: akima(xx) - half_max + + # --- Locate x_left by scanning left from the peak for sign change --- + x_left = x[0] # fallback in case no crossing is found + for i in range(i_max, 0, -1): + if y_fine[i - 1] - half_max < 0 < y_fine[i] - half_max or \ + y_fine[i - 1] - half_max > 0 > y_fine[i] - half_max: + # We found a bracket where f() changes sign + x_left = brentq(func, x_fine[i - 1], x_fine[i]) + break + + # --- Locate x_right by scanning right from the peak for sign change --- + x_right = x[-1] # fallback in case no crossing is found + for i in range(i_max, len(x_fine) - 1): + if y_fine[i] - half_max < 0 < y_fine[i + 1] - half_max or \ + y_fine[i] - half_max > 0 > y_fine[i + 1] - half_max: + # We found a bracket where f() changes sign + x_right = brentq(func, x_fine[i], x_fine[i + 1]) + break + + # The y-values at x_left and x_right are both half_max + return x_left, x_right, half_max def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, From 951d051fd69a931afb3ba0bbf46521b0c220eea3 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Wed, 1 Jan 2025 18:34:30 +0100 Subject: [PATCH 08/21] Create a cnum_check folder for the precise carrier number calculation with dynamics-pp --- .../postproc/calc_modes/dyna_run.py | 102 ++++++++++++++++++ .../postproc/utils/spectra_generate_pulse.py | 40 ++++++- 2 files changed, 139 insertions(+), 3 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 99dd6ba..7f04599 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -252,6 +252,108 @@ def extract_steady_drift_vel(self, dyna_pp_yaml_path): return steady_drift_vel, steady_conc + @staticmethod + def to_cdyna_h5(prefix, band_structure_ryd, snap_array, time_step_fs, + path='.', new=True, overwrite=False, num_runs=1): + """ + Write the dynamics data into the prefix_cdyna.h5 HDF5 file. + Follows the script format of the Perturbo dynamics-run calculation. + Currently implemented as a static method. The idea is that if we need to write an HDF5 cdyna + file, it is forcfully different from an existing one, therefore, one or more + attributes are different from the original cdyna file. + + Parameters + ---------- + + prefix : str + Prefix of the HDF5 file. Same as prefix in any Perturbo calculation. + + band_structure_ryd : np.ndarray + Band structure in Rydberg units. Shape: num_kpoints x num_bands + + snap_array : np.ndarray + Array of carrier occupations. Shape: (num_steps, (num_kpoints x num_bands)) + + time_step_fs : float or np.ndarray + Time step in fs. If float, the same time step is used for all runs. + If np.ndarray, the time steps for each run are specified. + + path : str + Path to the directory where the HDF5 file is saved. Default is the current directory. + + new : bool + Flag to indicate if the file is new. Default is True. + + overwrite : bool + Flag to indicate if the file is overwritten. Default is False. + + num_runs : int + Number of dynamics_run runs. Default is 1. + """ + + trun = TimingGroup('write cdyna') + trun.add('total', level=3).start() + + if not new: + raise NotImplementedError("Only new files are supported") + if num_runs != 1: + raise NotImplementedError("Only num_runs=1 is supported") + + if isinstance(time_step_fs, (float, np.floating)): + time_step_fs = np.array([time_step_fs]) + if time_step_fs.shape[0] != num_runs: + raise ValueError("The number of time steps must be equal to num_runs") + + if path != '.': + os.makedirs(path, exist_ok=True) + + new_cdyna_filename = os.path.join(path, f'{prefix}_cdyna.h5') + + if os.path.isfile(new_cdyna_filename): + if overwrite: + warnings.warn(f'File {new_cdyna_filename} already exists. Overwriting it.') + else: + raise FileExistsError(f'File {new_cdyna_filename} already exists. Set overwrite=True to overwrite it.') + + new_cdyna_file = open_hdf5(new_cdyna_filename, 'w') + + new_cdyna_file.create_dataset('band_structure_ryd', data=band_structure_ryd) + new_cdyna_file['band_structure_ryd'].attrs['ryd2ev'] = 13.605698066 + + new_cdyna_file.create_dataset('num_runs', data=num_runs) + + for irun in range(1, num_runs + 1): + + # In Perturbo, in dynamics_run_1, snap_t_ start from 0, where snap_t_0 is the initial state + # In dynamics_run_2 and onwards, snap_t_ start from + if irun == 1: + offset = 0 + + if snap_array.shape[0] == 1: + # duplicate the initial state + snap_array = np.vstack((snap_array, snap_array)) + num_steps = snap_array.shape[0] - 1 + + else: + offset = 1 + num_steps = snap_array.shape[0] + + dyn_str = f'dynamics_run_{irun}' + new_cdyna_file.create_group(dyn_str) + + + new_cdyna_file[dyn_str].create_dataset('num_steps', data=num_steps) + new_cdyna_file[dyn_str].create_dataset('time_step_fs', data=time_step_fs[irun - 1]) + + for itime in range(snap_array.shape[0]): + new_cdyna_file[dyn_str].create_dataset(f'snap_t_{itime + offset}', data=snap_array[itime, :, np.newaxis]) + + close_hdf5(new_cdyna_file) + + trun.timings['total'].stop() + print(f'{"Total time to write cdyna":>30}: {trun.timings["total"].total_runtime} s') + print() + class PumpPulse(): """ diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index d1c0e39..9dda86a 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -11,8 +11,6 @@ from perturbopy.io_utils.io import open_hdf5, close_hdf5 -from matplotlib.animation import FuncAnimation - from perturbopy.postproc import PumpPulse from .memory import get_size @@ -175,7 +173,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_time_window=50.0, pump_factor=0.3, finite_width=True, - animate=True): + animate=True, + cnum_check=True): """ Setup the Gaussian pump pulse excitation for electrons and holes. Write into the pump_pulse.h5 HDF5 file. @@ -227,6 +226,14 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, finite_width : bool, optional If True, the pulse is finite in time. If False, the pulse is a step function and occs_amplitude will be set as initial occupation. + + animate : bool, optional + If True, animate the pump pulse excitation. + + cnum_check : bool, optional + If True, create a cnum_check folder with a prefix_cdyna.h5 file with the total + occupation summed over time steps. Run dynamics-pp with Perturbo to get the accurate + total carrier number in the cnum_check folder (do not forget to link the epr file there). """ print(f"{'PUMP PULSE PARAMETERS':*^70}") @@ -452,4 +459,31 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_pulse = PumpPulse(pump_dict) + # Create the cnum_check folder + if cnum_check: + + if finite_width: + total_occ_elec = np.sum(elec_delta_occs_array, axis=1) + total_occ_hole = np.sum(hole_delta_occs_array, axis=1) + + else: + total_occ_elec = elec_occs_amplitude + total_occ_hole = hole_occs_amplitude + + elec_cnum_check_path = os.path.join(os.path.dirname(elec_pump_pulse_path), 'cnum_check') + hole_cnum_check_path = os.path.join(os.path.dirname(hole_pump_pulse_path), 'cnum_check') + elec_dyna_run.to_cdyna_h5(elec_dyna_run.prefix, + elec_dyna_run._energies, + total_occ_elec, + elec_dyna_run[1].time_step, + path=elec_cnum_check_path, + overwrite=True) + + hole_dyna_run.to_cdyna_h5(hole_dyna_run.prefix, + hole_dyna_run._energies, + total_occ_hole, + hole_dyna_run[1].time_step, + path=hole_cnum_check_path, + overwrite=True) + return pump_pulse From 89b2e576509531c49af9a5b6f01cebe72f3f089c Mon Sep 17 00:00:00 2001 From: imaliyov Date: Wed, 1 Jan 2025 18:37:09 +0100 Subject: [PATCH 09/21] flake8 fixes --- src/perturbopy/postproc/calc_modes/dyna_run.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 7f04599..59bf176 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -340,8 +340,6 @@ def to_cdyna_h5(prefix, band_structure_ryd, snap_array, time_step_fs, dyn_str = f'dynamics_run_{irun}' new_cdyna_file.create_group(dyn_str) - - new_cdyna_file[dyn_str].create_dataset('num_steps', data=num_steps) new_cdyna_file[dyn_str].create_dataset('time_step_fs', data=time_step_fs[irun - 1]) @@ -514,4 +512,3 @@ def plot_energy_profile(self, ax=None): ax.legend(loc='upper right') return ax - From 07a33ac379060744ff62a3c194901e452996571c Mon Sep 17 00:00:00 2001 From: imaliyov Date: Wed, 1 Jan 2025 19:15:04 +0100 Subject: [PATCH 10/21] plot_scale factor for plotting occupations --- .../postproc/calc_modes/dyna_run.py | 2 ++ .../postproc/utils/spectra_generate_pulse.py | 9 +++++-- .../postproc/utils/spectra_plots.py | 26 +++++++++++++------ 3 files changed, 27 insertions(+), 10 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 59bf176..5b1c0aa 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -486,6 +486,7 @@ def plot_time_profile(self, ax=None): ax.set_xlabel('Time (fs)', fontsize=24) ax.set_ylabel('Time Gaussian', fontsize=24) ax.legend() + ax.grid() return ax @@ -510,5 +511,6 @@ def plot_energy_profile(self, ax=None): ax.set_xlabel('Energy (eV)', fontsize=24) ax.set_ylabel('Energy Gaussian', fontsize=24) ax.legend(loc='upper right') + ax.grid() return ax diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 9dda86a..26724e0 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -174,6 +174,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_factor=0.3, finite_width=True, animate=True, + plot_scale=1.0, cnum_check=True): """ Setup the Gaussian pump pulse excitation for electrons and holes. @@ -230,6 +231,10 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, animate : bool, optional If True, animate the pump pulse excitation. + plot_scale : float + Scale factor for the scatter object sizes for carrier occupations. + Default is 1.0. + cnum_check : bool, optional If True, create a cnum_check folder with a prefix_cdyna.h5 file with the total occupation summed over time steps. Run dynamics-pp with Perturbo to get the accurate @@ -425,14 +430,14 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, print('\nPlotting...') spectra_plots.plot_occ_ampl(elec_occs_amplitude, elec_kpoint_array, elec_energy_array, hole_occs_amplitude, hole_kpoint_array, - hole_energy_array, pump_energy) + hole_energy_array, pump_energy, plot_scale=plot_scale) if animate and finite_width: print('\nAnimating...') spectra_plots.animate_pump_pulse(pump_time_step, elec_delta_occs_array, elec_kpoint_array, elec_energy_array, hole_delta_occs_array, hole_kpoint_array, hole_energy_array, - pump_energy) + pump_energy, plot_scale=plot_scale) # Setup the pump pulse object pump_dict = { diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index 7da1f39..12660b7 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -101,7 +101,7 @@ def find_fwhm(x, y, num_interp_points=2000): def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, - h_occs, hole_kpoint_array, hole_energy_array, pump_energy): + h_occs, hole_kpoint_array, hole_energy_array, pump_energy, plot_scale=1.0): """ Plot occupation amplitude. Currently hardcoded to the section kz = 0. @@ -127,6 +127,9 @@ def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, pump_energy: float Pump energy in eV. + + plot_scale : float + Scale factor for the scatter object sizes. Default is 1.0. """ # find where kz == 0 @@ -137,13 +140,13 @@ def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, idx = np.where((elec_kpoint_array[:, 1] == 0) & (elec_kpoint_array[:, 0] == 0)) for i in range(np.size(elec_energy_array, axis=1)): ax.scatter(elec_kpoint_array[idx, 2], elec_energy_array[idx, i], s=0.5, c='black', alpha=0.5) - ax.scatter(elec_kpoint_array[idx, 2], elec_energy_array[idx, i], s=e_occs[idx, i][0] * 1000 + 1e-10, c='red', alpha=0.5) + ax.scatter(elec_kpoint_array[idx, 2], elec_energy_array[idx, i], s=e_occs[idx, i][0] * 1e4 * plot_scale + 1e-10, c='red', alpha=0.5) # Plot hole occupations idx = np.where((hole_kpoint_array[:, 1] == 0) & (hole_kpoint_array[:, 0] == 0)) for i in range(np.size(hole_energy_array, axis=1)): ax.scatter(hole_kpoint_array[idx, 2], hole_energy_array[idx, i], s=0.5, c='black', alpha=0.5) - ax.scatter(hole_kpoint_array[idx, 2], hole_energy_array[idx, i], s=h_occs[idx, i][0] * 1000 + 1e-10, c='red', + ax.scatter(hole_kpoint_array[idx, 2], hole_energy_array[idx, i], s=h_occs[idx, i][0] * 1e4 * plot_scale + 1e-10, c='red', alpha=0.5) fsize = 16 @@ -157,7 +160,8 @@ def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, def animate_pump_pulse(time_step, elec_delta_occs_array, elec_kpoint_array, elec_energy_array, hole_delta_occs_array, hole_kpoint_array, hole_energy_array, - pump_energy): + pump_energy, + plot_scale=1.0): """ Animate the pump pulse excitation for electrons and holes. Defines fig and ax, initializes scatter objects for electron and hole occupations, and calls update_scatter. @@ -188,6 +192,9 @@ def animate_pump_pulse(time_step, pump_energy: float Pump energy in eV, used only in title. + + plot_scale : float + Scale factor for the scatter object sizes. Default is 1.0. """ fig, ax = plt.subplots(1, 1, figsize=(9, 6)) @@ -217,7 +224,7 @@ def animate_pump_pulse(time_step, ani = FuncAnimation(fig, update_scatter, frames=elec_delta_occs_array.shape[1], fargs=(ax, time_step, idx_elec, idx_hole, elec_scat_list, hole_scat_list, - elec_delta_occs_array, hole_delta_occs_array), + elec_delta_occs_array, hole_delta_occs_array, plot_scale), interval=100, repeat=False) # Save the animation to gif @@ -228,7 +235,7 @@ def animate_pump_pulse(time_step, def update_scatter(anim_time, ax, time_step, idx_elec, idx_hole, elec_scat_list, hole_scat_list, - elec_delta_occs_array, hole_delta_occs_array): + elec_delta_occs_array, hole_delta_occs_array, plot_scale=1.0): """ Animate the pump pulse excitation for electrons and holes. @@ -261,16 +268,19 @@ def update_scatter(anim_time, ax, time_step, idx_elec, idx_hole, hole_delta_occs_array : numpy.ndarray Array of hole occupation changes, similar to elec_delta_occs_array. + + plot_scale : float + Scale factor for the scatter object sizes. Default is 1.0. """ elec_num_bands = elec_delta_occs_array.shape[0] hole_num_bands = hole_delta_occs_array.shape[0] for i in range(elec_num_bands): - elec_scat_list[i].set_sizes(elec_delta_occs_array[i, anim_time, idx_elec].flatten() * 2e4 + 1e-10) + elec_scat_list[i].set_sizes(elec_delta_occs_array[i, anim_time, idx_elec].flatten() * 2e4 * plot_scale + 1e-10) for i in range(hole_num_bands): - hole_scat_list[i].set_sizes(hole_delta_occs_array[i, anim_time, idx_hole].flatten() * 2e4 + 1e-10) + hole_scat_list[i].set_sizes(hole_delta_occs_array[i, anim_time, idx_hole].flatten() * 2e4 * plot_scale + 1e-10) suffix = ax.get_title().split(';')[0] ax.set_title(f'{suffix}; Time: {anim_time * time_step:.2f} fs') From 3538de2a7f694b875c633cde35e4d010cc28ae80 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sat, 4 Jan 2025 12:22:12 +0100 Subject: [PATCH 11/21] Better scale factor for occupation plotting --- .../postproc/calc_modes/dyna_run.py | 10 ++++--- .../postproc/utils/spectra_generate_pulse.py | 10 ------- .../postproc/utils/spectra_plots.py | 28 +++++++++++-------- 3 files changed, 22 insertions(+), 26 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 5b1c0aa..c88a3b1 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -456,12 +456,12 @@ def __str__(self): """ text = 'Pump pulse excitation parameters:\n' - text += f"{'Pump energy':>30}: {self.pump_energy} {self.pump_energy_units}\n" - text += f"{'Pump duration FWHM':>30}: {self.pump_duration_fwhm} {self.pump_duration_fwhm_units}\n" - text += f"{'Energy broadening FWHM':>30}: {self.spectral_width_fwhm} {self.spectral_width_fwhm_units}\n" + text += f"{'Pump energy':>30}: {self.pump_energy:.4f} {self.pump_energy_units}\n" + text += f"{'Pump duration FWHM':>30}: {self.pump_duration_fwhm:.2f} {self.pump_duration_fwhm_units}\n" + text += f"{'Energy broadening FWHM':>30}: {self.spectral_width_fwhm:.4f} {self.spectral_width_fwhm_units}\n" text += f"{'Pump factor':>30}: {self.pump_factor}\n" text += f"{'Number of steps':>30}: {self.num_steps}\n" - text += f"{'Time step':>30}: {self.pump_time_step} {self.pump_time_step_units}\n" + text += f"{'Time step':>30}: {self.pump_time_step:.3f} {self.pump_time_step_units}\n" text += f"{'Hole':>30}: {self.hole}\n" return text @@ -485,6 +485,7 @@ def plot_time_profile(self, ax=None): ax.plot([time_left_FWHM, time_right_FWHM], [time_half_max, time_half_max], marker='o', color='tab:red', lw=3, label='FWHM') ax.set_xlabel('Time (fs)', fontsize=24) ax.set_ylabel('Time Gaussian', fontsize=24) + ax.set_title(f'Time profile (FWHM = {time_right_FWHM - time_left_FWHM:.4f} fs)') ax.legend() ax.grid() @@ -510,6 +511,7 @@ def plot_energy_profile(self, ax=None): ax.axvline(self.pump_energy, color='gray', lw=3, label='Pump energy') ax.set_xlabel('Energy (eV)', fontsize=24) ax.set_ylabel('Energy Gaussian', fontsize=24) + ax.set_title(f'Energy profile (FWHM = {energy_right_FWHM - energy_left_FWHM:.4f} eV)') ax.legend(loc='upper right') ax.grid() diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 26724e0..964f759 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -241,16 +241,6 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, total carrier number in the cnum_check folder (do not forget to link the epr file there). """ - print(f"{'PUMP PULSE PARAMETERS':*^70}") - print(f"{'Pump pulse energy (eV):':>40} {pump_energy:.3f}") - print(f"{'Pump pulse energy broadening (eV):':>40} {pump_spectral_width_fwhm:.4f}") - print(f"{'Finite width:':>40} {finite_width}") - if finite_width: - print(f"{'Pump pulse time step (fs):':>40} {pump_time_step:.3f}") - print(f"{'Pump pulse FWHM (fs):':>40} {pump_duration_fwhm:.3f}") - print(f"{'Pump pulse time window (fs):':>40} {pump_time_window:.3f}") - print("") - # Check electron and hole time steps elec_dyna_time_step = elec_dyna_run[1].time_step hole_dyna_time_step = hole_dyna_run[1].time_step diff --git a/src/perturbopy/postproc/utils/spectra_plots.py b/src/perturbopy/postproc/utils/spectra_plots.py index 12660b7..ec91871 100644 --- a/src/perturbopy/postproc/utils/spectra_plots.py +++ b/src/perturbopy/postproc/utils/spectra_plots.py @@ -101,7 +101,7 @@ def find_fwhm(x, y, num_interp_points=2000): def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, - h_occs, hole_kpoint_array, hole_energy_array, pump_energy, plot_scale=1.0): + h_occs, hole_kpoint_array, hole_energy_array, pump_energy, plot_scale=1e3): """ Plot occupation amplitude. Currently hardcoded to the section kz = 0. @@ -129,7 +129,7 @@ def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, Pump energy in eV. plot_scale : float - Scale factor for the scatter object sizes. Default is 1.0. + Scale factor for the scatter object sizes. """ # find where kz == 0 @@ -138,15 +138,17 @@ def plot_occ_ampl(e_occs, elec_kpoint_array, elec_energy_array, # Plot electron occupations idx = np.where((elec_kpoint_array[:, 1] == 0) & (elec_kpoint_array[:, 0] == 0)) + e_occs_max = np.max(e_occs[idx, :].ravel()) for i in range(np.size(elec_energy_array, axis=1)): ax.scatter(elec_kpoint_array[idx, 2], elec_energy_array[idx, i], s=0.5, c='black', alpha=0.5) - ax.scatter(elec_kpoint_array[idx, 2], elec_energy_array[idx, i], s=e_occs[idx, i][0] * 1e4 * plot_scale + 1e-10, c='red', alpha=0.5) + ax.scatter(elec_kpoint_array[idx, 2], elec_energy_array[idx, i], s=e_occs[idx, i][0] / e_occs_max * plot_scale + 1e-10, c='red', alpha=0.5) # Plot hole occupations idx = np.where((hole_kpoint_array[:, 1] == 0) & (hole_kpoint_array[:, 0] == 0)) + h_occs_max = np.max(h_occs[idx, :].ravel()) for i in range(np.size(hole_energy_array, axis=1)): ax.scatter(hole_kpoint_array[idx, 2], hole_energy_array[idx, i], s=0.5, c='black', alpha=0.5) - ax.scatter(hole_kpoint_array[idx, 2], hole_energy_array[idx, i], s=h_occs[idx, i][0] * 1e4 * plot_scale + 1e-10, c='red', + ax.scatter(hole_kpoint_array[idx, 2], hole_energy_array[idx, i], s=h_occs[idx, i][0] / h_occs_max * plot_scale + 1e-10, c='red', alpha=0.5) fsize = 16 @@ -161,7 +163,7 @@ def animate_pump_pulse(time_step, elec_delta_occs_array, elec_kpoint_array, elec_energy_array, hole_delta_occs_array, hole_kpoint_array, hole_energy_array, pump_energy, - plot_scale=1.0): + plot_scale=1e3): """ Animate the pump pulse excitation for electrons and holes. Defines fig and ax, initializes scatter objects for electron and hole occupations, and calls update_scatter. @@ -194,7 +196,7 @@ def animate_pump_pulse(time_step, Pump energy in eV, used only in title. plot_scale : float - Scale factor for the scatter object sizes. Default is 1.0. + Scale factor for the scatter object sizes. """ fig, ax = plt.subplots(1, 1, figsize=(9, 6)) @@ -235,7 +237,7 @@ def animate_pump_pulse(time_step, def update_scatter(anim_time, ax, time_step, idx_elec, idx_hole, elec_scat_list, hole_scat_list, - elec_delta_occs_array, hole_delta_occs_array, plot_scale=1.0): + elec_delta_occs_array, hole_delta_occs_array, plot_scale=1e3): """ Animate the pump pulse excitation for electrons and holes. @@ -270,17 +272,19 @@ def update_scatter(anim_time, ax, time_step, idx_elec, idx_hole, Array of hole occupation changes, similar to elec_delta_occs_array. plot_scale : float - Scale factor for the scatter object sizes. Default is 1.0. + Scale factor for the scatter object sizes. """ elec_num_bands = elec_delta_occs_array.shape[0] hole_num_bands = hole_delta_occs_array.shape[0] + e_occs_max = np.max(elec_delta_occs_array[:, :, idx_elec].ravel()) for i in range(elec_num_bands): - elec_scat_list[i].set_sizes(elec_delta_occs_array[i, anim_time, idx_elec].flatten() * 2e4 * plot_scale + 1e-10) + elec_scat_list[i].set_sizes(elec_delta_occs_array[i, anim_time, idx_elec].flatten() / e_occs_max * plot_scale + 1e-10) + h_occs_max = np.max(hole_delta_occs_array[:, :, idx_hole].ravel()) for i in range(hole_num_bands): - hole_scat_list[i].set_sizes(hole_delta_occs_array[i, anim_time, idx_hole].flatten() * 2e4 * plot_scale + 1e-10) + hole_scat_list[i].set_sizes(hole_delta_occs_array[i, anim_time, idx_hole].flatten() / h_occs_max * plot_scale + 1e-10) suffix = ax.get_title().split(';')[0] ax.set_title(f'{suffix}; Time: {anim_time * time_step:.2f} fs') @@ -378,7 +382,7 @@ def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, Pump energy in eV. pump_spectral_width_fwhm : float - Pump energy broadening in eV. + Pump energy broadening full width at half maximum (FWHM) in eV. scale : float Scale factor for the occupation amplitude. For plotting purposes. @@ -412,7 +416,7 @@ def plot_occs_on_bands(ax, kpath, bands, first_el_band_idx, hole_band = bands[i_hole_band] occ = gaussian(elec_band - hole_band, - pump_energy, pump_spectral_width_fwhm, hole_nband, elec_nband) + pump_energy, pump_spectral_width_fwhm) occs_amplitude_bands[i_elec_band - 1, :] += occ occs_amplitude_bands[i_hole_band - 1, :] += occ From ccc61a093162d2b4aaf6eb87cdbb35427bcf0207 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sat, 4 Jan 2025 13:20:47 +0100 Subject: [PATCH 12/21] pump pulse print format --- src/perturbopy/postproc/calc_modes/dyna_run.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index c88a3b1..ead37db 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -456,12 +456,12 @@ def __str__(self): """ text = 'Pump pulse excitation parameters:\n' - text += f"{'Pump energy':>30}: {self.pump_energy:.4f} {self.pump_energy_units}\n" - text += f"{'Pump duration FWHM':>30}: {self.pump_duration_fwhm:.2f} {self.pump_duration_fwhm_units}\n" - text += f"{'Energy broadening FWHM':>30}: {self.spectral_width_fwhm:.4f} {self.spectral_width_fwhm_units}\n" + text += f"{f'Pump energy ({self.pump_energy_units})':>30}: {self.pump_energy:.4f}\n" + text += f"{f'Energy broadening FWHM ({self.spectral_width_fwhm_units})':>30}: {self.spectral_width_fwhm:.4f}\n" + text += f"{f'Pulse duration FWHM ({self.pump_duration_fwhm_units})':>30}: {self.pump_duration_fwhm:.4f}\n" text += f"{'Pump factor':>30}: {self.pump_factor}\n" text += f"{'Number of steps':>30}: {self.num_steps}\n" - text += f"{'Time step':>30}: {self.pump_time_step:.3f} {self.pump_time_step_units}\n" + text += f"{f'Time step ({self.pump_time_step_units})':>30}: {self.pump_time_step}\n" text += f"{'Hole':>30}: {self.hole}\n" return text From a6be019c14f315baf3c90b03777b99fc1342d40d Mon Sep 17 00:00:00 2001 From: Ivan Maliyov Date: Sat, 4 Jan 2025 05:39:15 -0800 Subject: [PATCH 13/21] read_snaps option for DynaRun --- .../postproc/calc_modes/dyna_run.py | 40 ++++++++++++++----- 1 file changed, 31 insertions(+), 9 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index ead37db..b315351 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -19,6 +19,10 @@ class DynaRun(CalcMode): """ Class representation of a Perturbo dynamics-run calculation. + TODO: better handling with snaps. Make it as property. + Print a sensible message if snaps is None and access is attempted. + Implement a getter and setter for snaps linked to HDF5 file. + Attributes ---------- kpt : RecipPtDB @@ -40,19 +44,32 @@ class DynaRun(CalcMode): Python dictionary of DynaIndivRun objects containing results from each simulation """ - def __init__(self, cdyna_file, tet_file, pert_dict): + def __init__(self, cdyna_file, tet_file, pert_dict, read_snaps=True): """ Constructor method Parameters ---------- + cdyna_file : h5py.File + HDF5 file containing the results of the dynamics-run calculation. + The file is kept open after the object is created. + + tet_file : h5py.File + Tetrahedron HDF5 file obtained from the setup calculation. + pert_dict : dict - Dictionary containing the inputs and outputs from the dynamics-run calculation. + Dictionary containing the information from the dynamics-run YAML file. + + read_snaps : bool, optional + Flag to read the snapshots from the HDF5 file. + Reading the snapshots can be time-consuming and memory-intensive. """ self.timings = TimingGroup("dynamics-run") + self.read_snaps = read_snaps + super().__init__(pert_dict) if self.calc_mode != 'dynamics-run': @@ -102,10 +119,14 @@ def __init__(self, cdyna_file, tet_file, pert_dict): # a dynamics run must have at least one snap numk, numb = cdyna_file[dyn_str]['snap_t_1'][()].shape - snap_t = np.zeros((numb, numk, num_steps), dtype=np.float64) + if read_snaps: + print(f'Reading snapshots for {dyn_str}...') + snap_t = np.zeros((numb, numk, num_steps), dtype=np.float64) - for itime in range(num_steps): - snap_t[:, :, itime] = cdyna_file[dyn_str][f'snap_t_{itime + 1}'][()].T + for itime in range(num_steps): + snap_t[:, :, itime] = cdyna_file[dyn_str][f'snap_t_{itime + 1}'][()].T + else: + snap_t = None # Get E-field, which is only present if nonzero if "efield" in cdyna_file[dyn_str].keys(): @@ -116,7 +137,7 @@ def __init__(self, cdyna_file, tet_file, pert_dict): self._data[irun] = DynaIndivRun(num_steps, time_step, snap_t, time_units='fs', efield=efield) @classmethod - def from_hdf5_yaml(cls, cdyna_path, tet_path, yaml_path='pert_output.yml'): + def from_hdf5_yaml(cls, cdyna_path, tet_path, yaml_path='pert_output.yml', read_snaps=True): """ Class method to create a DynamicsRunCalcMode object from the HDF5 file and YAML file generated by a Perturbo calculation @@ -148,7 +169,7 @@ def from_hdf5_yaml(cls, cdyna_path, tet_path, yaml_path='pert_output.yml'): cdyna_file = open_hdf5(cdyna_path) tet_file = open_hdf5(tet_path) - return cls(cdyna_file, tet_file, yaml_dict) + return cls(cdyna_file, tet_file, yaml_dict, read_snaps=read_snaps) def close_hdf5_files(self): """ @@ -204,7 +225,7 @@ def __str__(self): Method to print the dynamics run information """ - title = f'dynamics-run: {self.prefix}' + title = f' dynamics-run: {self.prefix} ' text = f"{title:*^60}\n" text += f"This simulation has {self.num_runs} runs\n" @@ -213,6 +234,7 @@ def __str__(self): text += f"{'Number of steps':>30}: {dynamics_run.num_steps}\n" text += f"{'Time step (fs)':>30}: {dynamics_run.time_step}\n" text += f"{'Electric field (V/cm)':>30}: {dynamics_run.efield}\n\n" + text += f"{'Read snaps at init:':>30}: {self.read_snaps}\n" return text @@ -455,7 +477,7 @@ def __str__(self): Method to print the pump pulse excitation parameters. """ - text = 'Pump pulse excitation parameters:\n' + text = f'{" Pump pulse parameters ":*^60}\n' text += f"{f'Pump energy ({self.pump_energy_units})':>30}: {self.pump_energy:.4f}\n" text += f"{f'Energy broadening FWHM ({self.spectral_width_fwhm_units})':>30}: {self.spectral_width_fwhm:.4f}\n" text += f"{f'Pulse duration FWHM ({self.pump_duration_fwhm_units})':>30}: {self.pump_duration_fwhm:.4f}\n" From 56fa0c89b9893fd184d9e1078cc495ac47e92c10 Mon Sep 17 00:00:00 2001 From: Ivan Maliyov Date: Sat, 4 Jan 2025 06:59:18 -0800 Subject: [PATCH 14/21] print more info about DynaRun --- src/perturbopy/postproc/calc_modes/dyna_run.py | 14 +++++++++++--- src/perturbopy/postproc/utils/spectra_trans_abs.py | 4 ++-- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index b315351..06c1a65 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -227,14 +227,22 @@ def __str__(self): title = f' dynamics-run: {self.prefix} ' text = f"{title:*^60}\n" - text += f"This simulation has {self.num_runs} runs\n" + + text += f"{'Read snaps at init':>30}: {self.read_snaps}\n" + text += f"{'Number of k-points':>30}: {self.kpt.points_cart.shape[1]}\n" + nband = self._pert_dict['input parameters']['after conversion']['band_max'] - \ + self._pert_dict['input parameters']['after conversion']['band_min'] + 1 + text += f"{'Number of bands':>30}: {nband}\n" + text += f"{'Hole':>30}: {self._pert_dict['input parameters']['after conversion']['hole']}\n" + text += "\n" + + text += f"{'Number of runs':>30}: {self.num_runs}\n" for irun, dynamics_run in self._data.items(): text += f"{'Dynamics run':>30}: {irun}\n" text += f"{'Number of steps':>30}: {dynamics_run.num_steps}\n" text += f"{'Time step (fs)':>30}: {dynamics_run.time_step}\n" - text += f"{'Electric field (V/cm)':>30}: {dynamics_run.efield}\n\n" - text += f"{'Read snaps at init:':>30}: {self.read_snaps}\n" + text += f"{'Electric field (V/cm)':>30}: {dynamics_run.efield}\n" return text diff --git a/src/perturbopy/postproc/utils/spectra_trans_abs.py b/src/perturbopy/postproc/utils/spectra_trans_abs.py index 763847f..77b5c05 100644 --- a/src/perturbopy/postproc/utils/spectra_trans_abs.py +++ b/src/perturbopy/postproc/utils/spectra_trans_abs.py @@ -123,7 +123,7 @@ def compute_trans_abs(elec_dyna_run, # Energy grid for the transient absorption spectrum, affect the calculation time trans_abs_energy_grid = np.arange(bandgap, energy_max - energy_min, de_grid) num_energy_points = trans_abs_energy_grid.shape[0] - print(f"{'Number of energy points':<30}: {num_energy_points}\n") + print(f"{'Number of energy grid points':>30}: {num_energy_points}\n") # Find the interection between the electron and hole k-points print('Finding intersect k-points...') @@ -208,7 +208,7 @@ def compute_trans_abs(elec_dyna_run, print('Saving the data...') with trun.add('save data') as t: pump_energy = elec_dyna_run.pump_pulse.pump_energy - ending = f'_Epump_{pump_energy:.3f}' + ending = f'_Epump_{pump_energy:.4f}' # Total transient absorption np.save(f'trans_abs_dA{ending}', dA_elec + dA_hole) # Electron contribution From 1ff8ff1b02f32008eb22301d3148b73ccf12011e Mon Sep 17 00:00:00 2001 From: imaliyov Date: Wed, 8 Jan 2025 14:03:26 +0100 Subject: [PATCH 15/21] transient dipoles for pump pulse generation --- .../postproc/utils/spectra_generate_pulse.py | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 964f759..73f3104 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -175,7 +175,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, finite_width=True, animate=True, plot_scale=1.0, - cnum_check=True): + cnum_check=True, + tr_dipoles=None): """ Setup the Gaussian pump pulse excitation for electrons and holes. Write into the pump_pulse.h5 HDF5 file. @@ -239,6 +240,10 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, If True, create a cnum_check folder with a prefix_cdyna.h5 file with the total occupation summed over time steps. Run dynamics-pp with Perturbo to get the accurate total carrier number in the cnum_check folder (do not forget to link the epr file there). + + tr_dipoles: np.ndarray, optional + Transition dipoles squared, valence-to-conduction for each k-point and band. + Experimental feature. """ # Check electron and hole time steps @@ -286,6 +291,16 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, print(f"{'Band gap (Eg) (eV):':>40} {bandgap:.4f}") print("") + # Currently, the number of k-points of the transition dipoles + # must be the same as the number of k-points of the electron energy array + # In a future development, we will interpolate the transition dipoles on the fly + if tr_dipoles is not None: + elec_nbard_tr_dip, hole_nband_tr_dip, num_k_tr_dip = tr_dipoles.shape + + if elec_nbard_tr_dip != elec_nband or hole_nband_tr_dip != hole_nband \ + or num_k_tr_dip != elec_kpoint_array.shape[0]: + raise ValueError(f'Wrong shape of the transition dipoles array: {tr_dipoles.shape}') + elec_occs_amplitude = np.zeros_like(elec_energy_array) hole_occs_amplitude = np.zeros_like(hole_energy_array) @@ -317,6 +332,10 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, spectra_plots.gaussian(elec_energy_array[ekidx, iband] - hole_energy_array[hkidx, jband], pump_energy, pump_energy_broadening_sigma) + # Multiply by the transition dipoles squared if provided + if tr_dipoles is not None: + delta *= tr_dipoles[iband, jband, :] + # Only for the intersected k points, we add the delta occupation elec_occs_amplitude[ekidx, iband] += delta hole_occs_amplitude[hkidx, jband] += delta From 3494be476e28e6edb9db83ffa9d9cce30ece3bd6 Mon Sep 17 00:00:00 2001 From: Ivan Maliyov Date: Wed, 8 Jan 2025 03:50:49 -0800 Subject: [PATCH 16/21] energy_grid_max optional parameter for trans abs --- .../postproc/utils/spectra_trans_abs.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/src/perturbopy/postproc/utils/spectra_trans_abs.py b/src/perturbopy/postproc/utils/spectra_trans_abs.py index 77b5c05..f5f049a 100644 --- a/src/perturbopy/postproc/utils/spectra_trans_abs.py +++ b/src/perturbopy/postproc/utils/spectra_trans_abs.py @@ -4,6 +4,7 @@ """ import numpy as np +import warnings from .memory import get_size from .timing import TimingGroup @@ -30,6 +31,7 @@ def gaussian_delta(x, mu, sig): def compute_trans_abs(elec_dyna_run, hole_dyna_run, de_grid=0.02, + energy_grid_max=None, eta=0.02, save_npy=True, ): @@ -49,6 +51,10 @@ def compute_trans_abs(elec_dyna_run, de_grid : float Energy grid spacing for the transient absorption spectrum. Affects the calculation time. + energy_grid_max : float + Maximum probe energy for the spectrum. If None, the maximum band energy difference (condunction-valence) is used. + Minimum probe energy is always the bandgap. + eta : float Broadening parameter for the Gaussian delta functions applied to the energy grid. @@ -117,11 +123,18 @@ def compute_trans_abs(elec_dyna_run, CBM = np.min(elec_energy_array.ravel()) bandgap = CBM - VBM - energy_max = np.max(elec_energy_array.ravel()) - energy_min = np.min(hole_energy_array.ravel()) + energy_band_max = np.max(elec_energy_array.ravel()) + energy_band_min = np.min(hole_energy_array.ravel()) + + if energy_grid_max is None: + energy_grid_max = energy_band_max - energy_band_min + elif energy_grid_max < bandgap: + warnings.warn(f'The maximum probe energy specified ({energy_grid_max}) is smaller than the bandgap ({bandgap}).' + f' The maximum probe energy is set to max(cond)-min(val) {energy_band_max - energy_band_min}.') + energy_grid_max = energy_band_max - energy_band_min # Energy grid for the transient absorption spectrum, affect the calculation time - trans_abs_energy_grid = np.arange(bandgap, energy_max - energy_min, de_grid) + trans_abs_energy_grid = np.arange(bandgap, energy_grid_max, de_grid) num_energy_points = trans_abs_energy_grid.shape[0] print(f"{'Number of energy grid points':>30}: {num_energy_points}\n") From 4a08a0a950fba74508a7e83b82f589af1cf3e985 Mon Sep 17 00:00:00 2001 From: Ivan Maliyov Date: Sun, 4 May 2025 02:52:57 -0700 Subject: [PATCH 17/21] Include transition dipoles into the transient absorption calculation --- .../postproc/utils/spectra_generate_pulse.py | 15 +++++----- .../postproc/utils/spectra_trans_abs.py | 30 ++++++++++++++++--- 2 files changed, 34 insertions(+), 11 deletions(-) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 73f3104..dc304b8 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -176,7 +176,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, animate=True, plot_scale=1.0, cnum_check=True, - tr_dipoles=None): + tr_dipoles_sqr=None): """ Setup the Gaussian pump pulse excitation for electrons and holes. Write into the pump_pulse.h5 HDF5 file. @@ -241,8 +241,9 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, occupation summed over time steps. Run dynamics-pp with Perturbo to get the accurate total carrier number in the cnum_check folder (do not forget to link the epr file there). - tr_dipoles: np.ndarray, optional + tr_dipoles_sqr: np.ndarray, optional Transition dipoles squared, valence-to-conduction for each k-point and band. + Currently, the k-grid for dipoles must match the one for electrons. Experimental feature. """ @@ -294,12 +295,12 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, # Currently, the number of k-points of the transition dipoles # must be the same as the number of k-points of the electron energy array # In a future development, we will interpolate the transition dipoles on the fly - if tr_dipoles is not None: - elec_nbard_tr_dip, hole_nband_tr_dip, num_k_tr_dip = tr_dipoles.shape + if tr_dipoles_sqr is not None: + elec_nbard_tr_dip, hole_nband_tr_dip, num_k_tr_dip = tr_dipoles_sqr.shape if elec_nbard_tr_dip != elec_nband or hole_nband_tr_dip != hole_nband \ or num_k_tr_dip != elec_kpoint_array.shape[0]: - raise ValueError(f'Wrong shape of the transition dipoles array: {tr_dipoles.shape}') + raise ValueError(f'Wrong shape of the transition dipoles array: {tr_dipoles_sqr.shape}') elec_occs_amplitude = np.zeros_like(elec_energy_array) hole_occs_amplitude = np.zeros_like(hole_energy_array) @@ -333,8 +334,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_energy, pump_energy_broadening_sigma) # Multiply by the transition dipoles squared if provided - if tr_dipoles is not None: - delta *= tr_dipoles[iband, jband, :] + if tr_dipoles_sqr is not None: + delta *= tr_dipoles_sqr[iband, jband, :] # Only for the intersected k points, we add the delta occupation elec_occs_amplitude[ekidx, iband] += delta diff --git a/src/perturbopy/postproc/utils/spectra_trans_abs.py b/src/perturbopy/postproc/utils/spectra_trans_abs.py index f5f049a..699ec04 100644 --- a/src/perturbopy/postproc/utils/spectra_trans_abs.py +++ b/src/perturbopy/postproc/utils/spectra_trans_abs.py @@ -34,7 +34,7 @@ def compute_trans_abs(elec_dyna_run, energy_grid_max=None, eta=0.02, save_npy=True, - ): + tr_dipoles_sqr=None): """ Compute the transient absorption spectrum from the electron and hole dynamics simulations. The data is saved in the current directory as numpy binary files (.npy). @@ -75,6 +75,11 @@ def compute_trans_abs(elec_dyna_run, dA_hole : np.ndarray Hole contribution to the transient absorption spectrum. Shape (num_energy_points, num_steps). + + tr_dipoles_sqr : np.ndarray + Transition dipoles squared, valence-to-conduction for each k-point and band. + Currently, the k-grid for dipoles must match the one for electrons. + Experimental feature. """ trun = TimingGroup('trans. abs.') @@ -83,7 +88,7 @@ def compute_trans_abs(elec_dyna_run, # Check the consistency of the two DynaRun objects # 1. Check the k-grids. q-grids can be different - if not np.alltrue(elec_dyna_run.boltz_kdim == hole_dyna_run.boltz_kdim): + if not np.all(elec_dyna_run.boltz_kdim == hole_dyna_run.boltz_kdim): raise ValueError('The k-grids of the electron and hole simulations are different.') # 2. Check the simulation times @@ -150,6 +155,19 @@ def compute_trans_abs(elec_dyna_run, num_intersect_kpoints = ekidx.size + # Currently, the number of k-points of the transition dipoles + # must be the same as the number of k-points of the electron energy array + # In a future development, we will interpolate the transition dipoles on the fly + if tr_dipoles_sqr is not None: + elec_nband_tr_dip, hole_nband_tr_dip, num_k_tr_dip = tr_dipoles_sqr.shape + + if elec_nband_tr_dip != elec_nband or hole_nband_tr_dip != hole_nband \ + or num_k_tr_dip != elec_kpoint_array.shape[0]: + raise ValueError(f'Wrong shape of the transition dipoles array: {tr_dipoles_sqr.shape}') + else: + print('\nTransition dipoles are provided.') + print(f"{'Tr. dip. shape (num. elec. bands, num. hole bands, num. k points)':>30}: {tr_dipoles_sqr.shape}\n") + # Compute all the Gaussian deltas on the intersect grid print('Computing Gaussian deltas...') with trun.add('compute deltas') as t: @@ -161,12 +179,16 @@ def compute_trans_abs(elec_dyna_run, for jband in range(hole_nband): for ienergy in range(num_energy_points): # Factor of two from spin. - # TODO: add transient dipoles here - ALL_DELTAS[iband, jband, :, ienergy] = 2.0 * \ + ALL_DELTAS[iband, jband, :, ienergy] = \ gaussian_delta(elec_energy_array[ekidx, iband] - hole_energy_array[hkidx, jband], trans_abs_energy_grid[ienergy], eta) + if tr_dipoles_sqr is not None: + ALL_DELTAS[iband, jband, :, ienergy] *= 0.5 * tr_dipoles_sqr[iband, jband, :] + else: + ALL_DELTAS[iband, jband, :, ienergy] *= 2.0 + get_size(ALL_DELTAS, 'ALL_DELTAS', dump=True) print('') From 7b540b1af5f210adb01604a7a25ce5053c263ef9 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sat, 3 May 2025 15:35:59 +0200 Subject: [PATCH 18/21] fwhm_from_sigma function --- .../postproc/utils/spectra_generate_pulse.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index dc304b8..a50577e 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -34,6 +34,16 @@ def sigma_from_fwhm(fwhm): return sigma +def fwhm_from_sigma(sigma): + """ + Comupte FWHM from Gaussian sigma. + """ + + fwhm = 2.0 * np.sqrt(2.0 * np.log(2.0)) * sigma + + return fwhm + + def delta_occs_pulse_coef(t, dt, tw, sigma): """ Additional occupation due to the pulse excitation. From d372b54296292f85f6baa05d5440ec7658a13068 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Thu, 8 May 2025 20:32:34 +0200 Subject: [PATCH 19/21] Update for spectroscopy scripts: generation and postprocessing. * fixed bug with the k-grid indexing for transition dipoles * fixed bug with writing the cdyna HDF5 file for the carrier concentration check (cnum_check flag) * more strict requirements for the .to_cdyna method on the array of snapshots shape --- .../postproc/calc_modes/dyna_run.py | 15 +++++++++-- .../postproc/utils/spectra_generate_pulse.py | 25 ++++++++++++++----- .../postproc/utils/spectra_trans_abs.py | 2 +- 3 files changed, 33 insertions(+), 9 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 06c1a65..8864585 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -302,7 +302,7 @@ def to_cdyna_h5(prefix, band_structure_ryd, snap_array, time_step_fs, Band structure in Rydberg units. Shape: num_kpoints x num_bands snap_array : np.ndarray - Array of carrier occupations. Shape: (num_steps, (num_kpoints x num_bands)) + Array of carrier occupations. Shape: (num_steps, num_kpoints, num_bands) time_step_fs : float or np.ndarray Time step in fs. If float, the same time step is used for all runs. @@ -345,6 +345,17 @@ def to_cdyna_h5(prefix, band_structure_ryd, snap_array, time_step_fs, else: raise FileExistsError(f'File {new_cdyna_filename} already exists. Set overwrite=True to overwrite it.') + if snap_array.ndim != 3: + raise ValueError("snap_array must be 3D. Shape: (num_steps, num_kpoints, num_bands)") + else: + num_steps_read, num_kpoints_read, num_bands_read = snap_array.shape + print("Occupations shape to write:") + print(f"num_steps: {num_steps_read}, num_kpoints: {num_kpoints_read}, num_bands: {num_bands_read}") + + if num_kpoints_read < num_steps_read or num_kpoints_read < num_bands_read: + warnings.warn("The number of k-points of snap array to write is smaller " + "than the number of time steps or bands.") + new_cdyna_file = open_hdf5(new_cdyna_filename, 'w') new_cdyna_file.create_dataset('band_structure_ryd', data=band_structure_ryd) @@ -374,7 +385,7 @@ def to_cdyna_h5(prefix, band_structure_ryd, snap_array, time_step_fs, new_cdyna_file[dyn_str].create_dataset('time_step_fs', data=time_step_fs[irun - 1]) for itime in range(snap_array.shape[0]): - new_cdyna_file[dyn_str].create_dataset(f'snap_t_{itime + offset}', data=snap_array[itime, :, np.newaxis]) + new_cdyna_file[dyn_str].create_dataset(f'snap_t_{itime + offset}', data=snap_array[itime, :, :]) close_hdf5(new_cdyna_file) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index a50577e..7bcd5f9 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -248,8 +248,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, cnum_check : bool, optional If True, create a cnum_check folder with a prefix_cdyna.h5 file with the total - occupation summed over time steps. Run dynamics-pp with Perturbo to get the accurate - total carrier number in the cnum_check folder (do not forget to link the epr file there). + occupation summed over time steps. Run dynamics-pp with Perturbo (without dynamics-run) + to get the accurate total carrier number in the cnum_check folder (link the epr file there). tr_dipoles_sqr: np.ndarray, optional Transition dipoles squared, valence-to-conduction for each k-point and band. @@ -310,7 +310,8 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, if elec_nbard_tr_dip != elec_nband or hole_nband_tr_dip != hole_nband \ or num_k_tr_dip != elec_kpoint_array.shape[0]: - raise ValueError(f'Wrong shape of the transition dipoles array: {tr_dipoles_sqr.shape}') + raise ValueError(f'Wrong shape of the transition dipoles array: {tr_dipoles_sqr.shape}\n' + f'Expected shape: ({elec_nband}, {hole_nband}, {elec_kpoint_array.shape[0]})') elec_occs_amplitude = np.zeros_like(elec_energy_array) hole_occs_amplitude = np.zeros_like(hole_energy_array) @@ -344,8 +345,10 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, pump_energy, pump_energy_broadening_sigma) # Multiply by the transition dipoles squared if provided + # Transition dipoles are assumed to be on the same k-point grid + # as the electron energy array if tr_dipoles_sqr is not None: - delta *= tr_dipoles_sqr[iband, jband, :] + delta *= tr_dipoles_sqr[iband, jband, ekidx] # Only for the intersected k points, we add the delta occupation elec_occs_amplitude[ekidx, iband] += delta @@ -497,16 +500,26 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, elec_cnum_check_path = os.path.join(os.path.dirname(elec_pump_pulse_path), 'cnum_check') hole_cnum_check_path = os.path.join(os.path.dirname(hole_pump_pulse_path), 'cnum_check') + + # Format the occupations for cdyna.h5 writing + total_occ_elec_write = total_occ_elec.copy() + # Transpose: num_band, num_k -> num_k, num_band + total_occ_elec_write = total_occ_elec_write.T + # Add the dummy time axis: 1, num_k, num_band + total_occ_elec_write = np.expand_dims(total_occ_elec_write, axis=0) elec_dyna_run.to_cdyna_h5(elec_dyna_run.prefix, elec_dyna_run._energies, - total_occ_elec, + total_occ_elec_write, elec_dyna_run[1].time_step, path=elec_cnum_check_path, overwrite=True) + total_occ_hole_write = total_occ_hole.copy() + total_occ_hole_write = total_occ_hole_write.T + total_occ_hole_write = np.expand_dims(total_occ_hole_write, axis=0) hole_dyna_run.to_cdyna_h5(hole_dyna_run.prefix, hole_dyna_run._energies, - total_occ_hole, + total_occ_hole_write, hole_dyna_run[1].time_step, path=hole_cnum_check_path, overwrite=True) diff --git a/src/perturbopy/postproc/utils/spectra_trans_abs.py b/src/perturbopy/postproc/utils/spectra_trans_abs.py index 699ec04..7ec7192 100644 --- a/src/perturbopy/postproc/utils/spectra_trans_abs.py +++ b/src/perturbopy/postproc/utils/spectra_trans_abs.py @@ -185,7 +185,7 @@ def compute_trans_abs(elec_dyna_run, trans_abs_energy_grid[ienergy], eta) if tr_dipoles_sqr is not None: - ALL_DELTAS[iband, jband, :, ienergy] *= 0.5 * tr_dipoles_sqr[iband, jband, :] + ALL_DELTAS[iband, jband, :, ienergy] *= 0.5 * tr_dipoles_sqr[iband, jband, ekidx] else: ALL_DELTAS[iband, jband, :, ienergy] *= 2.0 From 4450390512d18d29d252a4b47f0da78a938bbfd0 Mon Sep 17 00:00:00 2001 From: imaliyov Date: Sun, 11 May 2025 15:50:32 +0200 Subject: [PATCH 20/21] Corrected the total occupation for holes for cnum_check --- src/perturbopy/postproc/utils/spectra_generate_pulse.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/perturbopy/postproc/utils/spectra_generate_pulse.py b/src/perturbopy/postproc/utils/spectra_generate_pulse.py index 7bcd5f9..265d6c0 100644 --- a/src/perturbopy/postproc/utils/spectra_generate_pulse.py +++ b/src/perturbopy/postproc/utils/spectra_generate_pulse.py @@ -517,6 +517,7 @@ def setup_pump_pulse(elec_pump_pulse_path, hole_pump_pulse_path, total_occ_hole_write = total_occ_hole.copy() total_occ_hole_write = total_occ_hole_write.T total_occ_hole_write = np.expand_dims(total_occ_hole_write, axis=0) + total_occ_hole_write = 1.0 - total_occ_hole_write hole_dyna_run.to_cdyna_h5(hole_dyna_run.prefix, hole_dyna_run._energies, total_occ_hole_write, From 617f3124bc68ee6bc8145e5e09afa8e4217f9879 Mon Sep 17 00:00:00 2001 From: Ivan Maliyov Date: Sun, 11 May 2025 08:25:01 -0700 Subject: [PATCH 21/21] Format output data of DynaRun --- src/perturbopy/postproc/calc_modes/dyna_run.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/perturbopy/postproc/calc_modes/dyna_run.py b/src/perturbopy/postproc/calc_modes/dyna_run.py index 8864585..06e55d2 100644 --- a/src/perturbopy/postproc/calc_modes/dyna_run.py +++ b/src/perturbopy/postproc/calc_modes/dyna_run.py @@ -496,11 +496,13 @@ def __str__(self): Method to print the pump pulse excitation parameters. """ - text = f'{" Pump pulse parameters ":*^60}\n' + carrier = 'hole' if self.hole else 'electron' + title = f'Pump pulse parameters ({carrier})' + text = f'{title:*^60}\n' text += f"{f'Pump energy ({self.pump_energy_units})':>30}: {self.pump_energy:.4f}\n" text += f"{f'Energy broadening FWHM ({self.spectral_width_fwhm_units})':>30}: {self.spectral_width_fwhm:.4f}\n" text += f"{f'Pulse duration FWHM ({self.pump_duration_fwhm_units})':>30}: {self.pump_duration_fwhm:.4f}\n" - text += f"{'Pump factor':>30}: {self.pump_factor}\n" + text += f"{'Pump factor':>30}: {self.pump_factor:.3e}\n" text += f"{'Number of steps':>30}: {self.num_steps}\n" text += f"{f'Time step ({self.pump_time_step_units})':>30}: {self.pump_time_step}\n" text += f"{'Hole':>30}: {self.hole}\n"