From 011d0ed828756e044c765a372163bc8ceb6d9eae Mon Sep 17 00:00:00 2001
From: Valentin Bruch <valentin.bruch@rwth-aachen.de>
Date: Fri, 22 Jul 2022 14:13:17 +0200
Subject: [PATCH] updated package
---
README.md | 7 +-
compact_rtrg.py | 3 -
kondo.py | 2 +
package/README.md | 5 +
package/pyproject.toml | 2 +-
package/setup.py | 3 +
package/src/frtrg_kondo/compact_rtrg.py | 126 ++--
package/src/frtrg_kondo/data_management.py | 222 +++++--
package/src/frtrg_kondo/gen_data.py | 53 +-
package/src/frtrg_kondo/kondo.py | 706 ++++++++++++++-------
package/src/frtrg_kondo/plot.py | 529 ++++++++++++---
package/src/frtrg_kondo/reservoirmatrix.py | 220 ++++---
package/src/frtrg_kondo/rtrg.py | 108 ++--
package/src/frtrg_kondo/rtrg_c.c | 260 +++++++-
package/src/frtrg_kondo/settings.py | 8 +-
plot.py | 2 +-
setup.py | 93 +--
17 files changed, 1641 insertions(+), 708 deletions(-)
diff --git a/README.md b/README.md
index 1ed682b..9537143 100644
--- a/README.md
+++ b/README.md
@@ -48,14 +48,13 @@ In `kondo.py` the framework of FRTRG is applied to the Kondo model.
## Versions
-baseversion 14:
+version 14.x:
* based on versions 8 and 13
-* using the same unitary transformation of the initial Hamiltonian as in version 8, but the improvements in the code from version 13
+* implements the RG for the system with and without unitary transformation (version 8 and 13)
* completely new user interface
baseversion 13:
-* best tested version so far
-* used to generate data for upcoming publication
+* used to generate data for <https://arxiv.org/abs/2206.06263>
baseversion 12:
* handling of frequency shifts in initial conditions for Floquet matrices is a bit strange
diff --git a/compact_rtrg.py b/compact_rtrg.py
index 0882a40..d79472c 100644
--- a/compact_rtrg.py
+++ b/compact_rtrg.py
@@ -22,7 +22,6 @@ class SymRGfunction(RGfunction):
Subclass of RGfunction for Floquet matrices representing functions which in
time-domain fulfill f(t + T/2) = ±f(t).
"""
-
def __init__(self, global_properties, values, symmetry=0, diag=None, offdiag=None, **kwargs):
self.global_properties = global_properties
self.symmetry = symmetry
@@ -84,8 +83,6 @@ class SymRGfunction(RGfunction):
self.submatrix01 = values[0::2,1::2] # (n+1)x n matrix for -
self.submatrix10 = values[1::2,0::2] # n x(n+1) matrix for -
-
-
@values.setter
def values(self, values):
values = np.asarray(values, np.complex128)
diff --git a/kondo.py b/kondo.py
index 5fe83d5..d6e31fe 100644
--- a/kondo.py
+++ b/kondo.py
@@ -63,6 +63,7 @@ from compact_rtrg import SymRGfunction
REF_TIME = time()
LAST_LOG_TIME = REF_TIME
+
def driving_voltage(tau, *fourier_coef):
"""
Generate function of time given Fourier coefficients.
@@ -77,6 +78,7 @@ def driving_voltage(tau, *fourier_coef):
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.
diff --git a/package/README.md b/package/README.md
index ddd54d5..bc8469a 100644
--- a/package/README.md
+++ b/package/README.md
@@ -18,3 +18,8 @@ python -m frtrg_kondo.gen_data --db_filename=/path/to/db.sqlite --filename=/path
# Generate 21 data points with vdc=0,1,...,20 using 4 parallel processes
OMP_NUM_TREADS=1 python -m frtrg_kondo.gen_data --db_filename=/path/to/db.sqlite --filename=/path/to/HDF5-file.h5 --nmax=14 --voltage_branches=3 --method=mu --vdc 0 20 --steps=21 --vac=5 --omega=10 --threads=4
```
+
+Build with:
+```sh
+python -m build -w
+```
diff --git a/package/pyproject.toml b/package/pyproject.toml
index 843ed73..fc9497e 100644
--- a/package/pyproject.toml
+++ b/package/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
[project]
name = "frtrg_kondo"
-version = "0.14.6"
+version = "0.14.8"
authors = [
{ name="Valentin Bruch", email="valentin.bruch@rwth-aachen.de" },
]
diff --git a/package/setup.py b/package/setup.py
index b707a46..5ca33fe 100644
--- a/package/setup.py
+++ b/package/setup.py
@@ -35,6 +35,9 @@ def main():
if 'DEBUG' in environ:
compiler_args += ['-DDEBUG']
+ if 'ANALYZE' in environ:
+ compiler_args += ['-DANALYZE']
+
module = Extension(
"frtrg_kondo.rtrg_c",
sources = ['src/frtrg_kondo/rtrg_c.c'],
diff --git a/package/src/frtrg_kondo/compact_rtrg.py b/package/src/frtrg_kondo/compact_rtrg.py
index a355832..4243e66 100644
--- a/package/src/frtrg_kondo/compact_rtrg.py
+++ b/package/src/frtrg_kondo/compact_rtrg.py
@@ -18,23 +18,11 @@ OVERWRITE_BOTH = bytes((3,))
class SymRGfunction(RGfunction):
- '''
+ """
Subclass of RGfunction for Floquet matrices representing functions which in
time-domain fulfill f(t + T/2) = ±f(t).
- '''
- # Flags:
- # 1 << 0 : matrix C contiguous
- # 1 << 1 : matrix F contiguous
- INVC_COUNTER = [0 for i in range(1 << 2)]
- # Flags:
- # 1 << 0 = 0x01 : symmetric
- # 1 << 1 = 0x02 : matrix 1 C contiguous
- # 1 << 2 = 0x04 : matrix 1 F contiguous
- # 1 << 3 = 0x08 : matrix 2 C contiguous
- # 1 << 4 = 0x10 : matrix 2 F contiguous
- MMC_COUNTER = [0 for i in range(1 << 5)]
-
- def __init__(self, global_properties, values, symmetry=0, **kwargs):
+ """
+ def __init__(self, global_properties, values, symmetry=0, diag=None, offdiag=None, **kwargs):
self.global_properties = global_properties
self.symmetry = symmetry
self.energy_shifted_copies = {}
@@ -54,15 +42,20 @@ class SymRGfunction(RGfunction):
self.submatrix00 = np.identity(self.nmax+1, dtype=np.complex128)
self.submatrix11 = np.identity(self.nmax, dtype=np.complex128)
elif values is not None:
- self.values = values
+ if diag is None or offdiag is None or settings.CHECK_SYMMETRIES:
+ self.values = values
+ else:
+ self.setValues(values, diag, offdiag)
@classmethod
- def fromRGfunction(cls, rgfunc):
+ def fromRGfunction(cls, rgfunc, diag=None, offdiag=None):
assert isinstance(rgfunc, RGfunction)
return SymRGfunction(
rgfunc.global_properties,
rgfunc.values,
rgfunc.symmetry,
+ diag,
+ offdiag,
)
@property
@@ -78,6 +71,18 @@ class SymRGfunction(RGfunction):
values[1::2,1::2] = self.submatrix11
return values
+ def setValues(self, values, diag:bool, offdiag:bool):
+ """
+ More efficient than setting values by self.values = values if
+ it is known which parts are non-zero.
+ """
+ if diag:
+ self.submatrix00 = values[0::2,0::2] # (n+1)x(n+1) matrix for +
+ self.submatrix11 = values[1::2,1::2] # n x n matrix for +
+ if offdiag:
+ self.submatrix01 = values[0::2,1::2] # (n+1)x n matrix for -
+ self.submatrix10 = values[1::2,0::2] # n x(n+1) matrix for -
+
@values.setter
def values(self, values):
values = np.asarray(values, np.complex128)
@@ -87,21 +92,21 @@ class SymRGfunction(RGfunction):
self.submatrix01 = values[0::2,1::2] # (n+1)x n matrix for -
self.submatrix10 = values[1::2,0::2] # n x(n+1) matrix for -
self.submatrix11 = values[1::2,1::2] # n x n matrix for +
- if np.abs(self.submatrix00).max() < 1e-15:
+ if np.all(np.abs(self.submatrix00) < 1e-15):
self.submatrix00 = None
- if np.abs(self.submatrix01).max() < 1e-15:
+ if np.all(np.abs(self.submatrix01) < 1e-15):
self.submatrix01 = None
- if np.abs(self.submatrix10).max() < 1e-15:
+ if np.all(np.abs(self.submatrix10) < 1e-15):
self.submatrix10 = None
- if np.abs(self.submatrix11).max() < 1e-15:
+ if np.all(np.abs(self.submatrix11) < 1e-15):
self.submatrix11 = None
assert (self.submatrix00 is None and self.submatrix11 is None) or (self.submatrix01 is None and self.submatrix10 is None)
self.energy_shifted_copies.clear()
def copy(self):
- '''
+ """
Copy only values, take a reference to global_properties.
- '''
+ """
return SymRGfunction(
self.global_properties,
values = None,
@@ -113,9 +118,9 @@ class SymRGfunction(RGfunction):
)
def reduced(self, shift=0):
- '''
+ """
Remove voltage-shifted copies
- '''
+ """
assert shift == 0
return self
@@ -123,7 +128,7 @@ class SymRGfunction(RGfunction):
return self
def floquetConjugate(self):
- '''
+ """
For a Floquet matrix A(E)_{nm} this returns the C-transform
C A(E)_{nm} C = A(-E*)_{-n,-m}
with the superoperator C defined by
@@ -134,7 +139,7 @@ class SymRGfunction(RGfunction):
This can only be evaluated if the energy of self lies on the
imaginary axis.
- '''
+ """
if self.symmetry == 1:
return self.copy()
elif self.symmetry == -1:
@@ -150,13 +155,13 @@ class SymRGfunction(RGfunction):
)
def __matmul__(self, other):
- '''
+ """
Convolution (or product in Floquet space) of two RG functions.
Other must be of type SymRGfunction.
Note: This is only approximately associative, as long the function
converges to 0 for |n| of order of nmax.
- '''
+ """
if not isinstance(other, SymRGfunction):
if isinstance(other, RGfunction):
return self.toRGfunction() @ other
@@ -171,30 +176,18 @@ class SymRGfunction(RGfunction):
assert (self.submatrix11 is not None and other.submatrix11 is not None)
res00 = rtrg_c.multiply_extended(other.submatrix00.T, self.submatrix00.T, self.padding//2, symmetry, self.clear_corners//2).T
res11 = rtrg_c.multiply_extended(other.submatrix11.T, self.submatrix11.T, self.padding//2, symmetry, self.clear_corners//2).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix00.T.flags.c_contiguous << 3) | (self.submatrix00.T.flags.f_contiguous << 4) | (other.submatrix00.T.flags.c_contiguous << 1) | (other.submatrix00.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix11.T.flags.c_contiguous << 3) | (self.submatrix11.T.flags.f_contiguous << 4) | (other.submatrix11.T.flags.c_contiguous << 1) | (other.submatrix11.T.flags.f_contiguous << 2)] += 1
elif (self.submatrix01 is not None and other.submatrix01 is not None):
assert (self.submatrix10 is not None and other.submatrix10 is not None)
res00 = rtrg_c.multiply_extended(other.submatrix10.T, self.submatrix01.T, self.padding//2, symmetry, self.clear_corners//2).T
res11 = rtrg_c.multiply_extended(other.submatrix01.T, self.submatrix10.T, self.padding//2, symmetry, self.clear_corners//2).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix01.flags.c_contiguous << 3) | (self.submatrix01.flags.f_contiguous << 4) | (other.submatrix10.flags.c_contiguous << 1) | (other.submatrix10.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix10.flags.c_contiguous << 3) | (self.submatrix10.flags.f_contiguous << 4) | (other.submatrix01.flags.c_contiguous << 1) | (other.submatrix01.flags.f_contiguous << 2)] += 1
elif (self.submatrix00 is not None and other.submatrix01 is not None):
assert (self.submatrix11 is not None and other.submatrix10 is not None)
res01 = rtrg_c.multiply_extended(other.submatrix01.T, self.submatrix00.T, self.padding//2, symmetry, self.clear_corners//2).T
res10 = rtrg_c.multiply_extended(other.submatrix10.T, self.submatrix11.T, self.padding//2, symmetry, self.clear_corners//2).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix00.T.flags.c_contiguous << 3) | (self.submatrix00.T.flags.f_contiguous << 4) | (other.submatrix01.T.flags.c_contiguous << 1) | (other.submatrix01.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix11.T.flags.c_contiguous << 3) | (self.submatrix11.T.flags.f_contiguous << 4) | (other.submatrix10.T.flags.c_contiguous << 1) | (other.submatrix10.T.flags.f_contiguous << 2)] += 1
elif (self.submatrix01 is not None and other.submatrix00 is not None):
assert (self.submatrix10 is not None and other.submatrix11 is not None)
res01 = rtrg_c.multiply_extended(other.submatrix11.T, self.submatrix01.T, self.padding//2, symmetry, self.clear_corners//2).T
res10 = rtrg_c.multiply_extended(other.submatrix00.T, self.submatrix10.T, self.padding//2, symmetry, self.clear_corners//2).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix01.T.flags.c_contiguous << 3) | (self.submatrix01.T.flags.f_contiguous << 4) | (other.submatrix11.T.flags.c_contiguous << 1) | (other.submatrix11.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix10.T.flags.c_contiguous << 3) | (self.submatrix10.T.flags.f_contiguous << 4) | (other.submatrix00.T.flags.c_contiguous << 1) | (other.submatrix00.T.flags.f_contiguous << 2)] += 1
return SymRGfunction(
self.global_properties,
values = None,
@@ -222,44 +215,32 @@ class SymRGfunction(RGfunction):
assert (self.submatrix11 is not None and other.submatrix11 is not None)
self.submatrix00 = rtrg_c.multiply_extended(other.submatrix00.T, self.submatrix00.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
self.submatrix11 = rtrg_c.multiply_extended(other.submatrix11.T, self.submatrix11.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix00.T.flags.c_contiguous << 3) | (self.submatrix00.T.flags.f_contiguous << 4) | (other.submatrix00.T.flags.c_contiguous << 1) | (other.submatrix00.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix11.T.flags.c_contiguous << 3) | (self.submatrix11.T.flags.f_contiguous << 4) | (other.submatrix11.T.flags.c_contiguous << 1) | (other.submatrix11.T.flags.f_contiguous << 2)] += 1
elif (self.submatrix01 is not None and other.submatrix01 is not None):
assert (self.submatrix10 is not None and other.submatrix10 is not None)
self.submatrix00 = rtrg_c.multiply_extended(other.submatrix10.T, self.submatrix01.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
self.submatrix11 = rtrg_c.multiply_extended(other.submatrix01.T, self.submatrix10.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix01.T.flags.c_contiguous << 3) | (self.submatrix01.T.flags.f_contiguous << 4) | (other.submatrix10.T.flags.c_contiguous << 1) | (other.submatrix10.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix10.T.flags.c_contiguous << 3) | (self.submatrix10.T.flags.f_contiguous << 4) | (other.submatrix01.T.flags.c_contiguous << 1) | (other.submatrix01.T.flags.f_contiguous << 2)] += 1
self.submatrix01 = None
self.submatrix10 = None
elif (self.submatrix00 is not None and other.submatrix01 is not None):
assert (self.submatrix11 is not None and other.submatrix10 is not None)
self.submatrix01 = rtrg_c.multiply_extended(other.submatrix01.T, self.submatrix00.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
self.submatrix10 = rtrg_c.multiply_extended(other.submatrix10.T, self.submatrix11.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix00.T.flags.c_contiguous << 3) | (self.submatrix00.T.flags.f_contiguous << 4) | (other.submatrix01.T.flags.c_contiguous << 1) | (other.submatrix01.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix11.T.flags.c_contiguous << 3) | (self.submatrix11.T.flags.f_contiguous << 4) | (other.submatrix10.T.flags.c_contiguous << 1) | (other.submatrix10.T.flags.f_contiguous << 2)] += 1
self.submatrix00 = None
self.submatrix11 = None
elif (self.submatrix01 is not None and other.submatrix00 is not None):
assert (self.submatrix10 is not None and other.submatrix11 is not None)
self.submatrix01 = rtrg_c.multiply_extended(other.submatrix11.T, self.submatrix01.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
self.submatrix10 = rtrg_c.multiply_extended(other.submatrix00.T, self.submatrix10.T, self.padding//2, symmetry, self.clear_corners//2, OVERWRITE_LEFT).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix01.T.flags.c_contiguous << 3) | (self.submatrix01.T.flags.f_contiguous << 4) | (other.submatrix11.T.flags.c_contiguous << 1) | (other.submatrix11.T.flags.f_contiguous << 2)] += 1
- SymRGfunction.MMC_COUNTER[(symmetry != 0) | (self.submatrix10.T.flags.c_contiguous << 3) | (self.submatrix10.T.flags.f_contiguous << 4) | (other.submatrix00.T.flags.c_contiguous << 1) | (other.submatrix00.T.flags.f_contiguous << 2)] += 1
return self
def __add__(self, other):
- '''
+ """
Add other to self. If other is a scalar or a scalar function of energy
represented by an array of values at self.energies, this treats other
as an identity (or diagonal) Floquet matrix.
Other must be a scalar or array of same shape as self.energies or an
RGfunction of the same shape and energies as self.
- '''
+ """
if isinstance(other, SymRGfunction):
assert self.global_properties is other.global_properties
symmetry = (self.symmetry == other.symmetry) * self.symmetry
@@ -352,9 +333,9 @@ class SymRGfunction(RGfunction):
raise TypeError("unsupported operand types for +: RGfunction and", type(other))
def __neg__(self):
- '''
+ """
Return a copy of self with inverted sign of self.values.
- '''
+ """
return SymRGfunction(
self.global_properties,
values = None,
@@ -423,12 +404,12 @@ class SymRGfunction(RGfunction):
return self
def __mul__(self, other):
- '''
+ """
Multiplication by scalar or SymRGfunction.
If other is a SymRGfunction, this calculates the matrix product
(equivalent to __matmul__). If other is a scalar, this multiplies
self.values by other.
- '''
+ """
if isinstance(other, RGfunction):
return self @ other
if isinstance(other, Number):
@@ -450,12 +431,12 @@ class SymRGfunction(RGfunction):
return NotImplemented
def __imul__(self, other):
- '''
+ """
In-place multiplication by scalar or SymRGfunction.
If other is a SymRGfunction, this calculates the matrix product
(equivalent to __matmul__). If other is a scalar, this multiplies
self.values by other.
- '''
+ """
if isinstance(other, SymRGfunction):
self @= other
elif isinstance(other, Number):
@@ -477,12 +458,12 @@ class SymRGfunction(RGfunction):
return self
def __rmul__(self, other):
- '''
+ """
Reverse multiplication by scalar or RGfunction.
If other is an RGfunction, this calculates the matrix product
(equivalent to __matmul__). If other is a scalar, this multiplies
self.values by other.
- '''
+ """
if isinstance(other, RGfunction):
return other @ self
if isinstance(other, Number):
@@ -496,17 +477,17 @@ class SymRGfunction(RGfunction):
return NotImplemented
def __truediv__(self, other):
- '''
+ """
Divide self by other, which must be a scalar.
- '''
+ """
if isinstance(other, Number):
return self * (1/other)
return NotImplemented
def __itruediv__(self, other):
- '''
+ """
Divide self in-place by other, which must be a scalar.
- '''
+ """
if isinstance(other, Number):
self *= (1/other)
return self
@@ -533,7 +514,7 @@ class SymRGfunction(RGfunction):
and self.symmetry == other.symmetry
def k2lambda(self, shift_matrix=None):
- '''
+ """
Assume that self is K_n^m(E) = K_n(E-mΩ).
Then calculate Λ_n^m(E) such that (approximately)
@@ -547,7 +528,7 @@ class SymRGfunction(RGfunction):
approximate solution because an exact solution does not exist.
TODO: implement direct energy shift?
- '''
+ """
assert shift_matrix is None
assert self.submatrix01 is None and self.submatrix10 is None
invert = -self
@@ -557,7 +538,7 @@ class SymRGfunction(RGfunction):
return invert.inverse()
def inverse(self):
- '''
+ """
For a given object self = A try to calculate B such that
A @ B = identity
with identity[n,k] = δ(n,0).
@@ -565,16 +546,13 @@ class SymRGfunction(RGfunction):
Some of the linear systems of equation which we need to solve here are
overdetermined. This means that we can in general only get an
approximate solution because an exact solution does not exist.
- '''
+ """
assert self.submatrix01 is None and self.submatrix10 is None
assert self.submatrix00 is not None and self.submatrix11 is not None
try:
res00 = rtrg_c.invert_extended(self.submatrix00.T, self.padding//2, round(settings.LAZY_INVERSE_FACTOR*self.padding/2)).T
res11 = rtrg_c.invert_extended(self.submatrix11.T, self.padding//2, round(settings.LAZY_INVERSE_FACTOR*self.padding/2)).T
- if settings.logger.level == settings.logging.DEBUG:
- SymRGfunction.INVC_COUNTER[self.submatrix00.T.flags.c_contiguous | (self.submatrix00.T.flags.f_contiguous << 1)] += 1
- SymRGfunction.INVC_COUNTER[self.submatrix11.T.flags.c_contiguous | (self.submatrix11.T.flags.f_contiguous << 1)] += 1
except:
settings.logger.exception("padded inversion failed in compact RTRG")
res00 = np.linalg.inv(self.submatrix00)
diff --git a/package/src/frtrg_kondo/data_management.py b/package/src/frtrg_kondo/data_management.py
index 54b7344..01fc352 100644
--- a/package/src/frtrg_kondo/data_management.py
+++ b/package/src/frtrg_kondo/data_management.py
@@ -113,6 +113,13 @@ class KondoExport:
solver_flags |= DataManager.SOLVER_FLAGS["simplified_initial_conditions"]
except AttributeError:
pass
+ try:
+ if self.kondo.truncation_order == 2:
+ solver_flags |= DataManager.SOLVER_FLAGS["second_order_rg_equations"]
+ elif self.kondo.truncation_order != 3:
+ settings.logger.warning("Invalid truncation order: %s"%self.kondo.truncation_order)
+ except AttributeError:
+ pass
for (key, value) in self.kondo.global_settings.items():
if value:
try:
@@ -294,25 +301,27 @@ class KondoExport:
If overwrite is False and a file would be overwritten, append a random
string to the end of the filename.
"""
- os.sync()
- while os.path.exists(filename + '.lock'):
+ while True:
try:
- settings.logger.warning('File %s is locked, waiting 0.5s'%filename)
- sleep(0.5)
- except KeyboardInterrupt:
- answer = input('Ignore lock file? Then type "yes": ')
- if answer.lower() == "yes":
- break
- answer = input('Save with filename extended by random string? (Yn): ')
- if answer.lower()[0] != "n":
- return self.save_h5(filename + random_string(8) + ".h5", include, overwrite)
- pathlib.Path(filename + '.lock').touch()
+ pathlib.Path(filename + '.lock').touch(exist_ok=False)
+ break
+ except FileExistsError:
+ try:
+ settings.logger.warning('File %s is locked, waiting ~0.5s'%filename)
+ sleep(0.4 + 0.2*random.random())
+ except KeyboardInterrupt:
+ answer = input('Ignore lock file? Then type "yes": ')
+ if answer.lower() == "yes":
+ break
+ answer = input('Save with filename extended by random string? (Yn): ')
+ if answer.lower()[0] != "n":
+ return self.save_h5(filename.removesuffix(".h5") + random_string(8) + ".h5", include, overwrite)
try:
file_exists = os.path.exists(filename)
h5file = None
while h5file is None:
try:
- h5file = tb.open_file(filename, "a")
+ h5file = tb.open_file(filename, "a", MAX_NUMEXPR_THREADS=1, MAX_BLOSC_THREADS=1)
except tb.exceptions.HDF5ExtError:
settings.logger.warning('Error opening file %s, waiting 0.5s'%filename)
sleep(0.5)
@@ -320,8 +329,9 @@ class KondoExport:
if file_exists:
try:
h5file.is_visible_node('/data/' + self.hash)
- settings.logger.warning("Hash exists in file %s!"%filename)
- return self.save_h5(filename + random_string(8) + ".h5", include, overwrite)
+ new_filename = filename.removesuffix(".h5") + random_string(8) + ".h5"
+ settings.logger.warning("Hash exists in file %s! Saving to %s"%(filename, new_filename))
+ return self.save_h5(new_filename, include, overwrite)
except tb.exceptions.NoSuchNodeError:
pass
metadata_table = h5file.get_node("/metadata/mdtable")
@@ -414,7 +424,12 @@ class KondoImport:
counter = 0
for row in metadatatable.where(f"hash == '{khash}'"):
metadata = {key:row[key] for key in metadatatable.colnames}
+ metadata.pop("idnum", None)
+ metadata["hash"] = row["hash"].decode()
+ metadata["method"] = KondoExport.METHOD_ENUM(row["method"])
+ metadata["solver_method"] = KondoExport.SOLVER_METHOD_ENUM(row["solver_method"])
item = cls(metadata, datanode, h5file)
+ item._rawmetadata = row
yield item
counter += 1
if counter == 1:
@@ -429,8 +444,13 @@ class KondoImport:
counter = 0
for row in metadatatable:
metadata = {key:row[key] for key in metadatatable.colnames}
- datanode = h5file.get_node('/data/' + row.hash)
+ metadata.pop("idnum", None)
+ metadata["hash"] = row["hash"].decode()
+ metadata["method"] = KondoExport.METHOD_ENUM(row["method"])
+ metadata["solver_method"] = KondoExport.SOLVER_METHOD_ENUM(row["solver_method"])
+ datanode = h5file.get_node("/data/" + metadata["hash"])
item = cls(metadata, datanode, h5file)
+ item._rawmetadata = row
yield item
counter += 1
if counter == 1:
@@ -438,6 +458,59 @@ class KondoImport:
else:
settings.logger.warning("h5file will not be closed automatically")
+ @property
+ def main_results(self):
+ """
+ dictionary of main results: DC current, DC conductance, AC current (absolute value and phase)
+ """
+ results = dict(
+ dc_current = np.nan,
+ dc_conductance = np.nan,
+ ac_current_abs = np.nan,
+ ac_current_phase = np.nan
+ )
+ nmax = self.nmax
+ if self.method in ('unknown', 'mu', 'J', 'mu-reference', 'J-reference', 'mu-extrap-voltage', 'J-extrap-voltage'):
+ try:
+ results['dc_current'] = self.gammaL[nmax, nmax].real
+ except:
+ pass
+ try:
+ results['dc_conductance'] = self.deltaGammaL[nmax, nmax].real
+ except:
+ pass
+ if nmax == 0:
+ results['ac_current_abs'] = 0
+ else:
+ try:
+ results['ac_current_abs'] = np.abs(self.gammaL[nmax-1, nmax])
+ results['ac_current_phase'] = np.angle(self.gammaL[nmax-1, nmax])
+ except:
+ pass
+ elif self.method in ("J-compact-1", "J-compact-2"):
+ results['dc_current'] = 0
+ if nmax % 2:
+ try:
+ results['dc_conductance'] = self.deltaGammaL.submatrix11[nmax//2, nmax//2].real
+ except:
+ pass
+ try:
+ results['ac_current_abs'] = np.abs(self.gammaL.submatrix01[nmax//2, nmax//2])
+ results['ac_current_phase'] = np.angle(self.gammaL.submatrix01[nmax//2, nmax//2])
+ except:
+ pass
+ else:
+ try:
+ results['dc_conductance'] = self.deltaGammaL.submatrix00[nmax//2, nmax//2].real
+ except:
+ pass
+ try:
+ results['ac_current_abs'] = np.abs(self.gammaL.submatrix10[nmax//2-1, nmax//2])
+ results['ac_current_phase'] = np.angle(self.gammaL.submatrix10[nmax//2-1, nmax//2])
+ except:
+ pass
+ return results
+
def __getitem__(self, name):
if name in self.metadata:
return self.metadata[name]
@@ -453,13 +526,12 @@ class KondoImport:
raise AttributeError("Unknown attribute name: %s"%name)
-
class DataManager:
- '''
+ """
Database structure
tables:
datapoints (single data point)
- '''
+ """
SOLVER_FLAGS = dict(
contains_flow = 0x001,
reduced = 0x002,
@@ -471,6 +543,7 @@ class DataManager:
extrapolate_voltage = 0x080,
use_cublas = 0x100,
use_reference_implementation = 0x200,
+ second_order_rg_equations = 0x400,
)
def __init__(self):
@@ -520,10 +593,57 @@ class DataManager:
self.table.create(bind=connection)
def insert_from_h5file(self, filename):
- raise NotImplementedError()
- basename = os.path.basename(filename)
+ """
+ Scan data in HDF5 file and insert datasets in database if they are not
+ included yet.
+ """
+ filename = os.path.abspath(filename)
dirname = os.path.dirname(filename)
- # TODO
+ basename = os.path.basename(filename)
+ datasets = []
+ skipped = 0
+ for dataset in KondoImport.read_all_from_h5(filename):
+ settings.logger.debug("Checking hash=" + dataset.hash)
+ candidates = self.df_table.loc[self.df_table.hash == dataset.hash]
+ if candidates.shape[0] > 0:
+ settings.logger.debug("Found %d times the same hash"%candidates.shape[0])
+ exists = False
+ for idx, candidate in candidates.iterrows():
+ if os.path.join(candidate.dirname, candidate.basename) == filename:
+ exists = True
+ break
+ if exists:
+ settings.logger.debug("Entry exists, skipping")
+ skipped += 1
+ continue
+ else:
+ settings.logger.debug("File seems new, continuing anyway")
+ metadata = dataset.metadata.copy()
+ energy = metadata.pop('energy')
+ metadata.update(
+ energy_re = energy.real,
+ energy_im = energy.imag,
+ timestamp = datetime.fromtimestamp(metadata.pop("timestamp")).isoformat().replace('T', ' '),
+ dirname = dirname,
+ basename = basename,
+ )
+ metadata.update(dataset.main_results)
+ datasets.append(metadata)
+ try:
+ if not dataset._owns_h5file:
+ dataset._h5file.close()
+ settings.logger.info("Closed HDF5 file")
+ except:
+ pass
+ settings.logger.info("Inserting %d new entries, ignoring %d"%(len(datasets), skipped))
+ new_frame = pd.DataFrame(datasets)
+ new_frame.to_sql(
+ 'datapoints',
+ self.engine,
+ if_exists = 'append',
+ index = False,
+ )
+ del self.df_table
def insert_in_db(self, filename : str, kondo : KondoExport):
"""
@@ -551,32 +671,6 @@ class DataManager:
except AttributeError:
pass
- def import_from_db(self, db_string, replace_base_path={}):
- """
- e.g. replace_base_path = {'/path/on/cluster/to/data':'/path/to/local/data'}
- """
- raise NotImplementedError()
- # TODO: rewrite
- import_engine = db.create_engine(db_string, future=True, echo=False)
- import_metadata = db.MetaData()
- import_table = db.Table('datapoints', import_metadata, autoload=True, autoload_with=import_engine)
- with import_engine.begin() as connection:
- import_df_table = pd.read_sql_table('datapoints', connection, index_col='id')
- valid_indices = []
- for idx in import_df_table.index:
- import_df_table.dirname[idx] = replace_all(import_df_table.dirname[idx], replace_base_path)
- # TODO: rewrite this
- selection = self.df_table.basename == import_df_table.basename[idx]
- if not any(self.df_table.dirname[selection] == import_df_table.dirname[idx]):
- valid_indices.append(idx)
- settings.logger.info('Importing %d entries'%len(valid_indices))
- import_df_table.loc[valid_indices].to_sql(
- 'datapoints',
- self.engine,
- if_exists='append',
- index=False,
- )
-
def save_h5(self, kondo : KondoExport, filename : str = None, include='all', overwrite=False):
"""
Save all data in given filename and keep metadata in database.
@@ -589,7 +683,7 @@ class DataManager:
self.insert_in_db(filename, kondo)
def cache_df_table(self, min_version=(0,5,-1)):
- settings.logger.debug('DataManager: cache df_table', flush=True)
+ settings.logger.debug('DataManager: cache df_table')
with self.engine.begin() as connection:
df_table = pd.read_sql_table('datapoints', connection, index_col='id')
selection = (df_table.solver_flags & DataManager.SOLVER_FLAGS['deleted']) == 0
@@ -621,21 +715,29 @@ class DataManager:
return kondo
def list(self, min_version=(14,0,-1,-1), **parameters):
- '''
+ """
Print and return DataFrame with selection of physical parameters.
- '''
+ """
selection = (self.df_table.version_major > min_version[0]) | (self.df_table.version_major == min_version[0]) & (self.df_table.version_minor >= min_version[1])
selection &= self.df_table.energy_re == 0
selection &= self.df_table.energy_im == 0
if len(min_version) > 2 and min_version[2] > 0:
selection &= self.df_table.git_commit_count >= min_version[2]
+ method = parameters.pop("method", None)
+ if method == "J":
+ selection &= self.df_table.method != "mu"
+ elif method is not None:
+ selection &= self.df_table.method == method
for key, value in parameters.items():
if value is None:
continue
try:
- selection &= self.df_table[key] == value
- except KeyError:
- settings.logger.warning("Unknown key: %s"%key)
+ selection &= np.isclose(self.df_table[key], value, rtol=1e-6, atol=1e-15)
+ except TypeError:
+ try:
+ selection &= self.df_table[key] == value
+ except KeyError:
+ settings.logger.warning("Unknown key: %s"%key)
if selection is True:
result = self.df_table
else:
@@ -643,9 +745,9 @@ class DataManager:
return result
def load_from_table(self, table, load_flow=False, load_old_files=True):
- '''
+ """
Extend table by adding a "solver" column.
- '''
+ """
solvers = []
reduced_table = table
for idx, row in table.iterrows():
@@ -662,18 +764,18 @@ class DataManager:
return reduced_table.assign(solver = solvers)
def list_kondo(self, **kwargs):
- '''
+ """
Returns a DataFrame with an extra column "solvers" with the filters
from kwargs applied (see documentation of DataManager.list for the
filters).
- '''
+ """
return self.load_from_table(self.list(**kwargs))
def clean_database(self):
- '''
+ """
Flag all database entries as 'deleted' for which no solver file can be found.
Delete duplicated entries.
- '''
+ """
raise NotImplementedError
with self.engine.begin() as connection:
full_df_table = pd.read_sql_table('datapoints', connection, index_col='id')
diff --git a/package/src/frtrg_kondo/gen_data.py b/package/src/frtrg_kondo/gen_data.py
index b046efc..3194d76 100644
--- a/package/src/frtrg_kondo/gen_data.py
+++ b/package/src/frtrg_kondo/gen_data.py
@@ -18,13 +18,13 @@ from frtrg_kondo.data_management import DataManager
from frtrg_kondo.kondo import Kondo
-def gen_option_iter(steps=None, scales=None, **options):
+def gen_option_iter(steps=None, scale=None, **options):
"""
Interpret given options to swipe over parameters.
Arguments:
steps: number of steps for each swipe dimension
- scales: spacing (linear or logarithmic) for each swipe dimension
+ scale: spacing (linear or logarithmic) for each swipe dimension
**options: arguments for Kondo(...) taken from an argument parser.
These are not the options for Kondo.run(...).
@@ -49,33 +49,33 @@ def gen_option_iter(steps=None, scales=None, **options):
options.pop(key)
steps = [max_length]
else:
- if scales is None:
- scales = len(steps) * ["linear"]
- elif isinstance(scales, str):
- scales = len(steps)*[scales]
- elif len(scales) == 1:
- scales *= len(steps)
+ if scale is None:
+ scale = len(steps) * ["linear"]
+ elif isinstance(scale, str):
+ scale = len(steps)*[scale]
+ elif len(scale) == 1:
+ scale *= len(steps)
for key, value in options.items():
if type(value) != list or key == "fourier_coef":
continue
if len(value) == 1:
options[key], = value
elif len(value) == 2:
- if scales[0] in ("lin", "linear"):
+ if scale[0] in ("lin", "linear"):
iter_variables[key] = (0, np.linspace(value[0], value[1], steps[0], dtype=type(value[0])))
- elif scales[0] in ("log", "logarithmic"):
+ elif scale[0] in ("log", "logarithmic"):
iter_variables[key] = (0, np.logspace(np.log10(value[0]), np.log10(value[1]), steps[0], dtype=type(value[0])))
else:
- raise ValueError("Unexpected value for parameters \"scales\": %s"%scales[0])
+ raise ValueError("Unexpected value for parameters \"scale\": %s"%scale[0])
elif len(value) == 3:
dim = round(value[2])
assert 0 <= dim < len(steps)
- if scales[dim] in ("lin", "linear"):
+ if scale[dim] in ("lin", "linear"):
iter_variables[key] = (dim, np.linspace(value[0], value[1], steps[dim], dtype=type(value[0])))
- elif scales[dim] == "log":
+ elif scale[dim] == "log":
iter_variables[key] = (dim, np.logspace(np.log10(value[0]), np.log10(value[1]), steps[dim], dtype=type(value[0])))
else:
- raise ValueError("Unexpected value for parameters \"scales\": %s"%scales[dim])
+ raise ValueError("Unexpected value for parameters \"scale\": %s"%scale[dim])
else:
raise ValueError("Array parameters must be of the form (start, stop, dim)")
for key in iter_variables.keys():
@@ -207,8 +207,11 @@ python -m frtrg_kondo.gen_data \\
help = "Set initial condition for gammaL to 0")
method_group.add_argument("--d", metavar='float', type=float, nargs='+', default=1e9,
help = "D (UV cutoff), units of Tkrg")
- method_group.add_argument("--resonant_dc_shift", metavar='int', type=int, default=0,
+ method_group.add_argument("--resonant_dc_shift", metavar='int', type=int, nargs='+', default=0,
help = "Describe DC voltage (partially) by shift in Floquet matrices. --vdc is the full voltage.")
+ method_group.add_argument("--truncation_order", metavar='int', type=int, nargs='+',
+ choices=(2,3), default=3,
+ help = "Truncation order of RG equations.")
# Numerical parameters concerning Floquet matrices
numerics_group = parser.add_argument_group(title="Numerical parameters")
@@ -218,7 +221,7 @@ python -m frtrg_kondo.gen_data \\
help = "Floquet matrix ppadding")
numerics_group.add_argument("--voltage_branches", metavar='int', type=int, required=True,
help = "Voltage branches")
- numerics_group.add_argument("--compact", metavar='{0,1,2}', type=int, default=0,
+ numerics_group.add_argument("--compact", metavar='{0,1,2}', type=int, nargs='+', default=0,
help = "compact FRTRG implementation (0, 1, or 2)")
numerics_group.add_argument("--lazy_inverse_factor", metavar='float', type=float,
default = settings.LAZY_INVERSE_FACTOR,
@@ -300,11 +303,14 @@ python -m frtrg_kondo.gen_data \\
# no parallelization
dm = DataManager()
for kondo_options in gen_option_iter(**options):
- vdc = kondo_options['vdc'] + kondo_options['omega']*kondo_options['resonant_dc_shift']
- settings.logger.info(f"Starting with Vdc={vdc}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']}")
+ if kondo_options["voltage_branches"] == 0:
+ vdc = kondo_options["resonant_dc_shift"] * kondo_options["omega"]
+ assert np.abs(kondo_options["vdc"] - vdc) < 1e-9
+ kondo_options["vdc"] = vdc
+ settings.logger.info(f"Starting with Vdc={kondo_options['vdc']}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']}")
kondo = Kondo(**kondo_options)
kondo.run(**solver_options)
- settings.logger.info(f"Saving Vdc={vdc}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']} to {filename}")
+ settings.logger.info(f"Saving Vdc={kondo_options['vdc']}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']} to {filename}")
dm.save_h5(kondo, filename, include)
else:
# generate data points in parallel
@@ -334,13 +340,16 @@ def save_data_mp(queue, lock, solver_options, filename, include='all', overwrite
kondo_options = queue.get()
if kondo_options is None:
break
- vdc = kondo_options['vdc'] + kondo_options['omega']*kondo_options['resonant_dc_shift']
- settings.logger.info(f"Vdc={vdc}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']}")
+ if kondo_options["voltage_branches"] == 0:
+ vdc = kondo_options["resonant_dc_shift"] * kondo_options["omega"]
+ assert np.abs(kondo_options["vdc"] - vdc) < 1e-9
+ kondo_options["vdc"] = vdc
+ settings.logger.info(f"Starting with Vdc={kondo_options['vdc']}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']}")
kondo = Kondo(**kondo_options)
kondo.run(**solver_options)
lock.acquire()
try:
- settings.logger.info(f"Saving Vdc={vdc}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']} to {filename}")
+ settings.logger.info(f"Saving Vdc={kondo_options['vdc']}, Vac={kondo_options['vac']}, Ω={kondo_options['omega']} to {filename}")
dm.save_h5(kondo, filename, include, overwrite)
finally:
lock.release()
diff --git a/package/src/frtrg_kondo/kondo.py b/package/src/frtrg_kondo/kondo.py
index f4b7f32..dfe3a6a 100644
--- a/package/src/frtrg_kondo/kondo.py
+++ b/package/src/frtrg_kondo/kondo.py
@@ -1,5 +1,6 @@
# Copyright 2022 Valentin Bruch <valentin.bruch@rwth-aachen.de>
# License: MIT
+# type: ignore
"""
Kondo FRTRG, main module for RG calculations
@@ -11,8 +12,22 @@ Example usage:
>>> nmax = 10
>>> vb = 3
>>> # Compute the RG flow in 2 different ways
->>> kondo1 = Kondo(unitary_transformation=1, omega=10, nmax=nmax, padding=8, vdc=4, vac=5, voltage_branches=vb)
->>> kondo2 = Kondo(unitary_transformation=0, omega=10, nmax=nmax, padding=0, vdc=4, vac=5, voltage_branches=vb)
+>>> 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
@@ -25,7 +40,8 @@ Example usage:
>>> np.abs(kondo1.z[vb,:,nmax] - kondo2.z[vb,:,nmax]).max()
1.421144031910071e-05
-Further information: https://vbruch.eu/frtrg.html
+Further information:
+https://arxiv.org/abs/2206.06263 and https://vbruch.eu/frtrg.html
"""
import hashlib
@@ -34,15 +50,20 @@ 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
from frtrg_kondo import settings
from frtrg_kondo.rtrg import GlobalRGproperties, RGfunction
-from frtrg_kondo.reservoirmatrix import ReservoirMatrix, einsum_34_12_43, einsum_34_12_43_double, product_combinations
+from frtrg_kondo.reservoirmatrix import ReservoirMatrix, \
+ einsum_34_12_43, \
+ einsum_34_12_43_double, \
+ product_combinations
from frtrg_kondo.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.
@@ -57,6 +78,7 @@ def driving_voltage(tau, *fourier_coef):
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.
@@ -68,7 +90,8 @@ def driving_voltage_integral(tau, *fourier_coef):
return Ω ∫dt' V(t') , t = tau*T
0
- fourier_coef should be in units of Ω, then the result has unit hbar/e = 1
+ 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):
@@ -78,7 +101,8 @@ def driving_voltage_integral(tau, *fourier_coef):
def gen_init_matrix(nmax, *fourier_coef, resonant_dc_shift=0, resolution=5000):
"""
- Generate Floquet matrix of the bare coupling without scalar coupling prefactor.
+ Generate Floquet matrix of the bare coupling without scalar
+ coupling prefactor.
fourier_coef must be in units of Ω:
fourier_coef[n] = V_{n+1}/Ω
@@ -86,7 +110,9 @@ def gen_init_matrix(nmax, *fourier_coef, resonant_dc_shift=0, resolution=5000):
"""
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])
+ 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
@@ -109,7 +135,58 @@ def gen_init_matrix(nmax, *fourier_coef, resonant_dc_shift=0, resolution=5000):
return init_matrix
-def solveTV0_scalar(d, tk=1, rtol=1e-8, atol=1e-10, full_output=False, **solveopts):
+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)
@@ -124,7 +201,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 imaginary axis, ODE of functions of variable lmbd = Λ"""
gamma, theta, j = values
dgamma = theta
dtheta = -4*j**2/(lmbd + gamma)
@@ -145,8 +222,16 @@ 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.
@@ -157,23 +242,31 @@ def solveTV0_Utransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveopt
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)
- gamma, z, j, scalar_solver = solveTV0_scalar(d, **solveopts)
+ 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)
- dj = dtheta*(1 + j/(1 + theta))/2
+ 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
@@ -181,10 +274,12 @@ def solveTV0_Utransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveopt
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.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_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)
@@ -218,8 +313,16 @@ 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
matrix μ used for the voltage shift.
@@ -228,12 +331,17 @@ def solveTV0_untransformed(d, properties, tk=1, rtol=1e-8, atol=1e-10, **solveop
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)
- gamma, z, j, scalar_solver = solveTV0_scalar(d, **solveopts)
+ 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,10 +352,14 @@ 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)
- dj = dtheta*(1 + j/(1 + theta))/2
+ 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))
+ 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))
@@ -285,38 +397,42 @@ 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())
- zj_square = np.einsum('kij,kj,klj->kil', eigvecs, (j_raw*z_raw)**2, eigvecs.conjugate())
-
- return gamma, z, j, zj_square
+ 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.
+ 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
+ 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 : matrix in reservoir space of RG functions, coupling vertex G2
- g3 : matrix in reservoir space of RG functions, coupling vertex G3
- current : matrix in reservoir space fo RG functions, representing the
- current vertex I^L
+ 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:
+ 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
@@ -325,7 +441,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,
@@ -341,41 +457,46 @@ class Kondo:
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.
+ """
+ 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
+ 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.
+ 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)
- 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!
+ 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==True and the
- symmetry V(t+(π/Ω)) = - V(t). Valid values are:
+ 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.
- '''
+ 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,
@@ -387,6 +508,8 @@ class Kondo:
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
@@ -403,13 +526,19 @@ class Kondo:
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 2*abs(resonant_dc_shift) > self.nmax:
- raise ValueError("Bad parameters: resonant_dc_shift must be <= 2*nmax")
- if compact:
+ 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
- assert fourier_coef is None or np.allclose(fourier_coef[1::2], 0)
- assert self.resonant_dc_shift == 0 # extending the implementation to allow for even resonant_dc_shift requires some checks, odd resonant_dc_shift require some rewriting.
+ 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
@@ -427,36 +556,41 @@ class Kondo:
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.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>'
+ """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.
- '''
+ """
+ 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_DC/Ω' : self.resonant_dc_shift,
'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),
- 'IR_cutoff' : getattr(self, 'ir_cutoff', 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.
@@ -464,14 +598,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)
@@ -492,7 +630,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
@@ -502,14 +643,29 @@ 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.
- 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] )
+ 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)
+ 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(
@@ -519,7 +675,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],
@@ -531,7 +687,7 @@ class Kondo:
self.yL = RGfunction(
self.global_properties,
np.zeros((2*self.nmax+1,2*self.nmax+1), dtype=np.complex128),
- symmetry=-symmetry
+ symmetry = -symmetry
)
### Full current, also includes AC current
@@ -544,15 +700,17 @@ class Kondo:
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.
+ **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γ, δΓ, δΓγ.
- '''
+ 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
@@ -560,7 +718,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)
@@ -578,6 +740,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.
@@ -592,11 +757,23 @@ 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)
- 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])
+ 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
@@ -613,14 +790,20 @@ 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
- 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])
+ 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
@@ -634,9 +817,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:
- current_entry = np.diag( 2*sqrtxx * jvalues[diag_idx] * (1 - jvalues[diag_idx] * zvalues[diag_idx] ) )
+ 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:
@@ -649,11 +840,23 @@ 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:
- i0 = j0 * (1 - j0*z0)
- 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])
+ 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
@@ -667,7 +870,9 @@ class Kondo:
## 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),
+ 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
)
@@ -680,14 +885,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])
@@ -708,9 +915,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))
@@ -732,7 +939,6 @@ class Kondo:
self.global_properties.energy = 1j*self.d
-
def initialize(self, **solveopts):
if self.unitary_transformation:
self.initialize_Utransformed(**solveopts)
@@ -746,27 +952,24 @@ 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:
- d = D = UV cutoff (on imaginary axis)
- vac = amplitude of sin(Ωt) driving of bias voltage.
- resonant_dc_shift >= 0, int: add DC voltage of (Ω resonant_dc_shift)
- xL = 1 - xR: asymmetry factor of coupling. Must fulfill 0 <= xL <= 1.
- 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γ, δΓ, δΓγ.
+ 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γ, δΓ, δΓγ.
+ 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
@@ -795,13 +998,13 @@ class Kondo:
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
@@ -822,10 +1025,9 @@ class Kondo:
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.logger.level == settings.logging.DEBUG:
- settings.logger.debug(" ", " ".join("%2s:%4d"%(hex(i)[2:], c) for i,c in enumerate(RGfunction.MM_COUNTER) if c))
- settings.logger.debug(" ", " ".join("%2s:%4d"%(hex(i)[2:], c) for i,c in enumerate(SymRGfunction.MMC_COUNTER) if c))
+ 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()
@@ -858,47 +1060,51 @@ class Kondo:
g2E1tr = g2E1.tr()
g2E1tr.symmetry = -symmetry
if self.compact and not isinstance(g2E1tr, SymRGfunction):
- g2E1tr = SymRGfunction.fromRGfunction(g2E1tr)
+ 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 Z
- # third term for g2E
- # -1/4 G34 ( Π G12 Z + Z G12 Π ) G43
+ ## 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
+ # 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
+ 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 = g2E1 + g2E2 - 0.25*g2_bracket_g2
+ ## 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 = 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, g2_bracket_g2
-
+ 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
-
- # third term for g3E
- # 1/2 G2_34 ( Π G2_12 Z + Z G2_12 Π ) G3_43
- # -> already calculated
+ del g2_pi
- self.g3E = g3E1 + g3E2 + 0.5 * g2_bracket_g3
+ 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, g2_bracket_g3
+ del g3E1, g3E2
# first terms of iE:
@@ -911,38 +1117,39 @@ class Kondo:
assert i_pi[1,1].symmetry == 1
iE1, iE2 = product_combinations( i_pi, self.g2 ) # costs: 12/7/4
- # third term for iE
- # 1/2 I34 ( Π G2_12 Z + Z G2_12 Π ) G43
- # The part in the brackets has been calculated before.
- iE3 = 0.5 * einsum_34_12_43( self.current, pi_g2_z, self.g2 ) # costs: 32/20/10
- del pi_g2_z
-
- self.currentE = iE1 + iE2 + iE3
+ 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
- del iE1, iE2, iE3
# 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
- deltaGammaLE1 = i_pi_g3_01 - i_pi_g3_10
+ self.deltaGammaLE = 1.5j*(i_pi_g3_01 - i_pi_g3_10)
# Second part: I_12 Z Π δΓ G^3_21
- 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
- if self.compact:
- assert isinstance(z_dgamma_pi, SymRGfunction)
- del z_reduced, pi_reduced, dGamma_reduced
- if settings.ENFORCE_SYMMETRIC:
- assert z_dgamma_pi.symmetry == -1
- deltaGammaLE2 = einsum_34_12_43( self.current, z_dgamma_pi, self.g3 ) # costs: 32/20/10
- self.deltaGammaLE = 1.5j*(deltaGammaLE1 + 0.5j*deltaGammaLE2)
+ 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
- del deltaGammaLE1, deltaGammaLE2, z_dgamma_pi
# RG equation for ΓL
@@ -953,36 +1160,36 @@ class Kondo:
if self.compact:
if not isinstance(self.deltaGammaE, SymRGfunction):
- self.deltaGammaE = SymRGfunction.fromRGfunction(self.deltaGammaE)
+ self.deltaGammaE = SymRGfunction.fromRGfunction(self.deltaGammaE, diag=False, offdiag=True)
if not isinstance(self.deltaGammaLE, SymRGfunction):
- self.deltaGammaLE = SymRGfunction.fromRGfunction(self.deltaGammaLE)
+ 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])
+ 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])
+ 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])
+ 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])
+ 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])
+ 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])
+ 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)
+ 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.
- '''
+ """
# Notation (mainly allow using objects without "self")
z = self.z
gamma = self.gamma
@@ -1027,24 +1234,45 @@ class Kondo:
# dE
self.zE = z @ (g2 @ pi @ g2).tr() @ z
- # 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)
+ 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)
- # 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)
+ else:
+ raise ValueError("Invalid truncation order: %s"%self.truncation_order)
# dδΓ ⎛ ⎞
# ——— = i ⎜ δ - δ ⎟ G² Π G²
@@ -1058,12 +1286,19 @@ class Kondo:
pi_reduced = pi.reduced_to_voltage_branches(1)
z_reduced = z.reduced_to_voltage_branches(1)
- # γ
- # 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()
+ 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Γ γ
@@ -1082,9 +1317,9 @@ class Kondo:
def check_symmetry(self):
- '''
+ """
Check if all symmetries are fulfilled
- '''
+ """
self.gamma.check_symmetry()
self.gammaL.check_symmetry()
self.deltaGamma.check_symmetry()
@@ -1097,14 +1332,14 @@ class Kondo:
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
@@ -1212,13 +1447,13 @@ class Kondo:
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,
@@ -1370,14 +1605,14 @@ class Kondo:
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,
@@ -1512,14 +1747,14 @@ class Kondo:
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:
@@ -1559,12 +1794,12 @@ class Kondo:
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()
@@ -1590,7 +1825,7 @@ class Kondo:
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.
@@ -1606,24 +1841,14 @@ class Kondo:
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 self.compact == 0:
- try:
- # Index points to G2[0,0][nmax, nmax, voltage_branches]
- idx = 2*self.gamma.values.size + self.nmax * (2*self.nmax+1) * (2*self.voltage_branches+1) + self.nmax * (2*self.voltage_branches+1) + self.voltage_branches
- except TypeError:
- # Index points to G2[0,0][nmax, nmax]
- idx = 2*self.gamma.values.size + self.nmax * (2*self.nmax+1) + self.nmax
- event = lambda t, y: abs(y[idx]) - g2max
- event.terminal = True
if init_values is None:
init_values = self.packFlattenedValues()
output = solve_ivp(
self.odeFunctionIm,
(eiminit, eimfinal),
init_values,
- events = event if self.compact == 0 else None,
t_eval = ((eimfinal,) if only_final else None),
**solveopts
)
@@ -1631,7 +1856,7 @@ class Kondo:
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.
@@ -1647,23 +1872,14 @@ class Kondo:
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
- try:
- # Index points to G2[0,0][nmax, nmax, voltage_branches]
- idx = 2*self.gamma.values.size + self.nmax * (2*self.nmax+1) * (2*self.voltage_branches+1) + self.nmax * (2*self.voltage_branches+1) + self.voltage_branches
- except TypeError:
- # Index points to G2[0,0][nmax, nmax]
- idx = 2*self.gamma.values.size + self.nmax * (2*self.nmax+1) + self.nmax
- event = lambda t, y: abs(y[idx]) - g2max
- event.terminal = True
if init_values is None:
init_values = self.packFlattenedValues()
output = solve_ivp(
self.odeFunctionRe,
(reEinit, reEfinal),
init_values,
- events = event,
t_eval = ((eimfinal,) if only_final else None),
**solveopts
)
@@ -1671,7 +1887,7 @@ class Kondo:
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
@@ -1679,7 +1895,7 @@ class Kondo:
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:
@@ -1693,10 +1909,10 @@ class Kondo:
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'])
diff --git a/package/src/frtrg_kondo/plot.py b/package/src/frtrg_kondo/plot.py
index becdbfe..bfa5df1 100644
--- a/package/src/frtrg_kondo/plot.py
+++ b/package/src/frtrg_kondo/plot.py
@@ -40,8 +40,12 @@ def main():
Parse command line arguments and call other functions
"""
parser = argparse.ArgumentParser(description=main.__doc__)
- valid_functions = {f.__name__:f for f in globals().values() if type(f) == type(main) and getattr(f, '__module__', '') == '__main__' and f.__name__[0] != '_'}
- parser.add_argument("functions", type=str, nargs='+', choices=valid_functions.keys(), help="functions to be called")
+ valid_functions = {f.__name__:f for f in globals().values()
+ if type(f) == type(main) \
+ and getattr(f, '__module__', '') == '__main__' \
+ and f.__name__[0] != '_'}
+ parser.add_argument("functions", type=str, nargs='+', metavar="function",
+ choices=valid_functions.keys(), help="functions to be called")
parser.add_argument("--db_filename", metavar='file', type=str,
help = "SQLite database file for saving metadata")
parser.add_argument("--log_level", metavar="str", type=str,
@@ -128,16 +132,32 @@ def plot_vdc0(dm, **parameters):
ax4.set_xlabel("Ω")
# DC conductance
ax2.set_title('DC conductance')
- plot = ax2.scatter(results.omega, results.vac, c=np.pi*results.dc_conductance, marker='x', cmap=plt.cm.viridis)
+ plot = ax2.scatter(
+ results.omega,
+ results.vac,
+ c=np.pi*results.dc_conductance,
+ marker='x',
+ cmap=plt.cm.viridis)
fig.colorbar(plot, ax=ax2)
# AC current (abs)
ax3.set_title('AC current')
- plot = ax3.scatter(results.omega, results.vac, c=results.ac_current_abs, marker='x', cmap=plt.cm.viridis)
+ plot = ax3.scatter(
+ results.omega,
+ results.vac,
+ c=results.ac_current_abs,
+ marker='x',
+ cmap=plt.cm.viridis)
fig.colorbar(plot, ax=ax3)
# AC current (phase)
ax4.set_title('AC phase')
phase_norm = plt.Normalize(-np.pi, np.pi)
- plot = ax4.scatter(results.omega, results.vac, c=results.ac_current_phase, marker='+', norm=phase_norm, cmap=plt.cm.hsv)
+ plot = ax4.scatter(
+ results.omega,
+ results.vac,
+ c=results.ac_current_phase,
+ marker='+',
+ norm=phase_norm,
+ cmap=plt.cm.hsv)
fig.colorbar(plot, ax=ax4)
def plot_overview(dm, omega, **parameters):
@@ -156,37 +176,106 @@ def plot_overview(dm, omega, **parameters):
ax4.set_xlabel("Vdc")
# DC current
- idc_min = min(results_J.size and results_J.dc_current.min(), results_mu.size and results_mu.dc_current.min())
- idc_max = max(results_J.size and results_J.dc_current.max(), results_mu.size and results_mu.dc_current.max())
+ idc_min = min(
+ results_J.size and results_J.dc_current.min(),
+ results_mu.size and results_mu.dc_current.min())
+ idc_max = max(
+ results_J.size and results_J.dc_current.max(),
+ results_mu.size and results_mu.dc_current.max())
idc_norm = plt.Normalize(idc_min, idc_max)
ax1.set_title('DC current')
- ax1.scatter(results_J.vdc, results_J.vac, c=results_J.dc_current, marker='x', norm=idc_norm, cmap=plt.cm.viridis)
- plot = ax1.scatter(results_mu.vdc, results_mu.vac, c=results_mu.dc_current, marker='+', norm=idc_norm, cmap=plt.cm.viridis)
+ ax1.scatter(
+ results_J.vdc,
+ results_J.vac,
+ c=results_J.dc_current,
+ marker='x',
+ norm=idc_norm,
+ cmap=plt.cm.viridis)
+ plot = ax1.scatter(
+ results_mu.vdc,
+ results_mu.vac,
+ c=results_mu.dc_current,
+ marker='+',
+ norm=idc_norm,
+ cmap=plt.cm.viridis)
fig.colorbar(plot, ax=ax1)
# DC conductance
- gmin = np.pi*min(results_J.dc_conductance.min() if results_J.size else 1, results_mu.dc_conductance.min() if results_mu.size else 1)
- gmax = np.pi*max(results_J.dc_conductance.max() if results_J.size else 0, results_mu.dc_conductance.max() if results_mu.size else 0)
+ gmin = np.pi*min(
+ results_J.dc_conductance.min() if results_J.size else 1,
+ results_mu.dc_conductance.min() if results_mu.size else 1)
+ gmax = np.pi*max(
+ results_J.dc_conductance.max() if results_J.size else 0,
+ results_mu.dc_conductance.max() if results_mu.size else 0)
gnorm = plt.Normalize(gmin, gmax)
ax2.set_title('DC conductance')
- ax2.scatter(results_J.vdc, results_J.vac, c=np.pi*results_J.dc_conductance, marker='x', norm=gnorm, cmap=plt.cm.viridis)
- plot = ax2.scatter(results_mu.vdc, results_mu.vac, c=np.pi*results_mu.dc_conductance, marker='+', norm=gnorm, cmap=plt.cm.viridis)
+ ax2.scatter(
+ results_J.vdc,
+ results_J.vac,
+ c=np.pi*results_J.dc_conductance,
+ marker='x',
+ norm=gnorm,
+ cmap=plt.cm.viridis)
+ plot = ax2.scatter(
+ results_mu.vdc,
+ results_mu.vac,
+ c=np.pi*results_mu.dc_conductance,
+ marker='+',
+ norm=gnorm,
+ cmap=plt.cm.viridis)
fig.colorbar(plot, ax=ax2)
# AC current (abs)
ax3.set_title('AC current')
- iac_min = min(results_J.size and results_J.ac_current_abs.min(), results_mu.size and results_mu.ac_current_abs.min())
- iac_max = max(results_J.size and results_J.ac_current_abs.max(), results_mu.size and results_mu.ac_current_abs.max())
+ iac_min = min(
+ results_J.size and results_J.ac_current_abs.min(),
+ results_mu.size and results_mu.ac_current_abs.min())
+ iac_max = max(
+ results_J.size and results_J.ac_current_abs.max(),
+ results_mu.size and results_mu.ac_current_abs.max())
iac_norm = plt.Normalize(iac_min, iac_max)
- ax3.scatter(results_J.vdc, results_J.vac, c=results_J.ac_current_abs, marker='x', norm=iac_norm, cmap=plt.cm.viridis)
- plot = ax3.scatter(results_mu.vdc, results_mu.vac, c=results_mu.ac_current_abs, marker='+', norm=iac_norm, cmap=plt.cm.viridis)
+ ax3.scatter(
+ results_J.vdc,
+ results_J.vac,
+ c=results_J.ac_current_abs,
+ marker='x',
+ norm=iac_norm,
+ cmap=plt.cm.viridis)
+ plot = ax3.scatter(
+ results_mu.vdc,
+ results_mu.vac,
+ c=results_mu.ac_current_abs,
+ marker='+',
+ norm=iac_norm,
+ cmap=plt.cm.viridis)
fig.colorbar(plot, ax=ax3)
# AC current (phase)
ax4.set_title('AC phase')
phase_norm = plt.Normalize(-np.pi, np.pi)
- ax4.scatter(results_J.vdc, results_J.vac, c=results_J.ac_current_phase, marker='x', norm=phase_norm, cmap=plt.cm.hsv)
- plot = ax4.scatter(results_mu.vdc, results_mu.vac, c=results_mu.ac_current_phase, marker='+', norm=phase_norm, cmap=plt.cm.hsv)
+ ax4.scatter(
+ results_J.vdc,
+ results_J.vac,
+ c=results_J.ac_current_phase,
+ marker='x',
+ norm=phase_norm,
+ cmap=plt.cm.hsv)
+ plot = ax4.scatter(
+ results_mu.vdc,
+ results_mu.vac,
+ c=results_mu.ac_current_phase,
+ marker='+',
+ norm=phase_norm,
+ cmap=plt.cm.hsv)
fig.colorbar(plot, ax=ax4)
-def plot_interpolate(dm, omega, dc_res=100, ac_res=100, vdc_min=0, vac_min=0, vdc_max=100, vac_max=50, **parameters):
+def plot_interpolate(
+ dm,
+ omega,
+ dc_res=100,
+ ac_res=100,
+ vdc_min=0,
+ vac_min=0,
+ vdc_max=100,
+ vac_max=50,
+ **parameters):
"""
Plot overview of dc and ac current and dc conductance for harmonic driving
at fixed frequency as function of Vdc and Vac.
@@ -208,10 +297,34 @@ def plot_interpolate(dm, omega, dc_res=100, ac_res=100, vdc_min=0, vac_min=0, vd
all_figs = []
all_axes = []
if show_J:
- gdc_J_tck = bisplrep(results.vac[j], results.vdc[j], results.dc_conductance[j], s=1e-5, kx=3, ky=3)
- idc_J_tck = bisplrep(results.vac[j], results.vdc[j], results.dc_current[j], s=1e-5, kx=3, ky=3)
- iac_J_tck = bisplrep(results.vac[j], results.vdc[j], results.ac_current_abs[j], s=1e-5, kx=3, ky=3)
- phase_J_tck = bisplrep(results.vac[j], results.vdc[j], results.ac_current_phase[j], s=1e-5, kx=3, ky=3)
+ gdc_J_tck = bisplrep(
+ results.vac[j],
+ results.vdc[j],
+ results.dc_conductance[j],
+ s=1e-5,
+ kx=3,
+ ky=3)
+ idc_J_tck = bisplrep(
+ results.vac[j],
+ results.vdc[j],
+ results.dc_current[j],
+ s=1e-5,
+ kx=3,
+ ky=3)
+ iac_J_tck = bisplrep(
+ results.vac[j],
+ results.vdc[j],
+ results.ac_current_abs[j],
+ s=1e-5,
+ kx=3,
+ ky=3)
+ phase_J_tck = bisplrep(
+ results.vac[j],
+ results.vdc[j],
+ results.ac_current_phase[j],
+ s=1e-5,
+ kx=3,
+ ky=3)
gdc_g_J_interp = bisplev(vac_arr, vdc_arr, gdc_J_tck)
gdc_i_J_interp = bisplev(vac_arr, vdc_arr, idc_J_tck, dy=1)
idc_J_interp = bisplev(vac_arr, vdc_arr, idc_J_tck)
@@ -226,10 +339,34 @@ def plot_interpolate(dm, omega, dc_res=100, ac_res=100, vdc_min=0, vac_min=0, vd
all_figs.append(fig_j)
all_axes.append((ax1_j, ax2_j, ax3_j, ax4_j))
if show_mu:
- gdc_mu_tck = bisplrep(results.vac[mu], results.vdc[mu], results.dc_conductance[mu], s=5e-6, kx=3, ky=3)
- idc_mu_tck = bisplrep(results.vac[mu], results.vdc[mu], results.dc_current[mu], s=5e-6, kx=3, ky=3)
- iac_mu_tck = bisplrep(results.vac[mu], results.vdc[mu], results.ac_current_abs[mu], s=5e-6, kx=3, ky=3)
- phase_mu_tck = bisplrep(results.vac[mu], results.vdc[mu], results.ac_current_phase[mu], s=5e-6, kx=3, ky=3)
+ gdc_mu_tck = bisplrep(
+ results.vac[mu],
+ results.vdc[mu],
+ results.dc_conductance[mu],
+ s=5e-6,
+ kx=3,
+ ky=3)
+ idc_mu_tck = bisplrep(
+ results.vac[mu],
+ results.vdc[mu],
+ results.dc_current[mu],
+ s=5e-6,
+ kx=3,
+ ky=3)
+ iac_mu_tck = bisplrep(
+ results.vac[mu],
+ results.vdc[mu],
+ results.ac_current_abs[mu],
+ s=5e-6,
+ kx=3,
+ ky=3)
+ phase_mu_tck = bisplrep(
+ results.vac[mu],
+ results.vdc[mu],
+ results.ac_current_phase[mu],
+ s=5e-6,
+ kx=3,
+ ky=3)
gdc_g_mu_interp = bisplev(vac_arr, vdc_arr, gdc_mu_tck)
gdc_i_mu_interp = bisplev(vac_arr, vdc_arr, idc_mu_tck, dy=1)
idc_mu_interp = bisplev(vac_arr, vdc_arr, idc_mu_tck)
@@ -311,26 +448,73 @@ def plot_interpolate(dm, omega, dc_res=100, ac_res=100, vdc_min=0, vac_min=0, vd
ax5.set_title("AC phase: mu vs. J")
check_norm = plt.Normalize(-0.03, 0.03)
if show_J and show_mu:
- ax1.imshow((idc_mu_interp - idc_J_interp)/idc_mu_interp, extent=extent, aspect='auto', cmap=plt.cm.seismic, norm=check_norm, origin='lower')
- ax2.imshow((gdc_g_mu_interp - gdc_g_J_interp)/gdc_g_J_interp, extent=extent, aspect='auto', cmap=plt.cm.seismic, norm=check_norm, origin='lower')
- ax4.imshow((iac_mu_interp - iac_J_interp)/iac_mu_interp, extent=extent, aspect='auto', cmap=plt.cm.seismic, norm=check_norm, origin='lower')
- ax5.imshow(phase_mu_interp - phase_J_interp, extent=extent, aspect='auto', cmap=plt.cm.seismic, norm=check_norm, origin='lower')
+ ax1.imshow(
+ (idc_mu_interp - idc_J_interp)/idc_mu_interp,
+ extent=extent,
+ aspect='auto',
+ cmap=plt.cm.seismic,
+ norm=check_norm,
+ origin='lower')
+ ax2.imshow(
+ (gdc_g_mu_interp - gdc_g_J_interp)/gdc_g_J_interp,
+ extent=extent,
+ aspect='auto',
+ cmap=plt.cm.seismic,
+ norm=check_norm,
+ origin='lower')
+ ax4.imshow(
+ (iac_mu_interp - iac_J_interp)/iac_mu_interp,
+ extent=extent,
+ aspect='auto',
+ cmap=plt.cm.seismic,
+ norm=check_norm,
+ origin='lower')
+ ax5.imshow(
+ phase_mu_interp - phase_J_interp,
+ extent=extent,
+ aspect='auto',
+ cmap=plt.cm.seismic,
+ norm=check_norm,
+ origin='lower')
if show_mu:
- img = ax3.imshow((gdc_g_mu_interp - gdc_i_mu_interp)/gdc_g_mu_interp, extent=extent, aspect='auto', cmap=plt.cm.seismic, norm=check_norm, origin='lower')
+ img = ax3.imshow(
+ (gdc_g_mu_interp - gdc_i_mu_interp)/gdc_g_mu_interp,
+ extent=extent,
+ aspect='auto',
+ cmap=plt.cm.seismic,
+ norm=check_norm,
+ origin='lower')
if show_J:
- img = ax6.imshow((gdc_g_J_interp - gdc_i_J_interp)/gdc_g_J_interp, extent=extent, aspect='auto', cmap=plt.cm.seismic, norm=check_norm, origin='lower')
+ img = ax6.imshow(
+ (gdc_g_J_interp - gdc_i_J_interp)/gdc_g_J_interp,
+ extent=extent,
+ aspect='auto',
+ cmap=plt.cm.seismic,
+ norm=check_norm,
+ origin='lower')
fig.colorbar(img, ax=[ax1,ax2,ax3,ax4,ax5,ax6])
-def plot_comparison(dm, omega, **trashoptions):
+def plot_comparison(dm, omega, d=1e9, **parameters):
"""
Compare current and dc conductance computed from both methods (J and mu) at
fixed frequency as function of Vdc and Vac.
Only data points which exist for both methods with equal Vdc, Vac, D and
frequency are considered.
"""
- results_J = dm.list(omega=omega, method='J')
- results_mu = dm.list(omega=omega, method='mu')
- merged = pd.merge(results_J, results_mu, how="inner", on=["vdc","vac","d","omega"], suffixes=("_J", "_mu"))
+ for name in ("omega", "method"):
+ parameters.pop(name, None)
+ results_J = dm.list(omega=omega, method='J', d=d, **parameters).sort_values(["vdc","vac"])
+ results_mu = dm.list(omega=omega, method='mu', d=d, **parameters).sort_values(["vdc","vac"])
+ results_J.vac = results_J.vac.round(6)
+ results_J.vdc = results_J.vdc.round(6)
+ results_mu.vac = results_mu.vac.round(6)
+ results_mu.vdc = results_mu.vdc.round(6)
+ merged = pd.merge(
+ results_J,
+ results_mu,
+ how="inner",
+ on=["vdc","vac","d","omega"],
+ suffixes=("_J", "_mu"))
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, sharey=True)
ax1.set_ylabel("Vac")
@@ -367,16 +551,27 @@ def plot_comparison(dm, omega, **trashoptions):
plot = ax4.scatter(merged.vdc, merged.vac, c=phase_diff, marker='o', norm=phase_norm, cmap=plt.cm.seismic)
fig.colorbar(plot, ax=ax4)
-def plot_comparison_relative(dm, omega, **trashoptions):
+def plot_comparison_relative(dm, omega, d=1e9, **parameters):
"""
Compare current and dc conductance computed from both methods (J and mu) at
fixed frequency as function of Vdc and Vac.
Only data points which exist for both methods with equal Vdc, Vac, D and
frequency are considered.
"""
- results_J = dm.list(omega=omega, method='J')
- results_mu = dm.list(omega=omega, method='mu')
- merged = pd.merge(results_J, results_mu, how="inner", on=["vdc","vac","d","omega"], suffixes=("_J", "_mu"))
+ for name in ("omega", "method"):
+ parameters.pop(name, None)
+ results_J = dm.list(omega=omega, method='J', d=d, **parameters).sort_values(["vdc","vac"])
+ results_mu = dm.list(omega=omega, method='mu', d=d, **parameters).sort_values(["vdc","vac"])
+ results_J.vac = results_J.vac.round(6)
+ results_J.vdc = results_J.vdc.round(6)
+ results_mu.vac = results_mu.vac.round(6)
+ results_mu.vdc = results_mu.vdc.round(6)
+ merged = pd.merge(
+ results_J,
+ results_mu,
+ how="inner",
+ on=["vdc","vac","d","omega"],
+ suffixes=("_J", "_mu"))
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, sharex=True, sharey=True)
ax1.set_ylabel("Vac")
@@ -531,13 +726,29 @@ def check_convergence(dm, **parameters):
if reference_values is not None:
x = (table.d.min(), table.d.max())
ax1.plot(x, (reference_values['idc'],reference_values['idc']), 'k:')
- ax1.fill_between(x, (reference_values['idc']-reference_values['idc_err'],reference_values['idc']-reference_values['idc_err']), (reference_values['idc']+reference_values['idc_err'],reference_values['idc']+reference_values['idc_err']), color='#80808080')
+ ax1.fill_between(
+ x,
+ (reference_values['idc']-reference_values['idc_err'], reference_values['idc']-reference_values['idc_err']),
+ (reference_values['idc']+reference_values['idc_err'], reference_values['idc']+reference_values['idc_err']),
+ color='#80808080')
ax2.plot(x, (reference_values['gdc'],reference_values['gdc']), 'k:')
- ax2.fill_between(x, (reference_values['gdc']-reference_values['gdc_err'],reference_values['gdc']-reference_values['gdc_err']), (reference_values['gdc']+reference_values['gdc_err'],reference_values['gdc']+reference_values['gdc_err']), color='#80808080')
+ ax2.fill_between(
+ x,
+ (reference_values['gdc']-reference_values['gdc_err'], reference_values['gdc']-reference_values['gdc_err']),
+ (reference_values['gdc']+reference_values['gdc_err'], reference_values['gdc']+reference_values['gdc_err']),
+ color='#80808080')
ax3.plot(x, (reference_values['iac'],reference_values['iac']), 'k:')
- ax3.fill_between(x, (reference_values['iac']-reference_values['iac_err'],reference_values['iac']-reference_values['iac_err']), (reference_values['iac']+reference_values['iac_err'],reference_values['iac']+reference_values['iac_err']), color='#80808080')
+ ax3.fill_between(
+ x,
+ (reference_values['iac']-reference_values['iac_err'], reference_values['iac']-reference_values['iac_err']),
+ (reference_values['iac']+reference_values['iac_err'], reference_values['iac']+reference_values['iac_err']),
+ color='#80808080')
ax4.plot(x, (reference_values['acphase'],reference_values['acphase']), 'k:')
- ax4.fill_between(x, (reference_values['acphase']-reference_values['acphase_err'],reference_values['acphase']-reference_values['acphase_err']), (reference_values['acphase']+reference_values['acphase_err'],reference_values['acphase']+reference_values['acphase_err']), color='#80808080')
+ ax4.fill_between(
+ x,
+ (reference_values['acphase']-reference_values['acphase_err'], reference_values['acphase']-reference_values['acphase_err']),
+ (reference_values['acphase']+reference_values['acphase_err'], reference_values['acphase']+reference_values['acphase_err']),
+ color='#80808080')
ax1.set_title("DC current")
ax1.scatter(d_j, table.dc_current[j], marker='x', s=s_j, c=c_j, linewidths=lw_j, norm=norm)
@@ -593,13 +804,29 @@ def check_convergence_nmax(dm, **parameters):
if reference_values is not None:
x = (table.nmax.min(), table.nmax.max())
ax1.plot(x, (reference_values['idc'],reference_values['idc']), 'k:')
- ax1.fill_between(x, (reference_values['idc']-reference_values['idc_err'],reference_values['idc']-reference_values['idc_err']), (reference_values['idc']+reference_values['idc_err'],reference_values['idc']+reference_values['idc_err']), color='#80808080')
+ ax1.fill_between(
+ x,
+ (reference_values['idc']-reference_values['idc_err'], reference_values['idc']-reference_values['idc_err']),
+ (reference_values['idc']+reference_values['idc_err'], reference_values['idc']+reference_values['idc_err']),
+ color='#80808080')
ax2.plot(x, (reference_values['gdc'],reference_values['gdc']), 'k:')
- ax2.fill_between(x, (reference_values['gdc']-reference_values['gdc_err'],reference_values['gdc']-reference_values['gdc_err']), (reference_values['gdc']+reference_values['gdc_err'],reference_values['gdc']+reference_values['gdc_err']), color='#80808080')
+ ax2.fill_between(
+ x,
+ (reference_values['gdc']-reference_values['gdc_err'], reference_values['gdc']-reference_values['gdc_err']),
+ (reference_values['gdc']+reference_values['gdc_err'], reference_values['gdc']+reference_values['gdc_err']),
+ color='#80808080')
ax3.plot(x, (reference_values['iac'],reference_values['iac']), 'k:')
- ax3.fill_between(x, (reference_values['iac']-reference_values['iac_err'],reference_values['iac']-reference_values['iac_err']), (reference_values['iac']+reference_values['iac_err'],reference_values['iac']+reference_values['iac_err']), color='#80808080')
+ ax3.fill_between(
+ x,
+ (reference_values['iac']-reference_values['iac_err'], reference_values['iac']-reference_values['iac_err']),
+ (reference_values['iac']+reference_values['iac_err'], reference_values['iac']+reference_values['iac_err']),
+ color='#80808080')
ax4.plot(x, (reference_values['acphase'],reference_values['acphase']), 'k:')
- ax4.fill_between(x, (reference_values['acphase']-reference_values['acphase_err'],reference_values['acphase']-reference_values['acphase_err']), (reference_values['acphase']+reference_values['acphase_err'],reference_values['acphase']+reference_values['acphase_err']), color='#80808080')
+ ax4.fill_between(
+ x,
+ (reference_values['acphase']-reference_values['acphase_err'], reference_values['acphase']-reference_values['acphase_err']),
+ (reference_values['acphase']+reference_values['acphase_err'], reference_values['acphase']+reference_values['acphase_err']),
+ color='#80808080')
ax1.set_title("DC current")
ax1.scatter(nmax_j, table.dc_current[j], marker='x', s=s_j, c=c_j, linewidths=lw_j, norm=norm)
@@ -621,9 +848,11 @@ def check_convergence_nmax(dm, **parameters):
ax4.scatter(nmax_mu, table.ac_current_phase[mu], marker='+', s=s_mu, c=c_mu, linewidths=lw_mu, norm=norm)
ax4.scatter(nmax_mod, table.ac_current_phase[mod], marker='_', s=s_mod, c=c_mod, linewidths=lw_mod, norm=norm)
+
def check_convergence_fit(dm, **parameters):
if not ("omega" in parameters and "vdc" in parameters and "vac" in parameters):
- settings.logger.warning("check_convergence expects specification of all physical parameters")
+ settings.logger.warning(
+ "check_convergence expects specification of all physical parameters")
table = dm.list(**parameters)
mod = (table.solver_flags & DataManager.SOLVER_FLAGS["simplified_initial_conditions"]) != 0
d_fit_max = 9e11
@@ -654,10 +883,24 @@ def check_convergence_fit(dm, **parameters):
for r in range(10):
suffix = '_%d'%r if r else ''
if r in table.resonant_dc_shift.values:
- mu_shift = (~mod) & (table.method == "mu") & (table.d <= d_fit_max) & (table.resonant_dc_shift == r)
- mu_mod_shift = mod & (table.method == "mu") & (table.d <= d_fit_max) & (table.resonant_dc_shift == r)
- j_shift = (~mod) & (table.method == "J") & (((table.padding > 0) & (table.d <= d_fit_max)) | (table.d <= (d_fit_max_j_nopadding if r == 0 else d_fit_max_j_shifted_nopadding))) & (table.resonant_dc_shift == r)
- j_mod_shift = mod & (table.method == "J") & (((table.padding > 0) & (table.d <= d_fit_max)) | (table.d <= (d_fit_max_j_nopadding if r == 0 else d_fit_max_j_shifted_nopadding))) & (table.resonant_dc_shift == r)
+ mu_shift = (~mod) \
+ & (table.method == "mu") \
+ & (table.d <= d_fit_max) \
+ & (table.resonant_dc_shift == r)
+ mu_mod_shift = mod \
+ & (table.method == "mu") \
+ & (table.d <= d_fit_max) \
+ & (table.resonant_dc_shift == r)
+ j_shift = (~mod) \
+ & (table.method == "J") \
+ & ( ((table.padding > 0) & (table.d <= d_fit_max)) \
+ | (table.d <= (d_fit_max_j_nopadding if r == 0 else d_fit_max_j_shifted_nopadding)) ) \
+ & (table.resonant_dc_shift == r)
+ j_mod_shift = mod \
+ & (table.method == "J") \
+ & ( ((table.padding > 0) & (table.d <= d_fit_max)) \
+ | (table.d <= (d_fit_max_j_nopadding if r == 0 else d_fit_max_j_shifted_nopadding)) ) \
+ & (table.resonant_dc_shift == r)
if mu_shift.any():
selections["mu" + suffix] = mu_shift
if mu_mod_shift.any():
@@ -688,24 +931,77 @@ def check_convergence_fit(dm, **parameters):
# Show reference values
try:
- reference_values = REFERENCE_VALUES[(parameters['omega'], parameters['vdc'], parameters['vac'])]
+ reference_values = REFERENCE_VALUES[(
+ parameters['omega'],
+ parameters['vdc'],
+ parameters['vac'],
+ )]
except KeyError:
reference_values = None
if reference_values is not None:
- ax1.plot((0,logd_inv_arr[-1]), (reference_values['idc'],reference_values['idc']), 'k:')
- ax1.fill_between((0,logd_inv_arr[-1]), (reference_values['idc']-reference_values['idc_err'],reference_values['idc']-reference_values['idc_err']), (reference_values['idc']+reference_values['idc_err'],reference_values['idc']+reference_values['idc_err']), color='#80808080')
- ax2.plot((0,logd_inv_arr[-1]), (reference_values['gdc'],reference_values['gdc']), 'k:')
- ax2.fill_between((0,logd_inv_arr[-1]), (reference_values['gdc']-reference_values['gdc_err'],reference_values['gdc']-reference_values['gdc_err']), (reference_values['gdc']+reference_values['gdc_err'],reference_values['gdc']+reference_values['gdc_err']), color='#80808080')
- ax3.plot((0,logd_inv_arr[-1]), (reference_values['iac'],reference_values['iac']), 'k:')
- ax3.fill_between((0,logd_inv_arr[-1]), (reference_values['iac']-reference_values['iac_err'],reference_values['iac']-reference_values['iac_err']), (reference_values['iac']+reference_values['iac_err'],reference_values['iac']+reference_values['iac_err']), color='#80808080')
- ax4.plot((0,logd_inv_arr[-1]), (reference_values['acphase'],reference_values['acphase']), 'k:')
- ax4.fill_between((0,logd_inv_arr[-1]), (reference_values['acphase']-reference_values['acphase_err'],reference_values['acphase']-reference_values['acphase_err']), (reference_values['acphase']+reference_values['acphase_err'],reference_values['acphase']+reference_values['acphase_err']), color='#80808080')
+ ax1.plot(
+ (0, logd_inv_arr[-1]),
+ (reference_values['idc'], reference_values['idc']),
+ 'k:')
+ ax1.fill_between(
+ (0, logd_inv_arr[-1]),
+ (reference_values['idc']-reference_values['idc_err'], reference_values['idc']-reference_values['idc_err']),
+ (reference_values['idc']+reference_values['idc_err'], reference_values['idc']+reference_values['idc_err']),
+ color='#80808080')
+ ax2.plot(
+ (0, logd_inv_arr[-1]),
+ (reference_values['gdc'], reference_values['gdc']),
+ 'k:')
+ ax2.fill_between(
+ (0, logd_inv_arr[-1]),
+ (reference_values['gdc']-reference_values['gdc_err'], reference_values['gdc']-reference_values['gdc_err']),
+ (reference_values['gdc']+reference_values['gdc_err'], reference_values['gdc']+reference_values['gdc_err']),
+ color='#80808080')
+ ax3.plot(
+ (0, logd_inv_arr[-1]),
+ (reference_values['iac'], reference_values['iac']),
+ 'k:')
+ ax3.fill_between(
+ (0, logd_inv_arr[-1]),
+ (reference_values['iac']-reference_values['iac_err'], reference_values['iac']-reference_values['iac_err']),
+ (reference_values['iac']+reference_values['iac_err'], reference_values['iac']+reference_values['iac_err']),
+ color='#80808080')
+ ax4.plot(
+ (0, logd_inv_arr[-1]),
+ (reference_values['acphase'], reference_values['acphase']),
+ 'k:')
+ ax4.fill_between(
+ (0, logd_inv_arr[-1]),
+ (reference_values['acphase']-reference_values['acphase_err'], reference_values['acphase']-reference_values['acphase_err']),
+ (reference_values['acphase']+reference_values['acphase_err'], reference_values['acphase']+reference_values['acphase_err']),
+ color='#80808080')
ax1.set_title("DC current")
- ax1.scatter(x_j, table.dc_current[j], marker='x', s=s_j, c=c_j, linewidths=lw_j, norm=norm)
- ax1.scatter(x_mu, table.dc_current[mu], marker='+', s=s_mu, c=c_mu, linewidths=lw_mu, norm=norm)
- ax1.scatter(x_mod, table.dc_current[mod], marker='*', s=s_mod, c=c_mod, linewidths=lw_mod, norm=norm)
+ ax1.scatter(
+ x_j,
+ table.dc_current[j],
+ marker = 'x',
+ s = s_j,
+ c = c_j,
+ linewidths = lw_j,
+ norm = norm)
+ ax1.scatter(
+ x_mu,
+ table.dc_current[mu],
+ marker = '+',
+ s = s_mu,
+ c = c_mu,
+ linewidths = lw_mu,
+ norm = norm)
+ ax1.scatter(
+ x_mod,
+ table.dc_current[mod],
+ marker = '*',
+ s = s_mod,
+ c = c_mod,
+ linewidths = lw_mod,
+ norm = norm)
bb = table.dc_current[~mod].mean()
aa = (x_shifted*table.dc_current[~mod]).sum()/s_xx
@@ -717,7 +1013,11 @@ def check_convergence_fit(dm, **parameters):
for name, selection in selections.items():
try:
- (a, b, c), covar = curve_fit(fit_func, logd_inv[selection], table.dc_current[selection], (bb-aa*x_mean, aa, 3))
+ (a, b, c), covar = curve_fit(
+ fit_func,
+ logd_inv[selection],
+ table.dc_current[selection],
+ (bb-aa*x_mean, aa, 3))
aerr, berr, cerr = covar.diagonal()**0.5
print(f"DC current ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
ax1.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
@@ -726,14 +1026,56 @@ def check_convergence_fit(dm, **parameters):
pass
ax2.set_title("DC conductance")
- ax2.scatter(x_j, table.dc_conductance[j], marker='x', s=s_j, c=c_j, linewidths=lw_j, norm=norm)
- ax2.scatter(x_mu, table.dc_conductance[mu], marker='+', s=s_mu, c=c_mu, linewidths=lw_mu, norm=norm)
- ax2.scatter(x_mod, table.dc_conductance[mod], marker='*', s=s_mod, c=c_mod, linewidths=lw_mod, norm=norm)
+ ax2.scatter(
+ x_j,
+ table.dc_conductance[j],
+ marker = 'x',
+ s = s_j,
+ c = c_j,
+ linewidths = lw_j,
+ norm = norm)
+ ax2.scatter(
+ x_mu,
+ table.dc_conductance[mu],
+ marker = '+',
+ s = s_mu,
+ c = c_mu,
+ linewidths = lw_mu,
+ norm = norm)
+ ax2.scatter(
+ x_mod,
+ table.dc_conductance[mod],
+ marker = '*',
+ s = s_mod,
+ c = c_mod,
+ linewidths = lw_mod,
+ norm = norm)
ax3.set_title("AC current")
- ax3.scatter(x_j, table.ac_current_abs[j], marker='x', s=s_j, c=c_j, linewidths=lw_j, norm=norm)
- ax3.scatter(x_mu, table.ac_current_abs[mu], marker='+', s=s_mu, c=c_mu, linewidths=lw_mu, norm=norm)
- ax3.scatter(x_mod, table.ac_current_abs[mod], marker='*', s=s_mod, c=c_mod, linewidths=lw_mod, norm=norm)
+ ax3.scatter(
+ x_j,
+ table.ac_current_abs[j],
+ marker = 'x',
+ s = s_j,
+ c = c_j,
+ linewidths = lw_j,
+ norm = norm)
+ ax3.scatter(
+ x_mu,
+ table.ac_current_abs[mu],
+ marker = '+',
+ s = s_mu,
+ c = c_mu,
+ linewidths = lw_mu,
+ norm = norm)
+ ax3.scatter(
+ x_mod,
+ table.ac_current_abs[mod],
+ marker = '*',
+ s = s_mod,
+ c = c_mod,
+ linewidths = lw_mod,
+ norm = norm)
b = table.ac_current_abs[~mod].mean()
a = (x_shifted*table.ac_current_abs[~mod]).sum()/s_xx
@@ -745,7 +1087,11 @@ def check_convergence_fit(dm, **parameters):
for name, selection in selections.items():
try:
- (a, b, c), covar = curve_fit(fit_func, logd_inv[selection], table.ac_current_abs[selection], (bb-aa*x_mean, aa, 3))
+ (a, b, c), covar = curve_fit(
+ fit_func,
+ logd_inv[selection],
+ table.ac_current_abs[selection],
+ (bb-aa*x_mean, aa, 3))
aerr, berr, cerr = covar.diagonal()**0.5
print(f"AC current ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
ax3.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
@@ -754,9 +1100,30 @@ def check_convergence_fit(dm, **parameters):
pass
ax4.set_title("AC phase")
- ax4.scatter(x_j, table.ac_current_phase[j], marker='x', s=s_j, c=c_j, linewidths=lw_j, norm=norm)
- ax4.scatter(x_mu, table.ac_current_phase[mu], marker='+', s=s_mu, c=c_mu, linewidths=lw_mu, norm=norm)
- ax4.scatter(x_mod, table.ac_current_phase[mod], marker='*', s=s_mod, c=c_mod, linewidths=lw_mod, norm=norm)
+ ax4.scatter(
+ x_j,
+ table.ac_current_phase[j],
+ marker='x',
+ s=s_j,
+ c=c_j,
+ linewidths=lw_j,
+ norm=norm)
+ ax4.scatter(
+ x_mu,
+ table.ac_current_phase[mu],
+ marker='+',
+ s=s_mu,
+ c=c_mu,
+ linewidths=lw_mu,
+ norm=norm)
+ ax4.scatter(
+ x_mod,
+ table.ac_current_phase[mod],
+ marker='*',
+ s=s_mod,
+ c=c_mod,
+ linewidths=lw_mod,
+ norm=norm)
b = table.ac_current_phase[~mod].mean()
a = (x_shifted*table.ac_current_phase[~mod]).sum()/s_xx
@@ -768,7 +1135,11 @@ def check_convergence_fit(dm, **parameters):
for name, selection in selections.items():
try:
- (a, b, c), covar = curve_fit(fit_func, logd_inv[selection], table.ac_current_phase[selection], (bb-aa*x_mean, aa, 3))
+ (a, b, c), covar = curve_fit(
+ fit_func,
+ logd_inv[selection],
+ table.ac_current_phase[selection],
+ (bb-aa*x_mean, aa, 3))
aerr, berr, cerr = covar.diagonal()**0.5
print(f"AC phase ({name:9}): a={a:.9}±{aerr:.9}, b={b:.9}±{berr:.9}, c={c:.9}±{cerr:.9}")
ax4.plot(logd_inv_arr, fit_func(logd_inv_arr, a, b, c), label=name)
diff --git a/package/src/frtrg_kondo/reservoirmatrix.py b/package/src/frtrg_kondo/reservoirmatrix.py
index c685c66..03c627d 100644
--- a/package/src/frtrg_kondo/reservoirmatrix.py
+++ b/package/src/frtrg_kondo/reservoirmatrix.py
@@ -16,10 +16,10 @@ from frtrg_kondo.rtrg import RGobj, RGfunction
class ReservoirMatrix(RGobj):
- '''
+ """
2x2 matrix of RGfunctions.
- This includes a system of book keeping for energy shifts by multiples of
- the voltage in products of ReservoirMatrices.
+ This includes a system of book keeping for energy shifts by
+ multiples of the voltage in products of ReservoirMatrices.
if symmetry != 0:
self.data[1,0] == symmetry * self.data[0,1].floquetConjugate()
if symmety != 0 and global_properties.symmetric:
@@ -27,14 +27,15 @@ class ReservoirMatrix(RGobj):
Multiplication operators for ReservoirMatrix objects are:
* for multiplication with scalars
- @ for multiplication with RGfunctions and normal matrix multiplication
- with ReservoirMatrices (matrix indices ij, jk → ik with sum over j)
+ @ for multiplication with RGfunctions and normal matrix
+ multiplication with ReservoirMatrices
+ (matrix indices ij, jk → ik with sum over j)
% for transpose matrix multiplication with ReservoirMatrices
(matrix indices jk, ij → ik with sum over j)
⎛ ⊤ ⊤⎞⊤
- i.e. A % B = ⎜A @ B ⎟ when there are no energy argument shifts.
+ i.e. A % B = ⎜A @ B ⎟ if there are no energy argument shifts.
⎝ ⎠
- '''
+ """
def __init__(self, global_properties, symmetry=0):
super().__init__(global_properties, symmetry)
@@ -49,7 +50,8 @@ class ReservoirMatrix(RGobj):
self.data.__setitem__(indices, value)
if isinstance(value, RGfunction):
assert value.global_properties is self.global_properties
- self.data.__getitem__(indices).voltage_shifts = self.voltage_shifts + indices[1] - indices[0]
+ self.data.__getitem__(indices).voltage_shifts = \
+ self.voltage_shifts + indices[1] - indices[0]
def __add__(self, other):
if isinstance(other, ReservoirMatrix):
@@ -61,7 +63,9 @@ class ReservoirMatrix(RGobj):
res.data = self.data + other.data
return res
else:
- raise NotImplementedError('Addition is not defined for types %s and %s'%(ReservoirMatrix, type(other)))
+ raise NotImplementedError(
+ 'Addition is not defined for types %s and %s'%(
+ ReservoirMatrix, type(other)))
def __iadd__(self, other):
if isinstance(other, ReservoirMatrix):
@@ -71,7 +75,9 @@ class ReservoirMatrix(RGobj):
if self.symmetry != other.symmetry:
self.symmetry = 0
else:
- raise NotImplementedError('Addition is not defined for types %s and %s'%(ReservoirMatrix, type(other)))
+ raise NotImplementedError(
+ 'Addition is not defined for types %s and %s'%(
+ ReservoirMatrix, type(other)))
return self
def __sub__(self, other):
@@ -84,7 +90,9 @@ class ReservoirMatrix(RGobj):
res.data = self.data - other.data
return res
else:
- raise NotImplementedError('Subtraction is not defined for types %s and %s'%(ReservoirMatrix, type(other)))
+ raise NotImplementedError(
+ 'Subtraction is not defined for types %s and %s'%(
+ ReservoirMatrix, type(other)))
def __isub__(self, other):
if isinstance(other, ReservoirMatrix):
@@ -94,7 +102,9 @@ class ReservoirMatrix(RGobj):
if self.symmetry != other.symmetry:
self.symmetry = 0
else:
- raise NotImplementedError('Subtraction is not defined for types %s and %s'%(ReservoirMatrix, type(other)))
+ raise NotImplementedError(
+ 'Subtraction is not defined for types %s and %s'%(
+ ReservoirMatrix, type(other)))
return self
def __neg__(self):
@@ -105,7 +115,7 @@ class ReservoirMatrix(RGobj):
def __mul__(self, other):
if isinstance(other, ReservoirMatrix) or isinstance(other, RGfunction):
- raise TypeError("Multiplication of reservoir matrices uses the @ symbol.")
+ raise TypeError('Multiplication of reservoir matrices uses the @ symbol.')
else:
res = ReservoirMatrix(self.global_properties)
if other.imag == 0:
@@ -119,7 +129,7 @@ class ReservoirMatrix(RGobj):
def __imul__(self, other):
if isinstance(other, ReservoirMatrix):
- raise TypeError("Multiplication of reservoir matrices uses the @ symbol.")
+ raise TypeError('Multiplication of reservoir matrices uses the @ symbol.')
else:
if other.imag != 0:
if other.real == 0:
@@ -131,7 +141,7 @@ class ReservoirMatrix(RGobj):
def __rmul__(self, other):
if isinstance(other, ReservoirMatrix):
- raise TypeError("Multiplication of reservoir matrices uses the @ symbol.")
+ raise TypeError('Multiplication of reservoir matrices uses the @ symbol.')
elif isinstance(other, RGfunction):
res = ReservoirMatrix(self.global_properties, other.symmetry * self.symmetry)
res.data = other * self.data
@@ -176,7 +186,9 @@ class ReservoirMatrix(RGobj):
# 3 multiplications with symmetry but xL != xR
# 2 multiplications with symmetry and xL == xR
assert self.global_properties is other.global_properties
- res = ReservoirMatrix(self.global_properties, self.symmetry * other.symmetry * (self.voltage_shifts == 0))
+ res = ReservoirMatrix(
+ self.global_properties,
+ self.symmetry * other.symmetry * (self.voltage_shifts == 0))
res.voltage_shifts = self.voltage_shifts + other.voltage_shifts
res[0,0] = self[0,0] @ other
if res.symmetry == 0 or settings.IGNORE_SYMMETRIES:
@@ -192,13 +204,17 @@ class ReservoirMatrix(RGobj):
res[1,0] = res.symmetry * res[0,1].floquetConjugate()
return res
else:
- raise TypeError('Math multiplication is not defined for types %s and %s'%(ReservoirMatrix, type(other)))
+ raise TypeError(
+ 'Math multiplication is not defined for types %s and %s'%(
+ ReservoirMatrix, type(other)))
def __rmatmul__(self, other):
if isinstance(other, RGfunction):
assert self.voltage_shifts == 0
assert self.global_properties is other.global_properties
- res = ReservoirMatrix(self.global_properties, self.symmetry * other.symmetry * (other.voltage_shifts == 0))
+ res = ReservoirMatrix(
+ self.global_properties,
+ self.symmetry * other.symmetry * (other.voltage_shifts == 0))
res.voltage_shifts = other.voltage_shifts
res[0,0] = other @ self[0,0]
if res.symmetry == 0 or settings.IGNORE_SYMMETRIES:
@@ -214,7 +230,9 @@ class ReservoirMatrix(RGobj):
res[1,0] = res.symmetry * res[0,1].floquetConjugate()
return res
else:
- raise TypeError('Math multiplication is not defined for types %s and %s'%(type(other), ReservoirMatrix))
+ raise TypeError(
+ 'Math multiplication is not defined for types %s and %s'%(
+ type(other), ReservoirMatrix))
def __imatmul__(self, other):
if isinstance(other, ReservoirMatrix):
@@ -244,15 +262,17 @@ class ReservoirMatrix(RGobj):
self.symmetry = 0
return self
else:
- raise TypeError('Math multiplication is not defined for types %s and %s'%(ReservoirMatrix, type(other)))
+ raise TypeError(
+ 'Math multiplication is not defined for types %s and %s'%(
+ ReservoirMatrix, type(other)))
def __mod__(self, other):
- '''
+ """
Transpose multiplication: Given A, B return C such that
C = A B
12 32 13
- '''
+ """
assert isinstance(other, ReservoirMatrix)
assert other.voltage_shifts == 0
res = ReservoirMatrix(self.global_properties, symmetry=0)
@@ -288,24 +308,8 @@ class ReservoirMatrix(RGobj):
return res
def __eq__(self, other):
- return self.global_properties is other.global_properties and np.all(self.data == other.data)
-
- def shift_energies(self, n=1, derivative=None, **kwargs):
- '''
- Apply RGfunction.shift_energies to every entry.
- '''
- raise DeprecationWarning("function ReservoirMatrix.shift_energies should not be used")
- shifted = ReservoirMatrix(self.global_properties, 0)
- if derivative is None:
- for i in range(2):
- for j in range(2):
- shifted[i,j] = self[i,j].shift_energies(n=n, derivative=None, **kwargs)
- else:
- assert isinstance(derivative, ReservoirMatrix)
- for i in range(2):
- for j in range(2):
- shifted[i,j] = self[i,j].shift_energies(n=n, derivative=derivative[i,j], **kwargs)
- return shifted
+ return self.global_properties is other.global_properties \
+ and np.all(self.data == other.data)
def to_numpy_array(self):
array = np.ndarray((2,2,*self[0,0].values.shape), dtype=np.complex128)
@@ -326,22 +330,28 @@ class ReservoirMatrix(RGobj):
def einsum_34_12_43(a:ReservoirMatrix, b:RGobj, c:ReservoirMatrix) -> RGobj:
- '''
- A_34 B_12 C_43
+ """
+ Compute
+ A_34 B_12 C_43
+ with (implicit) summation over indices 3 and 4.
+ B can be either a scalar or a reservoir matrix.
8 multiplications if b is a scalar,
32 multiplications if b is a reservoir matrix without symmetry,
20 multiplications if b is a reservoir matrix with symmetry and xL != xR,
10 multiplications if b is a reservoir matrix with symmetry and xL == xR.
- '''
+ """
assert isinstance(a, ReservoirMatrix)
+ assert isinstance(b, RGobj)
assert isinstance(c, ReservoirMatrix)
assert a.global_properties is b.global_properties is c.global_properties
assert c.voltage_shifts == 0
symmetry = a.symmetry * b.symmetry * c.symmetry
if symmetry == 0 or settings.IGNORE_SYMMETRIES:
if settings.ENFORCE_SYMMETRIC:
- raise RuntimeError("Unsymmetric einsum_34_12_43: %d %d %d"%(a.symmetry, b.symmetry, c.symmetry))
+ raise RuntimeError(
+ 'Unsymmetric einsum_34_12_43: %d %d %d'%(
+ a.symmetry, b.symmetry, c.symmetry))
return a[0,0] @ b @ c[0,0] \
+ a[0,1] @ b @ c[1,0] \
+ a[1,0] @ b @ c[0,1] \
@@ -349,16 +359,25 @@ def einsum_34_12_43(a:ReservoirMatrix, b:RGobj, c:ReservoirMatrix) -> RGobj:
if not isinstance(b, ReservoirMatrix):
res_01_10 = a[0,1] @ b @ c[1,0]
if a.global_properties.symmetric:
- res = 2*a[0,0] @ b @ c[0,0] + res_01_10 + symmetry * res_01_10.floquetConjugate()
+ res = 2*a[0,0] @ b @ c[0,0] \
+ + res_01_10 \
+ + symmetry * res_01_10.floquetConjugate()
else:
- res = a[0,0] @ b @ c[0,0] + a[1,1] @ b @ c[1,1] + res_01_10 + symmetry * res_01_10.floquetConjugate()
+ res = a[0,0] @ b @ c[0,0] \
+ + a[1,1] @ b @ c[1,1] \
+ + res_01_10 \
+ + symmetry * res_01_10.floquetConjugate()
res.symmetry = symmetry
return res
if a.global_properties.symmetric:
# xL = xR = 0.5
res_00_01 = a[0,1] @ b[0,0] @ c[1,0]
- res_00 = 2 * a[0,0] @ b[0,0] @ c[0,0] + res_00_01 + symmetry * res_00_01.floquetConjugate()
- res_01 = 2 * a[0,0] @ b[0,1] @ c[0,0] + a[0,1] @ b[0,1] @ c[1,0] + a[1,0] @ b[0,1] @ c[0,1]
+ res_00 = 2 * a[0,0] @ b[0,0] @ c[0,0] \
+ + res_00_01 \
+ + symmetry * res_00_01.floquetConjugate()
+ res_01 = 2 * a[0,0] @ b[0,1] @ c[0,0] \
+ + a[0,1] @ b[0,1] @ c[1,0] \
+ + a[1,0] @ b[0,1] @ c[0,1]
res = ReservoirMatrix(a.global_properties, symmetry)
res.voltage_shifts = a.voltage_shifts + b.voltage_shifts + c.voltage_shifts
res[0,0] = res_00
@@ -371,9 +390,18 @@ def einsum_34_12_43(a:ReservoirMatrix, b:RGobj, c:ReservoirMatrix) -> RGobj:
# xL != xR
res_00_01 = a[0,1] @ b[0,0] @ c[1,0]
res_11_01 = a[0,1] @ b[1,1] @ c[1,0]
- res_00 = a[0,0] @ b[0,0] @ c[0,0] + a[1,1] @ b[0,0] @ c[1,1] + res_00_01 + symmetry * res_00_01.floquetConjugate()
- res_11 = a[0,0] @ b[1,1] @ c[0,0] + a[1,1] @ b[1,1] @ c[1,1] + res_11_01 + symmetry * res_11_01.floquetConjugate()
- res_01 = a[1,1] @ b[0,1] @ c[1,1] + a[0,0] @ b[0,1] @ c[0,0] + a[0,1] @ b[0,1] @ c[1,0] + a[1,0] @ b[0,1] @ c[0,1]
+ res_00 = a[0,0] @ b[0,0] @ c[0,0] \
+ + a[1,1] @ b[0,0] @ c[1,1] \
+ + res_00_01 \
+ + symmetry * res_00_01.floquetConjugate()
+ res_11 = a[0,0] @ b[1,1] @ c[0,0] \
+ + a[1,1] @ b[1,1] @ c[1,1] \
+ + res_11_01 \
+ + symmetry * res_11_01.floquetConjugate()
+ res_01 = a[1,1] @ b[0,1] @ c[1,1] \
+ + a[0,0] @ b[0,1] @ c[0,0] \
+ + a[0,1] @ b[0,1] @ c[1,0] \
+ + a[1,0] @ b[0,1] @ c[0,1]
res = ReservoirMatrix(a.global_properties, symmetry)
res.voltage_shifts = a.voltage_shifts + b.voltage_shifts + c.voltage_shifts
res[0,0] = res_00
@@ -382,25 +410,36 @@ def einsum_34_12_43(a:ReservoirMatrix, b:RGobj, c:ReservoirMatrix) -> RGobj:
res[1,0] = symmetry * res_01.floquetConjugate()
return res
-def einsum_34_12_43_double(a:ReservoirMatrix, b:ReservoirMatrix, c:ReservoirMatrix, d:ReservoirMatrix) -> (ReservoirMatrix,ReservoirMatrix):
- '''
+
+def einsum_34_12_43_double(
+ a: ReservoirMatrix,
+ b: ReservoirMatrix,
+ c: ReservoirMatrix,
+ d: ReservoirMatrix
+ ) -> (ReservoirMatrix, ReservoirMatrix):
+ """
A_34 B_12 C_43 , A_34 B_12 D_43
48 multiplications if b is a reservoir matrix,
30 multiplications with symmetries if xL != xR,
15 multiplications with symmetries if xL == xR.
- '''
+ """
assert isinstance(a, ReservoirMatrix)
assert isinstance(b, ReservoirMatrix)
assert isinstance(c, ReservoirMatrix)
assert isinstance(d, ReservoirMatrix)
- assert a.global_properties is b.global_properties is c.global_properties is d.global_properties
+ assert a.global_properties \
+ is b.global_properties \
+ is c.global_properties \
+ is d.global_properties
assert c.voltage_shifts == d.voltage_shifts == 0
symmetry_c = a.symmetry * b.symmetry * c.symmetry
symmetry_d = a.symmetry * b.symmetry * d.symmetry
if symmetry_c == 0 or symmetry_d == 0 or settings.IGNORE_SYMMETRIES:
if settings.ENFORCE_SYMMETRIC:
- raise RuntimeError("Unsymmetric einsum_34_12_43_double: %d %d %d %d"%(a.symmetry, b.symmetry, c.symmetry, d.symmetry))
+ raise RuntimeError(
+ 'Unsymmetric einsum_34_12_43_double: %d %d %d %d'%(
+ a.symmetry, b.symmetry, c.symmetry, d.symmetry))
ab_00 = a[0,0] @ b
ab_01 = a[0,1] @ b
ab_10 = a[1,0] @ b
@@ -409,8 +448,6 @@ def einsum_34_12_43_double(a:ReservoirMatrix, b:ReservoirMatrix, c:ReservoirMatr
ab_00 @ c[0,0] + ab_01 @ c[1,0] + ab_10 @ c[0,1] + ab_11 @ c[1,1],
ab_00 @ d[0,0] + ab_01 @ d[1,0] + ab_10 @ d[0,1] + ab_11 @ d[1,1]
)
- if not isinstance(b, ReservoirMatrix):
- raise NotImplementedError
if a.global_properties.symmetric:
# xL == xR == 0.5
ab_00_00 = a[0,0] @ b[0,0]
@@ -418,17 +455,25 @@ def einsum_34_12_43_double(a:ReservoirMatrix, b:ReservoirMatrix, c:ReservoirMatr
ab_01_00 = a[0,1] @ b[0,0]
ab_01_01 = a[0,1] @ b[0,1]
ab_10_01 = a[1,0] @ b[0,1]
+
res_c_00_01 = ab_01_00 @ c[1,0]
- res_c_00 = 2 * ab_00_00 @ c[0,0] + res_c_00_01 + symmetry_c * res_c_00_01.floquetConjugate()
- res_c_01 = 2 * ab_00_01 @ c[0,0] + ab_01_01 @ c[1,0] + ab_10_01 @ c[0,1]
+ res_c_00 = 2 * ab_00_00 @ c[0,0] \
+ + res_c_00_01 \
+ + symmetry_c * res_c_00_01.floquetConjugate()
+ res_c_01 = 2 * ab_00_01 @ c[0,0] \
+ + ab_01_01 @ c[1,0] \
+ + ab_10_01 @ c[0,1]
res_c = ReservoirMatrix(a.global_properties, symmetry_c)
res_c.voltage_shifts = a.voltage_shifts + b.voltage_shifts + c.voltage_shifts
res_c[0,0] = res_c_00
res_c[1,1] = res_c_00.copy()
res_c[0,1] = res_c_01
res_c[1,0] = symmetry_c * res_c_01.floquetConjugate()
+
res_d_00_01 = ab_01_00 @ d[1,0]
- res_d_00 = 2 * ab_00_00 @ d[0,0] + res_d_00_01 + symmetry_d * res_d_00_01.floquetConjugate()
+ res_d_00 = 2 * ab_00_00 @ d[0,0] \
+ + res_d_00_01 \
+ + symmetry_d * res_d_00_01.floquetConjugate()
res_d_01 = 2 * ab_00_01 @ d[0,0] + ab_01_01 @ d[1,0] + ab_10_01 @ d[0,1]
res_d = ReservoirMatrix(a.global_properties, symmetry_d)
res_d.voltage_shifts = a.voltage_shifts + b.voltage_shifts + d.voltage_shifts
@@ -450,22 +495,42 @@ def einsum_34_12_43_double(a:ReservoirMatrix, b:ReservoirMatrix, c:ReservoirMatr
ab_01_11 = a[0,1] @ b[1,1]
ab_01_01 = a[0,1] @ b[0,1]
ab_10_01 = a[1,0] @ b[0,1]
+
res_c_00_01 = ab_01_00 @ c[1,0]
res_c_11_01 = ab_01_11 @ c[1,0]
- res_c_00 = ab_00_00 @ c[0,0] + ab_11_00 @ c[1,1] + res_c_00_01 + symmetry_c * res_c_00_01.floquetConjugate()
- res_c_11 = ab_00_11 @ c[0,0] + ab_11_11 @ c[1,1] + res_c_11_01 + symmetry_c * res_c_11_01.floquetConjugate()
- res_c_01 = ab_11_01 @ c[1,1] + ab_00_01 @ c[0,0] + ab_01_01 @ c[1,0] + ab_10_01 @ c[0,1]
+ res_c_00 = ab_00_00 @ c[0,0] \
+ + ab_11_00 @ c[1,1] \
+ + res_c_00_01 \
+ + symmetry_c * res_c_00_01.floquetConjugate()
+ res_c_11 = ab_00_11 @ c[0,0] \
+ + ab_11_11 @ c[1,1] \
+ + res_c_11_01 \
+ + symmetry_c * res_c_11_01.floquetConjugate()
+ res_c_01 = ab_11_01 @ c[1,1] \
+ + ab_00_01 @ c[0,0] \
+ + ab_01_01 @ c[1,0] \
+ + ab_10_01 @ c[0,1]
res_c = ReservoirMatrix(a.global_properties, symmetry_c)
res_c.voltage_shifts = a.voltage_shifts + b.voltage_shifts + c.voltage_shifts
res_c[0,0] = res_c_00
res_c[1,1] = res_c_11
res_c[0,1] = res_c_01
res_c[1,0] = symmetry_c * res_c_01.floquetConjugate()
+
res_d_00_01 = ab_01_00 @ d[1,0]
res_d_11_01 = ab_01_11 @ d[1,0]
- res_d_00 = ab_00_00 @ d[0,0] + ab_11_00 @ d[1,1] + res_d_00_01 + symmetry_d * res_d_00_01.floquetConjugate()
- res_d_11 = ab_00_11 @ d[0,0] + ab_11_11 @ d[1,1] + res_d_11_01 + symmetry_d * res_d_11_01.floquetConjugate()
- res_d_01 = ab_11_01 @ d[1,1] + ab_00_01 @ d[0,0] + ab_01_01 @ d[1,0] + ab_10_01 @ d[0,1]
+ res_d_00 = ab_00_00 @ d[0,0] \
+ + ab_11_00 @ d[1,1] \
+ + res_d_00_01 \
+ + symmetry_d * res_d_00_01.floquetConjugate()
+ res_d_11 = ab_00_11 @ d[0,0] \
+ + ab_11_11 @ d[1,1] \
+ + res_d_11_01 \
+ + symmetry_d * res_d_11_01.floquetConjugate()
+ res_d_01 = ab_11_01 @ d[1,1] \
+ + ab_00_01 @ d[0,0] \
+ + ab_01_01 @ d[1,0] \
+ + ab_10_01 @ d[0,1]
res_d = ReservoirMatrix(a.global_properties, symmetry_d)
res_d.voltage_shifts = a.voltage_shifts + b.voltage_shifts + d.voltage_shifts
res_d[0,0] = res_d_00
@@ -473,10 +538,13 @@ def einsum_34_12_43_double(a:ReservoirMatrix, b:ReservoirMatrix, c:ReservoirMatr
res_d[0,1] = res_d_01
res_d[1,0] = symmetry_d * res_d_01.floquetConjugate()
return res_c, res_d
- raise NotImplementedError
-def product_combinations(a:ReservoirMatrix, b:ReservoirMatrix) -> (ReservoirMatrix,ReservoirMatrix):
- '''
+
+def product_combinations(
+ a: ReservoirMatrix,
+ b: ReservoirMatrix
+ ) -> (ReservoirMatrix, ReservoirMatrix):
+ """
Equivalent to
lambda a, b: a @ b, a % b
@@ -486,17 +554,19 @@ def product_combinations(a:ReservoirMatrix, b:ReservoirMatrix) -> (ReservoirMatr
12 multiplications (instead of 16) without symmetry,
1 ReservoirMatrix multiplication with symmetry
- -> 4 multiplications with xL == xR,
- -> 7 multiplications with xL != xR,
- -> 8 multiplications without symmetry
- '''
+ Number of multiplications when using symmetries:
+ 4 multiplications with xL == xR,
+ 7 multiplications with xL != xR,
+ 8 multiplications without symmetry
+ """
assert isinstance(a, ReservoirMatrix)
assert isinstance(b, ReservoirMatrix)
assert a.global_properties is b.global_properties
symmetry = a.symmetry * b.symmetry
if symmetry == 0 or settings.IGNORE_SYMMETRIES:
if settings.ENFORCE_SYMMETRIC:
- raise RuntimeError("Unsymmetric product_combinations: %d %d"%(a.symmetry, b.symmetry))
+ raise RuntimeError(
+ 'Unsymmetric product_combinations: %d %d'%(a.symmetry, b.symmetry))
ab_0000 = a[0,0] @ b[0,0]
ab_0110 = a[0,1] @ b[1,0]
ab_1001 = a[1,0] @ b[0,1]
diff --git a/package/src/frtrg_kondo/rtrg.py b/package/src/frtrg_kondo/rtrg.py
index 2c4572b..39f6ce9 100644
--- a/package/src/frtrg_kondo/rtrg.py
+++ b/package/src/frtrg_kondo/rtrg.py
@@ -30,14 +30,14 @@ else:
class GlobalRGproperties:
- '''
+ """
Shared properties container object for RGfunction.
The properties stored in this object must be shared by all RGfunctions
which are used together in the same calculation. Usually all RGfunctions
own a reference to the same GlobalRGproperties object, which makes it
simple to adopt the energies of all these RGfunctions and to check whether
different RGfunctions are compatible.
- '''
+ """
DEFAULTS = dict(
nmax = 0,
padding = 0,
@@ -51,31 +51,31 @@ class GlobalRGproperties:
)
def __init__(self, **kwargs):
- '''
+ """
expected arguments:
energy
omega
nmax
voltage_branches = 0
vdc = 0 or mu = 0
- '''
+ """
settings.logger.debug('Created new GlobalRGproperties object')
self.__dict__.update(kwargs)
def __getattr__(self, name):
- '''
+ """
Return property of take value from DEFAULTS if property is not set.
- '''
+ """
try:
return self.DEFAULTS[name]
except KeyError:
raise AttributeError()
def shape(self):
- '''
+ """
Shape of RGfunction.values for all RGfunctions with these global
properties.
- '''
+ """
floquet_dim = 2*self.__dict__['nmax'] + 1
if self.__dict__['voltage_branches']:
return (2*self.__dict__['voltage_branches'] + 1, floquet_dim, floquet_dim)
@@ -83,10 +83,10 @@ class GlobalRGproperties:
return (floquet_dim, floquet_dim)
def copy(self):
- '''
+ """
Create a copy of self. It is assumed that all attributes of self are
implicitly copied on assignment.
- '''
+ """
return GlobalRGproperties(**self.__dict__)
@@ -100,18 +100,7 @@ class RGobj:
class RGfunction(RGobj):
- # Flags:
- # 1 << 0 : matrix C contiguous
- # 1 << 1 : matrix F contiguous
- INV_COUNTER = [0 for i in range(1 << 2)]
- # Flags:
- # 1 << 0 : symmetric
- # 1 << 1 : matrix 1 C contiguous
- # 1 << 2 : matrix 1 F contiguous
- # 1 << 3 : matrix 2 C contiguous
- # 1 << 4 : matrix 2 F contiguous
- MM_COUNTER = [0 for i in range(1 << 5)]
- '''
+ """
Matrix X_{nm}(E) = X_{n-m}(E+mΩ) with n,m = -nmax...nmax
self.values:
@@ -131,7 +120,7 @@ class RGfunction(RGobj):
self.symmetry:
Valid values are 0, -1, +1. If non-zero, it states the symmetry
X[::-1,::-1,::-1] == symmetry * X.conjugate() for this Floquet matrix.
- '''
+ """
def __init__(self, global_properties, values=None, voltage_shifts=0, symmetry=0, **kwargs):
super().__init__(global_properties, symmetry)
self.energy_shifted_copies = {}
@@ -164,13 +153,13 @@ class RGfunction(RGobj):
self.energy_shifted_copies.clear()
def copy(self):
- '''
+ """
Copy only values, take a reference to global_properties.
- '''
+ """
return RGfunction(self.global_properties, self._values.copy(), self.voltage_shifts, self.symmetry)
def floquetConjugate(self):
- '''
+ """
For a Floquet matrix A(E)_{nm} this returns the C-transform
C A(E)_{nm} C = A(-E*)_{-n,-m}
with the superoperator C defined by
@@ -181,7 +170,7 @@ class RGfunction(RGobj):
This can only be evaluated if the energy of self lies on the
imaginary axis.
- '''
+ """
if self.symmetry == 1:
return self.copy()
elif self.symmetry == -1:
@@ -196,13 +185,13 @@ class RGfunction(RGobj):
raise NotImplementedError()
def __matmul__(self, other):
- '''
+ """
Convolution (or product in Floquet space) of two RG functions.
Other must be of type RGfunction.
Note: This is only approximately associative, as long the function
converges to 0 for |n| of order of nmax.
- '''
+ """
if isinstance(other, SymRGfunction):
return NotImplemented
if not isinstance(other, RGfunction):
@@ -217,8 +206,6 @@ class RGfunction(RGobj):
right = other.shift_energies(self.voltage_shifts)
symmetry = self.symmetry * other.symmetry * (self.voltage_shifts == 0)
matrix = rtrg_c.multiply_extended(right._values.T, self._values.T, self.padding, symmetry*(self._values.ndim==2), self.clear_corners).T
- if settings.logger.level == settings.logging.DEBUG:
- RGfunction.MM_COUNTER[(symmetry != 0) | (self._values.T.flags.c_contiguous << 3) | (self._values.T.flags.f_contiguous << 4) | (right._values.T.flags.c_contiguous << 1) | (right._values.T.flags.f_contiguous << 2)] += 1
res = RGfunction( self.global_properties, matrix, self.voltage_shifts + other.voltage_shifts, symmetry )
return res
@@ -238,19 +225,17 @@ class RGfunction(RGobj):
right = other.shift_energies(self.voltage_shifts)
self.symmetry *= other.symmetry
self._values = rtrg_c.multiply_extended(right._values.T, self._values.T, self.padding, self.symmetry*(self._values.ndim==2), self.clear_corners, OVERWRITE_LEFT).T
- if settings.logger.level == settings.logging.DEBUG:
- RGfunction.MM_COUNTER[(symmetry != 0) | (self._values.T.flags.c_contiguous << 3) | (self._values.T.flags.f_contiguous << 4) | (right._values.T.flags.c_contiguous << 1) | (right._values.T.flags.f_contiguous << 2)] += 1
self.voltage_shifts += other.voltage_shifts
return self
def __add__(self, other):
- '''
+ """
Add other to self. If other is a scalar or a scalar function of energy
represented by an array of values at self.energies, this treats other
as an identity (or diagonal) Floquet matrix.
Other must be a scalar or array of same shape as self.energies or an
RGfunction of the same shape and energies as self.
- '''
+ """
if isinstance(other, RGfunction):
assert self.global_properties is other.global_properties
assert self.voltage_shifts == other.voltage_shifts
@@ -281,9 +266,9 @@ class RGfunction(RGobj):
raise TypeError("unsupported operand types for +: RGfunction and", type(other))
def __neg__(self):
- '''
+ """
Return a copy of self with inverted sign of self._values.
- '''
+ """
return RGfunction(self.global_properties, -self._values, self.voltage_shifts, self.symmetry)
def __iadd__(self, other):
@@ -318,12 +303,12 @@ class RGfunction(RGobj):
return self
def __mul__(self, other):
- '''
+ """
Multiplication by scalar or RGfunction.
If other is an RGfunction, this calculates the matrix product
(equivalent to __matmul__). If other is a scalar, this multiplies
self._values by other.
- '''
+ """
if isinstance(other, RGfunction):
return self @ other
if isinstance(other, Number):
@@ -337,12 +322,12 @@ class RGfunction(RGobj):
return NotImplemented
def __imul__(self, other):
- '''
+ """
In-place multiplication by scalar or RGfunction.
If other is an RGfunction, this calculates the matrix product
(equivalent to __matmul__). If other is a scalar, this multiplies
self._values by other.
- '''
+ """
if isinstance(other, RGfunction):
self @= other
self.energy_shifted_copies.clear()
@@ -360,12 +345,12 @@ class RGfunction(RGobj):
return self
def __rmul__(self, other):
- '''
+ """
Reverse multiplication by scalar or RGfunction.
If other is an RGfunction, this calculates the matrix product
(equivalent to __matmul__). If other is a scalar, this multiplies
self._values by other.
- '''
+ """
if isinstance(other, RGfunction):
return other @ self
if isinstance(other, Number):
@@ -379,9 +364,9 @@ class RGfunction(RGobj):
return NotImplemented
def __truediv__(self, other):
- '''
+ """
Divide self by other, which must be a scalar.
- '''
+ """
if isinstance(other, Number):
if other.imag == 0:
symmetry = self.symmetry
@@ -393,9 +378,9 @@ class RGfunction(RGobj):
return NotImplemented
def __itruediv__(self, other):
- '''
+ """
Divide self in-place by other, which must be a scalar.
- '''
+ """
if isinstance(other, Number):
if other.imag != 0:
if other.real == 0:
@@ -427,7 +412,7 @@ class RGfunction(RGobj):
return ( self.global_properties is other.global_properties ) and self.voltage_shifts == other.voltage_shifts and np.allclose(self._values, other._values) and self.symmetry == other.symmetry
def k2lambda(self, shift_matrix=None, calculate_energy_shifts=False):
- '''
+ """
shift is usually n*zinv*mu.
Assume that self is K_n^m(E) = K_n(E-mΩ).
Then calculate Λ_n^m(E) such that (approximately)
@@ -440,14 +425,19 @@ class RGfunction(RGobj):
Some of the linear systems of equation which we need to solve here are
overdetermined. This means that we can in general only get an
approximate solution because an exact solution does not exist.
- '''
+ """
if shift_matrix is not None:
shift_matrix_vb = shift_matrix._values.shape[0]//2
if self._values.ndim == 3 and calculate_energy_shifts:
- invert_array = np.ndarray((2*self.voltage_branches+5,2*self.nmax+1,2*self.nmax+1), dtype=np.complex128)
+ invert_array = np.ndarray(
+ (2*self.voltage_branches+5,2*self.nmax+1,2*self.nmax+1),
+ dtype=np.complex128)
invert_array[2:-2] = -self._values
for v in range(-self.voltage_branches, self.voltage_branches+1):
- invert_array[(v+2+self.voltage_branches, *np.diag_indices(2*self.nmax+1))] += self.energy + v*self.vdc + self.omega*np.arange(-self.nmax, self.nmax+1)
+ invert_array[(v+2+self.voltage_branches, *np.diag_indices(2*self.nmax+1))] += \
+ self.energy \
+ + v*self.vdc \
+ + self.omega*np.arange(-self.nmax, self.nmax+1)
if shift_matrix is not None:
invert_array[v+2+self.voltage_branches] += shift_matrix[v + shift_matrix_vb]
dim0 = 2*self.voltage_branches+1
@@ -466,8 +456,6 @@ class RGfunction(RGobj):
invert_array[:2] = invert_array[2]
invert_array[:2] += shift_matrix[shift_matrix_vb-2:shift_matrix_vb]
inverted = rtrg_c.invert_extended(invert_array.T, self.padding, round(settings.LAZY_INVERSE_FACTOR*self.padding)).T
- if settings.logger.level == settings.logging.DEBUG:
- RGfunction.INV_COUNTER[invert_array.T.flags.c_contiguous | (invert_array.T.flags.f_contiguous << 1)] += 1
symmetry = -1 if (self.symmetry == -1 and shift_matrix.symmetry == -1) else 0
res = RGfunction(self.global_properties, inverted[2:-2], voltage_shifts=self.voltage_shifts, symmetry=symmetry)
res.energy_shifted_copies[2] = RGfunction(self.global_properties, inverted[4:], self.voltage_shifts, symmetry=0)
@@ -492,9 +480,9 @@ class RGfunction(RGobj):
return invert.inverse()
def reduced(self, shift=0):
- '''
+ """
Remove voltage-shifted copies
- '''
+ """
if self._values.ndim == 2 and shift == 0:
return self
assert self._values.ndim == 3
@@ -510,17 +498,15 @@ class RGfunction(RGobj):
return RGfunction(self.global_properties, self._values[diff:-diff], symmetry=self.symmetry, voltage_shifts=self.voltage_shifts)
def inverse(self):
- '''
+ """
multiplicative inverse
- '''
+ """
assert self.voltage_shifts == 0
if self.padding == 0 or settings.LAZY_INVERSE_FACTOR*self.padding < 0.5:
res = np.linalg.inv(self._values)
else:
try:
res = rtrg_c.invert_extended(self._values.T, self.padding, round(settings.LAZY_INVERSE_FACTOR*self.padding)).T
- if settings.logger.level == settings.logging.DEBUG:
- RGfunction.INV_COUNTER[self._values.T.flags.c_contiguous | (self._values.T.flags.f_contiguous << 1)] += 1
except:
settings.logger.exception("padded inversion failed")
res = np.linalg.inv(self._values)
@@ -528,14 +514,14 @@ class RGfunction(RGobj):
def shift_energies(self, n=0):
# TODO: use symmetries
- '''
+ """
Shift energies by n*self.vdc. The calculated RGfunction is kept in cache.
Assumptions:
* On the last axis of self.energies is linear with values separated by self.vdc
* n is an integer
* derivative is a RGfunction or None
* if the length of the last axis of self.energies is < 2, then derivative must not be None
- '''
+ """
if n==0 or (self.vdc==0 and self.mu is None):
return self
try:
diff --git a/package/src/frtrg_kondo/rtrg_c.c b/package/src/frtrg_kondo/rtrg_c.c
index cb975c9..cac57d6 100644
--- a/package/src/frtrg_kondo/rtrg_c.c
+++ b/package/src/frtrg_kondo/rtrg_c.c
@@ -1,7 +1,7 @@
/*
MIT License
-Copyright (c) 2021 Valentin Bruch
+Copyright (c) 2021 Valentin Bruch <valentin.bruch@rwth-aachen.de>
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
@@ -44,6 +44,8 @@ SOFTWARE.
* is recommended.
* * DEBUG: print debugging information to stderr. This is neither complete
* nor really useful.
+ * * ANALYZE: collect some data for analyzing number of calls of specific
+ * functions and structure of input matrices for optimization.
* The following macros can be redefined to optimize performance, adapt to your
* BLAS and LAPACK installation, and adapt to the concrete mathematical problem.
* * TRIANGULAR_OPTIMIZE_THRESHOLD: Threshold for subdividing multiplication
@@ -110,6 +112,20 @@ static const complex_type zero = 0.;
static const complex_type one = 1.;
+#ifdef ANALYZE
+static int TOTAL_MATRIX_MULTIPLICATIONS = 0;
+#define MATRIX_MULTIPLICATIONS_SIZE 0x40
+static int MATRIX_MULTIPLICATIONS[MATRIX_MULTIPLICATIONS_SIZE];
+#define LEFT_F_CONTIGUOUS 0x1
+#define RIGHT_F_CONTIGUOUS 0x2
+#define LEFT_C_CONTIGUOUS 0x4
+#define RIGHT_C_CONTIGUOUS 0x8
+#define SYMMETRIC 0x10
+#define TWO_DIM 0x20
+#define STRINGIFY(x) #x
+#define TOSTRING(x) STRINGIFY(x)
+#endif
+
#define PY_SSIZE_T_CLEAN
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
@@ -214,7 +230,8 @@ static void extrapolate_top(
row2 = input[(chunk+2)*rows_in+2];
i = chunk + 1;
while(--i > 0)
- output[cutoff-chunk+(chunk-i)*(out_rows+1)] = extrapolate(i, row0, row1, row2);
+ output[cutoff-chunk+(chunk-i)*(out_rows+1)]
+ = extrapolate(i, row0, row1, row2);
}
}
#else /* PARALLEL_EXTRAPOLATION */
@@ -284,7 +301,8 @@ static void extrapolate_bottom(
row2 = input[-(chunk+2)*rows_in-2];
i = chunk + 1;
while(--i > 0)
- output[chunk - cutoff + (cutoff - 1 - chunk + i)*(out_rows+1)] = extrapolate(i, row0, row1, row2);
+ output[chunk - cutoff + (cutoff - 1 - chunk + i)*(out_rows+1)]
+ = extrapolate(i, row0, row1, row2);
}
}
#else /* PARALLEL_EXTRAPOLATION */
@@ -349,7 +367,10 @@ static void extrapolate_left(
jmax = cutoff - i;
j = -1;
while (++j <= jmax)
- output[out_rows * jmax + j] = extrapolate(i, input[i+j], input[i+j+rows_in+1], input[i+j+2*(rows_in+1)]);
+ output[out_rows * jmax + j] = extrapolate(i,
+ input[i + j],
+ input[i + j + rows_in + 1],
+ input[i + j + 2*rows_in + 2]);
}
}
#else /* PARALLEL_EXTRAPOLATION */
@@ -415,7 +436,10 @@ static void extrapolate_right(
jmax = cutoff - i;
j = -1;
while (++j <= jmax)
- output[(out_rows+1) * (i-1) + j] = extrapolate(i, input[j], input[j-rows_in-1], input[j-2*rows_in-2]);
+ output[(out_rows+1) * (i-1) + j] = extrapolate(i,
+ input[j],
+ input[j - rows_in - 1],
+ input[j - 2*rows_in - 2]);
}
}
#else /* PARALLEL_EXTRAPOLATION */
@@ -462,8 +486,8 @@ static void extrapolate_right(
* b00 b01 r20 r21
* b11 b12 r31
*
- * In this representation, the pointers handed to the following functions have
- * the meaning:
+ * In this representation, the pointers handed to the following functions
+ * have the meaning:
* extrapolate_top: output=t00, input=m00
* extrapolate_left: output=l00, input=m00
* extrapolate_bottom: output=b00, input=m45
@@ -776,7 +800,12 @@ static void multiply_UL_inplace(
#endif /* CBLAS */
/* Step 2: overwrite Ad with Ad Bd */
- multiply_UL_inplace(part2, a + (a_dim0+1)*part1, a_dim0, b + (b_dim0+1)*part1, b_dim0);
+ multiply_UL_inplace(
+ part2,
+ a + (a_dim0+1)*part1,
+ a_dim0,
+ b + (b_dim0+1)*part1,
+ b_dim0);
/* Step 3: overwrite Aa with Aa Ba */
multiply_UL_inplace(part1, a, a_dim0, b, b_dim0);
/* Step 4: add Ab Bc to Aa */
@@ -959,7 +988,12 @@ static void multiply_LU_inplace(
/* Step 2: overwrite Aa with Aa Ba */
multiply_LU_inplace(part1, a, a_dim0, b, b_dim0);
/* Step 3: overwrite Ad with Ad Bd */
- multiply_LU_inplace(part2, a + (a_dim0+1)*part1, a_dim0, b + (b_dim0+1)*part1, b_dim0);
+ multiply_LU_inplace(
+ part2,
+ a + (a_dim0+1)*part1,
+ a_dim0,
+ b + (b_dim0+1)*part1,
+ b_dim0);
/* Step 4: add Ac Bb to Ad */
#ifdef CBLAS
gemm(
@@ -1198,9 +1232,13 @@ static PyObject* extend_matrix(PyObject *self, PyObject *args)
return NULL;
if (PyArray_TYPE(input) != NPY_COMPLEX_TYPE)
- return PyErr_Format(PyExc_ValueError, "array of type complex128 required");
+ return PyErr_Format(
+ PyExc_ValueError,
+ "array of type complex128 required");
if (PyArray_NDIM(input) < 2)
- return PyErr_Format(PyExc_ValueError, "1st argument must have at least 2 dimensions.");
+ return PyErr_Format(
+ PyExc_ValueError,
+ "1st argument must have at least 2 dimensions.");
if (cutoff <= 0)
{
@@ -1209,7 +1247,9 @@ static PyObject* extend_matrix(PyObject *self, PyObject *args)
}
if ((PyArray_DIM(input, 0) < cutoff + 3) || (PyArray_DIM(input, 1) < cutoff + 3))
- return PyErr_Format(PyExc_ValueError, "Matrix is too small or cutoff too large.");
+ return PyErr_Format(
+ PyExc_ValueError,
+ "Matrix is too small or cutoff too large.");
PyArrayObject *finput = (PyArrayObject*) PyArray_FromArray(
input,
@@ -1223,7 +1263,10 @@ static PyObject* extend_matrix(PyObject *self, PyObject *args)
if (PyArray_NDIM(finput) == 2)
{
- npy_intp shape[] = {PyArray_DIM(finput, 0) + 2*cutoff, PyArray_DIM(finput, 1) + 2*cutoff};
+ npy_intp shape[] = {
+ PyArray_DIM(finput, 0) + 2*cutoff,
+ PyArray_DIM(finput, 1) + 2*cutoff
+ };
output = (PyArrayObject*) PyArray_ZEROS(2, shape, NPY_COMPLEX_TYPE, 1);
if (output)
extend_matrix_worker(
@@ -1346,7 +1389,11 @@ static PyArrayObject *invert_2d(
output = (PyArrayObject*) PyArray_EMPTY(2, PyArray_DIMS(finput), NPY_COMPLEX_TYPE, 1);
if (output)
for (int i=0; i<size; ++i)
- memcpy( PyArray_GETPTR2(output, 0, i), extended + cutoff + (i+cutoff)*extended_stride, size*sizeof(complex_type) );
+ memcpy(
+ PyArray_GETPTR2(output, 0, i),
+ extended + cutoff + (i+cutoff)*extended_stride,
+ size*sizeof(complex_type)
+ );
}
free( extended );
@@ -1417,7 +1464,11 @@ static PyArrayObject *invert_nd(
PyErr_WarnEx(PyExc_RuntimeWarning, "encountered singular matrix.", 1);
outptr = PyArray_DATA(output) + i*out_matrixstride;
for (j=0; j<size; ++j)
- memcpy( outptr + j*out_colstride, extended + cutoff + (j+cutoff)*extended_stride, size*sizeof(complex_type) );
+ memcpy(
+ outptr + j*out_colstride,
+ extended + cutoff + (j+cutoff)*extended_stride,
+ size*sizeof(complex_type)
+ );
}
free( extended );
@@ -1493,7 +1544,11 @@ static PyArrayObject *invert_nd_parallel(
PyErr_WarnEx(PyExc_RuntimeWarning, "encountered singular matrix.", 1);
outptr = PyArray_DATA(output) + i*out_matrixstride;
for (int j=0; j<size; ++j)
- memcpy( outptr + j*out_colstride, extended + cutoff + (j+cutoff)*extended_stride, size*sizeof(complex_type) );
+ memcpy(
+ outptr + j*out_colstride,
+ extended + cutoff + (j+cutoff)*extended_stride,
+ size*sizeof(complex_type)
+ );
free( extended );
}
if (fatal_error)
@@ -1705,12 +1760,20 @@ static void multiply_extended_worker(
#pragma omp for
#endif
for (int i=0; i<clear_corners; i++)
- memset( output + (ncolr-clear_corners+i)*out_dim0, 0, (i+1)*sizeof(complex_type) );
+ memset(
+ output + (ncolr-clear_corners+i)*out_dim0,
+ 0,
+ (i+1)*sizeof(complex_type)
+ );
#ifdef PARALLEL_EXTRAPOLATION
#pragma omp for
#endif
for (int i=0; i<clear_corners; i++)
- memset( output + i*(out_dim0+1) + nrowl - clear_corners, 0, (clear_corners-i)*sizeof(complex_type) );
+ memset(
+ output + i*(out_dim0+1) + nrowl - clear_corners,
+ 0,
+ (clear_corners-i)*sizeof(complex_type)
+ );
}
if (cutoff <= 0)
@@ -2200,12 +2263,20 @@ static void multiply_symmetric_nonsquare(
#pragma omp for
#endif
for (int i=0; i<clear_corners; i++)
- memset( left + (ncolr-clear_corners+i)*left_dim0, 0, (i+1)*sizeof(complex_type) );
+ memset(
+ left + (ncolr-clear_corners+i)*left_dim0,
+ 0,
+ (i+1)*sizeof(complex_type)
+ );
#ifdef PARALLEL_EXTRAPOLATION
#pragma omp for
#endif
for (int i=0; i<clear_corners; i++)
- memset( left + i*(left_dim0+1) + nrowl - clear_corners, 0, (clear_corners-i)*sizeof(complex_type) );
+ memset(
+ left + i*(left_dim0+1) + nrowl - clear_corners,
+ 0,
+ (clear_corners-i)*sizeof(complex_type)
+ );
}
}
@@ -2298,12 +2369,20 @@ static void multiply_symmetric_square(
#pragma omp for
#endif
for (int i=0; i<clear_corners; i++)
- memset( left + (ncolr-clear_corners+i)*left_dim0, 0, (i+1)*sizeof(complex_type) );
+ memset(
+ left + (ncolr-clear_corners+i)*left_dim0,
+ 0,
+ (i+1)*sizeof(complex_type)
+ );
#ifdef PARALLEL_EXTRAPOLATION
#pragma omp for
#endif
for (int i=0; i<clear_corners; i++)
- memset( left + i*(left_dim0+1) + nrowl - clear_corners, 0, (clear_corners-i)*sizeof(complex_type) );
+ memset(
+ left + i*(left_dim0+1) + nrowl - clear_corners,
+ 0,
+ (clear_corners-i)*sizeof(complex_type)
+ );
}
#ifdef DEBUG
@@ -2499,7 +2578,9 @@ static PyArrayObject* multiply_extended_2d_symmetric(
{
Py_DECREF(fright);
Py_DECREF(fleft);
- PyErr_SetString(PyExc_RuntimeError, "Failed to resize (enlargen) matrix.");
+ PyErr_SetString(
+ PyExc_RuntimeError,
+ "Failed to resize (enlargen) matrix.");
return NULL;
}
}
@@ -2544,7 +2625,9 @@ static PyArrayObject* multiply_extended_2d_symmetric(
if (!PyArray_Resize(fleft, &dims, 1, NPY_FORTRANORDER))
{
Py_DECREF(fleft);
- PyErr_SetString(PyExc_RuntimeError, "Failed to resize (truncate) matrix.");
+ PyErr_SetString(
+ PyExc_RuntimeError,
+ "Failed to resize (truncate) matrix.");
return NULL;
}
}
@@ -2564,6 +2647,8 @@ static PyArrayObject* multiply_extended_nd_symmetric(
const char flags
)
{
+ if (PyErr_WarnEx(NULL, "This function should probably not be used!", 1))
+ return NULL;
#ifdef DEBUG
fprintf(stderr, "Entering multiply_extended_nd_symmetric\n");
#endif
@@ -2640,7 +2725,11 @@ static PyArrayObject* multiply_extended_nd_symmetric(
{
if (ncolr != ncoll)
/* copy data */
- memcpy( left_data + i*left_matrixstride, PyArray_DATA(fleft) + i*PyArray_STRIDE(fleft, 2), PyArray_STRIDE(fleft, 2) );
+ memcpy(
+ left_data + i*left_matrixstride,
+ PyArray_DATA(fleft) + i*PyArray_STRIDE(fleft, 2),
+ PyArray_STRIDE(fleft, 2)
+ );
multiply_symmetric_worker(
nrowl,
ncoll,
@@ -2697,26 +2786,49 @@ static PyObject* multiply_extended(PyObject *self, PyObject *args)
char flags = 0;
/* Parse the arguments. symmetry and flags are optional. */
- if (!PyArg_ParseTuple(args, "O!O!i|iic", &PyArray_Type, &left, &PyArray_Type, &right, &cutoff, &symmetry, &clear_corners, &flags))
+ if (!PyArg_ParseTuple(
+ args,
+ "O!O!i|iic",
+ &PyArray_Type,
+ &left,
+ &PyArray_Type,
+ &right,
+ &cutoff,
+ &symmetry,
+ &clear_corners,
+ &flags))
return NULL;
- if (PyArray_TYPE(left) != NPY_COMPLEX_TYPE || PyArray_TYPE(right) != NPY_COMPLEX_TYPE)
+ if ( PyArray_TYPE(left) != NPY_COMPLEX_TYPE
+ || PyArray_TYPE(right) != NPY_COMPLEX_TYPE)
return PyErr_Format(PyExc_ValueError, "arrays of type complex128 required");
if ( PyArray_NDIM(left) != PyArray_NDIM(right) ||
PyArray_NDIM(left) < 2 ||
PyArray_DIM(left, 1) != PyArray_DIM(right, 0))
- return PyErr_Format(PyExc_ValueError, "Shape of the 2 matrices does not allow multiplication.");
+ return PyErr_Format(
+ PyExc_ValueError,
+ "Shape of the 2 matrices does not allow multiplication.");
const int min_size = cutoff + (clear_corners > 4 ? clear_corners-1 : 3);
if (cutoff > 0 && (
PyArray_DIM(left, 0) < min_size ||
PyArray_DIM(left, 1) < min_size ||
PyArray_DIM(right, 1) < min_size))
- return PyErr_Format(PyExc_ValueError, "Matrices is too small (%d, %d, %d) or cutoff too large (%d).", PyArray_DIM(left, 0), PyArray_DIM(left, 1), PyArray_DIM(right, 1), cutoff);
+ return PyErr_Format(
+ PyExc_ValueError,
+ "Matrices is too small (%d, %d, %d) or cutoff too large (%d).",
+ PyArray_DIM(left, 0),
+ PyArray_DIM(left, 1),
+ PyArray_DIM(right, 1),
+ cutoff);
if (PyArray_NDIM(left) > 2 &&
- memcmp(PyArray_DIMS(left) + 2, PyArray_DIMS(right) + 2, (PyArray_NDIM(left)-2)*sizeof(npy_intp)))
+ memcmp(
+ PyArray_DIMS(left) + 2,
+ PyArray_DIMS(right) + 2,
+ (PyArray_NDIM(left)-2)*sizeof(npy_intp))
+ )
return PyErr_Format(PyExc_ValueError, "dimensions of matrices do not match");
/* Get an F-contiguous version of right. */
@@ -2731,20 +2843,52 @@ static PyObject* multiply_extended(PyObject *self, PyObject *args)
return PyErr_Format(PyExc_RuntimeError, "Failed to create array");
if (symmetry && (
- abs(((int) PyArray_DIM(left, 0)) - ((int) PyArray_DIM(left, 1))) > 1 ||
- abs(((int) PyArray_DIM(left, 0)) - ((int) PyArray_DIM(right, 1))) > 1 ||
- abs(((int) PyArray_DIM(right, 0)) - ((int) PyArray_DIM(right, 1))) > 1 ))
+ abs(((int) PyArray_DIM(left, 0)) - ((int) PyArray_DIM(left, 1))) > 1
+ || abs(((int) PyArray_DIM(left, 0)) - ((int) PyArray_DIM(right, 1))) > 1
+ || abs(((int) PyArray_DIM(right, 0)) - ((int) PyArray_DIM(right, 1))) > 1 ))
symmetry = 0;
+#ifdef ANALYZE
+ int mm_flags = 0;
+ if (PyArray_FLAGS(left) & NPY_ARRAY_F_CONTIGUOUS)
+ mm_flags |= LEFT_F_CONTIGUOUS;
+ else if (PyArray_FLAGS(left) & NPY_ARRAY_C_CONTIGUOUS)
+ mm_flags |= LEFT_C_CONTIGUOUS;
+ if (PyArray_FLAGS(right) & NPY_ARRAY_F_CONTIGUOUS)
+ mm_flags |= RIGHT_F_CONTIGUOUS;
+ else if (PyArray_FLAGS(right) & NPY_ARRAY_C_CONTIGUOUS)
+ mm_flags |= RIGHT_C_CONTIGUOUS;
+ if (symmetry)
+ mm_flags |= SYMMETRIC;
+ if (PyArray_NDIM(left) == 2)
+ mm_flags |= TWO_DIM;
+ TOTAL_MATRIX_MULTIPLICATIONS++;
+ MATRIX_MULTIPLICATIONS[mm_flags]++;
+#endif
+
PyArrayObject *out;
if (symmetry)
{
/* In this case, functions decide dyamically whether fleft needs to be copied.
* Something like fleft will be created in these functions. */
if (PyArray_NDIM(left) == 2)
- out = multiply_extended_2d_symmetric(left, fright, cutoff, symmetry, clear_corners, flags);
+ out = multiply_extended_2d_symmetric(
+ left,
+ fright,
+ cutoff,
+ symmetry,
+ clear_corners,
+ flags
+ );
else
- out = multiply_extended_nd_symmetric(left, fright, cutoff, symmetry, clear_corners, flags);
+ out = multiply_extended_nd_symmetric(
+ left,
+ fright,
+ cutoff,
+ symmetry,
+ clear_corners,
+ flags
+ );
}
else
{
@@ -2779,6 +2923,24 @@ static PyObject* multiply_extended(PyObject *self, PyObject *args)
return (PyObject*) out;
}
+#ifdef ANALYZE
+static PyObject* reset_statistics(PyObject *self, PyObject *args)
+{
+ TOTAL_MATRIX_MULTIPLICATIONS = 0;
+ memset(&MATRIX_MULTIPLICATIONS, 0, MATRIX_MULTIPLICATIONS_SIZE*sizeof(int));
+ Py_RETURN_NONE;
+}
+static PyObject* get_statistics(PyObject *self, PyObject *args)
+{
+ int flags = -1;
+ if (!PyArg_ParseTuple(args, "|i", &flags))
+ return NULL;
+ if (flags >= 0 && flags < MATRIX_MULTIPLICATIONS_SIZE)
+ return PyLong_FromLong(MATRIX_MULTIPLICATIONS[flags]);
+ return PyLong_FromLong(TOTAL_MATRIX_MULTIPLICATIONS);
+}
+#endif
+
/**
* @file
@@ -2839,6 +3001,34 @@ static PyMethodDef FloquetMethods[] =
" inverse, numpy array of shape (n,n,...)"
)
},
+#ifdef ANALYZE
+ {
+ "reset_statistics",
+ reset_statistics,
+ METH_VARARGS,
+ PyDoc_STR(
+ "Reset statistics."
+ )
+ },
+ {
+ "get_statistics",
+ get_statistics,
+ METH_VARARGS,
+ PyDoc_STR(
+ "Get some statistics.\n"
+ "Argument: flags (int)\n"
+ "If no argument or an invalid argument is given,\n"
+ "the total number of matrix products is returned.\n"
+ "These flags can be combined:\n"
+ " left array F contiguous: " TOSTRING(LEFT_F_CONTIGUOUS) "\n"
+ " right array F contiguous: " TOSTRING(RIGHT_F_CONTIGUOUS) "\n"
+ " left array C contiguous: " TOSTRING(LEFT_C_CONTIGUOUS) "\n"
+ " right array C contiguous: " TOSTRING(RIGHT_C_CONTIGUOUS) "\n"
+ " use symmetry in multiply: " TOSTRING(SYMMETRIC) "\n"
+ " arrays are 2 dimensinoal: " TOSTRING(TWO_DIM) "\n"
+ )
+ },
+#endif
{NULL, NULL, 0, NULL}
};
diff --git a/package/src/frtrg_kondo/settings.py b/package/src/frtrg_kondo/settings.py
index c73ba91..c958099 100644
--- a/package/src/frtrg_kondo/settings.py
+++ b/package/src/frtrg_kondo/settings.py
@@ -32,7 +32,7 @@ except:
class GlobalFlags:
- '''
+ """
Define global settings that should be available in all modules.
BASEPATH
@@ -78,7 +78,7 @@ class GlobalFlags:
USE_REFERENCE_IMPLEMENTATION = 0
Use the (slower) reference implementation of RG equation.
Enabling this option also sets IGNORE_SYMMETRIES=1.
- '''
+ """
# Default values of settings. These can be overwritten directly by setting
# environment variables.
@@ -86,8 +86,8 @@ class GlobalFlags:
BASEPATH = os.path.expanduser("~/frtrg_data"),
DB_CONNECTION_STRING = "sqlite:///" + os.path.expanduser("~/frtrg_data/frtrg.sqlite"),
FILENAME = "frtrg-01.h5",
- VERSION = (14, 6, 80, 30318097),
- MIN_VERSION = (14,0),
+ VERSION = (14, 8, 99, 2444441),
+ MIN_VERSION = (14, 0),
LOG_TIME = 10, # in s
ENFORCE_SYMMETRIC = 0,
CHECK_SYMMETRIES = 0,
diff --git a/plot.py b/plot.py
index b5f23e1..505a285 100644
--- a/plot.py
+++ b/plot.py
@@ -494,7 +494,7 @@ def plot_interpolate(
origin='lower')
fig.colorbar(img, ax=[ax1,ax2,ax3,ax4,ax5,ax6])
-def plot_comparison(dm, omega, d=1e-9, **parameters):
+def plot_comparison(dm, omega, d=1e9, **parameters):
"""
Compare current and dc conductance computed from both methods (J and mu) at
fixed frequency as function of Vdc and Vac.
diff --git a/setup.py b/setup.py
index fb3f038..ed9aeb9 100644
--- a/setup.py
+++ b/setup.py
@@ -5,47 +5,52 @@ from setuptools import setup, Extension
import numpy as np
from os import environ
-compiler_args = ['-O3','-Wall','-Wextra','-std=c11']
-linker_args = []
-include_dirs = [np.get_include()]
-name = 'rtrg_c'
-sources = ['rtrg_c.c']
-libraries = ['lapack']
-
-if 'CBLAS' in environ:
- compiler_args += ['-DCBLAS']
- libraries += ['cblas']
- #libraries += ['mkl_rt']
-else:
- libraries += ['blas']
-
-if 'LAPACK_C' in environ:
- compiler_args += ['-DLAPACK_C']
-
-parallel_modifiers = ('PARALLEL', 'PARALLEL_EXTRAPOLATION', 'PARALLEL_EXTRA_DIMS')
-
-need_omp = False
-for modifier in parallel_modifiers:
- if modifier in environ:
- compiler_args += ['-D' + modifier]
- need_omp = True
-
-if need_omp:
- compiler_args += ['-fopenmp']
- linker_args += ['-fopenmp']
-
-if 'DEBUG' in environ:
- compiler_args += ['-DDEBUG']
-
-if 'ANALYZE' in environ:
- compiler_args += ['-DANALYZE']
-
-module = Extension(
- name,
- sources = sources,
- include_dirs = include_dirs,
- libraries = libraries,
- extra_compile_args = compiler_args,
- extra_link_args = linker_args
- )
-setup(ext_modules = [module])
+def main():
+ compiler_args = ['-O3','-Wall','-Wextra','-std=c11']
+ linker_args = []
+ include_dirs = [np.get_include()]
+ name = 'rtrg_c'
+ sources = ['rtrg_c.c']
+ libraries = ['lapack']
+
+ if 'CBLAS' in environ:
+ compiler_args += ['-DCBLAS']
+ libraries += ['cblas']
+ #libraries += ['mkl_rt']
+ else:
+ libraries += ['blas']
+
+ if 'LAPACK_C' in environ:
+ compiler_args += ['-DLAPACK_C']
+
+ parallel_modifiers = ('PARALLEL', 'PARALLEL_EXTRAPOLATION', 'PARALLEL_EXTRA_DIMS')
+
+ need_omp = False
+ for modifier in parallel_modifiers:
+ if modifier in environ:
+ compiler_args += ['-D' + modifier]
+ need_omp = True
+
+ if need_omp:
+ compiler_args += ['-fopenmp']
+ linker_args += ['-fopenmp']
+
+ if 'DEBUG' in environ:
+ compiler_args += ['-DDEBUG']
+
+ if 'ANALYZE' in environ:
+ compiler_args += ['-DANALYZE']
+
+ module = Extension(
+ name,
+ sources = sources,
+ include_dirs = include_dirs,
+ libraries = libraries,
+ extra_compile_args = compiler_args,
+ extra_link_args = linker_args
+ )
+ setup(ext_modules = [module])
+
+
+if __name__ == '__main__':
+ main()
--
GitLab