Spin dynamics calculations for Gd chains¶

# construct starting values for different helixes (single or double nodes)
x = np.arange(50)
x0 = 24.5

# one node
l = 16; phi = np.pi
m_helix1 = np.zeros((50,3))
m_helix1[:,0] = np.sin(np.pi*(x-x0)/l + phi)
m_helix1[:,2] = np.cos(np.pi*(x-x0)/l + phi)
plt.plot(x, m_helix1[:,2], color='C0', label='helix 1')

# two nodes
l = 11
m_helix2 = np.zeros((50,3))
m_helix2[:,0] = np.sin(np.pi*(x-x0)/l)
m_helix2[:,2] = np.cos(np.pi*(x-x0)/l)
plt.plot(x, m_helix2[:,2], color='C1', label='helix 2')

# three nodes
l = 7; phi = np.pi
m_helix3 = np.zeros((50,3))
m_helix3[:,0] = np.sin(np.pi*(x-x0)/l + phi)
m_helix3[:,2] = np.cos(np.pi*(x-x0)/l + phi)
plt.plot(x, m_helix3[:,2], color='C2', label='helix 3')

# four nodes
l = 5.8; phi = 0#np.pi
m_helix4 = np.zeros((50,3))
m_helix4[:,0] = np.sin(np.pi*(x-x0)/l + phi)
m_helix4[:,2] = np.cos(np.pi*(x-x0)/l + phi)
plt.plot(x, m_helix4[:,2], color='C3', label='helix 4')

# five nodes
l = 4.5; phi = np.pi
m_helix5 = np.zeros((50,3))
m_helix5[:,0] = np.sin(np.pi*(x-x0)/l + phi)
m_helix5[:,2] = np.cos(np.pi*(x-x0)/l + phi)
plt.plot(x, m_helix5[:,2], color='C4', label='helix 5')

plt.legend()

plt.show()

# this data is the contribution to the J and D with energy
J_D_data_short_EF = orm.load_node('02ea3f7d-fe61-46f3-ba78-287c93ef09f6').get_dict()
j = np.array(J_D_data_short_EF['J'])
d = np.array(J_D_data_short_EF['D'])
# since the above data does not integrate from -infty to EF+E, we need to shift the values to match the values for the original Fermi level
j -= (j[95]+j[96])/2 # substract value at E=0
d -= (d[95]+d[96])/2
j += 30.42 # meV (shift to correct value at EF)
d += 0.036 # meV

# energy range where J and D is computed
e = np.linspace(-2,2,192)


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

plt.subplot(1,2,1)
plt.plot(e[abs(e)<=1], j[abs(e)<=1], color='C0')
plt.axhline(0, color='C0')
plt.ylabel(r'$J(E)$ (meV)', color='C0')
plt.xlabel(r'$E-E_F$ (eV)')
plt.twinx()
plt.plot(e[abs(e)<=1], d[abs(e)<=1], color='C1')
plt.ylabel(r'$D(E)$ (meV)', color='C1')
plt.xlim(-1,1)

plt.subplot(1,2,2)
plt.plot(e[abs(e)<=0.65], d[abs(e)<=0.65] / j[abs(e)<=0.65], color='C0')
plt.axhline(0, color='grey', lw=1)
plt.ylabel(r'$D(E)\,/\, J(E)$', color='C0')
plt.xlabel(r'$E-E_F$ (eV)')

plt.twinx()
plt.plot(e[abs(e)<=0.65], np.arctan(d[abs(e)<=0.65] / j[abs(e)<=0.65])/np.pi*180, '--', color='C1')
plt.ylabel(r'$\theta$ (degrees)', color='C1')

plt.xlim(-0.65,0.65)

print('E, J, D:')
for E in [0, -0.2, -0.4, -0.6]:    
    plt.axvline(E, ls=':', color='grey', lw=1)
    ee = e[abs(e-E)==abs(e-E).min()][0]
    jj = j[abs(e-E)==abs(e-E).min()][0]
    dd = d[abs(e-E)==abs(e-E).min()][0]
    print(f'{ee:.2f} {jj:.3f} {dd:.3f} {dd/jj:.3f} {np.arctan(dd/jj)/np.pi*180:.3f}')


plt.tight_layout()
plt.show()

energies_all = []
plots = True

for E in [-0.6]:
    groupname0 =  'Gd_Nb110_chain_short' + f'/E={E:.1f}'
    print(groupname0)    
    energies_helix1 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix0')
        if calc.is_finished_ok:
            energies_helix1.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix1 = np.array(energies_helix1)
    
    
    energies_helix2 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix1')
        if calc.is_finished_ok:
            energies_helix2.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix2 = np.array(energies_helix2)
    
    
    
    energies_helix3 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix2')
        if calc.is_finished_ok:
            energies_helix3.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix3 = np.array(energies_helix3)


    energies_mz = []
    if plots:
        plt.figure()
    # for Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_mz')
        if calc.is_finished_ok:
            energies_mz.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_mz = np.array(energies_mz)
    
    energies_all.append([energies_helix1, energies_helix2, energies_helix3, energies_mz])

for (energies_helix1, energies_helix2, energies_helix3, energies_mz) in energies_all:

    E0 = energies_helix1[energies_helix1[:,0]==0.14,1]

    plt.figure()
    plt.title(f'short direction, $E=-0.6$ eV')
    plt.plot(energies_mz[energies_mz[:,0]>0.15,0], energies_mz[energies_mz[:,0]>0.15,1]-E0, 's-', label=r'$m_z$', ms=5)
    plt.plot(energies_helix1[:,0], energies_helix1[:,1]-E0, 'o-', label='helix 1 node', ms=4.5)
    plt.plot(energies_helix2[:,0], energies_helix2[:,1]-E0, 'o-', label='helix 2 nodes', ms=4.5)
    # plt.plot(energies_helix3[:,0], energies_helix3[:,1]-E0, 'o-', label='helix 3 nodes', ms=4.5)
    plt.legend(title='initial state')
    plt.ylabel(r'$E_\mathrm{tot}$ (eV / spin)')
    plt.xlabel(r'$B_z$ (T)')

    # plt.axvline(0.18, color='k', ls='--', lw=1)
    plt.axvline(0.140, color='k', ls='--', lw=1)
    plt.axvline(0.285, color='k', ls=':', lw=1)
    plt.axvline(0.355, color='k', ls='--', lw=1)
    
    plt.axvspan(0., 0.14, color='C2', alpha=0.3)
    plt.axvspan(0.14, 0.355, color='C1', alpha=0.3)
    plt.axvspan(0.355, 0.5, color='C0', alpha=0.3)

    plt.xlim(0, 0.5)

    plt.show()

init_spinview(vfr_frame_id='3', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.02', 'Gd_Nb110_chain_short/E=-0.6_helix1')
show_spins(calc, vfr_frame_id='3', scale_lattice=0.2)

init_spinview(vfr_frame_id='3_1', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.20', 'Gd_Nb110_chain_short/E=-0.6_helix0')
show_spins(calc, vfr_frame_id='3_1', scale_lattice=0.2)

init_spinview(vfr_frame_id='3_2', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.50', 'Gd_Nb110_chain_short/E=-0.6_mz')
show_spins(calc, vfr_frame_id='3_2', scale_lattice=0.2)

J_D_data_long_EF = orm.load_node('6c71914b-c28f-491e-9c6b-1c00716222ef').get_dict()
j = np.array(J_D_data_long_EF['J'])
d = np.array(J_D_data_long_EF['D'])

j -= (j[95]+j[96])/2
d -= (d[95]+d[96])/2
j += 1.18
d += 0.41

e = np.linspace(-2,2,192)

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

# plt.subplot(1,2,1)
# plt.plot(e[abs(e)<=1], j[abs(e)<=1], color='C0')
# plt.axhline(0, color='C0')
# plt.ylabel(r'$J(E)$ (meV)', color='C0')
# plt.xlabel(r'$E-E_F$ (eV)')
# plt.twinx()
# plt.plot(e[abs(e)<=1], d[abs(e)<=1], color='C1')
# plt.ylabel(r'$D(E)$ (meV)', color='C1')
# plt.xlim(-1,1)

# plt.subplot(1,2,2)
plt.plot(e[abs(e)<=0.65], d[abs(e)<=0.65] / j[abs(e)<=0.65], color='C0')
plt.axhline(0, color='grey', lw=1)
plt.ylabel(r'$D(E)\,/\, J(E)$', color='C0')
plt.xlabel(r'$E-E_F$ (eV)')

plt.twinx()
plt.plot(e[abs(e)<=0.65], np.arctan(d[abs(e)<=0.65] / j[abs(e)<=0.65])/np.pi*180, '--', color='C1')
plt.ylabel(r'$\theta$ (degrees)', color='C1')

plt.xlim(-0.4,0.4)

print('E, J, D, D/J, theta:')
for E in [-0.05, 0.0, 0.05]:
    plt.axvline(E, ls=':', color='grey', lw=1)
    ee = e[abs(e-E)==abs(e-E).min()][0]
    jj = j[abs(e-E)==abs(e-E).min()][0]
    dd = d[abs(e-E)==abs(e-E).min()][0]
    print(f'{ee:.2f} {jj:.3f} {dd:.3f} {dd/jj:.3f} {np.arctan(dd/jj)/np.pi*180:.3f}')

plt.savefig('D_J_ratio_long.svg')
plt.show()

energies_all_long = []
plots = True

for E in [-0.05, 0.0, 0.05]:
    groupname0 =  'Gd_Nb110_chain_long' + f'/E={E:.1f}'
    print(groupname0)
    
    energies_helix1 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix0')
        if calc.is_finished_ok:
            energies_helix1.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix1 = np.array(energies_helix1)
    
    
    energies_helix2 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix1')
        if calc.is_finished_ok:
            energies_helix2.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix2 = np.array(energies_helix2)
    
    
    
    energies_helix3 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix2')
        if calc.is_finished_ok:
            energies_helix3.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix3 = np.array(energies_helix3)
    
    
    
    energies_helix4 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix3')
        if calc.is_finished_ok:
            energies_helix4.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix4 = np.array(energies_helix4)
    
    
    
    energies_helix5 = []

    if plots:
        plt.figure()
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_helix4')
        if calc.is_finished_ok:
            energies_helix5.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_helix5 = np.array(energies_helix5)


    energies_mz = []
    if plots:
        plt.figure()
    # for Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
    for Bz in tqdm(np.linspace(0, 0.5, 51)):
        calc = get_calc_from_group(f'Bz={Bz:.2f}', groupname0 + '_mz')
        if calc.is_finished_ok:
            energies_mz.append([Bz, calc.outputs.energies.get_array('energies')[-1,1]])
            if Bz in [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 1, 2]:
                m = calc.outputs.magnetization.get_array('final')
                if plots:
                    plt.plot(m[:,2], '.--', label=f'{Bz:.2f}')

    if plots:
        plt.legend(title=r'$B_z$ (T)')
        plt.xlabel('site index')
        plt.ylabel('$m_z$')
        plt.show()

    energies_mz = np.array(energies_mz)
    
    energies_all_long.append([energies_helix1, energies_helix2, energies_helix3, energies_helix4, energies_helix5, energies_mz])

for ie, (energies_helix1, energies_helix2, energies_helix3, energies_helix4, energies_helix5, energies_mz) in enumerate([energies_all_long[0]]):
    E = [-0.05, 0.0, 0.05][ie]

    E0 = energies_helix1[energies_helix1[:,0]==0.14,1]

    plt.figure()
    plt.title(f'$E={E:.2f}$ eV')
    plt.plot(energies_helix2[energies_helix2[:,0]>0.1,0], energies_helix2[energies_helix2[:,0]>0.1,1]-E0, 's-', label='m_z', ms=5)
    plt.plot(energies_helix1[energies_helix1[:,0]>0,0], energies_helix1[energies_helix1[:,0]>0,1]-E0, 'o-', label='helix 1', ms=4.5)
    plt.plot(energies_helix2[energies_helix2[:,0]<=0.1,0], energies_helix2[energies_helix2[:,0]<=0.1,1]-E0, '.-', label='helix 2', ms=4.5)
    plt.legend(title='initial state')
    plt.ylabel(r'$E_\mathrm{tot}$ (eV / spin)')
    plt.xlabel(r'$B_z$ (T)')
    
    plt.axvspan(0., 0.105, color='C1', alpha=0.5)
    plt.axvspan(0.105, 0.5, color='C0', alpha=0.5)

    plt.xlim(0, 0.5)
    
    plt.savefig('E_tot_long_E=-50meV.svg')

    plt.show()

init_spinview(vfr_frame_id='4_1', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.02', 'Gd_Nb110_chain_long/E=-0.1_helix0')
show_spins(calc, vfr_frame_id='4_1', scale_lattice=0.2)

init_spinview(vfr_frame_id='4_2', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.50', 'Gd_Nb110_chain_long/E=-0.1_mz')
show_spins(calc, vfr_frame_id='4_2', scale_lattice=0.2)

for ie, (energies_helix1, energies_helix2, energies_helix3, energies_helix4, energies_helix5, energies_mz) in enumerate([energies_all_long[1]]):

    E = 0.0

    E0 = energies_helix1[energies_helix2[:,0]==0,1]

    plt.figure()
    plt.plot(energies_mz[energies_mz[:,0]>0.15,0], energies_mz[energies_mz[:,0]>0.15,1]-E0, 's-', label=r'$m \, || \, \hat{z}$', ms=5)
    plt.plot(energies_helix1[energies_helix1[:,0]>0.08,0], energies_helix1[energies_helix1[:,0]>0.08,1]-E0, 'o-', label='helix 1', ms=4.5)
    plt.plot(energies_helix2[energies_helix2[:,0]<0.43,0], energies_helix2[energies_helix2[:,0]<0.43,1]-E0, 'd-', label='helix 2', ms=4.5)
    plt.legend(title='initial state')
    plt.ylabel(r'$E_\mathrm{tot}$ (eV / spin)', fontsize=12)
    plt.xlabel(r'$B_z$ (T)', fontsize=12)
    
    plt.axvspan(0., 0.165, color='C2', alpha=0.5)
    plt.axvspan(0.165, 0.195, color='C1', alpha=0.5)
    plt.axvspan(0.195, 0.5, color='C0', alpha=0.5)

    plt.xlim(0, 0.5)
    
    
    plt.savefig(f'E_tot_phases_long_dimer_E={E:.2f}eV.svg')

    plt.show()

init_spinview(vfr_frame_id='5_1', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.02', 'Gd_Nb110_chain_long/E=0.0_helix1')
show_spins(calc, vfr_frame_id='5_1', scale_lattice=0.2)

init_spinview(vfr_frame_id='5_2', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.18', 'Gd_Nb110_chain_long/E=0.0_helix0')
show_spins(calc, vfr_frame_id='5_2', scale_lattice=0.2)

init_spinview(vfr_frame_id='5_3', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.30', 'Gd_Nb110_chain_long/E=0.0_mz')
show_spins(calc, vfr_frame_id='5_3', scale_lattice=0.2)

for ie, (energies_helix1, energies_helix2, energies_helix3, energies_helix4, energies_helix5, energies_mz) in enumerate([energies_all_long[2]]):
    E = 0.05

    E0 = energies_helix5[energies_helix5[:,0]==0,1]

    plt.figure()
    plt.title(f'$E={E:.2f}$ eV')
    plt.plot(energies_mz[energies_mz[:,0]>0.18,0], energies_mz[energies_mz[:,0]>0.18,1]-E0, 's-', label=r'$m_z$', ms=5)
    plt.plot(energies_helix3[energies_helix3[:,0]>0.21,0], energies_helix3[energies_helix3[:,0]>0.21,1]-E0, '.--', label='helix 3', ms=4.5)
    e5 = energies_helix5[energies_helix5[:,0]<0.35,:]
    e5 = e5[abs(e5[:,0]-0.15)>0.01]
    e5 = e5[abs(e5[:,0]-0.23)>0.01]
    plt.plot(e5[:,0], e5[:,1]-E0, '.--', label='helix 5', ms=4.5)
    plt.legend(title='initial state')
    plt.ylabel(r'$E_\mathrm{tot}$ (eV / spin)')
    plt.xlabel(r'$B_z$ (T)')
    
    plt.axvspan(0., 0.255, color='C2', alpha=0.5)
    plt.axvspan(0.255, 0.405, color='C1', alpha=0.5)
    plt.axvspan(0.405, 0.5, color='C0', alpha=0.5)

    plt.xlim(0, 0.5)
    
    plt.savefig('E_tot_long_E+50meV.svg')

    plt.show()

init_spinview(vfr_frame_id='6_1', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.02', 'Gd_Nb110_chain_long/E=0.1_helix4')
show_spins(calc, vfr_frame_id='6_1', scale_lattice=0.2)

init_spinview(vfr_frame_id='6_2', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.30', 'Gd_Nb110_chain_long/E=0.1_helix2')
show_spins(calc, vfr_frame_id='6_2', scale_lattice=0.2)

init_spinview(vfr_frame_id='6_3', height_px=200)

sleep(10) # wait 10 seconds for spin view frame to initialize
calc = get_calc_from_group('Bz=0.30', 'Gd_Nb110_chain_long/E=0.1_mz')
show_spins(calc, vfr_frame_id='6_3', scale_lattice=0.2)

Bcrit = np.array([
    [-0.05, 0.105],
    [0.0, 0.195],
    [0.05, 0.405]])

plt.figure()
plt.plot(1000*Bcrit[:,0], Bcrit[:,1], 'o--', lw=3, ms=10)
plt.ylabel(r'$B^c_z$ (T)', color='C0')
plt.xlabel(r'$E-E_F$ (meV)')
plt.xticks([-50, -25, 0, 25, 50])
plt.yticks([0.10, 0.15, 0.20, 0.25, 0.30, 0.35, 0.40])


plt.twinx()

J_D_data_long_EF = orm.load_node('6c71914b-c28f-491e-9c6b-1c00716222ef').get_dict()
j = np.array(J_D_data_long_EF['J'])
d = np.array(J_D_data_long_EF['D'])

j -= (j[95]+j[96])/2
d -= (d[95]+d[96])/2
j += 1.18
d += 0.41

e = np.linspace(-2,2,192)


plt.plot(1000*e[abs(e)<=0.08], np.arctan(d[abs(e)<=0.08] / j[abs(e)<=0.08])/np.pi*180, '-', color='C1', lw=3)
plt.xlim(-0.06*1000, 0.06*1000)
plt.ylim(10.8, 37.2)
plt.ylabel(r'$\theta = \arctan(D\,/\,J)$ (degrees)', color='C1')
plt.yticks([15, 20, 25, 30, 35])

plt.savefig('inset_Bcrit.svg')


plt.show()