%load_ext aiida
from aiida import load_profile
import aiida.orm as orm
import numpy as np
import sisl
from matplotlib.collections import LineCollection
import matplotlib.pyplot as plt
import matplotlib.transforms as mtransforms
from masci_tools.io.common_functions import interpolate_dos
from masci_tools.util.constants import RY_TO_EV
profile = 'BdG' # Modify to match name of profile archive was imported into
load_profile(profile)
(oldwidth, oldheight) = plt.rcParams['figure.figsize']
plt.rcParams['figure.figsize'] = (2*oldwidth, 2*oldheight)
plt.rcParams.update({'font.size': 32})
plt.style.use('seaborn-v0_8-deep')
axes_linewidth = 3
plot_linewidth = 3
plot_markersize = 15
NB The JuKKR data was taken from https://archive.materialscloud.org/record/2021.163 This has to be separately downloaded and imported into your AiiDA profile to generate the plot
# load nodes from database
# BdG DOS
BdGdos_1K = {
'0/2': orm.load_node('16d774e5-774d-4d77-bf2e-34c9e399d1ca'),
'1/2': orm.load_node('0d15f0ca-2701-47b5-997f-ecb15805308f'),
}
BdGdos_0_1K = {
'0/2': orm.load_node('1f92e219-03a0-45e7-903e-c25dce4f0373'),
'1/2': orm.load_node('79890f8e-d430-4d08-88a5-df269be0e12d')
}
with BdGdos_1K['0/2'].outputs.retrieved.open('complex.dos') as _f:
ef, d1_1 = interpolate_dos(_f)
with BdGdos_1K['1/2'].outputs.retrieved.open('complex.dos') as _f:
ef, d2_1 = interpolate_dos(_f)
with BdGdos_0_1K['0/2'].outputs.retrieved.open('complex.dos') as _f:
ef, d1_2 = interpolate_dos(_f)
with BdGdos_0_1K['1/2'].outputs.retrieved.open('complex.dos') as _f:
ef, d2_2 = interpolate_dos(_f)
def plot_BdG_dos(d1, d2, iatom, color='C0', label='', **kwargs):
xdata1 = (d1[iatom, :, 0]-ef) * RY_TO_EV * 1000 # in meV
ydata1 = d1[iatom,:, 1]/RY_TO_EV
xdata2 = (d2[iatom, :, 0]-ef) * RY_TO_EV * 1000 # in meV
ydata2 = d2[iatom,:, 1]/RY_TO_EV
return xdata1, ydata1, xdata2, ydata2
def BdG_dos():
T0_1K = plot_BdG_dos(d1_2, d2_2, 0, label='T=0.1K', color='C0', lw=4, )
T1K = plot_BdG_dos(d1_1, d2_1, 0, label='T=1K', color='C1', lw=4, ls='--')
return T0_1K, T1K
T0_1K, T1K = BdG_dos()
xdata = np.concatenate((T0_1K[0], T0_1K[2]))
ydata = np.concatenate((T0_1K[1], T0_1K[3]))
data = np.column_stack((xdata, ydata))
The experimental dI/dV data was extracted with permission from a PhD dissertation of Lucas Schneider
exp_data_Schneider = orm.load_node(uuid='7fbe01d5-ba2d-47cc-90d5-24281d0f4790')
data_exp = exp_data_Schneider.get_array('exp_dIdV')
scale_exp = data_exp[0,1] - np.min(data_exp[:,1])
scale_jukkr = ydata[113] - np.min(ydata)
ydata_exp = (data_exp[:,1] - np.min(data_exp[:,1])) * (scale_jukkr / scale_exp)
bulk_Nb_fixedlambda_fatbands1 = orm.load_node(uuid='89816871-fece-4a32-850f-23931fbc21eb')
bulk_Nb_fixedlambda_fatbands1_zoomed = orm.load_node(uuid='57b20f26-c9b5-43d6-bf39-8e583ce1081f')
bulk_Nb_fixedlambda_pdos_adaptive1 = orm.load_node(uuid='6404f989-edcd-4667-ad63-2afbbb492bc9')
ws = bulk_Nb_fixedlambda_fatbands1.outputs.weights.get_array('weights')
no = ws.shape[-1] // 2
ws_e = np.sum(ws[:,:,:no], axis=-1)
ws_zoomed = bulk_Nb_fixedlambda_fatbands1_zoomed.outputs.weights.get_array('weights')
no = ws_zoomed.shape[-1] // 2
ws_zoomed_e = np.sum(ws_zoomed[:,:,:no], axis=-1)
with bulk_Nb_fixedlambda_pdos_adaptive1.outputs.retrieved.open('aiida.PDOS.xml') as fh:
geom, E4x, PDOS4x = sisl.get_sile_class('PDOS.xml')(fh).read_data()
def plot_bands_aiida(node, include_bands=None, color='b', linestyle='-', marker=None, spin=0, plot_vlines=True, y_origin=None, k_offset=0, E_units='eV', figax=None, label=None, axis_linewidth=2.5, weights=None, colormap='coolwarm', colorbar_label=None, show_colorbar=True, linewidth=2, debug=False, code='SIESTA', **kwargs):
"""Plot bands stored in an AiiDA node."""
if figax is None:
fig, ax = plt.subplots()
else:
fig, ax = figax
if y_origin is None:
y_origin = 'fermi'
# Get bands data
if not y_origin.lower() == 'noshift':
if code == 'SIESTA':
Ef = node.outputs.output_parameters['E_Fermi']
else:
print(f'WARNING: code {code} not yet implemented. Not shifting bands wrt Fermi energy')
Ef = 0
else:
Ef = 0
bandsdata = node.outputs.bands._get_bandplot_data(True, y_origin=Ef, prettify_format='latex_seekpath', join_symbol='|', get_segments=True)
xdata = bandsdata['x']
ydata = bandsdata['y']
# Determine spin
spin_idx = bandsdata['band_type_idx'] == spin
ydata = ydata[:, spin_idx]
# Determine zero energy point
if y_origin.lower() == 'fermi':
ref_sub = 'F'
zero_point = 0
elif y_origin.lower() == 'vbm':
zero_point = np.max(ydata[ydata < 0])
ref_sub = 'VBM'
elif y_origin.lower() == 'cbm':
zero_point = np.min(ydata[ydata > 0])
ref_sub = 'CBM'
elif y_origin.lower() == 'noshift':
ref_sub = 'noshift'
zero_point = -Ef
else:
raise ValueError('Unrecognized value {0} for y_origin, possible values are: fermi, vbm, cbm, noshift'.format(y_origin))
# Determine conversion factor for units
conv_factor = sisl.unit.unit_convert('eV', E_units)
# Determine tick locations and plot vertical lines at segment boundaries if requested
tick_locs, tick_labels = zip(*bandsdata['labels'])
if plot_vlines:
for i in range(len(tick_locs)):
if '|' in tick_labels[i]:
p = tick_locs[i] + k_offset
ax.vlines(p, -100, +100, colors='k', linestyles='-', linewidth=axis_linewidth)
# Plot horizatonal line at zero energy point
ax.hlines(0, xdata[0], xdata[-1], colors='k', linestyles='-', linewidth=axis_linewidth)
# Determine plot style
if marker is not None:
fmt_str = f'{color}{marker}'
else:
fmt_str = f'{color}{linestyle}'
# Determine indices of bands to plot
if include_bands is None:
no = ydata.shape[1]
include_bands = np.arange(no)
# Construct the joined segments that make up the bands
segments = []
no = ydata.shape[1]
x = np.empty((0,))
y = np.empty((no,0))
iw = 0
for path in bandsdata['paths']:
if path['length'] == 1:
if np.isclose(path['x'][0], path['x'][-1]):
# This path is at a segment boundary so don't include
# Also add segment constructed so far to segments
if not x.size == 0:
if weights is not None:
segments.append((x+k_offset, y, weights[iw:(iw+x.size), :]))
iw += x.size
else:
segments.append((x+k_offset, y))
x = np.empty((0,))
y = np.empty((no,0))
continue
if x.size == 0:
x = np.concatenate((x, path['x']))
y = np.hstack((y, conv_factor*(np.array(path['values'])[spin_idx,:] - zero_point)))
else:
x = np.concatenate((x, path['x'][1:]))
y = np.hstack((y, conv_factor*(np.array(path['values'])[spin_idx,1:] - zero_point)))
if not x.size == 0:
# Add last segment to list
if weights is not None:
segments.append((x+k_offset, y, weights[iw:(iw+x.size), :]))
iw += x.size
else:
segments.append((x+k_offset, y))
# Plot the bands
first_line = True
if weights is not None:
# Plot bands with color that corresponds to supplied weights
norm = plt.Normalize(weights[:,include_bands].min(), weights[:,include_bands].max())
for ks, Es, ws in segments:
for iband in include_bands:
points = np.array([ks, Es[iband,:]]).T.reshape(-1, 1, 2)
ss = np.concatenate([points[:-1], points[1:]], axis=1)
lc = LineCollection(ss, cmap=colormap, norm=norm)
lc.set_array(ws[:,iband])
lc.set_linewidth(linewidth)
line = ax.add_collection(lc)
if first_line:
if show_colorbar:
cbar = fig.colorbar(line, label=colorbar_label)
else:
cbar = None
first_line = False
else:
# Plot regular bands
for ks, Es in segments:
for iband in include_bands:
if first_line:
ax.plot(ks, Es[iband,:], fmt_str, label=label, linewidth=linewidth, **kwargs)
first_line = False
else:
ax.plot(ks, Es[iband,:], fmt_str, linewidth=linewidth, **kwargs)
#if weights is not None:
# fig.colorbar(line, ax=ax)
# Do some plot layout
if ref_sub == 'noshift':
ax.set_ylabel(r'$E$ [{0}]'.format(E_units))
else:
ax.set_ylabel(r'$E-E_\mathrm{{{0}}}$ [{1}]'.format(ref_sub, E_units))
ax.set_xticks(ticks=np.array(tick_locs)+k_offset, labels=tick_labels)
ax.set_xlim(xdata[0]+k_offset, xdata[-1]+k_offset)
ax.set_ylim(-3, 3)
for spine in ax.spines.values():
spine.set_linewidth(axis_linewidth)
if debug:
return fig, ax, segments
else:
if weights is not None:
return fig, ax, cbar
else:
return fig, ax
iN = bulk_Nb_fixedlambda_fatbands1.inputs.bandskpoints.labels[2][0]
kN = bulk_Nb_fixedlambda_fatbands1.inputs.bandskpoints.get_kpoints(cartesian=True)[iN,:]
kN_dist = bulk_Nb_fixedlambda_fatbands1.outputs.bands._get_bandplot_data(True, prettify_format='latex_seekpath', join_symbol='|', get_segments=True)['labels'][2][0]
k_zoom_start = bulk_Nb_fixedlambda_fatbands1_zoomed.inputs.bandskpoints.get_kpoints(cartesian=True)[0,:]
k_zoom_end = bulk_Nb_fixedlambda_fatbands1_zoomed.inputs.bandskpoints.get_kpoints(cartesian=True)[-1,:]
offset = kN_dist + np.linalg.norm(k_zoom_start - kN)
k_zoom_start_dist = offset
k_zoom_end_dist = offset + np.linalg.norm(k_zoom_end - k_zoom_start)
(width, height) = plt.rcParams['figure.figsize']
fig, axs = plt.subplot_mosaic([['a)'], ['b)']],
layout='constrained', figsize=(width+2, 2*height))
axs['a)'].plot(1000*E4x, np.sum(np.sum(PDOS4x[:2,:,:], axis=0), axis=0), 'b-', label=r'SIESTA', linewidth=plot_linewidth)
axs['a)'].plot(xdata, ydata, 'r-', label='JuKKR', linewidth=plot_linewidth)
axs['a)'].plot(data_exp[:,0], ydata_exp, 'g.', label='Exp', markersize=plot_markersize)
axs['a)'].legend()
axs['a)'].set_xlabel(r'$E - E_F$ [meV]')
axs['a)'].set_ylabel(r'DOS [1/eV]')
axs['a)'].set_xlim((-4,4))
for axis in ['top','bottom','left','right']:
axs['a)'].spines[axis].set_linewidth(axes_linewidth)
_, _, cbar = plot_bands_aiida(bulk_Nb_fixedlambda_fatbands1, figax=(fig, axs['b)']), weights=ws_e, y_origin='noshift', axis_linewidth=axes_linewidth, linewidth=plot_linewidth, colormap='coolwarm')
axs['b)'].set_ylim((-1, 1))
axs['b)'].set_ylabel(r'$E - E_F$ [eV]')
cbar.set_ticks(ticks=[0, 1], labels=['h', 'e'])
cbar.set_label(r'$\sum_{\mu \nu \alpha} u_{i\nu}^{*\alpha} u_{i\mu}^{\alpha} S_{\nu\mu}$')
dk = k_zoom_end_dist - k_zoom_start_dist
xlim = [k_zoom_start_dist-dk/2, k_zoom_end_dist+dk/2]
ylim = [-0.02, 0.02]
axins = axs['b)'].inset_axes(
[0.65, 0.65, 0.34, 0.34], # [x, y, width, height] w.r.t. axes
xlim=xlim, ylim=ylim, # sets viewport & tells relation to main axes
)
_, _, _ = plot_bands_aiida(bulk_Nb_fixedlambda_fatbands1_zoomed, figax=(fig, axins), weights=ws_zoomed_e, y_origin='noshift', k_offset=offset, axis_linewidth=axes_linewidth, linewidth=plot_linewidth, colormap='coolwarm', show_colorbar=False)
axins.set_xticks([])
axins.set_yticks([])
axins.set_ylabel('')
axins.set_xlim(xlim)
axins.set_ylim(ylim)
axs['b)'].indicate_inset_zoom(axins, edgecolor="black", linewidth=axes_linewidth-1, alpha=1)
for label, ax in axs.items():
# label physical distance to the left and up:
trans = mtransforms.ScaledTranslation(-50/72, 20/72, fig.dpi_scale_trans)
ax.text(0.0, 1.0, label, transform=ax.transAxes + trans,
fontsize='large', va='bottom')
#fig.savefig('bulkNB-DOS-fatbands.pdf', bbox_inches='tight')
plt.show()