Select Git revision
final_plots.py
final_plots.py 60.29 KiB
#!/usr/bin/env python3
# Copyright 2022 Valentin Bruch <valentin.bruch@rwth-aachen.de>
# License: MIT
"""
Kondo FRTRG, generate high-quality plots for publication
"""
import os
import itertools
import scipy.constants as sc
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mplcolors
from matplotlib.widgets import Slider
import argparse
import numpy as np
from scipy.interpolate import bisplrep, bisplev, splrep, splev, BSpline, griddata, interp1d
from scipy.special import jn
from imageio import imwrite
from scipy.integrate import quad, quad_vec
from scipy.optimize import newton
import settings
from data_management import DataManager, KondoImport
from kondo import gen_init_matrix
from numbers import Number
from gen_pulse_data import fourier_coef_gauss_symmetric
from plot_pulses import load_all_pulses_full, integrate_ft
# In this program all energies are given in units of the RTRG Kondo
# temperature Tkrg, which is an integration constant of the E-flow RG
# equations. The more conventional definition of the Kondo temperature is
# G(V=Tk)=G(V=0)/2=e²/h. The ratio Tk/Tkrg is:
#TK_VOLTAGE = 3.44249 # for rtol=1e-10, atol=1e-12, voltage_branches=10
TK_VOLTAGE = 3.4425351 # for rtol=1e-9, atol=1e-11, voltage_branches=10
#TK_VOLTAGE = 3.44334 # with correction
TK_VOLTAGE_O3P = 3.458524 # for rtol=1e-9, atol=1e-11, voltage_branches=10
#TK_VOLTAGE_O3P = 3.4593 # with correction
TK_VOLTAGE_ORDER2 = 10.1368086 # for rtol=1e-9, atol=1e-11, voltage_branches=10
#TK_VOLTAGE_ORDER2 = 10.13754 # with correction
def save_overview(
omega = 16.5372,
vdc_res = 501,
vac_res = 501,
vdc_max = 165.372,
vac_max = 165.372,
method = "mu",
d = 1e9,
xL = 0.5,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
lazy_inverse_factor = 0.25,
voltage_branches = 4,
s_g = 1e-4,
s_idc = 5e-6,
s_iac = 5e-6,
filename = "figdata/vdc_vac_omega16.5372_interp.npz",
**kwargs
):
"""
3d plot of ac and dc differential conductance and current as function of Vac and Vdc.
"""
dm = DataManager()
data = dm.list(omega=omega, vdc=None, vac=None, method=method, d=d, xL=xL, solver_tol_abs=solver_tol_abs, solver_tol_rel=solver_tol_rel, voltage_branches=voltage_branches, lazy_inverse_factor=lazy_inverse_factor, **kwargs)
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"]
filters = (data.vac < 1.2*vac_max) \
& (data.vdc < 1.2*vdc_max) \
& np.isfinite(data.dc_conductance) \
& (data.solver_flags & bad_flags == 0)
data = data.loc[filters]
print("Interpolate Gdc", flush=True)
gdc_tck = bisplrep(data.vac, data.vdc, data.dc_conductance, s=s_g, kx=3, ky=3)
print("Interpolate Idc", flush=True)
idc_tck = bisplrep(data.vac, data.vdc, data.dc_current, s=s_idc, kx=3, ky=3)
print("Interpolate Iac", flush=True)
iac_tck = bisplrep(data.vac, data.vdc, data.ac_current_abs, s=s_iac, kx=3, ky=3)
print("done.")
vac_arr = np.linspace(0, min(vac_max, data.vac.max()), vac_res)
vdc_arr = np.linspace(0, min(vdc_max, data.vdc.max()), vdc_res)
gdc_g_interp = bisplev(vac_arr, vdc_arr, gdc_tck).T
idc_interp = bisplev(vac_arr, vdc_arr, idc_tck).T
iac_interp = bisplev(vac_arr, vdc_arr, iac_tck).T
gdc_i_interp = bisplev(vac_arr, vdc_arr, idc_tck, dy=1).T
gac_interp = 2*bisplev(vac_arr, vdc_arr, iac_tck, dx=1).T
np.savez(filename, vac=vac_arr/TK_VOLTAGE, vdc=vdc_arr/TK_VOLTAGE, gdc_g=np.pi*gdc_g_interp, gdc_i=np.pi*gdc_i_interp, gac=np.pi*gac_interp, idc=idc_interp/TK_VOLTAGE, iac=iac_interp/TK_VOLTAGE)
#vac_data, vdc_data = np.meshgrid(vac_arr, vdc_arr)
#with open(filename, 'w') as file:
# file.write("vdc,vac,gdc_g,gdc_i,gac,idc,iac\n")
# np.savetxt(file, np.array([vdc_data/TK_VOLTAGE, vac_data/TK_VOLTAGE, np.pi*gdc_g_interp, np.pi*gdc_i_interp, np.pi*gac_interp, idc_interp/TK_VOLTAGE, iac_interp/TK_VOLTAGE]).reshape((7,-1)).T, fmt="%.9e", delimiter=",")
# #for i in range(vdc_res):
# # file.write('\n')
# # np.savetxt(file, np.array([vdc_data[i]/omega, vac_data[i]/omega, np.pi*gdc_g_interp[i], np.pi*gdc_i_interp[i], np.pi*gac_interp[i], idc_interp[i], iac_interp[i]]).T)
def interp_vac0(dm, xL=0.5, include_Ga=True, truncation_order=3, integral_method=-15, **trashparams):
"""Return differential conductance G as function of Vdc at Vac=0"""
data = dm.list(vac=0, omega=0, include_Ga=include_Ga, d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, nmax=0, voltage_branches=4, xL=xL, truncation_order=truncation_order, integral_method=integral_method)
data = data.sort_values("vdc")
vdc = np.concatenate((-data.vdc[:0:-1], data.vdc))
gdc = np.concatenate((data.dc_conductance[:0:-1], data.dc_conductance))
return interp1d(vdc, gdc, kind="cubic")
def photon_assisted_tunneling(dm=None, omega=0, vdc=0, vac=0, nmax=50, interp=None, **parameters):
"""
Return differential conductance G calculated from photon
assisted tunneling at given parameters.
"""
result = np.zeros(np.broadcast_shapes(np.shape(omega), np.shape(vdc), np.shape(vac)), dtype=np.float64)
if interp is None:
interp = interp_vac0(dm, **parameters)
for n in range(-nmax, nmax+1):
try:
result += jn(n, vac/omega)**2 * interp(vdc+n*omega)
except ValueError:
pass
return result
def gen_photon_assisted_tunneling(dm, fourier_coef, omega, vdc, nmax=50, **parameters):
"""
Return differential conductance G calculated from generalized
photon assisted tunneling at given parameters.
check:
>>> gen_photon_assisted_tunneling(dm, (vac/2,0), omega, vdc, ...) \
>>> == photon_assisted_tunneling(dm, omega, vdc, vac, ...)
"""
result = np.zeros(np.broadcast_shapes(np.shape(omega), np.shape(vdc)), dtype=np.complex128)
interp = interp_vac0(dm, **parameters)
f0 = gen_init_matrix(nmax, *(f/omega for f in fourier_coef))
for n in range(-nmax, nmax+1):
try:
result += (abs(f0[nmax+n,nmax])**2 + abs(f0[nmax,nmax+n])**2)/2 * interp(vdc+n*omega)
except ValueError:
pass
assert (np.abs(result.imag) < 1e-9).all()
return result.real
def filter_grid_data(
data : pd.DataFrame,
omega = 16.5372,
vac_min = 0,
vac_max = 165.372,
vac_num = 101,
vdc_min = 0,
vdc_max = 165.372,
vdc_num = 101,
v_tol = 1e-6,
):
vac_step = (vac_max - vac_min) / (vac_num - 1)
vdc_step = (vdc_max - vdc_min) / (vdc_num - 1)
data = data.loc[(data.vac >= vac_min - v_tol) & (data.vdc >= vdc_min - v_tol) & (data.vac <= vac_max + v_tol) & (data.vdc <= vdc_max + v_tol) & np.isfinite(data.dc_conductance)]
grid_data = data.loc[(np.abs(((data.vac - vac_min + v_tol) % vac_step) - v_tol) < v_tol) & (np.abs(((data.vdc - vdc_min + v_tol) % vdc_step) - v_tol) < v_tol)]
data = grid_data.sort_values(["vac", "vdc"])
vac_arr = np.linspace(vac_min, vac_max, vac_num)
vdc_arr = np.linspace(vdc_min, vdc_max, vdc_num)
gdc_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
idc_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
iac_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
phase_arr = np.empty((vac_num, vdc_num), dtype=np.float64)
gdc_arr.fill(np.nan)
idc_arr.fill(np.nan)
iac_arr.fill(np.nan)
phase_arr.fill(np.nan)
lower_index = 0
for i, vac in enumerate(vac_arr):
upper_index = data.vac.searchsorted(vac + v_tol)
indices = vdc_arr.searchsorted(data.vdc[lower_index:upper_index] - v_tol)
indices = indices[np.abs(vdc_arr[indices] - data.vdc[lower_index:upper_index]) < v_tol]
gdc_arr[i, indices] = data.dc_conductance[lower_index:upper_index]
idc_arr[i, indices] = data.dc_current[lower_index:upper_index]
iac_arr[i, indices] = data.ac_current_abs[lower_index:upper_index]
phase_arr[i, indices] = data.ac_current_phase[lower_index:upper_index]
lower_index = upper_index
return *np.meshgrid(vdc_arr, vac_arr), gdc_arr, idc_arr, iac_arr, phase_arr
def export_matrix_pgfplots(filename, *arrays, header="", fmt="%.6g"):
array = np.array(arrays).T
with open(filename, "w") as file:
if header:
file.write(header)
for sector in array:
file.write("\n")
np.savetxt(file, sector, fmt=fmt)
def export_omega5_pgfplots(filename="figdata/omega5_interp.dat", dc_steps=3, ac_steps=2):
data = np.load("figdata/omega5_interp.npz")
omega = 16.5372
vdc = data["vdc"][::ac_steps,::dc_steps]/omega
vac = data["vac"][::ac_steps,::dc_steps]/omega
gdc = data["gdc_mu_o3a"][::ac_steps,::dc_steps]*np.pi
gac = data["gac_mu_o3a_ispline"][::ac_steps,::dc_steps]*np.pi
idc = data["idc_mu_o3a"][::ac_steps,::dc_steps]
iac = data["iac_mu_o3a"][::ac_steps,::dc_steps]
phase = data["phase_mu_o3a"][::ac_steps,::dc_steps]
export_matrix_pgfplots(filename, vdc, vac, gdc, gac, idc, iac, phase, header="vdc vac gdc gac idc iac phase")
def export_omega5():
"deprecated!"
omega = 16.5372
dm = DataManager()
# Full overview
data = dm.list(
omega = omega,
vdc = None,
vac = None,
method = "mu",
d = 1e9,
xL = 0.5,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
voltage_branches = 4,
truncation_order = 3,
include_Ga = True,
solve_integral_exactly = False,
)
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["deleted"]
data = data.loc[data.solver_flags & bad_flags == 0]
vdc, vac, gdc, idc, iac, phase = filter_grid_data(
data,
omega = omega,
vac_min = 0,
vac_max = 165.372,
vac_num = 51,
vdc_min = 0,
vdc_max = 165.372,
vdc_num = 51,
)
with open("figdata/vdc_vac_omega16.5372.dat", "w") as file:
file.write("vac vdc gdc idc iac")
for i in range(51):
file.write("\n")
np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
## Lines at constant Vac
#for i in range(11):
# with open(f"figdata/vdc_vac{i}_omega16.5372.dat", "w") as file:
# file.write("vac vdc gdc idc iac\n")
# np.savetxt(file, np.array([vac[5*i]/omega, vdc[5*i]/omega, np.pi*gdc[5*i], idc[5*i], iac[5*i]]).T)
## Lines at constant Vdc
#for i in range(11):
# with open(f"figdata/vac_vdc{i}_omega16.5372.dat", "w") as file:
# file.write("vac vdc gdc idc iac\n")
# np.savetxt(file, np.array([vac[:,5*i]/omega, vdc[:,5*i]/omega, np.pi*gdc[:,5*i], idc[:,5*i], iac[:,5*i]]).T)
# Zoom to primary Kondo peak
vdc, vac, gdc, idc, iac, phase = filter_grid_data(
data,
omega = omega,
vac_min = 0,
vac_max = 16.5372,
vac_num = 21,
vdc_min = 0,
vdc_max = 16.5372,
vdc_num = 21,
)
with open("figdata/vdc_vac_omega16.5372_zoom.dat", "w") as file:
file.write("vac vdc gdc idc iac")
for i in range(21):
file.write("\n")
np.savetxt(file, np.array([vac[i]/omega, vdc[i]/omega, np.pi*gdc[i], idc[i], iac[i]]).T)
def export_interp(
filename,
fixed_parameter,
fixed_value,
x_parameter,
y_parameter,
x_arr,
y_arr,
x_sort_parameter = None,
y_sort_parameter = None,
x_func = None,
y_func = None,
kxorder = 2,
kyorder = 2,
spline_s = 2e-4,
special_selection = None,
o2_scale_x = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
o2_scale_y = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
o2_scale_fixed = TK_VOLTAGE_ORDER2/TK_VOLTAGE,
**parameters
):
"""
Export interpolated data. Parameters:
filename: string, file path
fixed_parameter: name of fixed parameter (e.g. "omega")
fixed_value: value of fixed parameter (e.g. 16.5372)
x_parameter: name of x axis parameter (e.g. "vdc")
y_parameter: name of y axis parameter (e.g. "vac")
x_sort_parameter: name of x axis parameter for sorting
y_sort_parameter: name of y axis parameter for sorting
x_func: function to obtain x value from data
y_func: function to obtain y value from data
kxorder: order of spline interpolation in x direction
kyorder: order of spline interpolation in y direction
special_selection: function for filtering input data. Arguments:
data, order, method, padding, **kwargs
Data are stored with units e=hbar=kB=Tkrg=1.
All results are exported as 2d arrays of the same shape.
The parameters vdc and vac are included as 2d arrays.
Results have the naming convention
{observable}_{method}_{order}{padding}{suffix}
observable:
gdc DC differential conductance, may have suffix "_ispline"
gac AC differential conductance, must have suffix "_ispline"
idc average current, may have suffix "_spline"
iac oscillating current, may have suffix "_spline"
phase AC phase
gamma Γ(E=0) (not an observable)
order:
o2 2nd order RG equations
o3 3rd order RG equations without Ga and with approximated integral
o3a 3rd order RG equations with approximated integral
o3p 3rd order RG equations with full integral
padding:
"" no Floquet matrix extrapolation
"_p" Floquet matrix extrapolation to mitigate truncation effects
method:
mu no unitary transformation, oscillating chemical potentials
J include oscillating voltage in coupling by unitary transformation
suffix
"" interpolation is done using griddata
"_spline" interpolation is done using spline
"_ispline" data is derived from spline interpolation for the
current (dc or ac)
"""
assert isinstance(filename, str)
if x_func is None:
x_func = lambda data: data[x_parameter]
if y_func is None:
y_func = lambda data: data[y_parameter]
x_mesh, y_mesh = np.meshgrid(x_arr, y_arr)
mesh = (x_mesh, y_mesh)
dm = DataManager()
data = dm.list(**{fixed_parameter:fixed_value}, **parameters)
global_bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"]
data = data.loc[data.solver_flags & global_bad_flags == 0]
data = data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
required_flags = dict(
o2 = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
o3 = 0,
o3a = DataManager.SOLVER_FLAGS["include_Ga"],
o3p = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
)
bad_flags = dict(
o2 = DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
o3 = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["include_Ga"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
o3a = DataManager.SOLVER_FLAGS["second_order_rg_equations"] | DataManager.SOLVER_FLAGS["solve_integral_exactly"],
o3p = DataManager.SOLVER_FLAGS["second_order_rg_equations"],
)
o2_data = dm.list(**{fixed_parameter:fixed_value*o2_scale_fixed}, **parameters, min_version=(14,15,-1,-1))
o2_data = o2_data.loc[o2_data.solver_flags & (global_bad_flags | bad_flags["o2"] | required_flags["o2"]) == required_flags["o2"]]
o2_data = o2_data.sort_values([x_sort_parameter or x_parameter, y_sort_parameter or y_parameter])
def selection(order, method, padding=False, **kwargs):
thedata = o2_data if order == "o2" else data
selected = (thedata.solver_flags & required_flags[order] == required_flags[order]) \
& (thedata.solver_flags & bad_flags[order] == 0) \
& (thedata.method == method)
if padding:
selected &= (thedata["padding"] > 0) | (thedata["nmax"] == 0)
else:
selected &= thedata["padding"] == 0
for key, value in kwargs.items():
selected &= (thedata[key] == value)
if special_selection is not None:
selected &= special_selection(thedata, order, method, padding, **kwargs)
return thedata.loc[selected]
def gen_data(order, method, padding=False, s=spline_s, extend_vdc=10, extend_vac=10, **kwargs):
suffix = "_p" if padding else ""
reduced_data = selection(order, method, padding, **kwargs)
settings.logger.info(f"Starting with {order}{suffix} {method}, found {reduced_data.shape[0]} data points")
if extend_vdc or extend_vac:
extended_data = [reduced_data]
if extend_vdc:
dc_mirrored = reduced_data[reduced_data.vdc<extend_vdc].copy()
dc_mirrored.vdc *= -1
dc_mirrored.dc_current *= -1
extended_data.append(dc_mirrored)
del dc_mirrored
if extend_vac:
ac_mirrored = reduced_data[reduced_data.vac<extend_vac].copy()
ac_mirrored.vac *= -1
ac_mirrored.ac_current_phase += np.pi
extended_data.append(ac_mirrored)
if extend_vdc:
acdc_mirrored = ac_mirrored[ac_mirrored.vdc<extend_vdc].copy()
acdc_mirrored.vdc *= -1
acdc_mirrored.dc_current *= -1
extended_data.append(acdc_mirrored)
del acdc_mirrored
del ac_mirrored
reduced_data = pd.concat(extended_data)
results = {}
if reduced_data.shape[0] < 100:
return results
xy_data = (x_func(reduced_data), y_func(reduced_data))
if order == "o2":
xy_data = (xy_data[0]/o2_scale_x, xy_data[1]/o2_scale_y)
try:
results[f"gdc_{method}_{order}{suffix}"] = griddata(
xy_data,
reduced_data.dc_conductance,
mesh,
method="cubic")
#gdc_tck = bisplrep(reduced_data[y], reduced_data[x], reduced_data.dc_conductance, s=s, kx=kxorder, ky=kyorder)
#results[f"gdc_{method}_{order}{suffix}"] = bisplev(y_arr, x_arr, gdc_tck)
except:
settings.logger.exception(f"in griddata interpolation of Gdc for {order}{suffix} {method}")
try:
results[f"idc_{method}_{order}{suffix}"] = griddata(
xy_data,
reduced_data.dc_current,
mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of Idc for {order}{suffix} {method}")
try:
idc_tck = bisplrep(*xy_data[::-1], reduced_data.dc_current, s=s, kx=kxorder, ky=kyorder)
results[f"idc_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, idc_tck)
results[f"gdc_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, idc_tck, dy=1)
except:
settings.logger.exception(f"in spline interpolation of Idc for {order}{suffix} {method}")
try:
results[f"iac_{method}_{order}{suffix}"] = griddata(
xy_data,
2*reduced_data.ac_current_abs,
mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of Iac for {order}{suffix} {method}")
try:
iac_tck = bisplrep(*xy_data[::-1], 2*reduced_data.ac_current_abs, s=s, kx=kxorder, ky=kyorder)
results[f"iac_{method}_{order}{suffix}_spline"] = bisplev(y_arr, x_arr, iac_tck)
results[f"gac_{method}_{order}{suffix}_ispline"] = bisplev(y_arr, x_arr, iac_tck, dx=1)
except:
settings.logger.exception(f"in spline interpolation of Iac for {order}{suffix} {method}")
try:
sel = reduced_data.vac > 0
results[f"phase_{method}_{order}{suffix}"] = griddata(
(xy_data[0][sel], xy_data[1][sel]),
reduced_data.ac_current_phase[sel],
mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of phase for {order}{suffix} {method}")
try:
results[f"gamma_{method}_{order}{suffix}"] = griddata(
xy_data,
reduced_data.gamma,
mesh,
method="cubic")
except:
settings.logger.exception(f"in griddata interpolation of gamma for {order}{suffix} {method}")
return results
np.savez(
filename,
**{fixed_parameter:fixed_value, x_parameter:x_mesh, y_parameter:y_mesh},
**gen_data("o2", "mu"),
**gen_data("o3", "mu"),
**gen_data("o3a", "mu"),
**gen_data("o3p", "mu", integral_method=-1),
**gen_data("o2", "J", False),
**gen_data("o3", "J", False),
**gen_data("o3a", "J", False),
**gen_data("o3p", "J", False, integral_method=-1),
**gen_data("o2", "J", True),
**gen_data("o3", "J", True),
**gen_data("o3a", "J", True),
**gen_data("o3p", "J", True, integral_method=-1),
)
def export_omega5_interp(
filename = "figdata/omega5_interp.npz",
omega = 16.5372,
vdc_min = 0,
vdc_max = 165.372,
vac_min = 0,
vac_max = 165.372,
dc_res = 301,
ac_res = 201,
vdc_max_J = 82.686,
vac_max_J = 82.686,
voltage_branches = 4,
):
"""
Export interpolated data for fixed frequency omega.
"""
def special_selection(data, order, method, padding, **kwargs):
if method == "J":
return (data["vdc"] < vdc_max_J + 1e-3) & (data["vac"] < vac_max_J + 1e-3)
return True
yarr = np.linspace(vac_min, vac_max, ac_res)
yarr[-1] -= 1e-6
export_interp(filename,
fixed_parameter = "omega",
fixed_value = omega,
x_parameter = "vdc",
y_parameter = "vac",
x_arr = np.linspace(vdc_min, vdc_max, dc_res),
y_arr = yarr,
kyorder = 3,
special_selection = special_selection,
spline_s = 2e-5,
d = 1e9,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
voltage_branches = voltage_branches,
xL = 0.5)
def export_vdc0_interp(
filename = "figdata/vdc0_interp.npz",
omega_min = 0.1,
omega_max = 40,
vac_omega_min = 0.01,
vac_omega_max = 10,
omega_res = 201,
vac_res = 201,
korder = 2,
):
"""
Export interpolated data for fixed frequency omega.
"""
parameters = dict(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, voltage_branches=0, xL=0.5)
export_interp(filename,
fixed_parameter = "vdc",
fixed_value = 0,
x_parameter = "omega",
y_parameter = "vac_omega",
y_sort_parameter = "vac",
x_arr = np.linspace(omega_min, omega_max, omega_res),
y_arr = np.linspace(vac_omega_min, vac_omega_max, vac_res),
y_func = (lambda data: data["vac"]/data["omega"]),
korder = korder,
special_selection = (lambda data, *args, **kwargs: data["omega"] > 0),
**parameters)
def prepare_plotly_csv():
dm = DataManager()
# General overview
reduction_dict = dict(omega="omega", vdc="vdc", vac="vac", dc_conductance="g", dc_current="idc", ac_current_abs="iac", ac_current_phase="ac_phase")
data = dm.list(d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, truncation_order=3)
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["include_Ga"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["deleted"]
data = data.loc[data.solver_flags & bad_flags == 0]
data = data.rename(columns=reduction_dict)
data.g *= np.pi
data.omega /= TK_VOLTAGE
data.vdc /= TK_VOLTAGE
data.vac /= TK_VOLTAGE
data.idc /= TK_VOLTAGE
data.iac /= TK_VOLTAGE
data = data.sort_values(["vac","vdc","omega"])
#data.to_csv("html/full.csv", columns=reduction_dict.values())
# omega = 5Tk
# TODO: this is outdated!
for omega, name in zip((16.5372, 9.2159791, 5.8206184, 7.1271), ("omega5", "compare_bruhat18a", "compare_bruhat18b", "compare_kogan04")):
reduced = data.loc[np.isclose(data.omega*TK_VOLTAGE, omega) & (data.method == "mu") & (data.voltage_branches == 4)]
reduced.to_csv(f"html/{name}.csv", columns=reduction_dict.values())
# vdc = 0
reduced = data.loc[(data.vdc == 0) & (data.method != "mu")]
reduced.to_csv("html/vdc0.csv", columns=reduction_dict.values())
def prepare_RGconvergence():
dc_indices = np.array([0, 9, 21, 30, 90, 165])
ac_indices = np.array([32, 52, 80, 124, 172])
real_cmap = lambda x: plt.cm.seismic((x+1)/2)[...,:3]
imwrite(f"figdata/colorbar_seismic.png", np.array(0xffff*plt.cm.seismic(np.linspace(0, 1, 0x1000))[...,:3], dtype=np.uint16).reshape((1,0x1000,3)), format="PNG-FI")
def export_img(name, array):
imwrite(f"figdata/{name}.png", np.array(0xffff*real_cmap(array[::-1]), dtype=np.uint16), format="PNG-FI")
data = np.load("figdata/omega5_interp.npz")
data_vb7 = np.load("figdata/omega5_interp_vb7.npz")
# O2
omega = 16.5372
omega_o2 = 48.695054
vdc = data["vdc"][0] / omega
vac = data["vac"][:,0] / omega
for method, ref in (("o2","o3p"), ("o3","o3p"), ("o3a","o3p"), ("o3a_ispline","o3a")):
gdiff = data[f"gdc_mu_{method}"] - data[f"gdc_mu_{ref}"]
normalize = np.nanmax(np.abs(gdiff))
export_img(f"convergence_{method}_diff", gdiff/normalize)
print(f"Color scale for {method} difference: {np.pi*normalize}")
grel = gdiff/data[f"gdc_mu_{ref}"]
normalize = np.nanmax(np.abs(grel))
export_img(f"convergence_{method}_relative", grel/normalize)
print(f"Color scale for {method} relative: {normalize}")
for method in ("o2", "o3", "o3a", "o3p", "o3a_ispline"):
g = np.pi*data[f"gdc_mu_{method}"]
with open(f"figdata/convergence_mu_{method}_dc.dat", "w") as file:
file.write("vac " + " ".join("gdc%.2f"%vdc[i] for i in dc_indices) + "\n")
np.savetxt(file, np.array([vac, *(g[:,i] for i in dc_indices)]).T)
with open(f"figdata/convergence_mu_{method}_ac.dat", "w") as file:
file.write("vdc " + " ".join("gdc%.2f"%vac[i] for i in ac_indices) + "\n")
np.savetxt(file, np.array([vdc, *(g[i] for i in ac_indices)]).T)
print("DC voltage: " + " ".join("%.6g"%vdc[i] for i in dc_indices))
print("AC voltage: " + " ".join("%.6g"%vac[i] for i in ac_indices))
def adjust_argument_array(start, end, *funcs, initial_resolution=51, s1=0.01, maxit=5, scale_coef=-0.5, yscales=1):
x = np.linspace(start, end, initial_resolution)
ys = [func(x) for func in funcs]
if isinstance(yscales, Number):
yscales = len(funcs)*[yscales]
else:
try:
assert len(yscales) == len(funcs)
except:
yscales = len(funcs)*[1]
for _ in range(maxit):
y1s = [(y[1:] - y[:-1]) / (x[1:] - x[:-1]) for y in ys]
w1s = [(y1[1:] - y1[:-1]) / (1 + yscale*(y1[1:] + y1[:-1])) for y1,yscale in zip(y1s,yscales)]
if scale_coef != 0:
for i in range(len(funcs)):
w1s[i] *= (x[2:] - x[:-2])**scale_coef
w1max = max(w1.max() for w1 in w1s)
settings.logger.debug(f"adjust argument: {w1max} {s1}, iteration {_+1}")
if w1max <= s1:
return (x, *ys)
x_selection = np.zeros(x.shape, dtype=bool)
for w1 in w1s:
x_selection[1:-1] |= np.abs(w1) > s1
x_selection[:-2] |= x_selection[1:-1]
#x_selection[2:-1] |= np.abs(w2) > s2
new_x = (x[x_selection] + x[1:][x_selection[:-1]])/2
full_x = np.append(x, new_x)
sorting = np.argsort(full_x)
x = full_x[sorting]
for i, func in enumerate(funcs):
ys[i] = np.append(ys[i], func(new_x))[sorting]
return (x, *ys)
def export_gapproximations_pgfplots(
omega = 16.5372,
vdc = 165.372,
vac1 = 16.5372,
vac2 = 165.372,
):
parameters = dict(
d = 1e9,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
xL = 0.5,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"]
)
dm = DataManager()
data_omega = dm.list(vac=vac1, vdc=0.0, method="J", voltage_branches=0, **parameters).sort_values("omega")
data_omega.drop_duplicates("omega", inplace=True)
data_vac = dm.list(omega=omega, vdc=0.0, method="J", voltage_branches=0, **parameters).sort_values("vac")
data_vac.drop_duplicates("vac", inplace=True)
data_vdc = dm.list(vac=vac2, omega=omega, method="mu", voltage_branches=4, **parameters).sort_values("vdc")
data_vdc.drop_duplicates("vdc", inplace=True)
data_vac_vdc = dm.list(vdc=vdc, omega=omega, method="J", voltage_branches=0, resonant_dc_shift=10, **parameters)
# messy selection of data points for data_vac_vdc
sel = (data_vac_vdc.padding>10) & (data_vac_vdc.nmax>=40) & ((data_vac_vdc.version_minor==15) | (data_vac_vdc.vac>=165.372))
data_vac_vdc = data_vac_vdc[sel].sort_values("vac")
data_vac_vdc.drop_duplicates("vac", keep="last", inplace=True)
interp = interp_vac0(dm, voltage_branches=4, method="mu", **parameters)
def adiabatic(vac, vdc):
if np.all(vdc == 0):
result, error = quad_vec((lambda x: interp(vac*np.sin(x))), 0, np.pi/2, epsrel=1e-7, epsabs=1e-14)
result *= 2
error *= 2
assert error < 1e-6
return result
result, error = quad_vec((lambda x: interp(vdc+vac*np.sin(x))), -np.pi/2, np.pi/2, epsrel=1e-7, epsabs=1e-14)
assert error < 1e-6
return result
def kaminski_high_vac(vac):
"""
Comparison to Kaminski et al PRB 62.8154 (2000)
Low frequency (adiabatic) limit
Eq. (77) in Kaminski et al 2000
"""
g = 3*np.pi**2/(16*np.log(vac/TK_VOLTAGE)**2)
if isinstance(g, Number):
return g
g[vac <= TK_VOLTAGE] = np.nan
g[g > 1.2] = np.nan
return g
def kaminski_high_omega(vac, omega):
"""
Comparison to Kaminski et al PRB 62.8154 (2000)
High frequency limit
Eqs. (72) - (74) in Kaminski et al 2000
"""
# Eq. (72) in Kaminski et al 2000
decoherence_rate = (vac/TK_VOLTAGE)**2 / (np.pi*omega/TK_VOLTAGE*np.log(omega/TK_VOLTAGE)**2)
# Eq. (74) in Kaminski et al 2000
g = 1 - decoherence_rate
# Eq. (73) in Kaminski et al 2000
g[decoherence_rate>1] = 3*np.pi**2/(16*np.log(decoherence_rate[decoherence_rate>1])**2)
return g
def kaminski_low_energy(vac):
"""
Comparison to Kaminski et al PRB 62.8154 (2000)
low energies, near equilibrium
Eq. (76) in Kaminski et al 2000
"""
return 1 - 3/16 * (vac/TK_VOLTAGE)**2
common_adjust_kw = dict(s1=1e-4, initial_resolution=101, maxit=4, scale_coef=0.9)
# first plot: swipe omega at constant vac and vdc=0
g_pat_func = lambda log_omega: np.pi*photon_assisted_tunneling(omega=np.exp(log_omega), vdc=0, vac=vac1, interp=interp)
gdc_func = interp1d(np.log(data_omega.omega), np.pi*data_omega.dc_conductance, kind="cubic")
log_omega_arr, g_pat_arr, gdc_arr = adjust_argument_array(
np.log(data_omega.omega.iloc[0]), np.log(data_omega.omega.iloc[-1]),
g_pat_func, gdc_func,
**common_adjust_kw,
yscales=(10,10))
omega_arr = np.exp(log_omega_arr)
g_adiabatic_arr = np.ones_like(omega_arr) * adiabatic(vac1, 0)
g_adiabatic_arr[omega_arr > 2.5*TK_VOLTAGE] = np.nan
g_kaminski1_arr = np.ones_like(omega_arr) * kaminski_high_vac(vac1)
g_kaminski1_arr[omega_arr > 2.5*TK_VOLTAGE] = np.nan
g_kaminski2_arr = kaminski_high_omega(vac1, omega_arr)
g_kaminski2_arr[np.exp(vac1/(np.pi*omega_arr*TK_VOLTAGE)**0.5)*TK_VOLTAGE > omega_arr] = np.nan
np.savetxt("figdata/omega_vac5_vdc0.dat",
np.array([omega_arr/TK_VOLTAGE,
gdc_arr,
g_adiabatic_arr,
g_pat_arr,
g_kaminski1_arr,
g_kaminski2_arr
]).T,
header="omega gdc gdc_adiabatic gdc_pat gdc_kaminski1 gdc_kaminski2",
comments = "",
fmt = "%.6g")
# second plot: swipe vac at constant omega and vdc=0
g_pat_func = lambda log_vac: np.pi*photon_assisted_tunneling(omega=omega, vac=np.exp(log_vac), vdc=0, interp=interp)
gdc_func = interp1d(np.log(data_vac.vac), np.pi*data_vac.dc_conductance, kind="cubic")
log_vac_arr, g_pat_arr, gdc_arr = adjust_argument_array(
np.log(data_vac.vac.iloc[0]), np.log(data_vac.vac.iloc[-1]),
g_pat_func, gdc_func,
**common_adjust_kw,
yscales=(10,10))
vac_arr = np.exp(log_vac_arr)
g_adiabatic_arr = adiabatic(vac_arr, 0)
g_kaminski1_arr = kaminski_high_vac(vac_arr)
g_kaminski1_arr[g_kaminski1_arr>=1] = np.nan
g_kaminski2_arr = kaminski_high_omega(vac_arr, omega)
g_kaminski2_arr[np.exp(vac_arr/(np.pi*omega*TK_VOLTAGE)**0.5)*TK_VOLTAGE > omega] = np.nan
np.savetxt("figdata/vac_omega5_vdc0.dat",
np.array([vac_arr/TK_VOLTAGE,
gdc_arr,
g_adiabatic_arr,
g_pat_arr,
g_kaminski1_arr,
g_kaminski2_arr
]).T,
header="vac gdc gdc_adiabatic gdc_pat gdc_kaminski1 gdc_kaminski2",
comments = "")
# third plot: swipe vdc at constant omega and vac
g_pat_func = lambda vdc: np.pi*photon_assisted_tunneling(omega=omega, vac=vac2, vdc=vdc, interp=interp)
gdc_func = interp1d(data_vdc.vdc, np.pi*data_vdc.dc_conductance, kind="cubic")
g_adiabatic_func = lambda vdc: adiabatic(vac2, vdc)
vdc_arr, g_pat_arr, gdc_arr, g_adiabatic_arr = adjust_argument_array(
data_vdc.vdc.iloc[0], data_vdc.vdc.iloc[-1],
g_pat_func, gdc_func, g_adiabatic_func,
**common_adjust_kw)
np.savetxt("figdata/vdc_omega5_vac10.dat",
np.array([vdc_arr/omega,
gdc_arr,
g_adiabatic_arr,
g_pat_arr,
]).T,
header="vdc gdc gdc_adiabatic gdc_pat",
comments = "")
# fourth plot: swipe vac at constant omega and vdc
g_pat_func = lambda vac: np.pi*photon_assisted_tunneling(omega=omega, vac=vac, vdc=vdc, interp=interp)
gdc_func = interp1d(data_vac_vdc.vac, np.pi*data_vac_vdc.dc_conductance, kind="cubic")
g_adiabatic_func = lambda vac: adiabatic(vac, vdc)
vac_arr, g_pat_arr, gdc_arr, g_adiabatic_arr = adjust_argument_array(
data_vac_vdc.vac.iloc[0], data_vac_vdc.vac.iloc[-1],
g_pat_func, gdc_func, g_adiabatic_func,
**common_adjust_kw)
np.savetxt("figdata/vac_omega5_vdc10.dat",
np.array([vac_arr/omega,
gdc_arr,
g_adiabatic_arr,
g_pat_arr,
]).T,
header="vac gdc gdc_adiabatic gdc_pat",
comments = "",
fmt = "%.6g")
def export_asymmetry_pgfplots(
omega = 16.5372,
vac = 26.45952,
xL = (0.5, 0.4, 0.3, 0.2, 0.1, 0.001),
vdc_max = 82.686,
):
dm = DataManager()
data = dm.list(
omega = omega,
vac = vac,
d = 1e9,
method = "mu",
voltage_branches = 4,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"]
).sort_values("vdc")
data.drop_duplicates(("vdc", "xL"), inplace=True)
gdc_funcs = []
for x in xL:
sel_data = data[np.abs(data.xL - x) < 1e-6]
gdc_funcs.append(interp1d(sel_data.vdc, np.pi/(4*x*(1-x))*sel_data.dc_conductance, kind="cubic"))
vdc_arr, *gdc_arrs = adjust_argument_array(0, vdc_max, *gdc_funcs, s1=1e-3, yscales=[0.1 for x in xL], initial_resolution=101, maxit=3)
np.savetxt("figdata/asymmetry.dat",
np.array([vdc_arr/omega,
*gdc_arrs
]).T,
header="vdc " + " ".join(f"xL{x:.3g}" for x in xL),
comments = "",
fmt = "%.6g")
def export_kogan04_pgfplots():
"""
In the paper the numbers Tk≃300mK and Ω≃2Tk are given.
These correspond to 7.4182 and 6.8851, respectively.
"""
#omega = 3.1
#xL = 0.05
omega = 3.5
xL = 0.166667
#omega = 4
#xL = 0.1
dm = DataManager()
# adjusted by factor 1.4:
vac_omega = np.array([0.72881, 1.13091, 1.50788, 1.6838, 3.61891])
data = dm.list(
omega = omega,
d = 1e9,
method = "mu",
voltage_branches = 4,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"]
).sort_values(["xL", "vac", "vdc"])
gdc_funcs = []
gdc_funcs_asym = []
#prefactor_sym = 0.113*np.pi
#prefactor_asym = 0.62*np.pi
#g_shift_sym = 0.023
#g_shift_asym = 0.022
#prefactor_sym = 0.1*np.pi
#g_shift_sym = 0.015
prefactor_sym = 0.096*np.pi
g_shift_sym = 0.02
#prefactor_sym = 0.0915*np.pi
#g_shift_sym = 0.0245
prefactor_asym = prefactor_sym / (4*xL*(1-xL))
g_shift_asym = g_shift_sym
for v in vac_omega:
data_sel = data[(np.abs(data.xL-0.5) < 1e-6) & (np.abs(data.vac-v*omega)<1e-6)]
try:
gdc_funcs.append(interp1d(data_sel.vdc, g_shift_sym + prefactor_sym*data_sel.dc_conductance, kind="cubic", bounds_error=False))
except ValueError:
gdc_funcs.append(lambda x: np.nan*np.empty_like(x))
data_sel = data[(np.abs(data.xL-xL) < 1e-6) & (np.abs(data.vac-v*omega)<1e-6)]
try:
gdc_funcs_asym.append(interp1d(data_sel.vdc, g_shift_asym + prefactor_asym*data_sel.dc_conductance, kind="cubic", bounds_error=False))
except ValueError:
gdc_funcs_asym.append(lambda x: np.nan*np.empty_like(x))
vdc, *arrays = adjust_argument_array(0, 25.5, *gdc_funcs, *gdc_funcs_asym, s1=5e-4, maxit=4, initial_resolution=51)
vdc = np.append(-vdc[:0:-1], vdc)
np.savetxt(
"figdata/kogan04.dat",
np.array([vdc/omega, *(np.append(a[:0:-1], a) for a in arrays)]).T,
header = "vdc " + " ".join(f"g{v:.3g}" for v in vac_omega) + " " + " ".join(f"g{v:.3g}asym" for v in vac_omega),
fmt = "%.6g",
comments = "")
def export_bruhat18_pgfplots():
dm = DataManager()
omega1 = 4.14
omega2 = 2.61
prefactor = 0.274*np.pi
g_shift = 0.0642
xL = 0.333333
prefactor_asym = prefactor/(4*xL*(1-xL))
# including correction factors:
vac_omega1 = np.array([0.90611, 1.132637, 1.359164, 1.585692, 1.812219, 2.038747, 2.265274, 2.491802, 2.718329, 3.397911])
vac_omega2 = np.array([0.431208, 0.862416, 1.293624, 1.724832, 2.15604, 2.587248, 3.018457, 3.449665, 3.880873, 4.743289, 6.468121])
parameters = dict(
d = 1e9,
method = "mu",
voltage_branches = 4,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"],
)
gdc_funcs1 = []
gdc_funcs2 = []
gdc_funcs1_asym = []
# Collect data
data1 = dm.list(omega=omega1, **parameters).sort_values(["xL", "vac", "vdc"])
for v in vac_omega1:
data_sel = data1[(np.abs(data1.xL-0.5) < 1e-6) & (np.abs(data1.vac-v*omega1)<1e-6)]
try:
gdc_funcs1.append(interp1d(data_sel.vdc, g_shift + prefactor*data_sel.dc_conductance, kind="cubic", bounds_error=False))
except ValueError:
gdc_funcs1.append(lambda x: np.nan*np.empty_like(x))
data_sel = data1[(np.abs(data1.xL-0.333333) < 1e-6) & (np.abs(data1.vac-v*omega1)<1e-6)]
try:
gdc_funcs1_asym.append(interp1d(data_sel.vdc, g_shift + prefactor_asym*data_sel.dc_conductance, kind="cubic", bounds_error=False))
except ValueError:
gdc_funcs1_asym.append(lambda x: np.nan*np.empty_like(x))
data2 = dm.list(omega=omega2, xL=0.5, **parameters).sort_values(["vac", "vdc"])
for v in vac_omega2:
data_sel = data2[np.abs(data2.vac-v*omega2)<1e-6]
try:
gdc_funcs2.append(interp1d(data_sel.vdc, g_shift + prefactor*data_sel.dc_conductance, kind="cubic", bounds_error=False))
except ValueError:
gdc_funcs2.append(lambda x: np.nan*np.empty_like(x))
vdc1, *arrays1 = adjust_argument_array(0, 21.5, *gdc_funcs1, s1=5e-4, maxit=4, initial_resolution=51)
vdc1_asym, *arrays1_asym = adjust_argument_array(0, 21.5, *gdc_funcs1_asym, s1=5e-4, maxit=4, initial_resolution=51)
vdc2, *arrays2 = adjust_argument_array(0, 21.5, *gdc_funcs2, s1=5e-4, maxit=4, initial_resolution=51)
vdc1 = np.append(-vdc1[:0:-1], vdc1) / omega1
vdc1_asym = np.append(-vdc1_asym[:0:-1], vdc1_asym) / omega1
vdc2 = np.append(-vdc2[:0:-1], vdc2) / omega2
# collect experimental data
#exp_data1 = np.genfromtxt("exp_data/KondoAC_Freq19GHz.dat", names=True)
#vdc_omega_exp1 = exp_data1["Vsd_mV"] * sc.eV/(19e12*sc.h)
#g_exp1_80 = exp_data1["Vac80"]
exp_data2 = np.genfromtxt("exp_data/KondoAC_Freq12GHz.dat", names=True)
vdc_omega_exp2 = exp_data2["Vsd_mV"] * sc.eV/(12e12*sc.h)
g_exp2_20 = exp_data2["Vac20"]
# get background
#g_exp1_bg = g_exp1_80 - gdc_funcs1[0](vdc_omega_exp1*omega1)
g_exp2_bg = g_exp2_20 - gdc_funcs2[0](np.abs(vdc_omega_exp2)*omega2)
sel = ~np.isnan(g_exp2_bg)
g_exp2_bg_spl = splrep(vdc_omega_exp2[sel], g_exp2_bg[sel], s=1e-2)
g2_adiabatic_bg = np.array([quad_vec((lambda t: splev(vdc2+vac*np.cos(t), g_exp2_bg_spl, ext=3)), 0, np.pi)[0]/np.pi for vac in vac_omega2])
g2 = np.array([np.append(a[:0:-1], a) for a in arrays2])
np.savetxt("figdata/bruhat18_19GHz.dat",
np.array([vdc1, *(np.append(a[:0:-1], a) for a in arrays1)]).T,
header = "vdc " + " ".join(f"g{v:.3g}" for v in vac_omega1),
fmt = "%.6g",
comments = "")
np.savetxt("figdata/bruhat18_12GHz.dat",
np.array([vdc2, *g2, *g2_adiabatic_bg, *(g2 + g2_adiabatic_bg)]).T,
header = "vdc " + " ".join(f"g{v:.3g}" for v in vac_omega2) + " " + " ".join(f"bg{v:.3g}" for v in vac_omega2) + " " + " ".join(f"gbg{v:.3g}" for v in vac_omega2),
fmt = "%.6g",
comments = "")
np.savetxt("figdata/bruhat18_19GHz_asym.dat",
np.array([vdc1_asym, *(np.append(a[:0:-1], a) for a in arrays1_asym)]).T,
header = "vdc " + " ".join(f"g{v:.3g}" for v in vac_omega1),
fmt = "%.6g",
comments = "")
def export_pulse_current_pgfplots(omega=1.5, pulse_duration=0.01):
dm = DataManager()
i_coef = []
g_coef = []
for phase in (0.25, 0.5, 0.75, 1):
vdc, fourier_coef = fourier_coef_gauss_symmetric(10, omega, pulse_duration, None, phase)
kondo, = dm.load_all(
omega = omega,
d = 1e9,
method = "J",
voltage_branches = 0,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
has_fourier_coef = True,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"],
fourier_coef = fourier_coef
)
nmax = kondo.nmax
i_coef.append(kondo.gammaL[nmax:,nmax]/TK_VOLTAGE)
g_coef.append(np.pi*kondo.deltaGammaL[nmax:,nmax])
del kondo
t_prefactor = 2*np.pi/omega * TK_VOLTAGE
t = np.concatenate((np.linspace(-0.01, -0.05/t_prefactor, 10, endpoint=False),
np.linspace(-0.05/t_prefactor, 0.08/t_prefactor, 130, endpoint=False),
np.linspace(0.08/t_prefactor, 0.5/t_prefactor, 42, endpoint=False),
np.linspace(0.5/t_prefactor, 3.33/t_prefactor, 40, endpoint=False),
np.linspace(3.33/t_prefactor, 0.49, 41)))
i_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in i_coef]
g_arrs = [-coef[0].real + 2*sum((np.exp(-2j*np.pi*t*n)*c).real for n,c in enumerate(coef)) for coef in g_coef]
t *= t_prefactor
array = np.array([t, *i_arrs, *g_arrs]).T
np.savetxt("figdata/pulse_current_full.dat",
array,
header = "t " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)),
fmt = "%.6g",
comments = "")
np.savetxt("figdata/pulse_current_zoom.dat",
array[10:141],
header = "t " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)),
fmt = "%.6g",
comments = "")
np.savetxt("figdata/pulse_current_tail.dat",
array[140:223],
header = "t " + " ".join(f"i{i}" for i in range(1,5)) + " " + " ".join(f"g{i}" for i in range(1,5)),
fmt = "%.6g",
comments = "")
def export_pulse_charge_pgfplots(omega=1.5):
dm = DataManager()
duration_data = [0.04]
phase_data = [0]
charge_data = [0]
trec_data = [np.nan]
_h5files = set()
parameters = dict(
d = 1e9,
method = "J",
voltage_branches = 0,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"],
pulse_type = "gauss_symmetric",
)
for (row, kondo) in itertools.chain(
load_all_pulses_full(dm, omega=omega, **parameters),
load_all_pulses_full(dm, omega=1.9, **parameters),
):
ref_time = min(4*row.pulse_duration*row.omega/(2*np.pi), 0.15)
nmax = kondo.nmax
charge_data.append(integrate_ft(-ref_time, 0.5-ref_time, kondo.gammaL[nmax-1::-1,nmax])/row.omega*2*np.pi)
_h5files.add(kondo._h5file)
duration_data.append(row.pulse_duration)
phase_data.append(row.pulse_phase)
nmax = kondo.nmax
g_coef = np.pi*kondo.deltaGammaL[nmax::-1,nmax]
g_coef[1:] *= 2
aux = 2j*np.pi*np.arange(g_coef.size)
func = lambda t: (g_coef*np.exp(t*aux)).real.sum() - 0.8
func_d1 = lambda t: (aux*g_coef*np.exp(t*aux)).real.sum()
func_d2 = lambda t: (aux**2*g_coef*np.exp(t*aux)).real.sum()
#func = lambda t: sum((c*np.exp(t*2j*np.pi*n)).real for n,c in enumerate(g_coef)) - 1
trec = np.nan
for init in np.array([1.2, 2, 4, 8, 16]) * row.pulse_duration:
try:
trec = newton(func, init, fprime=func_d1, fprime2=func_d2)
assert trec > 0.8*row.pulse_duration and trec < 0.45
break
except (RuntimeError, AssertionError):
pass
trec = np.nan
trec_data.append(TK_VOLTAGE*2*np.pi/row.omega * trec)
#print(trec, func(trec), charge_data[-1], row.pulse_phase, row.pulse_duration)
for file in _h5files:
file.close()
phase_data = np.array(phase_data)
duration_data = np.array(duration_data)
charge_data = np.array(charge_data)
trec_data = np.array(trec_data)
order = np.argsort(duration_data)
duration_data = duration_data[order]
charge_data = charge_data[order]
phase_data = phase_data[order]
trec_data = trec_data[order]
for phase in (0.25, 0.5, 0.75, 1):
sel = np.abs(phase_data - phase) < 1e-5
np.savetxt(f"figdata/pulse_charge_duration_p{phase:.2g}.dat",
np.array([duration_data[sel]*TK_VOLTAGE, charge_data[sel], trec_data[sel]]).T,
header = "t q trec",
fmt = "%.6g",
comments = "")
sel = np.abs(duration_data - 0.04) < 1e-6
np.savetxt("figdata/pulse_charge_phase.dat",
np.array([phase_data[sel], charge_data[sel], trec_data[sel]]).T,
header = "phase q trec",
fmt = "%.6g",
comments = "")
q_duration_funcs = []
trec_duration_funcs = []
for phase in (0.25, 0.5, 0.75, 1):
sel = np.abs(phase_data - phase) < 1e-5
q_duration_funcs.append(interp1d(duration_data[sel], charge_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan))
sel *= np.isfinite(trec_data)
trec_duration_funcs.append(interp1d(duration_data[sel], trec_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan))
t, *arrs = adjust_argument_array(0.01, 0.2, *q_duration_funcs, *trec_duration_funcs, s1=1e-3, maxit=3)
np.savetxt("figdata/pulse_charge_duration_interp.dat",
np.array([t*TK_VOLTAGE, *arrs]).T,
header = "t " + " ".join(f"q{i}" for i in range(1,5)) + " " + " ".join(f"trec{i}" for i in range(1,5)),
fmt = "%.6g",
comments = "")
sel = np.abs(duration_data - 0.04) < 1e-5
q_duration_func = interp1d(phase_data[sel], charge_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan)
sel *= np.isfinite(trec_data)
trec_duration_func = interp1d(phase_data[sel], trec_data[sel], kind="cubic", bounds_error=False, fill_value=np.nan)
phases, *arrs = adjust_argument_array(0.0, 2.8, q_duration_func, trec_duration_func, s1=1e-3, maxit=3)
np.savetxt("figdata/pulse_charge_phase_interp.dat",
np.array([phases, *arrs]).T,
header = "phase q trec",
fmt = "%.6g",
comments = "")
def export_harmonic_modes(
#omega = 3.4425351,
omega = 68.850702,
#omega = 16.5372,
n_phase_crossings = 10,
):
label = "-omega20"
vac_omega = np.array([5, 10, 20, 40, 80])
dm = DataManager()
data = dm.list(
omega = omega,
d = 1e9,
vdc = 0,
method = "J",
voltage_branches = 0,
include_Ga = True,
integral_method = -15,
solver_tol_rel = 1e-8,
solver_tol_abs = 1e-10,
bad_flags = DataManager.SOLVER_FLAGS["simplified_initial_conditions"] \
| DataManager.SOLVER_FLAGS["second_order_rg_equations"] \
| DataManager.SOLVER_FLAGS["solve_integral_exactly"] \
| DataManager.SOLVER_FLAGS["extrapolate_voltage"] \
| DataManager.SOLVER_FLAGS["deleted"],
)
current = []
t = np.linspace(-0.5, 0.5, 201)
current = []
phase = []
valid_vac_omega = []
phase_crossing_times = []
phase_crossing_currents = []
for vac in vac_omega:
sel = np.abs(data.vac-vac*omega) < 1e-6
if sel.sum() != 1:
continue
row = data[sel].iloc[-1]
filename = os.path.join(row.dirname, row.basename)
kondo, = KondoImport.read_from_h5(filename, row.hash)
nmax = row.nmax
coef = kondo.gammaL[nmax::-1,nmax]
iarr = 2*sum((c*np.exp(2j*np.pi*t*n)).real for n,c in enumerate(coef)) - coef.real[0]
iarr /= TK_VOLTAGE
current.append(iarr)
pc_times = np.arcsin(2*np.pi/vac*np.arange(n_phase_crossings)-1)/(2*np.pi)
pc_currents = 2*sum((c*np.exp(2j*np.pi*pc_times*n)).real for n,c in enumerate(coef)) - coef.real[0]
pc_currents /= TK_VOLTAGE
phase_crossing_times.append(pc_times)
phase_crossing_currents.append(pc_currents)
phase.append(vac/(2*np.pi)*(1+np.sin(2*np.pi*t)))
valid_vac_omega.append(vac)
coef /= coef[1]
np.savetxt(
f"figdata/harmonic_modes_vac{vac:.3g}{label}.dat",
np.array([np.arange(1, coef.size, 2), np.abs(coef[1::2])]).T,
header = "n i",
comments = "",
fmt = "%.6g")
np.savetxt(
f"figdata/harmonic_current{label}.dat",
np.array([t, *current, *phase]).T,
header = "t " + " ".join(f"i{vac:.3g}" for vac in valid_vac_omega) + " " + " ".join(f"p{vac:.3g}" for vac in valid_vac_omega),
comments = "",
fmt = "%.6g")
np.savetxt(
f"figdata/phase_crossings{label}.dat",
np.array([np.arange(n_phase_crossings), *phase_crossing_times, *phase_crossing_currents]).T,
header = "phase " + " ".join(f"t{vac:.3g}" for vac in valid_vac_omega) + " " + " ".join(f"i{vac:.3g}" for vac in valid_vac_omega),
comments = "",
fmt = "%.6g")
# adiabatic result
n = np.arange(1, 101, 2)
adiabatic_modes = []
data = dm.list(vac=0, omega=0, include_Ga=True, d=1e9, solver_tol_rel=1e-8, solver_tol_abs=1e-10, nmax=0, voltage_branches=4, xL=0.5, truncation_order=3, integral_method=-15)
data = data.sort_values("vdc")
vdc = np.concatenate((-data.vdc[:0:-1], data.vdc))
idc = np.concatenate((data.dc_current[:0:-1], data.dc_current))
interp = interp1d(vdc, idc, kind="cubic")
for i, vac in enumerate(omega*vac_omega):
ni, nierr = quad_vec((lambda t: np.sin(n*t) * interp(vac*np.sin(t))), 0, np.pi, limit=200)
adiabatic_modes.append(np.abs(ni/ni[0]))
np.savetxt(
f"figdata/harmonic_modes_adiabatic{label}.dat",
np.array([n, *adiabatic_modes]).T,
header = "n " + " ".join(f"mode{vac:.0f}" for vac in omega*vac_omega),
comments = "",
fmt = "%.6g")
fig, (ax1, ax2) = plt.subplots(nrows=2)
for i, (c, p, pt, pc) in enumerate(zip(current, phase, phase_crossing_times, phase_crossing_currents)):
color = f"C{i+1}"
ax1.plot(t, c, color=color)
ax1.plot(pt, pc, "x", color=color)
ax2.plot(p, c, color=color)
ax2.plot(np.arange(n_phase_crossings), pc, "x", color=color)
plt.show()
def export_convergence_pgfplots(omega = 16.5372):
"""
Plot convergence of current as function of D (Λ₀) for a fine
grid of dc and ac voltage at fixed omega.
"""
dm = DataManager()
assert isinstance(parameters["omega"], Number)
vac_num = 26
vdc_num = 26
vdc_max = 82.686
vac_max = 82.686
bad_flags = 0x0c84
good_flags = 0x1008
exponent = 2 if good_flags & 8 else 3
d_num = 5
log10d_min = 5
log10d_max = 9
d_arr = np.logspace(log10d_min, log10d_max, d_num)
gdc = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
idc = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
iac = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
phase = np.empty((d_num, vac_num, vdc_num), dtype=np.float64)
for i, d in enumerate(d_arr):
data = dm.list(d=d, **parameters)
vdc, vac, gdc[i], idc[i], iac[i], phase[i] = filter_grid_data(data, vac_max=vac_max, vdc_max=vdc_max, vac_num=vac_num, vdc_num=vdc_num)
gdc_fit_a = np.ndarray((vac_num, vdc_num), dtype=np.float64)
gdc_fit_b = np.ndarray((vac_num, vdc_num), dtype=np.float64)
idc_fit_a = np.ndarray((vac_num, vdc_num), dtype=np.float64)
idc_fit_b = np.ndarray((vac_num, vdc_num), dtype=np.float64)
iac_fit_a = np.ndarray((vac_num, vdc_num), dtype=np.float64)
iac_fit_b = np.ndarray((vac_num, vdc_num), dtype=np.float64)
j = np.vectorize(lambda d: solveTV0_scalar(d, rtol=1e-9, atol=1e-11)[2])
logd_inv3 = j(d_arr)**exponent
#logd_inv3 = np.log(d_arr)**-exponent
fit_func = lambda x, y: x[0] + x[1] * logd_inv3 - y
for i in range(vac_num):
for j in range(vdc_num):
(gdc_fit_a[i,j], gdc_fit_b[i,j]), trash = leastsq(fit_func, (gdc[-1,i,j], 0.), args=(gdc[:,i,j]))
(idc_fit_a[i,j], idc_fit_b[i,j]), trash = leastsq(fit_func, (idc[-1,i,j], 0.), args=(idc[:,i,j]))
(iac_fit_a[i,j], iac_fit_b[i,j]), trash = leastsq(fit_func, (iac[-1,i,j], 0.), args=(iac[:,i,j]))
extent = (-0.5*vdc_max/(vdc_num-1), vdc_max*(1+0.5/(vdc_num-1)), -0.5*vac_max/(vac_num-1), vac_max*(1+0.5/(vac_num-1)))
idc_diff = (abs(idc_fit_a) > 1e-6) * (((idc_fit_a.reshape((1,vac_num,vdc_num)) + idc_fit_b.reshape((1,vac_num,vdc_num))*logd_inv3.reshape((d_num,1,1)) - idc)**2).sum(axis=0)/d_num)**0.5 / (idc_fit_b * logd_inv3[-1])
iac_diff = (abs(iac_fit_a) > 1e-6) * (((iac_fit_a.reshape((1,vac_num,vdc_num)) + iac_fit_b.reshape((1,vac_num,vdc_num))*logd_inv3.reshape((d_num,1,1)) - iac)**2).sum(axis=0)/d_num)**0.5 / (iac_fit_b * logd_inv3[-1])
gdc_diff = (((gdc_fit_a.reshape((1,vac_num,vdc_num)) + gdc_fit_b.reshape((1,vac_num,vdc_num))*logd_inv3.reshape((d_num,1,1)) - gdc)**2).sum(axis=0)/d_num)**0.5 / (gdc_fit_b * logd_inv3[-1])
img = ax1.imshow(idc_fit_a, extent=extent, origin="lower")
img = ax2.imshow((abs(idc_fit_a)>1e-6)*idc_fit_b*(-logd_inv3[4])/idc_fit_a, extent=extent, origin="lower")
img = ax3.imshow(-idc_diff, extent=extent, origin="lower")
img = ax4.imshow(iac_fit_a, extent=extent, origin="lower")
img = ax5.imshow((abs(iac_fit_a)>1e-6)*iac_fit_b*(-logd_inv3[4])/iac_fit_a, extent=extent, origin="lower")
img = ax6.imshow(-iac_diff, extent=extent, origin="lower")
if __name__ == "__main__":
from sys import argv
globals()[argv[1]]()