diff --git a/python/triqs_dftkit/wien2k/_dmftproj.py b/python/triqs_dftkit/wien2k/_dmftproj.py
index 57f9c4a..7e16a99 100644
--- a/python/triqs_dftkit/wien2k/_dmftproj.py
+++ b/python/triqs_dftkit/wien2k/_dmftproj.py
@@ -103,7 +103,9 @@ def read_two_complex(r):
 _COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
 
 _CUBIC = {
-    1: np.array([[2**-0.5, 0, -2**-0.5], [-1j*2**-0.5, 0, -1j*2**-0.5], [0, 1, 0]], dtype=complex),
+    1: np.array([[2 ** -0.5, 0, -(2 ** -0.5)],
+                 [-1j * 2 ** -0.5, 0, -1j * 2 ** -0.5],
+                 [0, 1, 0]], dtype=complex),
     2: np.array([
         [0, 0, 1, 0, 0],
         [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
@@ -225,6 +227,58 @@ def mixing_timeinv_op(l, P):
     return P @ tinv @ P.T
 
 
+def euler(R):
+    """Euler angles (a, b, c) of the proper rotation R, in the WIEN2k
+    convention of set_rotloc.f/SRC_lapwdm euler.f: R is the transpose of a
+    krotm matrix. a, c in [0, 2pi), b in [0, pi]."""
+    z = np.array([0.0, 0.0, 1.0])
+    y = np.array([0.0, 1.0, 0.0])
+    yyy = R[:, 1]
+    zz = R[:, 2]
+    yy = np.cross(z, zz)
+    y_norm = math.sqrt(yy @ yy)
+    if y_norm < 1e-10:
+        a = math.acos(max(-1.0, min(1.0, y @ yyy)))
+        if z @ zz > 0.0:
+            b = c = 0.0
+            if yyy[0] > 0.0:
+                a = 2 * math.pi - a
+        else:
+            c, a, b = a, 0.0, math.pi
+            if yyy[0] < 0.0:
+                c = 2 * math.pi - c
+        return a, b, c
+    yy = yy / y_norm
+    a = math.acos(max(-1.0, min(1.0, y @ yy)))
+    b = math.acos(max(-1.0, min(1.0, z @ zz)))
+    c = math.acos(max(-1.0, min(1.0, yy @ yyy)))
+    if yy[0] > 0.0:
+        a = 2 * math.pi - a
+    if np.cross(yy, yyy) @ zz < 0.0:
+        c = 2 * math.pi - c
+    return a, b, c
+
+
+def _rloc_spinor_s(ref, op, timeinv):
+    """2x2 spinor Rloc of an s shell under SP+SO (outputqmc.f, l=0): built
+    from the Euler angles of the composed proper rotation
+    krotm(op) krotm(ref), which fixes its SU(2) sign. For a magnetic op the
+    columns carry the time reversal -i sigma_y."""
+    a, b, g = euler((op['krotm'] @ ref['krotm']).T)
+    s = np.zeros((2, 2), dtype=complex)
+    if timeinv:
+        s[1, 0] = np.exp(1j * (a + g) / 2) * math.cos(b / 2)
+        s[0, 1] = -np.conj(s[1, 0])
+        s[1, 1] = -np.exp(-1j * (a - g) / 2) * math.sin(b / 2)
+        s[0, 0] = np.conj(s[1, 1])
+    else:
+        s[0, 0] = np.exp(1j * (a + g) / 2) * math.cos(b / 2)
+        s[1, 1] = np.conj(s[0, 0])
+        s[0, 1] = np.exp(-1j * (a - g) / 2) * math.sin(b / 2)
+        s[1, 0] = -np.conj(s[0, 1])
+    return s
+
+
 def rotloc_rotl_so(l, ref, ops, iatom, iref):
     """The composed 2(2l+1) Rloc rotation rotloc(iatom)%rotl(l) under SP+SO
     (setsym.f:496-528 + set_rotloc.f): the representative spinor rotloc
@@ -233,6 +287,12 @@ def rotloc_rotl_so(l, ref, ops, iatom, iref):
     Returns (rotl, timeinv); callers apply their own basis transform (the full
     transmat for a mixing basis, blkdiag(transmat, transmat) otherwise)."""
     d = 2 * l + 1
+    op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+    det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+            - op['krotm'][0, 1] * op['krotm'][1, 0])
+    timeinv = det2 < 0.0
+    if l == 0:
+        return _rloc_spinor_s(ref, op, timeinv), timeinv
     Dref = dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
     f = (ref['a'] + ref['g']) / 2.0
     spmt = np.zeros((2, 2), dtype=complex)
@@ -247,10 +307,6 @@ def rotloc_rotl_so(l, ref, ops, iatom, iref):
     rotl[:d, d:] = spmt[0, 1] * Dref
     rotl[d:, :d] = spmt[1, 0] * Dref
 
-    op = next(o for o in ops if o['perm'][iref - 1] == iatom)
-    det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
-            - op['krotm'][0, 1] * op['krotm'][1, 0])
-    timeinv = det2 < 0.0
     srot_phase = (op['g'] - op['a']) if timeinv else (op['a'] + op['g'])
     rotl_sym = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
     if timeinv:
@@ -263,6 +319,13 @@ def rotloc_rotl_so(l, ref, ops, iatom, iref):
     return rotl, timeinv
 
 
+def rotloc_orbital(l, ref, ops, iatom, iref):
+    """Orbital part of the full local rotation R[isym] * rotloc_ref."""
+    op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+    return dmat(l, op['a'], op['b'], op['g'], float(op['iprop'])) @ \
+        dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
+
+
 def mixing_rotrep(op, l, P, ti):
     """The full 2(2l+1) spinor representation D(R)_{new_i} = P spinrot P^dag of
     one symmetry operation in a spin-coupling (mixing) basis, with the spinor
@@ -625,18 +688,6 @@ def fmt_scalar(x):
     return '%.16f' % x if abs(x) < 1e5 else '%.16E' % x
 
 
-def require_identity_struct(case):
-    """Fail explicitly until representative struct rotations are ported."""
-    with open(case + '.struct') as f:
-        lines = f.readlines()
-    for i, line in enumerate(lines):
-        if 'LOCAL ROT MATRIX:' in line:
-            rows = [line.split(':', 1)[1], *lines[i + 1:i + 3]]
-            rotation = np.array([[to_float(x) for x in row.split()] for row in rows])
-            if not np.allclose(rotation, np.eye(3), atol=1e-7, rtol=0):
-                raise NotImplementedError('non-identity struct local rotations are not implemented yet')
-
-
 def validate_spin_windows(spins, info, windows, below=None):
     """Require valid selected windows and the same absolute bands in both spins."""
     if len(spins) != 2:
diff --git a/python/triqs_dftkit/wien2k/ctqmcout.py b/python/triqs_dftkit/wien2k/ctqmcout.py
index 22f6e16..0243a90 100644
--- a/python/triqs_dftkit/wien2k/ctqmcout.py
+++ b/python/triqs_dftkit/wien2k/ctqmcout.py
@@ -60,9 +60,9 @@
 from scipy.linalg.lapack import zheev
 from scipy.linalg.blas import zgemm
 
-from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_spin_windows, band_index_window, dmat, fmt_scalar, read_almblm,
+from ._dmftproj import (build_orbs as _build_orbs, validate_spin_windows, band_index_window, dmat, fmt_scalar, read_almblm,
                         read_dmftsym, read_fromfile, read_indmftpr, reptrans,
-                        rotloc_rotl_so, select_band_window, select_window,
+                        rotloc_orbital, rotloc_rotl_so, select_band_window, select_window,
                         tmat, write_row)
 
 
@@ -102,34 +102,23 @@ def _build_crorbs(info):
 
 # --- Rloc spinor representation (set_rotloc.f / setsym.f) ---------------------
 
-def _rloc_rotrep(op, l, transmat, ifSO):
-    """rotloc%rotrep(l)%mat, the Rloc rotation in the new basis for a non-mixing
-    shell, with identity struct local rotation (rotloc_ref Euler = 0).
-
-    Under SO it is the 2(2l+1) spinor rotation: rotloc%rotl =
-    blkdiag(ephase*D, conj(ephase)*D), rotrep = S rotl S^H, S =
-    blkdiag(transmat, transmat). Without SO srot%timeinv is always false, there
-    is no spin phase, and it reduces to the bare (2l+1) transmat D transmat^H
+def _rloc_rotrep(l, ref, ops, iatom, iref, transmat, ifSO):
+    """rotloc(iatom)%rotrep(l)%mat for a non-mixing shell: the composed Rloc
+    spinor rotation rotloc_ref then the symmetry op mapping iref onto iatom
+    (set_rotloc.f / setsym.f), put into the new basis with S = blkdiag(transmat,
+    transmat). Shares rotloc_rotl_so with the mixing path so the non-trivial
+    struct reference rotation rotloc_ref (e.g. Te, beta=pi/2) is not dropped.
+    Without SO it reduces to the bare (2l+1) D(op)D(ref) transmat sandwich
     (set_rotloc.f non-SO branch). Returns (rotrep, timeinv)."""
-    a, b, g, iprop = op['a'], op['b'], op['g'], op['iprop']
-    D = dmat(l, a, b, g, float(iprop))
     if not ifSO:
+        D = rotloc_orbital(l, ref, ops, iatom, iref)
         return transmat @ D @ np.conj(transmat.T), False
-    krotm = op['krotm']
-    det2 = krotm[0, 0] * krotm[1, 1] - krotm[0, 1] * krotm[1, 0]
-    timeinv = det2 < 0.0
-    phase = (g - a) if timeinv else (a + g)
-    if timeinv:
-        D = tmat(l) @ np.conj(D)       # setsym.f:320-326 orbital time reversal
-    ephase = np.exp(1j * phase / 2)
+    rotl, timeinv = rotloc_rotl_so(l, ref, ops, iatom, iref)
     d = 2 * l + 1
-    rotl = np.zeros((2 * d, 2 * d), dtype=complex)
-    rotl[:d, :d] = ephase * D
-    rotl[d:, d:] = np.conj(ephase) * D
     S = np.zeros((2 * d, 2 * d), dtype=complex)
     S[:d, :d] = transmat
     S[d:, d:] = transmat
-    rotrep = S @ rotl @ np.conj(S.T)
+    rotrep = (S @ rotl @ S.T) if timeinv else (S @ rotl @ np.conj(S.T))
     return rotrep, timeinv
 
 
@@ -148,7 +137,6 @@ def write_ctqmcout(case):
     """Read <case>.almblm{up,dn}, <case>.indmftpr, <case>.struct,
     <case>.dmftsym and write <case>.ctqmcout in the dmftproj format."""
     info = read_indmftpr(case + '.indmftpr')
-    require_identity_struct(case)
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
 
     up, dn = case + '.almblmup', case + '.almblmdn'
@@ -230,12 +218,6 @@ def write_ctqmcout(case):
     # rot_projectmat applies dmat(l, a, b, g, iprop) in the |lm> basis
     # (rot_projectmat.f:59-65). For atom 2 this is a non-trivial in-plane
     # rotation (a=270) that mixes the m-rows; omitting it flips the m=1,2 rows.
-    rotloc_op = {}
-    for icr, cr in enumerate(crorbs):
-        iref = sum(info['mult'][:cr['sort'] - 1]) + 1
-        rotloc_op[icr] = next(o for o in ops
-                              if o['perm'][iref - 1] == cr['atom'])
-
     # ---- raw correlated projector mat_rep, then Loewdin orthonormalize ----
     # Non-mixing: mat_rep[(icr, ik, is)] -> (2l+1, nbsel) per spin block.
     # Mixing: mat_rep[(icr, ik)] -> (2*(2l+1), nbsel), the stacked spinor block.
@@ -247,8 +229,8 @@ def write_ctqmcout(case):
         d = 2 * l + 1
         transmat = transmats[(l, sort)]
         ismix = mixing[(l, sort)]
-        op = rotloc_op[icr]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[sort - 1]
+        rot = rotloc_orbital(l, ref, ops, atom, sum(info['mult'][:sort - 1]) + 1)
         for ik in range(nk):
             incl, nb_bot, nb_top = windows[ik]
             if not incl:
@@ -328,8 +310,8 @@ def write_ctqmcout(case):
             rloc_blocks.append(_rloc_rotrep_mixing(
                 l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat))
         else:
-            op = next(o for o in ops if o['perm'][iref - 1] == cr['atom'])
-            rloc_blocks.append(_rloc_rotrep(op, l, transmat, ifSO))
+            rloc_blocks.append(_rloc_rotrep(
+                l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat, ifSO))
 
     # ---- write ----
     with open(case + '.ctqmcout', 'w') as f:
diff --git a/python/triqs_dftkit/wien2k/outband.py b/python/triqs_dftkit/wien2k/outband.py
index 65e5787..74814ad 100644
--- a/python/triqs_dftkit/wien2k/outband.py
+++ b/python/triqs_dftkit/wien2k/outband.py
@@ -51,8 +51,8 @@
 import numpy as np
 from scipy.linalg.blas import zgemm
 
-from ._dmftproj import (require_identity_struct, validate_spin_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
-                        read_indmftpr, reptrans, select_window, to_float,
+from ._dmftproj import (validate_spin_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
+                        read_indmftpr, reptrans, rotloc_orbital, select_window, to_float,
                         write_row)
 from .ctqmcout import _build_crorbs, _sqrt_inv
 from .parproj import (_build_matn_rep, _build_matn_rep_mixing, _build_orbs)
@@ -108,7 +108,6 @@ def write_outband(case):
     info, eferm = _read_indmftpr_band(case + '.indmftpr')
     if not info['so']:
         raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
-    require_identity_struct(case)
     if info['proj_mode'] != 0:
         raise NotImplementedError('this output currently supports projection mode 0')
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
@@ -158,12 +157,6 @@ def write_outband(case):
 
     # rot_projectmat local rotation: the op mapping the sort representative onto
     # this atom (set_projections.f via rot_projectmat). Same as ctqmcout/parproj.
-    rotloc_op = {}
-    for o in orbs:
-        iref = sum(info['mult'][:o['sort'] - 1]) + 1
-        rotloc_op[o['atom']] = next(op for op in ops
-                                    if op['perm'][iref - 1] == o['atom'])
-
     # ---- correlated projector mat_rep, then Loewdin orthonormalize ----
     mat_rep = {}
     for icr, cr in enumerate(crorbs):
@@ -173,8 +166,8 @@ def write_outband(case):
         d = 2 * l + 1
         transmat = transmats[(l, sort)]
         ismix = mixing[(l, sort)]
-        op = rotloc_op[atom]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[sort - 1]
+        rot = rotloc_orbital(l, ref, ops, atom, sum(info['mult'][:sort - 1]) + 1)
         for ik in range(nk):
             incl, nb_bot, nb_top = windows[ik]
             if not incl:
@@ -236,8 +229,9 @@ def write_outband(case):
     for o in orbs:
         l = o['l']
         transmat = transmats[(l, o['sort'])]
-        op = rotloc_op[o['atom']]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[o['sort'] - 1]
+        iref = sum(info['mult'][:o['sort'] - 1]) + 1
+        rot = rotloc_orbital(l, ref, ops, o['atom'], iref)
         if mixing[(l, o['sort'])]:
             matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
                 o, info, spins, windows, transmat, rot)
diff --git a/python/triqs_dftkit/wien2k/parproj.py b/python/triqs_dftkit/wien2k/parproj.py
index 22488a0..69322b4 100644
--- a/python/triqs_dftkit/wien2k/parproj.py
+++ b/python/triqs_dftkit/wien2k/parproj.py
@@ -57,9 +57,9 @@
 import os
 import numpy as np
 
-from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_spin_windows, band_index_window, dmat, mixing_rotrep, read_almblm, read_dmftsym,
-                        read_fromfile, read_indmftpr, reptrans,
-                        rotloc_rotl_so, select_band_window, select_window, tmat, write_row)
+from ._dmftproj import (build_orbs as _build_orbs, validate_spin_windows, band_index_window, dmat, mixing_rotrep, paramag_time_op,
+                        read_almblm, read_dmftsym, read_fromfile, read_indmftpr, reptrans,
+                        rotloc_orbital, rotloc_rotl_so, select_band_window, select_window, tmat, write_row)
 
 
 def _op_timeinv(op, ifSO=True):
@@ -117,15 +117,14 @@ def _srot_rotrep_nonmixing(op, l, transmat, ifSO=True):
 
 # --- rotloc rotrep (set_rotloc.f) under SP+SO, non-mixing --------------------
 
-def _rotloc_rotrep_nonso(orb, ops, info, transmat):
+def _rotloc_rotrep_nonso(orb, ops, info, ref, transmat):
     """rotloc(iatom)%rotrep(l)%mat without SO: the bare (2l+1) Rloc rotation in
     the new basis, transmat D(Rloc)_{lm} transmat^H with D(Rloc) the symmetry op
     mapping the sort representative onto this atom (set_rotloc.f non-SO branch).
     srot%timeinv is always false under non-SO, so timeinv=False."""
     l = orb['l']
     iref = sum(info['mult'][:orb['sort'] - 1]) + 1
-    op = next(o for o in ops if o['perm'][iref - 1] == orb['atom'])
-    D = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+    D = rotloc_orbital(l, ref, ops, orb['atom'], iref)
     return transmat @ D @ np.conj(transmat.T), False
 
 
@@ -389,6 +388,14 @@ def _symmetrize_densmat_nonso(orbs, info, ops, nsym, srot_rotreps, dens_raw):
     return out
 
 
+def _add_timeinv_nonso(blocks, l, transmat):
+    """add_timeinv for a paramagnetic non-SO shell (timeinv.f, nsp=1): the
+    single spin block D becomes (D + U conj(D) U^H)/2, U the time-reversal
+    operator in the shell basis."""
+    U = paramag_time_op(dict(l=l, P=transmat))
+    return [(D + U @ np.conj(D) @ np.conj(U.T)) / 2 for D in blocks]
+
+
 def _rotdens_densmat_nonso(blocks, rotrep_loc):
     """rotdens_mat for the non-SO non-mixing path (rot_dens.f:197-228):
     each (2l+1) spin block is rotated as conj(Rloc^T) D Rloc with the bare
@@ -422,7 +429,6 @@ def write_parproj(case):
     """Read <case>.almblm{up,dn}, <case>.indmftpr, <case>.struct,
     <case>.dmftsym and write <case>.parproj in the dmftproj format."""
     info = read_indmftpr(case + '.indmftpr')
-    require_identity_struct(case)
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
 
     spin_files = [case + '.almblmup', case + '.almblmdn']
@@ -476,12 +482,6 @@ def write_parproj(case):
 
     # rot_projectmat local rotation: the op mapping the representative atom of
     # the sort onto this atom (set_projections.f via rot_projectmat).
-    rotloc_op = {}
-    for o in orbs:
-        iref = sum(info['mult'][:o['sort'] - 1]) + 1
-        rotloc_op[o['atom']] = next(op for op in ops
-                                    if op['perm'][iref - 1] == o['atom'])
-
     # ---- projectors matn_rep per orbital (energy window for block A, full
     # band range for the density matrix) ----
     matn_reps = {}
@@ -489,8 +489,9 @@ def write_parproj(case):
     for o in orbs:
         l = o['l']
         transmat = transmats[(l, o['sort'])]
-        op = rotloc_op[o['atom']]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[o['sort'] - 1]
+        iref = sum(info['mult'][:o['sort'] - 1]) + 1
+        rot = rotloc_orbital(l, ref, ops, o['atom'], iref)
         if mixing[(l, o['sort'])]:
             matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
                 o, info, spins, windows, transmat, rot)
@@ -506,11 +507,11 @@ def write_parproj(case):
     rotloc_rotrep = {}
     for o in orbs:
         transmat = transmats[(o['l'], o['sort'])]
+        ref = rotloc_ref[o['sort'] - 1]
         if not ifSO:
             rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_nonso(
-                o, ops, info, transmat)
+                o, ops, info, ref, transmat)
         else:
-            ref = rotloc_ref[o['sort'] - 1]
             rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_so(
                 o, ops, info, ref, transmat, mixing[(o['l'], o['sort'])])
 
@@ -544,9 +545,12 @@ def write_parproj(case):
                                                  srot_rotreps, dens_raw)
             for o in nonmix_orbs:
                 key = (o['atom'], o['l'])
+                blocks = dens_sym[key]
+                if not ifSP:
+                    blocks = _add_timeinv_nonso(
+                        blocks, o['l'], transmats[(o['l'], o['sort'])])
                 rotrep_loc, _ = rotloc_rotrep[key]
-                densprint[key] = _rotdens_densmat_nonso(
-                    dens_sym[key], rotrep_loc)
+                densprint[key] = _rotdens_densmat_nonso(blocks, rotrep_loc)
 
     mixing_groups = {}
     for o in orbs:
diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt
index 9ecee85..f3dacea 100644
--- a/test/python/wien2k/CMakeLists.txt
+++ b/test/python/wien2k/CMakeLists.txt
@@ -45,6 +45,24 @@ add_test(NAME Py_wien2k_symqmc_mixing
 set_property(TEST Py_wien2k_symqmc_mixing APPEND PROPERTY ENVIRONMENT
              PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
 
+# Unitarity on a non-centrosymmetric spin-orbit cell (tellurium, p shell).
+# Native reference comparisons are covered by Py_wien2k_te_python below.
+add_test(NAME Py_wien2k_symqmc_noncentrosym
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symqmc_noncentrosym.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_symqmc_noncentrosym APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Projector and symmetry generators on tellurium against native dmftproj
+# references: the composed local frame of a non-centrosymmetric cell, SO s and p
+# shells, a magnetic operation mapping the atoms, and the non-SO paths. Reads its
+# inputs from the source directory.
+add_test(NAME Py_wien2k_te_python
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_te_python.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_te_python APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
 # Pure-Python case.oubwin band-window generator vs the dmftproj Fortran output.
 # Reads the existing compressed almblm inputs from the source fixtures and
 # compares both spin outputs against CaOs2.oubwinup in a scratch directory.
diff --git a/test/python/wien2k/Te.README.md b/test/python/wien2k/Te.README.md
new file mode 100644
index 0000000..fe65fc8
--- /dev/null
+++ b/test/python/wien2k/Te.README.md
@@ -0,0 +1,96 @@
+# Te projector fixture
+
+`Te.almblmup.gz` and `Te.almblmdn.gz` are the Wien2k `lapw2 -almd` projector
+inputs for trigonal tellurium at a single k-point. Te is non-centrosymmetric
+(P3_121) with three equivalent atoms, so these files exercise the composed
+local frame in `ctqmcout`, `parproj` and `outband`, which the centrosymmetric
+fixtures never reach. They pair with the shipped `Te.struct`, `Te.dmftsym` and
+`Te.indmftpr`, which are the files of the same calculation with trailing
+whitespace removed.
+
+## Producer provenance
+
+- Wien2k 19.1 (gfortran build), serial, no MPI.
+- Structure: `Te.struct`, P3_121, a = 9.41 bohr, c = 11.199 bohr, one
+  inequivalent Te with RMT = 2.10 bohr; spin-orbit with M || z, which leaves
+  the six operations of the space group; the three two-fold axes
+  perpendicular to z reverse M and enter combined with time reversal.
+- Density: PBE, spin-polarized scalar-relativistic self-consistent run on
+  12 irreducible k-points, converged to 1e-4 Ry in energy and 1e-3 e in
+  charge. The solution is non-magnetic (cell moment -0.001 μB).
+- Projector step on that density, spin-polarized with spin-orbit coupling, at
+  the single general k-point (3, 7, 11)/35 in reciprocal-lattice units with
+  weight 1. `Te.klist`:
+
+  ```
+           1         3         7        11        35  1.0 -7.0  1.5         1 k, div: (  1  1  1)
+  END
+  ```
+
+- Band range: 34 spin-orbit bands, 31 to 64 (Te 5s, 5p and the lowest
+  conduction bands). The Te 4d semicore bands are left out with `EMIN = -1.0`
+  in `Te.in2c` (NE stays 48, since the Fermi search counts all bands). The
+  spin-orbit basis is limited to scalar-relativistic states below 0.85 Ry
+  (`Te.inso` Emax; 1.9 Ry in the original projector run), which moves the
+  energies of bands 31 to 64 by at most 6.3e-5 Ry.
+- Fermi level: 0.26554 Ry from Fermi smearing (`TEMP 0.002` in `Te.in2c`),
+  because the tetrahedron method needs a mesh. The full-mesh run gave
+  0.24443 Ry. With the shipped `Te.indmftpr` (cubic p, -0.50 to 0.50 Ry) the
+  window holds bands 37 to 62, as in the full-mesh run.
+
+Commands, in a case directory named `Te`, after writing the klist and the
+changed `Te.in2c` and `Te.inso` lines above:
+
+```sh
+x lapw1 -up; x lapw1 -dn; x lapwso -up
+x lapw2 -up -c -so -almd; x lapw2 -dn -c -so -almd
+gzip -9 -n Te.almblmup Te.almblmdn
+```
+
+## Native references
+
+`wien2k_te_python.py` compares the Python generators against native dmftproj
+output. dmftproj derives the case name from the working directory, so each
+reference is produced in a directory named after it:
+
+| reference         | indmftpr          | dmftsym          | almblm         | dmftproj flags | files                                     |
+|-------------------|-------------------|------------------|----------------|----------------|-------------------------------------------|
+| `Te`              | `Te.indmftpr`     | `Te.dmftsym`     | up, dn         | `-sp -so`      | ctqmcout, symqmc, parproj                 |
+| `Te_sp`           | `Te_sp.indmftpr`  | `Te.dmftsym`     | up, dn         | `-sp -so`      | ctqmcout, symqmc, parproj                 |
+| `Te_sp_mag`       | `Te_sp.indmftpr`  | `Te_mag.dmftsym` | up, dn         | `-sp -so`      | ctqmcout, parproj                         |
+| `Te_noso`         | `Te_noso.indmftpr`| `Te.dmftsym`     | up, dn         | `-sp`          | ctqmcout                                  |
+| `Te_noso_paramag` | `Te_noso.indmftpr`| `Te.dmftsym`     | up as `.almblm`| none           | ctqmcout, symqmc, sympar, parproj         |
+
+`Te_mag.dmftsym` is `Te.dmftsym` with the operations reordered (4, 5, 6, 1, 2, 3)
+so that a two-fold axis, magnetic under spin-orbit, maps the representative
+onto each atom. `$WIENROOT/SRC_templates` must hold the full-precision
+`case.cf_p_cubic` from `fortran/dmftproj/SRC_templates`. From
+`test/python/wien2k`, with `DMFTPROJ` the built `fortran/dmftproj/dmftproj`:
+
+```sh
+src=$PWD; out=$(mktemp -d)
+ref() {  # ref <reference> <indmftpr> <dmftsym> <sp|pm> <flags> <files...>
+  name=$1 inp=$2 sym=$3 spin=$4 flags=$5; shift 5
+  mkdir "$out/$name" && cd "$out/$name"
+  cp "$src/$inp.indmftpr" "$name.indmftpr"
+  cp "$src/Te.struct" "$name.struct"
+  cp "$src/$sym.dmftsym" "$name.dmftsym"
+  if [ "$spin" = sp ]; then
+    for s in up dn; do gunzip -c "$src/Te.almblm$s.gz" > "$name.almblm$s"; done
+  else
+    gunzip -c "$src/Te.almblmup.gz" > "$name.almblm"
+  fi
+  "$DMFTPROJ" $flags > run.log
+  for f in "$@"; do gzip -9 -n -c "$name.$f" > "$name.$f.gz"; done
+  cd "$src"
+}
+ref Te              Te      Te     sp "-sp -so" ctqmcout symqmc parproj
+ref Te_sp           Te_sp   Te     sp "-sp -so" ctqmcout symqmc parproj
+ref Te_sp_mag       Te_sp   Te_mag sp "-sp -so" ctqmcout parproj
+ref Te_noso         Te_noso Te     sp "-sp"     ctqmcout
+ref Te_noso_paramag Te_noso Te     pm ""        ctqmcout symqmc sympar parproj
+for f in "$out"/*/*.gz; do cmp "$f" "$(basename "$f")"; done
+```
+
+The last line checks the regenerated references byte for byte against the
+committed ones.
diff --git a/test/python/wien2k/Te.almblmdn.gz b/test/python/wien2k/Te.almblmdn.gz
new file mode 100644
index 0000000..097989c
Binary files /dev/null and b/test/python/wien2k/Te.almblmdn.gz differ
diff --git a/test/python/wien2k/Te.almblmup.gz b/test/python/wien2k/Te.almblmup.gz
new file mode 100644
index 0000000..4f11866
Binary files /dev/null and b/test/python/wien2k/Te.almblmup.gz differ
diff --git a/test/python/wien2k/Te.ctqmcout.gz b/test/python/wien2k/Te.ctqmcout.gz
new file mode 100644
index 0000000..000d122
Binary files /dev/null and b/test/python/wien2k/Te.ctqmcout.gz differ
diff --git a/test/python/wien2k/Te.dmftsym b/test/python/wien2k/Te.dmftsym
new file mode 100644
index 0000000..a0234c9
--- /dev/null
+++ b/test/python/wien2k/Te.dmftsym
@@ -0,0 +1,49 @@
+           6
+           1           2           3
+           2           3           1
+           3           1           2
+           1           3           2
+           2           1           3
+           3           2           1
+
+Sym. op. :  1
+   0.0   0.0   0.0  1       - euler angles: a,b,c; iprop
+   1.0000000000000000        0.0000000000000000        0.0000000000000000
+   0.0000000000000000       0.99999999999999989        0.0000000000000000
+   0.0000000000000000        0.0000000000000000       0.99999999999999989
+
+Sym. op. :  2
+ 240.0   0.0   0.0  1       - euler angles: a,b,c; iprop
+ -0.50000000000000000      -0.86602540378443849        0.0000000000000000
+  0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000       0.99999999999999989
+
+Sym. op. :  3
+ 120.0   0.0   0.0  1       - euler angles: a,b,c; iprop
+ -0.50000000000000000       0.86602540378443849        0.0000000000000000
+ -0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000       0.99999999999999989
+
+Sym. op. :  4
+   0.0 180.0 120.0  1       - euler angles: a,b,c; iprop
+  0.50000000000000000       0.86602540378443849        0.0000000000000000
+  0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000      -0.99999999999999989
+
+Sym. op. :  5
+   0.0 180.0   0.0  1       - euler angles: a,b,c; iprop
+  -1.0000000000000000        0.0000000000000000        0.0000000000000000
+   0.0000000000000000       0.99999999999999989        0.0000000000000000
+   0.0000000000000000        0.0000000000000000      -0.99999999999999989
+
+Sym. op. :  6
+   0.0 180.0 240.0  1       - euler angles: a,b,c; iprop
+  0.50000000000000000      -0.86602540378443849        0.0000000000000000
+ -0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000      -0.99999999999999989
+ Global->local coordinates rotation matrices
+           1
+  -0.0000000000000000       0.50000000000000000       0.86602540378443860
+   0.0000000000000000      -0.86602540378443860       0.50000000000000000
+   1.0000000000000000        0.0000000000000000        0.0000000000000000
+   0.0  90.0 150.0  1       - euler angles: a,b,c; iprop
diff --git a/test/python/wien2k/Te.indmftpr b/test/python/wien2k/Te.indmftpr
new file mode 100644
index 0000000..f2fa44d
--- /dev/null
+++ b/test/python/wien2k/Te.indmftpr
@@ -0,0 +1,8 @@
+1
+3
+3
+cubic
+0 2 0 0
+0 0 0 0
+1
+-0.50 0.50
diff --git a/test/python/wien2k/Te.parproj.gz b/test/python/wien2k/Te.parproj.gz
new file mode 100644
index 0000000..73236f0
Binary files /dev/null and b/test/python/wien2k/Te.parproj.gz differ
diff --git a/test/python/wien2k/Te.struct b/test/python/wien2k/Te.struct
new file mode 100644
index 0000000..d0e338e
--- /dev/null
+++ b/test/python/wien2k/Te.struct
@@ -0,0 +1,37 @@
+ Tellur  (P3121)                         s-o calc. M||  0.00  0.00  1.00
+H                            1 NEW
+             RELA
+  9.410000  9.410000 11.199000  0.000000  0.000000  0.000000
+ATOM  -1: X=0.76468000 Y=0.76468000 Z=0.00000000
+          MULT= 3          ISPLIT= 8
+      -1: X=0.23532000 Y=0.00000000 Z=0.33333330
+      -1: X=0.00000000 Y=0.23532000 Z=0.66666670
+Te         NPT=  381  R0=.000005000 RMT=   2.10000   Z:  52.00000
+LOCAL ROT MATRIX:   -0.0000000 0.5000000 0.8660254
+                     0.0000000-0.8660254 0.5000000
+                     1.0000000 0.0000000 0.0000000
+   6      NUMBER OF SYMMETRY OPERATIONS
+ 1 0 0 0.00000000
+ 0 1 0 0.00000000
+ 0 0 1 0.00000000
+       1   A   1 so. oper.  type  orig. index
+ 0-1 0 0.00000000
+ 1-1 0 0.00000000
+ 0 0 1 0.33333330
+       2   A   3
+-1 1 0 0.00000000
+-1 0 0 0.00000000
+ 0 0 1 0.66666670
+       3   A   5
+ 0 1 0 0.00000000
+ 1 0 0 0.00000000
+ 0 0-1 0.00000000
+       4   B   2
+-1 0 0 0.00000000
+-1 1 0 0.00000000
+ 0 0-1 0.33333330
+       5   B   4
+ 1-1 0 0.00000000
+ 0-1 0 0.00000000
+ 0 0-1 0.66666670
+       6   B   6
diff --git a/test/python/wien2k/Te.symqmc.gz b/test/python/wien2k/Te.symqmc.gz
new file mode 100644
index 0000000..94785e7
Binary files /dev/null and b/test/python/wien2k/Te.symqmc.gz differ
diff --git a/test/python/wien2k/Te_mag.dmftsym b/test/python/wien2k/Te_mag.dmftsym
new file mode 100644
index 0000000..88c6ed8
--- /dev/null
+++ b/test/python/wien2k/Te_mag.dmftsym
@@ -0,0 +1,49 @@
+           6
+           1           3           2
+           2           1           3
+           3           2           1
+           1           2           3
+           2           3           1
+           3           1           2
+
+Sym. op. :  1
+   0.0 180.0 120.0  1       - euler angles: a,b,c; iprop
+  0.50000000000000000       0.86602540378443849        0.0000000000000000
+  0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000      -0.99999999999999989
+
+Sym. op. :  2
+   0.0 180.0   0.0  1       - euler angles: a,b,c; iprop
+  -1.0000000000000000        0.0000000000000000        0.0000000000000000
+   0.0000000000000000       0.99999999999999989        0.0000000000000000
+   0.0000000000000000        0.0000000000000000      -0.99999999999999989
+
+Sym. op. :  3
+   0.0 180.0 240.0  1       - euler angles: a,b,c; iprop
+  0.50000000000000000      -0.86602540378443849        0.0000000000000000
+ -0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000      -0.99999999999999989
+
+Sym. op. :  4
+   0.0   0.0   0.0  1       - euler angles: a,b,c; iprop
+   1.0000000000000000        0.0000000000000000        0.0000000000000000
+   0.0000000000000000       0.99999999999999989        0.0000000000000000
+   0.0000000000000000        0.0000000000000000       0.99999999999999989
+
+Sym. op. :  5
+ 240.0   0.0   0.0  1       - euler angles: a,b,c; iprop
+ -0.50000000000000000      -0.86602540378443849        0.0000000000000000
+  0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000       0.99999999999999989
+
+Sym. op. :  6
+ 120.0   0.0   0.0  1       - euler angles: a,b,c; iprop
+ -0.50000000000000000       0.86602540378443849        0.0000000000000000
+ -0.86602540378443882      -0.49999999999999994        0.0000000000000000
+   0.0000000000000000        0.0000000000000000       0.99999999999999989
+ Global->local coordinates rotation matrices
+           1
+  -0.0000000000000000       0.50000000000000000       0.86602540378443860
+   0.0000000000000000      -0.86602540378443860       0.50000000000000000
+   1.0000000000000000        0.0000000000000000        0.0000000000000000
+   0.0  90.0 150.0  1       - euler angles: a,b,c; iprop
diff --git a/test/python/wien2k/Te_noso.ctqmcout.gz b/test/python/wien2k/Te_noso.ctqmcout.gz
new file mode 100644
index 0000000..026b2e0
Binary files /dev/null and b/test/python/wien2k/Te_noso.ctqmcout.gz differ
diff --git a/test/python/wien2k/Te_noso.indmftpr b/test/python/wien2k/Te_noso.indmftpr
new file mode 100644
index 0000000..7d09fed
--- /dev/null
+++ b/test/python/wien2k/Te_noso.indmftpr
@@ -0,0 +1,8 @@
+1
+3
+3
+cubic
+1 2 0 0
+0 0 0 0
+0
+-0.50 0.50
diff --git a/test/python/wien2k/Te_noso_paramag.ctqmcout.gz b/test/python/wien2k/Te_noso_paramag.ctqmcout.gz
new file mode 100644
index 0000000..b29732b
Binary files /dev/null and b/test/python/wien2k/Te_noso_paramag.ctqmcout.gz differ
diff --git a/test/python/wien2k/Te_noso_paramag.parproj.gz b/test/python/wien2k/Te_noso_paramag.parproj.gz
new file mode 100644
index 0000000..fdee260
Binary files /dev/null and b/test/python/wien2k/Te_noso_paramag.parproj.gz differ
diff --git a/test/python/wien2k/Te_noso_paramag.sympar.gz b/test/python/wien2k/Te_noso_paramag.sympar.gz
new file mode 100644
index 0000000..1ad62c5
Binary files /dev/null and b/test/python/wien2k/Te_noso_paramag.sympar.gz differ
diff --git a/test/python/wien2k/Te_noso_paramag.symqmc.gz b/test/python/wien2k/Te_noso_paramag.symqmc.gz
new file mode 100644
index 0000000..1fd4266
Binary files /dev/null and b/test/python/wien2k/Te_noso_paramag.symqmc.gz differ
diff --git a/test/python/wien2k/Te_sp.ctqmcout.gz b/test/python/wien2k/Te_sp.ctqmcout.gz
new file mode 100644
index 0000000..59b2c1e
Binary files /dev/null and b/test/python/wien2k/Te_sp.ctqmcout.gz differ
diff --git a/test/python/wien2k/Te_sp.indmftpr b/test/python/wien2k/Te_sp.indmftpr
new file mode 100644
index 0000000..a798da4
--- /dev/null
+++ b/test/python/wien2k/Te_sp.indmftpr
@@ -0,0 +1,8 @@
+1
+3
+3
+cubic
+2 2 0 0
+0 0 0 0
+1
+-0.80 0.50
diff --git a/test/python/wien2k/Te_sp.parproj.gz b/test/python/wien2k/Te_sp.parproj.gz
new file mode 100644
index 0000000..a2c298b
Binary files /dev/null and b/test/python/wien2k/Te_sp.parproj.gz differ
diff --git a/test/python/wien2k/Te_sp.symqmc.gz b/test/python/wien2k/Te_sp.symqmc.gz
new file mode 100644
index 0000000..cf32d40
Binary files /dev/null and b/test/python/wien2k/Te_sp.symqmc.gz differ
diff --git a/test/python/wien2k/Te_sp_mag.ctqmcout.gz b/test/python/wien2k/Te_sp_mag.ctqmcout.gz
new file mode 100644
index 0000000..ec6f8cd
Binary files /dev/null and b/test/python/wien2k/Te_sp_mag.ctqmcout.gz differ
diff --git a/test/python/wien2k/Te_sp_mag.parproj.gz b/test/python/wien2k/Te_sp_mag.parproj.gz
new file mode 100644
index 0000000..cce0a61
Binary files /dev/null and b/test/python/wien2k/Te_sp_mag.parproj.gz differ
diff --git a/test/python/wien2k/wien2k_dmftproj_common_python.py b/test/python/wien2k/wien2k_dmftproj_common_python.py
index 50c6770..fad15cc 100644
--- a/test/python/wien2k/wien2k_dmftproj_common_python.py
+++ b/test/python/wien2k/wien2k_dmftproj_common_python.py
@@ -46,10 +46,27 @@ def _reader(text):
 assert exact.shape == (5, 5)
 # unitary
 assert np.max(np.abs(exact @ np.conj(exact.T) - np.eye(5))) < 1e-12
+# The cubic p harmonics are orthonormal.
+p1 = dp.reptrans('cubic', 1)
+assert p1.shape == (3, 3)
+assert np.max(np.abs(p1 @ np.conj(p1.T) - np.eye(3))) < 1e-12
 assert abs(abs(exact[1, 0]) - 2 ** -0.5) < 1e-15
 assert dp.reptrans('complex', 2).shape == (5, 5)
 assert np.allclose(dp.reptrans('complex', 2), np.eye(5))
 
+# The tellurium representative has a nontrivial local frame. Its composed
+# orbital rotation must stay unitary and must include more than the atom-mapping
+# symmetry operation alone.
+_, te_ops, te_refs = dp.read_dmftsym(os.path.join(HERE, 'Te.dmftsym'),
+                                    rotloc=True)
+te_rot = dp.rotloc_orbital(1, te_refs[0], te_ops, 1, 1)
+np.testing.assert_allclose(te_rot @ te_rot.conj().T, np.eye(3),
+                           atol=1e-12, rtol=0)
+mapping = next(op for op in te_ops if op['perm'][0] == 1)
+mapping_only = dp.dmat(1, mapping['a'], mapping['b'], mapping['g'],
+                       float(mapping['iprop']))
+assert np.max(np.abs(te_rot - mapping_only)) > 0.1
+
 
 # --- Wigner D and orbital time reversal --------------------------------------
 
diff --git a/test/python/wien2k/wien2k_symqmc_noncentrosym.py b/test/python/wien2k/wien2k_symqmc_noncentrosym.py
new file mode 100644
index 0000000..e7a3ecf
--- /dev/null
+++ b/test/python/wien2k/wien2k_symqmc_noncentrosym.py
@@ -0,0 +1,38 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Check unitarity of the cubic p transform and spinor symmetry matrices for
+# non-centrosymmetric tellurium (P3_1 21, point group 32). The comparison with
+# double-precision native dmftproj output is covered by wien2k_te_python.py.
+
+import os
+
+import numpy as np
+
+from triqs_dftkit.wien2k import symqmc
+from triqs_dftkit.wien2k._dmftproj import reptrans
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+P1 = reptrans('cubic', 1)
+assert np.max(np.abs(P1 @ np.conj(P1.T) - np.eye(3))) < 1e-12, 'cubic l=1 not unitary'
+
+_, ops = symqmc._read_dmftsym(os.path.join(HERE, 'Te.dmftsym'))
+shells, so = symqmc._read_correlated_shells(os.path.join(HERE, 'Te.indmftpr'),
+                                        os.path.join(HERE, 'Te.struct'))
+
+
+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
+
+
+for op in ops:
+    for sh in shells:
+        mat = symqmc._shell_matrix(op, sh, _timeinv(op, so), bool(so))
+        dev = np.max(np.abs(mat @ np.conj(mat.T) - np.eye(mat.shape[0])))
+        assert dev < 1e-7, ('not unitary', dev)
+
+print('wien2k_symqmc_noncentrosym: ok')
diff --git a/test/python/wien2k/wien2k_te_python.py b/test/python/wien2k/wien2k_te_python.py
new file mode 100644
index 0000000..94711f8
--- /dev/null
+++ b/test/python/wien2k/wien2k_te_python.py
@@ -0,0 +1,84 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.* generators on tellurium (P3_121), a
+# non-centrosymmetric cell with three equivalent atoms. The local frame of each
+# atom is the representative rotloc composed with the symmetry operation that
+# maps the representative onto it, which the centrosymmetric fixtures never
+# exercise. Inputs: Te.almblm{up,dn}.gz (one k-point, see Te.README.md) with
+# Te.struct and Te.dmftsym. References are native dmftproj output, double
+# precision, with the full-precision cubic templates.
+#
+#   Te              p correlated, SO (dmftproj -sp -so), window -0.5 0.5 Ry
+#   Te_sp           s and p correlated, SO, window -0.8 0.5 Ry
+#   Te_sp_mag       as Te_sp, with Te_mag.dmftsym: the operations of Te.dmftsym
+#                   reordered so that a two-fold axis, a magnetic operation under
+#                   SO, maps the representative onto each atom
+#   Te_noso         s included, p correlated, no SO (dmftproj -sp)
+#   Te_noso_paramag as Te_noso from Te.almblmup alone (dmftproj without flags)
+
+import gzip
+import os
+import shutil
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import ctqmcout, symqmc, sympar, parproj
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+WRITERS = dict(ctqmcout=ctqmcout.write_ctqmcout, symqmc=symqmc.write_symqmc,
+               sympar=sympar.write_sympar, parproj=parproj.write_parproj)
+
+
+def _gunzip(src, dst):
+    with gzip.open(src, 'rb') as fi, open(dst, 'wb') as fo:
+        shutil.copyfileobj(fi, fo)
+
+
+def _compare(got_path, ref_path):
+    got = open(got_path).read().split()
+    with gzip.open(ref_path, 'rt') as fh:
+        ref = fh.read().split()
+    assert len(got) == len(ref), (got_path, len(got), len(ref))
+    floats = []
+    for a, b in zip(got, ref):
+        if b.lstrip('-').isdigit():
+            assert int(a) == int(b), (got_path, a, b)
+        else:
+            floats.append((float(a), float(b)))
+    actual, expected = np.array(floats).T
+    np.testing.assert_allclose(actual, expected, atol=1e-11, rtol=0,
+                               equal_nan=False, err_msg=got_path)
+
+
+def _check(reference, indmftpr, exts, dmftsym='Te', paramagnetic=False):
+    tmp = tempfile.mkdtemp()
+    case = os.path.join(tmp, 'case')
+    try:
+        shutil.copy(os.path.join(HERE, f'{indmftpr}.indmftpr'), case + '.indmftpr')
+        shutil.copy(os.path.join(HERE, 'Te.struct'), case + '.struct')
+        shutil.copy(os.path.join(HERE, f'{dmftsym}.dmftsym'), case + '.dmftsym')
+        if paramagnetic:
+            _gunzip(os.path.join(HERE, 'Te.almblmup.gz'), case + '.almblm')
+        else:
+            for ext in ('almblmup', 'almblmdn'):
+                _gunzip(os.path.join(HERE, f'Te.{ext}.gz'), f'{case}.{ext}')
+        for ext in exts:
+            WRITERS[ext](case)
+            _compare(f'{case}.{ext}', os.path.join(HERE, f'{reference}.{ext}.gz'))
+    finally:
+        shutil.rmtree(tmp)
+
+
+_check('Te', 'Te', ('ctqmcout', 'symqmc', 'parproj'))
+_check('Te_sp', 'Te_sp', ('ctqmcout', 'symqmc', 'parproj'))
+_check('Te_sp_mag', 'Te_sp', ('ctqmcout', 'parproj'), dmftsym='Te_mag')
+_check('Te_noso', 'Te_noso', ('ctqmcout',))
+_check('Te_noso_paramag', 'Te_noso', ('ctqmcout', 'symqmc', 'sympar', 'parproj'),
+       paramagnetic=True)
+
+print('wien2k_te_python: ok')