diff --git a/.repo-hygiene-allow b/.repo-hygiene-allow
deleted file mode 100644
index 86f9197..0000000
--- a/.repo-hygiene-allow
+++ /dev/null
@@ -1 +0,0 @@
-test/python/wien2k/CaOs2_pcubic.symqmc.gz	Native dmftproj cubic-p regression output from the published CaOs2 inputs retained in this repository.
diff --git a/python/triqs_dftkit/wien2k/_dmftproj.py b/python/triqs_dftkit/wien2k/_dmftproj.py
index c4db637..5f373bf 100644
--- a/python/triqs_dftkit/wien2k/_dmftproj.py
+++ b/python/triqs_dftkit/wien2k/_dmftproj.py
@@ -175,6 +175,10 @@ def small_d(l, m, n, b):
 
 
 def dmat(l, a, b, c, det):
+    # det is the improper-rotation parity: pass op['iprop'] (+-1), not
+    # np.linalg.det(krotm). krotm in case.dmftsym is the proper part and its
+    # determinant is +1 even for improper ops, so it drops the (-1)^l parity
+    # factor. Harmless for even l, wrong sign for odd l (p, f).
     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):
@@ -198,6 +202,16 @@ def timeinv_orbital(l, mat):
     return tmat(l) @ np.conj(mat)
 
 
+def paramag_time_op(shell):
+    """Paramagnetic time-reversal operator P T P^T of one non-mixing shell in
+    its basis (outputqmc.f, .not.ifSP); the 1x1 identity for an s shell."""
+    l = shell['l']
+    if l == 0:
+        return np.ones((1, 1), dtype=complex)
+    P = shell['P']
+    return P @ tmat(l) @ P.T
+
+
 def mixing_timeinv_op(l, P):
     """The spinor time-reversal operator -i sigma_y (x) T in the mixing basis:
     tinv_{new} = P tinv_{lm} P^T (timeinv.f, ifmixing branch). Returned as the
@@ -304,9 +318,17 @@ def _nonmixing_matrix(op, shell, ti):
     return mat
 
 
-def shell_matrix(op, shell, ti):
-    """Spinor symmetry matrix for one shell under one operation."""
+def shell_matrix(op, shell, ti, ifSO=True):
+    """Spinor or orbital symmetry matrix for one shell under one operation."""
+    if shell['mixing'] and not ifSO:
+        raise NotImplementedError('spin-mixing bases without spin-orbit coupling are not implemented')
     l = shell['l']
+    if not ifSO:
+        if l == 0:
+            return np.array([[1.0 + 0j]], dtype=complex)
+        P = shell['P']
+        rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+        return P @ rotl @ np.conj(P.T)
     if l == 0:
         return _l0_matrix(op, ti)
     if shell['mixing']:
@@ -614,10 +636,10 @@ def require_identity_struct(case):
                 raise NotImplementedError('non-identity struct local rotations are not implemented yet')
 
 
-def validate_so_windows(spins, info, windows, below=None):
+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:
-        raise ValueError('spin-orbit projectors require both up and down almblm files')
+        raise ValueError('two-spin projectors require both up and down almblm files')
     bw = band_index_window(info, spins) if info['proj_mode'] else None
     for spin in spins[1:]:
         if spin['nk'] != spins[0]['nk']:
@@ -626,11 +648,11 @@ def validate_so_windows(spins, info, windows, below=None):
             selected = (select_window(kp['nbmin'], kp['nbmax'], kp['eband'], info['e_bot'], info['e_top'])
                         if bw is None else select_band_window(kp['nbmin'], kp['nbmax'], *bw))
             if selected != windows[ik]:
-                raise ValueError('spin-orbit run requires identical up/dn band windows')
+                raise ValueError('projector run requires identical up/dn band windows')
             if below is not None:
                 lower = (select_window(kp['nbmin'], kp['nbmax'], kp['eband'], -1e6, info['e_bot'])
                          if bw is None else select_band_window(kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
                 if lower != below[ik]:
-                    raise ValueError('spin-orbit run requires identical up/dn below-window bands')
+                    raise ValueError('projector run requires identical up/dn below-window bands')
     if not all(incl and lo <= hi for incl, lo, hi in windows):
         raise ValueError('projectors require a nonempty included window at every k-point')
diff --git a/python/triqs_dftkit/wien2k/ctqmcout.py b/python/triqs_dftkit/wien2k/ctqmcout.py
index aac711b..22f6e16 100644
--- a/python/triqs_dftkit/wien2k/ctqmcout.py
+++ b/python/triqs_dftkit/wien2k/ctqmcout.py
@@ -60,7 +60,7 @@
 from scipy.linalg.lapack import zheev
 from scipy.linalg.blas import zgemm
 
-from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_so_windows, band_index_window, dmat, fmt_scalar, read_almblm,
+from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, 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,
                         tmat, write_row)
@@ -102,18 +102,23 @@ def _build_crorbs(info):
 
 # --- Rloc spinor representation (set_rotloc.f / setsym.f) ---------------------
 
-def _rloc_rotrep(op, l, transmat):
-    """rotloc%rotrep(l)%mat, the 2(2l+1) spinor rotation in the new basis for a
-    non-mixing SO shell, with identity struct local rotation (rotloc_ref
-    Euler = 0). rotloc%rotl = blkdiag(ephase*D, conj(ephase)*D) with
-    D = D(R[isym])_{lm}, ephase = exp(i*phase/2); rotrep = S rotl S^H,
-    S = blkdiag(transmat, transmat). Returns (rotrep, timeinv)."""
+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
+    (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:
+        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)
-    D = dmat(l, a, b, g, float(iprop))
     if timeinv:
         D = tmat(l) @ np.conj(D)       # setsym.f:320-326 orbital time reversal
     ephase = np.exp(1j * phase / 2)
@@ -143,8 +148,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')
-    if not info['so']:
-        raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
     require_identity_struct(case)
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
 
@@ -194,7 +197,10 @@ def write_ctqmcout(case):
             win_below.append(select_band_window(
                 kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
 
-    validate_so_windows(spins, info, windows, win_below)
+    if ifSO or len(spins) == 2:
+        validate_spin_windows(spins, info, windows, win_below)
+    elif not all(window[0] and window[1] <= window[2] for window in windows):
+        raise ValueError('projectors require a nonempty included window at every k-point')
 
     # qbbot: point integration over bands below e_bot, is=1 only under SO
     qbbot = 0.0
@@ -215,6 +221,8 @@ def write_ctqmcout(case):
                 info['sorts'][cr['sort'] - 1]['sourcefile'], cr['l'])
         else:
             transmats[key], mixing[key] = reptrans(cr['basis'], cr['l']), False
+    if not ifSO and any(mixing.values()):
+        raise NotImplementedError('spin-mixing bases require spin-orbit coupling')
 
     # rotloc Euler angles per crorb: the symmetry op mapping the representative
     # atom of the sort onto this atom (set_rotloc.f). With identity struct
@@ -270,10 +278,17 @@ def write_ctqmcout(case):
                 for ispin in range(ns):
                     mat_rep[(icr, ik, ispin)] = transmat @ spin_blocks[ispin]
 
-    # Loewdin: per k stack all crorb rows. A mixing crorb contributes its full
-    # 2(2l+1) block once (is=1 only); a non-mixing one the two (2l+1) spin
-    # blocks (orthogonal_wannier_SO ndim layout).
-    for ik in range(nk):
+    # Loewdin orthonormalization. Under SO (orthogonal_wannier_SO) the up+dn
+    # blocks of each crorb are stacked together and a single overlap is
+    # orthonormalized per k. Without SO (orthogonal_wannier) spin is a good
+    # quantum number, so each spin is orthonormalized independently over the
+    # stack of (2l+1)-wide crorb blocks of that spin.
+    if ifSO:
+        ortho_groups = [[(ik, None)] for ik in range(nk)]
+    else:
+        ortho_groups = [[(ik, ispin)] for ik in range(nk) for ispin in range(ns)]
+    for group in ortho_groups:
+        (ik, gspin) = group[0]
         incl, nb_bot, nb_top = windows[ik]
         if not incl:
             continue
@@ -281,15 +296,18 @@ def write_ctqmcout(case):
         layout = []  # (key, nrows)
         for icr, cr in enumerate(crorbs):
             l = cr['l']
-            if mixing[(l, cr['sort'])]:
+            if ifSO and mixing[(l, cr['sort'])]:
                 blocks.append(mat_rep[(icr, ik)])
                 layout.append(((icr, ik), 2 * (2 * l + 1)))
-            else:
+            elif ifSO:
                 for ispin in range(ns):
                     blocks.append(mat_rep[(icr, ik, ispin)])
                     layout.append(((icr, ik, ispin), 2 * l + 1))
+            else:
+                blocks.append(mat_rep[(icr, ik, gspin)])
+                layout.append(((icr, ik, gspin), 2 * l + 1))
         D = np.vstack(blocks)                 # ndim x nbnd
-        # match the Fortran's exact BLAS calls (orthogonal_wannier_SO): the
+        # match the Fortran's exact BLAS calls (orthogonal_wannier[_SO]): the
         # near-singular O^{-1/2} amplifies any last-bit difference, so use the
         # identical ZGEMM trans flags rather than numpy's conj-transpose copies.
         O = zgemm(1.0, D, D, trans_b=2)        # ZGEMM('N','C'): D @ D^H
@@ -311,7 +329,7 @@ def write_ctqmcout(case):
                 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))
+            rloc_blocks.append(_rloc_rotrep(op, l, transmat, ifSO))
 
     # ---- write ----
     with open(case + '.ctqmcout', 'w') as f:
@@ -332,15 +350,18 @@ def write_ctqmcout(case):
             l = cr['l']
             size = 2 * (2 * l + 1) if ifSO else 2 * l + 1
             f.write('%6d %6d %6d %6d %6d %6d \n' %
-                    (cr['atom'], cr['sort'], l, size, 1, 1))
+                    (cr['atom'], cr['sort'], l, size, 1 if ifSO else 0, 1))
 
-        # Rloc block per crorb (non-mixing SP+SO whole shell)
+        # Rloc block per crorb. SO writes the 2*(2l+1) spinor whole shell;
+        # non-SO the bare (2l+1) shell. The time-reversal flag follows under
+        # SO always, under non-SO only when ifSP (outputqmc.f:347-351).
         for rotrep, timeinv in rloc_blocks:
             for m in range(rotrep.shape[0]):
                 write_row(f, rotrep[m, :].real)
             for m in range(rotrep.shape[0]):
                 write_row(f, rotrep[m, :].imag)
-            f.write('%6d\n' % (1 if timeinv else 0))
+            if ifSO or ifSP:
+                f.write('%6d\n' % (1 if timeinv else 0))
 
         # complex-harmonics -> basis transform block (crorb%first only). Mixing
         # writes the full 2(2l+1) transmat; non-mixing the spin block-diagonal.
@@ -352,14 +373,17 @@ def write_ctqmcout(case):
             d = 2 * l + 1
             if mixing[(l, cr['sort'])]:
                 spinrot = transmat
-            else:
+            elif ifSO:
                 spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
                 spinrot[:d, :d] = transmat
                 spinrot[d:, d:] = transmat
-            f.write('%6d %6d \n' % (1, 2 * d))
-            for m in range(2 * d):
+            else:
+                spinrot = transmat       # bare (2l+1) transform, non-SO
+            dim = spinrot.shape[0]
+            f.write('%6d %6d \n' % (1, dim))
+            for m in range(dim):
                 write_row(f, spinrot[m, :].real)
-            for m in range(2 * d):
+            for m in range(dim):
                 write_row(f, spinrot[m, :].imag)
 
         # number of bands per k (skip is=2 under SO)
diff --git a/python/triqs_dftkit/wien2k/outband.py b/python/triqs_dftkit/wien2k/outband.py
index e16f302..65e5787 100644
--- a/python/triqs_dftkit/wien2k/outband.py
+++ b/python/triqs_dftkit/wien2k/outband.py
@@ -51,7 +51,7 @@
 import numpy as np
 from scipy.linalg.blas import zgemm
 
-from ._dmftproj import (require_identity_struct, validate_so_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
+from ._dmftproj import (require_identity_struct, validate_spin_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
                         read_indmftpr, reptrans, select_window, to_float,
                         write_row)
 from .ctqmcout import _build_crorbs, _sqrt_inv
@@ -142,7 +142,7 @@ def write_outband(case):
         windows.append(select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
                                      info['e_bot'], info['e_top']))
 
-    validate_so_windows(spins, info, windows)
+    validate_spin_windows(spins, info, windows)
 
     transmats = {}
     mixing = {}
diff --git a/python/triqs_dftkit/wien2k/parproj.py b/python/triqs_dftkit/wien2k/parproj.py
index fd74118..22488a0 100644
--- a/python/triqs_dftkit/wien2k/parproj.py
+++ b/python/triqs_dftkit/wien2k/parproj.py
@@ -57,15 +57,15 @@
 import os
 import numpy as np
 
-from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_so_windows, band_index_window, dmat, mixing_rotrep, read_almblm, read_dmftsym,
+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)
 
 
-def _op_timeinv(op):
+def _op_timeinv(op, ifSO=True):
     det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
             - op['krotm'][0, 1] * op['krotm'][1, 0])
-    return bool(det2 < 0.0)
+    return bool(ifSO and det2 < 0.0)
 
 
 def _sqrtm_real_sym(O):
@@ -97,15 +97,15 @@ def _build_s12(l, isrt, n, spin):
 
 # --- srot rotrep (setsym.f) --------------------------------------------------
 
-def _srot_rotrep_nonmixing(op, l, transmat):
+def _srot_rotrep_nonmixing(op, l, transmat, ifSO=True):
     """srot(isym)%rotrep(l,isrt)%mat, the (2l+1) up/up block of the symmetry op
     in the new basis: transmat . D(R)_{lm} . transmat^H, with the orbital
     time-reversal operator applied for the magnetic (timeinv) operations.
 
     Returns (rotrep[2l+1,2l+1], timeinv, phase)."""
     a, b, g, iprop = op['a'], op['b'], op['g'], op['iprop']
-    timeinv = _op_timeinv(op)
-    phase = (g - a) if timeinv else (a + g)
+    timeinv = _op_timeinv(op, ifSO)
+    phase = ((g - a) if timeinv else (a + g)) if ifSO else 0.0
     rotl = dmat(l, a, b, g, float(iprop))
     rotrep = transmat @ rotl @ np.conj(transmat.T)
     if timeinv:
@@ -117,6 +117,18 @@ def _srot_rotrep_nonmixing(op, l, transmat):
 
 # --- rotloc rotrep (set_rotloc.f) under SP+SO, non-mixing --------------------
 
+def _rotloc_rotrep_nonso(orb, ops, info, 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']))
+    return transmat @ D @ np.conj(transmat.T), False
+
+
 def _rotloc_rotrep_so(orb, ops, info, ref, transmat, mixing):
     """rotloc(iatom)%rotrep(l)%mat, the full 2*(2l+1) Rloc spinor rotation in
     the new basis under SP+SO (set_rotloc.f), plus timeinv flag.
@@ -270,12 +282,14 @@ def _rotdens_densmat_mixing(blk, rotrep_loc, timeinv):
 # --- density matrix (density.f Theta path, non-mixing SP+SO) -----------------
 
 def _orbital_densmat_blocks(orb, info, spins, windows, matn_reps):
-    """The nsp=4 raw density-matrix blocks (up/up, dn/dn, up/dn, dn/up) for one
-    orbital, point-integrated with the geometric k-weight (density.f:889-907,
-    tetr=.FALSE. path). The window is the bands below e_bot (dmftproj.f:798)."""
+    """Raw density blocks below the window, integrated with geometric k-weights.
+
+    Spin-orbit coupling requires all four spin blocks; otherwise each input
+    spin contributes one diagonal block (density.f:67-72,889-907)."""
     l = orb['l']
     d = 2 * l + 1
-    blocks = [np.zeros((d, d), dtype=complex) for _ in range(4)]
+    nblocks = 4 if info['so'] else len(spins)
+    blocks = [np.zeros((d, d), dtype=complex) for _ in range(nblocks)]
     for ik in range(info['nk']):
         incl, nb_bot, nb_top = windows[ik]
         if not incl:
@@ -284,7 +298,7 @@ def _orbital_densmat_blocks(orb, info, spins, windows, matn_reps):
         off_bot = nb_bot - kp0['nbmin']
         off_top = nb_top - kp0['nbmin']
         n = spins[0]['nLO'][(l, orb['sort'])] + 2
-        for iss in range(1, 5):
+        for iss in range(1, nblocks + 1):
             if iss <= 2:
                 is_, is1 = iss, iss
             else:
@@ -344,6 +358,44 @@ def _symmetrize_densmat(orbs, ops, nsym, srot_rotreps, dens_raw):
     return out
 
 
+def _symmetrize_densmat_nonso(orbs, info, ops, nsym, srot_rotreps, dens_raw):
+    """symmetrize_mat for the non-SO non-mixing path (symmetrize_mat.f:196-262,
+    nsp=1 or 2). Each input spin has one diagonal block; srot%timeinv is
+    always false so the ephase phase factors are 1. Each block is summed over
+    symmetry ops as rotrep block rotrep^H, equivalent atoms scattering via the
+    permutation, divided by nsym. Returns the averaged blocks for each atom."""
+    out = {}
+    isort_groups = {}
+    for o in orbs:
+        isort_groups.setdefault((o['sort'], o['l']), []).append(o)
+    for (isrt, l), group in isort_groups.items():
+        d = 2 * l + 1
+        mult = len(group)
+        nblocks = len(dens_raw[(group[0]['atom'], l)])
+        sym = [[np.zeros((d, d), dtype=complex) for _ in range(nblocks)]
+               for _ in range(mult)]
+        for imult in range(mult):
+            iatom = group[imult]['atom']
+            for isym in range(nsym):
+                op = ops[isym]
+                rotrep, _, _ = srot_rotreps[(l, isrt, isym)]
+                jorb = op['perm'][iatom - 1] - iatom + imult
+                for iss in range(nblocks):
+                    blk = dens_raw[(iatom, l)][iss]
+                    sym[jorb][iss] += rotrep @ (blk @ np.conj(rotrep.T))
+        for imult in range(mult):
+            iatom = group[imult]['atom']
+            out[(iatom, l)] = [sym[imult][k] / nsym for k in range(nblocks)]
+    return out
+
+
+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
+    (2l+1) rotloc rotrep (timeinv always false). Returns one block per input spin."""
+    return [np.conj(rotrep_loc.T) @ (blk @ rotrep_loc) for blk in blocks]
+
+
 def _rotdens_densmat(orb, blocks, rotrep_loc, timeinv):
     """rotdens_mat for non-mixing SP+SO (rot_dens.f:151-193): assemble the four
     blocks into a 2*(2l+1) matrix, apply inverse(Rloc) D Rloc, return it as the
@@ -370,14 +422,15 @@ 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')
-    if not info['so']:
-        raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
     require_identity_struct(case)
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
 
     spin_files = [case + '.almblmup', case + '.almblmdn']
     if not all(os.path.exists(p) for p in spin_files):
-        raise ValueError('spin-orbit parproj requires both up and down almblm files')
+        spin_files = [case + '.almblm']
+    ifSO = bool(info['so'])
+    ifSP = len(spin_files) == 2
+    ns = len(spin_files)
 
     spins = [read_almblm(p, info, projectors=True) for p in spin_files]
     nk = spins[0]['nk']
@@ -404,7 +457,10 @@ def write_parproj(case):
             windows.append(select_band_window(kp['nbmin'], kp['nbmax'], *bw))
             below_windows.append(select_band_window(kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
 
-    validate_so_windows(spins, info, windows, below_windows)
+    if ifSO or ifSP:
+        validate_spin_windows(spins, info, windows, below_windows)
+    elif not all(window[0] and window[1] <= window[2] for window in windows):
+        raise ValueError('projectors require a nonempty included window at every k-point')
 
     transmats = {}
     mixing = {}
@@ -415,6 +471,8 @@ def write_parproj(case):
                 info['sorts'][o['sort'] - 1]['sourcefile'], o['l'])
         else:
             transmats[key], mixing[key] = reptrans(o['basis'], o['l']), False
+    if not ifSO and any(mixing.values()):
+        raise NotImplementedError('spin-mixing bases require spin-orbit coupling')
 
     # rot_projectmat local rotation: the op mapping the representative atom of
     # the sort onto this atom (set_projections.f via rot_projectmat).
@@ -448,12 +506,18 @@ def write_parproj(case):
     rotloc_rotrep = {}
     for o in orbs:
         transmat = transmats[(o['l'], o['sort'])]
-        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'])])
+        if not ifSO:
+            rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_nonso(
+                o, ops, info, 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'])])
 
     # ---- density matrices: raw -> symmetrize -> rotdens ----
-    # Non-mixing sorts use the 4-block path; mixing sorts the single-block path.
+    # Non-mixing sorts use the block path; mixing sorts the single-block path.
+    # Under SO each non-mixing density is the 2*(2l+1) 4-block matrix; without
+    # SO there is one (2l+1) diagonal block per input spin.
     nonmix_orbs = [o for o in orbs if not mixing[(o['l'], o['sort'])]]
     densprint = {}
 
@@ -463,16 +527,26 @@ def write_parproj(case):
             l, isrt = o['l'], o['sort']
             for isym in range(nsym):
                 srot_rotreps[(l, isrt, isym)] = _srot_rotrep_nonmixing(
-                    ops[isym], l, transmats[(l, isrt)])
+                    ops[isym], l, transmats[(l, isrt)], ifSO)
         dens_raw = {(o['atom'], o['l']): _orbital_densmat_blocks(
             o, info, spins, below_windows, matn_reps_full)
             for o in nonmix_orbs}
-        dens_sym = _symmetrize_densmat(nonmix_orbs, ops, nsym,
-                                       srot_rotreps, dens_raw)
-        for o in nonmix_orbs:
-            rotrep_loc, timeinv = rotloc_rotrep[(o['atom'], o['l'])]
-            densprint[(o['atom'], o['l'])] = _rotdens_densmat(
-                o, dens_sym[(o['atom'], o['l'])], rotrep_loc, timeinv)
+        if ifSO:
+            dens_sym = _symmetrize_densmat(nonmix_orbs, ops, nsym,
+                                           srot_rotreps, dens_raw)
+            for o in nonmix_orbs:
+                key = (o['atom'], o['l'])
+                rotrep_loc, timeinv = rotloc_rotrep[key]
+                densprint[key] = _rotdens_densmat(
+                    o, dens_sym[key], rotrep_loc, timeinv)
+        else:
+            dens_sym = _symmetrize_densmat_nonso(nonmix_orbs, info, ops, nsym,
+                                                 srot_rotreps, dens_raw)
+            for o in nonmix_orbs:
+                key = (o['atom'], o['l'])
+                rotrep_loc, _ = rotloc_rotrep[key]
+                densprint[key] = _rotdens_densmat_nonso(
+                    dens_sym[key], rotrep_loc)
 
     mixing_groups = {}
     for o in orbs:
@@ -505,7 +579,7 @@ def write_parproj(case):
             ismix = mixing[(l, o['sort'])]
 
             # (A) Theta projector (outputqmc.f:935-973). Mixing writes the full
-            # 2(2l+1) block from is=1 only; non-mixing the two (2l+1) spin blocks.
+            # 2(2l+1) block from is=1 only; non-mixing one (2l+1) block per input spin.
             for ik in range(nk):
                 incl, nb_bot, nb_top = windows[ik]
                 for ir in range(n):
@@ -516,26 +590,38 @@ def write_parproj(case):
                         for m in range(2 * (2 * l + 1)):
                             write_row(f, P[m, :].imag)
                         continue
-                    for ispin in range(2):
+                    for ispin in range(ns):
                         P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
                         for mi in range(2 * l + 1):
                             write_row(f, P[mi, :].real)
-                    for ispin in range(2):
+                    for ispin in range(ns):
                         P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
                         for mi in range(2 * l + 1):
                             write_row(f, P[mi, :].imag)
 
-            # (B) density matrix, SP+SO 2*(2l+1) (outputqmc.f:1012-1066).
+            # (B) density matrix (outputqmc.f:1012-1066). SO writes the single
+            # 2*(2l+1) block; mixing too; non-SO one (2l+1) block per input spin
+            # (each real then imag, per spin, outputqmc.f:1052-1064).
             dp = densprint[(atom, l)]
-            for m in range(2 * (2 * l + 1)):
-                write_row(f, dp[m, :].real)
-            for m in range(2 * (2 * l + 1)):
-                write_row(f, dp[m, :].imag)
-
-            # (C) Rloc rotrep, SP+SO 2*(2l+1) (outputqmc.f:1130-1177).
+            if ifSO or ismix:
+                for m in range(2 * (2 * l + 1)):
+                    write_row(f, dp[m, :].real)
+                for m in range(2 * (2 * l + 1)):
+                    write_row(f, dp[m, :].imag)
+            else:
+                for ispin in range(ns):
+                    for m in range(2 * l + 1):
+                        write_row(f, dp[ispin][m, :].real)
+                for ispin in range(ns):
+                    for m in range(2 * l + 1):
+                        write_row(f, dp[ispin][m, :].imag)
+
+            # (C) Rloc rotrep (outputqmc.f:1130-1177). SO/mixing 2*(2l+1);
+            # non-SO bare (2l+1). The time-reversal flag follows under ifSP.
             rotrep_loc, timeinv = rotloc_rotrep[(atom, l)]
-            for m in range(2 * (2 * l + 1)):
+            for m in range(rotrep_loc.shape[0]):
                 write_row(f, rotrep_loc[m, :].real)
-            for m in range(2 * (2 * l + 1)):
+            for m in range(rotrep_loc.shape[0]):
                 write_row(f, rotrep_loc[m, :].imag)
-            f.write('%6d\n' % (1 if timeinv else 0))
+            if ifSP:
+                f.write('%6d\n' % (1 if timeinv else 0))
diff --git a/python/triqs_dftkit/wien2k/sympar.py b/python/triqs_dftkit/wien2k/sympar.py
index 493fb2d..0a4b629 100644
--- a/python/triqs_dftkit/wien2k/sympar.py
+++ b/python/triqs_dftkit/wien2k/sympar.py
@@ -39,14 +39,15 @@
 no crorb%ifSOat/correp metadata. It is purely group-theoretic: header, perms,
 timeflag line, and the representation matrices.
 
-Only the spin-orbit case is implemented (indmftpr SO flag=1, which implies
-ifSP): 2*(2l+1)-wide spinor matrices with the timeflag line. The non-SO layouts
-(scalar l=0 entries, paramagnetic time-reversal tail) raise
-NotImplementedError. The reference case is CaOs2 (SP+SO).
+Spin-orbit runs write 2*(2l+1)-wide spinor matrices. Without spin-orbit
+coupling, the matrices have width 2l+1, with scalar s-shell entries and a
+paramagnetic time-reversal tail when the input is not spin-polarized.
 """
 
-from ._dmftproj import (read_dmftsym, read_fromfile, read_indmftpr, reptrans,
-                        shell_matrix, write_row)
+import os
+
+from ._dmftproj import (paramag_time_op, read_dmftsym, read_fromfile, read_indmftpr,
+                        reptrans, shell_matrix, write_row)
 
 _shell_matrix = shell_matrix
 
@@ -85,28 +86,36 @@ def _write_matrix(f, mat):
 def write_sympar(case):
     """Write <case>.sympar from <case>.dmftsym and <case>.indmftpr.
 
-    Spin-orbit only (indmftpr SO flag=1, implying ifSP): writes the timeflag
-    line and the 2*(2l+1)-wide spinor matrices."""
+    Detects ifSP/ifSO from the presence of <case>.almblm{up,dn} and the indmftpr
+    SO flag, exactly as oubwin/symqmc/ctqmcout. For SP+SO the file carries the
+    timeflag line and 2*(2l+1)-wide block-diag spinor matrices; the paramagnetic
+    time-reversal tail is written only when .not.ifSP."""
     info = read_indmftpr(case + '.indmftpr')
-    if not info['so']:
-        raise NotImplementedError('non-spin-orbit sympar is not implemented yet')
-    shells = _included_shells(info)
+    shells, so = _included_shells(info), info['so']
     nsym, ops = read_dmftsym(case + '.dmftsym')
     natom = len(ops[0]['perm'])
 
+    ifSP = os.path.exists(case + '.almblmup') and os.path.exists(case + '.almblmdn')
+    ifSO = bool(so)
+    ifSP = ifSP or ifSO                      # ifSO => ifSP
+
     timeflag = []
     for op in ops:
         det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
                 - op['krotm'][0, 1] * op['krotm'][1, 0])
-        timeflag.append(1 if det2 < 0.0 else 0)
+        timeflag.append(1 if (ifSO and det2 < 0.0) else 0)
 
     with open(case + '.sympar', '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')
-        f.write(''.join('%6d ' % t for t in timeflag) + '\n')
+        if ifSP:
+            f.write(''.join('%6d ' % t for t in timeflag) + '\n')
 
         for isym, op in enumerate(ops):
-            ti = bool(timeflag[isym])
             for sh in shells:
-                _write_matrix(f, _shell_matrix(op, sh, ti))
+                _write_matrix(f, _shell_matrix(op, sh, bool(timeflag[isym]), ifSO))
+
+        if not ifSP:
+            for sh in shells:
+                _write_matrix(f, paramag_time_op(sh))
diff --git a/python/triqs_dftkit/wien2k/symqmc.py b/python/triqs_dftkit/wien2k/symqmc.py
index 524d9e0..f8d4336 100644
--- a/python/triqs_dftkit/wien2k/symqmc.py
+++ b/python/triqs_dftkit/wien2k/symqmc.py
@@ -46,10 +46,10 @@
 through single-precision CMPLX differs from this output at about 1e-7.
 """
 
-import numpy as np
+import os
 
-from ._dmftproj import (read_dmftsym, read_fromfile, read_indmftpr, reptrans,
-                        shell_matrix)
+from ._dmftproj import (paramag_time_op, read_dmftsym, read_fromfile, read_indmftpr,
+                        reptrans, shell_matrix)
 
 _shell_matrix = shell_matrix
 
@@ -86,26 +86,39 @@ def _read_correlated_shells(indmftpr, struct):
 
 
 def write_symqmc(case):
-    """Write <case>.symqmc from <case>.dmftsym, <case>.indmftpr, <case>.struct."""
+    """Write <case>.symqmc from <case>.dmftsym, <case>.indmftpr, <case>.struct.
+
+    ifSO comes from the indmftpr SO flag; ifSP from the presence of
+    <case>.almblm{up,dn} (ifSO => ifSP). With SO the matrices are the 2*(2l+1)
+    block-diag spinor representation and a timeflag line is written; without SO
+    they are the bare (2l+1) representation, with a timeflag line (all zero)
+    only when ifSP. The paramagnetic time-reversal tail follows when .not.ifSP."""
     nsym, ops = read_dmftsym(case + '.dmftsym')
     info = read_indmftpr(case + '.indmftpr')
     shells, so = _correlated_shells(info), info['so']
     natom = len(ops[0]['perm'])
 
+    ifSO = bool(so)
+    ifSP = (os.path.exists(case + '.almblmup')
+            and os.path.exists(case + '.almblmdn')) or ifSO
+
     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)
+        timeinv.append(1 if (ifSO 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:
+        if ifSP:
             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])))
+                f.write(_format_matrix(_shell_matrix(op, sh, timeinv[isym], ifSO)))
+        if not ifSP:
+            for sh in shells:
+                f.write(_format_matrix(paramag_time_op(sh)))
 
 
 def _format_matrix(mat):
diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt
index f94fffa..a562874 100644
--- a/test/python/wien2k/CMakeLists.txt
+++ b/test/python/wien2k/CMakeLists.txt
@@ -91,3 +91,37 @@ add_test(NAME Py_wien2k_outband_python
          WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 set_property(TEST Py_wien2k_outband_python APPEND PROPERTY ENVIRONMENT
              PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Pure-Python generators on the non-spin-orbit path (ifSO=0, the common DMFT
+# case): dmftproj writes (2l+1) matrices, per-spin orthonormalization. Validated
+# against the native dmftproj run without -so. Reuses CaOs2.struct/dmftsym.
+file(COPY CaOs2_noso.indmftpr CaOs2_noso.almblmup.gz CaOs2_noso.almblmdn.gz
+          CaOs2_noso.ctqmcout.gz CaOs2_noso.symqmc.gz
+          CaOs2_noso_partial.indmftpr CaOs2_noso_partial.sympar.gz CaOs2_noso_partial.parproj.gz
+          CaOs2_noso_partial_paramag.parproj.gz
+          CaOs2_noso_paramag.ctqmcout.gz CaOs2_noso_paramag.symqmc.gz
+          CaOs2_noso_partial_paramag.sympar.gz
+          CaOs2_noso_s.indmftpr CaOs2_noso_s.symqmc.gz CaOs2_noso_s.sympar.gz
+          CaOs2_noso_s_paramag.ctqmcout.gz CaOs2_noso_s_paramag.symqmc.gz
+          CaOs2_noso_s_paramag.sympar.gz
+     DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+add_test(NAME Py_wien2k_noso_python
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_noso_python.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_noso_python APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Odd-l shells (p l=1, f l=3): the improper-rotation parity (-1)^l is invisible
+# for the even-l (d) fixtures, so these pin dmat's use of iprop. p uses the
+# standard almblm (wide window, full rank); f uses a high-band-cutoff almblm so
+# the high-energy window carries real l=3 character and the Loewdin overlap is
+# full rank. Both reproduce the native dmftproj to machine precision.
+file(COPY CaOs2_pshell.indmftpr CaOs2_pshell.ctqmcout.gz CaOs2_pshell.symqmc.gz
+          CaOs2_fshell.indmftpr CaOs2_fshell.almblmup.gz CaOs2_fshell.almblmdn.gz
+          CaOs2_fshell.ctqmcout.gz CaOs2_fshell.symqmc.gz
+     DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
+add_test(NAME Py_wien2k_oddl_python
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_oddl_python.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_oddl_python APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
diff --git a/test/python/wien2k/CaOs2_fshell.almblmdn.gz b/test/python/wien2k/CaOs2_fshell.almblmdn.gz
new file mode 100644
index 0000000..ff9c44c
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.almblmdn.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.almblmup.gz b/test/python/wien2k/CaOs2_fshell.almblmup.gz
new file mode 100644
index 0000000..0fbf89f
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.almblmup.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.ctqmcout.gz b/test/python/wien2k/CaOs2_fshell.ctqmcout.gz
new file mode 100644
index 0000000..2561edd
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_fshell.indmftpr b/test/python/wien2k/CaOs2_fshell.indmftpr
new file mode 100644
index 0000000..9624e85
--- /dev/null
+++ b/test/python/wien2k/CaOs2_fshell.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+complex
+0 0 0 2
+0 0 0 0
+1
+3.0 10.0
diff --git a/test/python/wien2k/CaOs2_fshell.symqmc.gz b/test/python/wien2k/CaOs2_fshell.symqmc.gz
new file mode 100644
index 0000000..541d819
Binary files /dev/null and b/test/python/wien2k/CaOs2_fshell.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.almblmdn.gz b/test/python/wien2k/CaOs2_noso.almblmdn.gz
new file mode 100644
index 0000000..10bba32
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.almblmdn.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.almblmup.gz b/test/python/wien2k/CaOs2_noso.almblmup.gz
new file mode 100644
index 0000000..8080ae5
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.almblmup.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.ctqmcout.gz b/test/python/wien2k/CaOs2_noso.ctqmcout.gz
new file mode 100644
index 0000000..7c8b27d
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_noso.indmftpr b/test/python/wien2k/CaOs2_noso.indmftpr
new file mode 100644
index 0000000..b54bf2f
--- /dev/null
+++ b/test/python/wien2k/CaOs2_noso.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+0
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_noso.symqmc.gz b/test/python/wien2k/CaOs2_noso.symqmc.gz
new file mode 100644
index 0000000..d18da41
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_paramag.ctqmcout.gz b/test/python/wien2k/CaOs2_noso_paramag.ctqmcout.gz
new file mode 100644
index 0000000..64957e9
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_paramag.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_paramag.symqmc.gz b/test/python/wien2k/CaOs2_noso_paramag.symqmc.gz
new file mode 100644
index 0000000..cdd0257
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_paramag.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_partial.indmftpr b/test/python/wien2k/CaOs2_noso_partial.indmftpr
new file mode 100644
index 0000000..19b63bb
--- /dev/null
+++ b/test/python/wien2k/CaOs2_noso_partial.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 1 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+0
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_noso_partial.parproj.gz b/test/python/wien2k/CaOs2_noso_partial.parproj.gz
new file mode 100644
index 0000000..10156b7
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_partial.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_partial.sympar.gz b/test/python/wien2k/CaOs2_noso_partial.sympar.gz
new file mode 100644
index 0000000..0547daa
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_partial.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_partial_paramag.parproj.gz b/test/python/wien2k/CaOs2_noso_partial_paramag.parproj.gz
new file mode 100644
index 0000000..bb206b3
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_partial_paramag.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_partial_paramag.sympar.gz b/test/python/wien2k/CaOs2_noso_partial_paramag.sympar.gz
new file mode 100644
index 0000000..30e68ba
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_partial_paramag.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_s.indmftpr b/test/python/wien2k/CaOs2_noso_s.indmftpr
new file mode 100644
index 0000000..b536c82
--- /dev/null
+++ b/test/python/wien2k/CaOs2_noso_s.indmftpr
@@ -0,0 +1,12 @@
+2
+1 2
+3
+complex
+2 0 0 0
+0 0 0 0
+0
+cubic
+0 0 2 0
+0 0 0 0
+0
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_noso_s.sympar.gz b/test/python/wien2k/CaOs2_noso_s.sympar.gz
new file mode 100644
index 0000000..379df05
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_s.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_s.symqmc.gz b/test/python/wien2k/CaOs2_noso_s.symqmc.gz
new file mode 100644
index 0000000..379df05
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_s.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_s_paramag.ctqmcout.gz b/test/python/wien2k/CaOs2_noso_s_paramag.ctqmcout.gz
new file mode 100644
index 0000000..b18c119
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_s_paramag.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_s_paramag.sympar.gz b/test/python/wien2k/CaOs2_noso_s_paramag.sympar.gz
new file mode 100644
index 0000000..eebfe48
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_s_paramag.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_noso_s_paramag.symqmc.gz b/test/python/wien2k/CaOs2_noso_s_paramag.symqmc.gz
new file mode 100644
index 0000000..eebfe48
Binary files /dev/null and b/test/python/wien2k/CaOs2_noso_s_paramag.symqmc.gz differ
diff --git a/test/python/wien2k/CaOs2_pshell.ctqmcout.gz b/test/python/wien2k/CaOs2_pshell.ctqmcout.gz
new file mode 100644
index 0000000..5669218
Binary files /dev/null and b/test/python/wien2k/CaOs2_pshell.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_pshell.indmftpr b/test/python/wien2k/CaOs2_pshell.indmftpr
new file mode 100644
index 0000000..db81ca3
--- /dev/null
+++ b/test/python/wien2k/CaOs2_pshell.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+complex
+0 2 0 0
+0 0 0 0
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_pshell.symqmc.gz b/test/python/wien2k/CaOs2_pshell.symqmc.gz
new file mode 100644
index 0000000..29e2b7e
Binary files /dev/null and b/test/python/wien2k/CaOs2_pshell.symqmc.gz differ
diff --git a/test/python/wien2k/wien2k_noso_python.py b/test/python/wien2k/wien2k_noso_python.py
new file mode 100644
index 0000000..b60ee4a
--- /dev/null
+++ b/test/python/wien2k/wien2k_noso_python.py
@@ -0,0 +1,153 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.* generators on the non-spin-orbit path (ifSO=0),
+# the common DMFT case (e.g. SrVO3). dmftproj writes (2l+1) matrices there, not
+# the 2*(2l+1) spinor matrices of the SO path: the orthonormalization is
+# per-spin (orthogonal_wannier), the symmetry and Rloc matrices are the bare
+# (2l+1) representation, there is one density block per input spin channel,
+# and srot%timeinv is always false (setsym.f sets it only under SP+SO).
+#
+# CaOs2_noso is the spin-polarized, non-SO CaOs2 (indmftpr SO flag 0, wide
+# window -2.0 3.0): two equivalent Os d shells, cubic harmonics. The wide window
+# keeps the Loewdin overlap full rank. ctqmcout and symqmc are checked against
+# the native dmftproj reference for this case.
+#
+# CaOs2_noso_partial adds an uncorrelated Ca d shell, so sympar and parproj have
+# an included-but-not-correlated shell to symmetrize. Same struct/dmftsym/almblm
+# as the correlated-only case. References built with the native
+# dmftproj run without -so.
+#
+# The _paramag references use only the up-channel almblm as case.almblm (native
+# dmftproj without -sp). CaOs2_noso_s correlates the Ca s shell as well.
+
+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__))
+
+
+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_lines):
+    got = open(got_path).read().split()
+    ref = ''.join(ref_lines).split()
+    assert len(got) == len(ref), (len(got), len(ref))
+    floats = []
+    for a, b in zip(got, ref):
+        if b.lstrip('-').isdigit():
+            assert int(a) == int(b), (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)
+
+
+def _stage(case, paramagnetic=False):
+    tmp = tempfile.mkdtemp()
+    shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+                os.path.join(tmp, f'{case}.indmftpr'))
+    for ext in ('struct', 'dmftsym'):
+        shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+                    os.path.join(tmp, f'{case}.{ext}'))
+    if paramagnetic:
+        _gunzip(os.path.join(HERE, 'CaOs2_noso.almblmup.gz'),
+                os.path.join(tmp, f'{case}.almblm'))
+    else:
+        for ext in ('almblmup', 'almblmdn'):
+            _gunzip(os.path.join(HERE, f'CaOs2_noso.{ext}.gz'),
+                    os.path.join(tmp, f'{case}.{ext}'))
+    return tmp
+
+
+def _check_output(case, ext, writer, paramagnetic=False):
+    tmp = _stage(case, paramagnetic)
+    try:
+        writer(os.path.join(tmp, case))
+        reference = case + '_paramag' if paramagnetic else case
+        with gzip.open(os.path.join(HERE, f'{reference}.{ext}.gz'), 'rt') as fh:
+            _compare(os.path.join(tmp, f'{case}.{ext}'), fh.readlines())
+    finally:
+        shutil.rmtree(tmp)
+
+
+# Correlated-only: the projector and correlated-shell symmetry files.
+_check_output('CaOs2_noso', 'ctqmcout', ctqmcout.write_ctqmcout)
+_check_output('CaOs2_noso', 'symqmc', symqmc.write_symqmc)
+
+# Partial: the included-shell symmetry and partial-projector files.
+_check_output('CaOs2_noso_partial', 'sympar', sympar.write_sympar)
+_check_output('CaOs2_noso_partial', 'parproj', parproj.write_parproj)
+
+# A single non-magnetic input carries one density block and no spin flag; the
+# symmetry files end with the paramagnetic time-reversal operator.
+_check_output('CaOs2_noso_partial', 'parproj', parproj.write_parproj,
+              paramagnetic=True)
+_check_output('CaOs2_noso_partial', 'sympar', sympar.write_sympar,
+              paramagnetic=True)
+_check_output('CaOs2_noso', 'ctqmcout', ctqmcout.write_ctqmcout,
+              paramagnetic=True)
+_check_output('CaOs2_noso', 'symqmc', symqmc.write_symqmc, paramagnetic=True)
+
+# Correlated Ca s shell next to Os d: the scalar s-shell matrices, with and
+# without spin polarization.
+for ext, writer in (('symqmc', symqmc.write_symqmc), ('sympar', sympar.write_sympar)):
+    _check_output('CaOs2_noso_s', ext, writer)
+for ext, writer in (('ctqmcout', ctqmcout.write_ctqmcout),
+                    ('symqmc', symqmc.write_symqmc), ('sympar', sympar.write_sympar)):
+    _check_output('CaOs2_noso_s', ext, writer, paramagnetic=True)
+
+# Shells with different l on one atom remain independent in the non-SO path.
+orbs = [dict(atom=1, sort=1, l=l) for l in (0, 1)]
+raw = {(1, l): [np.eye(2 * l + 1) * value for value in (1., 2.)]
+       for l in (0, 1)}
+representations = {(l, 1, 0): (np.eye(2 * l + 1), False, 0.) for l in (0, 1)}
+result = parproj._symmetrize_densmat_nonso(
+    orbs, {}, [dict(perm=[1])], 1, representations, raw)
+for shell in raw:
+    for actual, expected in zip(result[shell], raw[shell]):
+        np.testing.assert_allclose(actual, expected, atol=1e-14, rtol=0)
+
+# SO metadata requires both spin channels, including when a single-channel
+# almblm happens to be present. Never export a half-sized SO projector.
+tmp = _stage('CaOs2_noso')
+try:
+    case = os.path.join(tmp, 'CaOs2_noso')
+    os.rename(case + '.almblmup', case + '.almblm')
+    os.unlink(case + '.almblmdn')
+    with open(case + '.indmftpr') as fh:
+        lines = fh.read().splitlines()
+    lines[-2] = '1'
+    with open(case + '.indmftpr', 'w') as fh:
+        fh.write('\n'.join(lines) + '\n')
+    for writer in (ctqmcout.write_ctqmcout, parproj.write_parproj):
+        try:
+            writer(case)
+        except ValueError as error:
+            assert 'both up and down' in str(error)
+        else:
+            raise AssertionError('single-channel SO projectors were accepted')
+finally:
+    shutil.rmtree(tmp)
+
+for module in (symqmc, sympar):
+    try:
+        module._shell_matrix({}, dict(l=1, mixing=True, P=np.eye(6)), 0, False)
+    except NotImplementedError as error:
+        assert 'spin-mixing' in str(error)
+    else:
+        raise AssertionError('non-SO symmetry accepted a spinor basis')
+
+print('wien2k_noso_python: ok')
diff --git a/test/python/wien2k/wien2k_oddl_python.py b/test/python/wien2k/wien2k_oddl_python.py
new file mode 100644
index 0000000..74e65e1
--- /dev/null
+++ b/test/python/wien2k/wien2k_oddl_python.py
@@ -0,0 +1,83 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the generators on the odd-l shells (p, l=1; f, l=3). The d-only fixtures
+# (l=2) cannot catch the improper-rotation parity: the Wigner-D parity factor is
+# (-1)^l, which is +1 for even l and only flips sign for odd l. case.dmftsym
+# stores the proper part of each operation in krotm (det(krotm)=+1 even for
+# improper ops) and the true parity in iprop; dmat must use iprop. The p- and
+# f-shell symqmc here pin that sign.
+#
+# Both are the spin-orbit CaOs2 cell (two equivalent Os atoms) with the Os shell
+# promoted to a complex p or f shell. The Loewdin overlap must be full rank:
+#   p (l=1): the standard almblm, wide window -2..3 (6 spin-orbitals/atom).
+#   f (l=3): Os carries no f weight near E_F, so the projector is rank-deficient
+#     on the standard bands. CaOs2_fshell.almblm is recomputed with a high band
+#     cutoff (lapw1/lapwso EMAX raised on the converged density); the window
+#     3..10 then sits on high-energy bands with real l=3 character, full rank
+#     (overlap condition number ~1.3).
+
+import gzip
+import os
+import shutil
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import ctqmcout, symqmc
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+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_gz):
+    got = open(got_path).read().split()
+    with gzip.open(ref_gz, 'rt') as fh:
+        ref = fh.read().split()
+    assert len(got) == len(ref), (len(got), len(ref))
+    floats = []
+    for a, b in zip(got, ref):
+        if b.lstrip('-').isdigit():
+            assert int(a) == int(b), (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)
+
+
+def _check(case, almblm_case, exts):
+    tmp = tempfile.mkdtemp()
+    try:
+        shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+                    os.path.join(tmp, f'{case}.indmftpr'))
+        for ext in ('struct', 'dmftsym'):
+            shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
+                        os.path.join(tmp, f'{case}.{ext}'))
+        for spin in ('up', 'dn'):
+            _gunzip(os.path.join(HERE, f'{almblm_case}.almblm{spin}.gz'),
+                    os.path.join(tmp, f'{case}.almblm{spin}'))
+        case_path = os.path.join(tmp, case)
+        for ext, writer in exts:
+            writer(case_path)
+            _compare(case_path + '.' + ext,
+                     os.path.join(HERE, f'{case}.{ext}.gz'))
+    finally:
+        shutil.rmtree(tmp)
+
+
+# p shell: standard SO almblm, full rank on the wide window.
+_check('CaOs2_pshell', 'CaOs2',
+       [('ctqmcout', ctqmcout.write_ctqmcout), ('symqmc', symqmc.write_symqmc)])
+
+# f shell: high-band-cutoff almblm, full rank on the high-energy window.
+_check('CaOs2_fshell', 'CaOs2_fshell',
+       [('ctqmcout', ctqmcout.write_ctqmcout), ('symqmc', symqmc.write_symqmc)])
+
+print('wien2k_oddl_python: ok')
diff --git a/test/python/wien2k/wien2k_sympar_python.py b/test/python/wien2k/wien2k_sympar_python.py
index d79866b..732a38c 100644
--- a/test/python/wien2k/wien2k_sympar_python.py
+++ b/test/python/wien2k/wien2k_sympar_python.py
@@ -70,7 +70,7 @@ def _check(case):
 inversion = dict(a=0., b=0., c=0., iprop=-1, krotm=np.eye(3))
 for l in (1, 3):
     shell = dict(l=l, mixing=False, P=np.eye(2 * l + 1, dtype=complex))
-    np.testing.assert_allclose(sympar._shell_matrix(inversion, shell, False),
+    np.testing.assert_allclose(sympar._shell_matrix(inversion, shell, False, True),
                                -np.eye(2 * (2 * l + 1)), atol=1e-14, rtol=0)
 p = sympar.reptrans('cubic', 1)
 np.testing.assert_allclose(p @ p.conj().T, np.eye(3), atol=1e-11, rtol=0)
diff --git a/test/python/wien2k/wien2k_symqmc_mixing.py b/test/python/wien2k/wien2k_symqmc_mixing.py
index 9702446..6b8765a 100644
--- a/test/python/wien2k/wien2k_symqmc_mixing.py
+++ b/test/python/wien2k/wien2k_symqmc_mixing.py
@@ -82,7 +82,7 @@ def _check(indmftpr, ref_npy, extra=(), unit_tol=1e-7):
 def _check_unitary(shells, ops, so, unit_tol):
     for op in ops:
         for sh in shells:
-            mat = symqmc._shell_matrix(op, sh, _timeinv(op, so))
+            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 < unit_tol, ('not unitary', dev)