Select Git revision
compact_rtrg.py
kondo.py 85.02 KiB
# Copyright 2022 Valentin Bruch <valentin.bruch@rwth-aachen.de>
# License: MIT
# type: ignore
"""
Kondo FRTRG, main module for RG calculations
Floquet real-time renormalization group implementation for the
spin 1/2 isotropic Kondo model.
Example usage:
>>> from kondo import Kondo, np
>>> nmax = 10
>>> vb = 3
>>> # Compute the RG flow in 2 different ways
>>> kondo1 = Kondo(
... unitary_transformation=True,
... omega=10,
... nmax=nmax,
... padding=8,
... vdc=4,
... vac=5,
... voltage_branches=vb)
>>> kondo2 = Kondo(
... unitary_transformation=False,
... omega=10,
... nmax=nmax,
... padding=0,
... vdc=4,
... vac=5,
... voltage_branches=vb)
>>> solver1 = kondo1.run()
>>> solver2 = kondo2.run()
>>> # Check if the results agree
>>> np.abs(kondo1.gammaL[:,nmax] - kondo2.gammaL[:,nmax]).max()
2.4691660784226606e-05
>>> np.abs(kondo1.deltaGammaL[:,nmax] - kondo2.deltaGammaL[:,nmax]).max()
1.2758585339583961e-05
>>> np.abs(kondo1.gamma[vb,:,nmax] - kondo2.gamma[vb,:,nmax]).max()
0.00018744930313729924
>>> np.abs(kondo1.z[vb,:,nmax] - kondo2.z[vb,:,nmax]).max()
1.421144031910071e-05
Further information:
https://arxiv.org/abs/2206.06263 and https://vbruch.eu/frtrg.html
"""
import hashlib
import numpy as np
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
from compact_rtrg import SymRGfunction
# Log times (global variables)
REF_TIME = time()
LAST_LOG_TIME = REF_TIME
def driving_voltage(tau, *fourier_coef):
"""
Generate function of time given Fourier coefficients.
tau = t/T so that 0 <= tau <= 1.
fourier_coef[n] = V_{n+1}/Ω
The result is given in the same units as fourier_coef.
"""
res = np.zeros_like(tau)
for n, c in enumerate(fourier_coef, 1):
res += (c * np.exp(2j*np.pi*n*tau)).real
return 2*res
def driving_voltage_integral(tau, *fourier_coef):
"""
Compute time-integral given Fourier coefficients.
tau = t/T so that 0 <= tau <= 1.
fourier_coef[n] = V_{n+1}/Ω
t
return Ω ∫dt' V(t') , t = tau*T
0
fourier_coef should be in units of Ω, then the result has
unit hbar/e = 1
"""
res = np.zeros_like(tau)
for n, c in enumerate(fourier_coef, 1):
res += (c/n * (np.exp(2j*np.pi*n*tau) - 1)).imag
return 2*res
def gen_init_matrix(nmax, *fourier_coef, resonant_dc_shift=0, resolution=5000):
"""
Generate Floquet matrix of the bare coupling without scalar
coupling prefactor.
fourier_coef must be in units of Ω:
fourier_coef[n] = V_{n+1}/Ω
TODO: signs have been chosen such that the result looks correct
"""
if len(fourier_coef) == 1 or np.allclose(fourier_coef[1:], 0):
# Simple driving, only one frequency:
coef = bessel_jn(
np.arange(-2*nmax+resonant_dc_shift, 2*nmax+resonant_dc_shift+1),
-2*fourier_coef[0])
elif len(fourier_coef) == 0:
coef = np.zeros(4*nmax+1)
coef[2*nmax-resonant_dc_shift] = 1
else:
# Anharmonic driving: use FFT
assert resolution > 2*nmax + abs(resonant_dc_shift)
assert resolution % 2 == 0
fft_coef = fft(np.exp(
-1j*driving_voltage_integral(
np.linspace(0, 1, resolution, endpoint=False),
*fourier_coef
)
))
coef = np.ndarray(4*nmax+1, dtype=np.complex128)
coef[2*nmax-resonant_dc_shift:] = fft_coef[:2*nmax+resonant_dc_shift+1].conjugate() / resolution
coef[:2*nmax-resonant_dc_shift] = fft_coef[-2*nmax+resonant_dc_shift:].conjugate() / resolution
init_matrix = np.ndarray((2*nmax+1, 2*nmax+1), dtype=np.complex128)
for j in range(2*nmax+1):
init_matrix[:,j] = coef[2*nmax-j:4*nmax+1-j]
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)
The default value of tk is adjusted to yield comparible results to
the case of 3rd order truncation for d=1e9, rtol=1e-8, atol=1e-10,
voltage_branches=4
"""
# 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 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.
returns: (gamma, z, j, solver)
"""
# Initial conditions
jbar = 2/(np.pi*np.sqrt(3))
j0 = jbar/(1-jbar)**2
# Θ := d/dΛ Γ = 1/Z - 1
theta0 = (1-jbar)**-2 - 1
gamma0 = np.sqrt((1-jbar)/jbar)*np.exp(1/(2*jbar)) * tk
# Solve on imaginary axis
def ode_function_imaxis(lmbd, values):
"""RG eq. for Kondo model on imaginary axis, ODE of functions of variable lmbd = Λ"""
gamma, theta, j = values
dgamma = theta
dtheta = -4*j**2/(lmbd + gamma)
dj = dtheta*(1 + j/(1 + theta))/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_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.
returns: (gamma, z, j)
Return values are arrays representing the quantities at energies shifted
by Ω and μ as required for the initial conditions.
If properties.resonant_dc_shift ≠ 0, then a larger array of energies
is considered, equivalent to mapping nmax → nmax+resonant_dc_shift in
the shape of properties.energies.
"""
# 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)"""
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
# to know, Γ, Z, J
nmax_plus = nmax + abs(properties.resonant_dc_shift)
if vb:
energies_orig = \
properties.energy.real \
+ properties.vdc*np.arange(-vb, vb+1).reshape((2*vb+1, 1)) \
+ properties.omega*np.arange(-nmax_plus, nmax_plus+1).reshape((1, 2*nmax_plus+1))
else:
energies_orig = properties.energy.real \
+ properties.omega * np.arange(-nmax_plus, nmax_plus+1)
energies = energies_orig.flatten()
energies_unique, inverse_indices = np.unique(energies, return_inverse=True)
if energies_unique.size == 1:
y = scalar_solver.y[:,-1].reshape((3,1))
else:
split_idx = np.searchsorted(energies_unique, 0)
energies_left = energies_unique[:split_idx]
energies_right = energies_unique[split_idx:]
result_left = solve_ivp(
ode_function_imconst,
t_span = (0, energies_left[0]),
y0 = np.array(scalar_solver.y[:,-1], dtype=np.complex128),
t_eval = energies_left[::-1],
**solveopts
)
assert result_left.success
result_right = solve_ivp(
ode_function_imconst,
t_span = (0, energies_right[-1]),
y0 = np.array(scalar_solver.y[:,-1], dtype=np.complex128),
t_eval = energies_right,
**solveopts
)
assert result_right.success
y = np.concatenate((result_left.y[:,::-1], result_right.y), axis=1)
gamma = y[0][inverse_indices].reshape(energies_orig.shape)
z = ( 1/(1 + y[1]) )[inverse_indices].reshape(energies_orig.shape)
j = y[2][inverse_indices].reshape(energies_orig.shape)
return gamma, z, j
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
matrix μ used for the voltage shift.
returns: (gamma, z, j, zj_square)
Return values represent Floquet matrices of the quantities at energies
shifted by μ as required for the initial conditions.
"""
# 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)'
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))
assert np.allclose(shifts.imag, 0)
eigvals, eigvecs = np.linalg.eigh(shifts)
assert all(np.allclose(eigvecs[i] @ np.diag(eigvals[i]) @ eigvecs[i].T.conjugate(), shifts[i]) for i in range(2*vb+1))
energies = eigvals.flatten()
energies_unique, inverse_indices = np.unique(energies, return_inverse=True)
if energies_unique.size == 1:
y = scalar_solver.y[:,-1].reshape((3,1))
else:
split_idx = np.searchsorted(energies_unique, 0)
energies_left = energies_unique[:split_idx]
energies_right = energies_unique[split_idx:]
result_left = solve_ivp(
ode_function_imconst,
t_span = (0, energies_left[0]),
y0 = np.array(scalar_solver.y[:,-1], dtype=np.complex128),
t_eval = energies_left[::-1],
**solveopts
)
assert result_left.success
result_right = solve_ivp(
ode_function_imconst,
t_span = (0, energies_right[-1]),
y0 = np.array(scalar_solver.y[:,-1], dtype=np.complex128),
t_eval = energies_right,
**solveopts
)
assert result_right.success
y = np.concatenate((result_left.y[:,::-1], result_right.y), axis=1)
gamma_raw = y[0][inverse_indices].reshape(eigvals.shape)
z_raw = ( 1/(1 + y[1]) )[inverse_indices].reshape(eigvals.shape)
j_raw = y[2][inverse_indices].reshape(eigvals.shape)
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:
total_iterations: total number of calls to self.updateRGequations()
global_properties: Properties for the RG functions, energy, ...
Properties stored in global_properties can be accessed
directly from self, e.g. using self.omega or self.energy.
When setting initial values, the following properties are added:
xL : asymmetry factor of the coupling, defaults to 0.5
d : UV cutoff
vac : AC voltage amplitude relative to Kondo temperature Tk
ir_cutoff : IR cutoff, should be 0, but is adapted when RG flow
is interrupted earlier
z : RGfunction, Z = 1/( 1 + i dΓ/dE )
gamma : RGfunction, Γ as in Π = 1/(E+iΓ)
deltaGammaL : RGfunction, δΓL (conductivity)
deltaGamma : RGfunction, δΓ
g2 : Reservoir matrix, coupling vertex G2
g3 : Reservoir matrix, coupling vertex G3
current : matrix in reservoir space fo RG functions,
representing the current vertex I^L
When running self.updateRGequations() the following properties are
added:
pi : RGfunction, Π = 1/(E+iΓ)
zE : derivative of self.z with respect to E
gammaE : derivative of self.gamma with respect to E
deltaGammaE : derivative of self.deltaGamma with respect to E
deltaGammaLE : derivative of self.deltaGammaL with respect to E
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,
nmax = 0,
padding = 0,
vdc = 0,
vac = 0,
omega = 0,
d = 1e9,
fourier_coef = None,
voltage_branches = 0,
include_Ga = False,
solve_integral_exactly = False,
resonant_dc_shift : 'DC bias voltages, multiples of Ω, positive int' = 0,
xL : 'asymmetry factor' = 0.5,
compact = 0,
simplified_initial_conditions = False,
truncation_order = 3,
**rg_properties):
"""
Create Kondo object, initialize global properties shared by all
Floquet matrices.
Expected arguments:
omega : frequency Ω, in units of Tk
nmax : size of Floquet matrix = (2*nmax+1, 2*nmax+1)
padding : extrapolation to avoid Floquet matrix truncation
effects, valid values: 0 ≤ padding ≤ 2*nmax-2
vdc : DC voltage, in units of Tk, including voltage due to
resonant_dc_shift.
vac : AC voltage, in units of Tk
voltage_branches : keep copies of Floquet matrices with energies
shifted by n Vdc, n = ±1,...,±voltage_branches.
Must be 0 or >=2
resonant_dc_shift : Describe DC voltage vdc partially by shifts
in the Floquet in the initial conditions.
valid values: non-negative integers
d : UV cutoff (convergence parameter)
include_Ga : include vertex paramter Ga in RG equations
solve_integral_exactly : diagonalize χ to solve integral in
next-to-leading order terms exactly.
xL = 1 - xR : asymmetry factor of coupling. Must fulfill 0 ≤ xL ≤ 1.
clear_corners : improve convergence for large Floquet matrices
by setting parts of the matrices to 0 after each
multiplication. Handle with care!
valid values: padding + clear_corners <= 2*nmax + 1
compact : Use extra symmetry to improve efficiency for large
matrices. compact != 0 requires unitary_transformation
and the symmetry V(t+(π/Ω)) = - V(t). Valid values are:
0: don't use compact form.
1: compact form that ignores matrix elements which
are zero by symmetry
2: compact form which additionally uses symmetry of
the nonzero matrix elements. requires xL = 0.5.
simplified_initial_conditions : use simplified initial
conditions for the current kernel which lead to the same
result in the limit of large D.
truncation_order : Truncation order of RG equations, must be 2 or 3.
"""
self.global_properties = GlobalRGproperties(
nmax = nmax,
omega = omega,
vdc = 0,
mu = None,
voltage_branches = voltage_branches,
resonant_dc_shift = resonant_dc_shift,
padding = padding,
fourier_coef = fourier_coef,
energy = 0j,
**rg_properties)
assert truncation_order in (2, 3)
self.truncation_order = truncation_order
self.global_settings = settings.export()
self.unitary_transformation = unitary_transformation
self.simplified_initial_conditions = simplified_initial_conditions
self.include_Ga = include_Ga
self.solve_integral_exactly = solve_integral_exactly
if truncation_order == 2:
assert not include_Ga
assert not solve_integral_exactly
self.compact = compact
self.vdc = vdc
self.vac = vac
self.xL = xL
if xL != 0.5:
self.global_properties.symmetric = False
assert self.compact != 2
# Some checks of the input
if (vac or fourier_coef is not None) and self.nmax == 0:
raise ValueError("Bad parameters: driving != 0 requires nmax > 0")
if xL < 0 or xL > 1:
raise ValueError("Bad parameter: need 0 <= xL <= 1")
if resonant_dc_shift and abs(resonant_dc_shift) > self.nmax:
raise ValueError("Bad parameters: resonant_dc_shift must be <= nmax")
if compact or voltage_branches == 0:
assert unitary_transformation
assert self.vdc == omega*resonant_dc_shift
if compact:
assert voltage_branches == 0
assert fourier_coef is None or np.allclose(fourier_coef[1::2], 0)
# When resonant_dc_shift is finite, typically the driving does
# not obey symmetry required for compact calculation.
# This implementation cannot handle resonant_dc_shift in
# combination with "compact" matrices.
assert self.resonant_dc_shift == 0
assert d.imag == 0.
self.d = d
if self.unitary_transformation:
self.compact = compact
self.global_properties.vdc = vdc - resonant_dc_shift*omega
else:
assert resonant_dc_shift == 0
mu = np.zeros((2*nmax+1, 2*nmax+1), dtype=np.complex128)
mu[np.diag_indices(2*nmax+1)] = vdc
if fourier_coef is None:
mu[np.arange(2*nmax), np.arange(1, 2*nmax+1)] = vac/2
mu[np.arange(1, 2*nmax+1), np.arange(2*nmax)] = vac/2
else:
for i, f in enumerate(fourier_coef, 1):
mu[np.arange(2*nmax+1-i), np.arange(i, 2*nmax+1)] = f
mu[np.arange(i, 2*nmax+1), np.arange(2*nmax+1-i)] = f.conjugate()
self.global_properties.mu = RGfunction(
self.global_properties,
np.arange(-voltage_branches, voltage_branches+1)\
.reshape((2*voltage_branches+1,1,1)) \
* mu.reshape((1,2*nmax+1,2*nmax+1)),
symmetry = -1)
self.total_iterations = 0
def __getattr__(self, name):
"""self.<name> is defined as shortcut for self.global_properties.<name>"""
return getattr(self.global_properties, name)
def getParameters(self):
"""
Get most relevant parameters.
The returned dict can be used to label this object.
"""
return {
'method' : 'J' if self.unitary_transformation else 'mu',
'Ω' : self.omega,
'nmax' : self.nmax,
'padding' : self.padding,
'Vdc' : self.vdc,
'Vac' : self.vac,
'V_branches' : self.voltage_branches,
'Vac' : getattr(self, 'vac', None),
'D' : getattr(self, 'd', None),
'solveopts' : getattr(self, 'solveopts', None),
'xL' : getattr(self, 'xL', None),
}
def initialize_untransformed(self,
**solveopts : 'keyword arguments passed to solver',
):
"""
Arguments:
**solveopts: keyword arguments passed to the solver. Most relevant
are rtol and atol.
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,
truncation_order=self.truncation_order,
**solveopts)
# Create Γ and Z with the just calculated initial values.
self.gamma = RGfunction(self.global_properties, gamma0, symmetry=symmetry)
self.z = RGfunction(self.global_properties, z0, symmetry=symmetry)
# Create G2 from J: G2_{ij} = - 2 sqrt(x_i x_j) J
self.g2 = ReservoirMatrix(self.global_properties, symmetry=symmetry)
j_rgfunction = RGfunction(self.global_properties, j0, symmetry=symmetry)
self.g2[0,0] = -2*self.xL * j_rgfunction
self.g2[1,1] = -2*(1-self.xL) * j_rgfunction
j_rgfunction.symmetry = 0
self.g2[0,1] = -2*sqrtxx * j_rgfunction
self.g2[1,0] = -2*sqrtxx * j_rgfunction
## Initial conditions for G3
# 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)
self.g3[0,0] = 2*self.xL * g3_entry
self.g3[1,1] = 2*(1-self.xL) * g3_entry
g3_entry.symmetry = 0
self.g3[0,1] = 2*sqrtxx * g3_entry
self.g3[1,0] = 2*sqrtxx * g3_entry
## 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)
self.current[0,1] = RGfunction(
self.global_properties,
current_entry,
symmetry=0)
self.current[1,0] = RGfunction(
self.global_properties,
-current_entry,
symmetry=0)
## Initial conditions for voltage-variation of Γ: δΓ
self.deltaGamma = RGfunction(
self.global_properties,
np.zeros((3,2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
symmetry = symmetry
)
## Initial conditions for voltage-variation of current-Γ: δΓ_L
# 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],
symmetry = symmetry
)
### Derivative of full current
self.yL = RGfunction(
self.global_properties,
np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
symmetry = -symmetry
)
### Full current, also includes AC current
self.gammaL = self.mu.reduced(shift=1) @ self.deltaGammaL
if self.simplified_initial_conditions:
self.gammaL *= 0
self.global_properties.energy = 1j*self.d
def initialize_Utransformed(self,
**solveopts : 'keyword arguments passed to solver',
):
"""
Initialize all RG objects with unitary transformation.
Arguments:
**solveopts: keyword arguments passed to the solver. Most
relevant are rtol and atol.
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). Then initialize G3, Iγ, δΓ, δΓγ, Yγ.
"""
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 = 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)
zvalues = np.zeros(self.shape(), dtype=np.complex128)
jvalues = np.zeros(self.shape(), dtype=np.complex128)
# construct diagonal matrices
diag_idx = (..., *np.diag_indices(2*self.nmax+1))
if self.resonant_dc_shift:
gammavalues[diag_idx] = gamma0[...,self.resonant_dc_shift:-self.resonant_dc_shift]
zvalues[diag_idx] = z0[...,self.resonant_dc_shift:-self.resonant_dc_shift]
jvalues[diag_idx] = j0[...,self.resonant_dc_shift:-self.resonant_dc_shift]
else:
gammavalues[diag_idx] = gamma0
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.
self.gamma = RGclass(self.global_properties, gammavalues, symmetry=symmetry)
self.z = RGclass(self.global_properties, zvalues, symmetry=symmetry)
# Create G2 from J: G2_{ij} = - 2 sqrt(x_i x_j) J
self.g2 = ReservoirMatrix(self.global_properties, symmetry=symmetry)
j_rgfunction = RGclass(self.global_properties, jvalues, symmetry=symmetry)
self.g2[0,0] = -2*self.xL * j_rgfunction
self.g2[1,1] = -2*(1-self.xL) * j_rgfunction
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)
else:
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)
j_RL = RGfunction(self.global_properties, j_RL)
self.g2[0,1] = -2*sqrtxx * j_LR
self.g2[1,0] = -2*sqrtxx * j_RL
else:
assert self.compact == 0
j_rgfunction.symmetry = 0
self.g2[0,1] = -2*sqrtxx * j_rgfunction
self.g2[1,0] = -2*sqrtxx * j_rgfunction
## Initial conditions for G3
# 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 = np.zeros(self.shape(), dtype=np.complex128)
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*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)
g3_RL = RGfunction(self.global_properties, g3_RL)
self.g3[0,1] = 2*sqrtxx * g3_LR
self.g3[1,0] = 2*sqrtxx * g3_RL
else:
assert self.compact == 0
g3_entry.symmetry = 0
self.g3[0,1] = 2*sqrtxx * g3_entry
self.g3[1,0] = 2*sqrtxx * g3_entry
## Initial conditions for current I^{γ=L} = J0 (1 - Jtilde)
if self.voltage_branches:
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])
current_entry = RGclass(self.global_properties, current_entry, symmetry=-symmetry)
self.current = ReservoirMatrix(self.global_properties, symmetry=-symmetry)
if self.compact:
self.current[0,0] = RGclass(self.global_properties, None, symmetry=symmetry)
self.current[0,0].submatrix01 = np.zeros((self.nmax+1, self.nmax), np.complex128)
self.current[0,0].submatrix10 = np.zeros((self.nmax, self.nmax+1), np.complex128)
else:
self.current[0,0] = 0*current_entry
self.current[0,0].symmetry = -symmetry
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:
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
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)
i_RL = RGfunction(self.global_properties, i_RL)
self.current[0,1] = 2*sqrtxx * i_LR
self.current[1,0] = -2*sqrtxx * i_RL
else:
assert self.compact == 0
current_entry.symmetry = 0
self.current[0,1] = 2*sqrtxx * current_entry
self.current[1,0] = -2*sqrtxx * current_entry
## Initial conditions for voltage-variation of Γ: δΓ
self.deltaGamma = RGclass(
self.global_properties,
None if self.compact else np.zeros(
(3,2*self.nmax+1,2*self.nmax+1) if self.voltage_branches else self.shape(),
dtype=np.complex128),
symmetry = symmetry
)
## Initial conditions for voltage-variation of current-Γ: δΓ_L
self.deltaGammaL = RGclass(
self.global_properties,
None if self.compact else np.zeros((2*self.nmax+1, 2*self.nmax+1), dtype=np.complex128),
symmetry = symmetry
)
if self.resonant_dc_shift:
assert self.compact == 0
if self.voltage_branches:
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 \
* zj0[self.resonant_dc_shift:-self.resonant_dc_shift]**2
else:
if self.voltage_branches:
diag_values = 3*np.pi*sqrtxx**2 * zj0[self.voltage_branches]**2
else:
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])
self.deltaGammaL.submatrix11 = np.diag(diag_values[1::2])
else:
self.deltaGammaL.values[diag_idx] = diag_values
del diag_values
### Derivative of full current
self.yL = RGclass(
self.global_properties,
None if self.compact else np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
symmetry=-symmetry
)
### Full current, also includes AC current
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 * zj0[self.voltage_branches]**2
else:
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))
self.gammaL.values[idx] = gammaL_AC[...,1:]
idx = (np.arange(0, 2*self.nmax), np.arange(1, 2*self.nmax+1))
self.gammaL.values[idx] = gammaL_AC[...,:-1]
if self.simplified_initial_conditions:
self.gammaL *= 0
if self.compact:
self.deltaGamma.submatrix01 = np.zeros((self.nmax+1, self.nmax), dtype=np.complex128)
self.deltaGamma.submatrix10 = np.zeros((self.nmax, self.nmax+1), dtype=np.complex128)
self.yL.submatrix01 = np.zeros((self.nmax+1, self.nmax), dtype=np.complex128)
self.yL.submatrix10 = np.zeros((self.nmax, self.nmax+1), dtype=np.complex128)
self.gammaL.submatrix00 = None
self.gammaL.submatrix11 = None
self.gammaL.submatrix01 = np.zeros((self.nmax+1, self.nmax), dtype=np.complex128)
self.gammaL.submatrix10 = np.zeros((self.nmax, self.nmax+1), dtype=np.complex128)
self.global_properties.energy = 1j*self.d
def initialize(self, **solveopts):
if self.unitary_transformation:
self.initialize_Utransformed(**solveopts)
else:
self.initialize_untransformed(**solveopts)
def run(self,
ir_cutoff : 'IR cutoff of RG flow' = 0,
forget_flow : 'do not store RG flow' = True,
save_filename : 'save intermediate results: string containing %d for number of iterations' = '',
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:
ir_cutoff: Stop the RG flow at Λ = -iE = ir_cutoff (instead of
Λ=0). If the RG flow is interrupted earlier, ir_cutoff
will be adapted.
**solveopts: keyword arguments passed to the solver. Most
interesting are rtol and atol.
1. 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γ, δΓ, δΓγ, Yγ.
2. Solve RG equations from E=iD to E=0
for Γ, Z, G2, G3, Iγ, δΓ, δΓγ, Yγ.
Write parameters and solution for E=0 to self.<variables>
Return the ODE solver.
"""
self.initialize(**solveopts)
self.save_filename = save_filename
self.save_iterations = save_iterations
if ir_cutoff:
self.ir_cutoff = ir_cutoff
self.solveopts = solveopts
if ir_cutoff >= self.d:
return
### Solve RG ODE
output = self.solveOdeIm(self.d, ir_cutoff, only_final=forget_flow, **solveopts)
# Write final values to Floquet matrices in self.
try:
# Shift energy
self.global_properties.energy = self.energy.real + 1j*output.t[-1]
# Unpack values
self.unpackFlattenedValues(output.y[:,-1])
except:
settings.logger.exception("Failed to read solver results:")
return output
def updateRGequations(self):
"""
Calculates the energy derivatives using the RG equations.
The derivative of self.<name> is written to self.<name>E.
A human readable reference implementation for this function
is provided in updateRGequations_reference.
"""
if settings.ENFORCE_SYMMETRIC:
if hasattr(self, 'pi'):
assert self.pi.symmetry == -1
assert self.z.symmetry == 1
assert self.yL.symmetry == -1
assert self.gamma.symmetry == 1
assert self.gammaL.symmetry == 1
assert self.deltaGamma.symmetry == 1
assert self.deltaGammaL.symmetry == 1
assert self.g2[0,0].symmetry == 1
assert self.g2[1,1].symmetry == 1
assert self.g3[0,0].symmetry == -1
assert self.g3[1,1].symmetry == -1
assert self.current[0,0].symmetry == -1
assert self.current[1,1].symmetry == -1
# Print some log message to indicate progress
global LAST_LOG_TIME, REF_TIME
if settings.LOG_TIME > 0 and time() - LAST_LOG_TIME >= settings.LOG_TIME:
LAST_LOG_TIME = time()
settings.logger.info(
"%9.2fs: Λ = %.4e, iterations = %d"%(
process_time(), self.energy.imag, self.total_iterations))
if settings.USE_REFERENCE_IMPLEMENTATION:
return self.updateRGequations_reference()
symmetry = 0 if settings.IGNORE_SYMMETRIES else 1
if self.include_Ga:
raise NotImplementedError
if self.solve_integral_exactly:
raise NotImplementedError
# Denote costs in terms of Floquet matrix products as x/y/z where
# x is the number of multiplications for resonant_dc_shift != 0 without symmetry.
# y is the number of multiplications for xL != xR but resonant_dc_shift == 0,
# z is the number of multiplications for xL == xR,
# Total costs: (+ 2 inversions, optionally + 2 multiplications for pi.derivative)
# 178 / 116 / 67
## RG eq for Γ
zinv = self.z.inverse() # costs: 1 inversion
self.gammaE = -1j*( zinv - (SymRGfunction if self.compact else RGfunction)(self.global_properties, 'identity') )
del zinv
# Derivative of G2
# First calculate Π (at E and shifted by multiples of vdc or mu):
self.pi = (-1j*self.gamma).k2lambda(self.mu) # costs: 1 inversion
# first terms of g2E:
# 1/2 G13 Π G32 ,
# 1/2 G32 Π G13
g2_pi = self.g2 @ self.pi # costs: 4/3/2
g2E1, g2E2 = product_combinations( g2_pi, self.g2 ) # costs: 12/7/4
# g2E1.tr() = -i (d/dE)² Γ
g2E1tr = g2E1.tr()
g2E1tr.symmetry = -symmetry
if self.compact and not isinstance(g2E1tr, SymRGfunction):
g2E1tr = SymRGfunction.fromRGfunction(g2E1tr, diag=True, offdiag=False)
g2E1 *= 0.5
g2E2 *= 0.5
self.zE = self.z @ g2E1tr @ self.z # costs: 2/2/2
del g2E1tr
## RG eq for G2
self.g2E = g2E1 + g2E2
if self.truncation_order >= 3:
# third term for g2E
# -1/4 G34 ( Π G12 Z + Z G12 Π ) G43
# First the part inside the brackets:
pi_g2_z = self.pi @ self.g2 @ self.z + self.z @ g2_pi # costs: 12/9/6
g2_bracket_g2, g2_bracket_g3 = einsum_34_12_43_double(self.g2, pi_g2_z, self.g2, self.g3) # costs: 48/30/15
## RG eq for G2
self.g2E += (-0.25)*g2_bracket_g2
del g2_bracket_g2
self.g2E.symmetry = -symmetry
self.g2E[0,0].symmetry = -symmetry
self.g2E[1,1].symmetry = -symmetry
self.deltaGammaE = 1j * (g2E1[0,0] - g2E1[1,1] + g2E2[1,1] - g2E2[0,0]).reduced_to_voltage_branches(1)
self.deltaGammaE.symmetry = -symmetry
del g2E1, g2E2
# first terms of G3E:
# G2_13 Π G3_32 ,
# G2_32 Π G3_13
g3E1, g3E2 = product_combinations( g2_pi, self.g3 ) # costs: 12/7/4
del g2_pi
self.g3E = g3E1 + g3E2
if self.truncation_order >= 3:
# third term for g3E
# 1/2 G2_34 ( Π G2_12 Z + Z G2_12 Π ) G3_43
# -> already calculated
self.g3E += 0.5 * g2_bracket_g3
del g2_bracket_g3
self.g3E.symmetry = symmetry
self.g3E[0,0].symmetry = symmetry
self.g3E[1,1].symmetry = symmetry
del g3E1, g3E2
# first terms of iE:
# I13 Π G32 ,
# I32 Π G13
i_pi = self.current @ self.pi # costs: 4/3/2
i_pi.symmetry = 0
if settings.ENFORCE_SYMMETRIC:
assert i_pi[0,0].symmetry == 1
assert i_pi[1,1].symmetry == 1
iE1, iE2 = product_combinations( i_pi, self.g2 ) # costs: 12/7/4
self.currentE = iE1 + iE2
del iE1, iE2
if self.truncation_order >= 3:
# third term for iE
# 1/2 I34 ( Π G2_12 Z + Z G2_12 Π ) G43
# The part in the brackets has been calculated before.
self.currentE += 0.5 * einsum_34_12_43( self.current, pi_g2_z, self.g2 ) # costs: 32/20/10
del pi_g2_z
self.currentE.symmetry = symmetry
self.currentE[0,0].symmetry = symmetry
self.currentE[1,1].symmetry = symmetry
# RG equation for δΓ_L
# First part: (δ1L - δ2L) I_12 Π G^3_21
i_pi_g3_01 = i_pi[0,1] @ self.g3[1,0] # costs: 1
i_pi_g3_10 = i_pi[1,0] @ self.g3[0,1] # costs: 1
self.deltaGammaLE = 1.5j*(i_pi_g3_01 - i_pi_g3_10)
# Second part: I_12 Z Π δΓ G^3_21
if self.truncation_order >= 3:
z_reduced = self.z.reduced_to_voltage_branches(1)
dGamma_reduced = self.deltaGamma.reduced_to_voltage_branches(1)
pi_reduced = self.pi.reduced_to_voltage_branches(1)
z_dgamma_pi = z_reduced @ dGamma_reduced @ pi_reduced \
+ pi_reduced @ dGamma_reduced @ z_reduced # costs: 4/4/4
del z_reduced, pi_reduced, dGamma_reduced
if self.compact:
assert isinstance(z_dgamma_pi, SymRGfunction)
if settings.ENFORCE_SYMMETRIC:
assert z_dgamma_pi.symmetry == -1
self.deltaGammaLE += (-0.75)*einsum_34_12_43( self.current, z_dgamma_pi, self.g3 ) # costs: 32/20/10
del z_dgamma_pi
self.deltaGammaLE.symmetry = -symmetry
# RG equation for ΓL
self.gammaLE = self.yL
self.yLE = 1.5j*(i_pi[0,0] @ self.g3[0,0] + i_pi[1,1] @ self.g3[1,1] + i_pi_g3_01 + i_pi_g3_10) # costs: 2/2/2
self.yLE.symmetry = symmetry
del i_pi, i_pi_g3_10, i_pi_g3_01
if self.compact:
if not isinstance(self.deltaGammaE, SymRGfunction):
self.deltaGammaE = SymRGfunction.fromRGfunction(self.deltaGammaE, diag=False, offdiag=True)
if not isinstance(self.deltaGammaLE, SymRGfunction):
self.deltaGammaLE = SymRGfunction.fromRGfunction(self.deltaGammaLE, diag=True, offdiag=False)
if not isinstance(self.g2E[0,0], SymRGfunction):
self.g2E[0,0] = SymRGfunction.fromRGfunction(self.g2E[0,0], diag=True, offdiag=False)
if not isinstance(self.g2E[1,1], SymRGfunction):
self.g2E[1,1] = SymRGfunction.fromRGfunction(self.g2E[1,1], diag=True, offdiag=False)
if not isinstance(self.g3E[0,0], SymRGfunction):
self.g3E[0,0] = SymRGfunction.fromRGfunction(self.g3E[0,0], diag=True, offdiag=False)
if not isinstance(self.g3E[1,1], SymRGfunction):
self.g3E[1,1] = SymRGfunction.fromRGfunction(self.g3E[1,1], diag=True, offdiag=False)
if not isinstance(self.currentE[0,0], SymRGfunction):
self.currentE[0,0] = SymRGfunction.fromRGfunction(self.currentE[0,0], diag=False, offdiag=True)
if not isinstance(self.currentE[1,1], SymRGfunction):
self.currentE[1,1] = SymRGfunction.fromRGfunction(self.currentE[1,1], diag=False, offdiag=True)
if not isinstance(self.yLE, SymRGfunction):
self.yLE = SymRGfunction.fromRGfunction(self.yLE, diag=False, offdiag=True)
# Count calls to RG equations.
self.total_iterations += 1
def updateRGequations_reference(self):
"""
Reference implementation of updateRGequations without optimization.
This reference implementation serves as a check for the more efficient
function updateRGequations.
This function takes approximately twice as long as updateRGequations.
"""
if self.include_Ga:
raise NotImplementedError
if self.solve_integral_exactly:
raise NotImplementedError
# Notation (mainly allow using objects without "self")
z = self.z
gamma = self.gamma
deltaGamma = self.deltaGamma
deltaGammaL = self.deltaGammaL
g2 = self.g2
g3 = self.g3
il = self.current
yL = self.yL
# Identity matrix
identity = RGfunction(self.global_properties, 'identity')
# Resolvent
pi = self.pi = (-1j*gamma).k2lambda(self.mu)
# Compute the sum
# Σ A B C
# 1,2 12 21
einsum_34_x_43 = lambda a, b, c: \
a[0,0] @ b @ c[0,0] \
+ a[0,1] @ b @ c[1,0] \
+ a[1,0] @ b @ c[0,1] \
+ a[1,1] @ b @ c[1,1]
### RG equations in human readable form
# Note that the shifts in the energy arguments of all RGfunction
# objects is implicitly handled by the multiplication operators.
# The muliplication operators for ReservoirMatrix objects are:
# @ for normal matrix multiplication (matrix indices ij, jk → ik with sum over j)
# % for transpose matrix multiplication (matrix indices jk, ij → ik with sum over j)
# ⎛ ⊤ ⊤⎞⊤
# i.e. A % B = ⎜A @ B ⎟ when there are no energy argument shifts.
# ⎝ ⎠
# dΓ ⎛ 1 ⎞
# —— = -i ⎜ — - 1 ⎟
# dE ⎝ Z ⎠
self.gammaE = -1j*(z.inverse() - identity)
# dZ
# —— = Z tr( G² Π G² ) Z
# dE
self.zE = z @ (g2 @ pi @ g2).tr() @ z
if self.truncation_order == 3:
# bracket = Π G² Z + Z G² Π
bracket = pi @ g2 @ z + z @ g2 @ pi
# dG² 1 1 ⎛ ⊤ ⊤⎞⊤ 1 ⎛ ⎞
# ——— = — G² Π G² + — ⎜G² Π G² ⎟ - — G² ⎜ Π G² Z + Z G² Π ⎟ G²
# dE 2 2 ⎝ ⎠ 4 34 ⎝ ⎠ 43
self.g2E = .5 * g2 @ pi @ g2 + .5 * ((g2 @ pi) % g2) - .25 * einsum_34_x_43(g2, bracket, g2)
# dG³ ⎛ ⊤ ⊤⎞⊤ 1 ⎛ ⎞
# ——— = G² Π G³ + ⎜G² Π G³ ⎟ + — G² ⎜ Π G² Z + Z G² Π ⎟ G³
# dE ⎝ ⎠ 2 34 ⎝ ⎠ 43
self.g3E = g2 @ pi @ g3 + ((g2 @ pi) % g3) + .5 * einsum_34_x_43(g2, bracket, g3)
# γ
# dI γ ⎛ γ⊤ ⊤⎞⊤ 1 γ ⎛ ⎞
# ——— = I Π G² + ⎜I Π G² ⎟ + — I ⎜ Π G² Z + Z G² Π ⎟ G³
# dE ⎝ ⎠ 2 34 ⎝ ⎠ 43
self.currentE = il @ pi @ g2 + ((il @ pi) % g2) + .5 * einsum_34_x_43(il, bracket, g2)
elif self.tuncation_order == 2:
# dG² 1 1 ⎛ ⊤ ⊤⎞⊤
# ——— = — G² Π G² + — ⎜G² Π G² ⎟
# dE 2 2 ⎝ ⎠
self.g2E = .5 * g2 @ pi @ g2 + .5 * ((g2 @ pi) % g2)
# dG³ ⎛ ⊤ ⊤⎞⊤
# ——— = G² Π G³ + ⎜G² Π G³ ⎟
# dE ⎝ ⎠
self.g3E = g2 @ pi @ g3 + ((g2 @ pi) % g3)
# γ
# dI γ ⎛ γ⊤ ⊤⎞⊤
# ——— = I Π G² + ⎜I Π G² ⎟
# dE ⎝ ⎠
self.currentE = il @ pi @ g2 + ((il @ pi) % g2)
else:
raise ValueError("Invalid truncation order: %s"%self.truncation_order)
# dδΓ ⎛ ⎞
# ——— = i ⎜ δ - δ ⎟ G² Π G²
# dE ⎝ 1L 2L ⎠ 12 21
deltaGammaE = 1j * (g2[0,1] @ pi @ g2[1,0] - g2[1,0] @ pi @ g2[0,1])
# Reduction of voltage branches as required for the solver.
# In this step some information is thrown away that cannot affect the
# result of the physical observables.
self.deltaGammaE = deltaGammaE.reduced_to_voltage_branches(1)
pi_reduced = pi.reduced_to_voltage_branches(1)
z_reduced = z.reduced_to_voltage_branches(1)
if self.truncation_order == 3:
# γ
# dδΓ 3 ⎛ ⎞ γ 3 ⎛ γ⎛ ⎞ ⎞
# ———— = — i ⎜ δ - δ ⎟ I Π G³ - — tr ⎜I ⎜ Π δΓ Z + Z δΓ Π ⎟ G³⎟
# dE 2 ⎝ 1L 2L ⎠ 12 21 4 ⎝ ⎝ ⎠ ⎠
self.deltaGammaLE = 1.5j * (il[0,1] @ pi @ g3[1,0] - il[1,0] @ pi @ g3[0,1]) \
- 0.75 * (il @ (pi_reduced @ deltaGamma @ z_reduced + z_reduced @ deltaGamma @ pi_reduced) @ g3).tr()
elif self.truncation_order == 2:
# γ
# dδΓ 3 ⎛ ⎞ γ
# ———— = — i ⎜ δ - δ ⎟ I Π G³
# dE 2 ⎝ 1L 2L ⎠ 12 21
self.deltaGammaLE = 1.5j * (il[0,1] @ pi @ g3[1,0] - il[1,0] @ pi @ g3[0,1])
# γ
# dΓ γ
# ——— = Y
# dE
self.gammaLE = yL
# γ
# dY 3 ⎛ γ ⎞
# ——— = — i tr ⎜ I Π G³ ⎟
# dE 2 ⎝ ⎠
self.yLE = 1.5j * (il @ pi @ g3).tr()
# Count calls to RG equations.
self.total_iterations += 1
def check_symmetry(self):
"""
Check if all symmetries are fulfilled
"""
self.gamma.check_symmetry()
self.gammaL.check_symmetry()
self.deltaGamma.check_symmetry()
self.deltaGammaL.check_symmetry()
self.z.check_symmetry()
self.yL.check_symmetry()
self.g2.check_symmetry()
self.g3.check_symmetry()
self.current.check_symmetry()
def unpackFlattenedValues(self, flattened_values):
"""
Translate between 1d array used by the solver and Floquet matrices
used in RG equations. Given a 1d array, write the values of this array
to the Floquet matrices self.<values>.
Order of flattened_values:
Γ, Z, δΓ, *G2, *G3, *IL, δΓL, ΓL, YL
"""
if self.compact == 0:
s = self.yL.values.size
m = self.deltaGamma.values.size
l = self.z.values.size
assert flattened_values.size == 10*l+m+7*s
self.gamma.values = flattened_values[:l].reshape(self.gamma.values.shape)
self.z.values = flattened_values[l:2*l].reshape(self.z.values.shape)
self.deltaGamma.values = flattened_values[2*l:2*l+m].reshape(self.deltaGamma.values.shape)
for (g2i, flat) in zip(self.g2.data.flat, np.split(flattened_values[2*l+m:6*l+m], 4)):
g2i.values = flat.reshape(g2i.values.shape)
for (g3i, flat) in zip(self.g3.data.flat, np.split(flattened_values[6*l+m:10*l+m], 4)):
g3i.values = flat.reshape(g3i.values.shape)
for (ii, flat) in zip(self.current.data.flat, np.split(flattened_values[10*l+m:10*l+m+4*s], 4)):
ii.values = flat.reshape(ii.values.shape)
self.deltaGammaL.values = flattened_values[10*l+m+4*s:10*l+m+5*s].reshape(self.deltaGammaL.values.shape)
self.gammaL.values = flattened_values[10*l+m+5*s:10*l+m+6*s].reshape(self.gammaL.values.shape)
self.yL.values = flattened_values[10*l+m+6*s:10*l+m+7*s].reshape(self.yL.values.shape)
elif self.compact == 1:
nmax = self.nmax
f = (2*nmax+1)**2
o = (nmax+1)**2
i = nmax**2
m = nmax*(nmax+1)
shape_f = (2*nmax+1, 2*nmax+1)
shape00 = (nmax+1, nmax+1)
shape11 = (nmax, nmax)
shape01 = (nmax+1, nmax)
shape10 = (nmax, nmax+1)
assert flattened_values.size == 7*(o+i) + 10*m + 6*f
self.gamma.submatrix00 = flattened_values[:o].reshape(shape00)
self.z.submatrix00 = flattened_values[o:2*o].reshape(shape00)
self.deltaGammaL.submatrix00 = flattened_values[2*o:3*o].reshape(shape00)
self.g2[0,0].submatrix00 = flattened_values[3*o:4*o].reshape(shape00)
self.g2[1,1].submatrix00 = flattened_values[4*o:5*o].reshape(shape00)
self.g3[0,0].submatrix00 = flattened_values[5*o:6*o].reshape(shape00)
self.g3[1,1].submatrix00 = flattened_values[6*o:7*o].reshape(shape00)
self.gamma.submatrix11 = flattened_values[7*o:7*o+i].reshape(shape11)
self.z.submatrix11 = flattened_values[7*o+i:7*o+2*i].reshape(shape11)
self.deltaGammaL.submatrix11 = flattened_values[7*o+2*i:7*o+3*i].reshape(shape11)
self.g2[0,0].submatrix11 = flattened_values[7*o+3*i:7*o+4*i].reshape(shape11)
self.g2[1,1].submatrix11 = flattened_values[7*o+4*i:7*o+5*i].reshape(shape11)
self.g3[0,0].submatrix11 = flattened_values[7*o+5*i:7*o+6*i].reshape(shape11)
self.g3[1,1].submatrix11 = flattened_values[7*o+6*i:7*(o+i)].reshape(shape11)
self.gammaL.submatrix01 = flattened_values[7*(o+i):7*(o+i)+m].reshape(shape01)
self.gammaL.submatrix10 = flattened_values[7*(o+i)+m:7*(o+i)+2*m].reshape(shape10)
self.deltaGamma.submatrix01 = flattened_values[7*(o+i)+2*m:7*(o+i)+3*m].reshape(shape01)
self.deltaGamma.submatrix10 = flattened_values[7*(o+i)+3*m:7*(o+i)+4*m].reshape(shape10)
self.yL.submatrix01 = flattened_values[7*(o+i)+4*m:7*(o+i)+5*m].reshape(shape01)
self.yL.submatrix10 = flattened_values[7*(o+i)+5*m:7*(o+i)+6*m].reshape(shape10)
self.current[0,0].submatrix01 = flattened_values[7*(o+i)+6*m:7*(o+i+m)].reshape(shape01)
self.current[0,0].submatrix10 = flattened_values[7*(o+i+m):7*(o+i)+8*m].reshape(shape10)
self.current[1,1].submatrix01 = flattened_values[7*(o+i)+8*m:7*(o+i)+9*m].reshape(shape01)
self.current[1,1].submatrix10 = flattened_values[7*(o+i)+9*m:7*(o+i)+10*m].reshape(shape10)
self.g2[0,1].values = flattened_values[7*(o+i)+10*m:7*(o+i)+10*m+f].reshape(shape_f)
self.g2[1,0].values = flattened_values[7*(o+i)+10*m+f:7*(o+i)+10*m+2*f].reshape(shape_f)
self.g3[0,1].values = flattened_values[7*(o+i)+10*m+2*f:7*(o+i)+10*m+3*f].reshape(shape_f)
self.g3[1,0].values = flattened_values[7*(o+i)+10*m+3*f:7*(o+i)+10*m+4*f].reshape(shape_f)
self.current[0,1].values = flattened_values[7*(o+i)+10*m+4*f:7*(o+i)+10*m+5*f].reshape(shape_f)
self.current[1,0].values = flattened_values[7*(o+i)+10*m+5*f:].reshape(shape_f)
#self.g2[0,1].setValues(flattened_values[7*(o+i)+10*m:7*(o+i)+10*m+f].reshape(shape_f), True, True)
#self.g2[1,0].setValues(flattened_values[7*(o+i)+10*m+f:7*(o+i)+10*m+2*f].reshape(shape_f), True, True)
#self.g3[0,1].setValues(flattened_values[7*(o+i)+10*m+2*f:7*(o+i)+10*m+3*f].reshape(shape_f), True, True)
#self.g3[1,0].setValues(flattened_values[7*(o+i)+10*m+3*f:7*(o+i)+10*m+4*f].reshape(shape_f), True, True)
#self.current[0,1].setValues(flattened_values[7*(o+i)+10*m+4*f:7*(o+i)+10*m+5*f].reshape(shape_f), True, True)
#self.current[1,0].setValues(flattened_values[7*(o+i)+10*m+5*f:].reshape(shape_f), True, True)
elif self.compact == 2:
nmax = self.nmax
f = (2*nmax+1)**2
o = ((nmax+1)**2 + 1) // 2
i = (nmax**2 + 1) // 2
m = nmax*(nmax+1) // 2
assert flattened_values.size == 5*(o+i) + 8*m + 3*f
unpack01 = lambda flat, sym: np.concatenate((flat, sym*flat[::-1].conjugate()), axis=None).reshape((nmax+1,nmax))
unpack10 = lambda flat, sym: np.concatenate((flat, sym*flat[::-1].conjugate()), axis=None).reshape((nmax,nmax+1))
if nmax % 2:
unpack00 = lambda flat, sym: np.concatenate((flat, sym*flat[::-1].conjugate()), axis=None).reshape((nmax+1,nmax+1))
unpack11 = lambda flat, sym: np.concatenate((flat, sym*flat[-2::-1].conjugate()), axis=None).reshape((nmax,nmax))
else:
unpack00 = lambda flat, sym: np.concatenate((flat, sym*flat[-2::-1].conjugate()), axis=None).reshape((nmax+1,nmax+1))
unpack11 = lambda flat, sym: np.concatenate((flat, sym*flat[::-1].conjugate()), axis=None).reshape((nmax,nmax))
self.gamma.submatrix00 = unpack00(flattened_values[:o], 1)
self.gamma.submatrix11 = unpack11(flattened_values[5*o:5*o+i], 1)
self.z.submatrix00 = unpack00(flattened_values[o:2*o], 1)
self.z.submatrix11 = unpack11(flattened_values[5*o+i:5*o+2*i], 1)
self.deltaGammaL.submatrix00 = unpack00(flattened_values[2*o:3*o], 1)
self.deltaGammaL.submatrix11 = unpack11(flattened_values[5*o+2*i:5*o+3*i], 1)
self.g2[0,0].submatrix00 = unpack00(flattened_values[3*o:4*o], 1)
self.g2[0,0].submatrix11 = unpack11(flattened_values[5*o+3*i:5*o+4*i], 1)
self.g2[1,1] = self.g2[0,0].copy()
self.g3[0,0].submatrix00 = unpack00(flattened_values[4*o:5*o], -1)
self.g3[0,0].submatrix11 = unpack11(flattened_values[5*o+4*i:5*(o+i)], -1)
self.g3[1,1] = self.g3[0,0].copy()
self.gammaL.submatrix01 = unpack01(flattened_values[5*(o+i):5*(o+i)+m], 1)
self.gammaL.submatrix10 = unpack10(flattened_values[5*(o+i)+m:5*(o+i)+2*m], 1)
self.deltaGamma.submatrix01 = unpack01(flattened_values[5*(o+i)+2*m:5*(o+i)+3*m], 1)
self.deltaGamma.submatrix10 = unpack10(flattened_values[5*(o+i)+3*m:5*(o+i)+4*m], 1)
self.yL.submatrix01 = unpack01(flattened_values[5*(o+i)+4*m:5*(o+i+m)], -1)
self.yL.submatrix10 = unpack10(flattened_values[5*(o+i+m):5*(o+i)+6*m], -1)
self.current[0,0].submatrix01 = unpack01(flattened_values[5*(o+i)+6*m:5*(o+i)+7*m], -1)
self.current[0,0].submatrix10 = unpack10(flattened_values[5*(o+i)+7*m:5*(o+i)+8*m], -1)
self.current[1,1] = self.current[0,0].copy()
#self.g2[0,1].setValues(flattened_values[5*(o+i)+8*m:5*(o+i)+8*m+f].reshape((2*nmax+1, 2*nmax+1)), True, True)
#self.g3[0,1].setValues(flattened_values[5*(o+i)+8*m+f:5*(o+i)+8*m+2*f].reshape((2*nmax+1, 2*nmax+1)), True, True)
#self.current[0,1].setValues(flattened_values[5*(o+i)+8*m+2*f:].reshape((2*nmax+1, 2*nmax+1)), True, True)
self.g2[0,1].values = flattened_values[5*(o+i)+8*m:5*(o+i)+8*m+f].reshape((2*nmax+1, 2*nmax+1))
self.g2[1,0] = self.g2[0,1].floquetConjugate()
self.g3[0,1].values = flattened_values[5*(o+i)+8*m+f:5*(o+i)+8*m+2*f].reshape((2*nmax+1, 2*nmax+1))
self.g3[1,0] = -self.g3[0,1].floquetConjugate()
self.current[0,1].values = flattened_values[5*(o+i)+8*m+2*f:].reshape((2*nmax+1, 2*nmax+1))
self.current[1,0] = -self.current[0,1].floquetConjugate()
if settings.CHECK_SYMMETRIES:
self.check_symmetry()
def packFlattenedDerivatives(self):
"""
Pack Floquet matrices representing derivatives in one flattened
(1d) array for the solver.
Order of flattened_values:
Γ, Z, δΓ, *G2, *G3, *IL, δΓL, ΓL, YL
"""
if self.compact == 0:
return np.concatenate((
self.gammaE.values,
self.zE.values,
self.deltaGammaE.values,
*(g.values for g in self.g2E.data.flat),
*(g.values for g in self.g3E.data.flat),
*(i.values for i in self.currentE.data.flat),
self.deltaGammaLE.values,
self.gammaLE.values,
self.yLE.values,
), axis = None
)
elif self.compact == 1:
if settings.CHECK_SYMMETRIES:
assert self.gammaE.submatrix01 is None
assert self.gammaE.submatrix10 is None
assert self.zE.submatrix01 is None
assert self.zE.submatrix10 is None
assert self.deltaGammaLE.submatrix01 is None
assert self.deltaGammaLE.submatrix10 is None
assert self.gammaLE.submatrix00 is None
assert self.gammaLE.submatrix11 is None
assert self.deltaGammaE.submatrix00 is None
assert self.deltaGammaE.submatrix11 is None
assert self.yLE.submatrix00 is None
assert self.yLE.submatrix11 is None
assert self.g2E[0,0].submatrix01 is None
assert self.g2E[1,1].submatrix01 is None
assert self.g3E[0,0].submatrix01 is None
assert self.g3E[1,1].submatrix01 is None
assert self.g2E[0,0].submatrix10 is None
assert self.g2E[1,1].submatrix10 is None
assert self.g3E[0,0].submatrix10 is None
assert self.g3E[1,1].submatrix10 is None
assert self.currentE[0,0].submatrix00 is None
assert self.currentE[0,0].submatrix11 is None
assert self.currentE[1,1].submatrix00 is None
assert self.currentE[1,1].submatrix11 is None
if self.yLE.submatrix01 is None:
dummy = np.zeros(self.nmax*(self.nmax+1), dtype=np.complex128)
self.yLE.submatrix01 = dummy
self.yLE.submatrix10 = dummy
if self.deltaGammaE.submatrix01 is None:
self.deltaGammaE.submatrix01 = dummy
self.deltaGammaE.submatrix10 = dummy
if self.currentE[0,0].submatrix01 is None:
self.currentE[0,0].submatrix01 = dummy
self.currentE[0,0].submatrix10 = dummy
if self.currentE[1,1].submatrix01 is None:
self.currentE[1,1].submatrix01 = dummy
self.currentE[1,1].submatrix10 = dummy
return np.concatenate((
self.gammaE.submatrix00,
self.zE.submatrix00,
self.deltaGammaLE.submatrix00,
self.g2E[0,0].submatrix00,
self.g2E[1,1].submatrix00,
self.g3E[0,0].submatrix00,
self.g3E[1,1].submatrix00,
self.gammaE.submatrix11,
self.zE.submatrix11,
self.deltaGammaLE.submatrix11,
self.g2E[0,0].submatrix11,
self.g2E[1,1].submatrix11,
self.g3E[0,0].submatrix11,
self.g3E[1,1].submatrix11,
self.gammaLE.submatrix01,
self.gammaLE.submatrix10,
self.deltaGammaE.submatrix01,
self.deltaGammaE.submatrix10,
self.yLE.submatrix01,
self.yLE.submatrix10,
self.currentE[0,0].submatrix01,
self.currentE[0,0].submatrix10,
self.currentE[1,1].submatrix01,
self.currentE[1,1].submatrix10,
self.g2E[0,1].values,
self.g2E[1,0].values,
self.g3E[0,1].values,
self.g3E[1,0].values,
self.currentE[0,1].values,
self.currentE[1,0].values,
), axis = None
)
elif self.compact == 2:
if settings.CHECK_SYMMETRIES:
assert self.gammaE.symmetry == -1
assert self.gammaE.submatrix01 is None
assert self.gammaE.submatrix10 is None
assert self.zE.symmetry == -1
assert self.zE.submatrix01 is None
assert self.zE.submatrix10 is None
assert self.deltaGammaLE.symmetry == -1
assert self.deltaGammaLE.submatrix01 is None
assert self.deltaGammaLE.submatrix10 is None
assert self.gammaLE.symmetry == -1
assert self.gammaLE.submatrix00 is None
assert self.gammaLE.submatrix11 is None
assert self.deltaGammaE.symmetry == -1
assert self.deltaGammaE.submatrix00 is None
assert self.deltaGammaE.submatrix11 is None
assert self.yLE.symmetry == 1
assert self.yLE.submatrix00 is None
assert self.yLE.submatrix11 is None
assert self.g2E.symmetry == -1
assert self.g3E.symmetry == 1
assert self.currentE.symmetry == 1
assert self.g2E[0,0].submatrix01 is None
assert self.g3E[0,0].submatrix01 is None
assert self.g2E[0,0].submatrix10 is None
assert self.g3E[0,0].submatrix10 is None
assert self.currentE[0,0].submatrix00 is None
assert self.currentE[0,0].submatrix11 is None
if self.yLE.submatrix01 is None:
dummy = np.zeros(self.nmax*(self.nmax+1), dtype=np.complex128)
self.yLE.submatrix01 = dummy
self.yLE.submatrix10 = dummy
if self.deltaGammaE.submatrix01 is None:
self.deltaGammaE.submatrix01 = dummy
self.deltaGammaE.submatrix10 = dummy
if self.currentE[0,0].submatrix01 is None:
self.currentE[0,0].submatrix01 = dummy
self.currentE[0,0].submatrix10 = dummy
pack_flat = lambda x, sym: (x.reshape(-1)[:(x.size+1)//2] + sym*x.reshape(-1)[:x.size//2-1:-1].conjugate())/2
return np.concatenate((
pack_flat(self.gammaE.submatrix00, -1),
pack_flat(self.zE.submatrix00, -1),
pack_flat(self.deltaGammaLE.submatrix00, -1),
pack_flat(self.g2E[0,0].submatrix00, -1),
pack_flat(self.g3E[0,0].submatrix00, 1),
pack_flat(self.gammaE.submatrix11, -1),
pack_flat(self.zE.submatrix11, -1),
pack_flat(self.deltaGammaLE.submatrix11, -1),
pack_flat(self.g2E[0,0].submatrix11, -1),
pack_flat(self.g3E[0,0].submatrix11, 1),
pack_flat(self.gammaLE.submatrix01, -1),
pack_flat(self.gammaLE.submatrix10, -1),
pack_flat(self.deltaGammaE.submatrix01, -1),
pack_flat(self.deltaGammaE.submatrix10, -1),
pack_flat(self.yLE.submatrix01, 1),
pack_flat(self.yLE.submatrix10, 1),
pack_flat(self.currentE[0,0].submatrix01, 1),
pack_flat(self.currentE[0,0].submatrix10, 1),
self.g2E[0,1].values,
self.g3E[0,1].values,
self.currentE[0,1].values,
), axis=None)
def packFlattenedValues(self):
"""
Translate between 1d array used by the solver and Floquet matrices
used in RG equations. Collect all Floquet matrices in one flattened
(1d) array.
Order of flattened_values:
Γ, Z, δΓ, *G2, *G3, *IL, δΓL, ΓL, YL
"""
if self.compact == 0:
return np.concatenate((
self.gamma.values,
self.z.values,
self.deltaGamma.values,
*(g.values for g in self.g2.data.flat),
*(g.values for g in self.g3.data.flat),
*(i.values for i in self.current.data.flat),
self.deltaGammaL.values,
self.gammaL.values,
self.yL.values,
), axis = None
)
elif self.compact == 1:
if settings.CHECK_SYMMETRIES:
assert self.gamma.submatrix01 is None
assert self.gamma.submatrix10 is None
assert self.z.submatrix01 is None
assert self.z.submatrix10 is None
assert self.deltaGammaL.submatrix01 is None
assert self.deltaGammaL.submatrix10 is None
assert self.gammaL.submatrix00 is None
assert self.gammaL.submatrix11 is None
assert self.deltaGamma.submatrix00 is None
assert self.deltaGamma.submatrix11 is None
assert self.yL.submatrix00 is None
assert self.yL.submatrix11 is None
assert self.g2[0,0].submatrix01 is None
assert self.g2[1,1].submatrix01 is None
assert self.g3[0,0].submatrix01 is None
assert self.g3[1,1].submatrix01 is None
assert self.g2[0,0].submatrix10 is None
assert self.g2[1,1].submatrix10 is None
assert self.g3[0,0].submatrix10 is None
assert self.g3[1,1].submatrix10 is None
assert self.current[0,0].submatrix00 is None
assert self.current[0,0].submatrix11 is None
assert self.current[1,1].submatrix00 is None
assert self.current[1,1].submatrix11 is None
return np.concatenate((
self.gamma.submatrix00,
self.z.submatrix00,
self.deltaGammaL.submatrix00,
self.g2[0,0].submatrix00,
self.g2[1,1].submatrix00,
self.g3[0,0].submatrix00,
self.g3[1,1].submatrix00,
self.gamma.submatrix11,
self.z.submatrix11,
self.deltaGammaL.submatrix11,
self.g2[0,0].submatrix11,
self.g2[1,1].submatrix11,
self.g3[0,0].submatrix11,
self.g3[1,1].submatrix11,
self.gammaL.submatrix01,
self.gammaL.submatrix10,
self.deltaGamma.submatrix01,
self.deltaGamma.submatrix10,
self.yL.submatrix01,
self.yL.submatrix10,
self.current[0,0].submatrix01,
self.current[0,0].submatrix10,
self.current[1,1].submatrix01,
self.current[1,1].submatrix10,
self.g2[0,1].values,
self.g2[1,0].values,
self.g3[0,1].values,
self.g3[1,0].values,
self.current[0,1].values,
self.current[1,0].values,
), axis = None
)
elif self.compact == 2:
if settings.CHECK_SYMMETRIES:
assert self.gamma.symmetry == 1
assert self.gamma.submatrix01 is None
assert self.gamma.submatrix10 is None
assert self.z.symmetry == 1
assert self.z.submatrix01 is None
assert self.z.submatrix10 is None
assert self.deltaGammaL.symmetry == 1
assert self.deltaGammaL.submatrix01 is None
assert self.deltaGammaL.submatrix10 is None
assert self.gammaL.symmetry == 1
assert self.gammaL.submatrix00 is None
assert self.gammaL.submatrix11 is None
assert self.deltaGamma.symmetry == 1
assert self.deltaGamma.submatrix00 is None
assert self.deltaGamma.submatrix11 is None
assert self.yL.symmetry == -1
assert self.yL.submatrix00 is None
assert self.yL.submatrix11 is None
assert self.g2.symmetry == 1
assert self.g3.symmetry == -1
assert self.current.symmetry == -1
assert self.g2[0,0].submatrix01 is None
assert self.g3[0,0].submatrix01 is None
assert self.g2[0,0].submatrix10 is None
assert self.g3[0,0].submatrix10 is None
assert self.current[0,0].submatrix00 is None
assert self.current[0,0].submatrix11 is None
pack_flat = lambda x, sym: (x.reshape(-1)[:(x.size+1)//2] + sym*x.reshape(-1)[:x.size//2-1:-1].conjugate())/2
return np.concatenate((
pack_flat(self.gamma.submatrix00, 1),
pack_flat(self.z.submatrix00, 1),
pack_flat(self.deltaGammaL.submatrix00, 1),
pack_flat(self.g2[0,0].submatrix00, 1),
pack_flat(self.g3[0,0].submatrix00, -1),
pack_flat(self.gamma.submatrix11, 1),
pack_flat(self.z.submatrix11, 1),
pack_flat(self.deltaGammaL.submatrix11, 1),
pack_flat(self.g2[0,0].submatrix11, 1),
pack_flat(self.g3[0,0].submatrix11, -1),
pack_flat(self.gammaL.submatrix01, 1),
pack_flat(self.gammaL.submatrix10, 1),
pack_flat(self.deltaGamma.submatrix01, 1),
pack_flat(self.deltaGamma.submatrix10, 1),
pack_flat(self.yL.submatrix01, -1),
pack_flat(self.yL.submatrix10, -1),
pack_flat(self.current[0,0].submatrix01, -1),
pack_flat(self.current[0,0].submatrix10, -1),
self.g2[0,1].values,
self.g3[0,1].values,
self.current[0,1].values,
), axis=None)
def hash(self):
data = self.packFlattenedValues()
data.flags["WRITEABLE"] = False
return hashlib.sha1(data.data).hexdigest()
def odeFunctionIm(self, imenergy, flattened_values):
"""
ODE as given to the solver for solving the RG equations along the
imaginary axis. Given a flattened array containing all Floquet
matrices, evaluate the RG equations and return a flattened array of
all derivatives.
imenergy = Im(E) is used as function argument since the solver cannot
handle complex flow parameters.
"""
try:
if self.save_filename and self.save_iterations > 0 and self.iterations % self.save_iterations == 0:
try:
self.save_compact()
except:
settings.logger.exception("Failed to save intermediate result:")
except AttributeError:
pass
except:
settings.logger.exception("Failed trying to save intermediate result:")
self.global_properties.energy = self.energy.real + 1j*imenergy
self.unpackFlattenedValues(flattened_values)
# Evaluate RG equations
try:
self.updateRGequations()
except KeyboardInterrupt as error:
settings.logger.critical('Interrupted at Im(E) = %e'%imenergy)
try:
self.ir_cutoff = max(self.ir_cutoff, imenergy)
except AttributeError:
self.ir_cutoff = imenergy
raise error
except Exception as error:
settings.logger.exception('Unhandled error at Im(E) = %e'%imenergy)
try:
self.ir_cutoff = max(self.ir_cutoff, imenergy)
except AttributeError:
self.ir_cutoff = imenergy
raise error
# Pack values
# use d/d(Im E) = i d/dE
return 1j*self.packFlattenedDerivatives()
def odeFunctionRe(self, reenergy, flattened_values):
"""
ODE as given to the solver for solving the RG equations along the
real axis. Given a flattened array containing all Floquet matrices,
evaluate the RG equations and return a flattened array of all
derivatives.
"""
if self.save_filename and self.save_iterations > 0 and self.iterations % self.save_iterations == 0:
try:
self.save_compact()
except:
settings.logger.exception('Failed to save intermediate result:')
self.global_properties.energy = reenergy + 1j*self.energy.imag
self.unpackFlattenedValues(flattened_values)
# Evaluate RG equations
try:
self.updateRGequations()
except KeyboardInterrupt as error:
settings.logger.critical('Interrupted at Re(E) = %g, Im(E) = %g'%(reenergy, energy0.imag))
raise error
except:
settings.logger.exception('Unhandled error at Re(E) = %g, Im(E) = %g'%(reenergy, energy0.imag))
raise error
# Pack values
return self.packFlattenedDerivatives()
def solveOdeIm(self, eiminit, eimfinal, init_values=None, g2max=1e6, only_final=False, **solveopts):
"""
Solve the RG equations along the imaginary axis, starting from
E = Ereal + eiminit*1j and ending at E = Ereal + eimfinal*1j
where Ereal is the real part of the current energy.
Other arguments:
init_values : flattened array of initial values, by default taken from
self.packFlattenedValues()
g2max : Threshold for a trigger element in one of the Floquet matrices.
This element will detect whether a pole is reached by the RG
flow and will stop the flow at the pole.
This is only implemented for self.compact == 0
only_final : only save final result, do not save the RG flow.
This saves memory.
**solveopts : arguments that are directly passed on to the solver. Most
relevant are rtol and atol.
"""
assert np.allclose(self.energy.imag, eiminit)
if init_values is None:
init_values = self.packFlattenedValues()
output = solve_ivp(
self.odeFunctionIm,
(eiminit, eimfinal),
init_values,
t_eval = ((eimfinal,) if only_final else None),
**solveopts
)
return output
def solveOdeRe(self, reEinit, reEfinal, init_values=None, g2max=1e6, only_final=False, **solveopts):
"""
Solve the RG equations along the real axis, starting from
E = reEinit + 1j*Eimag and ending at E = reEfinal + 1j*Eimag.
where Eimag is the imaginary part of the current energy.
Other arguments:
init_values : flattened array of initial values, by default taken from
self.packFlattenedValues()
g2max : Threshold for a trigger element in one of the Floquet matrices.
This element will detect whether a pole is reached by the RG
flow and will stop the flow at the pole.
This is only implemented for self.compact == 0
only_final : only save final result, do not save the RG flow.
This saves memory.
**solveopts : arguments that are directly passed on to the solver. Most
relevant are rtol and atol.
"""
assert abs(self.energy.real - reEinit) < 1e-8
if init_values is None:
init_values = self.packFlattenedValues()
output = solve_ivp(
self.odeFunctionRe,
(reEinit, reEfinal),
init_values,
t_eval = ((eimfinal,) if only_final else None),
**solveopts
)
return output
def save_compact(self, values, compressed=False):
"""
Automatically save current state of the RG flow in the most compact
form. The file name will be
self.save_filename % self.iterations
or
self.save_filename.format(self.iterations).
In a computationally expensive RG flow this allows saving intermediate
steps of the RG flow.
"""
try:
filename = self.save_filename%self.iterations
except TypeError:
filename = self.save_filename.format(self.iterations)
(np.savez_compressed if compressed else np.savez)(
filename,
values = values,
energy = self.energy,
compact = self.compact,
)
def load_compact(self, filename):
"""
Load a file that was created with Kondo.save_compact. This overwrites
the current state of self with the values given in the file.
"""
data = np.load(filename)
assert data['compact'] == self.compact
self.unpackFlattenedValues(data['values'])
self.global_properties.energy = data['energy']