diff --git a/.repo-hygiene-allow b/.repo-hygiene-allow
new file mode 100644
index 0000000..86f9197
--- /dev/null
+++ b/.repo-hygiene-allow
@@ -0,0 +1 @@
+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/doc/ChangeLog.md b/doc/ChangeLog.md
index be9ffd0..e4a01d6 100644
--- a/doc/ChangeLog.md
+++ b/doc/ChangeLog.md
@@ -12,6 +12,7 @@ this section to the version number when the next release is cut.
* Run DMFT on the irreducible Brillouin zone: `convert_dft_input(use_ibz=...)` reduces the data to the irreducible k-points, sets `symm_op = 1` and writes the correlated-shell symmetry operations (real-harmonic rotations for `l = 0, 1, 2` and atom permutations, built from the space group in `vaspout.h5`) to `dft_symmcorr_input`, so that DFTTools and ModEST restore the full-BZ average by symmetrization. Enabled by default when VASP used symmetry and `vaspout.h5` holds the symmetry data, and only when it is known to work: correlated shells with `l <= 2`, a symmetry-closed orbital set with an orthonormal `TRANSFORM`, all equivalent atoms projected, no spin-orbit coupling, and the symmetrized IBZ occupation and local levels reproduce the full grid within `ibz_tol` (1e-3). Otherwise the converter falls back to the full grid with a note; `use_ibz=True` raises instead. The VASP `Driver` passes `use_ibz` through for CSC runs
### Wien2k
+* Add pure-Python band-mode projectors, projection modes 1/2, and custom `fromfile` bases, including spin-mixing bases, to the Wien2k projector generators
* Read the high-symmetry k-path labels from the end of `case.outband` and store them as `kpts_labels` / `kpts_labels_idx` in `dft_bands_input`, matching the VASP band conversion
### test
diff --git a/python/triqs_dftkit/wien2k/_dmftproj.py b/python/triqs_dftkit/wien2k/_dmftproj.py
new file mode 100644
index 0000000..c4db637
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/_dmftproj.py
@@ -0,0 +1,636 @@
+################################################################################
+#
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+#
+# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola, C. Martins
+#
+# 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 .
+#
+################################################################################
+
+"""Shared dmftproj machinery for the pure-Python case.* generators.
+
+The generators (symqmc, oubwin, ctqmcout, sympar, parproj, outband) reproduce
+different outputs of the dmftproj Fortran executable but draw on the same
+underlying pieces: the list-directed almblm reader, the Fortran real/complex
+token parsers, the Wigner-D rotation in the dmftproj convention, the angular
+basis transforms (transmat), the case.indmftpr / case.dmftsym parsers, the
+band-window selection and the gfortran-style real formatters.
+This module owns each of those exactly once.
+
+All generators use the exact double-precision cubic harmonics; with
+full-precision SRC_templates dmftproj reads the same coefficients, so the port
+reproduces its output to machine precision.
+"""
+
+import math
+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, *, lines=None):
+ if lines is None:
+ with open(path) as fh:
+ lines = fh.read().splitlines()
+ self._lines = lines
+ 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'))
+
+
+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 bases: transmat = , m = -l..l --------------------------
+# Standard cubic harmonics, Wien2k convention, exact double precision.
+
+_COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
+
+_CUBIC = {
+ 1: np.array([[2**-0.5, 0, -2**-0.5], [-1j*2**-0.5, 0, -1j*2**-0.5], [0, 1, 0]], dtype=complex),
+ 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 = : the exact double-precision cubic harmonics for a
+ cubic basis, the identity for a complex basis."""
+ if basis == 'cubic' and l in _CUBIC:
+ return _CUBIC[l]
+ if basis == 'cubic' and l != 0:
+ raise NotImplementedError('cubic f harmonics are not implemented yet')
+ return _COMPLEX[l]
+
+
+def read_fromfile(path, l):
+ """Parse a dmftproj fromfile basis (each line: the coefficients of a new
+ basis vector in {|m,up>, |m,dn>}, m = -l..l, real/imag interleaved; '*'
+ marks the end of an irep). Return (P, mixing) where P = is the
+ transform matrix (reptrans.transmat), full 2(2l+1) for a spin-mixing basis
+ or the (2l+1) up/up block otherwise."""
+ n = 2 * (2 * l + 1)
+ rows = []
+ for line in open(path):
+ line = line.rstrip('\n')
+ if not line.strip():
+ continue
+ body = line[1:]
+ vals = [to_float(x) for x in body.split()[:2 * n]]
+ if len(vals) != 2 * n:
+ raise ValueError(f'{path}: each basis row requires {2 * n} real coefficients')
+ rows.append([vals[2 * k] + 1j * vals[2 * k + 1] for k in range(n)])
+ if len(rows) == n:
+ break
+ R = np.array(rows) # R[i, :] = |new_i> in old basis
+ if R.shape != (n, n) or not np.all(np.isfinite(R)):
+ raise ValueError(f'{path}: expected a finite {n} by {n} basis matrix')
+ d = 2 * l + 1
+ up_up, up_dn = R[:d, :d], R[:d, d:]
+ dn_up, dn_dn = R[d:, :d], R[d:, d:]
+ mixing = not (np.array_equal(dn_dn, up_up) and
+ np.all(up_dn == 0) and np.all(dn_up == 0))
+ P = np.conj(R) if mixing else np.conj(up_up) # = conj()
+ return P, mixing
+
+
+# --- Wigner D matrix (dmftproj convention, setsym.f) -------------------------
+
+def small_d(l, m, n, b):
+ f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
+ (math.factorial(l + n) * math.factorial(l - n))
+ s = 0.0
+ for t in range(0, 2 * l + 1):
+ if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
+ f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
+ (math.factorial(l - m - t) * math.factorial(m + n + t) *
+ math.factorial(l - n - t) * math.factorial(t))
+ f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 \
+ else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
+ f4 = 1.0 if (2 * t + n + m) == 0 \
+ else math.cos(b / 2) ** (2 * t + n + m)
+ s += (-1) ** (l - m - t) * f2 * f3 * f4
+ return math.sqrt(f1) * s
+
+
+def dmat(l, a, b, c, det):
+ D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ for n in range(-l, l + 1):
+ v = np.exp(1j * n * a) * np.exp(1j * m * c) * small_d(l, m, n, b)
+ if det < -0.5:
+ v *= (-1) ** l
+ D[m + l, n + l] = v
+ return D
+
+
+def tmat(l):
+ """Complex-conjugation operator in the spherical-harmonic basis,
+ T[m, -m] = (-1)^m (timeinv.f)."""
+ T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
+ for m in range(-l, l + 1):
+ T[-m + l, m + l] = (-1) ** m
+ return T
+
+
+def timeinv_orbital(l, mat):
+ return tmat(l) @ np.conj(mat)
+
+
+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
+ bare operator; callers apply it as tinv @ conj(mat)."""
+ d = 2 * l + 1
+ tm = tmat(l)
+ tinv = np.zeros((2 * d, 2 * d), dtype=complex)
+ tinv[:d, d:] = -tm
+ tinv[d:, :d] = tm
+ return P @ tinv @ P.T
+
+
+def rotloc_rotl_so(l, ref, ops, iatom, iref):
+ """The composed 2(2l+1) Rloc rotation rotloc(iatom)%rotl(l) under SP+SO
+ (setsym.f:496-528 + set_rotloc.f): the representative spinor rotloc
+ spmt (x) D(rotloc_ref) composed with the first symmetry op R[isym] mapping
+ iref onto iatom, with the orbital time-reversal applied for the magnetic op.
+ Returns (rotl, timeinv); callers apply their own basis transform (the full
+ transmat for a mixing basis, blkdiag(transmat, transmat) otherwise)."""
+ d = 2 * l + 1
+ Dref = dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
+ f = (ref['a'] + ref['g']) / 2.0
+ spmt = np.zeros((2, 2), dtype=complex)
+ spmt[0, 0] = np.exp(1j * f) * math.cos(ref['b'] / 2.0)
+ spmt[1, 1] = np.conj(spmt[0, 0])
+ f = -(ref['a'] - ref['g']) / 2.0
+ spmt[0, 1] = np.exp(1j * f) * math.sin(ref['b'] / 2.0)
+ spmt[1, 0] = -np.conj(spmt[0, 1])
+ rotl = np.zeros((2 * d, 2 * d), dtype=complex)
+ rotl[:d, :d] = spmt[0, 0] * Dref
+ rotl[d:, d:] = spmt[1, 1] * Dref
+ rotl[:d, d:] = spmt[0, 1] * Dref
+ rotl[d:, :d] = spmt[1, 0] * Dref
+
+ op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+ det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+ - op['krotm'][0, 1] * op['krotm'][1, 0])
+ timeinv = det2 < 0.0
+ srot_phase = (op['g'] - op['a']) if timeinv else (op['a'] + op['g'])
+ rotl_sym = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ if timeinv:
+ rotl_sym = tmat(l) @ np.conj(rotl_sym)
+ ephase = np.exp(1j * srot_phase / 2.0)
+ tmp = np.zeros((2 * d, 2 * d), dtype=complex)
+ tmp[:d, :d] = ephase * rotl_sym
+ tmp[d:, d:] = np.conj(ephase) * rotl_sym
+ rotl = (tmp @ np.conj(rotl)) if timeinv else (tmp @ rotl)
+ return rotl, timeinv
+
+
+def mixing_rotrep(op, l, P, ti):
+ """The full 2(2l+1) spinor representation D(R)_{new_i} = P spinrot P^dag of
+ one symmetry operation in a spin-coupling (mixing) basis, with the spinor
+ time-reversal operator applied for the magnetic (timeinv) operations
+ (setsym.f spinrotmat + timeinv_op). This is srot%rotrep(l,isrt)%mat, shared
+ by the symqmc and sympar shell matrices."""
+ rotl = dmat(l, op['a'], op['b'], op['c'], float(op['iprop']))
+ phase = (op['c'] - op['a']) if ti else (op['a'] + op['c'])
+ e = np.exp(1j * phase / 2)
+ d = 2 * l + 1
+ spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
+ if ti: # beta = pi, block-antidiagonal
+ spinrot[:d, d:] = e * rotl
+ spinrot[d:, :d] = -np.conj(e) * rotl
+ else: # beta = 0, block-diagonal
+ spinrot[:d, :d] = e * rotl
+ spinrot[d:, d:] = np.conj(e) * rotl
+ rotrep = P @ spinrot @ np.conj(P.T)
+ if ti:
+ rotrep = mixing_timeinv_op(l, P) @ np.conj(rotrep)
+ return rotrep
+
+
+# --- shell symmetry matrices (shared by symqmc and sympar) -------------------
+
+def _phase(op, ti):
+ a, c = op['a'], op['c']
+ return (c - a) if ti else (a + c) # (g-a) on magnetic ops, else (a+g)
+
+
+def _l0_matrix(op, ti):
+ """s shell: 2x2 spin phase block diag(e, conj(e)) (outputqmc.f l==0)."""
+ e = np.exp(1j * _phase(op, ti) / 2)
+ return np.array([[e, 0], [0, np.conj(e)]], dtype=complex)
+
+
+def _nonmixing_matrix(op, shell, ti):
+ """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)
+ # Antiunitary basis changes conjugate the source coefficients.
+ rotrep = P @ rotl @ P.T
+ else:
+ 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):
+ """Spinor symmetry matrix for one shell under one operation."""
+ l = shell['l']
+ if l == 0:
+ return _l0_matrix(op, ti)
+ if shell['mixing']:
+ return mixing_rotrep(op, l, shell['P'], bool(ti))
+ return _nonmixing_matrix(op, shell, ti)
+
+
+def build_orbs(info):
+ """One entry per INCLUDED orbital (sort, l, atom order)."""
+ 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,
+ basis=sortinfo['basis']))
+ return orbs
+
+
+# --- case.indmftpr -----------------------------------------------------------
+
+def read_indmftpr(indmftpr, *, lines=None):
+ """Full structured parse of case.indmftpr.
+
+ Returns a dict with nsort, mult, lmax, the spin-orbit flag `so`, the energy
+ window (e_bot, e_top, proj_mode) and one `sorts` entry per atomic sort with
+ its basis name, optional fromfile sourcefile path, and the correlated /
+ included l lists (l_inc==2 / l_inc in {1,2}). Every generator derives its
+ shell or orbital list from this one parse. Optional lines replace the file
+ contents; relative basis paths still resolve beside indmftpr."""
+ if lines is None:
+ with open(indmftpr) as fh:
+ lines = fh.readlines()
+ raw = [l.split('!')[0].strip() for l in lines]
+ 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
+ sourcefile = None
+ if basis == 'fromfile':
+ sourcefile = os.path.join(os.path.dirname(indmftpr), raw[i])
+ i += 1
+ l_inc = [int(x) for x in raw[i].split()]
+ i += 1
+ ireps = [int(x) for x in raw[i].split()]
+ i += 1
+ correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
+ 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:
+ raise NotImplementedError('selection of individual irreducible representations is not implemented')
+ if ireps[l] > 0:
+ if raw[i].strip() != '1':
+ raise NotImplementedError('selection of individual irreducible representations is not implemented')
+ i += 1 # each correlated l has its own correps line
+ 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 # SO flag follows a correlated sort
+ i += 1
+ sorts.append(dict(basis=basis, sourcefile=sourcefile,
+ 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
+ if proj_mode not in (0, 1, 2) or not np.isfinite([e_bot, e_top]).all() or e_bot > e_top:
+ raise ValueError('projection window requires finite ordered bounds and mode 0, 1 or 2')
+ return dict(nsort=nsort, mult=mult, lmax=lmax, sorts=sorts, so=so,
+ e_bot=e_bot, e_top=e_top, proj_mode=proj_mode)
+
+
+# --- case.dmftsym ------------------------------------------------------------
+
+def read_dmftsym(path, rotloc=False):
+ """Symmetry operations from case.dmftsym. Each op carries (perm, a, b, g,
+ iprop, krotm); the third Euler angle is named both `c` and `g` for the
+ callers that use either name. With rotloc=True also parse the
+ 'Global->local' representative rotloc per sort and return it as a third
+ value (setsym.f:437-455)."""
+ 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, c=g, g=g, iprop=iprop,
+ krotm=krotm))
+ if not rotloc:
+ return nsym, ops
+
+ gl = next(i for i, l in enumerate(lines) if 'Global->local' in l)
+ rotloc_ref = []
+ j = gl + 1
+ while j < len(lines):
+ if lines[j].strip() == '' or not lines[j].split()[0].lstrip('-').isdigit():
+ j += 1
+ continue
+ # sort index line, then 3 krotm rows, then the Euler/iprop line.
+ if len(lines[j].split()) == 1:
+ krotm = np.array([[to_float(x) for x in lines[j + 1 + r].split()]
+ for r in range(3)])
+ ang = lines[j + 4].split()
+ a, b, g = (math.radians(float(x)) for x in ang[:3])
+ iprop = int(ang[3])
+ rotloc_ref.append(dict(krotm=krotm, a=a, b=b, g=g, iprop=iprop))
+ j += 5
+ else:
+ j += 1
+ return nsym, ops, rotloc_ref
+
+
+# --- almblm parsing ----------------------------------------------------------
+
+def read_almblm(path, info, projectors=False, *, lines=None):
+ """Read one spin's almblm file in the dmftproj.f read order.
+
+ With projectors=False keep only the header (elecn, nk, nloat, eferm), the
+ per-(sort,l) overlap block (u_dot_norm, nLO, ovl_LO_u, ovl_LO_udot) and the
+ per-k band data (nbmin/nbmax, Fermi-shifted eband, tetrahedron weights); the
+ per-(l,m) coefficient records are still consumed to keep the reader in sync
+ but their values are discarded (oubwin's window-only path).
+
+ With projectors=True also accumulate the per-(l,m) Alm/Blm/Clm coefficients
+ in the canonical Wien2k packing lm = l*l + (m+l), for every atom of every
+ sort (ctqmcout/parproj). Optional lines replace the file contents."""
+ nsort = info['nsort']
+ lmax = info['lmax']
+ mult = info['mult']
+ nlm = (lmax + 1) ** 2
+ natom = sum(mult)
+
+ r = Reader(path, lines=lines)
+ elecn = to_float(r.record()[0])
+ nk = int(r.record()[0])
+ r.record() # nloat
+ eferm = to_float(r.record()[0])
+
+ nLO = {}
+ u_dot_norm = {}
+ ovl_LO_u = {}
+ ovl_LO_udot = {}
+ kp = [None] * nk
+
+ for isrt in range(1, nsort + 1):
+ for l in range(lmax + 1):
+ u_dot_norm[(l, isrt)] = to_float(r.record()[0])
+ 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])
+ ovl_LO_udot[(ilo, l, isrt)] = to_float(toks[1])
+ 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, tetr=None)
+ if projectors:
+ kp[ik].update(
+ Alm=np.zeros((nlm, natom + 1, nb), dtype=complex),
+ Blm=np.zeros((nlm, natom + 1, nb), dtype=complex),
+ Clm=np.zeros((4, nlm, natom + 1, nb), dtype=complex))
+ eband = np.empty(nb)
+ tetr = np.empty(nb)
+ for off in range(nb):
+ toks = r.record()
+ tetr[off] = 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'] = tetr[0]
+ kp[ik]['tetr'] = tetr
+ 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)
+ for ilo in range(nLO[(l, isrt)]):
+ clm = read_complex(r)
+ if projectors:
+ kp[ik]['Clm'][ilo, lm, iatom, off] = clm
+ if projectors:
+ kp[ik]['Alm'][lm, iatom, off] = alm
+ kp[ik]['Blm'][lm, iatom, off] = blm
+ lm += 1
+
+ return dict(elecn=elecn, eferm=eferm, nk=nk, nlm=nlm, natom=natom,
+ nLO=nLO, u_dot_norm=u_dot_norm, ovl_LO_u=ovl_LO_u,
+ ovl_LO_udot=ovl_LO_udot, 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
+
+
+def select_band_window(nbmin, nbmax, b_bot, b_top):
+ """proj_mode 1/2 band-index selection for one k-point (set_projections.f
+ 70-88). b_bot/b_top are band indices, not energies: nb_bot = b_bot clamped
+ up to nbmin (strict INT(e1) > nbmin), nb_top = b_top clamped down to nbmax.
+ Returns (included, nb_bot, nb_top); included is False when the clamped
+ range is empty (nb_top < nb_bot), which contributes no bands."""
+ nb_bot = b_bot if b_bot > nbmin else nbmin
+ nb_top = b_top if b_top < nbmax else nbmax
+ return nb_bot <= nb_top, nb_bot, nb_top
+
+
+def band_index_window(info, spins):
+ """Resolve the (b_bot, b_top) band-index window for proj_mode 1 and 2.
+
+ proj_mode 2 (dmftproj.f:233-237): b_bot=INT(e_bot), b_top=INT(e_top) taken
+ directly from the indmftpr window line.
+
+ proj_mode 1 (dmftproj.f:704-722): e_bot/e_top are Fermi-shifted energies;
+ scan every spin and k-point for bands with e_bot < E <= e_top and take the
+ global min/max band index, seeded with b_bot=1000, b_top=1 so an empty scan
+ keeps that seed. `spins` is the list of read_almblm dicts (one per spin)."""
+ if info['proj_mode'] == 2:
+ return int(info['e_bot']), int(info['e_top'])
+ if info['proj_mode'] != 1:
+ raise ValueError('band_index_window is only valid for proj_mode 1 or 2')
+ e_bot, e_top = info['e_bot'], info['e_top']
+ b_bot, b_top = 1000, 1
+ for sp in spins:
+ for kp in sp['kp']:
+ for off, ib in enumerate(range(kp['nbmin'], kp['nbmax'] + 1)):
+ e = kp['eband'][off]
+ if e > e_bot and e <= e_top:
+ if ib > b_top:
+ b_top = ib
+ if ib < b_bot:
+ b_bot = ib
+ return b_bot, b_top
+
+
+# --- gfortran-style real formatters ------------------------------------------
+
+def fmt(x):
+ """One list-directed real, gfortran-style (leading sign space, ~17 sig)."""
+ return ' %.16E' % float(x)
+
+
+def write_row(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 require_identity_struct(case):
+ """Fail explicitly until representative struct rotations are ported."""
+ with open(case + '.struct') as f:
+ lines = f.readlines()
+ for i, line in enumerate(lines):
+ if 'LOCAL ROT MATRIX:' in line:
+ rows = [line.split(':', 1)[1], *lines[i + 1:i + 3]]
+ rotation = np.array([[to_float(x) for x in row.split()] for row in rows])
+ if not np.allclose(rotation, np.eye(3), atol=1e-7, rtol=0):
+ raise NotImplementedError('non-identity struct local rotations are not implemented yet')
+
+
+def validate_so_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')
+ bw = band_index_window(info, spins) if info['proj_mode'] else None
+ for spin in spins[1:]:
+ if spin['nk'] != spins[0]['nk']:
+ raise ValueError('up/down k-point counts differ')
+ for ik, kp in enumerate(spin['kp']):
+ 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')
+ 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')
+ 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 852389d..aac711b 100644
--- a/python/triqs_dftkit/wien2k/ctqmcout.py
+++ b/python/triqs_dftkit/wien2k/ctqmcout.py
@@ -49,293 +49,41 @@
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.
+ctqmcout layout is outputqmc.f:66-613. The cubic transform uses the exact
+analytic harmonics; against a dmftproj with double-precision casts and
+full-precision templates the port reproduces case.ctqmcout to machine
+precision on a full-rank window.
"""
-import math
import os
import numpy as np
+from scipy.linalg.lapack import zheev
+from scipy.linalg.blas import zgemm
-
-# --- 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
+from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_so_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)
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)
+ """O^{-1/2} of a Hermitian matrix, reproducing orthogonal.f sqrtm
+ (inv=.TRUE.: Z diag(w^{-1/2}) Z^H). Use the same LAPACK routine the Fortran
+ calls, ZHEEV('V','U'), via scipy so the eigenvectors match rather than the
+ divide-and-conquer ZHEEVD of numpy.linalg.eigh. dmftproj takes 1/sqrt of the
+ eigenvalue as a *complex* sqrt (W_comp = CMPLX(W,0)); reproduce that with a
+ complex power. The result D1 @ conj(Z).T matches the Fortran's ZGEMM('N','T').
+
+ Full-rank overlaps give a well-defined inverse square root. Singular
+ overlaps cannot define orthonormal correlated projectors and are rejected.
+ """
+ w, Z, info = zheev(O, compute_v=1, lower=0) # 'V', 'U'
+ if info != 0:
+ raise np.linalg.LinAlgError(f'Loewdin eigensolver failed (info={info})')
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
+ D1 = Z * (w.astype(complex) ** -0.5) # Z @ diag(w^{-1/2})
+ return zgemm(1.0, D1, np.conj(Z), trans_b=1) # ZGEMM('N','T'): D1 @ conj(Z)^T
# --- correlated / included orbital descriptors -------------------------------
@@ -352,17 +100,6 @@ def _build_crorbs(info):
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):
@@ -376,9 +113,9 @@ def _rloc_rotrep(op, l, transmat):
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))
+ D = dmat(l, a, b, g, float(iprop))
if timeinv:
- D = _tmat(l) @ np.conj(D) # setsym.f:320-326 orbital time reversal
+ 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)
@@ -391,39 +128,25 @@ def _rloc_rotrep(op, l, transmat):
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 _rloc_rotrep_mixing(l, ref, ops, iatom, iref, P):
+ """rotloc(iatom)%rotrep(l)%mat for a mixing SO shell: the composed Rloc
+ spinor rotation (rotloc_rotl_so) put into the new basis with the full
+ 2(2l+1) transmat P (set_rotloc.f mixing branch). Returns (rotrep, timeinv)."""
+ rotl, timeinv = rotloc_rotl_so(l, ref, ops, iatom, iref)
+ rotrep = (P @ rotl @ P.T) if timeinv else (P @ rotl @ np.conj(P.T))
+ return rotrep, timeinv
-def _fmt_scalar(x):
- x = float(x)
- return '%.16f' % x if abs(x) < 1e5 else '%.16E' % x
+# --- ctqmcout writer ---------------------------------------------------------
def write_ctqmcout(case):
"""Read .almblm{up,dn}, .indmftpr, .struct,
.dmftsym and write .ctqmcout in the dmftproj format."""
- info = _read_indmftpr(case + '.indmftpr')
+ 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')
+ raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
+ require_identity_struct(case)
+ nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
up, dn = case + '.almblmup', case + '.almblmdn'
if os.path.exists(up) and os.path.exists(dn):
@@ -432,54 +155,46 @@ def write_ctqmcout(case):
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]
+ spins = [read_almblm(p, info, projectors=True) 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]
+ # band-index window for proj_mode 1/2 (set_projections is called with band
+ # indices, not energies). proj_mode 1 scans all spins for the global window.
+ bw = band_index_window(info, spins) if info['proj_mode'] != 0 else None
+
+ # window in [e_bot, e_top] (proj_mode 0) or [b_bot, b_top] (mode 1/2)
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)
+ if info['proj_mode'] == 0:
+ windows.append(select_window(
+ kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top']))
+ else:
+ windows.append(select_band_window(
+ kp['nbmin'], kp['nbmax'], bw[0], bw[1]))
+
+ # window below e_bot (mode 0) or below b_bot (mode 1/2), for qbbot.
+ # Mode 1/2: set_projections(1, b_bot-1) -> select_band_window with top=b_bot-1.
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']))
+ if info['proj_mode'] == 0:
+ win_below.append(select_window(
+ kp['nbmin'], kp['nbmax'], kp['eband'], -1e6, info['e_bot']))
+ else:
+ win_below.append(select_band_window(
+ kp['nbmin'], kp['nbmax'], 1, bw[0] - 1))
- 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')
+ validate_so_windows(spins, info, windows, win_below)
# qbbot: point integration over bands below e_bot, is=1 only under SO
qbbot = 0.0
@@ -491,8 +206,15 @@ def write_ctqmcout(case):
if ifSO:
break
- transmats = {(cr['l'], cr['sort']): _reptrans(cr['basis'], cr['l'])
- for cr in crorbs}
+ transmats = {}
+ mixing = {}
+ for cr in crorbs:
+ key = (cr['l'], cr['sort'])
+ if cr['basis'] == 'fromfile':
+ transmats[key], mixing[key] = read_fromfile(
+ info['sorts'][cr['sort'] - 1]['sourcefile'], cr['l'])
+ else:
+ transmats[key], mixing[key] = reptrans(cr['basis'], cr['l']), False
# rotloc Euler angles per crorb: the symmetry op mapping the representative
# atom of the sort onto this atom (set_rotloc.f). With identity struct
@@ -507,27 +229,31 @@ def write_ctqmcout(case):
if o['perm'][iref - 1] == cr['atom'])
# ---- raw correlated projector mat_rep, then Loewdin orthonormalize ----
- # mat_rep[(icr, ik, is)] -> (2l+1, nbsel)
+ # Non-mixing: mat_rep[(icr, ik, is)] -> (2l+1, nbsel) per spin block.
+ # Mixing: mat_rep[(icr, ik)] -> (2*(2l+1), nbsel), the stacked spinor block.
mat_rep = {}
for icr, cr in enumerate(crorbs):
l = cr['l']
atom = cr['atom']
sort = cr['sort']
+ d = 2 * l + 1
transmat = transmats[(l, sort)]
+ ismix = mixing[(l, sort)]
op = rotloc_op[icr]
- rot = _dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ 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
+ nbsel = nb_top - nb_bot + 1
+ spin_blocks = []
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)
+ P = np.zeros((d, 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)):
@@ -536,29 +262,42 @@ def write_ctqmcout(case):
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
+ spin_blocks.append(rot @ P) # rot_projectmat per spin
+ if ismix:
+ stack = np.vstack(spin_blocks) # (2*(2l+1), nbsel), up then dn
+ mat_rep[(icr, ik)] = transmat @ stack
+ else:
+ for ispin in range(ns):
+ mat_rep[(icr, ik, ispin)] = transmat @ spin_blocks[ispin]
- # Loewdin: per k stack all crorb rows (is=1 block then is=2 block per crorb)
+ # 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):
incl, nb_bot, nb_top = windows[ik]
if not incl:
continue
blocks = []
- layout = [] # (icr, is, nrows)
+ layout = [] # (key, 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))
+ if mixing[(l, cr['sort'])]:
+ blocks.append(mat_rep[(icr, ik)])
+ layout.append(((icr, ik), 2 * (2 * l + 1)))
+ else:
+ for ispin in range(ns):
+ blocks.append(mat_rep[(icr, ik, ispin)])
+ layout.append(((icr, ik, ispin), 2 * l + 1))
D = np.vstack(blocks) # ndim x nbnd
- O = D @ np.conj(D.T)
+ # 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
S = _sqrt_inv(O)
- D_orth = S @ D
+ D_orth = zgemm(1.0, S, D) # ZGEMM('N','N'): O^{-1/2} @ D
row = 0
- for icr, ispin, nrows in layout:
- mat_rep[(icr, ik, ispin)] = D_orth[row:row + nrows, :]
+ for key, nrows in layout:
+ mat_rep[key] = D_orth[row:row + nrows, :]
row += nrows
# ---- Rloc rotrep per crorb ----
@@ -567,8 +306,12 @@ def write_ctqmcout(case):
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))
+ if mixing[(l, cr['sort'])]:
+ rloc_blocks.append(_rloc_rotrep_mixing(
+ l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat))
+ else:
+ op = next(o for o in ops if o['perm'][iref - 1] == cr['atom'])
+ rloc_blocks.append(_rloc_rotrep(op, l, transmat))
# ---- write ----
with open(case + '.ctqmcout', 'w') as f:
@@ -576,8 +319,8 @@ def write_ctqmcout(case):
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(' %s\n' % fmt_scalar(qbbot))
+ f.write(' %s\n' % fmt_scalar(elecn))
f.write('%6d\n' % norb)
for o in orbs:
@@ -594,26 +337,30 @@ def write_ctqmcout(case):
# 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)
+ write_row(f, rotrep[m, :].real)
for m in range(rotrep.shape[0]):
- _w(f, rotrep[m, :].imag)
+ write_row(f, rotrep[m, :].imag)
f.write('%6d\n' % (1 if timeinv else 0))
- # complex-harmonics -> basis transform block (crorb%first only)
+ # complex-harmonics -> basis transform block (crorb%first only). Mixing
+ # writes the full 2(2l+1) transmat; non-mixing the spin block-diagonal.
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
+ if mixing[(l, cr['sort'])]:
+ spinrot = transmat
+ else:
+ 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)
+ write_row(f, spinrot[m, :].real)
for m in range(2 * d):
- _w(f, spinrot[m, :].imag)
+ write_row(f, spinrot[m, :].imag)
# number of bands per k (skip is=2 under SO)
for ispin in range(ns):
@@ -623,22 +370,30 @@ def write_ctqmcout(case):
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)
+ # projector block: DO ik, DO icrorb. Mixing writes the full 2(2l+1)
+ # block from is=1; non-mixing the two (2l+1) spin blocks.
for ik in range(nk):
for icr, cr in enumerate(crorbs):
l = cr['l']
+ if mixing[(l, cr['sort'])]:
+ P = mat_rep[(icr, ik)]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
for ispin in range(ns):
P = mat_rep[(icr, ik, ispin)]
for mi in range(2 * l + 1):
- _w(f, P[mi, :].real)
+ write_row(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)
+ write_row(f, P[mi, :].imag)
# k-weights
for ik in range(nk):
- f.write(' %s\n' % _fmt_scalar(spins[0]['kp'][ik]['weight']))
+ f.write(' %s\n' % fmt_scalar(spins[0]['kp'][ik]['weight']))
# H(k) eigenvalues (skip is=2 under SO)
for ispin in range(ns):
@@ -648,4 +403,4 @@ def write_ctqmcout(case):
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']]))
+ 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
index a575536..a46e13c 100644
--- a/python/triqs_dftkit/wien2k/oubwin.py
+++ b/python/triqs_dftkit/wien2k/oubwin.py
@@ -31,10 +31,16 @@
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.
+by eferm, and apply the set_projections.f selection. For proj_mode==0 this is
+the energy-window rule (strict lower bound e_bot < E, inclusive upper bound
+E <= e_top, yielding a single contiguous [nb_bot, nb_top] per k). For
+proj_mode 1 and 2 the window is a pair of band indices (b_bot, b_top): every
+k-point is included with nb_bot/nb_top those indices clamped to the local
+nbmin/nbmax (set_projections.f:70-88). In mode 2 the indices come straight
+from the window line; in mode 1 they are the global min/max band index over
+all spins/k whose Fermi-shifted energy lies in (e_bot, e_top]
+(dmftproj.f:704-722). 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),
@@ -44,193 +50,28 @@
"""
import os
-import numpy as np
+from ._dmftproj import (band_index_window, read_almblm, read_indmftpr,
+ select_band_window, select_window)
-# --- 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)
+def _windows_from_spin(sp, info, band_window):
+ """Per-k (included, nb_bot, nb_top, weight) for one already-read spin dict.
+ proj_mode==0 uses the energy window; modes 1/2 use the band-index window
+ (b_bot, b_top) precomputed in band_window."""
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))
+ for kp in sp['kp']:
+ if info['proj_mode'] == 0:
+ incl, nb_bot, nb_top = select_window(
+ kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top'])
+ else:
+ incl, nb_bot, nb_top = select_band_window(
+ kp['nbmin'], kp['nbmax'], band_window[0], band_window[1])
+ out.append((incl, nb_bot, nb_top, kp['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
@@ -252,20 +93,25 @@ def write_oubwin(case):
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)
+ info = read_indmftpr(case + '.indmftpr')
+ ifso = bool(info['so'])
up, dn = case + '.almblmup', case + '.almblmdn'
spin_polarized = os.path.exists(up) and os.path.exists(dn)
+ if ifso and not spin_polarized:
+ raise ValueError('spin-orbit windows require both up and down almblm files')
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)
+ sp_up = read_almblm(up, info)
+ sp_dn = read_almblm(dn, info)
+ # mode 1 scans both spins for the global band-index window.
+ bw = (band_index_window(info, [sp_up, sp_dn])
+ if info['proj_mode'] != 0 else None)
+ win_up = _windows_from_spin(sp_up, info, bw)
+ win_dn = _windows_from_spin(sp_dn, info, bw)
if ifso:
+ if len(win_up) != len(win_dn):
+ raise ValueError('up/down k-point counts differ')
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(
@@ -274,6 +120,8 @@ def write_oubwin(case):
_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)
+ sp = read_almblm(case + '.almblm', info)
+ bw = (band_index_window(info, [sp])
+ if info['proj_mode'] != 0 else None)
+ win = _windows_from_spin(sp, info, bw)
_write_oubwin_file(case + '.oubwin', ifso, win)
diff --git a/python/triqs_dftkit/wien2k/outband.py b/python/triqs_dftkit/wien2k/outband.py
new file mode 100644
index 0000000..e16f302
--- /dev/null
+++ b/python/triqs_dftkit/wien2k/outband.py
@@ -0,0 +1,325 @@
+################################################################################
+#
+# 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-structure projector file
+(case.outband), replacing the -band path of the dmftproj Fortran executable
+(outband.f).
+
+case.outband feeds the converter method convert_bands_input: it is the
+band-structure analog of case.ctqmcout. The correlated-shell projectors are
+built exactly as in ctqmcout.py (raw Alm/Clm projector -> rot_projectmat local
+rotation -> cubic/fromfile transmat -> Loewdin orthonormalization
+orthogonal_wannier_SO), but evaluated on a band k-path and written in the
+outband layout. case.outband also carries the raw radial-normalized Theta
+projector pr_orb%matn_rep for every included shell (the parproj.py path, no
+Loewdin) and the k-path labels.
+
+Two band-mode differences from the ctqmcout path (dmftproj.f:557-581):
+
+- nkband, the number of k-points along the plotted k-path, is read from
+ .klist_band (read_k_list.f). It is written in the header and drives the
+ k-label block; the per-k projector/band data loop over the almblm nk, which is
+ the single representative k-point lapw2 -almd writes in band mode.
+- the Fermi energy is read from the LAST line of .indmftpr, not from the
+ almblm header; the band almblm keeps the eferm record as a placeholder that
+ the reader skips (READ(iualmblm,*) with no target).
+
+outband layout is outband.f:56-278: nkband, the per-spin per-k band counts, the
+Loewdin correlated projectors (DO ik, DO icrorb), H(k) eigenvalues, the
+included-orbital radial sizes, the raw Theta projectors, and the k-labels.
+"""
+
+import os
+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,
+ read_indmftpr, reptrans, select_window, to_float,
+ write_row)
+from .ctqmcout import _build_crorbs, _sqrt_inv
+from .parproj import (_build_matn_rep, _build_matn_rep_mixing, _build_orbs)
+
+
+def _read_indmftpr_band(indmftpr):
+ """case.indmftpr in band mode: the standard parse, plus the Fermi energy
+ appended as the last physical line after the (e_bot e_top proj_mode) window
+ line. The trailing scalar would be mistaken for the window line, so feed
+ read_indmftpr the lines without that record and return eferm separately."""
+ with open(indmftpr) as fh:
+ raw = [l.split('!')[0].rstrip('\n') for l in fh]
+ nonblank = [i for i, l in enumerate(raw) if l.strip() != '']
+ eferm = to_float(raw[nonblank[-1]].split()[0])
+ info = read_indmftpr(indmftpr, lines=raw[:nonblank[-1]])
+ return info, eferm
+
+
+def _read_almblm_band(path, info, eferm):
+ """Band-mode almblm read: identical to read_almblm but the header eferm
+ record is a placeholder the band path ignores; the Fermi shift uses the
+ eferm passed in (from case.indmftpr, dmftproj.f:576-581). read_almblm reads
+ the 4th record as eferm, so overwrite that record with the indmftpr value."""
+ with open(path) as fh:
+ text = fh.read().splitlines()
+ text[3] = ' %.16E' % eferm # Reader record 4 (physical line index 3)
+ return read_almblm(path, info, projectors=True, lines=text)
+
+
+def read_k_list(path):
+ """nkband and the k-path labels from .klist_band (read_k_list.f).
+
+ Every physical line up to the END marker is a k-point; a line whose first
+ character is not a blank carries a label, whose position is the k-point
+ index and whose name is the leading non-blank token. Returns
+ (nkband, [(pos, name), ...])."""
+ lines = open(path).read().splitlines()
+ nkband = 0
+ labels = []
+ for line in lines:
+ if line[:3] == 'END':
+ break
+ nkband += 1
+ if line[:1] != ' ':
+ labels.append((nkband, line.split()[0]))
+ return nkband, labels
+
+
+def write_outband(case):
+ """Read .almblm{up,dn}, .indmftpr, .struct,
+ .dmftsym and .klist_band and write .outband in the
+ dmftproj band format."""
+ info, eferm = _read_indmftpr_band(case + '.indmftpr')
+ if not info['so']:
+ raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
+ require_identity_struct(case)
+ if info['proj_mode'] != 0:
+ raise NotImplementedError('this output currently supports projection mode 0')
+ nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
+
+ 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'])
+ ns = 2 if ifSP else 1
+
+ spins = [_read_almblm_band(p, info, eferm) for p in spin_files]
+ nk = spins[0]['nk']
+ info['nk'] = nk
+
+ nkband, labels = read_k_list(case + '.klist_band')
+ if nkband != nk:
+ raise ValueError('klist_band and almblm k-point counts differ')
+
+ crorbs = _build_crorbs(info)
+ orbs = _build_orbs(info)
+ ncrorb = len(crorbs)
+ norb = len(orbs)
+
+ # proj_mode 0 energy window [e_bot, e_top] per k (band path: proj_mode 0).
+ 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']))
+
+ validate_so_windows(spins, info, windows)
+
+ transmats = {}
+ mixing = {}
+ for o in orbs:
+ key = (o['l'], o['sort'])
+ if key in transmats:
+ continue
+ if o['basis'] == 'fromfile':
+ transmats[key], mixing[key] = read_fromfile(
+ info['sorts'][o['sort'] - 1]['sourcefile'], o['l'])
+ else:
+ transmats[key], mixing[key] = reptrans(o['basis'], o['l']), False
+
+ # rot_projectmat local rotation: the op mapping the sort representative onto
+ # this atom (set_projections.f via rot_projectmat). Same as ctqmcout/parproj.
+ rotloc_op = {}
+ for o in orbs:
+ iref = sum(info['mult'][:o['sort'] - 1]) + 1
+ rotloc_op[o['atom']] = next(op for op in ops
+ if op['perm'][iref - 1] == o['atom'])
+
+ # ---- correlated projector mat_rep, then Loewdin orthonormalize ----
+ mat_rep = {}
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ atom = cr['atom']
+ sort = cr['sort']
+ d = 2 * l + 1
+ transmat = transmats[(l, sort)]
+ ismix = mixing[(l, sort)]
+ op = rotloc_op[atom]
+ rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ nbsel = nb_top - nb_bot + 1
+ spin_blocks = []
+ 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)]
+ P = np.zeros((d, nbsel), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l)
+ 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
+ spin_blocks.append(rot @ P)
+ if ismix:
+ stack = np.vstack(spin_blocks)
+ mat_rep[(icr, ik)] = transmat @ stack
+ else:
+ for ispin in range(ns):
+ mat_rep[(icr, ik, ispin)] = transmat @ spin_blocks[ispin]
+
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ blocks = []
+ layout = []
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ if mixing[(l, cr['sort'])]:
+ blocks.append(mat_rep[(icr, ik)])
+ layout.append(((icr, ik), 2 * (2 * l + 1)))
+ else:
+ for ispin in range(ns):
+ blocks.append(mat_rep[(icr, ik, ispin)])
+ layout.append(((icr, ik, ispin), 2 * l + 1))
+ D = np.vstack(blocks)
+ # 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
+ S = _sqrt_inv(O)
+ D_orth = zgemm(1.0, S, D) # ZGEMM('N','N'): O^{-1/2} @ D
+ row = 0
+ for key, nrows in layout:
+ mat_rep[key] = D_orth[row:row + nrows, :]
+ row += nrows
+
+ # ---- raw Theta projectors matn_rep per included orbital (parproj path) ----
+ matn_reps = {}
+ for o in orbs:
+ l = o['l']
+ transmat = transmats[(l, o['sort'])]
+ op = rotloc_op[o['atom']]
+ rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ if mixing[(l, o['sort'])]:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
+ o, info, spins, windows, transmat, rot)
+ else:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep(
+ o, info, spins, windows, transmat, rot)
+
+ # ---- write ----
+ with open(case + '.outband', 'w') as f:
+ f.write('%6d\n' % nkband)
+
+ # number of bands per k (skip is=2 under SP+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))
+
+ # correlated projector block: DO ik, DO icrorb (Loewdin orthonormalized).
+ for ik in range(nk):
+ for icr, cr in enumerate(crorbs):
+ l = cr['l']
+ if mixing[(l, cr['sort'])]:
+ P = mat_rep[(icr, ik)]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].real)
+ for ispin in range(ns):
+ P = mat_rep[(icr, ik, ispin)]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].imag)
+
+ # H(k) eigenvalues (skip is=2 under SO)
+ for ispin in range(ns):
+ if 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(' %.16E\n' % kp['eband'][ib - kp['nbmin']])
+
+ # included-orbital radial sizes norm_radf%n
+ for o in orbs:
+ n = spins[0]['nLO'][(o['l'], o['sort'])] + 2
+ f.write('%6d\n' % n)
+
+ # Theta projector block: DO iorb, DO ik, DO ir. Mixing writes the full
+ # 2(2l+1) block from is=1; non-mixing the two (2l+1) spin blocks.
+ for o in orbs:
+ l = o['l']
+ atom = o['atom']
+ n = spins[0]['nLO'][(l, o['sort'])] + 2
+ ismix = mixing[(l, o['sort'])]
+ for ik in range(nk):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ for ir in range(n):
+ if ismix:
+ P = matn_reps[(atom, l)][ik][:, :, ir]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ 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(ns):
+ P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
+ for mi in range(2 * l + 1):
+ write_row(f, P[mi, :].imag)
+
+ # k-labels
+ for i, (pos, name) in enumerate(labels, start=1):
+ f.write('%6d%6d%s\n' % (i, pos, name))
diff --git a/python/triqs_dftkit/wien2k/parproj.py b/python/triqs_dftkit/wien2k/parproj.py
index aefc927..fd74118 100644
--- a/python/triqs_dftkit/wien2k/parproj.py
+++ b/python/triqs_dftkit/wien2k/parproj.py
@@ -45,128 +45,27 @@
the per-orbital Rloc spinor rotrep (set_rotloc.f), and a time-reversal flag.
The reference case CaOs2 is SP+SO (ifSP=ifSO=1, ns=2), so every per-orbital
-matrix is the full 2*(2l+1)=10 spin+orbital block. The cubic transmat (cast to
-complex64 for the Os shells, as dmftproj does in set_ang_trans.f:146) and the
-rot_projectmat local rotation are applied to matn_rep before writing, exactly
-as in ctqmcout.
+matrix is the full 2*(2l+1)=10 spin+orbital block. The exact cubic transmat and
+the rot_projectmat local rotation are applied to matn_rep before writing,
+exactly as in ctqmcout.
outputqmc.f parproj writer: 897-1183. set_projections.f s12/projector:
163-197/273-622. density.f Theta path: 588-915. symmetrize_mat.f /
rot_dens.f / setsym.f / set_rotloc.f for the symmetrization and Rloc rotation.
"""
-import math
import os
import numpy as np
-from .ctqmcout import _select_window
+from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_so_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)
-# --- 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):
- 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 = ------------------
-
-_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 the cubic coefficients via a
- single-precision CMPLX cast (set_ang_trans.f:146); reproduce the float32
- truncation so the written numbers match bit-for-bit."""
- 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 parproj 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):
- """Complex-conjugation operator in the spherical-harmonic basis,
- T[m, -m] = (-1)^m (timeinv.f)."""
- T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
- for m in range(-l, l + 1):
- T[-m + l, m + l] = (-1) ** m
- return T
+def _op_timeinv(op):
+ det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+ - op['krotm'][0, 1] * op['krotm'][1, 0])
+ return bool(det2 < 0.0)
def _sqrtm_real_sym(O):
@@ -179,186 +78,6 @@ def _sqrtm_real_sym(O):
return (D1 @ Z.T).real
-# --- 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):
- """Symmetry operations and the per-sort representative Rloc rotation from
- case.dmftsym. The symop block carries (perm, a, b, g, iprop, krotm); the
- 'Global->local' tail carries the representative rotloc krotm + Euler angles
- per sort (setsym.f:437-455)."""
- 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))
- # Global->local representative rotloc per sort.
- gl = next(i for i, l in enumerate(lines) if 'Global->local' in l)
- rotloc_ref = []
- j = gl + 1
- while j < len(lines):
- if lines[j].strip() == '' or not lines[j].split()[0].lstrip('-').isdigit():
- j += 1
- continue
- # sort index line, then 3 krotm rows, then the Euler/iprop line.
- if len(lines[j].split()) == 1:
- krotm = np.array([[_to_float(x) for x in lines[j + 1 + r].split()]
- for r in range(3)])
- ang = lines[j + 4].split()
- a, b, g = (math.radians(float(x)) for x in ang[:3])
- iprop = int(ang[3])
- rotloc_ref.append(dict(krotm=krotm, a=a, b=b, g=g, iprop=iprop))
- j += 5
- else:
- j += 1
- return nsym, ops, rotloc_ref
-
-
-# --- almblm parsing ----------------------------------------------------------
-
-def _read_almblm(path, info):
- """Read one spin's almblm file. Keeps Alm, Blm and Clm (parproj needs the
- full radial vector), plus u_dot_norm, ovl_LO_u, ovl_LO_udot for the s12
- build and the per-band tetrahedron weights for the density matrix."""
- nsort = info['nsort']
- lmax = info['lmax']
- mult = info['mult']
- nlm = (lmax + 1) ** 2
- natom = sum(mult)
- nloat_max = 1 # max nLO over (l, sort); grown below
-
- r = _Reader(path)
- elecn = _to_float(r.record()[0])
- nk = int(r.record()[0])
- r.record() # nloat
- eferm = _to_float(r.record()[0])
-
- nLO = {}
- u_dot_norm = {}
- ovl_LO_u = {}
- ovl_LO_udot = {}
- kp = [None] * nk
-
- for isrt in range(1, nsort + 1):
- for l in range(lmax + 1):
- u_dot_norm[(l, isrt)] = _to_float(r.record()[0])
- 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
- nloat_max = max(nloat_max, n)
- for ilo in range(1, n + 1):
- toks = r.record()
- ovl_LO_u[(ilo, l, isrt)] = _to_float(toks[0])
- ovl_LO_udot[(ilo, l, isrt)] = _to_float(toks[1])
- 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,
- tetr=None,
- Alm=np.zeros((nlm, natom + 1, nb), dtype=complex),
- Blm=np.zeros((nlm, natom + 1, nb), dtype=complex),
- Clm=np.zeros((4, nlm, natom + 1, nb), dtype=complex))
- eband = np.empty(nb)
- tetr = np.empty(nb)
- for off in range(nb):
- toks = r.record()
- tetr[off] = _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'] = tetr[0]
- kp[ik]['tetr'] = tetr
- 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
- kp[ik]['Blm'][lm, iatom, off] = blm
- 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, u_dot_norm=u_dot_norm, ovl_LO_u=ovl_LO_u,
- ovl_LO_udot=ovl_LO_udot, kp=kp)
-
-
-# --- included orbital descriptors --------------------------------------------
-
-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,
- basis=sortinfo['basis']))
- return orbs
-
-
# --- s12 radial transform (set_projections.f:163-197) ------------------------
def _build_s12(l, isrt, n, spin):
@@ -385,82 +104,39 @@ def _srot_rotrep_nonmixing(op, l, transmat):
Returns (rotrep[2l+1,2l+1], timeinv, phase)."""
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
+ timeinv = _op_timeinv(op)
phase = (g - a) if timeinv else (a + g)
- rotl = _dmat(l, a, b, g, float(iprop))
+ rotl = dmat(l, a, b, g, float(iprop))
rotrep = transmat @ rotl @ np.conj(transmat.T)
if timeinv:
# timeinv_op in the new basis: reptrans T reptrans^T applied to conj.
- tmat = _tmat(l)
- tinv = transmat @ tmat @ transmat.T
+ tinv = transmat @ tmat(l) @ transmat.T
rotrep = tinv @ np.conj(rotrep)
return rotrep, timeinv, phase
# --- rotloc rotrep (set_rotloc.f) under SP+SO, non-mixing --------------------
-def _rotloc_rotrep_so(orb, ops, info, ref, transmat):
+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, non-mixing (set_rotloc.f), plus timeinv flag.
+ the new basis under SP+SO (set_rotloc.f), plus timeinv flag.
- ref is the representative-sort rotloc (krotm, a, b, g, iprop). set_rotloc
- runs the same composition for EVERY atom of the sort, including the
- representative: it finds the FIRST symmetry op R[isym] with
- perm(iref)==iatom and composes srot%rotl with the representative spinor
- rotloc. For the representative this op need not be the identity."""
+ ref is the representative-sort rotloc; rotloc_rotl_so composes the
+ representative spinor rotloc with the first symmetry op mapping the sort
+ representative onto this atom. The new-basis transform is the full transmat
+ for a mixing (spin-coupling) basis, and blkdiag(transmat, transmat) for a
+ spin-diagonal one."""
l = orb['l']
- isrt = orb['sort']
- iatom = orb['atom']
- iref = sum(info['mult'][:isrt - 1]) + 1
+ iref = sum(info['mult'][:orb['sort'] - 1]) + 1
d = 2 * l + 1
-
- # rotloc%rotl(2d) initial value (setsym.f:496-528): spmt (x) D(rotloc_ref).
- Dref = _dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
- f = (ref['a'] + ref['g']) / 2.0
- spmt = np.zeros((2, 2), dtype=complex)
- spmt[0, 0] = np.exp(1j * f) * math.cos(ref['b'] / 2.0)
- spmt[1, 1] = np.conj(spmt[0, 0])
- f = -(ref['a'] - ref['g']) / 2.0
- spmt[0, 1] = np.exp(1j * f) * math.sin(ref['b'] / 2.0)
- spmt[1, 0] = -np.conj(spmt[0, 1])
- rotl = np.zeros((2 * d, 2 * d), dtype=complex)
- rotl[:d, :d] = spmt[0, 0] * Dref
- rotl[d:, d:] = spmt[1, 1] * Dref
- rotl[:d, d:] = spmt[0, 1] * Dref
- rotl[d:, :d] = spmt[1, 0] * Dref
-
- # Compose with the first symmetry op mapping iref onto iatom (set_rotloc.f:60).
- isym = next(k for k, o in enumerate(ops)
- if o['perm'][iref - 1] == iatom)
- op = ops[isym]
- det2 = op['krotm'][0, 0] * op['krotm'][1, 1] - \
- op['krotm'][0, 1] * op['krotm'][1, 0]
- timeinv = det2 < 0.0
- srot_phase = (op['g'] - op['a']) if timeinv else (op['a'] + op['g'])
- # srot%rotl D(R[isym])_{lm}, with the timeinv operator applied if magnetic
- # (setsym.f:318-331 applies timeinv_op to srot%rotl before set_rotloc).
- rotl_sym = _dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
- if timeinv:
- rotl_sym = _tmat(l) @ np.conj(rotl_sym)
- ephase = np.exp(1j * srot_phase / 2.0)
- tmp = np.zeros((2 * d, 2 * d), dtype=complex)
- tmp[:d, :d] = ephase * rotl_sym
- tmp[d:, d:] = np.conj(ephase) * rotl_sym
- if timeinv:
- rotl = tmp @ np.conj(rotl)
+ rotl, timeinv = rotloc_rotl_so(l, ref, ops, orb['atom'], iref)
+ if mixing:
+ S = transmat
else:
- rotl = tmp @ rotl
-
- # rotrep = S rotl S^(H or T), S = blkdiag(transmat, transmat).
- S = np.zeros((2 * d, 2 * d), dtype=complex)
- S[:d, :d] = transmat
- S[d:, d:] = transmat
- if timeinv:
- rotrep = S @ rotl @ S.T
- else:
- rotrep = S @ rotl @ np.conj(S.T)
+ S = np.zeros((2 * d, 2 * d), dtype=complex)
+ S[:d, :d] = transmat
+ S[d:, d:] = transmat
+ rotrep = (S @ rotl @ S.T) if timeinv else (S @ rotl @ np.conj(S.T))
return rotrep, timeinv
@@ -502,6 +178,95 @@ def _build_matn_rep(orb, info, spins, windows, transmat, rot):
return out
+def _raw_matn(orb, sp, ik, nb_bot, nb_top, rot):
+ """The rot_projectmat'd (2l+1, nbsel, n) Theta projector in the |lm> basis
+ for one spin, before the angular transform (set_projections.f matn_rep)."""
+ l, atom, sort = orb['l'], orb['atom'], orb['sort']
+ kp = sp['kp'][ik]
+ off_bot = nb_bot - kp['nbmin']
+ off_top = nb_top - kp['nbmin']
+ nbsel = off_top - off_bot + 1
+ n = sp['nLO'][(l, sort)] + 2
+ s12 = _build_s12(l, sort, n, sp)
+ matn = np.zeros((2 * l + 1, nbsel, n), dtype=complex)
+ for mi, m in enumerate(range(-l, l + 1)):
+ lm = l * l + (m + l)
+ for j, off in enumerate(range(off_bot, off_top + 1)):
+ coeff = np.zeros(n, dtype=complex)
+ coeff[0] = kp['Alm'][lm, atom, off]
+ coeff[1] = kp['Blm'][lm, atom, off]
+ for ilo in range(n - 2):
+ coeff[2 + ilo] = kp['Clm'][ilo, lm, atom, off]
+ matn[mi, j, :] = coeff @ s12
+ for ir in range(n):
+ matn[:, :, ir] = rot @ matn[:, :, ir]
+ return matn
+
+
+def _build_matn_rep_mixing(orb, info, spins, windows, transmat, rot):
+ """matn_rep[ik] -> (2*(2l+1), nbsel, n): the mixing Theta projector. Each
+ spin's rot_projectmat'd |lm> block is stacked (up then dn) and multiplied by
+ the full 2(2l+1) transmat per radial channel (set_projections.f:515-585)."""
+ l = orb['l']
+ d = 2 * l + 1
+ out = {}
+ for ik in range(info['nk']):
+ incl, nb_bot, nb_top = windows[ik]
+ if not incl:
+ continue
+ up = _raw_matn(orb, spins[0], ik, nb_bot, nb_top, rot)
+ dn = _raw_matn(orb, spins[1], ik, nb_bot, nb_top, rot)
+ nbsel, n = up.shape[1], up.shape[2]
+ matn = np.zeros((2 * d, nbsel, n), dtype=complex)
+ matn[:d] = up
+ matn[d:] = dn
+ for ir in range(n):
+ matn[:, :, ir] = transmat @ matn[:, :, ir]
+ out[ik] = matn
+ return out
+
+
+def _orbital_densmat_mixing(orb, info, spins, windows, matn_rep):
+ """The single 2*(2l+1) raw density block for a mixing orbital, point-
+ integrated with the geometric k-weight (density.f:786-812): D = sum over k
+ and radial channels of matn_rep matn_rep^H, already in the new basis."""
+ l = orb['l']
+ d = 2 * l + 1
+ n = spins[0]['nLO'][(l, orb['sort'])] + 2
+ dens = np.zeros((2 * d, 2 * d), dtype=complex)
+ for ik in range(info['nk']):
+ incl, _, _ = windows[ik]
+ if not incl:
+ continue
+ weight = spins[0]['kp'][ik]['weight']
+ for i in range(n):
+ mat = matn_rep[(orb['atom'], orb['l'])][ik][:, :, i]
+ dens += (mat @ np.conj(mat.T)) * weight
+ return dens
+
+
+def _symmetrize_densmat_mixing(orbs, ops, nsym, rotreps, dens_raw):
+ """Scatter each mixing shell to its symmetry-permuted atom before averaging."""
+ out = {o['atom']: np.zeros_like(dens_raw[o['atom']]) for o in orbs}
+ for o in orbs:
+ atom = o['atom']
+ for isym, op in enumerate(ops):
+ rotrep, timeinv = rotreps[isym]
+ raw = dens_raw[atom]
+ tmp = np.conj(raw) if timeinv else raw
+ target = op['perm'][atom - 1]
+ out[target] += rotrep @ tmp @ rotrep.conj().T
+ return {atom: value / nsym for atom, value in out.items()}
+
+
+def _rotdens_densmat_mixing(blk, rotrep_loc, timeinv):
+ """rotdens_mat for a mixing orbital (rot_dens.f:119-141): inverse(Rloc) D
+ Rloc on the single 2*(2l+1) block."""
+ if timeinv:
+ return rotrep_loc.T @ np.conj(blk @ rotrep_loc)
+ return np.conj(rotrep_loc.T) @ (blk @ rotrep_loc)
+
+
# --- density matrix (density.f Theta path, non-mixing SP+SO) -----------------
def _orbital_densmat_blocks(orb, info, spins, windows, matn_reps):
@@ -529,8 +294,8 @@ def _orbital_densmat_blocks(orb, info, spins, windows, matn_reps):
weight = spins[isp]['kp'][ik]['weight']
acc = np.zeros((d, d), dtype=complex)
for i in range(n):
- mat = matn_reps[(orb['atom'], l)][(ik, isp)][:, :, i]
- cmat = matn_reps[(orb['atom'], l)][(ik, isp1)][:, :, i]
+ mat = matn_reps[(orb['atom'], orb['l'])][(ik, isp)][:, :, i]
+ cmat = matn_reps[(orb['atom'], orb['l'])][(ik, isp1)][:, :, i]
acc += mat @ np.conj(cmat.T)
blocks[iss - 1] += acc * weight
return blocks
@@ -601,47 +366,24 @@ def _rotdens_densmat(orb, blocks, rotrep_loc, timeinv):
# --- parproj writer ----------------------------------------------------------
-def _fmt(x):
- return ' %.16E' % float(x)
-
-
-def _w(f, arr):
- f.write(''.join(_fmt(x) for x in arr) + '\n')
-
-
def write_parproj(case):
"""Read .almblm{up,dn}, .indmftpr, .struct,
.dmftsym and write .parproj in the dmftproj format."""
- info = _read_indmftpr(case + '.indmftpr')
+ info = read_indmftpr(case + '.indmftpr')
if not info['so']:
- raise NotImplementedError('non-spin-orbit parproj is not implemented yet')
- if info['proj_mode'] != 0:
- raise NotImplementedError('parproj 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, rotloc_ref = _read_dmftsym(case + '.dmftsym')
+ 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')
- ns = 2
- spins = [_read_almblm(p, info) for p in spin_files]
+ spins = [read_almblm(p, info, projectors=True) for p in spin_files]
nk = spins[0]['nk']
info['nk'] = nk
orbs = _build_orbs(info)
- if any(o['basis'] == 'fromfile' for o in orbs):
- raise NotImplementedError(
- 'parproj for fromfile/mixing bases is not yet covered by a test '
- 'fixture; cubic and complex bases are supported')
- norb = len(orbs)
# Two band ranges (dmftproj.f:798,836): the energy window [e_bot, e_top]
# for the written Theta projector (block A), and the bands below e_bot
@@ -650,26 +392,29 @@ def write_parproj(case):
# density call is correlated-only), so densmat in outputqmc is that one.
windows = []
below_windows = []
+ bw = band_index_window(info, spins) if info['proj_mode'] else None
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']))
- below_windows.append(_select_window(kp['nbmin'], kp['nbmax'],
- kp['eband'], -1e6, info['e_bot']))
-
- for spin in spins[1:]:
- if spin['nk'] != nk:
- raise ValueError('up/down k-point counts differ')
- for ik, kp in enumerate(spin['kp']):
- for expected, lower, upper in ((windows[ik], info['e_bot'], info['e_top']),
- (below_windows[ik], -1e6, info['e_bot'])):
- if _select_window(kp['nbmin'], kp['nbmax'], kp['eband'], lower, upper) != expected:
- raise ValueError('spin-orbit run requires identical up/dn band windows')
- if not all(window[0] for window in windows):
- raise ValueError('parproj requires an included band window at every k-point')
-
- transmats = {(o['l'], o['sort']): _reptrans(o['basis'], o['l'])
- for o in orbs}
+ if bw is None:
+ windows.append(select_window(kp['nbmin'], kp['nbmax'], kp['eband'],
+ info['e_bot'], info['e_top']))
+ below_windows.append(select_window(kp['nbmin'], kp['nbmax'],
+ kp['eband'], -1e6, info['e_bot']))
+ else:
+ 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)
+
+ transmats = {}
+ mixing = {}
+ for o in orbs:
+ key = (o['l'], o['sort'])
+ if o['basis'] == 'fromfile':
+ transmats[key], mixing[key] = read_fromfile(
+ info['sorts'][o['sort'] - 1]['sourcefile'], o['l'])
+ else:
+ transmats[key], mixing[key] = reptrans(o['basis'], o['l']), False
# rot_projectmat local rotation: the op mapping the representative atom of
# the sort onto this atom (set_projections.f via rot_projectmat).
@@ -687,29 +432,17 @@ def write_parproj(case):
l = o['l']
transmat = transmats[(l, o['sort'])]
op = rotloc_op[o['atom']]
- rot = _dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
- matn_reps[(o['atom'], o['l'])] = _build_matn_rep(
- o, info, spins, windows, transmat, rot)
- matn_reps_full[(o['atom'], o['l'])] = _build_matn_rep(
- o, info, spins, below_windows, transmat, rot)
-
- # ---- srot rotrep per (sort, isym) for symmetrize_mat ----
- srot_rotreps = {}
- for isrt in range(1, info['nsort'] + 1):
- if not info['sorts'][isrt - 1]['included_ls']:
- continue
- for l in info['sorts'][isrt - 1]['included_ls']:
- transmat = transmats[(l, isrt)]
- for isym in range(nsym):
- srot_rotreps[(l, isrt, isym)] = _srot_rotrep_nonmixing(
- ops[isym], l, transmat)
-
- # ---- density matrices: raw -> symmetrize -> rotdens ----
- dens_raw = {}
- for o in orbs:
- dens_raw[(o['atom'], o['l'])] = _orbital_densmat_blocks(
- o, info, spins, below_windows, matn_reps_full)
- dens_sym = _symmetrize_densmat(orbs, ops, nsym, srot_rotreps, dens_raw)
+ rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+ if mixing[(l, o['sort'])]:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
+ o, info, spins, windows, transmat, rot)
+ matn_reps_full[(o['atom'], o['l'])] = _build_matn_rep_mixing(
+ o, info, spins, below_windows, transmat, rot)
+ else:
+ matn_reps[(o['atom'], o['l'])] = _build_matn_rep(
+ o, info, spins, windows, transmat, rot)
+ matn_reps_full[(o['atom'], o['l'])] = _build_matn_rep(
+ o, info, spins, below_windows, transmat, rot)
# ---- rotloc rotrep per orbital ----
rotloc_rotrep = {}
@@ -717,13 +450,47 @@ def write_parproj(case):
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)
+ 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.
+ nonmix_orbs = [o for o in orbs if not mixing[(o['l'], o['sort'])]]
densprint = {}
+
+ if nonmix_orbs:
+ srot_rotreps = {}
+ for o in nonmix_orbs:
+ 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)])
+ 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)
+
+ mixing_groups = {}
for o in 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)
+ key = (o['l'], o['sort'])
+ if mixing[key]:
+ mixing_groups.setdefault(key, []).append(o)
+ for (l, sort), group in mixing_groups.items():
+ rotreps = []
+ for op in ops:
+ ti = _op_timeinv(op)
+ rotreps.append((mixing_rotrep(op, l, transmats[(l, sort)], ti), ti))
+ raw = {o['atom']: _orbital_densmat_mixing(
+ o, info, spins, below_windows, matn_reps_full) for o in group}
+ sym = _symmetrize_densmat_mixing(group, ops, nsym, rotreps, raw)
+ for o in group:
+ key = (o['atom'], l)
+ rotrep_loc, timeinv = rotloc_rotrep[key]
+ densprint[key] = _rotdens_densmat_mixing(sym[o['atom']], rotrep_loc, timeinv)
# ---- write ----
with open(case + '.parproj', 'w') as f:
@@ -735,31 +502,40 @@ def write_parproj(case):
l = o['l']
atom = o['atom']
n = spins[0]['nLO'][(l, o['sort'])] + 2
+ ismix = mixing[(l, o['sort'])]
- # (A) Theta projector, non-mixing SP+SO (outputqmc.f:955-973).
+ # (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.
for ik in range(nk):
incl, nb_bot, nb_top = windows[ik]
for ir in range(n):
- for ispin in range(ns):
+ if ismix:
+ P = matn_reps[(atom, l)][ik][:, :, ir]
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].real)
+ for m in range(2 * (2 * l + 1)):
+ write_row(f, P[m, :].imag)
+ continue
+ for ispin in range(2):
P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
for mi in range(2 * l + 1):
- _w(f, P[mi, :].real)
- for ispin in range(ns):
+ write_row(f, P[mi, :].real)
+ for ispin in range(2):
P = matn_reps[(atom, l)][(ik, ispin)][:, :, ir]
for mi in range(2 * l + 1):
- _w(f, P[mi, :].imag)
+ write_row(f, P[mi, :].imag)
- # (B) density matrix, non-mixing SP+SO 2*(2l+1) (outputqmc.f:1033-1049).
+ # (B) density matrix, SP+SO 2*(2l+1) (outputqmc.f:1012-1066).
dp = densprint[(atom, l)]
for m in range(2 * (2 * l + 1)):
- _w(f, dp[m, :].real)
+ write_row(f, dp[m, :].real)
for m in range(2 * (2 * l + 1)):
- _w(f, dp[m, :].imag)
+ write_row(f, dp[m, :].imag)
- # (C) Rloc rotrep, non-mixing SP+SO (outputqmc.f:1148-1158).
+ # (C) Rloc rotrep, SP+SO 2*(2l+1) (outputqmc.f:1130-1177).
rotrep_loc, timeinv = rotloc_rotrep[(atom, l)]
for m in range(2 * (2 * l + 1)):
- _w(f, rotrep_loc[m, :].real)
+ write_row(f, rotrep_loc[m, :].real)
for m in range(2 * (2 * l + 1)):
- _w(f, rotrep_loc[m, :].imag)
+ write_row(f, rotrep_loc[m, :].imag)
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 8e37682..493fb2d 100644
--- a/python/triqs_dftkit/wien2k/sympar.py
+++ b/python/triqs_dftkit/wien2k/sympar.py
@@ -31,10 +31,8 @@
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).
+identity, a cubic shell uses the exact double-precision cubic-d harmonics, as
+ctqmcout does. Mixed within one file (Ca complex, Os cubic 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
@@ -42,216 +40,58 @@
timeflag line, and the representation matrices.
Only the spin-orbit case is implemented (indmftpr SO flag=1, which implies
-ifSP): the non-mixing 2*(2l+1)-wide block-diag spinor path 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).
+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).
"""
-import math
-import numpy as np
+from ._dmftproj import (read_dmftsym, read_fromfile, read_indmftpr, reptrans,
+ shell_matrix, write_row)
+_shell_matrix = shell_matrix
-# --- 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),
-}
+# --- included shells ---------------------------------------------------------
-
-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
+def _included_shells(info):
+ """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, a mixing
+ flag and the transform P: a spin-coupling fromfile basis gives the full
+ 2(2l+1) P (P spinrot P^dag matrix), otherwise the (2l+1) up/up block."""
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
+ for isort in range(info['nsort']):
+ s = info['sorts'][isort]
+ for l in s['included_ls']:
+ if s['basis'] == 'fromfile':
+ P, mixing = read_fromfile(s['sourcefile'], l)
+ else:
+ P, mixing = reptrans(s['basis'], l), False
+ for imu in range(1, info['mult'][isort] + 1):
+ atom = sum(info['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)
+ basis=s['basis'], P=P, mixing=mixing))
+ return shells
# --- 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)
+ write_row(f, mat[m, :].real)
for m in range(mat.shape[0]):
- _w(f, mat[m, :].imag)
+ write_row(f, mat[m, :].imag)
def write_sympar(case):
"""Write .sympar from .dmftsym and .indmftpr.
Spin-orbit only (indmftpr SO flag=1, implying ifSP): writes the timeflag
- line and 2*(2l+1)-wide block-diag spinor matrices."""
- shells, so = _read_indmftpr(case + '.indmftpr')
- if not so:
+ line and the 2*(2l+1)-wide spinor matrices."""
+ info = read_indmftpr(case + '.indmftpr')
+ if not info['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')
+ shells = _included_shells(info)
+ nsym, ops = read_dmftsym(case + '.dmftsym')
natom = len(ops[0]['perm'])
timeflag = []
diff --git a/python/triqs_dftkit/wien2k/symqmc.py b/python/triqs_dftkit/wien2k/symqmc.py
index bcc85ca..524d9e0 100644
--- a/python/triqs_dftkit/wien2k/symqmc.py
+++ b/python/triqs_dftkit/wien2k/symqmc.py
@@ -41,170 +41,55 @@
time-reversal operator -i sigma_y (x) T applied to the magnetic operations;
- l = 0 (s) shells: the 2x2 spin phase block.
-The mixing fromfile path is the one dft_tools #148 singles out. dmftproj reads
-both cubic and fromfile bases with single-precision CMPLX casts and a 250-column
-line cap (Fortran set_ang_trans.f), and evaluates spin phases in single precision
-(outputqmc.f, setsym.f). Its case.symqmc therefore differs at about 1e-7 from
-this double-precision generator, whose cubic d basis uses exact coefficients.
+All matrices are built in double precision from the exact cubic harmonics.
+A case.symqmc written by a dmftproj that casts the basis or the spin phases
+through single-precision CMPLX differs from this output at about 1e-7.
"""
-import math
-import os
import numpy as np
-# --- angular bases (transpose(P) = , m = -l..l) -----------------------
+from ._dmftproj import (read_dmftsym, read_fromfile, read_indmftpr, reptrans,
+ shell_matrix)
-_COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
+_shell_matrix = shell_matrix
-# standard cubic harmonics, Wien2k convention
-_CUBIC = {
- 1: np.array([[0, 1, 0], [-1j, 0, -1j], [1, 0, -1]], dtype=complex) / 1.0,
- 2: np.array([
- [0, 0, 1, 0, 0],
- [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
- [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5],
- [0, 2 ** -0.5, 0, -(2 ** -0.5), 0],
- [0, 2 ** -0.5, 0, 2 ** -0.5, 0],
- ], dtype=complex),
-}
+# --- correlated shells -------------------------------------------------------
-def _reptrans(basis, l):
- if basis == 'cubic' and l in _CUBIC:
- return _CUBIC[l]
- return _COMPLEX[l]
-
-
-# --- Wigner D matrix (dmftproj convention, setsym.f) -------------------------
-
-def _small_d(l, m, n, b):
- f1 = (math.factorial(l + m) * math.factorial(l - m)) / \
- (math.factorial(l + n) * math.factorial(l - n))
- s = 0.0
- for t in range(0, 2 * l + 1):
- if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0:
- f2 = (math.factorial(l + n) * math.factorial(l - n)) / \
- (math.factorial(l - m - t) * math.factorial(m + n + t) *
- math.factorial(l - n - t) * math.factorial(t))
- f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 else math.sin(b / 2) ** (2 * l - m - n - 2 * t)
- f4 = 1.0 if (2 * t + n + m) == 0 else math.cos(b / 2) ** (2 * t + n + m)
- s += (-1) ** (l - m - t) * f2 * f3 * f4
- return math.sqrt(f1) * s
-
-
-def _dmat(l, a, b, c, det):
- D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
- for m in range(-l, l + 1):
- for n in range(-l, l + 1):
- v = np.exp(1j * n * a) * np.exp(1j * m * c) * _small_d(l, m, n, b)
- if det < -0.5:
- v *= (-1) ** l
- D[m + l, n + l] = v
- return D
-
-
-def _tmat(l):
- """Complex-conjugation operator in the spherical-harmonic basis,
- T[m, -m] = (-1)^m (timeinv.f)."""
- T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex)
- for m in range(-l, l + 1):
- T[-m + l, m + l] = (-1) ** m
- return T
-
-
-def _timeinv_orbital(l, mat):
- return _tmat(l) @ np.conj(mat)
-
-
-# --- fromfile angular basis ---------------------------------------------------
-
-def _read_fromfile(path, l):
- """Parse a dmftproj fromfile basis (each line: the coefficients of a new
- basis vector in {|m,up>, |m,dn>}, m = -l..l, real/imag interleaved; '*'
- marks the end of an irep). Return (P, mixing) where P = is the
- transform matrix (reptrans.transmat), full 2(2l+1) for a spin-mixing basis
- or the (2l+1) up/up block otherwise."""
- n = 2 * (2 * l + 1)
- rows = []
- for line in open(path):
- line = line.rstrip('\n')
- if not line.strip():
- continue
- body = line[1:]
- vals = [float(x) for x in body.split()][:2 * n]
- rows.append([vals[2 * k] + 1j * vals[2 * k + 1] for k in range(n)])
- if len(rows) == n:
- break
- R = np.array(rows) # R[i, :] = |new_i> in old basis
- d = 2 * l + 1
- up_up, up_dn = R[:d, :d], R[:d, d:]
- dn_up, dn_dn = R[d:, :d], R[d:, d:]
- mixing = not (np.allclose(dn_dn, up_up) and
- np.allclose(up_dn, 0) and np.allclose(dn_up, 0))
- P = np.conj(R) if mixing else np.conj(up_up) # = conj()
- return P, mixing
-
+def _correlated_shells(info):
+ """One entry per correlated atom (l_inc==2), in sort order. Each shell
+ carries l, basis name, its transform matrix P = and a mixing flag
+ (True for a spin-coupling fromfile basis). symqmc builds its symmetry
+ matrices from the exact cubic harmonics."""
+ shells = []
+ for isort in range(info['nsort']):
+ s = info['sorts'][isort]
+ basis = s['basis']
+ for l in s['correlated_ls']:
+ if basis == 'fromfile':
+ P, mixing = read_fromfile(s['sourcefile'], l)
+ else:
+ P, mixing = reptrans(basis, l), False
+ for _ in range(info['mult'][isort]):
+ shells.append(dict(l=l, basis=basis, P=P, mixing=mixing))
+ return shells
-# --- input parsing -----------------------------------------------------------
-def _read_dmftsym(path):
- lines = open(path).read().split('\n')
- nsym = int(lines[0].split()[0])
- perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)]
- rest = lines[1 + nsym:]
- starts = [i for i, l in enumerate(rest) if 'Sym. op.' in l]
- ops = []
- for k, s in enumerate(starts):
- ang = rest[s + 1].split()
- a, b, c = (math.radians(float(x)) for x in ang[:3])
- krotm = np.array([[float(x) for x in rest[s + 2 + r].split()] for r in range(3)])
- ops.append(dict(perm=perms[k], a=a, b=b, c=c, krotm=krotm))
- return nsym, ops
+# The symqmc test reaches into these two names directly; keep them as the
+# module's parsing entry points.
+_read_dmftsym = read_dmftsym
def _read_correlated_shells(indmftpr, struct):
- """Return the list of correlated shells, one entry per correlated atom, plus
- the SO flag, from case.indmftpr and case.struct multiplicities. Each shell
- carries l, the basis name, its transform matrix P = and a mixing
- flag (True for a spin-coupling fromfile basis)."""
- raw = [l.split('!')[0].strip() for l in open(indmftpr)]
- raw = [l for l in raw if l != '']
- nsort = int(raw[0].split()[0])
- mult = [int(x) for x in raw[1].split()][:nsort]
- i = 3
- so = 0
- shells = []
- for isort in range(nsort):
- basis = raw[i].split()[0]
- i += 1
- sourcefile = None
- if basis == 'fromfile':
- sourcefile = os.path.join(os.path.dirname(indmftpr), raw[i])
- i += 1
- l_inc = [int(x) for x in raw[i].split()]
- i += 1
- ireps = [int(x) for x in raw[i].split()]
- i += 1
- correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2]
- if any(n > 0 for n in ireps):
- i += 1 # skip the correps line
- if correlated_ls:
- so = int(raw[i].split()[0]) # SO flag follows a correlated sort
- i += 1
- for l in correlated_ls:
- if basis == 'fromfile':
- P, mixing = _read_fromfile(sourcefile, l)
- else:
- P, mixing = _reptrans(basis, l), False
- for _ in range(mult[isort]):
- shells.append(dict(l=l, basis=basis, P=P, mixing=mixing))
- return shells, so
+ info = read_indmftpr(indmftpr)
+ return _correlated_shells(info), info['so']
def write_symqmc(case):
"""Write .symqmc from .dmftsym, .indmftpr, .struct."""
- nsym, ops = _read_dmftsym(case + '.dmftsym')
- shells, so = _read_correlated_shells(case + '.indmftpr', case + '.struct')
+ nsym, ops = read_dmftsym(case + '.dmftsym')
+ info = read_indmftpr(case + '.indmftpr')
+ shells, so = _correlated_shells(info), info['so']
natom = len(ops[0]['perm'])
timeinv = []
@@ -223,69 +108,6 @@ def write_symqmc(case):
f.write(_format_matrix(_shell_matrix(op, sh, timeinv[isym])))
-def _shell_matrix(op, shell, ti):
- """Spinor symmetry matrix for one correlated shell under one operation."""
- l = shell['l']
- if l == 0:
- return _l0_matrix(op, ti)
- if shell['mixing']:
- return _mixing_matrix(op, shell, ti)
- return _nonmixing_matrix(op, shell, ti)
-
-
-def _phase(op, ti):
- a, c = op['a'], op['c']
- return (c - a) if ti else (a + c) # (g-a) on magnetic ops, else (a+g)
-
-
-def _l0_matrix(op, ti):
- """s shell: 2x2 spin phase block diag(e, conj(e)) (outputqmc.f l==0)."""
- e = np.exp(1j * _phase(op, ti) / 2)
- return np.array([[e, 0], [0, np.conj(e)]], dtype=complex)
-
-
-def _nonmixing_matrix(op, shell, ti):
- """Spin-diagonal basis: the up/up block scaled by +-(a+g)/2, with the
- orbital time-reversal operator on the magnetic operations."""
- l, P = shell['l'], shell['P']
- rotl = _dmat(l, op['a'], op['b'], op['c'], np.linalg.det(op['krotm']))
- if ti:
- rotl = _timeinv_orbital(l, rotl)
- rotrep = P @ rotl @ np.conj(P.T)
- e = np.exp(1j * _phase(op, ti) / 2)
- d = 2 * l + 1
- mat = np.zeros((2 * d, 2 * d), dtype=complex)
- mat[:d, :d] = e * rotrep
- mat[d:, d:] = np.conj(e) * rotrep
- return mat
-
-
-def _mixing_matrix(op, shell, ti):
- """Spin-coupling basis: full 2(2l+1) spinor representation. The orbital
- rotation enters block-diagonal (beta=0) or block-antidiagonal (beta=pi, the
- magnetic operations); the spinor time-reversal -i sigma_y (x) T is then
- applied to the magnetic operations (setsym.f, timeinv.f)."""
- l, P = shell['l'], shell['P']
- rotl = _dmat(l, op['a'], op['b'], op['c'], np.linalg.det(op['krotm']))
- e = np.exp(1j * _phase(op, ti) / 2)
- d = 2 * l + 1
- spinrot = np.zeros((2 * d, 2 * d), dtype=complex)
- if ti: # beta = pi, block-antidiagonal
- spinrot[:d, d:] = e * rotl
- spinrot[d:, :d] = -np.conj(e) * rotl
- else: # beta = 0, block-diagonal
- spinrot[:d, :d] = e * rotl
- spinrot[d:, d:] = np.conj(e) * rotl
- rotrep = P @ spinrot @ np.conj(P.T)
- if ti:
- tm = _tmat(l)
- tinv = np.zeros((2 * d, 2 * d), dtype=complex)
- tinv[:d, d:] = -tm
- tinv[d:, :d] = tm
- rotrep = (P @ tinv @ P.T) @ np.conj(rotrep)
- return rotrep
-
-
def _format_matrix(mat):
out = []
for part in (mat.real, mat.imag):
diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt
index 1d7f454..f94fffa 100644
--- a/test/python/wien2k/CMakeLists.txt
+++ b/test/python/wien2k/CMakeLists.txt
@@ -37,9 +37,8 @@ add_test(NAME Py_wien2k_symqmc_python
set_property(TEST Py_wien2k_symqmc_python APPEND PROPERTY ENVIRONMENT
PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
-# Generator on the two paths the cubic test does not reach: the spin-mixing
-# fromfile (|j,m_j>) basis and an l=0 shell. Reuses CaOs2.dmftsym/struct; the
-# references are the dmftproj matrix data (single precision, ~26 kB total).
+# Native double-precision references for spin-mixing fromfile bases, an l=0
+# shell and cubic p harmonics, including magnetic operations.
add_test(NAME Py_wien2k_symqmc_mixing
COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symqmc_mixing.py
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
@@ -80,3 +79,15 @@ add_test(NAME Py_wien2k_kpath_labels
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
set_property(TEST Py_wien2k_kpath_labels APPEND PROPERTY ENVIRONMENT
PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+add_test(NAME Py_wien2k_dmftproj_common_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_dmftproj_common_python.py
+ WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_dmftproj_common_python APPEND PROPERTY ENVIRONMENT
+ PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})
+
+add_test(NAME Py_wien2k_outband_python
+ COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_outband_python.py
+ 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})
diff --git a/test/python/wien2k/CaOs2.README.md b/test/python/wien2k/CaOs2.README.md
index 8f6ed1b..53144ef 100644
--- a/test/python/wien2k/CaOs2.README.md
+++ b/test/python/wien2k/CaOs2.README.md
@@ -93,3 +93,17 @@ The Python generator uses exact cubic d coefficients and double-precision
phases; its generated matrices are checked for unitarity at `1e-12`.
It writes only into a temporary directory, so the converter tests always
read the original native fixture.
+
+## Double-precision generator references
+
+The separate generator comparisons in `wien2k_ctqmcout_python.py`,
+`wien2k_parproj_python.py`, `wien2k_sympar_python.py` and
+`wien2k_outband_python.py` use the precision-corrected native producer.
+Their cubic-d references require a full-precision `case.cf_d_eg_t2g` in
+`$WIENROOT/SRC_templates`. Replace the eight-digit `0.70710678` coefficients
+with `0.70710678118654752440` (preserving their signs and all other entries)
+before regenerating those references. The Python cubic basis uses the same
+double-precision coefficients. The eight-digit template can leave differences
+of about `7e-9`, exceeding the generator tests' absolute tolerance of `1e-11`.
+This requirement applies to the full-rank generator references; retain the
+original template when reproducing the historical converter fixture above.
diff --git a/test/python/wien2k/CaOs2_band.almblmdn.gz b/test/python/wien2k/CaOs2_band.almblmdn.gz
new file mode 100644
index 0000000..0329e74
Binary files /dev/null and b/test/python/wien2k/CaOs2_band.almblmdn.gz differ
diff --git a/test/python/wien2k/CaOs2_band.almblmup.gz b/test/python/wien2k/CaOs2_band.almblmup.gz
new file mode 100644
index 0000000..3eccf62
Binary files /dev/null and b/test/python/wien2k/CaOs2_band.almblmup.gz differ
diff --git a/test/python/wien2k/CaOs2_band.indmftpr b/test/python/wien2k/CaOs2_band.indmftpr
new file mode 100644
index 0000000..6aa6f67
--- /dev/null
+++ b/test/python/wien2k/CaOs2_band.indmftpr
@@ -0,0 +1,13 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+fromfile
+jbasis_d.dat
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0
+0.6529001225
diff --git a/test/python/wien2k/CaOs2_band.klist_band b/test/python/wien2k/CaOs2_band.klist_band
new file mode 100644
index 0000000..366e653
--- /dev/null
+++ b/test/python/wien2k/CaOs2_band.klist_band
@@ -0,0 +1,2 @@
+GAMMA 0 0 0 8 2.0
+END
diff --git a/test/python/wien2k/CaOs2_band.outband.gz b/test/python/wien2k/CaOs2_band.outband.gz
new file mode 100644
index 0000000..063276c
Binary files /dev/null and b/test/python/wien2k/CaOs2_band.outband.gz differ
diff --git a/test/python/wien2k/CaOs2_band_cubic.indmftpr b/test/python/wien2k/CaOs2_band_cubic.indmftpr
new file mode 100644
index 0000000..7452996
--- /dev/null
+++ b/test/python/wien2k/CaOs2_band_cubic.indmftpr
@@ -0,0 +1,12 @@
+2
+1 2
+3
+complex
+0 0 1 0
+0 0 0 0
+cubic
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0
+0.6529001225
diff --git a/test/python/wien2k/CaOs2_band_cubic.outband.gz b/test/python/wien2k/CaOs2_band_cubic.outband.gz
new file mode 100644
index 0000000..2201443
Binary files /dev/null and b/test/python/wien2k/CaOs2_band_cubic.outband.gz differ
diff --git a/test/python/wien2k/CaOs2_full.ctqmcout.gz b/test/python/wien2k/CaOs2_full.ctqmcout.gz
new file mode 100644
index 0000000..3cca2cc
Binary files /dev/null and b/test/python/wien2k/CaOs2_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_full.indmftpr b/test/python/wien2k/CaOs2_full.indmftpr
new file mode 100644
index 0000000..227d144
--- /dev/null
+++ b/test/python/wien2k/CaOs2_full.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
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_jbasis_full.ctqmcout.gz b/test/python/wien2k/CaOs2_jbasis_full.ctqmcout.gz
new file mode 100644
index 0000000..c3c7a7c
Binary files /dev/null and b/test/python/wien2k/CaOs2_jbasis_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_jbasis_full.indmftpr b/test/python/wien2k/CaOs2_jbasis_full.indmftpr
new file mode 100644
index 0000000..dab1b92
--- /dev/null
+++ b/test/python/wien2k/CaOs2_jbasis_full.indmftpr
@@ -0,0 +1,12 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+fromfile
+jbasis_d.dat
+0 0 2 0
+0 0 0 0
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_jbasis_full.parproj.gz b/test/python/wien2k/CaOs2_jbasis_full.parproj.gz
new file mode 100644
index 0000000..9a60325
Binary files /dev/null and b/test/python/wien2k/CaOs2_jbasis_full.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_jbasis_full.sympar.gz b/test/python/wien2k/CaOs2_jbasis_full.sympar.gz
new file mode 100644
index 0000000..a9abae8
Binary files /dev/null and b/test/python/wien2k/CaOs2_jbasis_full.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_mode1.indmftpr b/test/python/wien2k/CaOs2_mode1.indmftpr
new file mode 100644
index 0000000..1cabe50
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1.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
+1
+-0.15 0.30 1
diff --git a/test/python/wien2k/CaOs2_mode1.oubwindn b/test/python/wien2k/CaOs2_mode1.oubwindn
new file mode 100644
index 0000000..97fdd85
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1.oubwindn
@@ -0,0 +1,5 @@
+ 1
+ 1
+ 1
+ 60 76
+ 1.0000000000000000
diff --git a/test/python/wien2k/CaOs2_mode1.oubwinup b/test/python/wien2k/CaOs2_mode1.oubwinup
new file mode 100644
index 0000000..97fdd85
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1.oubwinup
@@ -0,0 +1,5 @@
+ 1
+ 1
+ 1
+ 60 76
+ 1.0000000000000000
diff --git a/test/python/wien2k/CaOs2_mode1_full.ctqmcout.gz b/test/python/wien2k/CaOs2_mode1_full.ctqmcout.gz
new file mode 100644
index 0000000..38b1ece
Binary files /dev/null and b/test/python/wien2k/CaOs2_mode1_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_mode1_full.indmftpr b/test/python/wien2k/CaOs2_mode1_full.indmftpr
new file mode 100644
index 0000000..13e6507
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode1_full.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
+1
+-2.0 3.0 1
diff --git a/test/python/wien2k/CaOs2_mode2_below.ctqmcout.gz b/test/python/wien2k/CaOs2_mode2_below.ctqmcout.gz
new file mode 100644
index 0000000..49d68ac
Binary files /dev/null and b/test/python/wien2k/CaOs2_mode2_below.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_mode2_below.indmftpr b/test/python/wien2k/CaOs2_mode2_below.indmftpr
new file mode 100644
index 0000000..0e008e1
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode2_below.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
+1
+0 92 2
diff --git a/test/python/wien2k/CaOs2_mode2_full.ctqmcout.gz b/test/python/wien2k/CaOs2_mode2_full.ctqmcout.gz
new file mode 100644
index 0000000..38b1ece
Binary files /dev/null and b/test/python/wien2k/CaOs2_mode2_full.ctqmcout.gz differ
diff --git a/test/python/wien2k/CaOs2_mode2_full.indmftpr b/test/python/wien2k/CaOs2_mode2_full.indmftpr
new file mode 100644
index 0000000..7732eb7
--- /dev/null
+++ b/test/python/wien2k/CaOs2_mode2_full.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
+1
+43 92 2
diff --git a/test/python/wien2k/CaOs2_partial.parproj.gz b/test/python/wien2k/CaOs2_partial.parproj.gz
index 7b893c4..6507818 100644
Binary files a/test/python/wien2k/CaOs2_partial.parproj.gz and b/test/python/wien2k/CaOs2_partial.parproj.gz differ
diff --git a/test/python/wien2k/CaOs2_partial.sympar.gz b/test/python/wien2k/CaOs2_partial.sympar.gz
index d553f10..e4e2f53 100644
Binary files a/test/python/wien2k/CaOs2_partial.sympar.gz and b/test/python/wien2k/CaOs2_partial.sympar.gz differ
diff --git a/test/python/wien2k/CaOs2_pcubic.indmftpr b/test/python/wien2k/CaOs2_pcubic.indmftpr
new file mode 100644
index 0000000..2cb0501
--- /dev/null
+++ b/test/python/wien2k/CaOs2_pcubic.indmftpr
@@ -0,0 +1,11 @@
+2
+1 2
+3
+complex
+0 0 0 0
+0 0 0 0
+cubic
+0 2 0 0
+0 0 0 0
+1
+-2.0 3.0
diff --git a/test/python/wien2k/CaOs2_pcubic.symqmc.gz b/test/python/wien2k/CaOs2_pcubic.symqmc.gz
new file mode 100644
index 0000000..ce1b94c
Binary files /dev/null and b/test/python/wien2k/CaOs2_pcubic.symqmc.gz differ
diff --git a/test/python/wien2k/wien2k_ctqmcout_python.py b/test/python/wien2k/wien2k_ctqmcout_python.py
index 83f5fe1..2e3a4c5 100644
--- a/test/python/wien2k/wien2k_ctqmcout_python.py
+++ b/test/python/wien2k/wien2k_ctqmcout_python.py
@@ -4,15 +4,27 @@
################################################################################
# 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.
+# reproduces the dmftproj Fortran correlated-shell projectors to machine
+# precision: regenerate case.ctqmcout from the almblm/struct/dmftsym + a wide
+# energy window and compare to the native 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).
+# The window (-2.0, 3.0) keeps 50 bands, more than the 20 correlated spin-orbitals
+# of the two SOC Os d shells, so the Loewdin overlap O is full rank. The
+# references are native dmftproj output with double-precision casts and a
+# full-precision cubic-d template.
+#
+# CaOs2_jbasis_full exercises the spin-mixing fromfile (|j, m_j>) path: the Os d
+# shell is a 2(2l+1)=10 spinor basis (jbasis_d.dat), so the projector, Rloc
+# rotrep and complex-harmonics transform blocks are the full mixing matrices.
+# Same wide window for the full-rank Loewdin.
+#
+# CaOs2_mode1_full and CaOs2_mode2_full exercise the band-index projection modes
+# (set_projections.f:70-88). Mode 1 ("-2.0 3.0 1") scans all spins/k for the
+# global band-index window; mode 2 ("43 92 2") takes the indices straight from
+# the window line. Both resolve to bands 43..92 (50 bands > 20 correlated
+# spin-orbitals), so the Loewdin overlap stays full rank and the projectors
+# match the native dmftproj to machine precision. CaOs2_mode2_below ("0 92 2")
+# puts b_bot below nbmin, so no bands lie below the window and qbbot is zero.
import gzip
import os
@@ -30,9 +42,9 @@ def _gunzip(src, dst):
shutil.copyfileobj(fi, fo)
-def _compare(got_path, ref_path):
+def _compare(got_path, ref_lines):
got = open(got_path).read().split()
- ref = open(ref_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):
@@ -41,59 +53,50 @@ def _compare(got_path, ref_path):
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)
+ np.testing.assert_allclose(a, b, atol=1e-11, 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}'),
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ if os.path.exists(os.path.join(HERE, 'jbasis_d.dat')):
+ shutil.copy(os.path.join(HERE, 'jbasis_d.dat'),
+ os.path.join(tmp, 'jbasis_d.dat'))
+ 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'{case}.{ext}.gz'),
+ _gunzip(os.path.join(HERE, f'CaOs2.{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')
+ with gzip.open(os.path.join(HERE, f'{case}.ctqmcout.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.ctqmcout'), fh.readlines())
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')
+_check('CaOs2_full')
+_check('CaOs2_jbasis_full')
+_check('CaOs2_mode1_full') # proj_mode 1: band-index window 43..92
+_check('CaOs2_mode2_full') # proj_mode 2: explicit band indices 43 92
+_check('CaOs2_mode2_below') # proj_mode 2 with b_bot below nbmin
-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)
+print('wien2k_ctqmcout_python: ok')
+
+# A singular overlap must fail clearly rather than writing NaN projectors.
try:
- ctqmcout._sqrt_inv(np.diag([0., 1.]).astype(complex))
+ ctqmcout._sqrt_inv(np.ones((2, 2), dtype=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')
+# Rank is invariant under a uniform projector scaling. Small, well-conditioned
+# overlaps must remain valid and normalize to the identity.
+for scale in (1., 1e-20):
+ overlap = np.diag([scale, 2 * scale]).astype(complex)
+ inverse_root = ctqmcout._sqrt_inv(overlap)
+ np.testing.assert_allclose(inverse_root @ overlap @ inverse_root.conj().T,
+ np.eye(2), atol=1e-14, rtol=0)
diff --git a/test/python/wien2k/wien2k_dmftproj_common_python.py b/test/python/wien2k/wien2k_dmftproj_common_python.py
new file mode 100644
index 0000000..50c6770
--- /dev/null
+++ b/test/python/wien2k/wien2k_dmftproj_common_python.py
@@ -0,0 +1,176 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Unit tests for the shared dmftproj machinery (triqs_dftkit.wien2k._dmftproj),
+# the module the five case.* generators are built on. The generators' golden
+# tests cover the assembled output; these exercise the isolated primitives
+# directly so a regression in a shared unit is caught at its source.
+
+import os
+import tempfile
+
+import numpy as np
+
+from triqs_dftkit.wien2k import _dmftproj as dp
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+
+
+# --- list-directed numeric parsing -------------------------------------------
+
+assert dp.to_float('1.0D-3') == 1.0e-3
+assert dp.to_float('4.57E-002') == 4.57e-2
+assert dp.to_float('-2.5d0') == -2.5
+
+assert dp.to_complex('(1.0,-2.0)') == complex(1.0, -2.0)
+
+
+def _reader(text):
+ fh = tempfile.NamedTemporaryFile('w', suffix='.tmp', delete=False)
+ fh.write(text)
+ fh.close()
+ return dp.Reader(fh.name)
+
+
+r = _reader('(1.0,2.0)\n(3.0,-4.0) (5.0,6.0)\n')
+assert dp.read_complex(r) == complex(1.0, 2.0)
+a, b = dp.read_two_complex(r)
+assert a == complex(3.0, -4.0) and b == complex(5.0, 6.0)
+
+
+# --- angular basis: cubic transform -------------------------------------------
+
+exact = dp.reptrans('cubic', 2)
+assert exact.shape == (5, 5)
+# unitary
+assert np.max(np.abs(exact @ np.conj(exact.T) - np.eye(5))) < 1e-12
+assert abs(abs(exact[1, 0]) - 2 ** -0.5) < 1e-15
+assert dp.reptrans('complex', 2).shape == (5, 5)
+assert np.allclose(dp.reptrans('complex', 2), np.eye(5))
+
+
+# --- Wigner D and orbital time reversal --------------------------------------
+
+for l in (1, 2):
+ D0 = dp.dmat(l, 0.0, 0.0, 0.0, 1.0)
+ assert np.max(np.abs(D0 - np.eye(2 * l + 1))) < 1e-12 # zero rotation
+ D = dp.dmat(l, 0.3, 0.7, 1.1, 1.0)
+ assert np.max(np.abs(D @ np.conj(D.T) - np.eye(2 * l + 1))) < 1e-12 # unitary
+ T = dp.tmat(l)
+ for m in range(-l, l + 1):
+ assert T[-m + l, m + l] == (-1) ** m
+
+
+# --- select_window: contiguous band range in (e1, e2] ------------------------
+
+eband = np.array([-1.0, -0.1, 0.05, 0.2, 0.5]) # bands 10..14, window (-0.15, 0.30]
+incl, lo, hi = dp.select_window(10, 14, eband, -0.15, 0.30)
+assert incl and lo == 11 and hi == 13
+
+
+# --- band-index window (proj_mode 1/2) ---------------------------------------
+
+# set_projections.f:70-88: indices clamped to nbmin/nbmax.
+incl, lo, hi = dp.select_band_window(1, 92, 60, 76)
+assert incl and lo == 60 and hi == 76
+incl, lo, hi = dp.select_band_window(50, 80, 43, 92) # clamp both ends
+assert incl and lo == 50 and hi == 80
+# the below-window (1, b_bot - 1) is empty once b_bot <= nbmin
+assert not dp.select_band_window(5, 92, 1, 4)[0]
+assert not dp.select_band_window(1, 92, 1, -1)[0]
+
+# proj_mode 2 takes b_bot/b_top straight from the (rounded) window line.
+m2 = {'proj_mode': 2, 'e_bot': 43.0, 'e_top': 92.0}
+assert dp.band_index_window(m2, []) == (43, 92)
+
+# proj_mode 1 scans the (e_bot, e_top] energies for the global band-index range.
+m1 = {'proj_mode': 1, 'e_bot': -0.15, 'e_top': 0.30}
+spins = [{'kp': [{'nbmin': 1, 'nbmax': 5,
+ 'eband': np.array([-1.0, -0.1, 0.05, 0.2, 0.5])}]}]
+assert dp.band_index_window(m1, spins) == (2, 4)
+
+
+# --- case.indmftpr / case.dmftsym structured parse ---------------------------
+
+info = dp.read_indmftpr(os.path.join(HERE, 'CaOs2.indmftpr'))
+assert info['nsort'] == 2 and info['mult'] == [1, 2] and info['lmax'] == 3
+assert info['so'] == 1
+assert info['sorts'][1]['correlated_ls'] == [2] # Os d correlated
+assert info['sorts'][0]['correlated_ls'] == [] # Ca nothing
+
+part = dp.read_indmftpr(os.path.join(HERE, 'CaOs2_partial.indmftpr'))
+assert part['sorts'][0]['included_ls'] == [2] # Ca d now an included shell
+assert part['sorts'][0]['correlated_ls'] == []
+assert part['sorts'][1]['correlated_ls'] == [2]
+
+nsym, ops = dp.read_dmftsym(os.path.join(HERE, 'CaOs2.dmftsym'))
+assert nsym == 16 and len(ops) == 16
+assert all(o['krotm'].shape == (3, 3) and len(o['perm']) == 3 for o in ops)
+
+with tempfile.TemporaryDirectory() as tmp:
+ path = os.path.join(tmp, 'noso.indmftpr')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n1\ncomplex\n0 2\n0 0\n0\n'
+ '-2.0 3.0 2 ! band window\n')
+ parsed = dp.read_indmftpr(path)
+ assert parsed['so'] == 0
+ assert (parsed['e_bot'], parsed['e_top'], parsed['proj_mode']) == (-2.0, 3.0, 2)
+
+ # Atomic sorts and symmetry operations are independent counts (P1 can have
+ # many inequivalent atoms but only the identity operation).
+ path = os.path.join(tmp, 'P1.dmftsym')
+ with open(path, 'w') as fh:
+ fh.write('1\n1 2\nSym. op. 1\n0 0 0 1\n'
+ '1 0 0\n0 1 0\n0 0 1\nGlobal->local\n'
+ '1\n1 0 0\n0 1 0\n0 0 1\n0 0 0 1\n'
+ '2\n1 0 0\n0 1 0\n0 0 1\n90 0 0 1\n')
+ _, _, refs = dp.read_dmftsym(path, rotloc=True)
+ assert len(refs) == 2
+ np.testing.assert_allclose(refs[1]['a'], np.pi / 2, atol=1e-14, rtol=0)
+
+ # A tiny but nonzero spin coupling must retain the full spinor basis.
+ theta = 1e-9
+ basis = np.eye(6)
+ basis[0, 0] = basis[3, 3] = np.cos(theta)
+ basis[0, 3], basis[3, 0] = np.sin(theta), -np.sin(theta)
+ path = os.path.join(tmp, 'basis.dat')
+ with open(path, 'w') as fh:
+ for row in basis:
+ fh.write(' ' + ' '.join(f'{value:.17E} 0' for value in row).replace('E', 'D') + '\n')
+ transform, mixing = dp.read_fromfile(path, 1)
+ assert mixing and transform.shape == (6, 6)
+ np.testing.assert_array_equal(transform, basis.conjugate())
+
+ # Each correlated l with an explicit irep owns its own correps record.
+ 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 = dp.read_indmftpr(path)
+ assert parsed['so'] == 1 and parsed['sorts'][0]['correlated_ls'] == [1, 2]
+
+ # Selecting individual irreps is not implemented: do not silently export
+ # the entire shell in place of the requested correlated subspace.
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\ncubic\n0 0 2\n0 0 2\n10\n1\n-2 3\n')
+ try:
+ dp.read_indmftpr(path)
+ except NotImplementedError as error:
+ assert 'irreducible' in str(error)
+ else:
+ raise AssertionError('split correlated shell was silently accepted')
+
+ # The input lacks the cross overlaps needed for two local orbitals;
+ # the native implementation rejects this case as well.
+ path = os.path.join(tmp, 'multi.almblm')
+ with open(path, 'w') as fh:
+ fh.write('1\n1\n2\n0\n1\n2\n')
+ try:
+ dp.read_almblm(path, dict(nsort=1, lmax=0, mult=[1]))
+ except NotImplementedError as error:
+ assert 'local orbital' in str(error)
+ else:
+ raise AssertionError('multiple local orbitals were silently accepted')
+
+print('wien2k_dmftproj_common: ok')
diff --git a/test/python/wien2k/wien2k_oubwin_python.py b/test/python/wien2k/wien2k_oubwin_python.py
index 26b427a..048f070 100644
--- a/test/python/wien2k/wien2k_oubwin_python.py
+++ b/test/python/wien2k/wien2k_oubwin_python.py
@@ -18,6 +18,7 @@
import numpy as np
from triqs_dftkit.wien2k import oubwin
+from triqs_dftkit.wien2k import _dmftproj as dp
HERE = os.path.dirname(os.path.abspath(__file__))
@@ -48,8 +49,8 @@ def _check(case):
(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.)
+ got = dp.select_window(nbmin, nbmin + len(energies) - 1,
+ np.array(energies), 0., 1.)
assert got == expected, (energies, got, expected)
# Three k-points with distinct weights, a Fermi shift of 10 Ry, and one
@@ -58,31 +59,17 @@ def _check(case):
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)
+ f.write('0\n3\n0\n10.0\n1.0\n0\n')
+ for ik, nbmin, energies, weights in (
+ (1, 10, [9., 10.25, 10.75, 12.], [.25, .80, .60, .25]),
+ (2, 20, [8., 9.], [.10, .10]),
+ (3, 30, [10.25, 10.75], [.75, .50])):
+ f.write(f'IK = {ik}\n0 0 0\n{ik} {nbmin} {nbmin + len(energies) - 1}\n')
+ for energy, weight in zip(energies, weights):
+ f.write(f'{weight} {energy}\n')
+ f.write('ATOM = 1\n1\n')
+ f.write('(0.,0.) (0.,0.)\n' * len(energies))
+ windows = oubwin._windows_from_spin(dp.read_almblm(almblm, dict(nsort=1, lmax=0, mult=[1])), dict(proj_mode=0, e_bot=0., e_top=1.), None)
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')
@@ -94,13 +81,4 @@ def _check(case):
_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_outband_python.py b/test/python/wien2k/wien2k_outband_python.py
new file mode 100644
index 0000000..7b70d51
--- /dev/null
+++ b/test/python/wien2k/wien2k_outband_python.py
@@ -0,0 +1,108 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+# Verify the pure-Python case.outband generator (triqs_dftkit.wien2k.outband)
+# reproduces the dmftproj Fortran band-structure projectors to machine
+# precision: regenerate case.outband from the band-mode almblm/struct/dmftsym,
+# the k-path (klist_band) and the band indmftpr, and compare to the native
+# dmftproj -band reference (full-precision cubic-d template).
+#
+# case.outband feeds convert_bands_input. It is the band analog of
+# case.ctqmcout: the same correlated-shell projector construction (raw Alm/Clm
+# projector -> rot_projectmat -> fromfile/cubic transmat -> Loewdin
+# orthonormalization) plus the raw Theta projectors of the parproj path, written
+# in the outband layout. The CaOs2_band fixture is SP+SO with the spin-mixing
+# fromfile (|j, m_j>) Os d shell (jbasis_d.dat); the wide window (-2.0, 3.0)
+# keeps 50 bands > 20 correlated spin-orbitals, so the Loewdin overlap is full
+# rank and the projectors match to machine precision.
+#
+# Band-mode specifics (dmftproj.f:557-581): nkband, the number of k-points along
+# the plotted path, is read from CaOs2_band.klist_band; the Fermi energy is the
+# last line of CaOs2_band.indmftpr (the band almblm keeps a placeholder eferm
+# record that the reader skips).
+#
+# CaOs2_band_cubic runs the same band inputs with a cubic (spin-diagonal) Os d
+# shell and an uncorrelated included Ca d shell, covering the non-mixing
+# projector path; its reference is native dmftproj -band output.
+
+import gzip
+import os
+import shutil
+import tempfile
+from unittest.mock import patch
+import numpy as np
+
+from triqs_dftkit.wien2k import outband
+
+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 a.lstrip('-').isdigit() and b.lstrip('-').isdigit():
+ assert int(a) == int(b), (a, b)
+ elif _isfloat(a) and _isfloat(b):
+ floats.append((float(a), float(b)))
+ else:
+ assert a == b, (a, b) # k-label tokens
+ a, b = np.array(floats).T
+ np.testing.assert_allclose(a, b, atol=1e-11, rtol=0, equal_nan=False)
+
+
+def _isfloat(tok):
+ try:
+ float(tok)
+ return True
+ except ValueError:
+ return False
+
+
+def _check(case, inputs=None):
+ inputs = inputs or case
+ tmp = tempfile.mkdtemp()
+ try:
+ shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
+ os.path.join(tmp, f'{case}.indmftpr'))
+ shutil.copy(os.path.join(HERE, f'{inputs}.klist_band'),
+ os.path.join(tmp, f'{case}.klist_band'))
+ shutil.copy(os.path.join(HERE, 'jbasis_d.dat'),
+ os.path.join(tmp, 'jbasis_d.dat'))
+ 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'{inputs}.{ext}.gz'),
+ os.path.join(tmp, f'{case}.{ext}'))
+ # Band parsing must also work without writable temporary storage.
+ with patch('tempfile.NamedTemporaryFile', side_effect=AssertionError('temporary file created')):
+ outband.write_outband(os.path.join(tmp, case))
+ with gzip.open(os.path.join(HERE, f'{case}.outband.gz'), 'rt') as fh:
+ _compare(os.path.join(tmp, f'{case}.outband'), fh.readlines())
+ with open(os.path.join(tmp, f'{case}.klist_band'), 'w') as fh:
+ fh.write('GAMMA 0 0 0 8 2.0\nX 4 0 0 8 2.0\nEND\n')
+ assert outband.read_k_list(os.path.join(tmp, f'{case}.klist_band')) == (2, [(1, 'GAMMA'), (2, 'X')])
+ try:
+ outband.write_outband(os.path.join(tmp, case))
+ except ValueError as exc:
+ assert 'k-point counts differ' in str(exc)
+ else:
+ raise AssertionError('inconsistent band-path dimensions must be rejected')
+ finally:
+ shutil.rmtree(tmp)
+
+
+_check('CaOs2_band')
+_check('CaOs2_band_cubic', inputs='CaOs2_band')
+
+print('wien2k_outband_python: ok')
diff --git a/test/python/wien2k/wien2k_parproj_python.py b/test/python/wien2k/wien2k_parproj_python.py
index 94a1ba0..3e42067 100644
--- a/test/python/wien2k/wien2k_parproj_python.py
+++ b/test/python/wien2k/wien2k_parproj_python.py
@@ -9,7 +9,8 @@
# an uncorrelated Ca d-shell. Inputs are the shared CaOs2 almblm/struct/dmftsym
# renamed to the CaOs2_partial case.
#
-# Floats compared at 1e-7 (dmftproj single-precision basis-transform floor).
+# Floats compared at 1e-11 against native dmftproj output (full-precision
+# cubic-d template).
import gzip
import os
@@ -38,14 +39,22 @@ def _compare(got_path, ref_lines):
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)
+ np.testing.assert_allclose(a, b, atol=1e-11, rtol=0, equal_nan=False)
-def _check(case):
+def _check(case, mode=0):
tmp = tempfile.mkdtemp()
try:
shutil.copy(os.path.join(HERE, f'{case}.indmftpr'),
os.path.join(tmp, f'{case}.indmftpr'))
+ if mode:
+ path = os.path.join(tmp, f'{case}.indmftpr')
+ with open(path) as fh:
+ lines = fh.read().splitlines()
+ lines[-1] = lines[-1] + ' 1' if mode == 1 else '60 76 2'
+ with open(path, 'w') as fh:
+ fh.write('\n'.join(lines) + '\n')
+ shutil.copy(os.path.join(HERE, 'jbasis_d.dat'), os.path.join(tmp, 'jbasis_d.dat'))
for ext in ('struct', 'dmftsym'):
shutil.copy(os.path.join(HERE, f'CaOs2.{ext}'),
os.path.join(tmp, f'{case}.{ext}'))
@@ -71,16 +80,18 @@ def _check(case):
np.testing.assert_allclose(actual, expected, atol=1e-14, rtol=0)
_check('CaOs2_partial')
-
-with tempfile.TemporaryDirectory() as tmp:
- path = os.path.join(tmp, 'P1.dmftsym')
- with open(path, 'w') as fh:
- fh.write('1\n1 2\nSym. op. 1\n0 0 0 1\n'
- '1 0 0\n0 1 0\n0 0 1\nGlobal->local\n'
- '1\n1 0 0\n0 1 0\n0 0 1\n0 0 0 1\n'
- '2\n1 0 0\n0 1 0\n0 0 1\n90 0 0 1\n')
- _, _, refs = parproj._read_dmftsym(path)
- assert len(refs) == 2
- np.testing.assert_allclose(refs[1]['a'], np.pi / 2, atol=1e-14, rtol=0)
+_check('CaOs2_jbasis_full')
+# All definitions select bands 60..76 here; the partial projectors and
+# below-window density must match the native energy-window reference.
+_check('CaOs2_partial', mode=1)
+_check('CaOs2_partial', mode=2)
print('wien2k_parproj_python: ok')
+
+# Equivalent atoms must share the average after identity and swap operations.
+raw_mixing = {1: np.diag([1., 3.]), 2: np.diag([5., 7.])}
+actual = parproj._symmetrize_densmat_mixing(
+ [dict(atom=1), dict(atom=2)], [dict(perm=[1, 2]), dict(perm=[2, 1])], 2,
+ [(np.eye(2), False)] * 2, raw_mixing)
+for matrix in actual.values():
+ np.testing.assert_allclose(matrix, np.diag([3., 5.]), atol=1e-14, rtol=0)
diff --git a/test/python/wien2k/wien2k_sympar_python.py b/test/python/wien2k/wien2k_sympar_python.py
index 0316fd5..d79866b 100644
--- a/test/python/wien2k/wien2k_sympar_python.py
+++ b/test/python/wien2k/wien2k_sympar_python.py
@@ -9,9 +9,8 @@
# 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.
+# Floats compared at 1e-11 against native dmftproj output (full-precision
+# cubic-d template).
import gzip
import os
@@ -41,12 +40,12 @@ def _compare(got_path, ref_lines):
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,
+ np.testing.assert_allclose(matrices, expected, atol=1e-11, 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)
+ atol=1e-11, rtol=0)
def _check(case):
@@ -70,19 +69,12 @@ def _check(case):
# 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))
+ 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.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)
+p = sympar.reptrans('cubic', 1)
+np.testing.assert_allclose(p @ p.conj().T, np.eye(3), atol=1e-11, 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')
diff --git a/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy b/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy
index b7e2e01..4cce82c 100644
Binary files a/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy and b/test/python/wien2k/wien2k_symqmc_jbasis.ref.npy differ
diff --git a/test/python/wien2k/wien2k_symqmc_l0.ref.npy b/test/python/wien2k/wien2k_symqmc_l0.ref.npy
index 679d0e0..ae517ee 100644
Binary files a/test/python/wien2k/wien2k_symqmc_l0.ref.npy and b/test/python/wien2k/wien2k_symqmc_l0.ref.npy differ
diff --git a/test/python/wien2k/wien2k_symqmc_mixing.py b/test/python/wien2k/wien2k_symqmc_mixing.py
index eb729a8..9702446 100644
--- a/test/python/wien2k/wien2k_symqmc_mixing.py
+++ b/test/python/wien2k/wien2k_symqmc_mixing.py
@@ -3,20 +3,20 @@
# (GPL-3.0-or-later)
################################################################################
-# Verify the pure-Python case.symqmc generator on the two paths the spin-diagonal
+# Verify the pure-Python case.symqmc generator on the paths the spin-diagonal
# cubic test (wien2k_symqmc_python) does not exercise:
#
# - a spin-mixing fromfile basis (the |j, m_j> basis), which uses the full
# 2(2l+1) spinor representation and the spinor time-reversal operator;
# - an l = 0 (s) correlated shell, the 2x2 spin phase block.
+# - cubic p harmonics, including the complex basis change of magnetic operations.
#
-# Both reuse the CaOs2 symmetry input (16 operations, 8 of them time-reversal).
-# The reference is the dmftproj Fortran case.symqmc for the same input, stored as
-# the matrix data in single precision: dmftproj reads the fromfile basis with a
-# single-precision CMPLX cast (set_ang_trans.f), so its output there carries a
-# ~1e-7 error. The generator is full double precision, so we compare at 1e-6 and
-# separately assert each Python spinor matrix is unitary.
+# These reuse the CaOs2 symmetry input (16 operations, 8 of them time-reversal).
+# The reference is the matrix data of the dmftproj case.symqmc for the same
+# input, stored in double precision and compared at 1e-11; each Python spinor
+# matrix is separately checked for unitarity.
+import gzip
import os
import shutil
import tempfile
@@ -36,7 +36,7 @@ def _matrix_data(symqmc_path, nsym, natom):
return np.array([float(x) for x in toks])
-def _run(indmftpr, extra=()):
+def _run(indmftpr, extra=(), ref_gz=None):
"""Generate case.symqmc in a scratch dir from the shared CaOs2 symmetry
input and the given indmftpr, returning (matrix data, shells, ops, so)."""
tmp = tempfile.mkdtemp()
@@ -51,6 +51,16 @@ def _run(indmftpr, extra=()):
shells, so = symqmc._read_correlated_shells(
case + '.indmftpr', case + '.struct')
natom = len(ops[0]['perm'])
+ if ref_gz is not None:
+ got = open(case + '.symqmc').read().split()
+ with gzip.open(os.path.join(HERE, ref_gz), 'rt') as fh:
+ ref = fh.read().split()
+ header = 2 + nsym * natom + nsym
+ assert len(got) == len(ref)
+ assert [int(x) for x in got[:header]] == [int(x) for x in ref[:header]]
+ np.testing.assert_allclose(np.array(got[header:], dtype=float),
+ np.array(ref[header:], dtype=float),
+ atol=1e-11, rtol=0, equal_nan=False)
return _matrix_data(case + '.symqmc', nsym, natom), shells, ops, so
finally:
shutil.rmtree(tmp)
@@ -65,7 +75,11 @@ def _check(indmftpr, ref_npy, extra=(), unit_tol=1e-7):
data, shells, ops, so = _run(indmftpr, extra)
ref = np.load(os.path.join(HERE, ref_npy))
assert data.shape == ref.shape, (data.shape, ref.shape)
- assert np.max(np.abs(data - ref)) < 1e-6, np.max(np.abs(data - ref))
+ assert np.max(np.abs(data - ref)) < 1e-11, np.max(np.abs(data - ref))
+ _check_unitary(shells, ops, so, unit_tol)
+
+
+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))
@@ -73,7 +87,7 @@ def _check(indmftpr, ref_npy, extra=(), unit_tol=1e-7):
assert dev < unit_tol, ('not unitary', dev)
-# spin-mixing |j, m_j> basis: the path dft_tools #148 singles out
+# spin-mixing |j, m_j> basis
data, shells, _, _ = _run('CaOs2_jbasis.indmftpr',
(('jbasis_d.dat', 'jbasis_d.dat'),))
assert shells[0]['mixing'] and shells[0]['P'].shape == (10, 10)
@@ -85,4 +99,9 @@ def _check(indmftpr, ref_npy, extra=(), unit_tol=1e-7):
assert shells[0]['l'] == 0
_check('CaOs2_l0.indmftpr', 'wien2k_symqmc_l0.ref.npy', unit_tol=1e-12)
+# Cubic p basis: full native header, all operations and both correlated atoms.
+_, shells, ops, so = _run('CaOs2_pcubic.indmftpr', ref_gz='CaOs2_pcubic.symqmc.gz')
+assert all(sh['l'] == 1 and not sh['mixing'] for sh in shells)
+_check_unitary(shells, ops, so, unit_tol=1e-12)
+
print('wien2k_symqmc_mixing: ok')