Run all cells to generate the figures used in the paper.
In [ ]:
import numpy as np
import holoviews as hv
import holoviews_rc
import kwant
from fun import *
import os
def ticks(plot, x=True, y=True):
hooks = [tick_marks]
if x:
xticks = [0, 1, 2]
else:
xticks = [(0,''), (1,''), (2,'')]
hooks.append(hide_x)
if y:
yticks = [0, 17, 35]
else:
yticks = [(0, ''), (17, ''), (35, '')]
hooks.append(hide_y)
return plot(plot={'Image': {'xticks': xticks, 'yticks': yticks},
'Overlay': {'final_hooks': hooks}})
def tick_marks(plot, element):
ax = plot.handles['axis']
fig = plot.state
ax.tick_params(which='major', color='k', size=3)
def hide_x(plot, element):
ax = plot.handles['axis']
ax.set_xlabel('')
def hide_y(plot, element):
ax = plot.handles['axis']
ax.set_ylabel('')
hv.notebook_extension()
%output size=100 dpi=250 css={'width': '3.4in'}
renderer = hv.Store.renderers['matplotlib'].instance(fig='pdf', size=100, dpi=250)
from holoviews.plotting.mpl import MPLPlot
MPLPlot.fig_inches = (3.4, None)
In [ ]:
import matplotlib.cm
import matplotlib.colors as mcolors
colors1 = matplotlib.cm.binary_r(np.linspace(0.5, 1, 128))
colors2 = matplotlib.cm.gist_heat_r(np.linspace(0, 0.8, 127))
colors = np.vstack((colors1, colors2))
mymap = mcolors.LinearSegmentedColormap.from_list('my_colormap', colors)
sc_on_side_alpha100 = create_holoviews('data/0_to_2T_1x1_angles_sc_on_side_mu_ranging_from_minus_2_to_plus_2_full_phase_diagram_with_correction_A_alpha100.h5')
sc_on_side_no_orb_alpha100 = create_holoviews('data/0_to_2T_1x1_angles_sc_on_side_mu_ranging_from_minus_2_to_plus_2_full_phase_diagram_with_correction_A_no_orbital_alpha100.h5')
sc_on_side = create_holoviews('data/0_to_2T_1x1_angles_sc_on_side_mu_ranging_from_minus_2_to_plus_2_full_phase_diagram_with_correction_A.h5')
sc_on_side_no_orb = create_holoviews('data/0_to_2T_1x1_angles_sc_on_side_mu_ranging_from_minus_2_to_plus_2_full_phase_diagram_with_correction_A_no_orbital.h5')
In [ ]:
%%opts Layout [vspace=0] Image (cmap=mymap clims=(-197, 197))
%%opts Layout [sublabel_position=(-0.4, 0.9) sublabel_format='({alpha})' sublabel_size=13]
%%opts Path (color='g')
im1 = sc_on_side_no_orb.Phase_diagram.Band_gap[0.5, 0]
im2 = sc_on_side.Phase_diagram.Band_gap[0.5, 0]
im1 = im1.relabel(r"$\bm{B} \parallel x, \; \bm{A} = 0$", depth=1)
im2 = im2.relabel(r"$\bm{B} \parallel x, \; \bm{A} \ne 0$", depth=1)
max1 = np.nanmax(im1.Im.Band_gap.data)
max2 = np.nanmax(im2.Im.Band_gap.data)
max_gap = np.max((max1, max2))
sc_on_side_hist = (ticks(im1, x=False).hist(bin_range=(0, max_gap)) +
ticks(im2).hist(bin_range=(0, max_gap)))
sc_on_side_hist.cols(1)
In [ ]:
# print the maximum band gaps
print("""The maximum band gap of the top plot is {:.4} meV.
The maximum band gap of the lower plot is {:.4} meV.""".format(max1, max2))
In [ ]:
%%opts Layout [vspace=0] Image (cmap=mymap clims=(-197, 197))
%%opts Layout [sublabel_position=(-0.4, 0.9) sublabel_format='({alpha})' sublabel_size=13]
%%opts Path (color='g')
im1_alpha100 = sc_on_side_no_orb_alpha100.Phase_diagram.Band_gap[0.5, 0]
im2_alpha100 = sc_on_side_alpha100.Phase_diagram.Band_gap[0.5, 0]
im1_alpha100 = im1_alpha100.relabel(r"$\bm{B} \parallel x, \; \bm{A} = 0$", depth=1)
im2_alpha100 = im2_alpha100.relabel(r"$\bm{B} \parallel x, \; \bm{A} \ne 0$", depth=1)
max1_alpha100 = np.nanmax(im1_alpha100.Im.Band_gap.data)
max2_alpha100 = np.nanmax(im2_alpha100.Im.Band_gap.data)
max_gap_alpha100 = np.max((max1_alpha100, max2_alpha100))
sc_on_side_hist_alpha100 = (ticks(im1_alpha100, x=False).hist(bin_range=(0, max_gap_alpha100)) +
ticks(im2_alpha100).hist(bin_range=(0, max_gap_alpha100)))
(sc_on_side_hist_alpha100).cols(1)
In [ ]:
# renderer.save(sc_on_side_hist, 'paper/figures/sc_on_side_hist', fmt='pdf')
In [ ]:
# print the maximum band gaps
print("""The maximum band gap of the top plot is {:.4} meV.
The maximum band gap of the lower plot is {:.4} meV.""".format(max1_alpha100, max2_alpha100))
In [ ]:
%%opts Layout [vspace=0] Image (clims=(0, 1.5))
%%opts Layout [sublabel_position=(-0.4, 0.9) sublabel_format='({alpha})' sublabel_size=13]
%%opts Path (color='g')
im1 = sc_on_side_no_orb.Phase_diagram.Inverse_decay_length[0.5, 0]
im2 = sc_on_side.Phase_diagram.Inverse_decay_length[0.5, 0]
im1 = im1.relabel(r"$\bm{B} \parallel x, \; \bm{A} = 0$", depth=1)
im2 = im2.relabel(r"$\bm{B} \parallel x, \; \bm{A} \ne 0$", depth=1)
dat1 = im1.Im.Inverse_decay_length.data
dat2 = im2.Im.Inverse_decay_length.data
dat1[dat1<0] = np.nan
dat2[dat2<0] = np.nan
sc_on_side_length = (ticks(im1, x=False).hist(bin_range=(0, 1)) +
ticks(im2).hist(bin_range=(0, 1)))
sc_on_side_length.cols(1)
In [ ]:
%%opts Layout [vspace=0] Image (clims=(0, 1.5))
%%opts Layout [sublabel_position=(-0.4, 0.9) sublabel_format='({alpha})' sublabel_size=13]
%%opts Path (color='g')
im1_alpha100 = sc_on_side_no_orb_alpha100.Phase_diagram.Inverse_decay_length[0.5, 0]
im2_alpha100 = sc_on_side_alpha100.Phase_diagram.Inverse_decay_length[0.5, 0]
im1_alpha100 = im1_alpha100.relabel(r"$\bm{B} \parallel x, \; \bm{A} = 0$", depth=1)
im2_alpha100 = im2_alpha100.relabel(r"$\bm{B} \parallel x, \; \bm{A} \ne 0$", depth=1)
dat1_alpha100 = im1_alpha100.Im.Inverse_decay_length.data
dat2_alpha100 = im2_alpha100.Im.Inverse_decay_length.data
dat1_alpha100[dat1_alpha100<0] = np.nan
dat2_alpha100[dat2_alpha100<0] = np.nan
sc_on_side_length = (ticks(im1_alpha100, x=False).hist(bin_range=(0, 1)) +
ticks(im2_alpha100).hist(bin_range=(0, 1)))
sc_on_side_length.cols(1)
In [ ]:
# renderer.save(sc_on_side_length, 'paper/figures/sc_on_side_length', fmt='pdf')
In [ ]:
# print the minimum decay lengths in nm
print("""The minimum decay length of the top plot is {:.3} nm.
The minimum decay length of the lower plot is {:.3} nm.""".format(1000 / np.nanmax(dat1), 1000 / np.nanmax(dat2)))
In [ ]:
# print the mode of the decay lengths
frequencies, edges = np.histogram(dat1[dat1>0].reshape(-1), bins=400)
max_mode1 = edges[np.argmax(frequencies)]
frequencies, edges = np.histogram(dat2[dat2>0].reshape(-1), bins=400)
max_mode2 = edges[np.argmax(frequencies)]
print("""The maximum mode of the top plot is {:.2} µm^-1.
The maximum mode of the lower plot is {:.2} µm^-1.
The ratio is {:.3}""".format(max_mode1, max_mode2, max_mode1 / max_mode2))
In [ ]:
p = make_params(mu=4.8, orbital=True, V=lambda x,y,z: 2/50 * z, t_interface=7*constants.t/8,
Delta=5, alpha=50, A_correction=False)
momenta = np.linspace(-0.6, 0.6, 200)
def bands(B):
p.B_x, p.B_y, p.B_z = B
bands_fun = kwant.physics.Bands(lead, args=[p])
_bands = np.array([bands_fun(k=k) for k in momenta])
return hv.Path((momenta, _bands), kdims=[r'$k$', r'$E$'])
E = (-1.5, 1.5)
k = (-0.65, 0.65)
lead = make_3d_wire_external_sc(a=constants.a, angle=0)
x1 = bands((0.5, 0, 0)).select(E=E, k=k)
y1 = bands((0, 0.5, 0)).select(E=E, k=k)
z1 = bands((0, 0, 0.5)).select(E=E, k=k)
lead = make_3d_wire_external_sc(a=constants.a)
x2 = bands((0.5, 0, 0)).select(E=E, k=k)
y2 = bands((0, 0.5, 0)).select(E=E, k=k)
z2 = bands((0, 0, 0.5)).select(E=E, k=k)
In [ ]:
%%output fig='svg'
%%opts Layout [vspace=0.1 hspace=0.1 sublabel_format='']
%%opts Path (color='k')
def labels(plot, x=False, y=False, label=''):
hooks = [tick_marks]
if not x:
hooks.append(hide_x)
if not y:
hooks.append(hide_y)
plot *= hv.HLine(0)(style=dict(lw=0.5, color='k', ls=(1, (3.0, 3.0))))
return plot.relabel(label)(plot={'Path': {'xticks': 0, 'yticks': 0},
'Overlay': {'final_hooks': hooks}})
opts = {'x': -0.62, 'y': 1.40, 'fontsize': 10, 'valign':'top', 'halign':'left'}
def rectangle(x=opts['x'], y=opts['y']-0.38, width=0.55, height=0.47):
box = np.array([(x,y), (x+width, y), (x+width, y+height), (x, y+height)])
return hv.Polygons([box])(style={'facecolor': '#F0F0F0'})
box2 = rectangle(width=0.55)
box3 = rectangle(width=0.80)
x1_txt = hv.Text(text="$\mathcal{P}$, $\mathcal{R}_x$, $\mathcal{C}'$", **opts) * box3
y1_txt = hv.Text(text="$\mathcal{P}$", **opts)
z1_txt = hv.Text(text="$\mathcal{P}$, $\mathcal{C}'$", **opts) * box2
x2_txt = hv.Text(text="$\mathcal{P}$, $\mathcal{R}_x$", **opts) * box2
y2_txt = hv.Text(text="$\mathcal{P}$", **opts)
z2_txt = hv.Text(text="$\mathcal{P}$", **opts)
gap_line = lambda x: hv.HLine(np.abs(np.array(x.data)[:, :, 1]).min())(style=dict(lw='0.5', c='r', ls=(1., (3., 3.))))
bands_layout = (labels(x1 * x1_txt * gap_line(x1), label=r"$\bm{B}\parallel \hat{x}$", y=True)+
labels((y1 * y1_txt),label=r"$\bm{B}\parallel \hat{y}$") +
labels((z1 * z1_txt * gap_line(z1)), label=r"$\bm{B}\parallel \hat{z}$") +
labels((x2 * x2_txt * gap_line(x2)), x=True, y=True) +
labels((y2 * y2_txt), x=True) +
labels((z2 * z2_txt), x=True)).cols(3)
bands_layout
In [ ]:
# renderer.save(bands_layout, 'paper/figures/bandstructure_annotated', fmt='pdf')
In [ ]:
orb = create_holoviews('data/0_to_2T_4x4_angles_misaligned_with_electric_field.h5')
no_orb = create_holoviews('data/0_to_2T_4x4_angles_misaligned_no_orbital_with_electric_field.h5')
In [ ]:
%%opts Path (color='g')
%%opts Image.d [colorbar=True cbar_ticks=np.linspace(0, 140, 5).tolist()]
%%opts Layout [vspace=0.20 hspace=0.15 sublabel_position=(-0.07, 0.79) sublabel_size=10 sublabel_format='({alpha})']
%%opts VLine (linewidth=0.5 color='k')
test = orb.Phase_diagram.Band_gap[0, 0.5]
comparing_phase_diagrams = (
ticks((no_orb.Phase_diagram.Band_gap * hv.VLine(1)).relabel(r"$\bm{B} \parallel \hat{x}, \; \bm{A} = 0$")[0.5, 0], x=False)
+ ticks(no_orb.Phase_diagram.Band_gap.relabel(label=r"$\bm{B} \parallel \hat{z}, \; \bm{A} = 0$")[0, 0.5], x=False, y=False)
+ ticks(orb.Phase_diagram.Band_gap.relabel(r"$\bm{B} \parallel \hat{x}, \; \bm{A} \ne 0$")[0.5, 0])
+ ticks(orb.Phase_diagram.Band_gap.relabel(label=r"$\bm{B} \parallel \hat{z}, \; \bm{A} \ne 0$", group='d', depth=2)[0, 0.5], y=False)).cols(2)
comparing_phase_diagrams
In [ ]:
# renderer.save(comparing_phase_diagrams, 'paper/figures/comparing_phase_diagrams', fmt='pdf')
In [ ]:
%%opts Path (color='g')
%%opts Image.d [colorbar=True cbar_ticks=np.linspace(0, 120, 5).tolist()]
%%opts Layout [vspace=0.20 hspace=0.15 sublabel_position=(-0.07, 0.79) sublabel_size=10 sublabel_format='({alpha})']
kys = no_orb.Phase_diagram.Band_gap.keys()
test = orb.Phase_diagram.Band_gap[nearest(kys, 0.05), 0.5]
misaligned = (
ticks(no_orb.Phase_diagram.Band_gap.relabel(label=r"$\bm{B} \parallel (10, 1, 0)^T, \; \bm{A} = 0$")[0.5, nearest(kys, 0.05)], x=False)
+ ticks(no_orb.Phase_diagram.Band_gap.relabel(label=r"$\bm{B} \parallel (0, 1, 10)^T, \; \bm{A} = 0$")[nearest(kys, 0.05), 0.5], x=False, y=False)
+ ticks(orb.Phase_diagram.Band_gap.relabel(label=r"$\bm{B} \parallel (10, 1, 0)^T, \; \bm{A} \ne 0$")[0.5, nearest(kys, 0.05)])
+ ticks(orb.Phase_diagram.Band_gap.relabel(label=r"$\bm{B} \parallel (0, 1, 10)^T, \; \bm{A} \ne 0$", group='d', depth=2)[nearest(kys, 0.05), 0.5], y=False)).cols(2)
misaligned
In [ ]:
# renderer.save(misaligned, 'paper/figures/misaligned', fmt='pdf')
In [ ]:
# import os
# from scripts.hpc05 import HPC05Client
# os.environ['SSH_AUTH_SOCK'] = os.path.join(os.path.expanduser('~'), 'ssh-agent.socket')
# cluster = HPC05Client()
In [ ]:
# v = cluster[:]
# v.use_dill()
# lview = cluster.load_balanced_view()
# len(v)
In [ ]:
# %%px
# import sys
# import os
# sys.path.append(os.path.join(os.path.expanduser('~'), 'orbitalfield'))
# import kwant
# import numpy as np
# from fun import *
# lead = make_3d_wire()
# p = make_params(orbital=False, B_x=1)
In [ ]:
lead = make_3d_wire()
p = make_params(orbital=False, B_x=1)
mus = np.linspace(0, 35, 2000)
if os.path.exists('data/gaps_plot.npy'):
gaps = np.load('data/gaps_plot.npy')
else:
print('Start cluster with the cells above.')
gaps = lview.map_async(lambda mu: find_gap(lead, p, ((1, 0, 0), mu, True), tol=1e-4), mus).result()
np.save('data/gaps_plot', gaps)
if os.path.exists('data/spectrum_ev_plot.npy'):
Es = np.load('data/spectrum_ev_plot.npy')
else:
Es = np.array([kwant.physics.Bands(lead, args=[p])(k=0) for p.mu in mus])
np.save('data/spectrum_ev_plot', Es)
In [ ]:
%%output fig='svg'
%%opts VLine (lw=0.5) HLine (lw=0.5, color='g')
%%opts Layout [vspace=.35 aspect_weight=1 sublabel_position=(-0.3, 0.9) sublabel_format='({alpha})' sublabel_size=13]
%%opts Overlay [yticks=3 aspect=1.5 vspace=0.]
E_dim = hv.Dimension(('E_k0', r'$E(k=0)$'), unit='meV')
spectrum = hv.Path((mus, Es), kdims=[dimensions.mu, E_dim])
ind_E = 100
idx = np.argsort(np.min(np.abs(Es), axis=1))
VPoints = hv.Points([(mus[ind_E], E) for E in Es[ind_E]])
p.mu = 0
phase_bounds = np.sort(find_phase_bounds(lead, p, (1, 0, 0), num_bands=40).real)[::2]
HPoints = hv.Points([(x, 0) for x in phase_bounds if x > 0])(style={'color': 'g'})
ev_plot = (spectrum * hv.VLine(mus[ind_E]) * VPoints * HPoints * hv.HLine(0))[:35, -10:10]
bool_array = np.array(np.digitize(mus, phase_bounds)%2, dtype=bool)
gaps_plot = (spectrum
* hv.Area((mus, np.array(gaps) * bool_array))(style={'facecolor': '#FF6700'})
* hv.Area((mus, np.array(gaps) * ~bool_array))(style={'facecolor': '#a9a9a9'})
* hv.HLine(0) * HPoints)
gaps_plot = gaps_plot.map(lambda x: x.clone(extents=(0, 0, 35, 0.2)), [hv.Element])
ev_problem = (ev_plot[:, -8:8](plot={'xticks':[(0, ''), (8, ''), (16, ''), (24, ''), (32, '')],
'final_hooks': [tick_marks, hide_x]}) +
gaps_plot(plot={'xticks': 5, 'final_hooks': [tick_marks]})).cols(1)
ev_problem
In [ ]:
# renderer.save(ev_problem, 'paper/figures/ev_problem', fmt='pdf')