From 9052192310b4e44ddf0bb11a6cabbc376f2a4e55 Mon Sep 17 00:00:00 2001 From: gbrunin Date: Thu, 14 Apr 2022 14:26:45 +0200 Subject: [PATCH 1/4] Updated the e-ph plotting script --- abipy/examples/plot/advanced/plot_eph_flow.py | 79 ++++++++++++------- 1 file changed, 49 insertions(+), 30 deletions(-) diff --git a/abipy/examples/plot/advanced/plot_eph_flow.py b/abipy/examples/plot/advanced/plot_eph_flow.py index cb768c324..6d0c58373 100644 --- a/abipy/examples/plot/advanced/plot_eph_flow.py +++ b/abipy/examples/plot/advanced/plot_eph_flow.py @@ -1,14 +1,29 @@ +#!/usr/bin/env python +r""" +Alas mobility +============= + +This example shows how to plot the electronic band structure, +phonon dispersion, electron linewidths and +the convergence of the electron mobility in the case of AlAs +(see example e-ph mobility flow). +Only electrons are considered, but the script can easily be +extended to holes (and combination of e/h as well). +""" + +#%% + import numpy as np +import matplotlib as mlp import matplotlib.pyplot as plt from abipy.abilab import abiopen -import os import glob from abipy.eph.rta import RtaRobot root = "../../flows/flow_eph_mob/" -fz = 13 -mz = 2 +fz = 13 # Fontsize for text, labels,... +mz = 2 # Markersize for dots (linewidths plot) fig, axes = plt.subplots(2, 2) plt.subplots_adjust(wspace=0.35, hspace=0.33) @@ -39,7 +54,7 @@ labels = [label for tick, label in enumerate(ebands.kpoints.names) if label != None] axes[0, 0].set_xticks(ticks, [], size=fz) for tick in ticks: - axes[0, 0].plot([tick, tick], [-100, 100], '-k', linewidth=1/2, alpha=0.6) + axes[0, 0].axvline(tick, color='k', linewidth=1/2, alpha=0.6) # Set limits of the axis axes[0, 0].set_xlim([ticks[0], ticks[-1]]) @@ -53,6 +68,7 @@ ### # Open the DDB file and get phonon dispersion +# Can be used to check the phonon dispersion obtained with the e-ph code with abiopen(root + "w1/outdata/out_DDB") as abifile: phbst, phdos = abifile.anaget_phbst_and_phdos_files(ndivsm=10) @@ -64,7 +80,7 @@ labels = [label for tick, label in enumerate(phbst.qpoints.names) if label != None] axes[1, 0].set_xticks(ticks, labels, size=fz-2) for tick in ticks: - axes[1, 0].plot([tick, tick], [-100, 100], '-k', linewidth=1/2, alpha=0.6) + axes[1, 0].axvline(tick, color='k', linewidth=1/2, alpha=0.6) # Set limits of the axis axes[1, 0].set_xlim([ticks[0], ticks[-1]]) @@ -80,21 +96,32 @@ # We get the linewidths for the densest mesh computed # We open the SIGEPH file corresponding to this task with abiopen(root + 'w3/t2/outdata/out_SIGEPH.nc') as abifile: - qplist = abifile.reader.read_allqps()[0] - eigs = qplist.get_field_itemp(field="qpe", itemp=0).real - lws = 2000*qplist.get_field_itemp(field="fan0", itemp=0).imag # Factor or 2x1000 to get from linewidths to 1/tau in meV - -# The zero of the energy axis is the CBM if there are holes, + # First within SERTA + qparray = abifile.get_qp_array(mode="ks+lifetimes", rta_type="serta") + qparray = qparray[np.nonzero(qparray)] + eigs = qparray.real + lws = 2000*qparray.imag # Factor or 2x1000 to get from linewidths to 1/tau in meV + + # Then within MRTA + qparray_mrta = abifile.get_qp_array(mode="ks+lifetimes", rta_type="mrta") + qparray_mrta = qparray_mrta[np.nonzero(qparray_mrta)] + eigs_mrta = qparray_mrta.real + lws_mrta = 2000*qparray_mrta.imag + +# The zero of the energy axis is the CBM zero = np.min(eigs) # CBM -# Plot 1/tau for electrons only : much easier -axes[0, 1].plot(eigs-zero, lws, 'ok', markersize=mz) -xlabel = r'$\varepsilon_{n\mathbf{k}} - \varepsilon_{\mathrm{CBM}}$ (eV)' +# Plot 1/tau +axes[0, 1].plot(eigs-zero, lws, 'ob', markersize=mz, label='SERTA') +axes[0, 1].plot(eigs_mrta-zero, lws_mrta, 'xr', markersize=mz, label='MRTA') + axes[0, 1].set_xticks([0, 0.25], ['0', '0.25']) -# Set the axis labels -axes[0, 1].set_xlabel(xlabel, labelpad=2, fontsize=fz) +# Set the axis labels and legend +axes[0, 1].set_xlabel(r'$\varepsilon_{n\mathbf{k}} - \varepsilon_{\mathrm{CBM}}$ (eV)', labelpad=2, fontsize=fz) axes[0, 1].set_ylabel(r'$\tau_{n\mathbf{k}}^{-1}$ (meV)', fontsize=fz) +axes[0, 1].legend(loc='best', labelcolor='linecolor') + ### ### Lower right : convergence of the mobility @@ -103,24 +130,16 @@ # First we find all the RTA.nc files abifiles = glob.glob(root+'w*/t*/outdata/*RTA.nc') -# We create a robot with these files and get the k-meshes and mobilities -# Alternatively we could sort the files and get the data ourselves +# We create a robot with these files and plot the mobilities robot = RtaRobot.from_files(abifiles) -figconv = robot.plot_mobility_kconv(eh=0, bte=["ibte"], show=False) +robot.plot_mobility_kconv(eh=0, bte=["ibte", "mrta", "serta"], ax=axes[1, 1], show=False, ax_grid=False) -# xdata : ngkpt, ydata: IBTE mobility -xdata = figconv.axes[0].lines[0].get_xdata() -ydata = figconv.axes[0].lines[0].get_ydata() - -# We plot the convergence on our axes -axes[1, 1].plot(xdata, ydata, '-ok') -# Set the xticks on the ngkpt values -axes[1, 1].set_xticks(xdata) -# Set the yticks at the minimum and maximum mobilities only -axes[1, 1].set_yticks([ydata[0], ydata[-1]], [int(np.round(ydata[0])), int(np.round(ydata[-1]))]) -# Set the labels +# Tune the plot +axes[1, 1].set_title(None) +axes[1, 1].legend(fontsize=fz, framealpha=0.5) axes[1, 1].set_xlabel('$N_k \\times$ $N_k \\times$ $N_k$ $\\mathbf{k}$-point grid', fontsize=fz) -axes[1, 1].set_ylabel(r'$\mu_e^\mathrm{IBTE}$'+' (cm$^2$/(V$\\cdot$s))', fontsize=fz) +axes[1, 1].set_ylabel(r'$\mu_e$'+' (cm$^2$/(V$\\cdot$s))', fontsize=fz) + # We save the figure in pdf format fig.savefig(pretty_formula + ".pdf", bbox_inches='tight') From 209efc9e7077aa6879d651c026353dc432891aa0 Mon Sep 17 00:00:00 2001 From: gbrunin Date: Thu, 14 Apr 2022 15:31:36 +0200 Subject: [PATCH 2/4] Make use of decorators for band structures --- abipy/examples/plot/advanced/plot_eph_flow.py | 51 ++++++------------- 1 file changed, 15 insertions(+), 36 deletions(-) diff --git a/abipy/examples/plot/advanced/plot_eph_flow.py b/abipy/examples/plot/advanced/plot_eph_flow.py index 6d0c58373..9ef2ade88 100644 --- a/abipy/examples/plot/advanced/plot_eph_flow.py +++ b/abipy/examples/plot/advanced/plot_eph_flow.py @@ -10,7 +10,6 @@ Only electrons are considered, but the script can easily be extended to holes (and combination of e/h as well). """ - #%% import numpy as np @@ -46,22 +45,12 @@ pretty_formula = "".join([i for i in formula if i != str(1) and i != " "]) # Remove '1's and spaces fig.suptitle(pretty_formula) -# Plot the bands -ebands.plot_ax(axes[0,0], e0="fermie", color="black") - -# Decorate the axes with the correct ticks and labels -ticks = [tick for tick, label in enumerate(ebands.kpoints.names) if label != None] -labels = [label for tick, label in enumerate(ebands.kpoints.names) if label != None] -axes[0, 0].set_xticks(ticks, [], size=fz) -for tick in ticks: - axes[0, 0].axvline(tick, color='k', linewidth=1/2, alpha=0.6) - -# Set limits of the axis -axes[0, 0].set_xlim([ticks[0], ticks[-1]]) -axes[0, 0].set_ylim([-3, 6]) +# Plot the bands (keep only the grid on the x-axis) +ebands.plot(ax=axes[0,0], show=False, linewidth=1.5, ylims=(-3, 6), ax_grid=False) # Set label axes[0, 0].set_ylabel('Energy (eV)', fontsize=fz) +axes[0, 0].set_xlabel(None) ### ### Lower left : phonon dispersion @@ -72,22 +61,12 @@ with abiopen(root + "w1/outdata/out_DDB") as abifile: phbst, phdos = abifile.anaget_phbst_and_phdos_files(ndivsm=10) -# Plot phonon dispersion -phbst.phbands.plot_ax(axes[1,0], branch=None, units="mev", color="black") - -# Decorate the axes with the correct ticks and labels -ticks = [tick for tick, label in enumerate(phbst.qpoints.names) if label != None] -labels = [label for tick, label in enumerate(phbst.qpoints.names) if label != None] -axes[1, 0].set_xticks(ticks, labels, size=fz-2) -for tick in ticks: - axes[1, 0].axvline(tick, color='k', linewidth=1/2, alpha=0.6) - -# Set limits of the axis -axes[1, 0].set_xlim([ticks[0], ticks[-1]]) -axes[1, 0].set_ylim([0, np.max(phbst.phbands.phfreqs)*1000+1]) +# Plot phonon dispersion (keep only the grid on the x-axis) +phbst.phbands.plot(ax=axes[1,0], show=False, units="mev", linewidth=1.5, ax_grid=False) # Set label axes[1, 0].set_ylabel('Frequency (meV)', fontsize=fz) +axes[1, 0].set_xlabel(None) ### ### Upper right : linewidths @@ -99,21 +78,18 @@ # First within SERTA qparray = abifile.get_qp_array(mode="ks+lifetimes", rta_type="serta") qparray = qparray[np.nonzero(qparray)] - eigs = qparray.real + eigs = qparray.real - np.min(qparray.real) # 0 to CBM lws = 2000*qparray.imag # Factor or 2x1000 to get from linewidths to 1/tau in meV # Then within MRTA qparray_mrta = abifile.get_qp_array(mode="ks+lifetimes", rta_type="mrta") qparray_mrta = qparray_mrta[np.nonzero(qparray_mrta)] - eigs_mrta = qparray_mrta.real + eigs_mrta = qparray_mrta.real - np.min(qparray_mrta.real) # 0 to CBM lws_mrta = 2000*qparray_mrta.imag -# The zero of the energy axis is the CBM -zero = np.min(eigs) # CBM - # Plot 1/tau -axes[0, 1].plot(eigs-zero, lws, 'ob', markersize=mz, label='SERTA') -axes[0, 1].plot(eigs_mrta-zero, lws_mrta, 'xr', markersize=mz, label='MRTA') +axes[0, 1].plot(eigs, lws, 'ob', markersize=mz, label='SERTA') +axes[0, 1].plot(eigs_mrta, lws_mrta, 'xr', markersize=mz, label='MRTA') axes[0, 1].set_xticks([0, 0.25], ['0', '0.25']) @@ -132,7 +108,7 @@ # We create a robot with these files and plot the mobilities robot = RtaRobot.from_files(abifiles) -robot.plot_mobility_kconv(eh=0, bte=["ibte", "mrta", "serta"], ax=axes[1, 1], show=False, ax_grid=False) +robot.plot_mobility_kconv(ax=axes[1, 1], eh=0, bte=["ibte", "mrta", "serta"], show=False, ax_grid=False) # Tune the plot axes[1, 1].set_title(None) @@ -140,6 +116,9 @@ axes[1, 1].set_xlabel('$N_k \\times$ $N_k \\times$ $N_k$ $\\mathbf{k}$-point grid', fontsize=fz) axes[1, 1].set_ylabel(r'$\mu_e$'+' (cm$^2$/(V$\\cdot$s))', fontsize=fz) - +# Reactivate the grid on the x-axis for the band structures +axes[0, 0].grid(True, axis='x') +axes[1, 0].grid(True, axis='x') + # We save the figure in pdf format fig.savefig(pretty_formula + ".pdf", bbox_inches='tight') From 3e58a885a315e99bba0f325669a236452e547eec Mon Sep 17 00:00:00 2001 From: gbrunin Date: Fri, 1 Dec 2023 15:14:30 +0100 Subject: [PATCH 3/4] Modified EventReport.stat to allow for paths that are not found on the local machine. --- abipy/flowtk/events.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/abipy/flowtk/events.py b/abipy/flowtk/events.py index 99e2d0aaa..80b94bdc3 100644 --- a/abipy/flowtk/events.py +++ b/abipy/flowtk/events.py @@ -273,7 +273,10 @@ def __init__(self, filename: str, events=None): events: List of Event objects """ self.filename = os.path.abspath(filename) - self.stat = os.stat(self.filename) + try: + self.stat = os.stat(self.filename) + except FileNotFoundError: + self.stat = None self.start_datetime, self.end_datetime = None, None self._events = [] From 9b11d60f2f99cc4ac44da5e0c8d5a035c455020d Mon Sep 17 00:00:00 2001 From: gbrunin Date: Fri, 1 Dec 2023 15:29:41 +0100 Subject: [PATCH 4/4] Clean up of old stuff. --- abipy/examples/plot/advanced/plot_eph_flow.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/abipy/examples/plot/advanced/plot_eph_flow.py b/abipy/examples/plot/advanced/plot_eph_flow.py index 212c3070f..b9a37ca9b 100755 --- a/abipy/examples/plot/advanced/plot_eph_flow.py +++ b/abipy/examples/plot/advanced/plot_eph_flow.py @@ -4,11 +4,11 @@ ============= This example shows how to plot the electronic band structure, -phonon dispersion, electron linewidths and -the convergence of the electron mobility in the case of AlAs +phonon dispersion, electron linewidths and +the convergence of the electron mobility in the case of AlAs (see example e-ph mobility flow). Only electrons are considered, but the script can easily be -extended to holes (and combination of e/h as well). +extended to holes (and combination of e/h as well). """ #%% @@ -80,13 +80,13 @@ qparray = qparray[np.nonzero(qparray)] eigs = qparray.real - np.min(qparray.real) # 0 to CBM lws = 2000*qparray.imag # Factor or 2x1000 to get from linewidths to 1/tau in meV - + # Then within MRTA qparray_mrta = abifile.get_qp_array(mode="ks+lifetimes", rta_type="mrta") qparray_mrta = qparray_mrta[np.nonzero(qparray_mrta)] eigs_mrta = qparray_mrta.real - np.min(qparray_mrta.real) # 0 to CBM lws_mrta = 2000*qparray_mrta.imag - + # Plot 1/tau axes[0, 1].plot(eigs, lws, 'ob', markersize=mz, label='SERTA') axes[0, 1].plot(eigs_mrta, lws_mrta, 'xr', markersize=mz, label='MRTA')