diff --git a/python/triqs_dftkit/wien2k/ctqmcout.py b/python/triqs_dftkit/wien2k/ctqmcout.py
new file mode 100644
index 0000000..852389d
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/ctqmcout.py
@@ -0,0 +1,651 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj correlated-shell projector file
+(case.ctqmcout), replacing the projector path of the dmftproj Fortran
+executable (outputqmc.f subroutine outqmc).
+
+Given the alm/blm coefficient files (case.almblmup / case.almblmdn, one per
+spin), the projector definition (case.indmftpr), the structure
+(case.struct) and the symmetry input (case.dmftsym), this rebuilds the
+correlated (Wannier) projector pr_crorb%mat_rep and writes case.ctqmcout in
+the format the converter reads.
+
+It follows the Fortran path field by field:
+
+- almblm read order (dmftproj.f:559-689): header, per-sort per-l overlap
+ block, then per-k banners + idum/nbmin/nbmax + rtetr/eband lines, then the
+ per-(l,m) (Alm,Blm) / Clm complex coefficients in the canonical Wien2k lm
+ packing lm = l*l + (m+l) + 1.
+- band-window selection (set_projections.f:89-129, proj_mode==0).
+- raw correlated projector P(m,ib) = Alm(lm) + sum_ilo Clm(ilo,lm)*ovl_LO_u
+ (set_projections.f:243-267); Blm and the radial s12 do not enter pr_crorb.
+- local rotation rot_projectmat (Wigner-D of rotloc, identity here) then the
+ angular-basis transform mat_rep = transmat . (Rloc . P) where transmat =
+ is the cubic-harmonics transform (set_projections.f:450-478).
+- Loewdin orthonormalization of the stacked correlated projectors,
+ P <- O^{-1/2} P with O = D D^H over all correlated rows
+ (orthogonal_wannier_SO, dmftproj.f:847). This step is what the ctqmcout
+ numbers carry; it is applied per k over the full ndim x nbands stack.
+- the Rloc spinor block rotloc%rotrep (outputqmc.f:303-313), built from the
+ symmetry operation mapping the representative atom onto each equivalent one
+ (set_rotloc.f / setsym.f).
+
+ctqmcout layout is outputqmc.f:66-613. The cubic transform coefficients are
+stored in dmftproj via a single-precision CMPLX cast (set_ang_trans.f:146);
+to match the committed file the transmat is cast to complex64 before use
+(e.g. 0.70710676908 for 1/sqrt2). This step supports spin-polarized SO
+calculations with identity struct local rotations; other paths are rejected.
+"""
+
+import math
+import os
+import numpy as np
+
+
+# --- list-directed token stream over a Fortran free-form text file -----------
+
+class _Reader:
+ def __init__(self, path):
+ with open(path) as fh:
+ self._lines = fh.read().splitlines()
+ self._i = 0
+
+ def skip(self):
+ self._i += 1
+
+ def record(self):
+ toks = self._lines[self._i].split()
+ self._i += 1
+ return toks
+
+
+def _to_float(tok):
+ return float(tok.replace('D', 'E').replace('d', 'e'))
+
+
+def _to_complex(s):
+ s = s.strip().lstrip('(').rstrip(')')
+ re_s, im_s = s.split(',')
+ return complex(_to_float(re_s), _to_float(im_s))
+
+
+def _read_complex(r):
+ """Read one list-directed complex `(re,im)` value, possibly split over
+ whitespace tokens."""
+ buf = list(r.record())
+ while ')' not in ''.join(buf):
+ buf += r.record()
+ return _to_complex(''.join(buf))
+
+
+def _read_two_complex(r):
+ """Read a record holding two list-directed complex values (Alm, Blm)."""
+ buf = list(r.record())
+ while ''.join(buf).count(')') < 2:
+ buf += r.record()
+ s = ''.join(buf)
+ cut = s.index(')') + 1
+ return _to_complex(s[:cut]), _to_complex(s[cut:])
+
+
+# --- angular basis: cubic d transform transmat = ------------------
+# Wien2k cubic d convention dz2, dx2-y2, dxy, dxz, dyz (set_ang_trans cubic);
+# identical to symqmc._CUBIC[2].
+
+_CUBIC = {
+ 2: np.array([
+ [0, 0, 1, 0, 0],
+ [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
+ [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5],
+ [0, 2 ** -0.5, 0, -(2 ** -0.5), 0],
+ [0, 2 ** -0.5, 0, 2 ** -0.5, 0],
+ ], dtype=complex),
+}
+
+
+def _reptrans(basis, l):
+ """transmat = . dmftproj stores cubic/fromfile coefficients via a
+ single-precision CMPLX cast (set_ang_trans.f:146); reproduce the float32
+ truncation of these coefficients. Other native phase operations still
+ differ at the single-precision floor."""
+ if basis == 'cubic' and l in _CUBIC:
+ return _CUBIC[l].astype(np.complex64).astype(np.complex128)
+ if basis == 'cubic' and l != 0:
+ raise NotImplementedError('cubic ctqmcout currently supports s and d shells')
+ return np.eye(2 * l + 1, dtype=complex)
+
+
+# --- Wigner D matrix (dmftproj convention, setsym.f dmat) --------------------
+
+def _small_d(l, m, n, b):
+ f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
+ (math.factorial(l + n) * math.factorial(l - n))
+ s = 0.0
+ for t in range(0, 2 * l + 1):
+ if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
+ f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
+ (math.factorial(l - m - t) * math.factorial(m + n + t) *
+ math.factorial(l - n - t) * math.factorial(t))
+ f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 \
+ else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
+ f4 = 1.0 if (2 * t + n + m) == 0 \
+ else math.cos(b / 2) ** (2 * t + n + m)
+ s += (-1) ** (l - m - t) * f2 * f3 * f4
+ return math.sqrt(f1) * s
+
+
+def _dmat(l, a, b, c, det):
+ D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ for n in range(-l, l + 1):
+ v = np.exp(1j * n * a) * np.exp(1j * m * c) * _small_d(l, m, n, b)
+ if det < -0.5:
+ v *= (-1) ** l
+ D[m + l, n + l] = v
+ return D
+
+
+def _tmat(l):
+ T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ T[-m + l, m + l] = (-1) ** m
+ return T
+
+
+def _sqrt_inv(O):
+ """O^{-1/2} of a Hermitian matrix (orthogonal.f sqrtm with inv=.TRUE.:
+ Z diag(w^{-1/2}) Z^H). dmftproj evaluates 1/sqrt of the eigenvalue as a
+ *complex* sqrt (W_comp = CMPLX(W,0)), so a (numerically) negative
+ eigenvalue contributes i/sqrt(|w|) rather than NaN; reproduce that with a
+ complex power. The final result D1 @ conj(Z).T matches sqrtm exactly."""
+ w, Z = np.linalg.eigh(O)
+ if w[0] <= np.finfo(float).eps * len(w) * abs(w[-1]):
+ raise ValueError('Loewdin overlap is rank deficient; widen the projection window')
+ D1 = Z * (w.astype(complex) ** -0.5) # Z @ diag(w^{-1/2})
+ return D1 @ np.conj(Z).T
+
+
+# --- inputs ------------------------------------------------------------------
+
+def _read_indmftpr(indmftpr):
+ raw = [l.split('!')[0].strip() for l in open(indmftpr)]
+ raw = [l for l in raw if l != '']
+ nsort = int(raw[0].split()[0])
+ mult = [int(x) for x in raw[1].split()][:nsort]
+ lmax = int(raw[2].split()[0])
+ i = 3
+ so = 0
+ sorts = []
+ for isort in range(nsort):
+ basis = raw[i].split()[0]
+ i += 1
+ if basis == 'fromfile':
+ i += 1
+ l_inc = [int(x) for x in raw[i].split()]
+ i += 1
+ ireps = [int(x) for x in raw[i].split()]
+ i += 1
+ correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
+ included_ls = [l for l in range(len(l_inc)) if l_inc[l] in (1, 2)]
+ for l in correlated_ls:
+ if ireps[l] > 1 or (ireps[l] == 1 and raw[i].strip() != '1'):
+ raise NotImplementedError('selection of individual irreducible representations is not implemented')
+ if ireps[l] > 0:
+ i += 1
+ if correlated_ls:
+ sort_so = int(raw[i].split()[0])
+ if sort_so not in (0, 1):
+ raise ValueError('spin-orbit flags must be 0 or 1')
+ if any(s['correlated_ls'] for s in sorts) and sort_so != so:
+ raise NotImplementedError('different spin-orbit flags across correlated sorts are not implemented')
+ so = sort_so
+ i += 1
+ sorts.append(dict(basis=basis, correlated_ls=correlated_ls,
+ included_ls=included_ls))
+ last = raw[-1].split()
+ e_bot, e_top = _to_float(last[0]), _to_float(last[1])
+ proj_mode = int(last[2]) if len(last) >= 3 else 0
+ return dict(nsort=nsort, mult=mult, lmax=lmax, sorts=sorts, so=so,
+ e_bot=e_bot, e_top=e_top, proj_mode=proj_mode)
+
+
+def _read_dmftsym(path):
+ lines = open(path).read().split('\n')
+ nsym = int(lines[0].split()[0])
+ perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)]
+ starts = [i for i, l in enumerate(lines) if 'Sym. op.' in l]
+ ops = []
+ for k, s in enumerate(starts[:nsym]):
+ toks = lines[s + 1].split()
+ a, b, g = (math.radians(float(x)) for x in toks[:3])
+ iprop = int(toks[3])
+ krotm = np.array([[_to_float(x) for x in lines[s + 2 + r].split()]
+ for r in range(3)])
+ ops.append(dict(perm=perms[k], a=a, b=b, g=g, iprop=iprop,
+ krotm=krotm))
+ return nsym, ops
+
+
+# --- almblm parsing ----------------------------------------------------------
+
+def _read_almblm(path, info):
+ nsort = info['nsort']
+ lmax = info['lmax']
+ mult = info['mult']
+ nlm = (lmax + 1) ** 2
+ natom = sum(mult)
+
+ r = _Reader(path)
+ elecn = _to_float(r.record()[0])
+ nk = int(r.record()[0])
+ r.record() # nloat
+ eferm = _to_float(r.record()[0])
+
+ nLO = {}
+ ovl_LO_u = {}
+ kp = [None] * nk
+
+ for isrt in range(1, nsort + 1):
+ for l in range(lmax + 1):
+ r.record() # u_dot_norm(l)
+ n = int(r.record()[0])
+ if n > 1:
+ raise NotImplementedError('more than one local orbital per angular momentum is not implemented')
+ nLO[(l, isrt)] = n
+ for ilo in range(1, n + 1):
+ toks = r.record()
+ ovl_LO_u[(ilo, l, isrt)] = _to_float(toks[0])
+ for ik in range(nk):
+ r.skip() # "IK = .." banner
+ r.skip() # 3-int line
+ head = r.record()
+ nbmin, nbmax = int(head[1]), int(head[2])
+ nb = nbmax - nbmin + 1
+ if kp[ik] is None:
+ kp[ik] = dict(
+ nbmin=nbmin, nbmax=nbmax, eband=None, weight=None,
+ Alm=np.zeros((nlm, natom + 1, nb), dtype=complex),
+ Clm=np.zeros((1, nlm, natom + 1, nb), dtype=complex))
+ eband = np.empty(nb)
+ weight = None
+ for off in range(nb):
+ toks = r.record()
+ if off == 0:
+ weight = _to_float(toks[0])
+ eband[off] = _to_float(toks[1])
+ eband = eband - eferm
+ if kp[ik]['eband'] is None:
+ kp[ik]['eband'] = eband
+ kp[ik]['weight'] = weight
+ for imu in range(1, mult[isrt - 1] + 1):
+ iatom = sum(mult[:isrt - 1]) + imu
+ r.skip() # banner
+ r.record() # idum
+ for off in range(nb):
+ lm = 0
+ for l in range(lmax + 1):
+ for m in range(-l, l + 1):
+ alm, _blm = _read_two_complex(r)
+ kp[ik]['Alm'][lm, iatom, off] = alm
+ for ilo in range(nLO[(l, isrt)]):
+ clm = _read_complex(r)
+ kp[ik]['Clm'][ilo, lm, iatom, off] = clm
+ lm += 1
+
+ return dict(elecn=elecn, eferm=eferm, nk=nk, nlm=nlm, natom=natom,
+ nLO=nLO, ovl_LO_u=ovl_LO_u, kp=kp)
+
+
+# --- band-window selection (set_projections.f, proj_mode==0) -----------------
+
+def _select_window(nbmin, nbmax, eband, e1, e2):
+ opened = False
+ included = False
+ nb_bot = nb_top = 0
+ for off, ib in enumerate(range(nbmin, nbmax + 1)):
+ e = eband[off]
+ if not opened and e > e1 and e <= e2:
+ opened = True
+ nb_bot = ib
+ elif opened and e > e2:
+ nb_top = ib - 1
+ included = True
+ break
+ elif ib == nbmax and e > e1 and e <= e2:
+ nb_top = ib
+ included = True
+ if not included:
+ nb_bot = nb_top = 0
+ return included, nb_bot, nb_top
+
+
+# --- correlated / included orbital descriptors -------------------------------
+
+def _build_crorbs(info):
+ crorbs = []
+ for isort in range(1, info['nsort'] + 1):
+ sortinfo = info['sorts'][isort - 1]
+ for l in sortinfo['correlated_ls']:
+ for imu in range(1, info['mult'][isort - 1] + 1):
+ atom = sum(info['mult'][:isort - 1]) + imu
+ crorbs.append(dict(atom=atom, sort=isort, l=l,
+ basis=sortinfo['basis'], first=(imu == 1)))
+ return crorbs
+
+
+def _build_orbs(info):
+ orbs = []
+ for isort in range(1, info['nsort'] + 1):
+ sortinfo = info['sorts'][isort - 1]
+ for l in sortinfo['included_ls']:
+ for imu in range(1, info['mult'][isort - 1] + 1):
+ atom = sum(info['mult'][:isort - 1]) + imu
+ orbs.append(dict(atom=atom, sort=isort, l=l))
+ return orbs
+
+
+# --- 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)."""
+ a, b, g, iprop = op['a'], op['b'], op['g'], op['iprop']
+ 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)
+ 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)
+ return rotrep, timeinv
+
+
+# --- ctqmcout writer ---------------------------------------------------------
+
+def _fmt(x):
+ """One list-directed real, gfortran-style (leading sign space, ~17 sig)."""
+ return ' %.16E' % float(x)
+
+
+def _w(f, arr):
+ f.write(''.join(_fmt(x) for x in arr) + '\n')
+
+
+def _fmt_scalar(x):
+ x = float(x)
+ return '%.16f' % x if abs(x) < 1e5 else '%.16E' % x
+
+
+def write_ctqmcout(case):
+ """Read .almblm{up,dn}, .indmftpr, .struct,
+ .dmftsym and write .ctqmcout in the dmftproj format."""
+ info = _read_indmftpr(case + '.indmftpr')
+ if not info['so']:
+ raise NotImplementedError('non-spin-orbit ctqmcout is not implemented yet')
+ if info['proj_mode'] != 0:
+ raise NotImplementedError('ctqmcout currently supports projection mode 0')
+ with open(case + '.struct') as f:
+ struct_lines = f.readlines()
+ for i, line in enumerate(struct_lines):
+ if 'LOCAL ROT MATRIX:' in line:
+ rows = [line.split(':', 1)[1], *struct_lines[i + 1:i + 3]]
+ rotation = np.array([[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')
+ nsym, ops = _read_dmftsym(case + '.dmftsym')
+
+ up, dn = case + '.almblmup', case + '.almblmdn'
+ if os.path.exists(up) and os.path.exists(dn):
+ spin_files = [up, dn]
+ else:
+ spin_files = [case + '.almblm']
+ ifSP = len(spin_files) == 2
+ ifSO = bool(info['so'])
+ if not ifSP:
+ raise ValueError('spin-orbit ctqmcout requires both up and down almblm files')
+ ns = 2 if ifSP else 1
+
+ spins = [_read_almblm(p, info) for p in spin_files]
+ nk = spins[0]['nk']
+ elecn = spins[0]['elecn']
+
+ crorbs = _build_crorbs(info)
+ if any(cr['basis'] == 'fromfile' for cr in crorbs):
+ raise NotImplementedError(
+ 'ctqmcout for fromfile/mixing correlated bases is not yet covered '
+ 'by a test fixture; cubic and complex bases are supported')
+ orbs = _build_orbs(info)
+ ncrorb = len(crorbs)
+ norb = len(orbs)
+
+ # window in [e_bot, e_top]
+ windows = []
+ for ik in range(nk):
+ kp = spins[0]['kp'][ik]
+ windows.append(_select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top']))
+
+ for spin in spins[1:]:
+ if spin['nk'] != nk:
+ raise ValueError('up/down k-point counts differ')
+ for ik, window in enumerate(windows):
+ kp = spin['kp'][ik]
+ other = _select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top'])
+ if other != window:
+ raise ValueError('spin-orbit run requires identical up/dn band windows')
+ if not all(window[0] for window in windows):
+ raise ValueError('ctqmcout requires an included band window at every k-point')
+
+ # window below e_bot (for qbbot)
+ win_below = []
+ for ik in range(nk):
+ kp = spins[0]['kp'][ik]
+ win_below.append(_select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
+ -1e6, info['e_bot']))
+
+ for spin in spins[1:]:
+ for ik, kp in enumerate(spin['kp']):
+ lower = _select_window(kp['nbmin'], kp['nbmax'], kp['eband'], -1e6, info['e_bot'])
+ if lower != win_below[ik]:
+ raise ValueError('spin-orbit run requires identical up/dn below-window bands')
+
+ # qbbot: point integration over bands below e_bot, is=1 only under SO
+ qbbot = 0.0
+ for ispin in range(ns):
+ for ik in range(nk):
+ incl, nb_bot, nb_top = win_below[ik]
+ if incl:
+ qbbot += (nb_top - nb_bot + 1) * spins[ispin]['kp'][ik]['weight']
+ if ifSO:
+ break
+
+ transmats = {(cr['l'], cr['sort']): _reptrans(cr['basis'], cr['l'])
+ for cr in crorbs}
+
+ # rotloc Euler angles per crorb: the symmetry op mapping the representative
+ # atom of the sort onto this atom (set_rotloc.f). With identity struct
+ # local rotation, rotloc%(a,b,g,iprop) are that op's Euler angles, and
+ # 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 ----
+ # mat_rep[(icr, ik, is)] -> (2l+1, nbsel)
+ mat_rep = {}
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ atom = cr['atom']
+ sort = cr['sort']
+ transmat = transmats[(l, sort)]
+ op = rotloc_op[icr]
+ rot = _dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ for ispin in range(ns):
+ sp = spins[ispin]
+ kp = sp['kp'][ik]
+ off_bot = nb_bot - kp['nbmin']
+ off_top = nb_top - kp['nbmin']
+ nlo = sp['nLO'][(l, sort)]
+ nbsel = off_top - off_bot + 1
+ P = np.zeros((2 * l + 1, nbsel), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l) # 0-based packed index
+ for j, off in enumerate(range(off_bot, off_top + 1)):
+ val = kp['Alm'][lm, atom, off]
+ for ilo in range(nlo):
+ val += kp['Clm'][ilo, lm, atom, off] * \
+ sp['ovl_LO_u'][(ilo + 1, l, sort)]
+ P[mi, j] = val
+ # rot_projectmat (local rotation) then the basis transform
+ P = transmat @ (rot @ P)
+ mat_rep[(icr, ik, ispin)] = P
+
+ # Loewdin: per k stack all crorb rows (is=1 block then is=2 block per crorb)
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ blocks = []
+ layout = [] # (icr, is, nrows)
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ for ispin in range(ns):
+ blocks.append(mat_rep[(icr, ik, ispin)])
+ layout.append((icr, ispin, 2 * l + 1))
+ D = np.vstack(blocks) # ndim x nbnd
+ O = D @ np.conj(D.T)
+ S = _sqrt_inv(O)
+ D_orth = S @ D
+ row = 0
+ for icr, ispin, nrows in layout:
+ mat_rep[(icr, ik, ispin)] = D_orth[row:row + nrows, :]
+ row += nrows
+
+ # ---- Rloc rotrep per crorb ----
+ rloc_blocks = []
+ for cr in crorbs:
+ l = cr['l']
+ transmat = transmats[(l, cr['sort'])]
+ iref = sum(info['mult'][:cr['sort'] - 1]) + 1
+ op = next(o for o in ops if o['perm'][iref - 1] == cr['atom'])
+ rloc_blocks.append(_rloc_rotrep(op, l, transmat))
+
+ # ---- write ----
+ with open(case + '.ctqmcout', 'w') as f:
+ f.write('13.605698\n')
+ f.write('%6d\n' % nk)
+ f.write('%6d\n' % (1 if ifSP else 0))
+ f.write('%6d\n' % (1 if ifSO else 0))
+ f.write(' %s\n' % _fmt_scalar(qbbot))
+ f.write(' %s\n' % _fmt_scalar(elecn))
+
+ f.write('%6d\n' % norb)
+ for o in orbs:
+ dim = 2 * (2 * o['l'] + 1) if ifSO else 2 * o['l'] + 1
+ f.write('%6d %6d %6d %6d \n' % (o['atom'], o['sort'], o['l'], dim))
+
+ f.write('%6d\n' % ncrorb)
+ for cr in crorbs:
+ 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))
+
+ # Rloc block per crorb (non-mixing SP+SO whole shell)
+ for rotrep, timeinv in rloc_blocks:
+ for m in range(rotrep.shape[0]):
+ _w(f, rotrep[m, :].real)
+ for m in range(rotrep.shape[0]):
+ _w(f, rotrep[m, :].imag)
+ f.write('%6d\n' % (1 if timeinv else 0))
+
+ # complex-harmonics -> basis transform block (crorb%first only)
+ for cr in crorbs:
+ if not cr['first']:
+ continue
+ l = cr['l']
+ transmat = transmats[(l, cr['sort'])]
+ d = 2 * l + 1
+ 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):
+ _w(f, spinrot[m, :].real)
+ for m in range(2 * d):
+ _w(f, spinrot[m, :].imag)
+
+ # number of bands per k (skip is=2 under SO)
+ for ispin in range(ns):
+ if ifSP and ifSO and ispin == 1:
+ continue
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ f.write('%6d\n' % abs(nb_top - nb_bot + 1))
+
+ # projector block: DO ik, DO icrorb (non-mixing SO whole-shell)
+ for ik in range(nk):
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ _w(f, P[mi, :].real)
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ _w(f, P[mi, :].imag)
+
+ # k-weights
+ for ik in range(nk):
+ f.write(' %s\n' % _fmt_scalar(spins[0]['kp'][ik]['weight']))
+
+ # H(k) eigenvalues (skip is=2 under SO)
+ for ispin in range(ns):
+ if ifSP and ifSO and ispin == 1:
+ continue
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ kp = spins[ispin]['kp'][ik]
+ for ib in range(nb_bot, nb_top + 1):
+ f.write(' %s\n' % _fmt_scalar(kp['eband'][ib - kp['nbmin']]))
diff --git a/python/triqs_dftkit/wien2k/oubwin.py b/python/triqs_dftkit/wien2k/oubwin.py
new file mode 100644
index 0000000..a575536
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/oubwin.py
@@ -0,0 +1,279 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by L. Pourovskii, V. Vildosola, C. Martins, M. Aichhorn
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj band-window file (case.oubwin),
+replacing the corresponding output of the dmftproj Fortran executable.
+
+Given the alm/blm coefficient file (case.almblm, or the spin pair
+case.almblmup / case.almblmdn) and the projector definition (case.indmftpr),
+this selects, per k-point and spin, the contiguous band range whose Fermi-
+shifted Kohn-Sham eigenvalues fall in the energy window, and writes the
+result in the case.oubwin format read by the charge self-consistency.
+
+It reproduces the Fortran path: read elecn/nk/nloat/eferm from the almblm
+header (outbwin.f / dmftproj.f), parse the energy window (e_bot, e_top,
+proj_mode) from the last line of case.indmftpr, shift every band eigenvalue
+by eferm, and apply the proj_mode==0 selection of set_projections.f
+(strict lower bound e_bot < E, inclusive upper bound E <= e_top, yielding a
+single contiguous [nb_bot, nb_top] per k). The weight written per included
+k-point is the tetrahedron weight of the lowest band nbmin.
+
+Spin-polarization is detected from the input files: case.almblmup /
+case.almblmdn present => two spin files (case.oubwinup / case.oubwindn),
+otherwise the single case.oubwin. The spin-orbit flag (ifSO) comes from
+case.indmftpr; with -so the up and dn windows must coincide and dmftproj
+aborts otherwise, a check this generator also performs.
+"""
+
+import os
+import numpy as np
+
+
+# --- list-directed token stream over a Fortran free-form text file -----------
+
+class _Reader:
+ """Record-oriented reader mirroring Fortran list-directed and READ('()')
+ semantics on a free-form almblm file.
+
+ `record()` returns the whitespace tokens of the next physical line (the
+ common case: every almblm logical record fits on one line). `skip()`
+ consumes one physical line unconditionally, matching READ('()')."""
+
+ def __init__(self, path):
+ with open(path) as fh:
+ self._lines = fh.read().splitlines()
+ self._i = 0
+
+ def skip(self):
+ self._i += 1
+
+ def record(self):
+ toks = self._lines[self._i].split()
+ self._i += 1
+ return toks
+
+
+def _to_float(tok):
+ """Parse a Fortran real token, accepting D/d exponents (1.0D-3) and the
+ list-directed 4.57E-002 form. Python float already handles E/e and the
+ embedded sign in the exponent; only D/d needs translation."""
+ return float(tok.replace('D', 'E').replace('d', 'e'))
+
+
+# --- inputs ------------------------------------------------------------------
+
+def _read_window(indmftpr):
+ """Energy/band window from the last non-empty line of case.indmftpr.
+
+ Try list-directed `e_bot e_top proj_mode`; on failure fall back to
+ `e_bot e_top` with proj_mode=0 (dmftproj.f:182-197)."""
+ lines = [l.split('!')[0].strip() for l in open(indmftpr)]
+ lines = [l for l in lines if l]
+ toks = lines[-1].split()
+ e_bot, e_top = _to_float(toks[0]), _to_float(toks[1])
+ proj_mode = int(toks[2]) if len(toks) >= 3 else 0
+ return e_bot, e_top, proj_mode
+
+
+def _read_indmftpr_geometry(indmftpr):
+ """Return (nsort, lmax) from case.indmftpr. nsort is its first value and
+ lmax its third (dmftproj.f:95). Both are needed to walk the per-sort
+ almblm header up to the k-loop."""
+ raw = [l.split('!')[0].strip() for l in open(indmftpr)]
+ raw = [l for l in raw if l != '']
+ nsort = int(raw[0].split()[0])
+ lmax = int(raw[2].split()[0])
+ return nsort, lmax
+
+
+def _read_so_flag(indmftpr):
+ """Spin-orbit flag from case.indmftpr: the single integer line that
+ follows a correlated sort's irep specification (the value dmftproj stores
+ as ifSOflag). Returns 0 or 1."""
+ raw = [l.split('!')[0].strip() for l in open(indmftpr)]
+ raw = [l for l in raw if l != '']
+ nsort = int(raw[0].split()[0])
+ i = 3
+ so = 0
+ seen_correlated = False
+ for _ in range(nsort):
+ basis = raw[i].split()[0]
+ i += 1
+ if basis == 'fromfile':
+ i += 1
+ l_inc = [int(x) for x in raw[i].split()]
+ i += 1
+ ireps = [int(x) for x in raw[i].split()]
+ i += 1
+ correlated = any(v == 2 for v in l_inc)
+ for l, included in enumerate(l_inc):
+ if included == 2 and ireps[l] > 0:
+ i += 1
+ if correlated:
+ sort_so = int(raw[i].split()[0])
+ if sort_so not in (0, 1):
+ raise ValueError('spin-orbit flags must be 0 or 1')
+ if seen_correlated and sort_so != so:
+ raise NotImplementedError('different spin-orbit flags across correlated sorts are not implemented')
+ so = sort_so
+ seen_correlated = True
+ i += 1
+ return so
+
+
+# --- almblm parsing ----------------------------------------------------------
+
+def _read_almblm_windows(path, nsort, lmax):
+ """Read case.almblm and return (eferm, [(nbmin, nbmax, eband, weight) per
+ k-point]) for the FIRST atomic sort. nbmin/nbmax/eband/weight are identical
+ across sorts at a given k, so the first sort's k-loop suffices for window
+ selection.
+
+ Follows the dmftproj.f read order field by field: header, then the first
+ sort's per-l header (u_dot_norm, nLO, nLO overlap lines), then its k-loop
+ (two skipped banner lines, the idum/nbmin/nbmax line, the per-band
+ rtetr/eband lines)."""
+ r = _Reader(path)
+ r.record() # elecn
+ nk = int(r.record()[0])
+ r.record() # nloat
+ eferm = _to_float(r.record()[0]) # non-band path: eferm is header line 4
+
+ # First sort's per-l header: u_dot_norm, nLO, then nLO overlap lines.
+ for _ in range(lmax + 1):
+ r.record() # u_dot_norm(l)
+ nlo = int(r.record()[0])
+ for _ in range(nlo):
+ r.record() # ovl_LO_u, ovl_LO_udot
+
+ kpoints = []
+ for _ in range(nk):
+ r.skip() # "IK = .. jatom=.." banner
+ r.skip() # 3-int line
+ idum_nb = r.record()
+ nbmin, nbmax = int(idum_nb[1]), int(idum_nb[2])
+ eband = np.empty(nbmax - nbmin + 1)
+ weight = None
+ for off in range(nbmax - nbmin + 1):
+ toks = r.record()
+ rtetr = _to_float(toks[0])
+ eband[off] = _to_float(toks[1])
+ if off == 0: # tetrahedron weight of the lowest band
+ weight = rtetr
+ kpoints.append((nbmin, nbmax, eband, weight))
+ return eferm, kpoints
+
+
+# --- band-window selection (set_projections.f, proj_mode==0) -----------------
+
+def _select_window(nbmin, nbmax, eband, eferm, e_bot, e_top):
+ """proj_mode==0 contiguous band selection for one k-point.
+
+ Scan absolute band index ib = nbmin..nbmax over the Fermi-shifted energy
+ E = eband - eferm. The first ib with e_bot < E <= e_top opens the window;
+ the first ib above it with E > e_top closes it at ib-1 (scan exits). A
+ window already open on reaching nbmax closes at nbmax. As in Fortran,
+ opening at nbmax alone does not validate the window. Returns
+ (included, nb_bot, nb_top), with zero bounds for an invalid window."""
+ opened = False
+ included = False
+ nb_bot = nb_top = 0
+ for off, ib in enumerate(range(nbmin, nbmax + 1)):
+ e = eband[off] - eferm
+ if not opened and e_bot < e <= e_top:
+ opened = True
+ nb_bot = ib
+ elif opened and e > e_top:
+ nb_top = ib - 1
+ included = True
+ break
+ elif ib == nbmax and e_bot < e <= e_top:
+ nb_top = ib
+ included = True
+ if not included:
+ nb_bot = nb_top = 0
+ return included, nb_bot, nb_top
+
+
+def _windows_for_spin(almblm, indmftpr, nsort, lmax, e_bot, e_top, proj_mode):
+ """Per-k (included, nb_bot, nb_top, weight) for one spin's almblm file."""
+ if proj_mode != 0:
+ raise NotImplementedError(
+ 'oubwin generation is implemented for the energy-window '
+ 'projection mode (proj_mode==0); band-index modes 1 and 2 are '
+ 'not yet covered by a test fixture')
+ eferm, kpoints = _read_almblm_windows(almblm, nsort, lmax)
+ out = []
+ for nbmin, nbmax, eband, weight in kpoints:
+ incl, nb_bot, nb_top = _select_window(
+ nbmin, nbmax, eband, eferm, e_bot, e_top)
+ out.append((incl, nb_bot, nb_top, weight))
+ return out
+
+
+# --- output ------------------------------------------------------------------
+
+def _write_oubwin_file(path, ifso, windows):
+ """Write one case.oubwin spin file (outbwin.f). Integer fields are i6;
+ the per-k weight is a list-directed real. Not-included k-points emit only
+ the flag line."""
+ with open(path, 'w') as f:
+ f.write('%6d\n' % len(windows))
+ f.write('%6d\n' % (1 if ifso else 0))
+ for incl, nb_bot, nb_top, weight in windows:
+ f.write('%6d\n' % (1 if incl else 0))
+ if incl:
+ f.write('%6d%6d\n' % (nb_bot, nb_top))
+ f.write(' %.16f \n' % weight)
+
+
+def write_oubwin(case):
+ """Read .almblm{up,dn} (or .almblm) and write the matching
+ .oubwin{up,dn} (or .oubwin) in the dmftproj format.
+
+ Spin-polarization is inferred from the input files. The spin-orbit flag
+ written into every file comes from case.indmftpr; under -sp+-so the up/dn
+ windows must coincide (dmftproj aborts otherwise)."""
+ indmftpr = case + '.indmftpr'
+ nsort, lmax = _read_indmftpr_geometry(indmftpr)
+ ifso = _read_so_flag(indmftpr)
+ e_bot, e_top, proj_mode = _read_window(indmftpr)
+
+ up, dn = case + '.almblmup', case + '.almblmdn'
+ spin_polarized = os.path.exists(up) and os.path.exists(dn)
+
+ if spin_polarized:
+ win_up = _windows_for_spin(
+ up, indmftpr, nsort, lmax, e_bot, e_top, proj_mode)
+ win_dn = _windows_for_spin(
+ dn, indmftpr, nsort, lmax, e_bot, e_top, proj_mode)
+ if ifso:
+ for u, d in zip(win_up, win_dn):
+ if u[0] != d[0] or u[1] != d[1] or u[2] != d[2]:
+ raise ValueError(
+ 'spin-orbit run requires identical up/dn band '
+ 'windows at every k-point')
+ _write_oubwin_file(case + '.oubwinup', ifso, win_up)
+ _write_oubwin_file(case + '.oubwindn', ifso, win_dn)
+ else:
+ win = _windows_for_spin(
+ case + '.almblm', indmftpr, nsort, lmax, e_bot, e_top, proj_mode)
+ _write_oubwin_file(case + '.oubwin', ifso, win)
diff --git a/python/triqs_dftkit/wien2k/sympar.py b/python/triqs_dftkit/wien2k/sympar.py
new file mode 100644
index 0000000..2e2c980
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/sympar.py
@@ -0,0 +1,298 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
+#
+# TRIQS is free software: you can redistribute it and/or modify it under the
+# terms of the GNU General Public License as published by the Free Software
+# Foundation, either version 3 of the License, or (at your option) any later
+# version.
+#
+# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# TRIQS. If not, see .
+#
+################################################################################
+
+"""Pure-Python generation of the dmftproj included-shell symmetry file
+(case.sympar), replacing the corresponding output of the dmftproj Fortran
+executable (outputqmc.f, unit ousympar).
+
+case.sympar is the sibling of case.symqmc. Both write the same per-operation
+spinor symmetry matrix srot(isym)%rotrep(l,isrt)%mat (setsym.f), built
+identically for every shell as transmat . D(R)_lm . transmat^dag. They differ
+only in the shell list: symqmc covers the CORRELATED shells (l_inc==2), sympar
+covers ALL INCLUDED shells (l_inc in {1,2}), iterating orb(1:norb) instead of
+crorb(1:ncrorb). The orb ordering is (sort, l, atom) per dmftproj.f:357-387.
+
+The per-shell basis transmat is per sort: a complex shell uses the exact
+identity, a cubic shell uses the single-precision CMPLX cast of the cubic-d
+harmonics (set_ang_trans.f:146), giving 0.70710676908 for 1/sqrt2 — the same
+truncation ctqmcout carries. Mixed within one file (Ca complex exact, Os cubic
+truncated for CaOs2).
+
+sympar drops three things symqmc/ctqmcout carry: no orbital-description block
+at the top, no ifsplit/irep sub-block selection (always the full matrix), and
+no crorb%ifSOat/correp metadata. It is purely group-theoretic: header, perms,
+optional timeflag line (ifSP), the representation matrices, and — only when
+.not.ifSP — a paramagnetic time-reversal operator per orbital.
+
+This reference case (CaOs2) is SP+SO: indmftpr SO flag=1 sets ifSO, and
+ifSO=>ifSP, so it exercises the non-mixing 2*(2l+1)=10-wide block-diag spinor
+path and writes the timeflag line; the paramagnetic tail is absent.
+"""
+
+import math
+import os
+import numpy as np
+
+
+# --- angular bases -----------------------------------------------------------
+# transpose(P) = , m = -l..l. dmftproj stores cubic/fromfile coefficients
+# via a single-precision CMPLX cast (set_ang_trans.f:146); reproduce the
+# float32 truncation so the written numbers match (0.70710676908 for 1/sqrt2).
+# A complex basis is the exact identity (no cast).
+
+_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),
+ 2: np.array([
+ [0, 0, 1, 0, 0],
+ [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
+ [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5],
+ [0, 2 ** -0.5, 0, -(2 ** -0.5), 0],
+ [0, 2 ** -0.5, 0, 2 ** -0.5, 0],
+ ], dtype=complex),
+}
+
+
+def _reptrans(basis, l):
+ if basis == 'cubic' and l in _CUBIC:
+ return _CUBIC[l].astype(np.complex64).astype(np.complex128)
+ if basis == 'cubic' and l != 0:
+ raise NotImplementedError('cubic f-shell sympar is not implemented yet')
+ return np.eye(2 * l + 1, dtype=complex)
+
+
+# --- Wigner D matrix (dmftproj convention, setsym.f) -------------------------
+
+def _small_d(l, m, n, b):
+ f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
+ (math.factorial(l + n) * math.factorial(l - n))
+ s = 0.0
+ for t in range(0, 2 * l + 1):
+ if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
+ f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
+ (math.factorial(l - m - t) * math.factorial(m + n + t) *
+ math.factorial(l - n - t) * math.factorial(t))
+ f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 \
+ else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
+ f4 = 1.0 if (2 * t + n + m) == 0 \
+ else math.cos(b / 2) ** (2 * t + n + m)
+ s += (-1) ** (l - m - t) * f2 * f3 * f4
+ return math.sqrt(f1) * s
+
+
+def _dmat(l, a, b, c, det):
+ D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ for n in range(-l, l + 1):
+ v = np.exp(1j * n * a) * np.exp(1j * m * c) * _small_d(l, m, n, b)
+ if det < -0.5:
+ v *= (-1) ** l
+ D[m + l, n + l] = v
+ return D
+
+
+def _tmat(l):
+ """T[m, -m] = (-1)^m, complex conjugation in the |lm> basis (timeinv.f)."""
+ T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ T[-m + l, m + l] = (-1) ** m
+ return T
+
+
+def _timeinv_orbital(l, mat):
+ return _tmat(l) @ np.conj(mat)
+
+
+# --- inputs ------------------------------------------------------------------
+
+def _read_indmftpr(indmftpr):
+ """Return one entry per INCLUDED shell (l_inc in {1,2}), one per atom of its
+ sort, in orb order (sort, l, atom). Each shell carries l, basis name and the
+ transform P. Also returns the SO flag."""
+ raw = [l.split('!')[0].strip() for l in open(indmftpr)]
+ raw = [l for l in raw if l != '']
+ nsort = int(raw[0].split()[0])
+ mult = [int(x) for x in raw[1].split()][:nsort]
+ i = 3
+ so = 0
+ seen_correlated = False
+ shells = []
+ for isort in range(nsort):
+ basis = raw[i].split()[0]
+ i += 1
+ if basis == 'fromfile':
+ i += 1
+ l_inc = [int(x) for x in raw[i].split()]
+ i += 1
+ ireps = [int(x) for x in raw[i].split()]
+ i += 1
+ correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
+ included_ls = [l for l in range(len(l_inc)) if l_inc[l] in (1, 2)]
+ for l in correlated_ls:
+ if ireps[l] > 1 or (ireps[l] == 1 and raw[i].strip() != '1'):
+ raise NotImplementedError('selection of individual irreducible representations is not implemented')
+ if ireps[l] > 0:
+ i += 1
+ if correlated_ls:
+ sort_so = int(raw[i].split()[0])
+ if sort_so not in (0, 1):
+ raise ValueError('spin-orbit flags must be 0 or 1')
+ if seen_correlated and sort_so != so:
+ raise NotImplementedError('different spin-orbit flags across correlated sorts are not implemented')
+ so = sort_so
+ seen_correlated = True
+ i += 1
+ for l in included_ls:
+ P = _reptrans(basis, l)
+ for imu in range(1, mult[isort] + 1):
+ atom = sum(mult[:isort]) + imu
+ shells.append(dict(l=l, sort=isort + 1, atom=atom,
+ basis=basis, P=P))
+ return shells, so
+
+
+def _read_dmftsym(path):
+ lines = open(path).read().split('\n')
+ nsym = int(lines[0].split()[0])
+ perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)]
+ starts = [i for i, l in enumerate(lines) if 'Sym. op.' in l]
+ ops = []
+ for k, s in enumerate(starts[:nsym]):
+ toks = lines[s + 1].split()
+ a, b, c = (math.radians(float(x)) for x in toks[:3])
+ krotm = np.array([[float(x.replace('D', 'E').replace('d', 'e'))
+ for x in lines[s + 2 + r].split()] for r in range(3)])
+ ops.append(dict(perm=perms[k], a=a, b=b, c=c, iprop=int(toks[3]),
+ krotm=krotm))
+ return nsym, ops
+
+
+# --- matrix construction (shared with symqmc) --------------------------------
+
+def _phase(op, ti):
+ a, c = op['a'], op['c']
+ return (c - a) if ti else (a + c) # (g-a) on magnetic ops, else (a+g)
+
+
+def _l0_matrix(op, ti):
+ e = np.exp(1j * _phase(op, ti) / 2)
+ return np.array([[e, 0], [0, np.conj(e)]], dtype=complex)
+
+
+def _nonmixing_matrix(op, shell, ti):
+ """Non-mixing SP+SO whole shell: the up/up block scaled by +-(a+g)/2,
+ block-diagonal over spin, with the orbital time-reversal operator on the
+ magnetic operations (setsym.f, outputqmc.f:1292-1339)."""
+ l, P = shell['l'], shell['P']
+ rotl = _dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+ if ti:
+ rotl = _timeinv_orbital(l, rotl)
+ rotrep = P @ rotl @ np.conj(P.T)
+ e = np.exp(1j * _phase(op, ti) / 2)
+ d = 2 * l + 1
+ mat = np.zeros((2 * d, 2 * d), dtype=complex)
+ mat[:d, :d] = e * rotrep
+ mat[d:, d:] = np.conj(e) * rotrep
+ return mat
+
+
+def _shell_matrix(op, shell, ti):
+ l = shell['l']
+ if l == 0:
+ return _l0_matrix(op, ti)
+ return _nonmixing_matrix(op, shell, ti)
+
+
+# --- output ------------------------------------------------------------------
+
+def _fmt(x):
+ """One gfortran list-directed real (leading sign space, ~17 sig figs)."""
+ return ' %.16E' % float(x)
+
+
+def _w(f, arr):
+ f.write(''.join(_fmt(x) for x in arr) + '\n')
+
+
+def _write_matrix(f, mat):
+ for m in range(mat.shape[0]):
+ _w(f, mat[m, :].real)
+ for m in range(mat.shape[0]):
+ _w(f, mat[m, :].imag)
+
+
+def write_sympar(case):
+ """Write .sympar from .dmftsym, .indmftpr, .struct.
+
+ Detects ifSP/ifSO from the presence of .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."""
+ shells, so = _read_indmftpr(case + '.indmftpr')
+ if not so:
+ raise NotImplementedError('non-spin-orbit sympar is not implemented yet')
+ if any(sh['basis'] == 'fromfile' for sh in shells):
+ raise NotImplementedError(
+ 'sympar for fromfile/mixing bases is not yet covered by a test '
+ 'fixture; cubic and complex bases are supported')
+ 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 (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')
+ if ifSP:
+ f.write(''.join('%6d ' % t for t in timeflag) + '\n')
+
+ for isym, op in enumerate(ops):
+ for sh in shells:
+ l = sh['l']
+ ti = bool(timeflag[isym])
+ if l == 0 and not (ifSP and ifSO):
+ f.write(_fmt(1.0) + '\n')
+ f.write(_fmt(0.0) + '\n')
+ continue
+ _write_matrix(f, _shell_matrix(op, sh, ti))
+
+ if not ifSP:
+ for sh in shells:
+ l = sh['l']
+ if l == 0:
+ f.write(_fmt(1.0) + '\n')
+ f.write(_fmt(0.0) + '\n')
+ continue
+ tm = _tmat(l)
+ op = sh['P'] @ tm @ sh['P'].T
+ ident = np.eye(2 * l + 1, dtype=complex)
+ time_op = op @ np.conj(ident)
+ _write_matrix(f, time_op)
diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt
index c43cbed..fac3a52 100644
--- a/test/python/wien2k/CMakeLists.txt
+++ b/test/python/wien2k/CMakeLists.txt
@@ -42,6 +42,29 @@ 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})
+# 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.
+add_test(NAME Py_wien2k_oubwin_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_oubwin_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_oubwin_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Pure-Python correlated-shell projector generator vs the native fixture.
+add_test(NAME Py_wien2k_ctqmcout_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_ctqmcout_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_ctqmcout_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+# Included-shell symmetry generator, including an uncorrelated Ca d shell.
+add_test(NAME Py_wien2k_sympar_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_sympar_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_sympar_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
add_test(NAME Py_wien2k_kpath_labels
COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_kpath_labels.py
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
diff --git a/test/python/wien2k/CaOs2_partial.indmftpr b/test/python/wien2k/CaOs2_partial.indmftpr
new file mode 100644
index 0000000..afcca02
--- /dev/null
+++ b/test/python/wien2k/CaOs2_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
+1
+-0.15 0.30
diff --git a/test/python/wien2k/CaOs2_partial.sympar.gz b/test/python/wien2k/CaOs2_partial.sympar.gz
new file mode 100644
index 0000000..d553f10
Binary files /dev/null and b/test/python/wien2k/CaOs2_partial.sympar.gz differ
diff --git a/test/python/wien2k/wien2k_ctqmcout_python.py b/test/python/wien2k/wien2k_ctqmcout_python.py
new file mode 100644
index 0000000..83f5fe1
--- /dev/null
+++ b/test/python/wien2k/wien2k_ctqmcout_python.py
@@ -0,0 +1,99 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.ctqmcout generator (triqs_dftkit.wien2k.ctqmcout)
+# reproduces the dmftproj Fortran correlated-shell projectors: regenerate
+# case.ctqmcout from case.almblm{up,dn} + case.indmftpr + case.struct +
+# case.dmftsym and compare to the committed dmftproj reference.
+#
+# CaOs2 is spin-orbit + spin-polarized, cubic, two equivalent correlated d atoms.
+# Integer fields must match exactly; floats are compared at 1e-7: dmftproj stores
+# the cubic basis transform through a single-precision CMPLX cast, so its
+# projectors carry ~1e-7 noise that the double-precision generator reproduces to
+# the same floor (cf. the symqmc mixing test).
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import ctqmcout
+
+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_path):
+ got = open(got_path).read().split()
+ ref = open(ref_path).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)))
+ a, b = np.array(floats).T
+ np.testing.assert_allclose(a, b, atol=1e-7, rtol=0, equal_nan=False)
+
+
+def _check(case):
+ tmp = tempfile.mkdtemp()
+ try:
+ for ext in ('indmftpr', 'struct', 'dmftsym'):
+ shutil.copy(os.path.join(HERE, f'{case}.{ext}'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ for ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'{case}.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ ctqmcout.write_ctqmcout(os.path.join(tmp, case))
+ _compare(os.path.join(tmp, f'{case}.ctqmcout'),
+ os.path.join(HERE, f'{case}.ctqmcout'))
+ # A mismatched SO spin window must fail before replacing the output.
+ dn = os.path.join(tmp, f'{case}.almblmdn')
+ with open(dn) as f:
+ lines = f.readlines()
+ lines[3] = '100.0\n' # Shift all down-spin bands outside the window.
+ with open(dn, 'w') as f:
+ f.writelines(lines)
+ try:
+ ctqmcout.write_ctqmcout(os.path.join(tmp, case))
+ except ValueError as error:
+ assert 'identical up/dn band windows' in str(error), error
+ else:
+ raise AssertionError('mismatched SO spin windows were accepted')
+ finally:
+ shutil.rmtree(tmp)
+
+
+assert ctqmcout._select_window(10, 13, [-3., -2., -1., 0.5], 0., 1.) == (False, 0, 0)
+assert ctqmcout._select_window(5, 5, [0.5], 0., 1.) == (False, 0, 0)
+_check('CaOs2')
+
+for scale in (1., 1e-20):
+ overlap = np.diag([scale, 2 * scale]).astype(complex)
+ inverse = ctqmcout._sqrt_inv(overlap)
+ np.testing.assert_allclose(inverse @ overlap @ inverse.conj().T,
+ np.eye(2), atol=1e-14, rtol=0)
+try:
+ ctqmcout._sqrt_inv(np.diag([0., 1.]).astype(complex))
+except ValueError as error:
+ assert 'rank deficient' in str(error)
+else:
+ raise AssertionError('singular overlap was accepted')
+
+with tempfile.TemporaryDirectory() as tmp:
+ path = os.path.join(tmp, 'multi.indmftpr')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\ncomplex\n0 2 2\n0 1 1\n1\n1\n1\n-2 3\n')
+ parsed = ctqmcout._read_indmftpr(path)
+ assert parsed['so'] == 1 and parsed['sorts'][0]['correlated_ls'] == [1, 2]
+
+print('wien2k_ctqmcout_python: ok')
diff --git a/test/python/wien2k/wien2k_oubwin_python.py b/test/python/wien2k/wien2k_oubwin_python.py
new file mode 100644
index 0000000..26b427a
--- /dev/null
+++ b/test/python/wien2k/wien2k_oubwin_python.py
@@ -0,0 +1,106 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.oubwin generator (triqs_dftkit.wien2k.oubwin)
+# reproduces the dmftproj Fortran output: regenerate case.oubwinup / case.oubwindn
+# from case.almblm{up,dn} + case.indmftpr and compare to the native reference,
+# ignoring only trailing whitespace removed from the shipped fixture.
+#
+# CaOs2 is spin-orbit + spin-polarized, so both outputs match CaOs2.oubwinup.
+# Reuse the compressed almblm inputs shipped with the SOC fixture.
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import oubwin
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+def _check(case):
+ tmp = tempfile.mkdtemp()
+ try:
+ for ext in ('almblmup', 'almblmdn'):
+ with gzip.open(os.path.join(HERE, f'{case}.{ext}.gz'), 'rb') as src, \
+ open(os.path.join(tmp, f'{case}.{ext}'), 'wb') as dst:
+ shutil.copyfileobj(src, dst)
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'), tmp)
+ oubwin.write_oubwin(os.path.join(tmp, case))
+ with open(os.path.join(HERE, f'{case}.oubwinup')) as f:
+ ref = [line.rstrip() for line in f]
+ for ext in ('oubwinup', 'oubwindn'):
+ with open(os.path.join(tmp, f'{case}.{ext}')) as f:
+ got = [line.rstrip() for line in f]
+ assert got == ref, f'{case}.{ext} differs from dmftproj reference'
+ finally:
+ shutil.rmtree(tmp)
+
+
+# Independent expected band ranges for set_projections.f's branch ordering.
+for nbmin, energies, expected in (
+ (10, [-3., -2., -1., 0.5], (False, 0, 0)), # opens only at nbmax
+ (5, [0.5], (False, 0, 0)), # nbmin == nbmax
+ (10, [0., 0.5, 1., 2.], (True, 11, 12)), # strict lower, inclusive upper
+ (10, [0.5, 1.], (True, 10, 11)), # open window reaches nbmax
+ (10, [-1., 0., 2.], (False, 0, 0))): # no in-window band
+ got = oubwin._select_window(nbmin, nbmin + len(energies) - 1,
+ np.array(energies), 0., 0., 1.)
+ assert got == expected, (energies, got, expected)
+
+# Three k-points with distinct weights, a Fermi shift of 10 Ry, and one
+# excluded point. The lowest-band weight is used even if that band is outside
+# the window; an excluded point emits no band range or weight record.
+with tempfile.TemporaryDirectory() as tmp:
+ almblm = os.path.join(tmp, 'case.almblm')
+ with open(almblm, 'w') as f:
+ f.write('''0
+3
+0
+10.0
+1.0
+0
+IK = 1
+0 0 0
+1 10 13
+0.25 9.0
+0.80 10.25
+0.60 10.75
+0.25 12.0
+IK = 2
+0 0 0
+2 20 21
+0.10 8.0
+0.10 9.0
+IK = 3
+0 0 0
+3 30 31
+0.75 10.25
+0.50 10.75
+''')
+ windows = oubwin._windows_for_spin(almblm, None, 1, 0, 0., 1., 0)
+ assert windows == [(True, 11, 12, 0.25), (False, 0, 0, 0.10),
+ (True, 30, 31, 0.75)], windows
+ output = os.path.join(tmp, 'case.oubwin')
+ oubwin._write_oubwin_file(output, 0, windows)
+ with open(output) as f:
+ assert [line.split() for line in f] == [
+ ['3'], ['0'], ['1'], ['11', '12'], ['0.2500000000000000'],
+ ['0'], ['1'], ['30', '31'], ['0.7500000000000000']]
+
+_check('CaOs2')
+
+with tempfile.TemporaryDirectory() as tmp:
+ path = os.path.join(tmp, 'multi.indmftpr')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\ncomplex\n0 2 2\n0 1 1\n1\n1\n0\n-2 3\n')
+ assert oubwin._read_so_flag(path) == 0
+ with open(path, 'a') as fh:
+ fh.write('-2 3 ! default energy-window mode\n')
+ assert oubwin._read_window(path) == (-2., 3., 0)
+
+print('wien2k_oubwin_python: ok')
diff --git a/test/python/wien2k/wien2k_sympar_python.py b/test/python/wien2k/wien2k_sympar_python.py
new file mode 100644
index 0000000..0316fd5
--- /dev/null
+++ b/test/python/wien2k/wien2k_sympar_python.py
@@ -0,0 +1,88 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.sympar generator (triqs_dftkit.wien2k.sympar)
+# reproduces the dmftproj Fortran output. sympar holds the symmetry matrices for
+# all INCLUDED orbitals (correlated + partial), so the CaOs2_partial fixture adds
+# an uncorrelated Ca d-shell to the correlated Os d-shells. Inputs are the shared
+# CaOs2 almblm/struct/dmftsym renamed to the CaOs2_partial case.
+#
+# Floats compared at 1e-7: dmftproj stores the cubic basis transform through a
+# single-precision CMPLX cast, so its matrices carry ~1e-7 noise the double-
+# precision generator reproduces to the same floor.
+
+import gzip
+import os
+import shutil
+import tempfile
+import numpy as np
+
+from triqs_dftkit.wien2k import sympar
+
+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))
+ # 16 operations, three atoms, time flags, then three 10x10 shell matrices.
+ header = 2 + 16 * 3 + 16
+ np.testing.assert_array_equal(np.array(got[:header], int),
+ np.array(ref[:header], int))
+ data = np.array(got[header:], float).reshape(16, 3, 2, 10, 10)
+ reference = np.array(ref[header:], float).reshape(data.shape)
+ matrices = data[:, :, 0] + 1j * data[:, :, 1]
+ expected = reference[:, :, 0] + 1j * reference[:, :, 1]
+ np.testing.assert_allclose(matrices, expected, atol=1e-7, rtol=0,
+ equal_nan=False)
+ for operation in matrices:
+ for matrix in operation:
+ np.testing.assert_allclose(matrix @ matrix.conj().T, np.eye(10),
+ atol=1e-7, rtol=0)
+
+
+def _check(case):
+ 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 ext in ('almblmup', 'almblmdn'):
+ _gunzip(os.path.join(HERE, f'CaOs2.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ sympar.write_sympar(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.sympar.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.sympar'), fh.readlines())
+ finally:
+ shutil.rmtree(tmp)
+
+
+# Inversion has proper part I but parity -1: odd-l matrices must be -I.
+inversion = dict(a=0., b=0., c=0., iprop=-1, krotm=np.eye(3))
+for l in (1, 3):
+ shell = dict(l=l, P=np.eye(2 * l + 1, dtype=complex))
+ np.testing.assert_allclose(sympar._shell_matrix(inversion, shell, False),
+ -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-7, rtol=0)
+
+_check('CaOs2_partial')
+
+with tempfile.TemporaryDirectory() as tmp:
+ path = os.path.join(tmp, 'multi.indmftpr')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\ncomplex\n0 2 2\n0 1 1\n1\n1\n1\n-2 3\n')
+ shells, so = sympar._read_indmftpr(path)
+ assert so == 1 and [shell['l'] for shell in shells] == [1, 2]
+
+print('wien2k_sympar_python: ok')