From 2b42dfe92fb7c16f8b066e5583721970a3c7c5b8 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Wed, 13 May 2026 11:50:15 +0200 Subject: [PATCH 01/11] added excph_resonant_raman.py --- .../exciton_phonon/excph_resonant_raman.py | 140 ++++++++++++++++++ 1 file changed, 140 insertions(+) create mode 100644 yambopy/exciton_phonon/excph_resonant_raman.py diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py new file mode 100644 index 00000000..a31a374d --- /dev/null +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -0,0 +1,140 @@ +## +# Authors: PMI, MN +## + +import numpy as np +from tqdm import tqdm +from yambopy.units import ha2ev + + +def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, + elec_dipoles, eph_g, n_val, cell_vol, + broad=0.1, ph_freq_threshold=5.0): + """ + Independent-particle one-phonon resonant Raman tensor at q=0 (Stokes). + + Computes the Stokes branch (phonon emission) of the third-order + light-matter / electron-phonon perturbation expression at the + independent-particle level. Resonant and anti-resonant time orderings + are both included; the anti-Stokes channel (phonon absorption) is not. + Both the electron-scattering and hole-scattering vertices are summed + with the standard relative minus sign. + + Mirrors `compute_Raman_oneph_ip` from + https://github.com/muralidhar-nalabothula/PhdScripts/blob/main/exph/raman.py + + The function takes all physical ingredients as input arrays and only + evaluates the formula; it does not load yambo databases. + + Returns + ------- + raman_tensor : (nfreqs, nmodes, 3, 3) complex ndarray + Raman tensor R^nu_{alpha,beta}(omega_L). For incoming polarization + e_in and outgoing polarization e_out the (un-normalized) intensity + is |sum_{ab} e_out[a] R[w, m, a, b] e_in[b]|^2. + + Parameters + ---------- + laser_energies : (nfreqs,) float ndarray + Incoming laser energies in eV. + ph_energies : (nmodes,) float ndarray + Phonon energies at q=0 in eV. + el_energies : (nk, nb) float ndarray + Single-particle (KS or QP) energies in eV. The first ``n_val`` bands + are valence, the remaining ``nc = nb - n_val`` are conduction + (yambo convention). + elec_dipoles : (3, nk, nc, nv) complex ndarray + Velocity-gauge electronic dipoles in atomic units, + as stored in yambo's ``ndb.dipoles``. The conduction index runs over + the ``nc`` states above ``n_val``; the valence index over the + ``n_val`` states below. + eph_g : (nmodes, nk, nb, nb) complex ndarray + Electron-phonon matrix elements at q=0 in Hartree, with index order + (mode, k, final_band, initial_band) and the standard + ``1/sqrt(2*omega_nu)`` normalization (yambo / LetzElPhC convention). + n_val : int + Number of valence bands. + cell_vol : float + Unit-cell volume in bohr^3. + broad : float, optional + Lorentzian broadening Gamma in eV. Default: 0.1. + ph_freq_threshold : float, optional + Acoustic-mode cutoff in cm^-1; modes with |omega_nu| below this + contribute zero Raman tensor. Default: 5. + """ + # cm^-1 -> Ha (Murali's constant: 1 cm^-1 = 0.12398e-3 eV) + cm1_to_eV = 0.12398e-3 + cm1_to_Ha = cm1_to_eV / ha2ev + + # Shape checks + assert el_energies.ndim == 2, "el_energies must be (nk, nb)" + nk, nb = el_energies.shape + nv = int(n_val) + nc = nb - nv + assert 0 < nv < nb, "n_val must satisfy 0 < n_val < nb" + assert elec_dipoles.shape == (3, nk, nc, nv), \ + f"elec_dipoles shape {elec_dipoles.shape} != (3, {nk}, {nc}, {nv})" + nmodes = ph_energies.shape[0] + assert eph_g.shape == (nmodes, nk, nb, nb), \ + f"eph_g shape {eph_g.shape} != ({nmodes}, {nk}, {nb}, {nb})" + assert laser_energies.ndim == 1 + assert cell_vol > 0. + + # Unit conversion to Hartree + laser_Ha = np.asarray(laser_energies) / ha2ev + ph_Ha = np.asarray(ph_energies) / ha2ev + El_Ha = np.asarray(el_energies) / ha2ev + # Murali's convention: half of input broadening enters the denominator + broad_Ha = (broad / ha2ev) / 2.0 + ph_thresh_Ha = ph_freq_threshold * cm1_to_Ha + + #TOCHECK sign on broadening + # Vertical c-v gaps minus i*Gamma : delta[k, c, v] = eps_c(k) - eps_v(k) - i*Gamma + delta = (El_Ha[:, nv:, None] - El_Ha[:, None, :nv]) - 1j * broad_Ha # (nk, nc, nv) + + # d_abs = d* (incoming-photon vertex), d_emi = d (outgoing-photon vertex) + d_abs = elec_dipoles.conj() + d_emi = elec_dipoles + + # c-c and v-v blocks of the e-ph matrix elements: index order [final, initial] + gcc = eph_g[:, :, nv:, nv:] # (nmodes, nk, nc, nc) + gvv = eph_g[:, :, :nv, :nv] # (nmodes, nk, nv, nv) + + # Output + nfreqs = laser_Ha.size + Ram = np.zeros((nfreqs, nmodes, 3, 3), dtype=np.complex128) + + # Murali's overall prefactor: 1 / N_k / sqrt(V) + ram_fac = 1.0 / nk / np.sqrt(cell_vol) + + for w_idx in tqdm(range(nfreqs), desc="IP one-phonon Raman"): + wL = laser_Ha[w_idx] + + # First-vertex dressed dipoles (no phonon shift in denominator) + dipS_res = d_abs / (wL - delta)[None, ...] # (3, nk, nc, nv) + dipS_ares = d_emi / (wL + delta)[None, ...] + + for m in range(nmodes): + wph = ph_Ha[m] + if abs(wph) <= ph_thresh_Ha: + continue + + # Second-vertex dressed dipoles (Stokes => -wph in the denominator) + dipSp_res = d_emi / (wL - delta - wph)[None, ...] + dipSp_ares = d_abs / (wL + delta - wph)[None, ...] + + gcc_m = gcc[m][None, ...] # (1, nk, nc, nc), broadcasts over polarization + gvv_m = gvv[m][None, ...] # (1, nk, nv, nv) + + # Electron-channel minus hole-channel (matmul on the last two axes) + tmp_res = gcc_m.conj() @ dipSp_res - dipSp_res @ gvv_m.conj() + tmp_ares = gcc_m @ dipSp_ares - dipSp_ares @ gvv_m + + # Sum over (k, c, v) -> Raman tensor element (alpha, beta) + Ram[w_idx, m] = np.einsum('akcv,bkcv->ab', dipS_res, tmp_res ) \ + + np.einsum('akcv,bkcv->ab', dipS_ares, tmp_ares) + + # Outgoing-photon radiation factor and BZ/volume normalization + Ram[w_idx, m] *= np.sqrt(abs(wL - wph) / wL) * ram_fac + + return Ram From 99a8090e617af6ed3d30fc34c4c9fab6b7203671 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Fri, 15 May 2026 14:35:33 +0200 Subject: [PATCH 02/11] removed shape checks --- .../exciton_phonon/excph_resonant_raman.py | 28 +++++++++---------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index a31a374d..cd26e18f 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -67,22 +67,22 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, cm1_to_Ha = cm1_to_eV / ha2ev # Shape checks - assert el_energies.ndim == 2, "el_energies must be (nk, nb)" - nk, nb = el_energies.shape - nv = int(n_val) - nc = nb - nv - assert 0 < nv < nb, "n_val must satisfy 0 < n_val < nb" - assert elec_dipoles.shape == (3, nk, nc, nv), \ - f"elec_dipoles shape {elec_dipoles.shape} != (3, {nk}, {nc}, {nv})" - nmodes = ph_energies.shape[0] - assert eph_g.shape == (nmodes, nk, nb, nb), \ - f"eph_g shape {eph_g.shape} != ({nmodes}, {nk}, {nb}, {nb})" - assert laser_energies.ndim == 1 - assert cell_vol > 0. + #assert el_energies.ndim == 2, "el_energies must be (nk, nb)" + #nk, nb = el_energies.shape + #nv = int(n_val) + #nc = nb - nv + #assert 0 < nv < nb, "n_val must satisfy 0 < n_val < nb" + #assert elec_dipoles.shape == (3, nk, nc, nv), \ + # f"elec_dipoles shape {elec_dipoles.shape} != (3, {nk}, {nc}, {nv})" + #nmodes = ph_energies.shape[0] + #assert eph_g.shape == (nmodes, nk, nb, nb), \ + # f"eph_g shape {eph_g.shape} != ({nmodes}, {nk}, {nb}, {nb})" + #assert laser_energies.ndim == 1 + #assert cell_vol > 0. # Unit conversion to Hartree - laser_Ha = np.asarray(laser_energies) / ha2ev - ph_Ha = np.asarray(ph_energies) / ha2ev + #laser_Ha = np.asarray(laser_energies) / ha2ev + #ph_Ha = np.asarray(ph_energies) / ha2ev El_Ha = np.asarray(el_energies) / ha2ev # Murali's convention: half of input broadening enters the denominator broad_Ha = (broad / ha2ev) / 2.0 From 575295c78764824058e691addc426b2db929ce40 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Fri, 15 May 2026 14:42:54 +0200 Subject: [PATCH 03/11] minor fixes --- yambopy/exciton_phonon/excph_resonant_raman.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index cd26e18f..9ef9ca57 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -68,13 +68,13 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, # Shape checks #assert el_energies.ndim == 2, "el_energies must be (nk, nb)" - #nk, nb = el_energies.shape - #nv = int(n_val) - #nc = nb - nv + nk, nb = el_energies.shape + nv = int(n_val) + nc = nb - nv #assert 0 < nv < nb, "n_val must satisfy 0 < n_val < nb" #assert elec_dipoles.shape == (3, nk, nc, nv), \ # f"elec_dipoles shape {elec_dipoles.shape} != (3, {nk}, {nc}, {nv})" - #nmodes = ph_energies.shape[0] + nmodes = ph_energies.shape[0] #assert eph_g.shape == (nmodes, nk, nb, nb), \ # f"eph_g shape {eph_g.shape} != ({nmodes}, {nk}, {nb}, {nb})" #assert laser_energies.ndim == 1 From a83eaed8012ea7f22563f4523244d951e064e8bf Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Fri, 15 May 2026 14:42:54 +0200 Subject: [PATCH 04/11] minor fixes 1 --- yambopy/exciton_phonon/excph_resonant_raman.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index cd26e18f..9ef9ca57 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -68,13 +68,13 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, # Shape checks #assert el_energies.ndim == 2, "el_energies must be (nk, nb)" - #nk, nb = el_energies.shape - #nv = int(n_val) - #nc = nb - nv + nk, nb = el_energies.shape + nv = int(n_val) + nc = nb - nv #assert 0 < nv < nb, "n_val must satisfy 0 < n_val < nb" #assert elec_dipoles.shape == (3, nk, nc, nv), \ # f"elec_dipoles shape {elec_dipoles.shape} != (3, {nk}, {nc}, {nv})" - #nmodes = ph_energies.shape[0] + nmodes = ph_energies.shape[0] #assert eph_g.shape == (nmodes, nk, nb, nb), \ # f"eph_g shape {eph_g.shape} != ({nmodes}, {nk}, {nb}, {nb})" #assert laser_energies.ndim == 1 From dd214a77a3fb00645cbd6a2ee316cba7b66daea3 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Fri, 15 May 2026 15:16:14 +0200 Subject: [PATCH 05/11] idk --- .../exciton_phonon/excph_resonant_raman.py | 203 ++++++++++++------ 1 file changed, 142 insertions(+), 61 deletions(-) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index 9ef9ca57..9660ead1 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -9,7 +9,8 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, elec_dipoles, eph_g, n_val, cell_vol, - broad=0.1, ph_freq_threshold=5.0): + broad=0.1, ph_freq_threshold=5.0, + eph_bands=None): """ Independent-particle one-phonon resonant Raman tensor at q=0 (Stokes). @@ -20,11 +21,17 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, Both the electron-scattering and hole-scattering vertices are summed with the standard relative minus sign. + The electron-phonon matrix elements (``eph_g``) may cover a band + window that is *smaller than* (or otherwise different from) the one + spanned by ``el_energies`` / ``elec_dipoles``. Use the ``eph_bands`` + argument to specify, for each band axis of ``eph_g``, which band of + ``el_energies`` it corresponds to. Bands in ``el_energies`` that are + not listed in ``eph_bands`` contribute to the optical (dipole) sums + but not to the phonon-mediated scattering. + Mirrors `compute_Raman_oneph_ip` from https://github.com/muralidhar-nalabothula/PhdScripts/blob/main/exph/raman.py - - The function takes all physical ingredients as input arrays and only - evaluates the formula; it does not load yambo databases. + (which is the special case ``eph_bands = np.arange(nb)``). Returns ------- @@ -40,20 +47,20 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, ph_energies : (nmodes,) float ndarray Phonon energies at q=0 in eV. el_energies : (nk, nb) float ndarray - Single-particle (KS or QP) energies in eV. The first ``n_val`` bands - are valence, the remaining ``nc = nb - n_val`` are conduction - (yambo convention). + Single-particle (KS or QP) energies in eV. The first ``n_val`` + bands are valence, the remaining ``nc = nb - n_val`` are + conduction (yambo convention). elec_dipoles : (3, nk, nc, nv) complex ndarray - Velocity-gauge electronic dipoles in atomic units, - as stored in yambo's ``ndb.dipoles``. The conduction index runs over - the ``nc`` states above ``n_val``; the valence index over the - ``n_val`` states below. - eph_g : (nmodes, nk, nb, nb) complex ndarray - Electron-phonon matrix elements at q=0 in Hartree, with index order - (mode, k, final_band, initial_band) and the standard - ``1/sqrt(2*omega_nu)`` normalization (yambo / LetzElPhC convention). + Velocity-gauge electronic dipoles in atomic + units, as stored in yambo's ``ndb.dipoles``. + eph_g : (nmodes, nk, nb_eph, nb_eph) complex ndarray + Electron-phonon matrix elements at q=0 in Hartree, with index + order (mode, k, final_band, initial_band) and the standard + ``1/sqrt(2*omega_nu)`` normalization. The matrix is square in + its last two axes; ``nb_eph`` can be less than or equal to + ``nb``. n_val : int - Number of valence bands. + Number of valence bands inside the ``el_energies`` window. cell_vol : float Unit-cell volume in bohr^3. broad : float, optional @@ -61,34 +68,71 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, ph_freq_threshold : float, optional Acoustic-mode cutoff in cm^-1; modes with |omega_nu| below this contribute zero Raman tensor. Default: 5. + eph_bands : array-like of int, optional + Length ``nb_eph``. ``eph_bands[i]`` is the index (within the + band axis of ``el_energies``, in [0, nb)) of the band that the + i-th band axis of ``eph_g`` corresponds to. Indices in + ``[0, n_val)`` are interpreted as valence, indices in + ``[n_val, nb)`` as conduction. If ``None`` (default), ``eph_g`` + must cover all ``nb`` bands of ``el_energies`` in their natural + order. """ # cm^-1 -> Ha (Murali's constant: 1 cm^-1 = 0.12398e-3 eV) cm1_to_eV = 0.12398e-3 cm1_to_Ha = cm1_to_eV / ha2ev - # Shape checks - #assert el_energies.ndim == 2, "el_energies must be (nk, nb)" + # ---- Input validation (raise so it survives `python -O`) --------------- + if el_energies.ndim != 2: + raise ValueError("el_energies must be (nk, nb)") nk, nb = el_energies.shape nv = int(n_val) nc = nb - nv - #assert 0 < nv < nb, "n_val must satisfy 0 < n_val < nb" - #assert elec_dipoles.shape == (3, nk, nc, nv), \ - # f"elec_dipoles shape {elec_dipoles.shape} != (3, {nk}, {nc}, {nv})" + if not (0 < nv < nb): + raise ValueError("n_val must satisfy 0 < n_val < nb; got n_val=%d, nb=%d" + % (nv, nb)) + if tuple(elec_dipoles.shape) != (3, nk, nc, nv): + raise ValueError("elec_dipoles shape %s != (3, %d, %d, %d)" + % (tuple(elec_dipoles.shape), nk, nc, nv)) nmodes = ph_energies.shape[0] - #assert eph_g.shape == (nmodes, nk, nb, nb), \ - # f"eph_g shape {eph_g.shape} != ({nmodes}, {nk}, {nb}, {nb})" - #assert laser_energies.ndim == 1 - #assert cell_vol > 0. - - # Unit conversion to Hartree - #laser_Ha = np.asarray(laser_energies) / ha2ev - #ph_Ha = np.asarray(ph_energies) / ha2ev - El_Ha = np.asarray(el_energies) / ha2ev + if eph_g.ndim != 4 or eph_g.shape[0] != nmodes or eph_g.shape[1] != nk: + raise ValueError("eph_g shape %s incompatible with (nmodes=%d, nk=%d, nb_eph, nb_eph)" + % (tuple(eph_g.shape), nmodes, nk)) + if eph_g.shape[2] != eph_g.shape[3]: + raise ValueError("eph_g must be square in its last two axes; got %s" + % (tuple(eph_g.shape),)) + nb_eph = eph_g.shape[2] + + if eph_bands is None: + if nb_eph != nb: + raise ValueError( + "eph_g spans nb_eph=%d bands but el_energies spans nb=%d; " + "pass `eph_bands` (length nb_eph) to specify the mapping." + % (nb_eph, nb)) + eph_bands = np.arange(nb) + else: + eph_bands = np.asarray(eph_bands, dtype=int) + if eph_bands.shape != (nb_eph,): + raise ValueError("eph_bands must have shape (%d,); got %s" + % (nb_eph, tuple(eph_bands.shape))) + if (eph_bands < 0).any() or (eph_bands >= nb).any(): + raise ValueError("eph_bands must contain band indices in [0, %d); got min=%d, max=%d" + % (nb, int(eph_bands.min()), int(eph_bands.max()))) + if len(np.unique(eph_bands)) != nb_eph: + raise ValueError("eph_bands must not contain duplicate band indices") + + if laser_energies.ndim != 1: + raise ValueError("laser_energies must be 1D") + if cell_vol <= 0.: + raise ValueError("cell_vol must be positive") + + # ---- Unit conversion --------------------------------------------------- + laser_Ha = np.asarray(laser_energies) / ha2ev + ph_Ha = np.asarray(ph_energies) / ha2ev + El_Ha = np.asarray(el_energies) / ha2ev # Murali's convention: half of input broadening enters the denominator - broad_Ha = (broad / ha2ev) / 2.0 + broad_Ha = (broad / ha2ev) / 2.0 ph_thresh_Ha = ph_freq_threshold * cm1_to_Ha - #TOCHECK sign on broadening # Vertical c-v gaps minus i*Gamma : delta[k, c, v] = eps_c(k) - eps_v(k) - i*Gamma delta = (El_Ha[:, nv:, None] - El_Ha[:, None, :nv]) - 1j * broad_Ha # (nk, nc, nv) @@ -96,45 +140,82 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, d_abs = elec_dipoles.conj() d_emi = elec_dipoles - # c-c and v-v blocks of the e-ph matrix elements: index order [final, initial] - gcc = eph_g[:, :, nv:, nv:] # (nmodes, nk, nc, nc) - gvv = eph_g[:, :, :nv, :nv] # (nmodes, nk, nv, nv) - - # Output - nfreqs = laser_Ha.size - Ram = np.zeros((nfreqs, nmodes, 3, 3), dtype=np.complex128) - - # Murali's overall prefactor: 1 / N_k / sqrt(V) + # ---- Identify elph valence / conduction subsets ------------------------ + v_mask = (eph_bands < nv) # (nb_eph,) bool + c_mask = ~v_mask + eph_v_local = eph_bands[v_mask] # indices in [0, nv) + eph_c_local = eph_bands[c_mask] - nv # indices in [0, nc) + n_eph_v = int(v_mask.sum()) + n_eph_c = int(c_mask.sum()) + + # cc and vv blocks of eph_g, with the same ordering as eph_c_local / + # eph_v_local (so the contractions below stay consistent). + gcc = eph_g[:, :, c_mask, :][:, :, :, c_mask] # (nmodes, nk, n_eph_c, n_eph_c) + gvv = eph_g[:, :, v_mask, :][:, :, :, v_mask] # (nmodes, nk, n_eph_v, n_eph_v) + + # ---- Sliced dipoles / gaps per channel --------------------------------- + # Electron channel: c restricted to elph_c, v spans full dipole-v. + if n_eph_c > 0: + delta_e = delta[:, eph_c_local, :] # (nk, n_eph_c, nv) + d_abs_e = d_abs[:, :, eph_c_local, :] # (3, nk, n_eph_c, nv) + d_emi_e = d_emi[:, :, eph_c_local, :] + # Hole channel: c spans full dipole-c, v restricted to elph_v. + if n_eph_v > 0: + delta_h = delta[:, :, eph_v_local] # (nk, nc, n_eph_v) + d_abs_h = d_abs[:, :, :, eph_v_local] # (3, nk, nc, n_eph_v) + d_emi_h = d_emi[:, :, :, eph_v_local] + + # ---- Output ------------------------------------------------------------ + nfreqs = laser_Ha.size + Ram = np.zeros((nfreqs, nmodes, 3, 3), dtype=np.complex128) ram_fac = 1.0 / nk / np.sqrt(cell_vol) for w_idx in tqdm(range(nfreqs), desc="IP one-phonon Raman"): wL = laser_Ha[w_idx] - # First-vertex dressed dipoles (no phonon shift in denominator) - dipS_res = d_abs / (wL - delta)[None, ...] # (3, nk, nc, nv) - dipS_ares = d_emi / (wL + delta)[None, ...] + # First-vertex dressed dipoles (no phonon shift in the denominator) + if n_eph_c > 0: + dipS_res_e = d_abs_e / (wL - delta_e)[None, ...] # (3, nk, n_eph_c, nv) + dipS_ares_e = d_emi_e / (wL + delta_e)[None, ...] + if n_eph_v > 0: + dipS_res_h = d_abs_h / (wL - delta_h)[None, ...] # (3, nk, nc, n_eph_v) + dipS_ares_h = d_emi_h / (wL + delta_h)[None, ...] for m in range(nmodes): wph = ph_Ha[m] if abs(wph) <= ph_thresh_Ha: continue - # Second-vertex dressed dipoles (Stokes => -wph in the denominator) - dipSp_res = d_emi / (wL - delta - wph)[None, ...] - dipSp_ares = d_abs / (wL + delta - wph)[None, ...] - - gcc_m = gcc[m][None, ...] # (1, nk, nc, nc), broadcasts over polarization - gvv_m = gvv[m][None, ...] # (1, nk, nv, nv) - - # Electron-channel minus hole-channel (matmul on the last two axes) - tmp_res = gcc_m.conj() @ dipSp_res - dipSp_res @ gvv_m.conj() - tmp_ares = gcc_m @ dipSp_ares - dipSp_ares @ gvv_m - - # Sum over (k, c, v) -> Raman tensor element (alpha, beta) - Ram[w_idx, m] = np.einsum('akcv,bkcv->ab', dipS_res, tmp_res ) \ - + np.einsum('akcv,bkcv->ab', dipS_ares, tmp_ares) - - # Outgoing-photon radiation factor and BZ/volume normalization - Ram[w_idx, m] *= np.sqrt(abs(wL - wph) / wL) * ram_fac + R_e = 0.0 + R_h = 0.0 + + # ---- Electron channel (g_cc) ------------------------------- + if n_eph_c > 0: + dipSp_res_e = d_emi_e / (wL - delta_e - wph)[None, ...] + dipSp_ares_e = d_abs_e / (wL + delta_e - wph)[None, ...] + gcc_m = gcc[m] # (nk, n_eph_c, n_eph_c) + # tmp[a, k, c_final, v] = sum_{c_init} g[k, c_final, c_init].conj + # * dipSp[a, k, c_init, v] + tmp_e_res = np.einsum('kij,lkjv->lkiv', gcc_m.conj(), dipSp_res_e) + tmp_e_ares = np.einsum('kij,lkjv->lkiv', gcc_m, dipSp_ares_e) + # Outer contraction over (k, c_final, v). + R_e = (np.einsum('akcv,bkcv->ab', dipS_res_e, tmp_e_res) + + np.einsum('akcv,bkcv->ab', dipS_ares_e, tmp_e_ares)) + + # ---- Hole channel (g_vv) ----------------------------------- + if n_eph_v > 0: + dipSp_res_h = d_emi_h / (wL - delta_h - wph)[None, ...] + dipSp_ares_h = d_abs_h / (wL + delta_h - wph)[None, ...] + gvv_m = gvv[m] # (nk, n_eph_v, n_eph_v) + # tmp[a, k, c, v_init] = sum_{v_final} dipSp[a, k, c, v_final] + # * g[k, v_final, v_init].conj + tmp_h_res = np.einsum('lkcv,kvw->lkcw', dipSp_res_h, gvv_m.conj()) + tmp_h_ares = np.einsum('lkcv,kvw->lkcw', dipSp_ares_h, gvv_m) + R_h = (np.einsum('akcv,bkcv->ab', dipS_res_h, tmp_h_res) + + np.einsum('akcv,bkcv->ab', dipS_ares_h, tmp_h_ares)) + + # Electron minus hole, with outgoing-photon radiation factor + # and BZ / volume normalisation. + Ram[w_idx, m] = (R_e - R_h) * np.sqrt(abs(wL - wph) / wL) * ram_fac return Ram From 46485f66e1c9e3589f3800054fcced694c9078d0 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Mon, 18 May 2026 15:27:50 +0200 Subject: [PATCH 06/11] rewriting again --- tutorial/exciton-phonon/ip_raman.py | 116 ++++++++++ .../exciton_phonon/excph_resonant_raman.py | 218 +++++------------- 2 files changed, 169 insertions(+), 165 deletions(-) create mode 100644 tutorial/exciton-phonon/ip_raman.py diff --git a/tutorial/exciton-phonon/ip_raman.py b/tutorial/exciton-phonon/ip_raman.py new file mode 100644 index 00000000..01f5dc6a --- /dev/null +++ b/tutorial/exciton-phonon/ip_raman.py @@ -0,0 +1,116 @@ +""" +IP one-phonon Stokes Raman with yambopy. + +The driver does all the band-window bookkeeping: + - loads lattice, electrons, dipoles, ndb.elph + - BZ-expands the dipoles via YamboDipolesDB(expand=True) and reorders + them into the LetzElPhC k-grid + - slices electron energies, dipoles and gkkp to the user-chosen + `raman_bands` window (Fortran 1-indexed, inclusive) + - passes ready-to-use arrays to ip_resonant_raman_oneph + +The Raman function itself takes only the energies, dipoles, gkkp, +cell volume and broadening and evaluates the formula. + +Assumes yambo wrote velocity-gauge dipoles (DIP_v in ndb.dipoles). +""" + +import numpy as np + +from yambopy import (YamboLatticeDB, YamboElectronsDB, + YamboDipolesDB, LetzElphElectronPhononDB) +from yambopy.exciton_phonon.excph_resonant_raman import ip_resonant_raman_oneph + + +# ---- User inputs ---------------------------------------------------------- +folder = '.' +SAVE_dir = folder + '/SAVE' +elph_path = folder + '/ndb.elph' +dipole_path = folder + '/dipoles/ndb.dipoles' + +omega_range = (1.0, 6.0, 501) # laser energies (eV): min, max, npts +broading = 0.10 # Lorentzian broadening (eV) +ph_fre_th = 5.0 # acoustic-mode cutoff (cm^-1) +raman_bands = (47, 56) # 1-indexed Fortran inclusive + + +# ---- Databases ------------------------------------------------------------ +lattice = YamboLatticeDB.from_db_file(filename='%s/ns.db1' % SAVE_dir) +electrons = YamboElectronsDB.from_db_file(folder=SAVE_dir, + filename='ns.db1', Expand=True) +dipdb = YamboDipolesDB.from_db_file(lattice, + filename=dipole_path, + dip_type='v', + project=False, + expand=True) +elph = LetzElphElectronPhononDB(elph_path, read_all=True) + + +# ---- Yambo BZ -> elph BZ permutation -------------------------------------- +diff = elph.kpoints[:, None, :] - lattice.red_kpoints[None, :, :] +diff = diff - np.round(diff) +perm = np.argmin(np.linalg.norm(diff, axis=-1), axis=1) + + +# ---- Resolve band ranges (1-indexed inclusive, Fortran) ------------------- +b_in_F, b_out_F = raman_bands +nb = b_out_F - b_in_F + 1 +n_val = int(electrons.nbandsv) - (b_in_F - 1) +nc = nb - n_val +cell_vol = float(lattice.lat_vol) + + +# ---- eph_g(q=0) sliced to the raman window -------------------------------- +elph_b_in_F = int(min(elph.bands)) +off = b_in_F - elph_b_in_F + +iq0 = 0 +g_q0 = elph.gkkp[iq0, :, :, 0, :, :] # (nk, nmodes, nb_i, nb_f) +g_q0 = np.swapaxes(g_q0, -1, -2) +eph_g = np.transpose(g_q0, (1, 0, 2, 3)) / 2.0 # (nmodes, nk, nb_eph, nb_eph), Ha +eph_g = eph_g[:, :, off : off + nb, off : off + nb] # (nmodes, nk, nb, nb) + +ph_eV = elph.ph_energies[iq0] # (nmodes,), eV + + +# ---- KS energies in elph k-order, sliced to the raman window -------------- +el_eV = electrons.eigenvalues[0][perm, b_in_F - 1 : b_out_F] # (nk, nb) + + +# ---- Velocity dipoles: cv block, reorder, slice to raman cv block -------- +nv_dip = int(dipdb.nbandsv) +nc_dip = int(dipdb.nbandsc) +elec_dipoles = dipdb.dipoles[:, :, nv_dip : nv_dip + nc_dip, :nv_dip] # (nk_y, 3, nc_dip, nv_dip) +elec_dipoles = elec_dipoles[perm] # -> elph k-order +elec_dipoles = elec_dipoles[:, :, :nc, -n_val:] # raman cv block +elec_dipoles = np.transpose(elec_dipoles, (1, 0, 2, 3)) # (3, nk, nc, n_val) + + +# ---- Raman calculation ---------------------------------------------------- +laser_eV = np.linspace(omega_range[0], omega_range[1], int(omega_range[2])) + +R = ip_resonant_raman_oneph( + laser_energies = laser_eV, + ph_energies = ph_eV, + el_energies = el_eV, + elec_dipoles = elec_dipoles, + eph_g = eph_g, + cell_vol = cell_vol, + broad = broading, + ph_freq_threshold = ph_fre_th, +) + + +# ---- Save ----------------------------------------------------------------- +raman_inte = np.sum(np.abs(R)**2, axis=(1, 2, 3)) + +np.savetxt('raman_intensity.dat', + np.column_stack([laser_eV, raman_inte]), + header='laser_energy(eV) raman_intensity(arb.u.)') + +np.savez('Raman_tensors.npz', + laser_energies_eV = laser_eV, + ph_freq_cm = ph_eV * 8065.54, + raman_tensor = R) + +print('Wrote raman_intensity.dat and Raman_tensors.npz') diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index 9660ead1..a05ebfea 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -7,38 +7,20 @@ from yambopy.units import ha2ev +@citation("PHYSICAL REVIEW B 113, 085201 (2026)") def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, - elec_dipoles, eph_g, n_val, cell_vol, - broad=0.1, ph_freq_threshold=5.0, - eph_bands=None): + elec_dipoles, eph_g, cell_vol, + broad=0.1, ph_freq_threshold=5.0): """ - Independent-particle one-phonon resonant Raman tensor at q=0 (Stokes). - - Computes the Stokes branch (phonon emission) of the third-order - light-matter / electron-phonon perturbation expression at the - independent-particle level. Resonant and anti-resonant time orderings - are both included; the anti-Stokes channel (phonon absorption) is not. - Both the electron-scattering and hole-scattering vertices are summed - with the standard relative minus sign. - - The electron-phonon matrix elements (``eph_g``) may cover a band - window that is *smaller than* (or otherwise different from) the one - spanned by ``el_energies`` / ``elec_dipoles``. Use the ``eph_bands`` - argument to specify, for each band axis of ``eph_g``, which band of - ``el_energies`` it corresponds to. Bands in ``el_energies`` that are - not listed in ``eph_bands`` contribute to the optical (dipole) sums - but not to the phonon-mediated scattering. + 1-phonon Raman tensor at independent-particle (IP) level (Stokes). Mirrors `compute_Raman_oneph_ip` from https://github.com/muralidhar-nalabothula/PhdScripts/blob/main/exph/raman.py - (which is the special case ``eph_bands = np.arange(nb)``). - Returns - ------- - raman_tensor : (nfreqs, nmodes, 3, 3) complex ndarray - Raman tensor R^nu_{alpha,beta}(omega_L). For incoming polarization - e_in and outgoing polarization e_out the (un-normalized) intensity - is |sum_{ab} e_out[a] R[w, m, a, b] e_in[b]|^2. + All band-window bookkeeping (which bands, IBZ -> BZ expansion, etc.) + is the caller's responsibility. The function assumes that every input + is already sliced and aligned onto the same nb-band window in the + order (valence, then conduction) along the band axis. Parameters ---------- @@ -47,125 +29,55 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, ph_energies : (nmodes,) float ndarray Phonon energies at q=0 in eV. el_energies : (nk, nb) float ndarray - Single-particle (KS or QP) energies in eV. The first ``n_val`` - bands are valence, the remaining ``nc = nb - n_val`` are - conduction (yambo convention). + Single-particle energies in eV, valence-first then conduction. elec_dipoles : (3, nk, nc, nv) complex ndarray - Velocity-gauge electronic dipoles in atomic - units, as stored in yambo's ``ndb.dipoles``. - eph_g : (nmodes, nk, nb_eph, nb_eph) complex ndarray + Velocity-gauge dipoles (yambo convention, a.u.). + nv + nc must equal nb. + eph_g : (nmodes, nk, nb, nb) complex ndarray Electron-phonon matrix elements at q=0 in Hartree, with index order (mode, k, final_band, initial_band) and the standard - ``1/sqrt(2*omega_nu)`` normalization. The matrix is square in - its last two axes; ``nb_eph`` can be less than or equal to - ``nb``. - n_val : int - Number of valence bands inside the ``el_energies`` window. + 1/sqrt(2*omega_nu) normalization. cell_vol : float Unit-cell volume in bohr^3. broad : float, optional Lorentzian broadening Gamma in eV. Default: 0.1. ph_freq_threshold : float, optional - Acoustic-mode cutoff in cm^-1; modes with |omega_nu| below this - contribute zero Raman tensor. Default: 5. - eph_bands : array-like of int, optional - Length ``nb_eph``. ``eph_bands[i]`` is the index (within the - band axis of ``el_energies``, in [0, nb)) of the band that the - i-th band axis of ``eph_g`` corresponds to. Indices in - ``[0, n_val)`` are interpreted as valence, indices in - ``[n_val, nb)`` as conduction. If ``None`` (default), ``eph_g`` - must cover all ``nb`` bands of ``el_energies`` in their natural - order. + Acoustic-mode cutoff in cm^-1. Default: 5. + + Returns + ------- + raman_tensor : (nfreqs, nmodes, 3, 3) complex ndarray """ - # cm^-1 -> Ha (Murali's constant: 1 cm^-1 = 0.12398e-3 eV) - cm1_to_eV = 0.12398e-3 - cm1_to_Ha = cm1_to_eV / ha2ev - - # ---- Input validation (raise so it survives `python -O`) --------------- - if el_energies.ndim != 2: - raise ValueError("el_energies must be (nk, nb)") - nk, nb = el_energies.shape - nv = int(n_val) - nc = nb - nv - if not (0 < nv < nb): - raise ValueError("n_val must satisfy 0 < n_val < nb; got n_val=%d, nb=%d" - % (nv, nb)) - if tuple(elec_dipoles.shape) != (3, nk, nc, nv): + cm1_to_Ha = 0.12398e-3 / ha2ev + + nk, nb = el_energies.shape + _, _, nc, nv = elec_dipoles.shape + nmodes = ph_energies.shape[0] + if nv + nc != nb: + raise ValueError("nv + nc != nb (%d + %d != %d)" % (nv, nc, nb)) + if elec_dipoles.shape != (3, nk, nc, nv): raise ValueError("elec_dipoles shape %s != (3, %d, %d, %d)" - % (tuple(elec_dipoles.shape), nk, nc, nv)) - nmodes = ph_energies.shape[0] - if eph_g.ndim != 4 or eph_g.shape[0] != nmodes or eph_g.shape[1] != nk: - raise ValueError("eph_g shape %s incompatible with (nmodes=%d, nk=%d, nb_eph, nb_eph)" - % (tuple(eph_g.shape), nmodes, nk)) - if eph_g.shape[2] != eph_g.shape[3]: - raise ValueError("eph_g must be square in its last two axes; got %s" - % (tuple(eph_g.shape),)) - nb_eph = eph_g.shape[2] - - if eph_bands is None: - if nb_eph != nb: - raise ValueError( - "eph_g spans nb_eph=%d bands but el_energies spans nb=%d; " - "pass `eph_bands` (length nb_eph) to specify the mapping." - % (nb_eph, nb)) - eph_bands = np.arange(nb) - else: - eph_bands = np.asarray(eph_bands, dtype=int) - if eph_bands.shape != (nb_eph,): - raise ValueError("eph_bands must have shape (%d,); got %s" - % (nb_eph, tuple(eph_bands.shape))) - if (eph_bands < 0).any() or (eph_bands >= nb).any(): - raise ValueError("eph_bands must contain band indices in [0, %d); got min=%d, max=%d" - % (nb, int(eph_bands.min()), int(eph_bands.max()))) - if len(np.unique(eph_bands)) != nb_eph: - raise ValueError("eph_bands must not contain duplicate band indices") - - if laser_energies.ndim != 1: - raise ValueError("laser_energies must be 1D") - if cell_vol <= 0.: - raise ValueError("cell_vol must be positive") - - # ---- Unit conversion --------------------------------------------------- + % (elec_dipoles.shape, nk, nc, nv)) + if eph_g.shape != (nmodes, nk, nb, nb): + raise ValueError("eph_g shape %s != (%d, %d, %d, %d)" + % (eph_g.shape, nmodes, nk, nb, nb)) + + # Unit conversion (eV -> Ha) laser_Ha = np.asarray(laser_energies) / ha2ev ph_Ha = np.asarray(ph_energies) / ha2ev El_Ha = np.asarray(el_energies) / ha2ev - # Murali's convention: half of input broadening enters the denominator broad_Ha = (broad / ha2ev) / 2.0 ph_thresh_Ha = ph_freq_threshold * cm1_to_Ha - # Vertical c-v gaps minus i*Gamma : delta[k, c, v] = eps_c(k) - eps_v(k) - i*Gamma + # Vertical c-v gaps minus i*Gamma delta = (El_Ha[:, nv:, None] - El_Ha[:, None, :nv]) - 1j * broad_Ha # (nk, nc, nv) - # d_abs = d* (incoming-photon vertex), d_emi = d (outgoing-photon vertex) - d_abs = elec_dipoles.conj() - d_emi = elec_dipoles - - # ---- Identify elph valence / conduction subsets ------------------------ - v_mask = (eph_bands < nv) # (nb_eph,) bool - c_mask = ~v_mask - eph_v_local = eph_bands[v_mask] # indices in [0, nv) - eph_c_local = eph_bands[c_mask] - nv # indices in [0, nc) - n_eph_v = int(v_mask.sum()) - n_eph_c = int(c_mask.sum()) - - # cc and vv blocks of eph_g, with the same ordering as eph_c_local / - # eph_v_local (so the contractions below stay consistent). - gcc = eph_g[:, :, c_mask, :][:, :, :, c_mask] # (nmodes, nk, n_eph_c, n_eph_c) - gvv = eph_g[:, :, v_mask, :][:, :, :, v_mask] # (nmodes, nk, n_eph_v, n_eph_v) - - # ---- Sliced dipoles / gaps per channel --------------------------------- - # Electron channel: c restricted to elph_c, v spans full dipole-v. - if n_eph_c > 0: - delta_e = delta[:, eph_c_local, :] # (nk, n_eph_c, nv) - d_abs_e = d_abs[:, :, eph_c_local, :] # (3, nk, n_eph_c, nv) - d_emi_e = d_emi[:, :, eph_c_local, :] - # Hole channel: c spans full dipole-c, v restricted to elph_v. - if n_eph_v > 0: - delta_h = delta[:, :, eph_v_local] # (nk, nc, n_eph_v) - d_abs_h = d_abs[:, :, :, eph_v_local] # (3, nk, nc, n_eph_v) - d_emi_h = d_emi[:, :, :, eph_v_local] - - # ---- Output ------------------------------------------------------------ + d_abs = elec_dipoles.conj() # incoming-photon vertex + d_emi = elec_dipoles # outgoing-photon vertex + + gcc = eph_g[:, :, nv:, nv:] # (nmodes, nk, nc, nc) + gvv = eph_g[:, :, :nv, :nv] # (nmodes, nk, nv, nv) + nfreqs = laser_Ha.size Ram = np.zeros((nfreqs, nmodes, 3, 3), dtype=np.complex128) ram_fac = 1.0 / nk / np.sqrt(cell_vol) @@ -173,49 +85,25 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, for w_idx in tqdm(range(nfreqs), desc="IP one-phonon Raman"): wL = laser_Ha[w_idx] - # First-vertex dressed dipoles (no phonon shift in the denominator) - if n_eph_c > 0: - dipS_res_e = d_abs_e / (wL - delta_e)[None, ...] # (3, nk, n_eph_c, nv) - dipS_ares_e = d_emi_e / (wL + delta_e)[None, ...] - if n_eph_v > 0: - dipS_res_h = d_abs_h / (wL - delta_h)[None, ...] # (3, nk, nc, n_eph_v) - dipS_ares_h = d_emi_h / (wL + delta_h)[None, ...] + dipS_res = d_abs / (wL - delta)[None, ...] + dipS_ares = d_emi / (wL + delta)[None, ...] for m in range(nmodes): wph = ph_Ha[m] if abs(wph) <= ph_thresh_Ha: continue - R_e = 0.0 - R_h = 0.0 - - # ---- Electron channel (g_cc) ------------------------------- - if n_eph_c > 0: - dipSp_res_e = d_emi_e / (wL - delta_e - wph)[None, ...] - dipSp_ares_e = d_abs_e / (wL + delta_e - wph)[None, ...] - gcc_m = gcc[m] # (nk, n_eph_c, n_eph_c) - # tmp[a, k, c_final, v] = sum_{c_init} g[k, c_final, c_init].conj - # * dipSp[a, k, c_init, v] - tmp_e_res = np.einsum('kij,lkjv->lkiv', gcc_m.conj(), dipSp_res_e) - tmp_e_ares = np.einsum('kij,lkjv->lkiv', gcc_m, dipSp_ares_e) - # Outer contraction over (k, c_final, v). - R_e = (np.einsum('akcv,bkcv->ab', dipS_res_e, tmp_e_res) - + np.einsum('akcv,bkcv->ab', dipS_ares_e, tmp_e_ares)) - - # ---- Hole channel (g_vv) ----------------------------------- - if n_eph_v > 0: - dipSp_res_h = d_emi_h / (wL - delta_h - wph)[None, ...] - dipSp_ares_h = d_abs_h / (wL + delta_h - wph)[None, ...] - gvv_m = gvv[m] # (nk, n_eph_v, n_eph_v) - # tmp[a, k, c, v_init] = sum_{v_final} dipSp[a, k, c, v_final] - # * g[k, v_final, v_init].conj - tmp_h_res = np.einsum('lkcv,kvw->lkcw', dipSp_res_h, gvv_m.conj()) - tmp_h_ares = np.einsum('lkcv,kvw->lkcw', dipSp_ares_h, gvv_m) - R_h = (np.einsum('akcv,bkcv->ab', dipS_res_h, tmp_h_res) - + np.einsum('akcv,bkcv->ab', dipS_ares_h, tmp_h_ares)) - - # Electron minus hole, with outgoing-photon radiation factor - # and BZ / volume normalisation. - Ram[w_idx, m] = (R_e - R_h) * np.sqrt(abs(wL - wph) / wL) * ram_fac + dipSp_res = d_emi / (wL - delta - wph)[None, ...] + dipSp_ares = d_abs / (wL + delta - wph)[None, ...] + + gcc_m = gcc[m][None, ...] + gvv_m = gvv[m][None, ...] + + tmp_res = gcc_m.conj() @ dipSp_res - dipSp_res @ gvv_m.conj() + tmp_ares = gcc_m @ dipSp_ares - dipSp_ares @ gvv_m + + Ram[w_idx, m] = (np.einsum('akcv,bkcv->ab', dipS_res, tmp_res) + + np.einsum('akcv,bkcv->ab', dipS_ares, tmp_ares)) + Ram[w_idx, m] *= np.sqrt(abs(wL - wph) / wL) * ram_fac return Ram From 5678f1504a68d34ce2ff77546f1bcebea2e458a3 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Mon, 18 May 2026 15:59:02 +0200 Subject: [PATCH 07/11] fixed citation --- yambopy/exciton_phonon/excph_resonant_raman.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index a05ebfea..3d7486cb 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -7,12 +7,14 @@ from yambopy.units import ha2ev -@citation("PHYSICAL REVIEW B 113, 085201 (2026)") +#@citation("PHYSICAL REVIEW B 113, 085201 (2026)") def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, elec_dipoles, eph_g, cell_vol, broad=0.1, ph_freq_threshold=5.0): """ - 1-phonon Raman tensor at independent-particle (IP) level (Stokes). + ! 1-phonon Raman tensor at independent-particle (IP) level + ! Phys. Rev. B 113, 085201 – Published 2 February, 2026 + ! DOI: 10.1103/ty8m-mgml Mirrors `compute_Raman_oneph_ip` from https://github.com/muralidhar-nalabothula/PhdScripts/blob/main/exph/raman.py @@ -20,7 +22,7 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, All band-window bookkeeping (which bands, IBZ -> BZ expansion, etc.) is the caller's responsibility. The function assumes that every input is already sliced and aligned onto the same nb-band window in the - order (valence, then conduction) along the band axis. + order (valence, then conduction) along the band axis. Parameters ---------- From c1dbd964a08399d0e52d8d79f06b71c7dbb4e09d Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Wed, 20 May 2026 15:48:53 +0200 Subject: [PATCH 08/11] added fixing laser energy to 1ph ip raman --- tutorial/exciton-phonon/ip_raman.py | 10 +- .../exciton_phonon/excph_resonant_raman.py | 319 +++++++++++++++++- 2 files changed, 321 insertions(+), 8 deletions(-) diff --git a/tutorial/exciton-phonon/ip_raman.py b/tutorial/exciton-phonon/ip_raman.py index 01f5dc6a..edc10040 100644 --- a/tutorial/exciton-phonon/ip_raman.py +++ b/tutorial/exciton-phonon/ip_raman.py @@ -7,7 +7,7 @@ them into the LetzElPhC k-grid - slices electron energies, dipoles and gkkp to the user-chosen `raman_bands` window (Fortran 1-indexed, inclusive) - - passes ready-to-use arrays to ip_resonant_raman_oneph + - passes ready-to-use arrays to ip_resonant_raman_tensor_oneph The Raman function itself takes only the energies, dipoles, gkkp, cell volume and broadening and evaluates the formula. @@ -19,7 +19,7 @@ from yambopy import (YamboLatticeDB, YamboElectronsDB, YamboDipolesDB, LetzElphElectronPhononDB) -from yambopy.exciton_phonon.excph_resonant_raman import ip_resonant_raman_oneph +from yambopy.exciton_phonon.excph_resonant_raman import ip_resonant_raman_tensor_oneph # ---- User inputs ---------------------------------------------------------- @@ -31,7 +31,7 @@ omega_range = (1.0, 6.0, 501) # laser energies (eV): min, max, npts broading = 0.10 # Lorentzian broadening (eV) ph_fre_th = 5.0 # acoustic-mode cutoff (cm^-1) -raman_bands = (47, 56) # 1-indexed Fortran inclusive +raman_bands = (44, 56) # 1-indexed Fortran inclusive # ---- Databases ------------------------------------------------------------ @@ -52,7 +52,7 @@ perm = np.argmin(np.linalg.norm(diff, axis=-1), axis=1) -# ---- Resolve band ranges (1-indexed inclusive, Fortran) ------------------- +# ---- Resolve band range (1-indexed inclusive, Fortran) ------------------- b_in_F, b_out_F = raman_bands nb = b_out_F - b_in_F + 1 n_val = int(electrons.nbandsv) - (b_in_F - 1) @@ -89,7 +89,7 @@ # ---- Raman calculation ---------------------------------------------------- laser_eV = np.linspace(omega_range[0], omega_range[1], int(omega_range[2])) -R = ip_resonant_raman_oneph( +R = ip_resonant_raman_tensor_oneph( laser_energies = laser_eV, ph_energies = ph_eV, el_energies = el_eV, diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index 3d7486cb..e7d2245e 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -4,17 +4,19 @@ import numpy as np from tqdm import tqdm +from numba import njit, prange from yambopy.units import ha2ev #@citation("PHYSICAL REVIEW B 113, 085201 (2026)") -def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, - elec_dipoles, eph_g, cell_vol, - broad=0.1, ph_freq_threshold=5.0): +def ip_resonant_raman_tensor_oneph(laser_energies, ph_energies, el_energies, + elec_dipoles, eph_g, cell_vol, + broad=0.1, ph_freq_threshold=5.0): """ ! 1-phonon Raman tensor at independent-particle (IP) level ! Phys. Rev. B 113, 085201 – Published 2 February, 2026 ! DOI: 10.1103/ty8m-mgml + Implements Eq. (D2) of the paper above. Mirrors `compute_Raman_oneph_ip` from https://github.com/muralidhar-nalabothula/PhdScripts/blob/main/exph/raman.py @@ -109,3 +111,314 @@ def ip_resonant_raman_oneph(laser_energies, ph_energies, el_energies, Ram[w_idx, m] *= np.sqrt(abs(wL - wph) / wL) * ram_fac return Ram + + +@njit(cache=True, nogil=True, parallel=True) +def _exc_raman_oneph_kernel(laser_Ha, ph_Ha, BS_energies, exc_dip_absorp, + exc_ph, ram_fac, ph_thresh_Ha): + """Numba kernel for `exc_resonant_raman_oneph`. + + Mirrors `compute_Raman_oneph_exc_numba` from Murali's PhdScripts. + All inputs in atomic units. + + Parameters + ---------- + laser_Ha : (nfreqs,) float – laser energies (Ha) + ph_Ha : (nmodes,) float – phonon energies (Ha) + BS_energies : (nexc,) complex – exciton energies minus i*Gamma/2 (Ha) + exc_dip_absorp : (npol, nexc) complex – absorption-form excitonic dipoles + (= conj of yambopy emission output) + exc_ph : (nmodes, nexc_in, nexc_out) complex + yambopy convention (phonon absorption, [mode, init, fin]). + The numba code conjugates internally to obtain the Stokes + (phonon emission) branch. + ram_fac : float – 1 / N_k / sqrt(cell_vol) + ph_thresh_Ha : float – acoustic-mode cutoff (Ha) + """ + nfreqs = len(laser_Ha) + nmodes, nexc, _ = exc_ph.shape + npol = exc_dip_absorp.shape[0] + Ram = np.zeros((nfreqs, nmodes, 3, 3), dtype=exc_dip_absorp.dtype) + + dipS_ares_base = np.conj(exc_dip_absorp) # = emission-form dipole + dipSp_res_base = np.conj(exc_dip_absorp) + + for i_ome in prange(nfreqs): + wL = laser_Ha[i_ome] + + dipS_res = exc_dip_absorp / (wL - BS_energies) + dipS_res_conj = np.conj(dipS_res) + dipS_ares = dipS_ares_base / (wL + BS_energies) + + for i in range(nmodes): + freq = ph_Ha[i] + if np.abs(freq) <= ph_thresh_Ha: + continue + + dipSp_res = dipSp_res_base / (wL - BS_energies - freq) + dipSp_ares = exc_dip_absorp / (wL + BS_energies - freq) + + ex_ph_T = exc_ph[i].T + + # Resonant term (Stokes): conj of full product gives g_emi + dipSp_res_T_conj = np.conj(dipSp_res.T) + term1 = dipS_res_conj @ ex_ph_T @ dipSp_res_T_conj + term1 = np.conj(term1) + + # Anti-resonant term + dipSp_ares_T = dipSp_ares.T + term2 = dipS_ares @ ex_ph_T @ dipSp_ares_T + + scale = np.sqrt(np.abs(wL - freq) / wL) * ram_fac + Ram[i_ome, i, :npol, :npol] = (term1 + term2) * scale + + return Ram + + +def exc_resonant_raman_oneph(laser_energies, ph_energies, exc_energies, + exc_dipoles, exc_ph_mat_el, n_kpts, cell_vol, + broad=0.1, ph_freq_threshold=5.0, + precision='d'): + """ + 1-phonon Raman tensor including excitonic effects (Stokes). + + Numba-accelerated. Mirrors `compute_Raman_oneph_exc` from + https://github.com/muralidhar-nalabothula/PhdScripts/blob/main/exph/raman.py + See also Sven Reichardt et al., Sci. Adv. 6, eabb5915 (2020). + + All band/exciton/k-point bookkeeping is the caller's responsibility. + The function assumes the inputs are already sliced consistently + (same nexc on every exciton-indexed axis) and that + `exc_dipoles` is in yambo's emission convention, + `exc_ph_mat_el` is in yambopy's raw (phonon-absorption) convention + with shape (nmodes, nexc_in, nexc_out). The kernel handles the + conjugation required to convert to phonon-emission (Stokes) Raman. + + Parameters + ---------- + laser_energies : (nfreqs,) float ndarray + Incoming laser energies in eV. + ph_energies : (nmodes,) float ndarray + Phonon energies at q=0 in eV. + exc_energies : (nexc,) float ndarray + Exciton energies at Q=0 in eV. + exc_dipoles : (3, nexc) complex ndarray + Velocity-gauge exciton dipoles in emission convention + (a.u.), as returned by yambopy `exc_dipoles_pol`. + exc_ph_mat_el : (nmodes, nexc, nexc) complex ndarray + Exciton-phonon matrix elements at q=0 in Hartree, in yambopy's + raw convention (phonon absorption, [mode, initial, final]). + This is `exciton_phonon_matelem(...)[0]` (q=0 slice) after + slicing nexc_in == nexc_out == nexc. + n_kpts : int + Number of k-points in the BZ (for the 1/N_k prefactor). + cell_vol : float + Unit-cell volume in bohr^3. + broad : float, optional + Lorentzian broadening Gamma in eV. Default: 0.1. + ph_freq_threshold : float, optional + Acoustic-mode cutoff in cm^-1. Default: 5. + precision : {'d', 's'}, optional + Floating-point precision used inside the numba kernel. + 'd' = double / complex128 (default), 's' = single / complex64. + + Returns + ------- + raman_tensor : (nfreqs, nmodes, 3, 3) complex ndarray + """ + cm1_to_Ha = 0.12398e-3 / ha2ev + + nfreqs = np.atleast_1d(laser_energies).size + nmodes = ph_energies.shape[0] + nexc = exc_energies.shape[0] + if exc_dipoles.shape[0] != 3 or exc_dipoles.shape[1] != nexc: + raise ValueError("exc_dipoles shape %s != (3, %d)" + % (tuple(exc_dipoles.shape), nexc)) + if exc_ph_mat_el.shape != (nmodes, nexc, nexc): + raise ValueError("exc_ph_mat_el shape %s != (%d, %d, %d)" + % (tuple(exc_ph_mat_el.shape), nmodes, nexc, nexc)) + if cell_vol <= 0.: + raise ValueError("cell_vol must be positive") + + prec = str(precision).strip().lower() + if prec in ('d', 'double'): + f_type, c_type = np.float64, np.complex128 + elif prec in ('s', 'single'): + f_type, c_type = np.float32, np.complex64 + else: + raise ValueError("precision must be 'd' or 's'") + + # eV -> Ha + laser_Ha = np.atleast_1d(np.asarray(laser_energies)) / ha2ev + ph_Ha = np.asarray(ph_energies) / ha2ev + broad_Ha = (broad / ha2ev) / 2.0 + ph_thresh_Ha = ph_freq_threshold * cm1_to_Ha + ram_fac = 1.0 / float(n_kpts) / np.sqrt(cell_vol) + + BS_energies = np.asarray(exc_energies) / ha2ev - 1j * broad_Ha + dip_absorp = np.conj(exc_dipoles) # yambopy emission -> absorption + + # Contiguous arrays in the chosen precision (required by the @njit kernel). + laser_c = np.ascontiguousarray(laser_Ha, dtype=f_type) + ph_c = np.ascontiguousarray(ph_Ha, dtype=f_type) + BS_c = np.ascontiguousarray(BS_energies, dtype=c_type) + dip_absorp_c = np.ascontiguousarray(dip_absorp, dtype=c_type) + exc_ph_c = np.ascontiguousarray(exc_ph_mat_el, dtype=c_type) + + return _exc_raman_oneph_kernel(laser_c, ph_c, BS_c, dip_absorp_c, + exc_ph_c, ram_fac, ph_thresh_Ha) + + +# 1 eV in cm^-1 +_EV_TO_CM1 = 8065.54 + + +def ip_raman_spectrum_oneph(laser_energy, ph_energies, + el_energies=None, elec_dipoles=None, eph_g=None, + cell_vol=None, + raman_tensor=None, + broad=0.1, ph_freq_threshold=5.0, + energy_grid=None, energy_units='cm-1', + spectrum_broad=0.1, broad_type='lorentzian', + pol_in=None, pol_out=None): + """ + 1-phonon IP Raman spectrum at a single laser energy. + + Implements Eq. (D1) of Phys. Rev. B 113, 085201 (2026): + + dσ/dΩ ∝ Σ_μν |M^{μν}(ω_L, ω_λ)|² · δ(E − ℏω_λ) + + where M^{μν}(ω_L, ω_λ) is the per-mode 3×3 IP Raman tensor (D2), + computed by `ip_resonant_raman_tensor_oneph`. The δ-function is + replaced by a Lorentzian (default) or Gaussian of FWHM + `spectrum_broad`. + + If a `raman_tensor` is given, it is used directly and no IP-Raman + tensor calculation is performed; the function then only needs + `laser_energy` (to set ω_L for the prefactor *already inside* the + tensor — i.e., we just trust it) and `ph_energies` (to know where + to put each peak). + + Parameters + ---------- + laser_energy : float + Single laser energy ω_L in eV. + ph_energies : (nmodes,) float ndarray + Phonon energies at q=0 in eV (used to place each peak on the + spectrum and, if ``raman_tensor is None``, passed through to + the tensor function). + el_energies, elec_dipoles, eph_g, cell_vol : optional + Required only when ``raman_tensor is None``. Same meaning as in + `ip_resonant_raman_tensor_oneph`. + raman_tensor : (nmodes, 3, 3) or (1, nmodes, 3, 3) complex ndarray, optional + Pre-computed Raman tensor (as returned by + `ip_resonant_raman_tensor_oneph` at a single laser energy). + If given, the tensor is reused and the IP-Raman calculation is + skipped. + broad : float, optional + Lorentzian broadening (eV) passed to the tensor function + (electronic Γ). Default: 0.1. + ph_freq_threshold : float, optional + Acoustic-mode cutoff in cm^-1 passed to the tensor function. + Default: 5. NOTE: acoustic modes are NOT removed from the + output spectrum — they appear with whatever (typically zero or + noisy) intensity comes out of the tensor calculation, so the + user can inspect them. + energy_grid : 1D ndarray, optional + Energy grid on which the spectrum is evaluated, in + `energy_units`. If None, an evenly-spaced grid is built + covering the phonon energy range with a step ≈ + `spectrum_broad/4`. + energy_units : {'cm-1', 'eV'}, optional + Units of the output `energy_grid` and of `spectrum_broad`. + Default: 'cm-1'. + spectrum_broad : float, optional + FWHM of the broadening function applied to each phonon + δ-peak, in `energy_units`. Default: 0.1. + broad_type : {'gaussian', 'lorentzian'}, optional + Shape of the broadening function. Default: 'gaussian'. + pol_in, pol_out : array-like of shape (3,), optional + Incoming / outgoing photon polarisation 3-vectors (Cartesian, + possibly complex). If BOTH are given, the per-mode intensity + is |Σ_{μν} e_out[μ] R[λ,μ,ν] e_in[ν]|² (vectors normalised + internally). If either is None (default), the unpolarised D1 + sum is used: I_λ = Σ_{μν} |R[λ,μ,ν]|². + + Returns + ------- + energy_grid : (nE,) float ndarray + Energy axis in `energy_units`. + intensity : (nE,) float ndarray + Raman intensity at each energy. + """ + # ---- Compute or accept the Raman tensor ------------------------------- + nmodes = int(ph_energies.shape[0]) + if raman_tensor is None: + for name, arr in (('el_energies', el_energies), + ('elec_dipoles', elec_dipoles), + ('eph_g', eph_g), + ('cell_vol', cell_vol)): + if arr is None: + raise ValueError("`%s` must be provided when `raman_tensor` is None" % name) + R = ip_resonant_raman_tensor_oneph( + np.array([float(laser_energy)]), + ph_energies, el_energies, elec_dipoles, eph_g, cell_vol, + broad=broad, ph_freq_threshold=ph_freq_threshold) + R = R[0] # (nmodes, 3, 3) + else: + R = np.asarray(raman_tensor) + if R.ndim == 4 and R.shape[0] == 1: + R = R[0] + if R.shape != (nmodes, 3, 3): + raise ValueError("raman_tensor shape %s != (%d, 3, 3)" % + (tuple(R.shape), nmodes)) + + # ---- Per-mode intensity ---------------------------------------------- + if (pol_in is None) or (pol_out is None): + # Eq. (D1) sum over polarisation indices + I_mode = np.sum(np.abs(R)**2, axis=(1, 2)) # (nmodes,) + else: + e_in = np.asarray(pol_in, dtype=complex) + e_out = np.asarray(pol_out, dtype=complex) + e_in = e_in / np.linalg.norm(e_in) + e_out = e_out / np.linalg.norm(e_out) + amp = np.einsum('m, lmn, n -> l', e_out, R, e_in) + I_mode = np.abs(amp)**2 # (nmodes,) + + # ---- Energy axis and broadening unit conversions --------------------- + units = str(energy_units).strip().lower() + if units in ('cm-1', 'cm^-1', 'cm', 'wavenumber'): + ph_axis = np.asarray(ph_energies) * _EV_TO_CM1 + unit_lbl = 'cm-1' + elif units in ('ev',): + ph_axis = np.asarray(ph_energies, dtype=float) + unit_lbl = 'eV' + else: + raise ValueError("energy_units must be 'cm-1' or 'eV'") + + fwhm = float(spectrum_broad) + if fwhm <= 0.: + raise ValueError("spectrum_broad (FWHM) must be > 0") + + if energy_grid is None: + emin = min(0.0, float(np.min(ph_axis)) - 5.0 * fwhm) + emax = float(np.max(ph_axis)) + 5.0 * fwhm + step = max(fwhm / 4.0, (emax - emin) / 5000.0) + energy_grid = np.arange(emin, emax + step, step) + energy_grid = np.asarray(energy_grid, dtype=float) + + # ---- Place each mode's intensity on the energy grid ------------------ + diff = energy_grid[:, None] - ph_axis[None, :] # (nE, nmodes) + btype = str(broad_type).strip().lower() + if btype == 'gaussian': + sigma = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0))) + profile = np.exp(-0.5 * (diff / sigma)**2) / (sigma * np.sqrt(2.0 * np.pi)) + elif btype == 'lorentzian': + hwhm = 0.5 * fwhm + profile = (hwhm / np.pi) / (diff**2 + hwhm**2) + else: + raise ValueError("broad_type must be 'gaussian' or 'lorentzian'") + + intensity = np.einsum('em, m -> e', profile, I_mode) + + return energy_grid, intensity From dfde3dc084edeebd2e7f97c011602a9aec1d6f5c Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Tue, 26 May 2026 16:35:25 +0200 Subject: [PATCH 09/11] added laser eng for excph_resonant_raman.py --- .../exciton_phonon/excph_resonant_raman.py | 152 +++++++++++++++++- 1 file changed, 150 insertions(+), 2 deletions(-) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index e7d2245e..4db60078 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -116,7 +116,7 @@ def ip_resonant_raman_tensor_oneph(laser_energies, ph_energies, el_energies, @njit(cache=True, nogil=True, parallel=True) def _exc_raman_oneph_kernel(laser_Ha, ph_Ha, BS_energies, exc_dip_absorp, exc_ph, ram_fac, ph_thresh_Ha): - """Numba kernel for `exc_resonant_raman_oneph`. + """Numba kernel for `exc_resonant_raman_tensor_oneph`. Mirrors `compute_Raman_oneph_exc_numba` from Murali's PhdScripts. All inputs in atomic units. @@ -175,7 +175,7 @@ def _exc_raman_oneph_kernel(laser_Ha, ph_Ha, BS_energies, exc_dip_absorp, return Ram -def exc_resonant_raman_oneph(laser_energies, ph_energies, exc_energies, +def exc_resonant_raman_tensor_oneph(laser_energies, ph_energies, exc_energies, exc_dipoles, exc_ph_mat_el, n_kpts, cell_vol, broad=0.1, ph_freq_threshold=5.0, precision='d'): @@ -422,3 +422,151 @@ def ip_raman_spectrum_oneph(laser_energy, ph_energies, intensity = np.einsum('em, m -> e', profile, I_mode) return energy_grid, intensity + + +def exc_raman_spectrum_oneph(laser_energy, ph_energies, + exc_energies=None, exc_dipoles=None, + exc_ph_mat_el=None, n_kpts=None, cell_vol=None, + raman_tensor=None, + broad=0.1, ph_freq_threshold=5.0, + energy_grid=None, energy_units='cm-1', + spectrum_broad=0.1, broad_type='lorentzian', + pol_in=None, pol_out=None, + precision='d'): + """ + 1-phonon excitonic Raman spectrum at a single laser energy (Stokes). + + Excitonic analogue of `ip_raman_spectrum_oneph`. Implements the + "place each phonon peak on the energy axis with a broadening" step + on top of the excitonic Raman tensor returned by + `exc_resonant_raman_tensor_oneph`: + + I(E) = Σ_λ |M^{μν}(ω_L, ω_λ)|² · δ(E − ℏω_λ) + + The δ-function is replaced by a Lorentzian (default) or Gaussian of + FWHM `spectrum_broad`. + + If `raman_tensor` is given, it is used directly and the excitonic + tensor calculation is skipped. + + Parameters + ---------- + laser_energy : float + Single laser energy ω_L in eV. + ph_energies : (nmodes,) float ndarray + Phonon energies at q=0 in eV (used to place each peak and, + if ``raman_tensor is None``, forwarded to the tensor function). + exc_energies, exc_dipoles, exc_ph_mat_el, n_kpts, cell_vol : optional + Required only when ``raman_tensor is None``. Same meaning as in + `exc_resonant_raman_tensor_oneph`. + raman_tensor : (nmodes, 3, 3) or (1, nmodes, 3, 3) complex ndarray, optional + Pre-computed excitonic Raman tensor at the single laser energy + of interest (as returned by `exc_resonant_raman_tensor_oneph`). + broad : float, optional + Lorentzian broadening (eV) for the electronic / excitonic Γ. + Default: 0.1. + ph_freq_threshold : float, optional + Acoustic-mode cutoff in cm^-1 passed to the tensor function. + Acoustic modes are KEPT in the output spectrum so the user can + inspect their (typically tiny / noise-level) contribution. + Default: 5. + energy_grid : 1D ndarray, optional + Energy grid on which the spectrum is evaluated, in + `energy_units`. Auto-built around the phonon range when None. + energy_units : {'cm-1', 'eV'}, optional + Units of the output `energy_grid` and of `spectrum_broad`. + Default: 'cm-1'. + spectrum_broad : float, optional + FWHM of the broadening function applied to each phonon peak, + in `energy_units`. Default: 0.1. + broad_type : {'lorentzian', 'gaussian'}, optional + Shape of the broadening function. Default: 'lorentzian'. + pol_in, pol_out : array-like of shape (3,), optional + Incoming / outgoing photon polarisation Cartesian 3-vectors + (complex allowed). If both are given, the per-mode intensity is + ``|Σ_{μν} e_out[μ] R[λ,μ,ν] e_in[ν]|²`` (vectors normalised + internally). If either is None (default), the unpolarised D1 + sum is used: ``I_λ = Σ_{μν} |R[λ,μ,ν]|²``. + precision : {'d', 's'}, optional + Forwarded to `exc_resonant_raman_tensor_oneph` when the tensor + is computed inside this function. Default: 'd'. + + Returns + ------- + energy_grid : (nE,) float ndarray + Energy axis in `energy_units`. + intensity : (nE,) float ndarray + Raman intensity at each energy. + """ + # ---- Compute or accept the Raman tensor ------------------------------- + nmodes = int(ph_energies.shape[0]) + if raman_tensor is None: + for name, arr in (('exc_energies', exc_energies), + ('exc_dipoles', exc_dipoles), + ('exc_ph_mat_el', exc_ph_mat_el), + ('n_kpts', n_kpts), + ('cell_vol', cell_vol)): + if arr is None: + raise ValueError("`%s` must be provided when `raman_tensor` is None" + % name) + R = exc_resonant_raman_tensor_oneph( + np.array([float(laser_energy)]), + ph_energies, exc_energies, exc_dipoles, exc_ph_mat_el, + n_kpts, cell_vol, + broad=broad, ph_freq_threshold=ph_freq_threshold, + precision=precision) + R = R[0] # (nmodes, 3, 3) + else: + R = np.asarray(raman_tensor) + if R.ndim == 4 and R.shape[0] == 1: + R = R[0] + if R.shape != (nmodes, 3, 3): + raise ValueError("raman_tensor shape %s != (%d, 3, 3)" % + (tuple(R.shape), nmodes)) + + # ---- Per-mode intensity ---------------------------------------------- + if (pol_in is None) or (pol_out is None): + I_mode = np.sum(np.abs(R)**2, axis=(1, 2)) # (nmodes,) + else: + e_in = np.asarray(pol_in, dtype=complex) + e_out = np.asarray(pol_out, dtype=complex) + e_in = e_in / np.linalg.norm(e_in) + e_out = e_out / np.linalg.norm(e_out) + amp = np.einsum('m, lmn, n -> l', e_out, R, e_in) + I_mode = np.abs(amp)**2 # (nmodes,) + + # ---- Energy axis and broadening unit conversions --------------------- + units = str(energy_units).strip().lower() + if units in ('cm-1', 'cm^-1', 'cm', 'wavenumber'): + ph_axis = np.asarray(ph_energies) * _EV_TO_CM1 + elif units in ('ev',): + ph_axis = np.asarray(ph_energies, dtype=float) + else: + raise ValueError("energy_units must be 'cm-1' or 'eV'") + + fwhm = float(spectrum_broad) + if fwhm <= 0.: + raise ValueError("spectrum_broad (FWHM) must be > 0") + + if energy_grid is None: + emin = min(0.0, float(np.min(ph_axis)) - 5.0 * fwhm) + emax = float(np.max(ph_axis)) + 5.0 * fwhm + step = max(fwhm / 4.0, (emax - emin) / 5000.0) + energy_grid = np.arange(emin, emax + step, step) + energy_grid = np.asarray(energy_grid, dtype=float) + + # ---- Place each mode's intensity on the energy grid ------------------ + diff = energy_grid[:, None] - ph_axis[None, :] # (nE, nmodes) + btype = str(broad_type).strip().lower() + if btype == 'gaussian': + sigma = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0))) + profile = np.exp(-0.5 * (diff / sigma)**2) / (sigma * np.sqrt(2.0 * np.pi)) + elif btype == 'lorentzian': + hwhm = 0.5 * fwhm + profile = (hwhm / np.pi) / (diff**2 + hwhm**2) + else: + raise ValueError("broad_type must be 'gaussian' or 'lorentzian'") + + intensity = np.einsum('em, m -> e', profile, I_mode) + + return energy_grid, intensity From 731a87b055d7e6338286071146138b3baaa53123 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Wed, 27 May 2026 19:47:32 +0200 Subject: [PATCH 10/11] added diagnostic function --- .../exciton_phonon/excph_resonant_raman.py | 383 ++++++++++++++++++ 1 file changed, 383 insertions(+) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index 4db60078..db31d36d 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -570,3 +570,386 @@ def exc_raman_spectrum_oneph(laser_energy, ph_energies, intensity = np.einsum('em, m -> e', profile, I_mode) return energy_grid, intensity + + +# ============================================================================= +# Diagnostic helpers — inspect every intermediate of the excitonic Raman +# tensor at a SINGLE laser energy. Pure numpy (no numba) so the code is easy +# to read; use for debugging / spectrum inspection. +# ============================================================================= + +def exc_raman_components_oneph(laser_energy, ph_energies, exc_energies, + exc_dipoles, exc_ph_mat_el, n_kpts, cell_vol, + broad=0.1, ph_freq_threshold=5.0, + save_per_pair=False): + """ + Diagnostic version of `exc_resonant_raman_tensor_oneph`. Computes the + Raman tensor at a SINGLE laser energy in pure numpy and returns every + intermediate quantity that goes into it. + + The returned dict is intended to be fed to `save_raman_components`, + which writes it both as a single .npz archive and as a folder of + column-format .dat files suitable for gnuplot. + + Parameters + ---------- + laser_energy : float + Single laser energy omega_L in eV. + ph_energies, exc_energies, exc_dipoles, exc_ph_mat_el, n_kpts, cell_vol : + Same as `exc_resonant_raman_tensor_oneph`. + broad, ph_freq_threshold : + Same as `exc_resonant_raman_tensor_oneph`. + save_per_pair : bool, optional + If True, include the full complex + (nmodes, 3, 3, nexc, nexc) per-pair amplitudes in the dict. + Memory ~ nmodes * 9 * nexc^2 * 16 bytes. Off by default. + + Returns + ------- + dict + See README produced by `save_raman_components` for the + complete list of keys. + """ + cm1_to_Ha = 0.12398e-3 / ha2ev + + # ----- Unit conversion (eV -> Ha) ------------------------------------ + wL = float(laser_energy) / ha2ev + ph_Ha = np.asarray(ph_energies, dtype=float) / ha2ev + exc_Ha = np.asarray(exc_energies, dtype=float) / ha2ev + broad_Ha = (broad / ha2ev) / 2.0 + ph_thresh_Ha = ph_freq_threshold * cm1_to_Ha + ram_fac = 1.0 / float(n_kpts) / np.sqrt(cell_vol) + + nmodes = ph_Ha.shape[0] + nexc = exc_Ha.shape[0] + + BS_energies = exc_Ha - 1j * broad_Ha # E_lam - i*Gamma/2 + + D_emi = np.asarray(exc_dipoles, dtype=complex) + D_abs = np.conj(D_emi) + + inv_denom_res_v1 = 1.0 / (wL - BS_energies) # (nexc,) + inv_denom_ares_v1 = 1.0 / (wL + BS_energies) + inv_denom_res_v2 = 1.0 / (wL - BS_energies[None, :] - ph_Ha[:, None]) # (nmodes, nexc) + inv_denom_ares_v2 = 1.0 / (wL + BS_energies[None, :] - ph_Ha[:, None]) + + dipS_res = D_abs * inv_denom_res_v1[None, :] # (3, nexc) + dipS_ares = D_emi * inv_denom_ares_v1[None, :] + dipSp_res = D_emi[None, :, :] * inv_denom_res_v2[:, None, :] # (nmodes, 3, nexc) + dipSp_ares = D_abs[None, :, :] * inv_denom_ares_v2[:, None, :] + + exc_ph_abs = np.asarray(exc_ph_mat_el, dtype=complex) # yambopy raw (mode, init, fin) + exc_ph_emi = np.conj(exc_ph_abs) # Stokes / emission + + radiation_factor = np.sqrt(np.abs(wL - ph_Ha) / wL) + active_modes = np.abs(ph_Ha) > ph_thresh_Ha + prefactor = radiation_factor * ram_fac # (nmodes,) real + + term_res = np.einsum('al, mLl, mbL -> mab', + dipS_res, exc_ph_emi, dipSp_res, optimize=True) + term_ares = np.einsum('al, mLl, mbL -> mab', + dipS_ares, exc_ph_abs, dipSp_ares, optimize=True) + + term_res = term_res * prefactor[:, None, None] + term_ares = term_ares * prefactor[:, None, None] + term_res[~active_modes] = 0.0 + term_ares[~active_modes] = 0.0 + raman_tensor = term_res + term_ares + + # Per-pair factorized intensity (sum over polarizations); cheap & always returned. + A_res = np.sum(np.abs(dipS_res)**2, axis=0) + A_ares = np.sum(np.abs(dipS_ares)**2, axis=0) + B_res = np.sum(np.abs(dipSp_res)**2, axis=1) + B_ares = np.sum(np.abs(dipSp_ares)**2, axis=1) + g_sq_pair = (np.abs(exc_ph_abs)**2).transpose(0, 2, 1) + pref_sq = (prefactor**2)[:, None, None] + + pair_intensity_res = (A_res[None, :, None] * g_sq_pair * B_res[:, None, :]) * pref_sq + pair_intensity_ares = (A_ares[None, :, None] * g_sq_pair * B_ares[:, None, :]) * pref_sq + pair_intensity_res[~active_modes] = 0.0 + pair_intensity_ares[~active_modes] = 0.0 + + components = { + 'laser_energy_eV' : float(laser_energy), + 'broad_eV' : float(broad), + 'ph_freq_threshold_cm-1' : float(ph_freq_threshold), + 'n_kpts' : int(n_kpts), + 'cell_vol_bohr3' : float(cell_vol), + 'ram_fac' : float(ram_fac), + 'exc_energies_eV' : exc_Ha * ha2ev, + 'exc_energies_Ha' : exc_Ha, + 'ph_energies_eV' : ph_Ha * ha2ev, + 'ph_energies_Ha' : ph_Ha, + 'active_modes' : active_modes.astype(np.int64), + 'exc_dip_emi' : D_emi, + 'exc_dip_absorp' : D_abs, + 'exc_ph_abs' : exc_ph_abs, + 'exc_ph_emi' : exc_ph_emi, + 'inv_denom_res_v1' : inv_denom_res_v1, + 'inv_denom_ares_v1' : inv_denom_ares_v1, + 'inv_denom_res_v2' : inv_denom_res_v2, + 'inv_denom_ares_v2' : inv_denom_ares_v2, + 'dipS_res' : dipS_res, + 'dipS_ares' : dipS_ares, + 'dipSp_res' : dipSp_res, + 'dipSp_ares' : dipSp_ares, + 'radiation_factor' : radiation_factor, + 'pair_intensity_res' : pair_intensity_res, + 'pair_intensity_ares' : pair_intensity_ares, + 'term_res' : term_res, + 'term_ares' : term_ares, + 'raman_tensor' : raman_tensor, + } + + if save_per_pair: + per_pair_res = np.einsum('al, mLl, mbL -> mablL', + dipS_res, exc_ph_emi, dipSp_res, optimize=True) + per_pair_ares = np.einsum('al, mLl, mbL -> mablL', + dipS_ares, exc_ph_abs, dipSp_ares, optimize=True) + per_pair_res = per_pair_res * prefactor[:, None, None, None, None] + per_pair_ares = per_pair_ares * prefactor[:, None, None, None, None] + per_pair_res[~active_modes] = 0.0 + per_pair_ares[~active_modes] = 0.0 + components['per_pair_res'] = per_pair_res + components['per_pair_ares'] = per_pair_ares + + return components + + +def save_raman_components(components, out_dir='raman_components', + top_n_pairs=20, heatmap_modes='auto'): + """ + Persist a diagnostic dict from `exc_raman_components_oneph` to disk. + + Writes one compressed .npz archive containing every array, plus a + set of column-format .dat files for gnuplot, plus a README and a + sample plot.gp. + + Parameters + ---------- + components : dict + As returned by `exc_raman_components_oneph`. + out_dir : str + Folder to create / overwrite. + top_n_pairs : int + Number of dominant (lambda1, lambda2) pairs listed per active mode. + heatmap_modes : {'auto', 'all', None} or list of int + Which modes get a full nexc x nexc pair-intensity heatmap file. + 'auto' = top 3 by |R|^2 (default). 'all' = every active mode. + None = skip. Explicit list = those mode indices. + """ + import os + + os.makedirs(out_dir, exist_ok=True) + + # 1. npz with everything + array_components = {} + for k, v in components.items(): + try: + array_components[k] = np.asarray(v) + except (TypeError, ValueError): + pass + np.savez_compressed(os.path.join(out_dir, 'components.npz'), + **array_components) + + wL_eV = float(components['laser_energy_eV']) + ph_eV = np.asarray(components['ph_energies_eV']) + exc_eV = np.asarray(components['exc_energies_eV']) + active = np.asarray(components['active_modes']).astype(bool) + rad_fac = np.asarray(components['radiation_factor']) + nmodes = ph_eV.shape[0] + nexc = exc_eV.shape[0] + ph_cm = ph_eV * _EV_TO_CM1 + + term_res = np.asarray(components['term_res']) + term_ares = np.asarray(components['term_ares']) + raman_t = np.asarray(components['raman_tensor']) + + mod_res_sum = np.sum(np.abs(term_res)**2, axis=(1, 2)) + mod_ares_sum = np.sum(np.abs(term_ares)**2, axis=(1, 2)) + mod_total = np.sum(np.abs(raman_t)**2, axis=(1, 2)) + + # 2. summary_modes.dat + with open(os.path.join(out_dir, 'summary_modes.dat'), 'w') as f: + f.write('# Excitonic Raman diagnostics at laser_energy = %.6f eV\n' % wL_eV) + f.write('# 1:mode_idx 2:ph_freq_cm-1 3:ph_freq_eV ' + '4:|R_res|^2_sum 5:|R_ares|^2_sum 6:|R_total|^2_sum ' + '7:radiation_factor 8:active(1)/skipped(0)\n') + for m in range(nmodes): + f.write('%6d %14.4f %14.6e %16.6e %16.6e %16.6e %12.4e %d\n' + % (m, ph_cm[m], ph_eV[m], + mod_res_sum[m], mod_ares_sum[m], mod_total[m], + rad_fac[m], int(active[m]))) + + # 3. summary_excitons.dat + D_emi = np.asarray(components['exc_dip_emi']) + inv_d_res_v1 = np.asarray(components['inv_denom_res_v1']) + inv_d_ares_v1 = np.asarray(components['inv_denom_ares_v1']) + Dx, Dy, Dz = np.abs(D_emi[0]), np.abs(D_emi[1]), np.abs(D_emi[2]) + D_sq = Dx**2 + Dy**2 + Dz**2 + + with open(os.path.join(out_dir, 'summary_excitons.dat'), 'w') as f: + f.write('# Excitonic Raman diagnostics at laser_energy = %.6f eV\n' % wL_eV) + f.write('# 1:exc_idx 2:E_eV 3:|D_x| 4:|D_y| 5:|D_z| 6:|D|^2 ' + '7:|1/denom_res_v1| 8:|1/denom_ares_v1|\n') + for l in range(nexc): + f.write('%6d %14.6f %14.6e %14.6e %14.6e %14.6e %14.6e %14.6e\n' + % (l, exc_eV[l], Dx[l], Dy[l], Dz[l], D_sq[l], + np.abs(inv_d_res_v1[l]), np.abs(inv_d_ares_v1[l]))) + + # 4. dominant pairs per mode + pair_res = np.asarray(components['pair_intensity_res']) + pair_ares = np.asarray(components['pair_intensity_ares']) + pair_total = pair_res + pair_ares + exc_ph_abs = np.asarray(components['exc_ph_abs']) + g_abs_sq = np.abs(exc_ph_abs)**2 + + pairs_dir = os.path.join(out_dir, 'pairs_per_mode') + os.makedirs(pairs_dir, exist_ok=True) + for m in range(nmodes): + if not active[m]: + continue + flat = pair_total[m].ravel() + if not np.any(flat > 0): + continue + n_top = min(top_n_pairs, flat.size) + top_idx = np.argsort(flat)[::-1][:n_top] + L1, L2 = np.unravel_index(top_idx, pair_total[m].shape) + filename = os.path.join(pairs_dir, 'pairs_mode_%03d.dat' % m) + with open(filename, 'w') as f: + f.write('# Top %d pair contributions, mode %d ' + '(omega = %.2f cm-1, %.4f eV)\n' + % (n_top, m, ph_cm[m], ph_eV[m])) + f.write('# laser_energy = %.6f eV\n' % wL_eV) + f.write('# lambda1 = first-vertex exciton; lambda2 = second-vertex exciton\n') + f.write('# 1:rank 2:lambda1 3:lambda2 4:E_lambda1_eV 5:E_lambda2_eV ' + '6:|g(lambda2,lambda1)|^2 7:|c_res|^2 8:|c_ares|^2 9:|c_total|^2\n') + for rank, (l1, l2) in enumerate(zip(L1, L2), start=1): + f.write('%6d %6d %6d %14.6f %14.6f %14.6e ' + '%14.6e %14.6e %14.6e\n' + % (rank, l1, l2, + exc_eV[l1], exc_eV[l2], + g_abs_sq[m, l2, l1], + pair_res[m, l1, l2], + pair_ares[m, l1, l2], + pair_total[m, l1, l2])) + + # 5. heatmap files + if heatmap_modes is None: + chosen = [] + elif isinstance(heatmap_modes, str): + opt = heatmap_modes.strip().lower() + if opt == 'all': + chosen = [m for m in range(nmodes) if active[m]] + elif opt == 'auto': + ranked = np.argsort(mod_total)[::-1] + chosen = [int(m) for m in ranked if active[m]][:3] + else: + raise ValueError("heatmap_modes string must be 'auto' or 'all'") + else: + chosen = [int(m) for m in heatmap_modes] + if chosen: + heatmap_dir = os.path.join(out_dir, 'heatmaps') + os.makedirs(heatmap_dir, exist_ok=True) + for m in chosen: + filename = os.path.join(heatmap_dir, 'pair_intensity_mode_%03d.dat' % m) + with open(filename, 'w') as f: + f.write('# Pair intensity heatmap, mode %d ' + '(omega = %.2f cm-1, %.4f eV)\n' + % (m, ph_cm[m], ph_eV[m])) + f.write('# laser_energy = %.6f eV\n' % wL_eV) + f.write('# gnuplot: splot "%s" u 1:2:5 with image\n' + % os.path.basename(filename)) + f.write('# 1:lambda1 2:lambda2 3:|c_res|^2 4:|c_ares|^2 5:|c_total|^2\n') + for l1 in range(nexc): + for l2 in range(nexc): + f.write('%6d %6d %14.6e %14.6e %14.6e\n' + % (l1, l2, + pair_res[m, l1, l2], + pair_ares[m, l1, l2], + pair_total[m, l1, l2])) + f.write('\n') + + # 6. README and plot.gp + with open(os.path.join(out_dir, 'README.txt'), 'w') as f: + f.write( + 'Excitonic Raman diagnostic dump\n' + 'Generated by save_raman_components(...).\n\n' + 'Laser energy : %.6f eV\n' + 'Broadening : %.4f eV (electronic Gamma, full width)\n' + 'Nmodes : %d Nexc : %d\n\n' + 'Files\n-----\n' + 'components.npz\n' + ' Compressed numpy archive with every array of the diagnostic\n' + ' dict (use np.load(...) in Python).\n\n' + 'summary_modes.dat\n' + ' One row per phonon mode. Columns:\n' + ' 1 mode_idx | 2 ph_freq (cm-1) | 3 ph_freq (eV)\n' + ' 4 |R_res|^2_sum_ab | 5 |R_ares|^2_sum_ab | 6 |R_total|^2_sum_ab\n' + ' 7 radiation_factor sqrt(|wL-wph|/wL) | 8 active(1)/skipped(0)\n\n' + 'summary_excitons.dat\n' + ' One row per exciton state. Columns:\n' + ' 1 exc_idx | 2 E (eV) | 3 |Dx| | 4 |Dy| | 5 |Dz| | 6 |D|^2\n' + ' 7 |1/(wL - E + iG/2)| | 8 |1/(wL + E - iG/2)|\n\n' + 'pairs_per_mode/pairs_mode_NNN.dat\n' + ' Top-N most contributing (lambda1, lambda2) pairs per active\n' + ' mode. Ranked by |c_res|^2 + |c_ares|^2.\n' + ' Columns: rank, lambda1, lambda2, E1, E2, |g|^2, |c_res|^2,\n' + ' |c_ares|^2, |c_total|^2\n\n' + 'heatmaps/pair_intensity_mode_NNN.dat\n' + ' Full (nexc x nexc) pair-intensity table for selected modes,\n' + ' in gnuplot matrix-with-coordinates format (blank line between\n' + ' rows). Columns: lambda1, lambda2, |c_res|^2, |c_ares|^2,\n' + ' |c_total|^2\n\n' + 'Conventions\n-----------\n' + ' lambda1 = first-vertex exciton (couples to incoming photon)\n' + ' lambda2 = second-vertex exciton (couples to outgoing photon)\n' + ' exc_ph[m, init, fin] = yambopy raw absorption-form matrix\n' + ' element. The Raman formula uses exc_ph[m, lambda2, lambda1].\n\n' + 'See plot.gp in this folder for example gnuplot commands.\n' + % (wL_eV, float(components['broad_eV']), nmodes, nexc)) + + with open(os.path.join(out_dir, 'plot.gp'), 'w') as f: + f.write( + '# Quick gnuplot recipe for the Raman diagnostic dump.\n' + '# Run inside this folder: gnuplot plot.gp\n\n' + 'set terminal pdfcairo size 9in,7in enhanced font "Helvetica,11"\n' + 'set output "diagnostic.pdf"\n' + 'set grid\n\n' + 'set multiplot layout 2,2 title ' + '"Excitonic Raman diagnostic at omega_L = %.3f eV"\n\n' + '# (1) Mode-resolved spectrum\n' + 'set title "Per-mode |R|^2 (sum over polarisations)"\n' + 'set xlabel "Raman shift (cm^{-1})"\n' + 'set ylabel "|R|^2 (arb.u.)"\n' + 'plot "summary_modes.dat" u 2:6 w impulses lw 2 t "total", \\\n' + ' "" u 2:4 w impulses lw 1 lt 3 t "resonant only"\n\n' + '# (2) Exciton oscillator strengths\n' + 'set title "Exciton |D|^2 vs E"\n' + 'set xlabel "Exciton energy (eV)"\n' + 'set ylabel "|D|^2"\n' + 'plot "summary_excitons.dat" u 2:6 w impulses lw 1.5 notitle\n\n' + '# (3) Resonance pattern\n' + 'set title "Resonance: |1/(omega_L - E +/- iGamma/2)|"\n' + 'set xlabel "Exciton energy (eV)"\n' + 'set ylabel "|1/denom|"\n' + 'plot "summary_excitons.dat" u 2:7 w l lw 1.5 t "resonant", \\\n' + ' "" u 2:8 w l lw 1 t "anti-resonant"\n\n' + '# (4) Top-pair contributions for the brightest mode\n' + 'set title "Top 20 (lambda1, lambda2) contributions"\n' + 'set xlabel "Rank"\n' + 'set ylabel "|c_total|^2"\n' + 'set logscale y\n' + 'plot for [f in system("ls pairs_per_mode/pairs_mode_*.dat | head -1")] \\\n' + ' f u 1:9 w impulses lw 2 t f\n\n' + 'unset multiplot\n' + 'unset output\n\n' + '# Pair-intensity heatmap (if heatmaps/ exists), uncomment:\n' + '# set terminal pdfcairo size 6in,5in\n' + '# set output "heatmap.pdf"\n' + '# set title "|c_total|^2(lambda1, lambda2)"\n' + '# set xlabel "lambda1"; set ylabel "lambda2"\n' + '# set palette defined (0 "white", 0.5 "orange", 1 "red")\n' + '# set logscale cb\n' + '# splot "heatmaps/pair_intensity_mode_012.dat" u 1:2:5 w image\n' + % wL_eV) From 80906521c192bc3435d3ef73ee8e28c47221d847 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Thu, 28 May 2026 15:39:38 +0200 Subject: [PATCH 11/11] added 1/energies and dip*g*dip heatmaps --- .../exciton_phonon/excph_resonant_raman.py | 153 ++++++++++++++++++ 1 file changed, 153 insertions(+) diff --git a/yambopy/exciton_phonon/excph_resonant_raman.py b/yambopy/exciton_phonon/excph_resonant_raman.py index db31d36d..c20bd4a7 100644 --- a/yambopy/exciton_phonon/excph_resonant_raman.py +++ b/yambopy/exciton_phonon/excph_resonant_raman.py @@ -953,3 +953,156 @@ def save_raman_components(components, out_dir='raman_components', '# set logscale cb\n' '# splot "heatmaps/pair_intensity_mode_012.dat" u 1:2:5 w image\n' % wL_eV) + + +def save_resonance_heatmap(components, mode_idx, out_path=None): + """ + Write a (nexc x nexc) heatmap of the bare resonance denominators for a + single phonon mode — no dipoles, no phonon matrix elements. + + Each cell (l1, l2) contains: + + res_map[l1, l2] = |1/(wL - E_l1 + iG/2)| * |1/(wL - E_l2 - wph + iG/2)| + ares_map[l1, l2] = |1/(wL + E_l1 - iG/2)| * |1/(wL + E_l2 - wph - iG/2)| + + This shows the pure double-resonance landscape: which (l1, l2) pairs are + simultaneously resonant at both photon vertices, irrespective of whether + they are optically bright or phonon-connected. + + Parameters + ---------- + components : dict + As returned by `exc_raman_components_oneph`. + mode_idx : int + Index of the phonon mode to plot. + out_path : str or None + File path for the output .dat file. If None, defaults to + 'resonance_heatmap_mode_.dat' in the current directory. + + Returns + ------- + res_map : (nexc, nexc) float ndarray + ares_map : (nexc, nexc) float ndarray + """ + import os + + inv_v1_res = np.abs(np.asarray(components['inv_denom_res_v1'])) # (nexc,) + inv_v1_ares = np.abs(np.asarray(components['inv_denom_ares_v1'])) + inv_v2_res = np.abs(np.asarray(components['inv_denom_res_v2'])) # (nmodes, nexc) + inv_v2_ares = np.abs(np.asarray(components['inv_denom_ares_v2'])) + + m = int(mode_idx) + # outer product: axis-0 = l1 (first vertex), axis-1 = l2 (second vertex) + res_map = np.outer(inv_v1_res, inv_v2_res[m]) # (nexc, nexc) + ares_map = np.outer(inv_v1_ares, inv_v2_ares[m]) + + ph_eV = np.asarray(components['ph_energies_eV']) + wL_eV = float(components['laser_energy_eV']) + exc_eV = np.asarray(components['exc_energies_eV']) + nexc = exc_eV.shape[0] + + if out_path is None: + out_path = 'resonance_heatmap_mode_%03d.dat' % m + + with open(out_path, 'w') as f: + f.write('# Bare resonance heatmap — mode %d ' + '(omega = %.2f cm-1, %.4f eV)\n' + % (m, ph_eV[m] * _EV_TO_CM1, ph_eV[m])) + f.write('# laser_energy = %.6f eV\n' % wL_eV) + f.write('# No dipoles, no phonon matrix elements.\n') + f.write('# res_map[l1,l2] = |1/(wL-E_l1+iG/2)| * |1/(wL-E_l2-wph+iG/2)|\n') + f.write('# ares_map[l1,l2] = |1/(wL+E_l1-iG/2)| * |1/(wL+E_l2-wph-iG/2)|\n') + f.write('# gnuplot: splot "%s" u 1:2:3 with image (res)\n' + % os.path.basename(out_path)) + f.write('# gnuplot: splot "%s" u 1:2:4 with image (ares)\n' + % os.path.basename(out_path)) + f.write('# 1:l1 2:l2 3:res_map 4:ares_map 5:E_l1_eV 6:E_l2_eV\n') + for l1 in range(nexc): + for l2 in range(nexc): + f.write('%6d %6d %14.6e %14.6e %14.6f %14.6f\n' + % (l1, l2, + res_map[l1, l2], ares_map[l1, l2], + exc_eV[l1], exc_eV[l2])) + f.write('\n') + + return res_map, ares_map + + +def save_dipole_g_dipole_heatmap(components, mode_idx, out_path=None): + """ + Write a (nexc x nexc) heatmap of the bare dipole-phonon-dipole product for + a single phonon mode — no energy denominators. + + Each cell (l1, l2) contains: + + dip_g_dip[l1, l2] = |D|²[l1] * |g[m, l2, l1]|² * |D|²[l2] + + where |D|²[l] = Σ_a |D_emi[a, l]|² is the total oscillator strength of + exciton l (summed over Cartesian components), and |g[m, l2, l1]|² is the + phonon matrix element squared for mode m connecting l1 -> l2. + + This shows which pairs are simultaneously optically bright and + phonon-connected, regardless of their resonance with the laser. + Comparing this with the full pair-intensity heatmap reveals how much of + the structure comes from the resonance denominators vs. the intrinsic + optical / phonon coupling. + + Parameters + ---------- + components : dict + As returned by `exc_raman_components_oneph`. + mode_idx : int + Index of the phonon mode to plot. + out_path : str or None + File path for the output .dat file. If None, defaults to + 'dip_g_dip_heatmap_mode_.dat' in the current directory. + + Returns + ------- + dip_g_dip : (nexc, nexc) float ndarray + """ + import os + + D_emi = np.asarray(components['exc_dip_emi']) # (3, nexc) + exc_ph = np.asarray(components['exc_ph_abs']) # (nmodes, nexc, nexc) [mode, init, fin] + exc_eV = np.asarray(components['exc_energies_eV']) + ph_eV = np.asarray(components['ph_energies_eV']) + wL_eV = float(components['laser_energy_eV']) + + m = int(mode_idx) + nexc = exc_eV.shape[0] + + D_sq = np.sum(np.abs(D_emi)**2, axis=0) # (nexc,) |D|² per exciton + g_sq = np.abs(exc_ph[m])**2 # (nexc, nexc) [init, fin] + # g_sq[l1, l2] = |g[m, init=l1, fin=l2]|² + # We want dip_g_dip[l1, l2] = D_sq[l1] * g_sq[l2, l1] * D_sq[l2] + # (l1 is first-vertex, l2 is second-vertex, phonon connects l1->l2 + # in the yambopy raw convention [mode, init, fin] where init=l2, fin=l1 + # matches exc_ph_emi[m, l2, l1] used in term_res) + dip_g_dip = D_sq[:, None] * g_sq.T * D_sq[None, :] # (nexc, nexc) + + if out_path is None: + out_path = 'dip_g_dip_heatmap_mode_%03d.dat' % m + + with open(out_path, 'w') as f: + f.write('# Dipole-phonon-dipole heatmap — mode %d ' + '(omega = %.2f cm-1, %.4f eV)\n' + % (m, ph_eV[m] * _EV_TO_CM1, ph_eV[m])) + f.write('# laser_energy = %.6f eV (no denominators)\n' % wL_eV) + f.write('# dip_g_dip[l1,l2] = |D|^2[l1] * |g[m,l2,l1]|^2 * |D|^2[l2]\n') + f.write('# gnuplot: splot "%s" u 1:2:3 with image\n' + % os.path.basename(out_path)) + f.write('# 1:l1 2:l2 3:dip_g_dip 4:|D|^2_l1 5:|g|^2 6:|D|^2_l2 ' + '7:E_l1_eV 8:E_l2_eV\n') + for l1 in range(nexc): + for l2 in range(nexc): + f.write('%6d %6d %14.6e %14.6e %14.6e %14.6e %14.6f %14.6f\n' + % (l1, l2, + dip_g_dip[l1, l2], + D_sq[l1], + g_sq[l2, l1], + D_sq[l2], + exc_eV[l1], exc_eV[l2])) + f.write('\n') + + return dip_g_dip