diff --git a/python/triqs_dftkit/wien2k/symqmc.py b/python/triqs_dftkit/wien2k/symqmc.py
new file mode 100644
index 0000000..bcc85ca
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/symqmc.py
@@ -0,0 +1,294 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see <http://www.gnu.org/licenses/>.
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj correlated-shell symmetry file
+(case.symqmc), replacing the corresponding output of the dmftproj Fortran
+executable.
+
+Given the dmftproj symmetry input (case.dmftsym), the projector definition
+(case.indmftpr) and the structure (case.struct), this builds the spinor symmetry
+matrices and writes them in the case.symqmc format the converter reads. It
+reproduces the Fortran construction: the Wigner D matrix (dmat), the orbital
+time-reversal operator for the magnetic (SP+SO) operations, the basis transform
+to the chosen angular harmonics, and the spin-1/2 phase blocks.
+
+Covers all dmftproj spin-orbit cases:
+
+- non-mixing spin-diagonal bases (complex, cubic, and a fromfile basis whose
+  spin-up and spin-down blocks coincide): the spin-reduced up/up block scaled
+  by the +-(a+g)/2 phase, with the orbital time-reversal operator on the
+  magnetic operations;
+- mixing bases (a fromfile basis that couples spin, e.g. the |j, m_j> basis):
+  the full 2(2l+1) spinor representation P spinrot P^dag, with the spinor
+  time-reversal operator -i sigma_y (x) T applied to the magnetic operations;
+- l = 0 (s) shells: the 2x2 spin phase block.
+
+The mixing fromfile path is the one dft_tools #148 singles out. dmftproj reads
+both cubic and fromfile bases with single-precision CMPLX casts and a 250-column
+line cap (Fortran set_ang_trans.f), and evaluates spin phases in single precision
+(outputqmc.f, setsym.f). Its case.symqmc therefore differs at about 1e-7 from
+this double-precision generator, whose cubic d basis uses exact coefficients.
+"""
+
+import math
+import os
+import numpy as np
+
+# --- angular bases (transpose(P) = <new|m>, m = -l..l) -----------------------
+
+_COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
+
+# standard cubic harmonics, Wien2k convention
+_CUBIC = {
+    1: np.array([[0, 1, 0], [-1j, 0, -1j], [1, 0, -1]], dtype=complex) / 1.0,
+    2: np.array([
+        [0, 0, 1, 0, 0],
+        [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
+        [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5],
+        [0, 2 ** -0.5, 0, -(2 ** -0.5), 0],
+        [0, 2 ** -0.5, 0, 2 ** -0.5, 0],
+    ], dtype=complex),
+}
+
+
+def _reptrans(basis, l):
+    if basis == 'cubic' and l in _CUBIC:
+        return _CUBIC[l]
+    return _COMPLEX[l]
+
+
+# --- Wigner D matrix (dmftproj convention, setsym.f) -------------------------
+
+def _small_d(l, m, n, b):
+    f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
+         (math.factorial(l + n) * math.factorial(l - n))
+    s = 0.0
+    for t in range(0, 2 * l + 1):
+        if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
+            f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
+                 (math.factorial(l - m - t) * math.factorial(m + n + t) *
+                  math.factorial(l - n - t) * math.factorial(t))
+            f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
+            f4 = 1.0 if (2 * t + n + m) == 0 else math.cos(b / 2) ** (2 * t + n + m)
+            s += (-1) ** (l - m - t) * f2 * f3 * f4
+    return math.sqrt(f1) * s
+
+
+def _dmat(l, a, b, c, det):
+    D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+    for m in range(-l, l + 1):
+        for n in range(-l, l + 1):
+            v = np.exp(1j * n * a) * np.exp(1j * m * c) * _small_d(l, m, n, b)
+            if det < -0.5:
+                v *= (-1) ** l
+            D[m + l, n + l] = v
+    return D
+
+
+def _tmat(l):
+    """Complex-conjugation operator in the spherical-harmonic basis,
+    T[m, -m] = (-1)^m (timeinv.f)."""
+    T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+    for m in range(-l, l + 1):
+        T[-m + l, m + l] = (-1) ** m
+    return T
+
+
+def _timeinv_orbital(l, mat):
+    return _tmat(l) @ np.conj(mat)
+
+
+# --- fromfile angular basis ---------------------------------------------------
+
+def _read_fromfile(path, l):
+    """Parse a dmftproj fromfile basis (each line: the coefficients of a new
+    basis vector in {|m,up>, |m,dn>}, m = -l..l, real/imag interleaved; '*'
+    marks the end of an irep). Return (P, mixing) where P = <new|m> is the
+    transform matrix (reptrans.transmat), full 2(2l+1) for a spin-mixing basis
+    or the (2l+1) up/up block otherwise."""
+    n = 2 * (2 * l + 1)
+    rows = []
+    for line in open(path):
+        line = line.rstrip('\n')
+        if not line.strip():
+            continue
+        body = line[1:]
+        vals = [float(x) for x in body.split()][:2 * n]
+        rows.append([vals[2 * k] + 1j * vals[2 * k + 1] for k in range(n)])
+        if len(rows) == n:
+            break
+    R = np.array(rows)                       # R[i, :] = |new_i> in old basis
+    d = 2 * l + 1
+    up_up, up_dn = R[:d, :d], R[:d, d:]
+    dn_up, dn_dn = R[d:, :d], R[d:, d:]
+    mixing = not (np.allclose(dn_dn, up_up) and
+                  np.allclose(up_dn, 0) and np.allclose(dn_up, 0))
+    P = np.conj(R) if mixing else np.conj(up_up)   # <new|m> = conj(<m|new>)
+    return P, mixing
+
+
+# --- input parsing -----------------------------------------------------------
+
+def _read_dmftsym(path):
+    lines = open(path).read().split('\n')
+    nsym = int(lines[0].split()[0])
+    perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)]
+    rest = lines[1 + nsym:]
+    starts = [i for i, l in enumerate(rest) if 'Sym. op.' in l]
+    ops = []
+    for k, s in enumerate(starts):
+        ang = rest[s + 1].split()
+        a, b, c = (math.radians(float(x)) for x in ang[:3])
+        krotm = np.array([[float(x) for x in rest[s + 2 + r].split()] for r in range(3)])
+        ops.append(dict(perm=perms[k], a=a, b=b, c=c, krotm=krotm))
+    return nsym, ops
+
+
+def _read_correlated_shells(indmftpr, struct):
+    """Return the list of correlated shells, one entry per correlated atom, plus
+    the SO flag, from case.indmftpr and case.struct multiplicities. Each shell
+    carries l, the basis name, its transform matrix P = <new|m> and a mixing
+    flag (True for a spin-coupling fromfile basis)."""
+    raw = [l.split('!')[0].strip() for l in open(indmftpr)]
+    raw = [l for l in raw if l != '']
+    nsort = int(raw[0].split()[0])
+    mult = [int(x) for x in raw[1].split()][:nsort]
+    i = 3
+    so = 0
+    shells = []
+    for isort in range(nsort):
+        basis = raw[i].split()[0]
+        i += 1
+        sourcefile = None
+        if basis == 'fromfile':
+            sourcefile = os.path.join(os.path.dirname(indmftpr), raw[i])
+            i += 1
+        l_inc = [int(x) for x in raw[i].split()]
+        i += 1
+        ireps = [int(x) for x in raw[i].split()]
+        i += 1
+        correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
+        if any(n > 0 for n in ireps):
+            i += 1                       # skip the correps line
+        if correlated_ls:
+            so = int(raw[i].split()[0])  # SO flag follows a correlated sort
+            i += 1
+            for l in correlated_ls:
+                if basis == 'fromfile':
+                    P, mixing = _read_fromfile(sourcefile, l)
+                else:
+                    P, mixing = _reptrans(basis, l), False
+                for _ in range(mult[isort]):
+                    shells.append(dict(l=l, basis=basis, P=P, mixing=mixing))
+    return shells, so
+
+
+def write_symqmc(case):
+    """Write <case>.symqmc from <case>.dmftsym, <case>.indmftpr, <case>.struct."""
+    nsym, ops = _read_dmftsym(case + '.dmftsym')
+    shells, so = _read_correlated_shells(case + '.indmftpr', case + '.struct')
+    natom = len(ops[0]['perm'])
+
+    timeinv = []
+    for op in ops:
+        det2 = op['krotm'][0, 0] * op['krotm'][1, 1] - op['krotm'][0, 1] * op['krotm'][1, 0]
+        timeinv.append(1 if (so and det2 < 0.0) else 0)
+
+    with open(case + '.symqmc', 'w') as f:
+        f.write('%6d %6d\n' % (nsym, natom))
+        for op in ops:
+            f.write(''.join('%6d ' % p for p in op['perm']) + '\n')
+        if so:
+            f.write(''.join('%6d ' % t for t in timeinv) + '\n')
+        for isym, op in enumerate(ops):
+            for sh in shells:
+                f.write(_format_matrix(_shell_matrix(op, sh, timeinv[isym])))
+
+
+def _shell_matrix(op, shell, ti):
+    """Spinor symmetry matrix for one correlated shell under one operation."""
+    l = shell['l']
+    if l == 0:
+        return _l0_matrix(op, ti)
+    if shell['mixing']:
+        return _mixing_matrix(op, shell, ti)
+    return _nonmixing_matrix(op, shell, ti)
+
+
+def _phase(op, ti):
+    a, c = op['a'], op['c']
+    return (c - a) if ti else (a + c)   # (g-a) on magnetic ops, else (a+g)
+
+
+def _l0_matrix(op, ti):
+    """s shell: 2x2 spin phase block diag(e, conj(e)) (outputqmc.f l==0)."""
+    e = np.exp(1j * _phase(op, ti) / 2)
+    return np.array([[e, 0], [0, np.conj(e)]], dtype=complex)
+
+
+def _nonmixing_matrix(op, shell, ti):
+    """Spin-diagonal basis: the up/up block scaled by +-(a+g)/2, with the
+    orbital time-reversal operator on the magnetic operations."""
+    l, P = shell['l'], shell['P']
+    rotl = _dmat(l, op['a'], op['b'], op['c'], np.linalg.det(op['krotm']))
+    if ti:
+        rotl = _timeinv_orbital(l, rotl)
+    rotrep = P @ rotl @ np.conj(P.T)
+    e = np.exp(1j * _phase(op, ti) / 2)
+    d = 2 * l + 1
+    mat = np.zeros((2 * d, 2 * d), dtype=complex)
+    mat[:d, :d] = e * rotrep
+    mat[d:, d:] = np.conj(e) * rotrep
+    return mat
+
+
+def _mixing_matrix(op, shell, ti):
+    """Spin-coupling basis: full 2(2l+1) spinor representation. The orbital
+    rotation enters block-diagonal (beta=0) or block-antidiagonal (beta=pi, the
+    magnetic operations); the spinor time-reversal -i sigma_y (x) T is then
+    applied to the magnetic operations (setsym.f, timeinv.f)."""
+    l, P = shell['l'], shell['P']
+    rotl = _dmat(l, op['a'], op['b'], op['c'], np.linalg.det(op['krotm']))
+    e = np.exp(1j * _phase(op, ti) / 2)
+    d = 2 * l + 1
+    spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
+    if ti:                                      # beta = pi, block-antidiagonal
+        spinrot[:d, d:] = e * rotl
+        spinrot[d:, :d] = -np.conj(e) * rotl
+    else:                                       # beta = 0, block-diagonal
+        spinrot[:d, :d] = e * rotl
+        spinrot[d:, d:] = np.conj(e) * rotl
+    rotrep = P @ spinrot @ np.conj(P.T)
+    if ti:
+        tm = _tmat(l)
+        tinv = np.zeros((2 * d, 2 * d), dtype=complex)
+        tinv[:d, d:] = -tm
+        tinv[d:, :d] = tm
+        rotrep = (P @ tinv @ P.T) @ np.conj(rotrep)
+    return rotrep
+
+
+def _format_matrix(mat):
+    out = []
+    for part in (mat.real, mat.imag):
+        for row in part:
+            out.append(''.join('  %.14E' % x for x in row) + '\n')
+    return ''.join(out)
diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt
index 780662c..c43cbed 100644
--- a/test/python/wien2k/CMakeLists.txt
+++ b/test/python/wien2k/CMakeLists.txt
@@ -25,6 +25,23 @@ add_test(NAME Py_wien2k_soc_convert
 set_property(TEST Py_wien2k_soc_convert APPEND PROPERTY ENVIRONMENT
              PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
 
+# Pure-Python case.symqmc generator vs the shipped dmftproj Fortran output.
+# The generator tests read source fixtures and write only to scratch directories.
+add_test(NAME Py_wien2k_symqmc_python
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symqmc_python.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_symqmc_python APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Generator on the two paths the cubic test does not reach: the spin-mixing
+# fromfile (|j,m_j>) basis and an l=0 shell. Reuses CaOs2.dmftsym/struct; the
+# references are the dmftproj matrix data (single precision, ~26 kB total).
+add_test(NAME Py_wien2k_symqmc_mixing
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symqmc_mixing.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_symqmc_mixing APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
 add_test(NAME Py_wien2k_kpath_labels
          COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_kpath_labels.py
          WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
diff --git a/test/python/wien2k/CaOs2.README.md b/test/python/wien2k/CaOs2.README.md
index e6d4d10..8f6ed1b 100644
--- a/test/python/wien2k/CaOs2.README.md
+++ b/test/python/wien2k/CaOs2.README.md
@@ -79,3 +79,17 @@ for conversion of the same fixture bytes. It does not use
 `triqs.utility.h5diff`, which compares arrays at `1e-6` regardless of its
 `precision` argument. The physical overlap check uses `1e-6` because the
 retained Wien2k inputs and native text output have finite precision.
+
+The Python `symqmc` generator is checked separately against the shipped
+`CaOs2.symqmc`, with exact header comparison and `atol=1e-7, rtol=0` for
+complex matrix entries. The maximum complex residual is about `8.6e-8`
+(`6.9e-8` for individual real or imaginary components). In
+`set_ang_trans.f`, the cubic reader converts the template's `0.70710678`
+through `CMPLX` without `KIND=8`, rounding it to `0.7071067690849304`.
+This gives an identity-operation diagonal of `0.9999999657714582`;
+the magnetic operations apply the rounded basis twice. The spin phases in
+`outputqmc.f` are also evaluated through single-precision `CMPLX` and `EXP`.
+The Python generator uses exact cubic d coefficients and double-precision
+phases; its generated matrices are checked for unitarity at `1e-12`.
+It writes only into a temporary directory, so the converter tests always
+read the original native fixture.
diff --git a/test/python/wien2k/CaOs2_jbasis.indmftpr b/test/python/wien2k/CaOs2_jbasis.indmftpr
new file mode 100644
index 0000000..5589430
--- /dev/null
+++ b/test/python/wien2k/CaOs2_jbasis.indmftpr
@@ -0,0 +1,12 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+fromfile
+jbasis_d.dat
+0 0 2 0
+0 0 0 0
+1
+-0.15 0.30
diff --git a/test/python/wien2k/CaOs2_l0.indmftpr b/test/python/wien2k/CaOs2_l0.indmftpr
new file mode 100644
index 0000000..15ca9ef
--- /dev/null
+++ b/test/python/wien2k/CaOs2_l0.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+complex
+2 0 0 0
+0 0 0 0
+1
+-0.15 0.30
diff --git a/test/python/wien2k/jbasis_d.dat b/test/python/wien2k/jbasis_d.dat
new file mode 100644
index 0000000..95c4843
--- /dev/null
+++ b/test/python/wien2k/jbasis_d.dat
@@ -0,0 +1,10 @@
+ -0.89442719 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.44721360 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
+ 0.00000000 0.00000000 -0.77459667 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.63245553 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
+ 0.00000000 0.00000000 0.00000000 0.00000000 -0.63245553 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.77459667 0.00000000 0.00000000 0.00000000
+*0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 -0.44721360 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.89442719 0.00000000
+ 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
+ 0.44721360 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.89442719 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
+ 0.00000000 0.00000000 0.63245553 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.77459667 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
+ 0.00000000 0.00000000 0.00000000 0.00000000 0.77459667 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.63245553 0.00000000 0.00000000 0.00000000
+ 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.89442719 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.44721360 0.00000000
+*0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 1.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
diff --git a/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy b/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy
new file mode 100644
index 0000000..b7e2e01
Binary files /dev/null and b/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy differ
diff --git a/test/python/wien2k/wien2k_symqmc_l0.ref.npy b/test/python/wien2k/wien2k_symqmc_l0.ref.npy
new file mode 100644
index 0000000..679d0e0
Binary files /dev/null and b/test/python/wien2k/wien2k_symqmc_l0.ref.npy differ
diff --git a/test/python/wien2k/wien2k_symqmc_mixing.py b/test/python/wien2k/wien2k_symqmc_mixing.py
new file mode 100644
index 0000000..eb729a8
--- /dev/null
+++ b/test/python/wien2k/wien2k_symqmc_mixing.py
@@ -0,0 +1,88 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.symqmc generator on the two paths the spin-diagonal
+# cubic test (wien2k_symqmc_python) does not exercise:
+#
+#   - a spin-mixing fromfile basis (the |j, m_j> basis), which uses the full
+#     2(2l+1) spinor representation and the spinor time-reversal operator;
+#   - an l = 0 (s) correlated shell, the 2x2 spin phase block.
+#
+# Both reuse the CaOs2 symmetry input (16 operations, 8 of them time-reversal).
+# The reference is the dmftproj Fortran case.symqmc for the same input, stored as
+# the matrix data in single precision: dmftproj reads the fromfile basis with a
+# single-precision CMPLX cast (set_ang_trans.f), so its output there carries a
+# ~1e-7 error. The generator is full double precision, so we compare at 1e-6 and
+# separately assert each Python spinor matrix is unitary.
+
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import symqmc
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _matrix_data(symqmc_path, nsym, natom):
+    """Flat array of the matrix blocks, skipping header/perm/timeflag lines."""
+    toks = iter(open(symqmc_path).read().split())
+    assert int(next(toks)) == nsym and int(next(toks)) == natom
+    for _ in range(nsym * natom + nsym):     # perm rows + the time-reversal flags
+        next(toks)
+    return np.array([float(x) for x in toks])
+
+
+def _run(indmftpr, extra=()):
+    """Generate case.symqmc in a scratch dir from the shared CaOs2 symmetry
+    input and the given indmftpr, returning (matrix data, shells, ops, so)."""
+    tmp = tempfile.mkdtemp()
+    try:
+        for src, dst in (('CaOs2.dmftsym', 'case.dmftsym'),
+                         ('CaOs2.struct', 'case.struct'),
+                         (indmftpr, 'case.indmftpr')) + extra:
+            shutil.copy(os.path.join(HERE, src), os.path.join(tmp, dst))
+        case = os.path.join(tmp, 'case')
+        symqmc.write_symqmc(case)
+        nsym, ops = symqmc._read_dmftsym(case + '.dmftsym')
+        shells, so = symqmc._read_correlated_shells(
+            case + '.indmftpr', case + '.struct')
+        natom = len(ops[0]['perm'])
+        return _matrix_data(case + '.symqmc', nsym, natom), shells, ops, so
+    finally:
+        shutil.rmtree(tmp)
+
+
+def _timeinv(op, so):
+    k = op['krotm']
+    return 1 if (so and k[0, 0] * k[1, 1] - k[0, 1] * k[1, 0] < 0.0) else 0
+
+
+def _check(indmftpr, ref_npy, extra=(), unit_tol=1e-7):
+    data, shells, ops, so = _run(indmftpr, extra)
+    ref = np.load(os.path.join(HERE, ref_npy))
+    assert data.shape == ref.shape, (data.shape, ref.shape)
+    assert np.max(np.abs(data - ref)) < 1e-6, np.max(np.abs(data - ref))
+    for op in ops:
+        for sh in shells:
+            mat = symqmc._shell_matrix(op, sh, _timeinv(op, so))
+            dev = np.max(np.abs(mat @ np.conj(mat.T) - np.eye(mat.shape[0])))
+            assert dev < unit_tol, ('not unitary', dev)
+
+
+# spin-mixing |j, m_j> basis: the path dft_tools #148 singles out
+data, shells, _, _ = _run('CaOs2_jbasis.indmftpr',
+                          (('jbasis_d.dat', 'jbasis_d.dat'),))
+assert shells[0]['mixing'] and shells[0]['P'].shape == (10, 10)
+_check('CaOs2_jbasis.indmftpr', 'wien2k_symqmc_jbasis.ref.npy',
+       (('jbasis_d.dat', 'jbasis_d.dat'),))
+
+# l = 0 (s) correlated shell
+_, shells, _, _ = _run('CaOs2_l0.indmftpr')
+assert shells[0]['l'] == 0
+_check('CaOs2_l0.indmftpr', 'wien2k_symqmc_l0.ref.npy', unit_tol=1e-12)
+
+print('wien2k_symqmc_mixing: ok')
diff --git a/test/python/wien2k/wien2k_symqmc_python.py b/test/python/wien2k/wien2k_symqmc_python.py
new file mode 100644
index 0000000..4ef4e53
--- /dev/null
+++ b/test/python/wien2k/wien2k_symqmc_python.py
@@ -0,0 +1,48 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.symqmc generator (triqs_dftkit.wien2k.symqmc)
+# reproduces the shipped dmftproj Fortran output within its precision limits.
+# Generate in a scratch directory to preserve the SOC converter's fixture.
+
+from pathlib import Path
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k.symqmc import write_symqmc
+
+HERE = Path(__file__).resolve().parent
+
+
+def read_symqmc(path):
+    with open(path) as f:
+        tokens = f.read().split()
+    # This fixture has 16 operations, 3 atoms and two 10-dimensional shells.
+    header_size = 2 + 16 * 3 + 16
+    header = np.array(tokens[:header_size], dtype=int)
+    blocks = np.array(tokens[header_size:], dtype=float).reshape(16, 2, 2, 10, 10)
+    return header, blocks[:, :, 0] + 1j * blocks[:, :, 1]
+
+
+with tempfile.TemporaryDirectory() as tmp:
+    case = Path(tmp) / 'CaOs2'
+    for suffix in ('dmftsym', 'indmftpr', 'struct'):
+        shutil.copyfile(HERE / f'CaOs2.{suffix}', f'{case}.{suffix}')
+    write_symqmc(str(case))
+    header, matrices = read_symqmc(f'{case}.symqmc')
+
+ref_header, ref_matrices = read_symqmc(HERE / 'CaOs2.symqmc')
+np.testing.assert_array_equal(header, ref_header)
+# set_ang_trans.f reads 0.70710678 through single-precision CMPLX; magnetic
+# operations use the rounded basis twice. outputqmc.f also evaluates the spin
+# phase in single precision. The complex matrix residual is about 8.6e-8.
+np.testing.assert_allclose(matrices, ref_matrices, atol=1e-7, rtol=0)
+for operation in matrices:
+    for matrix in operation:
+        np.testing.assert_allclose(matrix @ matrix.conj().T, np.eye(10),
+                                   atol=1e-12, rtol=0)
+
+print('wien2k_symqmc_python: ok')