Figure 4 band structure Nb110, DOS layer resolved and SC-DOS layer resolved¶

In [2]:

%load_ext autoreload
%autoreload 2
%matplotlib inline

from aiida import load_profile
load_profile()

from aiida.orm import load_node
import numpy as np
import matplotlib.pyplot as plt
from aiida_kkr.tools import plot_kkr
from aiida_kkr.calculations import KkrCalculation
from masci_tools.io.common_functions import interpolate_dos
from masci_tools.util.constants import RY_TO_EV

In [3]:

# load nodes from database

# DOS for unrelaxed 21 layer
dos_PBE21 = load_node('088a45d5-cbb9-4c13-b9ce-0c5b15b374ac')

# BdG DOS: EF+/-10meV, 48*4*2=384 energy points
BdGDOS_part1 = load_node('d8568999-a49d-46b9-89a2-ca54fb4a8b3b')
BdGDOS_part2 = load_node('c7bee213-6c4d-46f4-893f-aa96e0386851')

# band structure workchain
bs_calc = load_node(uuid = 'fc083491-5b28-4c29-b44c-fd5a26d33666')
bs_calc = bs_calc.get_outgoing(node_class=KkrCalculation).first().node

In [4]:

size_axis_labels = 20
size_tick_labels = 18
size_legend = 20
size_plabel = 24

# single column figure: width = 10
# double column = width=20

In [13]:

def load_data(bs_calc):
    from tqdm import tqdm
    from masci_tools.io.common_functions import search_string
    from masci_tools.util.constants import RY_TO_EV

    if 'saved_data_all_layers.npy' in bs_calc.outputs.retrieved.list_object_names():
        with bs_calc.outputs.retrieved.open('saved_data_all_layers.npy', 'rb') as _f:
            d2 = np.load(_f)
    else:
        for iatom in tqdm(range(1,28)):
            with bs_calc.outputs.retrieved.open(('qdos.'+f'{iatom:3}' + '.1.dat').replace(' ', '0'), 'r') as _f:
                tmp = np.loadtxt(_f)
            if iatom==1:
                d2 = tmp.copy()
                d_surf = tmp.copy()
            else:
                d2[:,4:] += tmp[:,4:]
            if iatom<=5 or iatom>=23:
                # 1-5, 23-27 surface and vacuum layers
                d_surf[:,4:] += tmp[:,4:]
        with bs_calc.outputs.retrieved.open('saved_data_all_layers.npy', 'wb') as _f:
            np.save(_f, d2)

    with bs_calc.outputs.retrieved.open('qdos.014.1.dat', 'r') as _f:
        d_mid = np.loadtxt(_f)
        
    e = np.sort(np.array(list(set(list(d2[:,0])))))
    ne = len(e)
    nk = len(d2)//ne

    with bs_calc.outputs.retrieved.open('output.0.txt', 'r') as _f:
        txt = _f.readlines()
    ef = float(txt[search_string('Fermi energy =', txt)].split()[-2])
    e = (e-ef)*RY_TO_EV

    d_surf = abs(d_surf[:,4]).reshape(ne, nk); d_surf = d_surf/d_surf.max()
    d_mid = abs(d_mid[:,4]).reshape(ne, nk); d_mid = d_mid/d_mid.max()
    d2 = abs(d2[:,4]).reshape(ne, nk); d_all = d2/d2.max()
        
    return d_all, d_surf, d_mid, e, nk


d_all, d_surf, d_mid, e, nk = load_data(bs_calc)

100%|██████████| 27/27 [01:08<00:00,  2.55s/it]

In [14]:

def plot_bandstruc(fig, ax, d_all, d_surf, d_mid, e, nk):
    from matplotlib.colors import LinearSegmentedColormap
    cmap_RedBlackBlue = LinearSegmentedColormap.from_list('RedBlackBlue', ['black', 'black', 'red'])

    def change_tranparency(pc, pn, fig, nmin = 0.):
        fig.canvas.draw()
        colors = pc.get_facecolors()
        n = pn.get_array()
        n = (n-n.min()) / (n.max()-n.min())
        n = (n-nmin); n[n<0] = 0.
        colors[:,3] = n
        def alpha_to_white(color):
            white = np.array([1,1,1])
            alpha = color[-1]
            color = color[:-1]
            return alpha*color + (1-alpha)*white
        colors = np.array([alpha_to_white(color) for color in colors])
        pc.set_facecolor(colors)
        fig.canvas.draw()

    n = np.log10(d_all);  n = (n-n.min()) / (n.max()-n.min())
    pn = ax.pcolormesh(range(nk), e, n, cmap='binary', shading='gouraud', rasterized=True)

    pc = ax.pcolormesh(range(nk), e, (d_surf/d_all), cmap=cmap_RedBlackBlue, shading='gouraud', rasterized=True, vmin=0, vmax=0.8)
    #plt.colorbar()

    change_tranparency(pc, pn, fig, nmin=0.5)

In [15]:

def dos_labels(ax, ylbl = '', xlbl = '', xticks = True, yticks = True, BdG=False):
    if not BdG:
        plt.xlim(-5,3)
    plt.ylabel(ylbl, fontsize=size_axis_labels)# labelpad=20)
    plt.xlabel(xlbl, fontsize=size_axis_labels)
    if not BdG:
        ax.get_legend().remove()
    if not yticks:
        ax.yaxis.set_ticklabels([])
    if not xticks:
        ax.xaxis.set_ticklabels([])
    if not BdG:
        plt.axvline(0, color='grey', ls=':', lw=2)
    else:
        plt.axvline(-1.7, color='C3', lw=2, ls='--')
        plt.axvline(1.7, color='C3', lw=2, ls='--')
    plt.xticks(fontsize=size_tick_labels)
    plt.yticks(fontsize=size_tick_labels)

In [16]:

with BdGDOS_part1.outputs.retrieved.open('complex.dos') as _f:
    ef, d1 = interpolate_dos(_f)
with BdGDOS_part2.outputs.retrieved.open('complex.dos') as _f:
    ef, d2 = interpolate_dos(_f)
    
def plot_BdG_dos(d1, d2, iatom, color='C0'):
    e = (d1[iatom, :, 0]-ef) * RY_TO_EV * 1000 # in meV
    plt.plot(e, d1[iatom,:, 1]/RY_TO_EV, color=color, lw=3)
    e = (d2[iatom, :, 0]-ef) * RY_TO_EV * 1000 # in meV
    plt.plot(e, d2[iatom,:, 1]/RY_TO_EV, color=color, lw=3)
    plt.ylim(0)
    plt.xlim(-7,7)

In [19]:

fig = plt.figure(figsize=(20, 10))

# band structure
axbs = plt.subplot2grid((3, 4), (0, 0), rowspan=3, colspan=2)
plot_bandstruc(fig, axbs, d_all, d_surf, d_mid, e, nk)
# add legend etc.
klbl = []
labels = '$\overline{S}$ $\overline{P}$ $\overline{N}$ $\overline{\Gamma}$ $\overline{H}$ $\overline{S}$ $\overline{\Gamma}$ '.split()
for i, ik in enumerate(bs_calc.inputs.kpoints.labels):
    klbl.append([int(ik[0]), labels[i]])
plt.xticks([int(i) for i in np.array(klbl)[:,0]], np.array(klbl)[:,1], fontsize=size_axis_labels)
plt.yticks([-4, -3, -2, -1, 0, 1, 2], fontsize=size_tick_labels)
plt.ylabel('$E - E_{\mathrm{F}}\,(\mathrm{eV})$', fontsize=size_axis_labels)
plt.axhline(0, ls=':', color='w')

###############################################################
### normal state DOS

axdos1 = plt.subplot2grid((3, 4), (0, 2), rowspan=1, colspan=1)
plot_kkr(dos_PBE21, l_channels=False, lw=3, iatom=3, silent=True, noshow=True, nofig=True, label='surface\nlayer', ptitle = '')
dos_labels(axdos1)

axdos2 = plt.subplot2grid((3, 4), (1, 2), rowspan=1, colspan=1)
plot_kkr(dos_PBE21, l_channels=False, lw=3, iatom=4, silent=True, noshow=True, nofig=True, label='sub-surface\nlayer', ptitle = '')
dos_labels(axdos2, ylbl = 'DOS (1/eV)')

axdos3 = plt.subplot2grid((3, 4), (2, 2), rowspan=1, colspan=1)
plot_kkr(dos_PBE21, l_channels=False, lw=3, iatom=14, silent=True, noshow=True, nofig=True, label='center\nlayer', ptitle = '')
dos_labels(axdos3, xlbl = '$E-E_{\mathrm{F}}\,(\mathrm{eV})$')

###############################################################
### BdG-DOS

axscdos1 = plt.subplot2grid((3, 4), (0, 3), rowspan=1, colspan=1)
plot_BdG_dos(d1, d2, 3)
dos_labels(axdos3, BdG=True)

axscdos2 = plt.subplot2grid((3, 4), (1, 3), rowspan=1, colspan=1)
plot_BdG_dos(d1, d2, 4)
dos_labels(axdos3, BdG=True)

axscdos3 = plt.subplot2grid((3, 4), (2, 3), rowspan=1, colspan=1)
plot_BdG_dos(d1, d2, 14)
dos_labels(axdos3, xlbl = '$E-E_{\mathrm{F}}\,(\mathrm{meV})$', BdG=True)

plt.tight_layout()

# panel labels
plt.annotate(xycoords='figure fraction', xy = (0.006, 0.95), text='(a)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.555, 0.94), text='(b)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.555, 0.62), text='(c)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.555, 0.31), text='(d)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.805, 0.94), text='(e)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.805, 0.62), text='(f)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.805, 0.31), text='(g)', fontsize=size_plabel, weight="bold")

plt.savefig('Fig4.png', dpi=150)
plt.savefig('Fig4.pdf')
plt.show()

No description has been provided for this image

In [ ]: