Skip to content
Snippets Groups Projects
Commit fa6d96c8 authored by Valentin Bruch's avatar Valentin Bruch
Browse files

initial conditions for 2nd order truncation

parent 6477161d
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
kondo.py
+ 139
− 43
View file @ fa6d96c8
......@@ -34,6 +34,7 @@ from time import time, process_time
from scipy.special import jn as bessel_jn
from scipy.fftpack import fft
from scipy.integrate import solve_ivp
from scipy.optimize import newton
import settings
from rtrg import GlobalRGproperties, RGfunction
from reservoirmatrix import ReservoirMatrix, einsum_34_12_43, einsum_34_12_43_double, product_combinations
......@@ -109,6 +110,41 @@ def gen_init_matrix(nmax, *fourier_coef, resonant_dc_shift=0, resolution=5000):
return init_matrix
def solveTV0_scalar_order2(d, tk=0.32633176486110027, rtol=1e-8, atol=1e-10, full_output=False, **solveopts):
"""
Solve the ODE for second order truncated RG equations in
equilibrium at T=V=0 for scalars from 0 to d.
returns: (gamma, z, j, solver)
"""
# Initial conditions
j0 = 2/(np.pi*3**.5)
z0 = 1/(2*j0+1)
theta0 = 1/z0 - 1
gamma0 = tk * np.exp(1/(2*j0)) / j0
def ode_function_imaxis(lmbd, values):
"""RG eq. for Kondo model on the imaginary axis, ODE of functions of variable lmbd = Λ"""
gamma, theta, j = values
dgamma = theta
dtheta = -4*j**2/(lmbd + gamma)
dj = dtheta/2
return np.array([dgamma, dtheta, dj])
result = solve_ivp(
ode_function_imaxis,
(0, d),
np.array([gamma0, theta0, j0]),
t_eval = None if full_output else (d,),
rtol = rtol,
atol = atol,
**solveopts)
assert result.success
gamma, theta, j = result.y[:, -1]
z = 1/(1+theta)
return (gamma, z, j, result)
def solveTV0_scalar(d, tk=1, rtol=1e-8, atol=1e-10, full_output=False, **solveopts):
"""
Solve the ODE in equilibrium at T=V=0 for scalars from 0 to d.
......@@ -124,7 +160,7 @@ def solveTV0_scalar(d, tk=1, rtol=1e-8, atol=1e-10, full_output=False, **solveop
# Solve on imaginary axis
def ode_function_imaxis(lmbd, values):
'RG eq. for Kondo model on the imaginary axis, ODE of functions of variable lmbd = Λ'
"""RG eq. for Kondo model on the imaginary axis, ODE of functions of variable lmbd = Λ"""
gamma, theta, j = values
dgamma = theta
dtheta = -4*j**2/(lmbd + gamma)
......@@ -145,7 +181,8 @@ def solveTV0_scalar(d, tk=1, rtol=1e-8, atol=1e-10, full_output=False, **solveop
z = 1/(1+theta)
return (gamma, z, j, result)
def solveTV0_Utransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveopts):
def solveTV0_Utransformed(d, properties, tk=1, truncation_order=3, rtol=1e-8, atol=1e-10, **solveopts):
'''
Solve the ODE in equilibrium at T=V=0 to obtain initial conditions for Γ, Z and J.
Here all time-dependence is assumed to be contained in J.
......@@ -160,20 +197,28 @@ def solveTV0_Utransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveopt
'''
# Check input
assert isinstance(properties, GlobalRGproperties)
assert truncation_order in (2, 3)
# Solve equilibrium RG equations from 0 to d.
solveopts.update(rtol=rtol, atol=atol)
if truncation_order == 3:
gamma, z, j, scalar_solver = solveTV0_scalar(d, **solveopts)
elif truncation_order == 2:
gamma, z, j, scalar_solver = solveTV0_scalar_order2(d, **solveopts)
# Solve for constant imaginary part, go to required points in complex plane.
nmax = properties.nmax
vb = properties.voltage_branches
def ode_function_imconst(rE, values):
'RG eq. for Kondo model for constant Im(E), ODE of functions of rE = Re(E)'
"""RG eq. for Kondo model for constant Im(E), ODE of functions of rE = Re(E)"""
gamma, theta, j = values
dgamma = theta
dtheta = -4*j**2/(d - 1j*rE + gamma)
if truncation_order == 3:
dj = dtheta*(1 + j/(1 + theta))/2
elif truncation_order == 2:
dj = dtheta/2
return -1j*np.array([dgamma, dtheta, dj])
# Define flattened array of real parts of energy values, for which we want
......@@ -218,7 +263,8 @@ def solveTV0_Utransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveopt
return gamma, z, j
def solveTV0_untransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveopts):
def solveTV0_untransformed(d, properties, tk=1, truncation_order=3, rtol=1e-8, atol=1e-10, **solveopts):
'''
Solve the ODE in equilibrium at T=V=0 to obtain initial conditions for Γ, Z and J.
Here all time-dependence is assumed to be included in the Floquet
......@@ -231,9 +277,14 @@ def solveTV0_untransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveop
'''
# Check input
assert isinstance(properties, GlobalRGproperties)
assert truncation_order in (2, 3)
# Solve equilibrium RG equations from 0 to d.
solveopts.update(rtol=rtol, atol=atol)
if truncation_order == 3:
gamma, z, j, scalar_solver = solveTV0_scalar(d, **solveopts)
elif truncation_order == 2:
gamma, z, j, scalar_solver = solveTV0_scalar_order2(d, **solveopts)
# Solve for constant imaginary part, go to required points in complex plane.
nmax = properties.nmax
......@@ -244,7 +295,10 @@ def solveTV0_untransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveop
gamma, theta, j = values
dgamma = theta
dtheta = -4*j**2/(d - 1j*rE + gamma)
if truncation_order == 3:
dj = dtheta*(1 + j/(1 + theta))/2
elif truncation_order == 2:
dj = dtheta/2
return -1j*np.array([dgamma, dtheta, dj])
shifts = properties.mu.values + properties.omega*np.diag(np.arange(-nmax, nmax+1)).reshape((1,2*nmax+1,2*nmax+1))
......@@ -285,14 +339,17 @@ def solveTV0_untransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveop
gamma = np.einsum('kij,kj,klj->kil', eigvecs, gamma_raw, eigvecs.conjugate())
z = np.einsum('kij,kj,klj->kil', eigvecs, z_raw, eigvecs.conjugate())
j = np.einsum('kij,kj,klj->kil', eigvecs, j_raw, eigvecs.conjugate())
if truncation_order == 3:
zj_square = np.einsum('kij,kj,klj->kil', eigvecs, (j_raw*z_raw)**2, eigvecs.conjugate())
return gamma, z, j, zj_square
else:
j_square = np.einsum('kij,kj,klj->kil', eigvecs, j_raw**2, eigvecs.conjugate())
return gamma, z, j, j_square
class Kondo:
'''
"""
Kondo model with RG flow equations and routines for initial conditions.
Always accessible properties:
......@@ -325,7 +382,7 @@ class Kondo:
g2E : derivative of self.g2 with respect to E
g3E : derivative of self.g3 with respect to E
currentE : derivative of self.current with respect to E
'''
"""
def __init__(self,
unitary_transformation = True,
......@@ -444,9 +501,9 @@ class Kondo:
return getattr(self.global_properties, name)
def getParameters(self):
'''
"""
Get most relevant parameters. The returned dict can be used to label this object.
'''
"""
return {
'Ω' : self.omega,
'nmax' : self.nmax,
......@@ -465,7 +522,7 @@ class Kondo:
def initialize_untransformed(self,
**solveopts : 'keyword arguments passed to solver',
):
'''
"""
Arguments:
**solveopts: keyword arguments passed to the solver. Most relevant
are rtol and atol.
......@@ -473,14 +530,18 @@ class Kondo:
Get initial conditions for Γ, Z and G2 by numerically solving the
equilibrium RG equations from E=0 to E=iD and for all required Re(E).
Initialize G3, Iγ, δΓ, δΓγ.
'''
"""
sqrtxx = np.sqrt(self.xL*(1-self.xL))
symmetry = 0 if settings.IGNORE_SYMMETRIES else 1
#### Initial conditions from exact results at T=V=0
# Get Γ, Z and J (G2) for T=V=0.
gamma0, z0, j0, zj0_square = solveTV0_untransformed(d=self.d, properties=self.global_properties, **solveopts)
gamma0, z0, j0, zj0_square = solveTV0_untransformed(
d=self.d,
properties=self.global_properties,
truncation_order=self.truncation_order,
**solveopts)
# Create Γ and Z with the just calculated initial values.
self.gamma = RGfunction(self.global_properties, gamma0, symmetry=symmetry)
......@@ -501,7 +562,10 @@ class Kondo:
# G3 ~ Jtilde^2 with Jtilde = Z J
# Every entry of G3 will be of the following form (up to prefactors):
self.g3 = ReservoirMatrix(self.global_properties, symmetry=-symmetry)
g3_entry = RGfunction(self.global_properties, 1j*np.pi * zj0_square, symmetry=-symmetry)
g3_entry = RGfunction(
self.global_properties,
1j*np.pi * zj0_square,
symmetry=-symmetry)
self.g3[0,0] = 2*self.xL * g3_entry
self.g3[1,1] = 2*(1-self.xL) * g3_entry
g3_entry.symmetry = 0
......@@ -511,9 +575,12 @@ class Kondo:
## Initial conditions for current I^{γ=L} = J0 (1 - Jtilde)
# Note that j0[self.voltage_branches] and z0[self.voltage_branches] are diagonal Floquet matrices.
if self.truncation_order >= 3:
current_entry = 2*sqrtxx * j0[self.voltage_branches] * ( \
np.identity(2*self.nmax+1, dtype=np.complex128) \
- j0[self.voltage_branches] * z0[self.voltage_branches] )
else:
current_entry = 2*sqrtxx * j0[self.voltage_branches]
self.current = ReservoirMatrix(self.global_properties, symmetry=-symmetry)
self.current[0,0] = RGfunction(self.global_properties, np.zeros_like(current_entry), symmetry=-symmetry)
self.current[1,1] = RGfunction(self.global_properties, np.zeros_like(current_entry), symmetry=-symmetry)
......@@ -528,7 +595,7 @@ class Kondo:
)
## Initial conditions for voltage-variation of current-Γ: δΓ_L
# Note that j0[self.voltage_branches] and z0[self.voltage_branches] are diagonal Floquet matrices.
# Note that zj0_square[self.voltage_branches] is a diagonal Floquet matrix.
self.deltaGammaL = RGfunction(
self.global_properties,
3*np.pi*sqrtxx**2 * zj0_square[self.voltage_branches],
......@@ -553,7 +620,7 @@ class Kondo:
def initialize_Utransformed(self,
**solveopts : 'keyword arguments passed to solver',
):
'''
"""
Arguments:
**solveopts: keyword arguments passed to the solver. Most relevant
are rtol and atol.
......@@ -561,7 +628,7 @@ class Kondo:
Get initial conditions for Γ, Z and G2 by numerically solving the
equilibrium RG equations from E=0 to E=iD and for all required Re(E).
Initialize G3, Iγ, δΓ, δΓγ.
'''
"""
sqrtxx = np.sqrt(self.xL*(1-self.xL))
symmetry = 0 if settings.IGNORE_SYMMETRIES else 1
......@@ -569,7 +636,11 @@ class Kondo:
#### Initial conditions from exact results at T=V=0
# Get Γ, Z and J (G2) for T=V=0.
gamma0, z0, j0 = solveTV0_Utransformed(d=self.d, properties=self.global_properties, **solveopts)
gamma0, z0, j0 = solveTV0_Utransformed(
d = self.d,
properties = self.global_properties,
truncation_order = self.truncation_order,
**solveopts)
# Write T=V=0 results to Floquet index n=0.
gammavalues = np.zeros(self.shape(), dtype=np.complex128)
......@@ -587,6 +658,9 @@ class Kondo:
zvalues[diag_idx] = z0
jvalues[diag_idx] = j0
zj0 = z0 * j0 if self.truncation_order >= 3 else j0
zjvalues = zvalues * jvalues if self.truncation_order >= 3 else jvalues
RGclass = SymRGfunction if self.compact else RGfunction
# Create Γ and Z with the just calculated initial values.
......@@ -601,9 +675,15 @@ class Kondo:
if self.vac or self.resonant_dc_shift or self.fourier_coef is not None:
# Coefficients are given by the Bessel function of the first kind.
if self.fourier_coef is not None:
init_matrix = gen_init_matrix(self.nmax, *(f/self.omega for f in self.fourier_coef), resonant_dc_shift=self.resonant_dc_shift)
init_matrix = gen_init_matrix(
self.nmax,
*(f/self.omega for f in self.fourier_coef),
resonant_dc_shift = self.resonant_dc_shift)
else:
init_matrix = gen_init_matrix(self.nmax, self.vac/(2*self.omega), resonant_dc_shift=self.resonant_dc_shift)
init_matrix = gen_init_matrix(
self.nmax,
self.vac/(2*self.omega),
resonant_dc_shift = self.resonant_dc_shift)
j_LR = np.einsum('ij,...j->...ij', init_matrix, j0[...,2*self.resonant_dc_shift:])
j_RL = np.einsum('ji,...j->...ij', init_matrix.conjugate(), j0[...,:j0.shape[-1]-2*self.resonant_dc_shift])
j_LR = RGfunction(self.global_properties, j_LR)
......@@ -622,12 +702,12 @@ class Kondo:
# Every entry of G3 will be of the following form (up to prefactors):
self.g3 = ReservoirMatrix(self.global_properties, symmetry=-symmetry)
g3_entry = np.zeros(self.shape(), dtype=np.complex128)
g3_entry[diag_idx] = 1j*np.pi * (jvalues[diag_idx]*zvalues[diag_idx])**2
g3_entry[diag_idx] = 1j*np.pi * zjvalues[diag_idx]**2
g3_entry = RGclass(self.global_properties, g3_entry, symmetry=-symmetry)
self.g3[0,0] = 2*self.xL * g3_entry
self.g3[1,1] = 2*(1-self.xL) * g3_entry
if self.vac or self.resonant_dc_shift or self.fourier_coef is not None:
g30 = 1j*np.pi*(z0*j0)**2
g30 = 1j*np.pi*zj0**2
g3_LR = np.einsum('ij,...j->...ij', init_matrix, g30[...,2*self.resonant_dc_shift:])
g3_RL = np.einsum('ji,...j->...ij', init_matrix.conjugate(), g30[...,:g30.shape[-1]-2*self.resonant_dc_shift])
g3_LR = RGfunction(self.global_properties, g3_LR)
......@@ -643,9 +723,17 @@ class Kondo:
## Initial conditions for current I^{γ=L} = J0 (1 - Jtilde)
if self.voltage_branches:
current_entry = np.diag( 2*sqrtxx * jvalues[self.voltage_branches][diag_idx] * (1 - jvalues[self.voltage_branches][diag_idx] * zvalues[self.voltage_branches][diag_idx] ) )
if self.truncation_order >= 3:
current_entry = np.diag( \
2*sqrtxx * jvalues[self.voltage_branches][diag_idx] \
* (1 - zjvalues[self.voltage_branches][diag_idx] ) )
else:
current_entry = np.diag(2*sqrtxx * jvalues[self.voltage_branches][diag_idx])
else:
if self.truncation_order >= 3:
current_entry = np.diag(2*sqrtxx * jvalues[diag_idx] * (1 - zjvalues[diag_idx]))
else:
current_entry = np.diag( 2*sqrtxx * jvalues[diag_idx] * (1 - jvalues[diag_idx] * zvalues[diag_idx] ) )
current_entry = np.diag(2*sqrtxx * jvalues[diag_idx])
current_entry = RGclass(self.global_properties, current_entry, symmetry=-symmetry)
self.current = ReservoirMatrix(self.global_properties, symmetry=-symmetry)
if self.compact:
......@@ -658,9 +746,15 @@ class Kondo:
self.current[1,1] = self.current[0,0].copy()
if self.vac or self.resonant_dc_shift or self.fourier_coef is not None:
if self.voltage_branches:
i0 = j0[self.voltage_branches] * (1 - j0[self.voltage_branches]*z0[self.voltage_branches])
if self.truncation_order >= 3:
i0 = j0[self.voltage_branches] * (1 - zj0[self.voltage_branches])
else:
i0 = j0[self.voltage_branches]
else:
if self.truncation_order >= 3:
i0 = j0 * (1 - zj0)
else:
i0 = j0 * (1 - j0*z0)
i0 = j0
i_LR = np.einsum('ij,...j->...ij', init_matrix, i0[2*self.resonant_dc_shift:])
i_RL = np.einsum('ji,...j->...ij', init_matrix.conjugate(), i0[:i0.size-2*self.resonant_dc_shift])
i_LR = RGfunction(self.global_properties, i_LR)
......@@ -689,14 +783,16 @@ class Kondo:
if self.resonant_dc_shift:
assert self.compact == 0
if self.voltage_branches:
self.deltaGammaL.values[diag_idx] = 3*np.pi*sqrtxx**2 * (j0[self.voltage_branches,self.resonant_dc_shift:-self.resonant_dc_shift]*z0[self.voltage_branches,self.resonant_dc_shift:-self.resonant_dc_shift])**2
self.deltaGammaL.values[diag_idx] = 3*np.pi*sqrtxx**2 \
* zj0[self.voltage_branches,self.resonant_dc_shift:-self.resonant_dc_shift]**2
else:
self.deltaGammaL.values[diag_idx] = 3*np.pi*sqrtxx**2 * (j0[self.resonant_dc_shift:-self.resonant_dc_shift]*z0[self.resonant_dc_shift:-self.resonant_dc_shift])**2
self.deltaGammaL.values[diag_idx] = 3*np.pi*sqrtxx**2 \
* zj0[self.resonant_dc_shift:-self.resonant_dc_shift]**2
else:
if self.voltage_branches:
diag_values = 3*np.pi*sqrtxx**2 * (j0[self.voltage_branches]*z0[self.voltage_branches])**2
diag_values = 3*np.pi*sqrtxx**2 * zj0[self.voltage_branches]**2
else:
diag_values = 3*np.pi*sqrtxx**2 * (j0*z0)**2
diag_values = 3*np.pi*sqrtxx**2 * zj0**2
if self.compact:
assert diag_values.dtype == np.complex128
self.deltaGammaL.submatrix00 = np.diag(diag_values[0::2])
......@@ -717,9 +813,9 @@ class Kondo:
self.gammaL = self.vdc * self.deltaGammaL.reduced()
if self.vac and self.fourier_coef is None:
if self.voltage_branches:
gammaL_AC = 3*np.pi*sqrtxx**2 * self.vac/2 * (j0[self.voltage_branches]*z0[self.voltage_branches])**2
gammaL_AC = 3*np.pi*sqrtxx**2 * self.vac/2 * zj0[self.voltage_branches]**2
else:
gammaL_AC = 3*np.pi*sqrtxx**2 * self.vac/2 * (j0*z0)**2
gammaL_AC = 3*np.pi*sqrtxx**2 * self.vac/2 * zj0**2
if self.resonant_dc_shift:
gammaL_AC = gammaL_AC[...,self.resonant_dc_shift:-self.resonant_dc_shift]
idx = (np.arange(1, 2*self.nmax+1), np.arange(2*self.nmax))
......@@ -754,7 +850,7 @@ class Kondo:
save_iterations : 'number of iterations after which intermediate result should be saved' = 0,
**solveopts : 'keyword arguments passed to solver',
):
'''
"""
Initialize and solve the RG equations.
Arguments:
......@@ -774,7 +870,7 @@ class Kondo:
for Γ, Z, G2, G3, Iγ, δΓ, δΓγ.
Write parameters and solution for E=0 to self.<variables>
Return the ODE solver.
'''
"""
self.initialize(**solveopts)
self.save_filename = save_filename
......
settings.py
+ 2
− 2
View file @ fa6d96c8
......@@ -86,7 +86,7 @@ class GlobalFlags:
BASEPATH = os.path.abspath("data"),
DB_CONNECTION_STRING = "sqlite:///" + os.path.join(os.path.abspath("data"), "frtrg.sqlite"),
FILENAME = "frtrg-04.h5",
VERSION = (14,7,-1,-1),
VERSION = (14, 8, -1, -1),
MIN_VERSION = (14, 0),
LOG_TIME = 10, # in s
ENFORCE_SYMMETRIC = 0,
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment