From af29a2bdaf686689be24f427d8a5920bb57f9f50 Mon Sep 17 00:00:00 2001 From: Kazuyoshi Yoshimi Date: Sat, 18 Jul 2026 14:37:06 +0900 Subject: [PATCH 1/3] Support complex-Hermitian Green's functions in the multi-orbital solver Allow complex off-diagonal G_ab(tau) / spectra (e.g. from spin-orbit coupling), so the semi-positive-definiteness of the full complex-Hermitian spectral matrix is enforced during continuation rather than dropping the imaginary part. - predict_Gtau: keep the dtype of rho_l so the complex off-diagonal G_ab(tau) is no longer truncated to its real part. - write_Gtau / write_rhol / write_rho: also write the imaginary part (0 for real data, backward compatible for readers of the first columns). Requires admmsolver with complex-Hermitian L1 / SPD support. Co-Authored-By: Claude Opus 4.8 (1M context) Claude-Session: https://claude.ai/code/session_01TMeDqXABvd876kMmVY7K34 --- src/spmac/solver/multiple.py | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/spmac/solver/multiple.py b/src/spmac/solver/multiple.py index 7af60a1..d280de3 100644 --- a/src/spmac/solver/multiple.py +++ b/src/spmac/solver/multiple.py @@ -242,9 +242,11 @@ def predict_Gtau(self, rho_l, idx=None): r_l = rho_l.reshape((-1, self.nflavor, self.nflavor)) if self.use_sparse_ir: sampler = sparse_ir.TauSampling(self.basis, self.ts[idx]) - ret = np.zeros((len(idx), self.nflavor, self.nflavor)) + # Keep the dtype of rho_l: complex off-diagonal G_ab(tau) must not be + # truncated to its real part (needed for complex-Hermitian / SOC data). + ret = np.zeros((len(idx), self.nflavor, self.nflavor), dtype=r_l.dtype) for i, j in itertools.product(range(self.nflavor), repeat=2): - ret[:, i, j] = np.real(sampler.evaluate(self.basis.s * r_l[:, i, j])) + ret[:, i, j] = sampler.evaluate(self.basis.s * r_l[:, i, j]) return ret else: return np.einsum("lab, l, li -> iab", r_l, self.s, self.u[:,idx]) @@ -277,7 +279,8 @@ def write_Gtau(self, ts, Gts, outdir: pathlib.Path): f.write(str(ts[i])) for gt in Gts: y = gt[i, ifl, jfl] - f.write(f" {np.real(y)}") + # real then imaginary part (imag is 0 for real data) + f.write(f" {np.real(y)} {np.imag(y)}") f.write("\n") def write_rhol( @@ -304,7 +307,7 @@ def write_rhol( if loglambda is not None: f.write(f"# log_lambda = {loglambda}\n") for i, r in enumerate(r_l[:, ifl, jfl]): - f.write(f"{i} {np.real(r)}\n") + f.write(f"{i} {np.real(r)} {np.imag(r)}\n") def write_rho( self, From dbb519507f6b3f2ee7a9b00ef180e6f4a7b53795 Mon Sep 17 00:00:00 2001 From: Yuichi Motoyama Date: Thu, 6 Aug 2026 17:47:14 +0900 Subject: [PATCH 2/3] Support complex G(tau) input and complex spectrum output Allow optional column_imag when loading G(tau), and write Re/Im for single- and multi-orbital spectrum/Gtau outputs so complex data is preserved end-to-end. Co-authored-by: Cursor --- README.md | 5 ++++- src/spmac/main.py | 30 ++++++++++++++++++++++++------ src/spmac/solver/multiple.py | 9 ++++++--- src/spmac/solver/single.py | 15 +++++++++------ 4 files changed, 43 insertions(+), 16 deletions(-) diff --git a/README.md b/README.md index 2585cd0..d2d333a 100644 --- a/README.md +++ b/README.md @@ -25,7 +25,10 @@ The names of parameters are case-insensitive. - Filename storing G(τ) - When `num_flavor > 1`, "{filein_g}.{a}_{b}" are read for G_{ab}, where a,b = 0,1,...,num_flavor-1 - column: int - - Index of column storing G(τ) (0-origin) + - Index of column storing Re G(τ) (0-origin) +- column_imag: int (optional) + - Index of column storing Im G(τ) (0-origin) + - If omitted, Im G(τ) is treated as 0 (real input, backward compatible) - beta: float - Inverse temperature β - max_omega: float diff --git a/src/spmac/main.py b/src/spmac/main.py index 69d7c23..eb853b5 100644 --- a/src/spmac/main.py +++ b/src/spmac/main.py @@ -49,9 +49,24 @@ def optimize( return loglambdas[i] +def _load_gtau_from_file( + filename: str, column: int, column_imag: Optional[int] = None +) -> NDArray: + """Load G(τ) from a text file. + + If ``column_imag`` is given, return complex G = Re + i Im. + Otherwise return the real column only (imaginary part is treated as 0). + """ + data = np.loadtxt(filename) + g = data[:, column] + if column_imag is not None: + g = g + 1j * data[:, column_imag] + return g + + def run( - params: Dict[str, Any], Gtau: NDArray[np.float64] -) -> Tuple[SolverBase, NDArray[np.float64], float]: + params: Dict[str, Any], Gtau: NDArray +) -> Tuple[SolverBase, NDArray, float]: params = dict_with_lowerkey(params) verbose = params.get("verbose", True) n_trials = params.get("num_trials", 10) @@ -103,18 +118,21 @@ def main(): nflavor = params.get("num_flavor", 1) column: int = params["column"] + column_imag: Optional[int] = params.get("column_imag") + print(f"column: {column}, column_imag: {column_imag}") if nflavor == 1: - Gtau = np.loadtxt(params["filein_g"])[:, column] + Gtau = _load_gtau_from_file(params["filein_g"], column, column_imag) else: - Gtau = np.zeros((1, nflavor, nflavor)) + dtype = np.complex128 if column_imag is not None else np.float64 + Gtau = np.zeros((1, nflavor, nflavor), dtype=dtype) ntau = -1 for i in range(nflavor): for j in range(nflavor): filename = params["filein_g"] + f".{i}_{j}" - Gt = np.loadtxt(filename)[:, column] + Gt = _load_gtau_from_file(filename, column, column_imag) if ntau < 0: ntau = len(Gt) - Gtau = np.zeros((ntau, nflavor, nflavor)) + Gtau = np.zeros((ntau, nflavor, nflavor), dtype=dtype) Gtau[:, i, j] = Gt[:] solver, rho_l, best_loglambda = run(params, Gtau) diff --git a/src/spmac/solver/multiple.py b/src/spmac/solver/multiple.py index d280de3..4c90e1f 100644 --- a/src/spmac/solver/multiple.py +++ b/src/spmac/solver/multiple.py @@ -229,8 +229,10 @@ def callback(): print(" Not converged") return opt - def predict_rho(self, rho_l) -> NDArray[np.float64]: - r_l = rho_l.reshape((-1, self.nflavor, self.nflavor)) + def predict_rho(self, rho_l) -> NDArray: + r_l = np.asarray(rho_l).reshape((-1, self.nflavor, self.nflavor)) + # Keep the dtype of rho_l: complex off-diagonal rho_ab(omega) must not be + # truncated to its real part (needed for complex-Hermitian / SOC data). if self.use_sparse_ir: return np.einsum("lab,lw->wab", r_l, self.basis.v(self.ws)) else: @@ -338,4 +340,5 @@ def write_rho( if loglambda is not None: f.write(f"# log_lambda = {loglambda}\n") for w, r in zip(ws, rs[:, ifl, jfl]): - f.write(f"{w} {np.real(r)/dw}\n") + # real then imaginary part (imag is 0 for real data) + f.write(f"{w} {np.real(r)/dw} {np.imag(r)/dw}\n") diff --git a/src/spmac/solver/single.py b/src/spmac/solver/single.py index 122f908..51ede83 100644 --- a/src/spmac/solver/single.py +++ b/src/spmac/solver/single.py @@ -232,17 +232,19 @@ def callback(): print(" Not converged") return opt - def predict_rho(self, rho_l: NDArray[np.float64]) -> NDArray[np.float64]: + def predict_rho(self, rho_l) -> NDArray: + # Keep the dtype of rho_l so complex spectra are not truncated to real. if self.use_sparse_ir: return rho_l @ self.basis.v(self.ws) else: return rho_l @ self.v - def predict_Gtau(self, rho_l: NDArray[np.float64], idx=None) -> NDArray[np.float64]: + def predict_Gtau(self, rho_l, idx=None) -> NDArray: if idx is None: idx = np.arange(self.ntau) if self.use_sparse_ir: sampler = sparse_ir.TauSampling(self.basis, self.ts[idx]) + # Keep the dtype of rho_l so complex G(tau) is not truncated to real. return sampler.evaluate(self.basis.s * rho_l) else: return rho_l @ self.u[:, idx] @@ -268,7 +270,8 @@ def write_Gtau(self, ts, Gts, outdir: pathlib.Path): f.write(str(ts[i])) for gt in Gts: y = gt[i] - f.write(f" {np.real(y)}") + # real then imaginary part (imag is 0 for real data) + f.write(f" {np.real(y)} {np.imag(y)}") f.write("\n") def write_rhol( @@ -292,7 +295,7 @@ def write_rhol( if loglambda is not None: f.write(f"# log_lambda = {loglambda}\n") for i, r in enumerate(rs): - f.write(f"{i} {np.real(r)}\n") + f.write(f"{i} {np.real(r)} {np.imag(r)}\n") def write_rho( self, @@ -321,5 +324,5 @@ def write_rho( if loglambda is not None: f.write(f"# log_lambda = {loglambda}\n") for w, r in zip(ws, rs): - f.write(f"{w} {np.real(r)/dw}\n") - # f.write(f"{w} {np.real(r)}\n") + # real then imaginary part (imag is 0 for real data) + f.write(f"{w} {np.real(r)/dw} {np.imag(r)/dw}\n") From 87f3ae1b3254d440b10d550e459a4b90653929f6 Mon Sep 17 00:00:00 2001 From: Yuichi Motoyama Date: Thu, 6 Aug 2026 21:53:54 +0900 Subject: [PATCH 3/3] use admmsolver with complex-hermition-support --- pyproject.toml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pyproject.toml b/pyproject.toml index 7bd815c..8eaffde 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -44,3 +44,6 @@ build-backend = "hatchling.build" [tool.mypy] files = "src" + +[tool.uv.sources] +admmsolver = { git = "https://github.com/SpM-lab/admmsolver", branch = "complex-hermitian-support" }