From 6adc23a835e0d306813aa7c9ac5bcf40a99ad852 Mon Sep 17 00:00:00 2001 From: Fulvio Paleari Date: Sun, 25 Jan 2026 21:12:53 +0100 Subject: [PATCH 1/3] fixed swapped array dimensions --- yambopy/kpoints.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/yambopy/kpoints.py b/yambopy/kpoints.py index 2f361957..3e9059cf 100644 --- a/yambopy/kpoints.py +++ b/yambopy/kpoints.py @@ -294,7 +294,7 @@ def generate_kpoint_grid(nk1,nk2,nk3,sym_and_trev,IBZ=True,eps=1.0e-5): Nsym = len(sym_red) # Generate full regular grid in crystal coordinates - xkg = regular_grid(nk1,nk2,nk3) + xkg = regular_grid(nk1,nk2,nk3) # [nk,3] if IBZ: # Now we have to start checking for equivalent points: @@ -305,7 +305,7 @@ def generate_kpoint_grid(nk1,nk2,nk3,sym_and_trev,IBZ=True,eps=1.0e-5): if equiv[ik] == ik: for i_s in range(Nsym): # Apply symmetry operations - xkr = np.dot(sym_red[i_s,:,:], xkg[:,ik]) + xkr = np.dot(sym_red[i_s,:,:], xkg[ik,:]) # Bring back in 1st BZ xkr -= np.round(xkr) # Take opposite if symmetry is composed with TR @@ -343,8 +343,8 @@ def generate_kpoint_grid(nk1,nk2,nk3,sym_and_trev,IBZ=True,eps=1.0e-5): # Filter unique k-points unique_kpoints = (equiv == np.arange(Nk)) # Bring back unique points into first BZ - xk = xkg[:,unique_kpoints] - np.round(xkg[:,unique_kpoints]) + xk = xkg[unique_kpoints,:] - np.round(xkg[unique_kpoints,:]) wk = wkk[unique_kpoints] / np.sum(wkk[unique_kpoints]) # Normalize weights nks = len(xk[0]) - return nks, xk.T, wk + return nks, xk, wk From 94f0f705bcf99788ad0c7e08351cbdf144ff831f Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Sun, 15 Feb 2026 20:11:43 +0100 Subject: [PATCH 2/3] elphondb.py 3d plots for the gkkp added --- yambopy/dbs/elphondb.py | 163 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 162 insertions(+), 1 deletion(-) diff --git a/yambopy/dbs/elphondb.py b/yambopy/dbs/elphondb.py index 02da6a18..a12851c1 100644 --- a/yambopy/dbs/elphondb.py +++ b/yambopy/dbs/elphondb.py @@ -64,6 +64,7 @@ def __init__(self,lattice,filename='ndb.elph_gkkp',folder_gkkp='SAVE',save='SAVE self.alat = lattice.alat self.rlat = lattice.rlat self.car_kpoints = lattice.car_kpoints + self.red_kpoints = lattice.red_kpoints # Check if databases exist. Exit only if header is missing. try: database = Dataset(filename) @@ -232,7 +233,7 @@ def get_gkkp_mixed(self): self.gkkp_mixed = np.real(self.gkkp)*np.real(self.gkkp_bare)+np.imag(self.gkkp)*np.imag(self.gkkp_bare) @add_fig_kwargs - def plot_elph(self,data,kcoords=None,plt_show=False,plt_cbar=False,**kwargs): + def plot_elph_old(self,data,kcoords=None,plt_show=False,plt_cbar=False,**kwargs): """ 2D scatterplot in the BZ: @@ -280,6 +281,166 @@ def plot_elph(self,data,kcoords=None,plt_show=False,plt_cbar=False,**kwargs): if plt_show: plt.show() else: print("Plot ready.\nYou can customise adding savefig, title, labels, text, show, etc...") + + def plot_elph(self, data, kcoords=None, plt_show=False, plt_cbar=False, + threeD=False, axis='z', tol=1e-6, ncols=3, **kwargs): + """ + 2D scatterplot in the BZ: + + (i) in k-space of the quantity A_{k}(iq,inu,ib1,ib2). + (ii) in q-space of the quantity A_{q}(ik,inu,ib1,ib2). + + Any real quantity which is a function of only the k-grid or q-grid may be supplied. + + The indices iq/ik,inu,ib1,ib2 are user-specified. + + - kcoords refers to the k/q-grid in Cartesian coordinates (i.e., yelph.car_qpoints and similar). + If None is specified, a k-space, fixed-q plot is assumed. + + - if plt_show plot is shown + - if plt_cbar colorbar is shown + - kwargs example: marker='H', s=300, cmap='viridis', etc. + + NB: So far requires a 2D system. + If threeD=True, plots BZ planes at constant component along `axis` + (x/y/z) in a subplot grid. + """ + import numpy as np + import matplotlib.pyplot as plt + import math + + # --- select k-points as in original code --- + if kcoords is None: + kpts = self.car_kpoints # Assume k-space plot + else: + kpts = kcoords # Plot on momentum map supplied by user + + # --- input check as in original code --- + if len(data) != len(kpts): + raise ValueError('Something wrong in data dimensions (%d data vs %d kpts)' % (len(data), len(kpts))) + + # ------------------------------------------------------------------------- + # Helper(s) for 3D layer selection + # ------------------------------------------------------------------------- + def _axis_to_index(ax): + ax = ax.lower() + if ax == 'x': return 0 + if ax == 'y': return 1 + if ax == 'z': return 2 + raise ValueError("axis must be one of 'x', 'y', 'z'") + + def _unique_with_tol(vals, tol_): + vals = np.asarray(vals) + decimals = max(0, int(np.ceil(-np.log10(tol_)))) + vals_r = np.round(vals, decimals=decimals) + uniq = np.unique(vals_r) + return uniq, vals_r, decimals + + ax_idx = _axis_to_index(axis) + + # ------------------------------------------------------------------------- + # 2D path: keep original structure/behavior + # ------------------------------------------------------------------------- + if not threeD: + # Global plot stuff (original) + self.fig, self.ax = plt.subplots(1, 1) + self.ax.add_patch(BZ_Wigner_Seitz(self.lattice)) + + if plt_cbar: + if 'cmap' in kwargs.keys(): color_map = plt.get_cmap(kwargs['cmap']) + else: color_map = plt.get_cmap('viridis') + + lim = 1.05 * np.linalg.norm(self.rlat[0]) + self.ax.set_xlim(-lim, lim) + self.ax.set_ylim(-lim, lim) + + # Reproduce plot also in adjacent BZs (original) + BZs = shifted_grids_2D(kpts, self.rlat) + for kpts_s in BZs: + plot = self.ax.scatter(kpts_s[:, 0], kpts_s[:, 1], c=data, **kwargs) + + if plt_cbar: + self.fig.colorbar(plot) + + plt.gca().set_aspect('equal') + + if plt_show: + plt.show() + else: + print("Plot ready.\nYou can customise adding savefig, title, labels, text, show, etc...") + return + + # ------------------------------------------------------------------------- + # 3D path: make a subplot grid of constant-axis slices (projections) + # ------------------------------------------------------------------------- + uniq_layers, kpts_r, decimals = _unique_with_tol(kpts[:, ax_idx], tol) + + # which two coordinates to plot + all_idx = [0, 1, 2] + all_idx.remove(ax_idx) + xidx, yidx = all_idx[0], all_idx[1] + axis_names = {0: 'x', 1: 'y', 2: 'z'} + proj_label = f"{axis_names[xidx]}{axis_names[yidx]}" + + nlay = len(uniq_layers) + nrows = int(math.ceil(nlay / ncols)) + + self.fig, axes = plt.subplots(nrows, ncols, figsize=(5.5 * ncols, 5.5 * nrows)) + axes = np.atleast_1d(axes).ravel() + + # Color normalization across all layers (so colors mean the same everywhere) + # Matplotlib will auto-scale by default; this keeps it global if you pass vmin/vmax. + if ('vmin' not in kwargs) and ('vmax' not in kwargs): + kwargs = dict(kwargs) # avoid mutating caller dict + kwargs['vmin'] = np.nanmin(data) + kwargs['vmax'] = np.nanmax(data) + + lim = 1.05 * np.linalg.norm(self.rlat[0]) + + last_plot = None + for i, layer_val in enumerate(uniq_layers): + ax = axes[i] + + # BZ border (same as original) + ax.add_patch(BZ_Wigner_Seitz(self.lattice)) + + ax.set_xlim(-lim, lim) + ax.set_ylim(-lim, lim) + ax.set_aspect('equal') + + ax.set_title(f"{axis} ≈ {layer_val:g} (rounded {decimals} dp)") + + # select points in this layer + mask = (kpts_r == layer_val) + k_layer = kpts[mask] + d_layer = np.asarray(data)[mask] + + # reproduce plot also in adjacent BZs (needs 2D shifts on the projected plane) + BZs = shifted_grids_2D(k_layer[:, [xidx, yidx]], self.rlat) + + # shifted_grids_2D expects 2D kpts; we use only the projected coords + for kpts_s in BZs: + last_plot = ax.scatter(kpts_s[:, 0], kpts_s[:, 1], c=d_layer, **kwargs) + + ax.set_xlabel(axis_names[xidx]) + ax.set_ylabel(axis_names[yidx]) + + # Turn off unused axes + for j in range(nlay, len(axes)): + axes[j].axis('off') + + # One shared colorbar for the whole figure + if plt_cbar and last_plot is not None: + self.fig.colorbar(last_plot, ax=axes[:nlay], shrink=0.9) + + self.fig.suptitle(f"3D slices: projection on {proj_label}-plane, grouped by {axis}", y=0.995) + + self.fig.tight_layout() + + if plt_show: + plt.show() + else: + print("3D subplot grid ready.\nYou can customise adding savefig, title, labels, text, show, etc...") def __str__(self,verbose=False): From 82d726103a9d3ef26165b288d6ee4fe29a7e3b80 Mon Sep 17 00:00:00 2001 From: Petru Milev Date: Thu, 19 Feb 2026 15:25:36 +0100 Subject: [PATCH 3/3] 3d plots and (no)split ibz/bz added to elphondb and excphondb --- yambopy/dbs/elphondb.py | 5 +- yambopy/dbs/excphondb.py | 183 +++++++++++++++++++++++++++++++++++++-- 2 files changed, 179 insertions(+), 9 deletions(-) diff --git a/yambopy/dbs/elphondb.py b/yambopy/dbs/elphondb.py index a12851c1..423c84f9 100644 --- a/yambopy/dbs/elphondb.py +++ b/yambopy/dbs/elphondb.py @@ -5,7 +5,7 @@ # from yambopy import * from netCDF4 import Dataset -from math import sqrt +import math import numpy as np from yambopy.tools.string import marquee import os @@ -305,9 +305,6 @@ def plot_elph(self, data, kcoords=None, plt_show=False, plt_cbar=False, If threeD=True, plots BZ planes at constant component along `axis` (x/y/z) in a subplot grid. """ - import numpy as np - import matplotlib.pyplot as plt - import math # --- select k-points as in original code --- if kcoords is None: diff --git a/yambopy/dbs/excphondb.py b/yambopy/dbs/excphondb.py index 1c09f663..a4aba7a1 100644 --- a/yambopy/dbs/excphondb.py +++ b/yambopy/dbs/excphondb.py @@ -5,7 +5,7 @@ # from yambopy import * from netCDF4 import Dataset -from math import sqrt +import math import numpy as np from yambopy.tools.string import marquee import os @@ -68,8 +68,9 @@ def __init__(self,lattice,save_excph="./",read_all=True): self.nexc_o = database.variables['EXCITON_SUM'][1].astype(int) self.nmodes = database.variables['PHONON_MODES'][0].astype(int) self.nqpoints = database.variables['HEAD_R_LATT'][3].astype(int) - self.type_exc_i = database.variables['L_kind_in'][...].tostring().decode().strip() - self.type_exc_o = database.variables['L_kind_out'][...].tostring().decode().strip() + self.type_exc_i = database.variables['L_kind_in'][...].tobytes().decode().strip() + self.type_exc_o = database.variables['L_kind_out'][...].tobytes().decode().strip() + database.close() #Check how many databases are present @@ -83,6 +84,7 @@ def __init__(self,lattice,save_excph="./",read_all=True): self.lattice = lattice self.alat = lattice.alat self.rlat = lattice.rlat + self.red_kpoints = lattice.red_kpoints # Keep reading if read_all: self.read_full_DB() @@ -145,9 +147,179 @@ def read_excph(self): self.excph = excph_full self.excph_sq = excph_sq_full + + def plot_excph(self, data, plt_show=False, plt_cbar=False, + threeD=False, axis='z', tol=1e-6, ncols=3, title=None, split_bz = True, **kwargs): + """ + 2D scatterplot in the q-BZ of the quantity A_{iq}(ib2,ib1,inu). + + Any real quantity which is a function of only the q-grid may be supplied. + The indices iq,inu,ib1,ib2 are user-specified. + + - if plt_show plot is shown + - if plt_cbar colorbar is shown + - kwargs example: marker='H', s=300, cmap='viridis', etc. + + NB: So far requires a 2D system. + If threeD=True, plots BZ planes at constant component along `axis` + (x/y/z) in a subplot grid. + """ + + qpts = self.car_qpoints + #qpts_reciprocal = self.red_qpoints + kpts_red = self.red_kpoints + + # Input check + if len(data) != len(qpts): + raise ValueError('Something wrong in data dimensions (%d data vs %d qpts)' % (len(data), len(qpts))) + + # ------------------------------------------------------------------------- + # Helpers for 3D layer selection + # ------------------------------------------------------------------------- + def _axis_to_index(ax): + ax = ax.lower() + if ax == 'x': return 0 + if ax == 'y': return 1 + if ax == 'z': return 2 + raise ValueError("axis must be one of 'x', 'y', 'z'") + """ + def _unique_with_tol(vals, tol_): + eps = 1e-6 + vals = np.asarray(vals) + decimals = max(0, int(np.ceil(-np.log10(tol_)))) + vals_r = np.round(vals, decimals=decimals) + if not split_bz: + vals_r=vals_r-np.round(vals_r + eps) + uniq = np.unique(vals_r) + else: + uniq = np.unique(vals_r) + return uniq, vals_r, decimals + + """ + def _unique_with_tol(vals, tol_, vals_red): + eps = 1e-5 + vals = np.asarray(vals) + decimals = max(0, int(np.ceil(-np.log10(tol_)))) + vals_r = np.round(vals, decimals=decimals) + if not split_bz: + vals_red = np.asarray(vals_red) + mask = np.round( vals_red + eps) == 1 + vals_r[mask] = -vals_r[mask] + uniq = np.unique(vals_r) + else: + uniq = np.unique(vals_r) + return uniq, vals_r, decimals + ax_idx = _axis_to_index(axis) + + # ------------------------------------------------------------------------- + # 2D path: keep original structure/behavior (plus optional title) + # ------------------------------------------------------------------------- + if not threeD: + # Global plot stuff + self.fig, self.ax = plt.subplots(1, 1) + self.ax.add_patch(BZ_Wigner_Seitz(self.lattice)) + + if title is not None: + self.ax.set_title(title) + + if plt_cbar: + if 'cmap' in kwargs.keys(): color_map = plt.get_cmap(kwargs['cmap']) + else: color_map = plt.get_cmap('viridis') + + lim = 1.05 * np.linalg.norm(self.rlat[0]) + self.ax.set_xlim(-lim, lim) + self.ax.set_ylim(-lim, lim) + + # Reproduce plot also in adjacent BZs + BZs = shifted_grids_2D(qpts, self.rlat) + for qpts_s in BZs: + plot = self.ax.scatter(qpts_s[:, 0], qpts_s[:, 1], c=data, **kwargs) + + if plt_cbar: + self.cbar = self.fig.colorbar(plot) + + plt.gca().set_aspect('equal') + + if plt_show: + plt.show() + else: + print("Plot ready.\nYou can customise adding savefig, title, labels, text, show, etc...") + return + + # ------------------------------------------------------------------------- + # 3D path: make a subplot grid of constant-axis slices (projections) + # ------------------------------------------------------------------------- + uniq_layers, qpts_r, decimals = _unique_with_tol(qpts[:, ax_idx], tol, vals_red = kpts_red[:, ax_idx]) + + # which two coordinates to plot + all_idx = [0, 1, 2] + all_idx.remove(ax_idx) + xidx, yidx = all_idx[0], all_idx[1] + axis_names = {0: 'x', 1: 'y', 2: 'z'} + proj_label = f"{axis_names[xidx]}{axis_names[yidx]}" + + nlay = len(uniq_layers) + nrows = int(math.ceil(nlay / ncols)) + + self.fig, axes = plt.subplots(nrows, ncols, figsize=(5.5 * ncols, 5.5 * nrows)) + axes = np.atleast_1d(axes).ravel() + + # Global color range across all layers (unless user provided vmin/vmax) + if ('vmin' not in kwargs) and ('vmax' not in kwargs): + kwargs = dict(kwargs) # avoid mutating caller dict + kwargs['vmin'] = np.nanmin(data) + kwargs['vmax'] = np.nanmax(data) + + lim = 1.05 * np.linalg.norm(self.rlat[0]) + + last_plot = None + for i, layer_val in enumerate(uniq_layers): + ax = axes[i] + + ax.add_patch(BZ_Wigner_Seitz(self.lattice)) + ax.set_xlim(-lim, lim) + ax.set_ylim(-lim, lim) + ax.set_aspect('equal') + + ax.set_title(f"{axis} ≈ {layer_val:g} (rounded {decimals} dp)") + + mask = (qpts_r == layer_val) + q_layer = qpts[mask] + d_layer = np.asarray(data)[mask] + + # Reproduce plot also in adjacent BZs (on projected plane) + BZs = shifted_grids_2D(q_layer[:, [xidx, yidx]], self.rlat) + for qpts_s in BZs: + last_plot = ax.scatter(qpts_s[:, 0], qpts_s[:, 1], c=d_layer, **kwargs) + + ax.set_xlabel(axis_names[xidx]) + ax.set_ylabel(axis_names[yidx]) + + # Turn off unused axes + for j in range(nlay, len(axes)): + axes[j].axis('off') + + # Shared colorbar + if plt_cbar and last_plot is not None: + self.cbar = self.fig.colorbar(last_plot, ax=axes[:nlay], shrink=0.9) + + # Figure title + if title is None: + title = f"3D slices in q-BZ: projection on {proj_label}-plane, grouped by {axis}" + self.fig.suptitle(title, y=0.995) + + self.fig.tight_layout() + + if plt_show: + plt.show() + else: + print("3D subplot grid ready.\nYou can customise adding savefig, title, labels, text, show, etc...") + + + @add_fig_kwargs - def plot_excph(self,data,plt_show=False,plt_cbar=False,**kwargs): + def plot_excph_old(self,data,plt_show=False,plt_cbar=False,**kwargs): """ 2D scatterplot in the q-BZ of the quantity A_{iq}(ib2,ib1,inu). @@ -183,7 +355,8 @@ def plot_excph(self,data,plt_show=False,plt_cbar=False,**kwargs): BZs = shifted_grids_2D(qpts,self.rlat) for qpts_s in BZs: plot=self.ax.scatter(qpts_s[:,0],qpts_s[:,1],c=data,**kwargs) - if plt_cbar: self.cbar = self.fig.colorbar(plot) + #if plt_cbar: self.cbar = self.fig.colorbar(plot) + if plt_cbar: self.cbar = self.fig.colorbar(plot) plt.gca().set_aspect('equal')