README.md

@@ -34,7 +34,8 @@ installed package directory, or opened with:
 ## Bibliography
 
 P. Virtanen, "Numerical saddle point finding for diffusive superconductor actions",
-https://arxiv.org/abs/2406.19894 (2024).
+[Phys. Rev. B 111, 024510 (2025)](https://doi.org/10.1103/PhysRevB.111.024510),
+[arXiv:2406.19894](https://arxiv.org/abs/2406.19894).
 
 ## Acknowledgments
 

pyproject.toml

@@ -44,4 +44,4 @@ filterwarnings = [
   'ignore:.*Broyden reset',
 ]
 addopts = ['--doctest-modules', '--log-level=ERROR']
-testpaths = ['usadelndsoc', 'tests']
+testpaths = ['src/usadelndsoc', 'tests']

src/tools/version_utils.py

+#!/usr/bin/env python3
+import os
+import subprocess
+import argparse
+
+
+def get_version_info(source_root, version):
+    if source_root and os.path.exists(os.path.join(source_root, ".git")):
+        rev, count = git_version(version, source_root)
+    else:
+        rev = count = None
+
+    is_released = not version.endswith(".dev0")
+    if not is_released and rev is not None:
+        version += f"+{count}.{rev}"
+
+    return version, rev, count, is_released
+
+
+def print_version(source_root, version):
+    cnt = """\
+# THIS FILE IS GENERATED DURING THE BUILD
+
+short_version = '%(version)s'
+version = '%(version)s'
+full_version = %(full_version)r
+git_revision = %(git_revision)r
+commit_count = %(commit_count)r
+release = %(isrelease)s
+
+if not release:
+    version = full_version
+"""
+    full_version, rev, count, is_released = get_version_info(source_root, version)
+
+    print(
+        cnt
+        % {
+            "version": version,
+            "full_version": full_version,
+            "git_revision": rev,
+            "commit_count": count,
+            "isrelease": is_released,
+        }
+    )
+
+
+# Return the git revision as a string
+def git_version(version, cwd):
+    def _minimal_ext_cmd(cmd):
+        # construct minimal environment
+        env = {}
+        for k in ["SYSTEMROOT", "PATH"]:
+            v = os.environ.get(k)
+            if v is not None:
+                env[k] = v
+        env["LANGUAGE"] = "C"
+        env["LANG"] = "C"
+        env["LC_ALL"] = "C"
+        out = subprocess.run(
+            cmd,
+            capture_output=True,
+            check=True,
+            encoding="utf-8",
+            errors="strict",
+            env=env,
+        )
+        return out.stdout.strip()
+
+    try:
+        git_dir = os.path.join(cwd, ".git")
+        out = _minimal_ext_cmd(["git", "--git-dir", git_dir, "rev-parse", "HEAD"])
+        rev = out[:7]
+
+        out = _minimal_ext_cmd(
+            [
+                "git",
+                "--git-dir",
+                git_dir,
+                "rev-list",
+                "HEAD",
+                "--count",
+            ]
+        )
+        count = 0 if not out else int(out)
+    except (OSError, UnicodeError):
+        rev = count = None
+
+    return rev, count
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "--source-root", type=str, default=None, help="Project source root"
+    )
+    parser.add_argument(
+        "--version", type=str, default="0.0.0.dev0", help="Project version base part"
+    )
+    args = parser.parse_args()
+
+    print_version(args.source_root, args.version)

src/usadelndsoc/__init__.py

@@ -3,6 +3,11 @@
 import logging
 import contextlib
 
+try:
+    from .version import version as __version__
+except ImportError:
+    __version__ = "Unknown"
+
 
 logger = logging.getLogger(__name__)
 """

src/usadelndsoc/meson.build

@@ -22,3 +22,16 @@ py3.install_sources(
   subdir : 'usadelndsoc',
   pure : false,
 )
+
+version_utils = meson.current_source_dir() / '../tools/version_utils.py'
+version_py = custom_target(
+  'generate-version',
+  build_always_stale: true,
+  build_by_default: true,
+  output: 'version.py',
+  input: version_utils,
+  capture: true,
+  command: [py3, version_utils, '--version', meson.project_version(), '--source-root', meson.project_source_root()],
+  install: true,
+  install_dir: py3.get_install_dir(subdir : 'usadelndsoc', pure: false)
+)

tests/test_solver.py

@@ -56,6 +56,23 @@ def test_example_sns_1d():
     assert_allclose(rF.squeeze(), F, rtol=0, atol=5e-3 * abs(F).max())
 
 
+def check_symmetry(res):
+    """Check symmetry of Q. This holds only when the matsubara frequency
+    is real, otherwise the symmetry is not local in frequency.
+
+    """
+    if not np.allclose(res.omega.imag, 0, atol=1e-13, rtol=0):
+        return
+
+    Phi = res._Phi
+    g = Phi[..., 0, :, :]
+    gt = Phi[..., 1, :, :]
+
+    s2 = np.array([[0, -1j], [1j, 0]])
+    gtr = s2 @ g.conjugate() @ s2
+    assert_allclose(gt, gtr, atol=1e-10, rtol=1e-5)
+
+
 @pytest.mark.skipif(usadel1 is None, reason="usadel1 not installed")
 def test_example_sns_1d_J_usadel1():
     Lx = 10
@@ -67,6 +84,8 @@ def test_example_sns_1d_J_usadel1():
     success, _ = sol.self_consistency(T=T, T_c0=0, tol=1e-10)
     I0 = sol.result.J_tot.A
 
+    check_symmetry(sol.result)
+
     tau3 = np.diag([1, 1, -1, -1])
     I0 = tr(I0[:-1, 0, 0] @ tau3)
 
@@ -108,6 +127,8 @@ def test_example_sss_1d_J_usadel1():
     success, _ = sol.self_consistency(T=T, T_c0=T_c0, workers=1)
     I0 = sol.result.J_tot.A
 
+    check_symmetry(sol.result)
+
     tau3 = np.diag([1, 1, -1, -1])
     I0 = tr(I0[:-1, 0, 0] @ tau3)
 
@@ -159,6 +180,8 @@ def test_example_sns_usadel1(omega):
     sol.solve(omega=omega)
     res = sol.result
 
+    check_symmetry(res)
+
     g = usadel1.Geometry(1, 2)
     g.t_type = [usadel1.NODE_CLEAN_S_TERMINAL, usadel1.NODE_CLEAN_S_TERMINAL]
     g.w_type = usadel1.WIRE_TYPE_N
@@ -331,6 +354,8 @@ def test_soc_analytic(T, n, alpha_soc):
     for ww, aa in zip(w, a):
         sol0.solve(omega=ww)
         res0 = sol0.result
+        check_symmetry(res0)
+
         sol.omega = ww
         J = sol._core.grad_A(res0.Phi.A.copy()).transpose(0, 1, 2, 4, 3) * (1 / 16)
         Jtot += -1j * pi * aa * J
@@ -346,7 +371,7 @@ def test_soc_analytic(T, n, alpha_soc):
 
         Jy = -1j * tr(J @ t3)[n // 2, n // 2, 1]
 
-        assert_allclose(Jy, Jy_an, rtol=1e-2)
+        assert_allclose(Jy, Jy_an, rtol=1e-2, atol=1e-8)
 
         Fs.append(res0.F.A.copy())
         Gs.append(res0.G.A.copy())
@@ -376,6 +401,8 @@ def test_selfcons_h():
         Delta = solver.result.Delta.A
         ac = solver.result.ac.real
 
+        check_symmetry(solver.result)
+
         if acval is None:
             acval = ac
 
@@ -413,6 +440,7 @@ def test_selfcons_Delta_min(h):
     for Delta in Deltas:
         solver.Delta[...] = Delta * np.eye(2)
         solver.matsubara_sums(T, T_c0, workers=1)
+        check_symmetry(solver.result)
         acs.append(float(solver.result.ac_tot.real))
         print(Delta, acs[-1])
         assert acs[-1] >= acs[0], (Delta, acs)