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 [21]:

# load nodes from database

# normal state DOS and band structure
dos = load_node('ec14c966-249c-4d40-989f-e45a8b918752')
bs = load_node('902f690c-f850-4216-b92a-0dbf837d7948')
bs = bs.get_outgoing(node_class=KkrCalculation).first().node

# BdG DOS
BdGdos_1K = {
    '0/2': load_node('16d774e5-774d-4d77-bf2e-34c9e399d1ca'),
    '1/2': load_node('0d15f0ca-2701-47b5-997f-ecb15805308f'),
}

BdGdos_0_1K = {
    '0/2': load_node('1f92e219-03a0-45e7-903e-c25dce4f0373'),
    '1/2': load_node('79890f8e-d430-4d08-88a5-df269be0e12d')
}

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 [5]:

def plot_dos_and_bandstruc(dos_calc, bs_calc, clim=None, ylim=None):
    plt.figure(figsize=(16,9))
    
    if not bs_calc.is_finished_ok:
        print('bandstructure calculation not finished ok (yet)')
    else:
        #plot_kkr(bs_calc, silent=True, nofig=True, noshow=True)
        #if clim is None:
        #    clim = (-6,1)
        #plt.clim(clim[0], clim[1])
        
        if clim is None:
            clim = (-6,1)# (0, 30) logscale=False, 
        plot_kkr(bs_calc.get_outgoing(node_class=KkrCalculation).first().node, silent=True, cmap='binary',
         clrbar=False, 
         noshow=True, nofig=True, ptitle='', clims = clim)
        
        if ylim is not None:
            plt.ylim(ylim)
    
    if not dos_calc.is_finished_ok:
        print('DOS calculation not finished ok (yet)')
    else:
        axBand = plt.gca()
        ylim_bands = axBand.get_ylim()

        # plot DOS on right hand side of bandstructure plot
        from mpl_toolkits.axes_grid1 import make_axes_locatable
        divider = make_axes_locatable(axBand)
        axDOS = divider.append_axes("right", 2.0, pad=0.1, sharey=axBand)
        plt.title('')
        plot_kkr(dos_calc, silent=True, nofig=True, switch_xy=True, sum_spins=True, noshow=True)
        plt.axhline(0, color='k', ls='dashed', lw=1)
        plt.title('')
        plt.ylim(ylim_bands)
        axDOS.yaxis.tick_right()
        axDOS.yaxis.set_label_position("right")

In [6]:

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):
    e = (d1[iatom, :, 0]-ef) * RY_TO_EV * 1000 # in meV
    plt.plot(e, d1[iatom,:, 1]/RY_TO_EV, color=color, label=label, **kwargs)
    e = (d2[iatom, :, 0]-ef) * RY_TO_EV * 1000 # in meV
    plt.plot(e, d2[iatom,:, 1]/RY_TO_EV, color=color, **kwargs)
    plt.ylim(0)
    plt.xlim(-7,7)
    

def BdG_dos():
    plot_BdG_dos(d1_2, d2_2, 0, label='T=0.1K', color='C0', lw=4, )
    plot_BdG_dos(d1_1, d2_1, 0, label='T=1K', color='C1', lw=4, ls='--')
    plt.axvline(-1.7, color='r', lw=3, ls=':', label='$\Delta^{\mathrm{exp.}}$', zorder=200)
    plt.axvline(1.7, color='r', lw=3, ls=':', zorder=200)
    plt.annotate(text='', xy=(-1.3, 2.2), xytext=(-0.7, 3),
                 arrowprops=dict(width=3, headlength=15, headwidth=10, color='k'))
    plt.legend(loc=2, fontsize=size_legend)
    plt.ylim(0); plt.xlim(-7, 7)
    plt.xlabel('$E-E_F$ (meV)', fontsize=size_axis_labels)
    plt.ylabel('DOS (1/eV)', fontsize=size_axis_labels)
    plt.xticks(fontsize=size_tick_labels)
    plt.yticks(fontsize=size_tick_labels)

In [81]:

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

# (a) normal state band structure
axbs = plt.subplot2grid((2, 3), (0, 0), rowspan=1, colspan=2)
plot_kkr(bs, silent=True, cmap='binary',
         clrbar=False, noshow=True, nofig=True, ptitle='', clims = [-4,3])
plt.ylim(-5,3)
plt.xticks(fontsize=size_tick_labels)
plt.yticks(fontsize=size_tick_labels)

# add legend etc.
klbl = []
labels = '$\Gamma$ H N $\Gamma$ P H|P P N'.split()
for i, ik in enumerate(bs.inputs.kpoints.labels):
    if i != 6:
        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)


# (b) normal state DOS with orbital resolution
axDOS = plt.subplot2grid((2, 3), (0, 2), rowspan=1, colspan=1)
plot_kkr(dos, silent=True, nofig=True, switch_xy=True, sum_spins=True, noshow=True, ptitle='', lw=3)
plt.axhline(0, color='k', ls='dashed', lw=1)
plt.ylabel('')
axDOS.yaxis.tick_right()
axDOS.yaxis.set_label_position("right")
plt.ylim(-5,3)
plt.xticks(fontsize=size_tick_labels)
plt.yticks(fontsize=size_tick_labels)
plt.xlim(0, 4)
axDOS.lines.remove(axDOS.lines[-2])
axDOS.lines.remove(axDOS.lines[-2])
L = plt.legend(loc=4, fontsize=size_axis_labels, numpoints=1, handlelength=0.8, ncol=2, handletextpad=0.5, labelspacing=0.1, columnspacing=0.5)
L.get_texts()[0].set_text('tot')
L.get_texts()[1].set_text('s')
L.get_texts()[2].set_text('p')
L.get_texts()[3].set_text('d')
plt.xlabel('DOS (1/eV)', fontsize=size_axis_labels)



# (c) BdG-DOS
axBdGdos = plt.subplot2grid((2, 3), (1, 0), rowspan=1, colspan=3)
BdG_dos()


plt.tight_layout()
plt.annotate(xycoords='figure fraction', xy = (0.02, 0.94), text='(a)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.85, 0.94), text='(b)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.02, 0.47), text='(c)', fontsize=size_plabel, weight="bold")


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

plt.show()

loading data

No description has been provided for this image

In [82]:

# parsed DOS files for BdG calcualtions

out = { # dos_zoom4
    'dos_data': load_node('da770436-9672-4505-ac5d-9029baf36b81'),
    'dos_data_interpol': load_node('565dce98-5514-4d3f-903d-d31b25be49ef'),
    'dos_data_eh': load_node('f68011e4-a181-4176-b5dd-2f6a3cac9c55'),
    'dos_data_interpol_eh': load_node('1a6e3d00-9f82-4b35-9e17-cdd00532254d'),
    'dos_data_he': load_node('04e5a81e-da4b-47ac-ab1e-eb4d77ad8e19'),
    'dos_data_interpol_he': load_node('62b66902-5573-4dfa-b4c4-9cf4ecd4b727'),
    'dos_data_hole': load_node('d4055aec-245d-4d8b-9f05-c5fefdc38681'),
    'dos_data_interpol_hole': load_node('7eb1b563-e073-4e33-84a6-5cd701058c2e')
}

xlbl, x, xunit = out['dos_data_interpol'].get_x(); x = x[0]
ylbl, y, yunit = out['dos_data_interpol'].get_y()[0] # total 
ylbl_eh, y_eh, yunit = out['dos_data_interpol_he'].get_y()[0] # total 
ylbl_hh, y_hh, yunit = out['dos_data_interpol_hole'].get_y()[0] # total

In [117]:

ret = load_node('da977eb2-d247-4633-af4d-49d061c58d65').outputs.retrieved
with ret.open('den_lm_ir.001.1.txt') as f:
    txt = f.readline()
    lmsize = int(txt.split()[1])
    txt = f.readline()
    irmd = int(txt.split()[-1])
    d = np.loadtxt(f)
    r = 0.52918 * d[:irmd]
    daver = d[irmd:irmd+lmsize] * 1e+5
    d = np.abs(d[irmd+lmsize:].reshape(irmd, lmsize))

In [118]:

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

plt.subplot(2,1,1)

plt.plot(x*1000, y[0], lw=4, label='normal')
plt.xlim(-10, 0.0)
plt.ylim(0)
plt.xlabel('$E-E_F$ (meV)', fontsize=size_axis_labels)
plt.ylabel('DOS (' + yunit + ')', fontsize=size_axis_labels)

plt.plot(x*1000, -1/np.pi * y_eh[0], lw=4, label='anomalous')
plt.legend(loc=2, fontsize=size_legend)

plt.xticks(fontsize=size_tick_labels)
plt.yticks(fontsize=size_tick_labels)


plt.subplot(2,1,2)

plt.bar(x=range(9), height=daver[:9])
plt.xticks(range(9), '$s$ $p_x$ $p_z$ $p_y$ $d_{xy}$ $d_{yz}$ $d_{z^2}$ $d_{xz}$ $d_{x^2-y^2}$'.split(), fontsize=size_axis_labels)
plt.title('$\overline{\chi}_L = \int \chi_L(r)\,r^2\mathrm{d}r \, / \, V$', fontsize=size_axis_labels)
plt.ylabel('$\overline{\chi}_L$ (a.u.)', fontsize=size_axis_labels)
# plt.xlabel('orbital', fontsize='xx-large')

plt.xticks(fontsize=size_axis_labels)
plt.yticks(fontsize=size_tick_labels)

im = plt.imread('t2g_eg_splitting.png')
newax = fig.add_axes([0.12,0.04,0.38,0.38], anchor='NE', zorder=1)
newax.imshow(im)
newax.axis('off')

plt.tight_layout()

plt.annotate(xycoords='figure fraction', xy = (0.02, 0.97), text='(a)', fontsize=size_plabel, weight="bold")
plt.annotate(xycoords='figure fraction', xy = (0.02, 0.47), text='(b)', fontsize=size_plabel, weight="bold")

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

plt.show()

<ipython-input-118-099429145b43>:34: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  plt.tight_layout()

Figure 2 and 3 - bulk bandstruc, DOS and BdG-DOS¶

Figure 2¶

Figure 3¶