From 086ef98a326351993d4a6085323a69a9c0aff6c3 Mon Sep 17 00:00:00 2001
From: Pauli Virtanen <pauli.t.virtanen@jyu.fi>
Date: Mon, 29 Apr 2024 16:13:06 +0300
Subject: [PATCH] Have interpolators only in Result(), not in Solver()
---
examples/cpr_sns.py | 2 +-
usadelndsoc/solver.py | 221 ++++++++++++++++++++++++++++++++----------
2 files changed, 172 insertions(+), 51 deletions(-)
diff --git a/examples/cpr_sns.py b/examples/cpr_sns.py
index f985c55..0479b7c 100644
--- a/examples/cpr_sns.py
+++ b/examples/cpr_sns.py
@@ -66,7 +66,7 @@ def j(T, h, phi, n=15, eta=0.1, alpha_soc=0.1, L=10, W=2):
res = sol.solve_many(omega=w)
t3 = np.diag([1, 1, -1, -1])
- J = tr(res.J @ t3)
+ J = tr(res.J.A @ t3)
J = -1j * pi * (J * a[:, None, None, None]).sum(axis=0)
Jx = (J[:, :, 0] + J[:, :, 2]).real / 2
diff --git a/usadelndsoc/solver.py b/usadelndsoc/solver.py
index a667c3d..3bafa3a 100644
--- a/usadelndsoc/solver.py
+++ b/usadelndsoc/solver.py
@@ -119,20 +119,11 @@ def _core_property(name, doc=None):
return property(fget, fset, doc=doc)
-def _array_property(name, doc=None, cls=None):
+def _array_property(name, doc=None):
def fget(self):
- v = getattr(self, name)
- if cls is not None:
- return cls(self, v)
- else:
- return v
+ return getattr(self, name)
fget.__name__ = name
- if doc is None:
- doc = ""
- doc = textwrap.dedent(doc)
- if cls is not None:
- doc += f"\n\n:Type: :any:`{cls.__name__}`"
return property(fget, doc=doc)
@@ -234,7 +225,7 @@ class Current:
--------
interpolate_J
"""
- return interpolate_S(
+ return interpolate_J(
self.__parent.x, self.__parent.y, self.__data, x, y, method
)
@@ -279,7 +270,6 @@ class Solver:
Phi = _array_property(
"_Phi",
doc=r"""Current saddle-point solution in Riccati parametrization. Array of complex, shape (nx, ny, 2, 2, 2).""",
- cls=Density,
)
J_tot = _array_property(
"_J_tot",
@@ -287,14 +277,12 @@ class Solver:
Matsubara-summed matrix current between cells, :math:`J = -i \pi T\sum_{\omega_n} \mathcal{J}(\omega_n)`.
Array of complex, shape (nx, ny, 4(direction), 4, 4).
""",
- cls=Current,
)
S_tot = _array_property(
"_S_tot",
doc=r"""Matsubara-summed density, :math:`S = \pi T \sum_{\omega_n} [Q(\omega_n) - \tau_3]`.
Array of complex, shape (nx, ny, 4, 4).
""",
- cls=Density,
)
def reset(self):
@@ -694,7 +682,9 @@ class Solver:
Phi00,
**kw,
)
- return Result(self, omega=omega, Phi=self._Phi)
+ return Result(
+ self, omega=omega, Phi=self._Phi, J_tot=self._J_tot, S_tot=self._S_tot
+ )
def _get_preconditioner(self, J, preconditioner=None):
if preconditioner is None:
@@ -754,7 +744,7 @@ class Solver:
omega0 = np.asarray(omega)
omega = omega0.ravel()
- Phi = np.zeros(omega.shape + self.Phi.A.shape, dtype=self.Phi.A.dtype)
+ Phi = np.zeros(omega.shape + self._Phi.shape, dtype=self._Phi.dtype)
js = np.lexsort((-abs(omega), np.sign(omega.real)))
@@ -778,7 +768,15 @@ class Solver:
Phi = Phi.reshape(omega0.shape + Phi.shape[1:])
- return Result(self, omega=omega0, Phi=Phi)
+ return Result(self, omega=omega0, Phi=Phi, J_tot=self._J_tot, S_tot=self._S_tot)
+
+ def result(self):
+ """
+ Return current values as :any:`Result`.
+ """
+ return Result(
+ self, omega=self.omega, Phi=self._Phi, J_tot=self._J_tot, S_tot=self._S_tot
+ )
def self_consistency(self, T, T_c0, continuation=False, **kw):
r"""
@@ -848,15 +846,15 @@ class Result:
"""
- def __init__(self, parent, omega, Phi, core=None):
+ def __init__(self, parent, omega, Phi, J_tot=None, S_tot=None, core=None):
self._core = core if core is not None else parent._core
self.shape = parent.shape
- self.Lx = parent.Lx
- self.Ly = parent.Ly
self.omega = omega
- self.Omega = parent.Omega
- self.Phi = Phi
+ self._Omega = parent.Omega
+ self._Phi = Phi
self.mask = parent.mask
+ self._J_tot = J_tot
+ self._S_tot = S_tot
def _density(self, v):
if np.asarray(self.omega).ndim == 0:
@@ -872,6 +870,21 @@ class Result:
r = np.moveaxis(v, 0, 3)
return Current(self, r)
+ Lx = _core_property("Lx")
+ Ly = _core_property("Ly")
+ x = _core_property("x")
+ y = _core_property("y")
+
+ @property
+ def Phi(self):
+ """
+ Solution Riccati parameters.
+
+ :Type: :any:`Density`
+ :Shape: (nx, ny, 2, 2, 2)
+ """
+ return Density(self, self._Phi)
+
@property
def Delta(self):
"""
@@ -880,7 +893,17 @@ class Result:
:Type: :any:`Density`
:Shape: (nx, ny, 2, 2)
"""
- return Density(self, _DeltaArray(self.Omega).A)
+ return Density(self, _DeltaArray(self._Omega).A)
+
+ @property
+ def Omega(self):
+ """
+ Local energy matrix.
+
+ :Type: :any:`Density`
+ :Shape: (nx, ny, 4, 4)
+ """
+ return Density(self, self._Omega)
@property
def G(self):
@@ -890,8 +913,8 @@ class Result:
:Type: :any:`Density`
:Shape: (nx, ny, [nw,] 2, 2)
"""
- g = self.Phi[..., 0, :, :]
- gt = self.Phi[..., 1, :, :]
+ g = self._Phi[..., 0, :, :]
+ gt = self._Phi[..., 1, :, :]
I = np.eye(2)
s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
@@ -906,8 +929,8 @@ class Result:
:Type: :any:`Density`
:Shape: (nx, ny, [nw,] 2, 2)
"""
- g = self.Phi[..., 0, :, :]
- gt = self.Phi[..., 1, :, :]
+ g = self._Phi[..., 0, :, :]
+ gt = self._Phi[..., 1, :, :]
I = np.eye(2)
s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
@@ -922,8 +945,8 @@ class Result:
:Type: :any:`Density`
:Shape: (nx, ny, [nw,] 2, 2)
"""
- g = self.Phi[..., 0, :, :]
- gt = self.Phi[..., 1, :, :]
+ g = self._Phi[..., 0, :, :]
+ gt = self._Phi[..., 1, :, :]
I = np.eye(2)
s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
@@ -938,8 +961,8 @@ class Result:
:Type: :any:`Density`
:Shape: (nx, ny, [nw,] 2, 2)
"""
- g = self.Phi[..., 0, :, :]
- gt = self.Phi[..., 1, :, :]
+ g = self._Phi[..., 0, :, :]
+ gt = self._Phi[..., 1, :, :]
I = np.eye(2)
s = _bcast_left(np.sign(self.omega), g.shape)
with np.errstate(invalid="ignore"):
@@ -991,10 +1014,11 @@ class Result:
:Shape: (nx, ny, [nw,] 4, 4)
"""
if np.asarray(self.omega).ndim == 0:
- r = self._get_S(self.omega, self.Phi)
+ r = self._get_S(self.omega, self._Phi)
else:
r = np.array(
- [self._get_S(w, P) for w, P in zip(self.omega, self.Phi)], dtype=complex
+ [self._get_S(w, P) for w, P in zip(self.omega, self._Phi)],
+ dtype=complex,
)
return self._density(r)
@@ -1014,13 +1038,52 @@ class Result:
:Shape: (nx, ny, 4(direction), [nw,] 4, 4)
"""
if np.asarray(self.omega).ndim == 0:
- r = self._get_J(self.omega, self.Phi)
+ r = self._get_J(self.omega, self._Phi)
else:
r = np.array(
- [self._get_J(w, P) for w, P in zip(self.omega, self.Phi)], dtype=complex
+ [self._get_J(w, P) for w, P in zip(self.omega, self._Phi)],
+ dtype=complex,
)
return self._current(r)
+ def _S_s(self, S, msummed=False):
+ Sx = tr(np.kron(S_z, S_x) @ S).imag
+ Sy = tr(np.kron(S_z, S_y) @ S).imag
+ Sz = tr(np.kron(S_z, S_z) @ S).imag
+ ss = np.array([Sx, Sy, Sz])
+ if msummed or np.asarray(self.omega).ndim == 0:
+ r = np.array([jx, jy, jz]).transpose(1, 2, 0)
+ else:
+ r = np.array([jx, jy, jz]).transpose(1, 2, 3, 0)
+ return Density(self, r)
+
+ def _J_c(self, J):
+ t3 = np.diag([1, 1, -1, -1])
+ r = tr(J @ t3).real
+ return Current(self, r)
+
+ def _J_s(self, J, msummed=False):
+ jx = tr(J @ np.kron(S_0, S_x)).real
+ jy = tr(J @ np.kron(S_0, S_y)).real
+ jz = tr(J @ np.kron(S_0, S_z)).real
+ if msummed or np.asarray(self.omega).ndim == 0:
+ r = np.array([jx, jy, jz]).transpose(1, 2, 3, 0)
+ else:
+ r = np.array([jx, jy, jz]).transpose(1, 2, 3, 4, 0)
+ return Current(self, r)
+
+ @property
+ def S_s(self):
+ """
+ Spectral spin density.
+
+ Direction index has same meaning as for :any:`J`.
+
+ :Type: :any:`Current`
+ :Shape: (nx, ny, [nw,] 3(spin-direction))
+ """
+ return self._S_s(self.S.A)
+
@property
def J_c(self):
"""
@@ -1029,11 +1092,9 @@ class Result:
Direction index has same meaning as for :any:`J`.
:Type: :any:`Current`
- :Shape: ([nw,] nx, ny, 4(direction))
+ :Shape: (nx, ny, [nw,] 4(direction))
"""
- t3 = np.diag([1, 1, -1, -1])
- r = tr(self.J.A @ t3)
- return Current(self, r)
+ return self._J_c(self.J.A)
@property
def J_s(self):
@@ -1047,14 +1108,74 @@ class Result:
:Type: :any:`Current`
:Shape: (nx, ny, 4(direction), [nw,] 3(spin-direction))
"""
- jx = tr(self.J.A @ np.kron(s0, sx))
- jy = tr(self.J.A @ np.kron(s0, sy))
- jz = tr(self.J.A @ np.kron(s0, sz))
- if np.asarray(self.omega).ndim == 0:
- r = np.array([jx, jy, jz]).transpose(1, 2, 3, 0)
- else:
- r = np.array([jx, jy, jz]).transpose(1, 2, 3, 4, 0)
- return Current(self, r)
+ return self._J_s(self.J.A)
+
+ @property
+ def J_tot(self):
+ r"""
+ Matsubara-summed matrix current between cells, :math:`J = -i \pi T\sum_{\omega_n} \mathcal{J}(\omega_n)`.
+
+ :Type: :any:`Current`
+ :Shape: (nx, ny, 4(direction), 4, 4)
+ """
+ if self._J_tot is None:
+ raise ValueError("J_tot not computed")
+ return Current(self, self._J_tot)
+
+ @property
+ def J_tot_c(self):
+ """
+ Charge current.
+
+ Direction index has same meaning as for :any:`J`.
+
+ :Type: :any:`Current`
+ :Shape: (nx, ny, 4(direction))
+ """
+ if self._J_tot is None:
+ raise ValueError("J_tot not computed")
+ return self._J_c(self._J_tot)
+
+ @property
+ def J_tot_s(self):
+ """
+ Spin current.
+
+ The current direction index has same meaning as for :any:`J`.
+
+ The spin direction index means (x, y, z).
+
+ :Type: :any:`Current`
+ :Shape: (nx, ny, 4(direction), 3(spin-direction))
+ """
+ if self._J_tot is None:
+ raise ValueError("J_tot not computed")
+ return self._J_s(self._J_tot, msummed=True)
+
+ @property
+ def S_tot(self):
+ """Matsubara-summed density, :math:`S = \pi T \sum_{\omega_n} [Q(\omega_n) - \tau_3]`.
+
+ :Type: :any:`Density`
+ :Shape: (nx, ny, 4, 4)
+ """
+ if self._S_tot is None:
+ raise ValueError("S_tot not computed")
+ return Density(self, self._S_tot)
+
+ @property
+ def S_s(self):
+ """
+ Spin density.
+
+ Direction index has same meaning as for :any:`J`.
+
+ :Type: :any:`Current`
+ :Shape: (nx, ny, 3(spin-direction))
+ """
+ if self._S_tot is None:
+ raise ValueError("S_tot not computed")
+ return self._S_s(self._S_tot, msummed=True)
class UmfpackILU(LinearOperator):
@@ -1674,10 +1795,10 @@ def cpr(
state.dss.append(float(np.squeeze(ds_cur)))
state.Deltas.append(Delta.copy())
state.phis.append(float(np.squeeze(phi)))
- state.Js.append(solver.J_tot.A.copy())
+ state.Js.append(solver.J_tot.copy())
t3 = np.diag([1, 1, -1, -1])
- J_avg = np.mean(tr(solver.J_tot.A[:, :, 0].real @ t3))
+ J_avg = np.mean(tr(solver.J_tot[:, :, 0].real @ t3))
_log.info(f"CPR #{j}: phi = {phi}, J = {J_avg}")
if filename is not None:
--
GitLab