From 41c5efd3a4084c54c472af996f762161ba5e65dd Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Tue, 16 Jun 2026 15:11:08 +0200 Subject: [PATCH 1/2] [wien2k] pure-Python case.ctqmcout correlated-shell projectors Read the almblm coefficients, assemble and locally rotate correlated projectors, apply the cubic/complex basis, and Loewdin-orthonormalize the stacked shells before writing ctqmcout. Part of #7, following #8. Keep the native separate window-open/valid flags. Validate matching SO spin windows and use each spin's own band offsets. Reject unsupported non-SO, band-index, cubic p/f and non-identity struct rotation paths instead of producing misleading files; subsequent port stages add them. Reuse the current full-rank CaOs2 inputs and #8 fixture handling. Compare integer fields exactly and finite float fields at atol=1e-7, rtol=0. Add regressions for final-band-only windows and mismatched SO spins. Local validation gives max native residual 1.3e-8 and full 20x50 projector orthonormality error 3.2e-15. The old rank-deficient-window claims no longer apply to this fixture. --- python/triqs_dftkit/wien2k/ctqmcout.py | 633 +++++++++++++++++++ test/python/wien2k/CMakeLists.txt | 7 + test/python/wien2k/wien2k_ctqmcout_python.py | 80 +++ 3 files changed, 720 insertions(+) create mode 100644 python/triqs_dftkit/wien2k/ctqmcout.py create mode 100644 test/python/wien2k/wien2k_ctqmcout_python.py diff --git a/python/triqs_dftkit/wien2k/ctqmcout.py b/python/triqs_dftkit/wien2k/ctqmcout.py new file mode 100644 index 0000000..839303f --- /dev/null +++ b/python/triqs_dftkit/wien2k/ctqmcout.py @@ -0,0 +1,633 @@ +################################################################################ +# +# 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) + 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)] + if any(n > 0 for n in ireps): + i += 1 + if correlated_ls: + so = int(raw[i].split()[0]) + 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]) + 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'])) + + # 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/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt index 05031b5..6b9c2ce 100644 --- a/test/python/wien2k/CMakeLists.txt +++ b/test/python/wien2k/CMakeLists.txt @@ -51,6 +51,13 @@ add_test(NAME Py_wien2k_oubwin_python 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}) + 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/wien2k_ctqmcout_python.py b/test/python/wien2k/wien2k_ctqmcout_python.py new file mode 100644 index 0000000..a85b5a1 --- /dev/null +++ b/test/python/wien2k/wien2k_ctqmcout_python.py @@ -0,0 +1,80 @@ +################################################################################ +# 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') + +print('wien2k_ctqmcout_python: ok') From 6d3188cf728cbbb1a715888a67958ceafeee25c3 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Mon, 12 Oct 2026 19:37:09 +0200 Subject: [PATCH 2/2] Validate correlated-shell records and projector overlap rank --- python/triqs_dftkit/wien2k/ctqmcout.py | 24 +++++++++++++++++--- test/python/wien2k/wien2k_ctqmcout_python.py | 19 ++++++++++++++++ 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/python/triqs_dftkit/wien2k/ctqmcout.py b/python/triqs_dftkit/wien2k/ctqmcout.py index 839303f..852389d 100644 --- a/python/triqs_dftkit/wien2k/ctqmcout.py +++ b/python/triqs_dftkit/wien2k/ctqmcout.py @@ -178,6 +178,8 @@ def _sqrt_inv(O): 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 @@ -204,10 +206,18 @@ def _read_indmftpr(indmftpr): 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)] - if any(n > 0 for n in ireps): - i += 1 + 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: - so = int(raw[i].split()[0]) + 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)) @@ -258,6 +268,8 @@ def _read_almblm(path, info): 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() @@ -463,6 +475,12 @@ def write_ctqmcout(case): 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): diff --git a/test/python/wien2k/wien2k_ctqmcout_python.py b/test/python/wien2k/wien2k_ctqmcout_python.py index a85b5a1..83f5fe1 100644 --- a/test/python/wien2k/wien2k_ctqmcout_python.py +++ b/test/python/wien2k/wien2k_ctqmcout_python.py @@ -77,4 +77,23 @@ def _check(case): 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')