From 3b45609cd2eb2be58cd26408da8245fb6771372b Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Thu, 11 Jun 2026 14:49:26 +0800 Subject: [PATCH 1/3] Update new script for RT-TDDFT absorption spectrum --- tools/rt-tddft-tools/dipole.py | 171 ------ tools/rt-tddft-tools/plot_absorption.py | 684 ++++++++++++++---------- 2 files changed, 406 insertions(+), 449 deletions(-) delete mode 100644 tools/rt-tddft-tools/dipole.py diff --git a/tools/rt-tddft-tools/dipole.py b/tools/rt-tddft-tools/dipole.py deleted file mode 100644 index 91fcbdb53f8..00000000000 --- a/tools/rt-tddft-tools/dipole.py +++ /dev/null @@ -1,171 +0,0 @@ -import numpy as np -from os import PathLike -from matplotlib.figure import Figure -from matplotlib import axes -import scipy.constants as sc - -Freq2eV = sc.h/sc.eV*1e15 # 1/fs to eV - - -class Dipole: - """Parse Dipole Data""" - - def __init__(self, dipolefile: PathLike, dt: float) -> None: - self.dipolefile = dipolefile - self._indices, self.dipole = self.read(self.dipolefile) - self.dt = dt # in fs - - @classmethod - def read(cls, filename: PathLike): - """Read dipole data file - - :params filename: string of dipole data file - """ - data = np.loadtxt(filename, dtype=float) - indices = data[:, 0] - dipole = data[:, 1:].transpose() - for i in range(3): - dipole[i,:] = dipole[i,:] - dipole[i,0] - - return indices, dipole - - @property - def dipole_data(self): - return self.dipole - - def plot_dipole(self, fig: Figure, ax: axes.Axes, directions: list = [0, 1, 2], time_range: list = []): - """Plot dipole data in x,y,z directions - - :params fig: (matplotlib.figure.Figure) - :params ax: (matplotlib.axes.Axes) - :params directions: (list) 0->X, 1->Y, 2->Z - :params time_range: (list) range of time (in unit fs) to plot - """ - - labels = {0: 'X', 1: 'Y', 2: 'Z'} - - for direc in directions: - ax.plot(self._indices*self.dt, self.dipole_data[direc], label=labels[direc]) - - ax.set_xlabel('Times (fs)') - ax.set_ylabel('Dipole') - ax.legend() - if time_range: - ax.set_xlim(time_range) - - return fig, ax - -class Efield: - """Parse Efield Data""" - - def __init__(self, Efile: list[list[PathLike]], dt: float, nstep:int) -> None: - self.Efile = Efile - self.dt=dt - self.efield=self.reade(self.Efile, nstep) - - @classmethod - def reade(cls, Efile: list[list[PathLike]], nstep: int): - """Read dipole data file - - :params Efile: string type 2D list of Efield data file, Efile[i][j] is the jth Efield data file in ith direction - :params nstep: number of steps in simulation - """ - Efield = np.zeros([3,nstep]) - for i in range(len(Efile)): - for file in Efile[i]: - Edata = np.loadtxt(file, dtype=float) - Efield[i,0:(len(Edata))] += Edata[:, 1] - return Efield - - @property - def efield_data(self): - return self.efield - -class Absorption(Dipole, Efield): - """Calculate Absorption Spectrum under light field""" - def __init__(self, dipolefile: PathLike,Efile: list[list[PathLike]], dt: float) -> None: - Dipole.__init__(self, dipolefile, dt) - Efield.__init__(self, Efile, dt, len(self._indices)) - - def padding(self, data: np.ndarray): - """Zero padding for FFT - - :params data: (np.ndarray) data to be padded - """ - #mask part - Ndim=len(self._indices)*10 - index=np.linspace(0,Ndim-1,Ndim) - t=self._indices*self.dt - b=5 - mask=np.exp(-b*t/t[-1]) - #padding part - padding=np.zeros(Ndim-len(data)) - data_pass=np.concatenate((data*mask, padding)) - return data_pass - - def alpha(self,dirc: int): - """Calculate alpha - - :params dirc: (int) 0->X, 1->Y, 2->Z - """ - #FFT part - dipole=self.padding(self.dipole_data[dirc]) - efield=self.padding(self.efield_data[dirc]) - dipole_fft = np.fft.fft(dipole) - efield_fft = np.fft.fft(efield) - alpha = np.abs((dipole_fft/efield_fft).imag) - return alpha - - def plot_abs(self, fig: Figure, ax: axes.Axes, directions: list = [0, 1, 2], x_range: list = [], unit: str = 'eV'): - """Plot Absportion Spectrum under Delta light field in x,y,z directions - - :params fig: (matplotlib.figure.Figure) - :params ax: (matplotlib.axes.Axes) - :params directions: (list) 0->X, 1->Y, 2->Z - :params x_range: (list) range of energies (in unit eV) to plot - :params unit: (str) - """ - - assert unit in ['eV', 'nm'] - labels = {0: 'X', 1: 'Y', 2: 'Z'} - Ndim=len(self._indices)*10 - index=np.linspace(0,Ndim-1,Ndim) - energies = Freq2eV*index/self.dt/len(index) - x_data = energies if unit == 'eV' else sc.nu2lambda( - sc.eV/sc.h*energies)*1e9 - - #plot the adsorption spectra and output the data - adsorption_spectra_data = x_data[:, np.newaxis] - for direc in directions: - alpha = self.alpha(direc) - ax.plot(x_data, alpha, label=labels[direc]) - adsorption_spectra_data = np.concatenate((adsorption_spectra_data, alpha[:, np.newaxis]),axis=1) - np.savetxt('absorpation_spectra.dat', adsorption_spectra_data) - - xlabel = 'Energy (eV)' if unit == 'eV' else 'Wave Length (nm)' - ax.set_xlabel(xlabel) - ax.set_ylabel('Absportion') - ax.legend() - if x_range: - ax.set_xlim(x_range) - lim_range=index[(x_data>x_range[0])&(x_data eV +# Conversion: V/Å -> a.u. V_ANGSTROM_TO_AU_EFIELD = (sc.e / (4 * sc.pi * sc.epsilon_0 * sc.physical_constants['Bohr radius'][0]**2) * 1e-10) -# Conversion: fs → a.u. of time, 1 fs = 41.34137 a.u. AU_TIME = sc.physical_constants['atomic unit of time'][0] -FS_TO_AU_TIME = 1e-15 / AU_TIME +FS_TO_AU_TIME = 1e-15 / AU_TIME # Conversion: fs -> a.u. of time + # ====================== -# Utility Functions +# Data Loading Utilities # ====================== -def validate_inputs(direc: List[int], efield_path: List[str]) -> None: - if len(direc) != len(efield_path): - raise ValueError(f"Number of --direc ({len(direc)}) != --efield_path ({len(efield_path)})") - if not all(d in (0, 1, 2) for d in direc): - raise ValueError("--direc values must be 0 (x), 1 (y), or 2 (z)") - -def group_files_by_direction(direc: List[int], paths: List[str]) -> List[List[str]]: - """Return [[x_files], [y_files], [z_files]]""" - groups = [[], [], []] - for d, p in zip(direc, paths): - groups[d].append(p) - return groups - -def zero_pad_and_mask(signal: np.ndarray, dt: float, pad_factor: int = 10, decay_beta: float = 5.0) -> np.ndarray: - """Apply exponential decay mask and zero-padding for FFT.""" - n = len(signal) - t = np.arange(n) * dt - mask = np.exp(-decay_beta * t / t[-1]) if t[-1] != 0 else np.ones_like(t) - padded = np.zeros(n * pad_factor) - padded[:n] = signal * mask - return padded +def load_signal(filename: str, step_start: int, step_end: int, remove_dc: bool) -> np.ndarray: + """Load polarization response data (dipole/current). Strict checking is enforced.""" + try: + data = np.loadtxt(filename) + except OSError as e: + raise RuntimeError(f"Error loading signal file {filename}: {e}") -def energy_from_fft(N: int, dt: float) -> np.ndarray: - """Return positive-frequency energy array in eV.""" - freqs = np.fft.fftfreq(N, d=dt)[:N//2] # in 1/fs - return freqs * FREQ2EV + available_rows = data.shape[0] -# ====================== -# Data Loaders -# ====================== + # Enforce strict validation: signal file length must cover the requested step range + if available_rows < step_end: + raise ValueError( + f"Execution halted: Signal file '{filename}' contains only {available_rows} rows, " + f"which is less than the requested --step_end ({step_end}). " + f"Please reduce --step_end or ensure the simulation completed fully." + ) -def load_dipole(filename: str, step_start: int, step_end: int, remove_dc: bool) -> np.ndarray: - """Load current/dipole. Optionally remove DC component.""" - data = np.loadtxt(filename) - end = min(step_end or len(data), data.shape[0]) - sliced = data[step_start:end, 1:].T # shape: (3, N) + # Slice data within the targeted window and extract direction components + sliced = data[step_start:step_end, 1:].T # shape: (3, n_steps) if remove_dc: - # Remove DC (zero-frequency) component by subtracting mean for i in range(3): sliced[i] -= np.mean(sliced[i]) else: - # Subtract initial value for i in range(3): sliced[i] -= sliced[i, 0] return sliced + def load_efield(file_groups: List[List[str]], step_start: int, step_end: int, remove_dc: bool) -> np.ndarray: - """Load and sum E-field files per direction. Convert V/Å → a.u. Optionally remove DC component.""" + """Load and sum E-field profiles. Elastic zero-padding is applied at the tail if mismatched.""" n_steps = step_end - step_start efield = np.zeros((3, n_steps)) + for i, files in enumerate(file_groups): for f in files: - arr = np.loadtxt(f) - efield[i] += arr[step_start:step_end, 1] / V_ANGSTROM_TO_AU_EFIELD # V/Å → a.u. - # Remove DC (zero-frequency) component by subtracting mean + try: + arr = np.loadtxt(f) + except OSError as e: + raise RuntimeError(f"Error loading efield file {f}: {e}") + + available_rows = arr.shape[0] + start_idx = min(step_start, available_rows) + end_idx = min(step_end, available_rows) + + sliced = arr[start_idx:end_idx, 1] + + # Elastic padding mechanism: append trailing zeros if row count is slightly less than step_end + if len(sliced) < n_steps: + padded = np.zeros(n_steps) + padded[:len(sliced)] = sliced + sliced = padded + + efield[i] += sliced / V_ANGSTROM_TO_AU_EFIELD + if remove_dc and n_steps > 0: efield[i] -= np.mean(efield[i]) return efield + +def validate_inputs(direc: List[int], efield_path: List[str]) -> None: + if len(direc) != len(efield_path): + raise ValueError(f"Input mismatch: Number of --direc ({len(direc)}) != --efield_path ({len(efield_path)})") + if not all(d in (0, 1, 2) for d in direc): + raise ValueError("Invalid direction: --direc values must be chosen from 0 (x), 1 (y), or 2 (z)") + + +def group_files_by_direction(direc: List[int], paths: List[str]) -> List[List[str]]: + groups = [[], [], []] + for d, p in zip(direc, paths): + groups[d].append(p) + return groups + + +def zero_pad_and_mask(signal: np.ndarray, dt: float, pad_factor: int = 10, decay_beta: float = 5.0) -> np.ndarray: + n = len(signal) + t = np.arange(n) * dt + mask = np.exp(-decay_beta * t / t[-1]) if t[-1] != 0 else np.ones_like(t) + padded = np.zeros(n * pad_factor) + padded[:n] = signal * mask + return padded + + # ====================== -# Core Calculator +# Core Calculator Class # ====================== class AbsorptionSpectrum: - def __init__( - self, - dipole: np.ndarray, - efield: np.ndarray, - dt: float, - system_type: str, - pad_factor: int = 10, - decay_beta: float = 5.0, - ): + def __init__(self, dipole: np.ndarray, efield: np.ndarray, dt: float, + response_source: str, spectrum_type: str, pad_factor: int = 10, + decay_beta: float = 5.0, volume: float = 1.0): self.dipole = dipole self.efield = efield self.dt = dt self.dt_au = self.dt * FS_TO_AU_TIME - self.system_type = system_type + self.response_source = response_source + self.spectrum_type = spectrum_type self.pad_factor = pad_factor self.decay_beta = decay_beta + self.volume = volume self.N_orig = dipole.shape[1] self.N_pad = self.N_orig * pad_factor - def compute_alpha(self, direction: int) -> np.ndarray: - """Compute alpha(ω) for a given direction (0,1,2). Returns full FFT array.""" + def compute_complex_alpha(self, direction: int) -> Tuple[np.ndarray, np.ndarray]: + """Compute the complex response function.""" d_pad = zero_pad_and_mask(self.dipole[direction], self.dt_au, self.pad_factor, self.decay_beta) e_pad = zero_pad_and_mask(self.efield[direction], self.dt_au, self.pad_factor, self.decay_beta) - d_fft = np.fft.fft(d_pad) - e_fft = np.fft.fft(e_pad) - - ratio = np.divide(d_fft, e_fft, out=np.zeros_like(d_fft), where=e_fft != 0) - - if self.system_type == "dipole_epsilon": - return np.abs(ratio.imag) - elif self.system_type == "dipole_sigma": - # Full frequency array (length = N_pad) - freqs_full = np.fft.fftfreq(self.N_pad, d=self.dt_au) # f in a.u.^{-1} - omega = 2 * np.pi * freqs_full # ω in a.u.^{-1} - alpha = np.multiply(ratio.imag, omega) - return np.abs(alpha) - elif self.system_type == "current_epsilon": - # Full frequency array (length = N_pad) - freqs_full = np.fft.fftfreq(self.N_pad, d=self.dt_au) # f in a.u.^{-1} - omega = 2 * np.pi * freqs_full # ω in a.u.^{-1} - alpha = np.divide(ratio.real, omega, out=np.zeros_like(ratio.real), where=omega != 0) - return np.abs(alpha) - elif self.system_type == "current_sigma": - return np.abs(ratio.real) - else: - raise ValueError("system_type must be 'dipole_epsilon', 'dipole_sigma', 'current_epsilon', or 'current_sigma'") - - def get_positive_alpha(self, direction: int) -> np.ndarray: - """Return alpha for positive frequencies only.""" - return self.compute_alpha(direction)[:self.N_pad // 2] + # Convert FFT output to physics convention (e^{iwt}) by applying complex conjugate + d_fft = np.conj(np.fft.fft(d_pad)) + e_fft = np.conj(np.fft.fft(e_pad)) + + # Safe complex division using boolean mask to completely bypass complex ufunc mask bugs + ratio = np.zeros_like(d_fft, dtype=complex) + non_zero_mask = (e_fft != 0) + ratio[non_zero_mask] = d_fft[non_zero_mask] / e_fft[non_zero_mask] + + freqs_full = np.fft.fftfreq(self.N_pad, d=self.dt_au) + omega = 2 * np.pi * freqs_full + + if self.response_source == "dipole": + if self.spectrum_type == "epsilon": + alpha_real = 1 + (4 * np.pi * ratio.real) / self.volume + alpha_imag = (4 * np.pi * ratio.imag) / self.volume + return alpha_real, alpha_imag + elif self.spectrum_type == "sigma": + alpha_real = (omega * ratio.imag) / self.volume + alpha_imag = (-omega * ratio.real) / self.volume + return alpha_real, alpha_imag + + elif self.response_source == "current": + if self.spectrum_type == "sigma": + return ratio.real, ratio.imag + elif self.spectrum_type == "epsilon": + alpha_real = np.zeros_like(ratio.real) + alpha_imag = np.zeros_like(ratio.real) + np.divide(4 * np.pi * ratio.imag, omega, out=alpha_real, where=omega != 0) + alpha_real = 1 - alpha_real + np.divide(4 * np.pi * ratio.real, omega, out=alpha_imag, where=omega != 0) + return alpha_real, alpha_imag + + raise ValueError("Invalid setup: Check your response_source and spectrum_type parameters.") + + def get_positive_alpha(self, direction: int) -> Tuple[np.ndarray, np.ndarray]: + alpha_real, alpha_imag = self.compute_complex_alpha(direction) + half = self.N_pad // 2 + return alpha_real[:half], alpha_imag[:half] def get_energy_axis(self) -> np.ndarray: - return energy_from_fft(self.N_pad, self.dt) + freqs = np.fft.fftfreq(self.N_pad, d=self.dt)[:self.N_pad//2] + return freqs * FREQ2EV + + +# ========================== +# Specialized Plotting Tools +# ========================== + +def calculate_y_limits(x_vals: np.ndarray, y_arrays: List[np.ndarray], x_range: Tuple[float, float], + percentile: float = 99.5, ignore_zero_energy: float = 0.2) -> Optional[Tuple[float, float]]: + """ + Calculate automatic Y-axis limits using percentile filtering. + A low-energy threshold is utilized to bypass divergences or numerical artifacts near zero frequency. + """ + safe_x_min = max(x_range[0], ignore_zero_energy) + mask = (x_vals >= safe_x_min) & (x_vals <= x_range[1]) + + # Fallback to the original requested range if the exclusion window covers the entire domain + if not np.any(mask): + mask = (x_vals >= x_range[0]) & (x_vals <= x_range[1]) + if not np.any(mask): + return (0, 1) + + valid_y = np.concatenate([arr[mask] for arr in y_arrays]) + valid_y = valid_y[np.isfinite(valid_y)] + + if len(valid_y) == 0: + return (0, 1) + + y_low = np.percentile(valid_y, 100 - percentile) + y_high = np.percentile(valid_y, percentile) + margin = 0.1 * (y_high - y_low) if y_high != y_low else 0.1 * max(abs(y_high), 1.0) + + return y_low - margin, y_high + margin -# ====================== -# Plotting Utilities -# ====================== def plot_time_series(ax, time: np.ndarray, data: np.ndarray, directions: List[int], labels: List[str], ylabel: str): for d in directions: @@ -170,275 +261,312 @@ def plot_time_series(ax, time: np.ndarray, data: np.ndarray, directions: List[in ax.set_xlabel("Time (fs)") ax.set_ylabel(ylabel) ax.legend() + ax.grid(True, linestyle=':', alpha=0.6) -def plot_absorption(ax, x_vals: np.ndarray, alphas: List[np.ndarray], directions: List[int], xlabel: str, x_range: Optional[Tuple] = None): - labels = {0: "$x$", 1: "$y$", 2: "$z$"} - for d, alpha in zip(directions, alphas): - ax.plot(x_vals, alpha, label=labels[d]) - ax.set_xlabel(xlabel) - ax.set_ylabel("Intensity") - ax.legend() - if x_range: - ax.set_xlim(x_range) - mask = (x_vals >= min(x_range)) & (x_vals <= max(x_range)) - if np.any(mask): - ax.set_ylim(0, 1.2 * max(max(a[mask]) for a in alphas if len(a) == len(x_vals))) -def plot_fft(ax, energies: np.ndarray, signals: List[np.ndarray], directions: List[int], data_type: str): - """Plot both real and imaginary parts of FFT for given signals.""" +def plot_fft(ax, energies: np.ndarray, signals: List[np.ndarray], directions: List[int], x_range: Tuple[float, float]): labels = {0: "$x$", 1: "$y$", 2: "$z$"} for d, sig in zip(directions, signals): - sig_pos = sig[:len(energies)] # positive frequency part + sig_pos = sig[:len(energies)] ax.plot(energies, sig_pos.real, '-', label=f'Re[{labels[d]}]') ax.plot(energies, sig_pos.imag, '--', label=f'Im[{labels[d]}]') ax.set_xlabel("Energy (eV)") ax.set_ylabel("FFT Amplitude") - ax.set_xlim(0, 20) + ax.set_xlim(x_range) ax.legend() + ax.grid(True, linestyle=':', alpha=0.6) + + +def generate_triple_plots(x_vals: np.ndarray, alphas_r: List[np.ndarray], alphas_i: List[np.ndarray], + directions: List[int], xlabel: str, x_range: Tuple[float, float], + ylim: Optional[Tuple[float, float]], title_base: str, filename_base: str): + """Generate and save three standalone layout variants: Mixed components, Real-only, and Imag-only.""" + axis_labels = {0: "$x$", 1: "$y$", 2: "$z$"} + is_energy = "Energy" in xlabel + zero_cutoff = 0.2 if is_energy else 0.0 + + plot_configurations = [ + {"suffix": "", "title_append": "", "arrays": alphas_r + alphas_i, "plot_r": True, "plot_i": True}, + {"suffix": "_real", "title_append": " (Real Part)", "arrays": alphas_r, "plot_r": True, "plot_i": False}, + {"suffix": "_imag", "title_append": " (Imaginary Part)", "arrays": alphas_i, "plot_r": False, "plot_i": True}, + ] + + for config in plot_configurations: + fig, ax = plt.subplots(figsize=(10, 6)) + for idx, d in enumerate(directions): + base_label = axis_labels[d] + if config["plot_r"]: + label_str = f"Re[{base_label}]" if config["plot_i"] else base_label + ax.plot(x_vals, alphas_r[idx], '-', label=label_str) + if config["plot_i"]: + label_str = f"Im[{base_label}]" if config["plot_r"] else base_label + style = '--' if config["plot_r"] else '-' + ax.plot(x_vals, alphas_i[idx], style, label=label_str) + + ax.set_xlabel(xlabel) + ax.set_ylabel("Intensity") + ax.set_xlim(x_range) -# ====================== -# Kramers-Kronig Relations -# ====================== - -def compute_response_function(dipole_fft: np.ndarray, efield_fft: np.ndarray) -> np.ndarray: - """Compute complex response χ(ω) = D(ω) / E(ω)""" - return np.divide(dipole_fft, efield_fft, out=np.zeros_like(dipole_fft), where=efield_fft != 0) + current_ylim = ylim or calculate_y_limits(x_vals, config["arrays"], x_range, ignore_zero_energy=zero_cutoff) + if current_ylim: + ax.set_ylim(current_ylim) -def auto_ylim(data: np.ndarray, x: np.ndarray, x_min: float = 0.5, x_max: float = 20.0, percentile: float = 99.5) -> Tuple[float, float]: - """ - Automatically determine y-limits by: - 1. Restricting to energy range [x_min, x_max] - 2. Using percentile to ignore outliers - """ - mask = (x >= x_min) & (x <= x_max) - if not np.any(mask): - return (np.nanmin(data), np.nanmax(data)) - y_valid = data[mask] - y_valid = y_valid[np.isfinite(y_valid)] - if len(y_valid) == 0: - return (0, 1) + ax.legend() + ax.set_title(f"{title_base}{config['title_append']}", fontsize=18, y=1.02) + ax.grid(True, linestyle=':', alpha=0.6) + fig.savefig(f"{filename_base}{config['suffix']}.png", dpi=300, bbox_inches="tight") + plt.close(fig) - # Remove extreme outliers via percentile - y_low = np.percentile(y_valid, 100 - percentile) - y_high = np.percentile(y_valid, percentile) - margin = 0.1 * (y_high - y_low) if y_high != y_low else 0.1 * max(abs(y_high), 1) - return y_low - margin, y_high + margin -def kk_hilbert_transform(signal_pos: np.ndarray, omega: np.ndarray, is_odd: bool = True) -> np.ndarray: - """ - Compute the Hilbert transform of a function defined on positive frequencies. +def export_spectrum_data(filename: str, energies: np.ndarray, wavelengths: np.ndarray, + data_arrays: List[np.ndarray], component_labels: List[str]): + """Export calculated physical properties to a structured text file.""" + data_out = np.column_stack([energies, wavelengths] + data_arrays) + header_str = " ".join(["Energy(eV)", "Wavelength(nm)"] + component_labels) + np.savetxt(filename, data_out, header=header_str) - Parameters: - signal_pos: values of the function on omega > 0 (array of length N) - omega: positive frequency grid (must be uniformly spaced, starting near 0) - is_odd: if True, extend as odd function (for χ''); if False, extend as even (for χ') - Returns: - ht_pos: Hilbert transform on the same positive omega grid +def apply_fermi_dirac_taper(energies: np.ndarray, signal: np.ndarray, cutoff_ev: float = 40.0, smear_ev: float = 2.0, baseline: float = 0.0) -> np.ndarray: + """ + Apply a smooth Fermi-Dirac-like tapering window in the frequency domain. + Forces the optical response to smoothly decay to a specified vacuum baseline at high energies. """ - from scipy.interpolate import interp1d - from scipy.signal import hilbert as hilbert_scipy + energy_abs = np.abs(energies) + taper_window = 1.0 / (1.0 + np.exp((energy_abs - cutoff_ev) / smear_ev)) + return (signal - baseline) * taper_window + baseline + +def kk_hilbert_transform(signal_pos: np.ndarray, omega: np.ndarray, is_odd: bool = True) -> np.ndarray: N = len(omega) if N < 2: return np.zeros_like(omega) - domega = omega[1] - omega[0] - omega_max = omega[-1] - - # Full frequency grid (symmetric around 0) - omega_full = np.linspace(-omega_max, omega_max, 2 * N - 1) - - # Interpolator for positive side - interp = interp1d(omega, signal_pos, kind='linear', bounds_error=False, fill_value=0.0) - - # Extend to negative frequencies + # Construct the full symmetric frequency axis array assuming uniform spacing + omega_full = np.concatenate((-omega[-1:0:-1], omega)) signal_full = np.zeros_like(omega_full) - pos_mask = omega_full >= 0 - signal_full[pos_mask] = interp(omega_full[pos_mask]) + + # Reconstruct the full signal utilizing symmetry properties + signal_full[N - 1:] = signal_pos if is_odd: - signal_full[~pos_mask] = -interp(-omega_full[~pos_mask]) # odd extension + signal_full[:N - 1] = -signal_pos[1:][::-1] else: - signal_full[~pos_mask] = interp(-omega_full[~pos_mask]) # even extension + signal_full[:N - 1] = signal_pos[1:][::-1] - # Compute analytic signal → Hilbert transform is imaginary part - analytic = hilbert_scipy(signal_full) - ht_full = analytic.imag # this is H[signal_full] + # Extend signal with zero-padding to completely eliminate periodic boundary wrap-around + pad_length = len(signal_full) * 10 + analytic = hilbert_scipy(signal_full, N=pad_length) - # Extract positive frequencies (center at index N-1) - ht_pos = ht_full[N - 1:] + # Extract imaginary part corresponding to the Hilbert transform and truncate padding + hilbert_trans = analytic.imag[:len(signal_full)] + return hilbert_trans[N - 1:] - return ht_pos # ====================== -# Main Entry Point +# Main Execution Entry # ====================== def main(): - parser = argparse.ArgumentParser(description="Plot RT-TDDFT absorption spectrum from ABACUS output.") - parser.add_argument("--td_dt", type=float, default=0.00484, help="Time step in fs") - parser.add_argument("--direc", type=int, nargs="+", required=True, help="Field directions: 0=x,1=y,2=z") - parser.add_argument("--efield_path", type=str, nargs="+", required=True, help="E-field file paths") - parser.add_argument("--dipolefile", type=str, default="./OUT.ABACUS/dipole_s1.txt") - # parser.add_argument("--dipolefile", type=str, default="./OUT.ABACUS/dipole_s1.txt") - parser.add_argument("--material_name", type=str, - default=r"CH$_3$COOH", help="Name for plot title") + parser = argparse.ArgumentParser(description="Script for processing ABACUS RT-TDDFT absorption spectra.") + + # Operational configuration switches + parser.add_argument("--response_source", choices=["dipole", "current"], default="dipole", + help="Select polarization physical driving mechanism source. Choices: ['dipole', 'current'].") + parser.add_argument("--spectrum_type", choices=["epsilon", "sigma"], default="epsilon", + help="Select optical calculation properties parameter. Choices: ['epsilon', 'sigma'].") + + # Input file specifications + parser.add_argument("--signal_path", type=str, default=None, + help="File path to polarization input data. Configured dynamically matching response source if omitted.") + parser.add_argument("--efield_path", type=str, nargs="+", required=True, + help="File paths to electric field source files.") + + # Signal windowing parameters + parser.add_argument("--td_dt", type=float, default=0.002, help="Simulation numerical time step in fs.") + parser.add_argument("--direc", type=int, nargs="+", required=True, + help="Field simulation spatial directions list: 0=x, 1=y, 2=z.") + parser.add_argument("--material_name", type=str, default=r"C$_6$H$_6$", + help="Target description context used inside labels.") parser.add_argument("--step_start", type=int, default=0) - parser.add_argument("--step_end", type=int, default=8000) - parser.add_argument("--system_type", choices=["dipole_epsilon", "dipole_sigma", "current_epsilon", "current_sigma"], - default="dipole_sigma", help="Type of system and response for absorption calculation") - parser.add_argument("--pad_factor", type=int, default=10, help="Zero-padding factor for FFT") - parser.add_argument("--decay_beta", type=float, default=5, help="Exponential decay parameter for windowing") + parser.add_argument("--step_end", type=int, default=5000) + parser.add_argument("--pad_factor", type=int, default=10, + help="Zero padding coefficient indicator for FFT scaling.") + parser.add_argument("--decay_beta", type=float, default=5.0, + help="Damping window filter boundary argument parameter.") parser.add_argument("--remove_dc", type=bool, default=False, - help="Whether to remove DC component from time signals") + help="Enable filtering out low-frequency background signals.") + parser.add_argument("--volume", type=float, default=1.0, + help="Cell volume in atomic units (Bohr^3). Only applied to dipole response source.") + + # Visual canvas bounding limits properties + parser.add_argument("--energy_range", type=float, nargs=2, + default=[0.0, 20.0], help="X-axis display window bounds inside Energy plots.") + parser.add_argument("--wavelength_range", type=float, nargs=2, + default=[80.0, 200.0], help="X-axis display window bounds inside Wavelength plots.") + parser.add_argument("--ylim_energy", type=float, nargs=2, default=None, + help="Y-axis bounding box limits overrides inside Energy plots.") + parser.add_argument("--ylim_wavelength", type=float, nargs=2, default=None, + help="Y-axis bounding box limits overrides inside Wavelength plots.") args = parser.parse_args() + # Intelligent fallback tracking path assignments matching runtime selections + if args.signal_path is None: + if args.response_source == "dipole": + args.signal_path = "./OUT.ABACUS/dipole_s1.txt" + else: + args.signal_path = "./OUT.ABACUS/current_tot.txt" + + print(f"--> Selected physical response: [{args.response_source}]") + print(f"--> Selected spectrum variant: [{args.spectrum_type}]") + print(f"--> Target parsing data path: {args.signal_path}") + validate_inputs(args.direc, args.efield_path) Efile_groups = group_files_by_direction(args.direc, args.efield_path) - # Load data - dipole_data = load_dipole(args.dipolefile, args.step_start, args.step_end, args.remove_dc) + # Load arrays via distinct specific mechanisms + signal_data = load_signal(args.signal_path, args.step_start, args.step_end, args.remove_dc) efield_data = load_efield(Efile_groups, args.step_start, args.step_end, args.remove_dc) - # Compute absorption + # Initialize calculation core workspace absorber = AbsorptionSpectrum( - dipole=dipole_data, - efield=efield_data, - dt=args.td_dt, - system_type=args.system_type, - pad_factor=args.pad_factor, - decay_beta=args.decay_beta, + dipole=signal_data, efield=efield_data, dt=args.td_dt, + response_source=args.response_source, spectrum_type=args.spectrum_type, + pad_factor=args.pad_factor, decay_beta=args.decay_beta, volume=args.volume ) energies = absorber.get_energy_axis() - with np.errstate(divide='ignore'): wavelengths = np.where(energies > 0, (sc.h * sc.c / sc.eV * 1e9) / energies, np.inf) - directions = sorted(set(args.direc)) # unique directions to plot + directions = sorted(set(args.direc)) - # --- Save data --- - alphas = [absorber.get_positive_alpha(d) for d in range(3)] - data_out = np.column_stack([energies, wavelengths] + alphas) - np.savetxt("abs_dat.txt", data_out, header="Energy(eV) Wavelength(nm) alpha_X alpha_Y alpha_Z") + # Compile array structures across active components + alphas_real = [] + alphas_imag = [] + for d in range(3): + r, i = absorber.get_positive_alpha(d) + alphas_real.append(r) + alphas_imag.append(i) - # --- Plotting --- - time = np.arange(args.step_end - args.step_start) * args.td_dt + # Export structured raw metrics matrices to external text files + export_spectrum_data( + f"{args.spectrum_type}_dat_real.txt", energies, wavelengths, alphas_real, + ["alpha_real_X", "alpha_real_Y", "alpha_real_Z"] + ) + export_spectrum_data( + f"{args.spectrum_type}_dat_imag.txt", energies, wavelengths, alphas_imag, + ["alpha_imag_X", "alpha_imag_Y", "alpha_imag_Z"] + ) - # 1. Dipole vs time - fig, ax = plt.subplots(figsize=(10, 6)) - plot_time_series(ax, time, dipole_data, directions, ["$J_x$", "$J_y$", "$J_z$"], "Current (a.u.)") - ax.set_xlim(0, args.step_end * args.td_dt) - fig.savefig("dipole.png", dpi=300, bbox_inches="tight") + time = np.arange(args.step_end - args.step_start) * args.td_dt - # 2. Absorption spectrum - fig, ax = plt.subplots(figsize=(10, 6)) - alphas_to_plot = [absorber.get_positive_alpha(d) for d in directions] - xmax = 20 - # plot_absorption(ax, energies, alphas_to_plot, directions, x_range=(5, xmax)) - plot_absorption(ax, energies, alphas_to_plot, directions, xlabel="Energy (eV)", x_range=(4, 20)) - # Automatically set y-limits - all_alpha = np.concatenate([a for a in alphas_to_plot if len(a) == len(energies)]) - if len(all_alpha) > 0: - ymin, ymax = auto_ylim(all_alpha, energies, x_min=0.5, x_max=xmax) - ax.set_ylim(0, ymax) - ax.set_ylim(0, 100) - title = rf"{args.material_name} Absorption Spectrum" - ax.set_title(title, fontsize=20, y=1.02) - fig.savefig("abs.png", dpi=300, bbox_inches="tight") + # Adjust terminology and identifiers map dynamically based on driving response origin + if args.response_source == "dipole": + time_ylabel, time_legends, file_prefix = "Dipole (a.u.)", ["$P_x$", "$P_y$", "$P_z$"], "dipole" + else: + time_ylabel, time_legends, file_prefix = "Current (a.u.)", ["$J_x$", "$J_y$", "$J_z$"], "current" + # Figure 1: Time domain response analysis fig, ax = plt.subplots(figsize=(10, 6)) - alphas_to_plot = [absorber.get_positive_alpha(d) for d in directions] + plot_time_series(ax, time, signal_data, directions, time_legends, time_ylabel) + ax.set_xlim(0, args.step_end * args.td_dt) + fig.savefig(f"{file_prefix}.png", dpi=300, bbox_inches="tight") + plt.close(fig) - # Set wavelength range for plotting - wl_min, wl_max = 80, 200 - plot_absorption(ax, wavelengths, alphas_to_plot, directions, xlabel="Wavelength (nm)", x_range=(wl_min, wl_max)) + # Figure 2 & 3 series: Generate combinations of graphs automatically across domains + title_label = r"$\epsilon$" if args.spectrum_type == "epsilon" else r"$\sigma$" + title_string = rf"{args.material_name} {title_label} Spectrum" - all_alpha = np.concatenate([a for a in alphas_to_plot if len(a) == len(wavelengths)]) - if len(all_alpha) > 0: - ymin, ymax = auto_ylim(all_alpha, wavelengths, x_min=wl_min, x_max=wl_max) - ax.set_ylim(0, ymax) + # Render all energy component variants + generate_triple_plots(energies, [alphas_real[d] for d in directions], [alphas_imag[d] for d in directions], + directions, "Energy (eV)", tuple(args.energy_range), getattr(args, "ylim_energy", None), + title_string, args.spectrum_type) - title = rf"{args.material_name} Absorption Spectrum" - ax.set_title(title, fontsize=20, y=1.02) - fig.savefig("abs_wavelength.png", dpi=300, bbox_inches="tight") + # Render all wavelength component variants + generate_triple_plots(wavelengths, [alphas_real[d] for d in directions], [alphas_imag[d] for d in directions], + directions, "Wavelength (nm)", tuple( + args.wavelength_range), getattr(args, "ylim_wavelength", None), + title_string, f"{args.spectrum_type}_wavelength") - # 3. E-field Fourier + # Figure 4: Incident field fourier response mapping fig, ax = plt.subplots(figsize=(10, 6)) efield_ffts = [np.fft.fft(zero_pad_and_mask(efield_data[d], args.td_dt, args.pad_factor)) for d in directions] - plot_fft(ax, energies, efield_ffts, directions, "efield") + plot_fft(ax, energies, efield_ffts, directions, tuple(args.energy_range)) fig.savefig("efield_fourier.png", dpi=300, bbox_inches="tight") + plt.close(fig) - # 4. Dipole Fourier + # Figure 5: Polarization domain tracking fourier mapping fig, ax = plt.subplots(figsize=(10, 6)) - dipole_ffts = [np.fft.fft(zero_pad_and_mask(dipole_data[d], args.td_dt, args.pad_factor)) for d in directions] - plot_fft(ax, energies, dipole_ffts, directions, "dipole") - fig.savefig("dipole_fourier.png", dpi=300, bbox_inches="tight") + signal_ffts = [np.fft.fft(zero_pad_and_mask(signal_data[d], args.td_dt, args.pad_factor)) for d in directions] + plot_fft(ax, energies, signal_ffts, directions, tuple(args.energy_range)) + fig.savefig(f"{file_prefix}_fourier.png", dpi=300, bbox_inches="tight") + plt.close(fig) - # 5. E-field vs time + # Figure 6: Incident electric field input time tracking fig, ax = plt.subplots(figsize=(10, 6)) plot_time_series(ax, time, efield_data, directions, ["$E_x$", "$E_y$", "$E_z$"], "Electric Field (a.u.)") - ax.set_xlim(0, 1) + ax.set_xlim(0, 1.0) fig.savefig("efield_time.png", dpi=300, bbox_inches="tight") + plt.close(fig) - # 6. Kramers-Kronig Validation + # Figure 7: Kramers-Kronig mathematical consistency checks validation workspace n_dirs = len(directions) fig_kk, axes_kk = plt.subplots(n_dirs, 2, figsize=(14, 5 * n_dirs), squeeze=False) - for idx, d in enumerate(directions): - # --- Get complex response --- - d_pad = zero_pad_and_mask(dipole_data[d], args.td_dt, args.pad_factor, args.decay_beta) - e_pad = zero_pad_and_mask(efield_data[d], args.td_dt, args.pad_factor, args.decay_beta) - d_fft = np.fft.fft(d_pad) - e_fft = np.fft.fft(e_pad) - chi = compute_response_function(d_fft, e_fft) - chi_pos = chi[:len(energies)] + # Initialize a dedicated calculator to strictly enforce dielectric function (epsilon) evaluation + kk_absorber = AbsorptionSpectrum( + dipole=signal_data, efield=efield_data, dt=args.td_dt, + response_source=args.response_source, spectrum_type="epsilon", + pad_factor=args.pad_factor, decay_beta=args.decay_beta, volume=args.volume + ) - chi_real = chi_pos.real - chi_imag = chi_pos.imag - omega = energies # in eV + for idx, d in enumerate(directions): + eps_real, eps_imag = kk_absorber.get_positive_alpha(d) - # --- Reconstruct real part from imag --- - chi_real_KK = kk_hilbert_transform(chi_imag, omega, is_odd=True) # χ'' is odd → use odd extension + # Isolate the high-frequency numerical smearing strictly to the KK verification step + eps_real_smooth = apply_fermi_dirac_taper(energies, eps_real, baseline=1.0) + eps_imag_smooth = apply_fermi_dirac_taper(energies, eps_imag, baseline=0.0) - # --- Reconstruct imag part from real --- - chi_imag_KK = -kk_hilbert_transform(chi_real, omega, is_odd=False) # χ' is even → use even extension + # Perform Kramers-Kronig relations analysis + eps_real_KK = 1.0 - kk_hilbert_transform(eps_imag_smooth, energies, is_odd=True) + eps_imag_KK = kk_hilbert_transform(eps_real_smooth - 1.0, energies, is_odd=False) - # --- Plot: Real part comparison --- + # Real domain validation subplot plotting ax1 = axes_kk[idx, 0] - ax1.plot(omega, chi_real, 'b-', lw=1.5, label=r"$\chi'(\omega)$ (FFT)") - ax1.plot(omega, chi_real_KK, 'r--', lw=1.5, label=r"$\mathcal{H}[\chi''](\omega)$") + ax1.plot(energies, eps_real, 'b-', label=r"$\mathrm{Re}[\epsilon(\omega)]$ (TDDFT)") + ax1.plot(energies, eps_real_KK, 'r--', label=r"$1 - \mathcal{H}[\mathrm{Im}(\epsilon)](\omega)$") ax1.set_xlabel("Energy (eV)") - ax1.set_ylabel(r"$\chi'(\omega)$") - ax1.set_title(f"Direction {['$x$','$y$','$z$'][d]}: Real part KK check", fontsize=18, y=1.02) + ax1.set_ylabel(r"Re[$\epsilon(\omega)$]") + ax1.set_xlim(tuple(args.energy_range)) + + ylim_real = calculate_y_limits(energies, [eps_real, eps_real_KK], tuple(args.energy_range)) + if ylim_real: + ax1.set_ylim(ylim_real) + ax1.legend() ax1.grid(True, linestyle=':', alpha=0.6) - ax1.set_xlim(0, min(20, omega[-1])) - # Automatically set y-limits - ymin, ymax = auto_ylim(chi_real, omega, x_min=0.5, x_max=20.0) - ymin_kk, ymax_kk = auto_ylim(chi_real_KK, omega, x_min=0.5, x_max=20.0) - ax1.set_ylim(-max(abs(ymin), abs(ymin_kk)), max(ymax, ymax_kk)) - # --- Plot: Imaginary part comparison --- + # Imaginary domain validation subplot plotting ax2 = axes_kk[idx, 1] - ax2.plot(omega, chi_imag, 'g-', lw=1.5, label=r"$\chi''(\omega)$ (FFT)") - ax2.plot(omega, chi_imag_KK, 'm--', lw=1.5, label=r"$-\mathcal{H}[\chi'](\omega)$") + ax2.plot(energies, eps_imag, 'g-', label=r"$\mathrm{Im}[\epsilon(\omega)]$ (TDDFT)") + ax2.plot(energies, eps_imag_KK, 'm--', label=r"$\mathcal{H}[\mathrm{Re}(\epsilon) - 1](\omega)$") ax2.set_xlabel("Energy (eV)") - ax2.set_ylabel(r"$\chi''(\omega)$") - ax2.set_title(f"Direction {['$x$','$y$','$z$'][d]}: Imaginary part KK check", fontsize=18, y=1.02) + ax2.set_ylabel(r"Im[$\epsilon(\omega)$]") + ax2.set_xlim(tuple(args.energy_range)) + + ylim_imag = calculate_y_limits(energies, [eps_imag, eps_imag_KK], tuple(args.energy_range)) + if ylim_imag: + ax2.set_ylim(ylim_imag) + ax2.legend() ax2.grid(True, linestyle=':', alpha=0.6) - ax2.set_xlim(0, min(20, omega[-1])) - # Automatically set y-limits - ymin2, ymax2 = auto_ylim(chi_imag, omega, x_min=0.5, x_max=20.0) - ymin2_kk, ymax2_kk = auto_ylim(chi_imag_KK, omega, x_min=0.5, x_max=20.0) - ax2.set_ylim(-max(abs(ymin2), abs(ymin2_kk)), max(ymax2, ymax2_kk)) fig_kk.tight_layout() fig_kk.savefig("kk_check.png", dpi=300, bbox_inches="tight") + plt.close(fig_kk) + + print(">>> Success: All data files and multi-component plots saved successfully.") - print("Plots and data saved successfully.") if __name__ == "__main__": main() From c44de6e50c53bc1d5e4ed7b03ae65dc5514bfa85 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Thu, 11 Jun 2026 15:36:45 +0800 Subject: [PATCH 2/3] Add more error handling --- tools/rt-tddft-tools/plot_absorption.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/tools/rt-tddft-tools/plot_absorption.py b/tools/rt-tddft-tools/plot_absorption.py index 81e071957e0..18ee23f8912 100644 --- a/tools/rt-tddft-tools/plot_absorption.py +++ b/tools/rt-tddft-tools/plot_absorption.py @@ -112,7 +112,11 @@ def load_efield(file_groups: List[List[str]], step_start: int, step_end: int, re raise RuntimeError(f"Error loading efield file {f}: {e}") available_rows = arr.shape[0] - start_idx = min(step_start, available_rows) + if available_rows <= step_start: + raise ValueError( + f"Execution halted: E-field file '{f}' contains only {available_rows} rows, which is <= --step_start ({step_start})." + ) + start_idx = step_start end_idx = min(step_end, available_rows) sliced = arr[start_idx:end_idx, 1] @@ -170,6 +174,8 @@ def __init__(self, dipole: np.ndarray, efield: np.ndarray, dt: float, self.pad_factor = pad_factor self.decay_beta = decay_beta self.volume = volume + if self.response_source == "dipole" and self.volume <= 0: + raise ValueError("volume must be > 0 when --response_source is 'dipole'") self.N_orig = dipole.shape[1] self.N_pad = self.N_orig * pad_factor @@ -409,6 +415,11 @@ def main(): args = parser.parse_args() + if args.step_start < 0: + raise ValueError(f"--step_start must be >= 0 (got {args.step_start})") + if args.step_end <= args.step_start: + raise ValueError(f"--step_end ({args.step_end}) must be > --step_start ({args.step_start})") + # Intelligent fallback tracking path assignments matching runtime selections if args.signal_path is None: if args.response_source == "dipole": From 5d058185451ebd29006d5909bedb686f06b70b20 Mon Sep 17 00:00:00 2001 From: AsTonyshment Date: Thu, 11 Jun 2026 16:37:34 +0800 Subject: [PATCH 3/3] Refine efield read-in logic --- tools/rt-tddft-tools/plot_absorption.py | 52 ++++++++++++++++--------- 1 file changed, 33 insertions(+), 19 deletions(-) diff --git a/tools/rt-tddft-tools/plot_absorption.py b/tools/rt-tddft-tools/plot_absorption.py index 18ee23f8912..1abf72ad20a 100644 --- a/tools/rt-tddft-tools/plot_absorption.py +++ b/tools/rt-tddft-tools/plot_absorption.py @@ -99,8 +99,8 @@ def load_signal(filename: str, step_start: int, step_end: int, remove_dc: bool) return sliced -def load_efield(file_groups: List[List[str]], step_start: int, step_end: int, remove_dc: bool) -> np.ndarray: - """Load and sum E-field profiles. Elastic zero-padding is applied at the tail if mismatched.""" +def load_efield(file_groups: List[List[str]], step_start: int, step_end: int, dt: float, remove_dc: bool) -> np.ndarray: + """Load and sum E-field profiles, mapping them exactly to absolute simulation time.""" n_steps = step_end - step_start efield = np.zeros((3, n_steps)) @@ -111,26 +111,39 @@ def load_efield(file_groups: List[List[str]], step_start: int, step_end: int, re except OSError as e: raise RuntimeError(f"Error loading efield file {f}: {e}") - available_rows = arr.shape[0] - if available_rows <= step_start: - raise ValueError( - f"Execution halted: E-field file '{f}' contains only {available_rows} rows, which is <= --step_start ({step_start})." - ) - start_idx = step_start - end_idx = min(step_end, available_rows) + if arr.size == 0: + continue - sliced = arr[start_idx:end_idx, 1] + # Ensure 2D shape for single-row files to prevent index out of bounds + if arr.ndim == 1: + arr = arr.reshape(1, -1) - # Elastic padding mechanism: append trailing zeros if row count is slightly less than step_end - if len(sliced) < n_steps: - padded = np.zeros(n_steps) - padded[:len(sliced)] = sliced - sliced = padded + times_fs = arr[:, 0] + fields = arr[:, 1] - efield[i] += sliced / V_ANGSTROM_TO_AU_EFIELD + # Map absolute physical time to a 0-based global step index + # E.g., t = 0.006 fs with dt = 0.002 fs -> global_step = 2 + abs_step_indices = np.round(times_fs / dt).astype(int) - 1 + + # Calculate relative indices within the current analysis window [step_start, step_end) + rel_indices = abs_step_indices - step_start + + # Create boolean mask to filter out data outside the requested window + valid_mask = (rel_indices >= 0) & (rel_indices < n_steps) + + valid_rel_indices = rel_indices[valid_mask] + valid_fields = fields[valid_mask] + + # Map valid fields into the pre-allocated zero array. + # Time slots before e-field starts or after it ends remain naturally zero. + efield[i, valid_rel_indices] += valid_fields / V_ANGSTROM_TO_AU_EFIELD + + if remove_dc and n_steps > 0: + for i in range(3): + # Only subtract DC component if the channel contains a non-zero signal + if np.any(efield[i] != 0): + efield[i] -= np.mean(efield[i]) - if remove_dc and n_steps > 0: - efield[i] -= np.mean(efield[i]) return efield @@ -436,7 +449,8 @@ def main(): # Load arrays via distinct specific mechanisms signal_data = load_signal(args.signal_path, args.step_start, args.step_end, args.remove_dc) - efield_data = load_efield(Efile_groups, args.step_start, args.step_end, args.remove_dc) + # Load electric fields utilizing absolute time step mapping + efield_data = load_efield(Efile_groups, args.step_start, args.step_end, args.td_dt, args.remove_dc) # Initialize calculation core workspace absorber = AbsorptionSpectrum(