[wien2k] fix dmftproj port for non-centrosymmetric cells

python/triqs_dftkit/wien2k/_dmftproj.py

@@ -103,7 +103,9 @@ def read_two_complex(r):
 _COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)}
 
 _CUBIC = {
-    1: np.array([[2**-0.5, 0, -2**-0.5], [-1j*2**-0.5, 0, -1j*2**-0.5], [0, 1, 0]], dtype=complex),
+    1: np.array([[2 ** -0.5, 0, -(2 ** -0.5)],
+                 [-1j * 2 ** -0.5, 0, -1j * 2 ** -0.5],
+                 [0, 1, 0]], dtype=complex),
     2: np.array([
         [0, 0, 1, 0, 0],
         [2 ** -0.5, 0, 0, 0, 2 ** -0.5],
@@ -225,6 +227,58 @@ def mixing_timeinv_op(l, P):
     return P @ tinv @ P.T
 
 
+def euler(R):
+    """Euler angles (a, b, c) of the proper rotation R, in the WIEN2k
+    convention of set_rotloc.f/SRC_lapwdm euler.f: R is the transpose of a
+    krotm matrix. a, c in [0, 2pi), b in [0, pi]."""
+    z = np.array([0.0, 0.0, 1.0])
+    y = np.array([0.0, 1.0, 0.0])
+    yyy = R[:, 1]
+    zz = R[:, 2]
+    yy = np.cross(z, zz)
+    y_norm = math.sqrt(yy @ yy)
+    if y_norm < 1e-10:
+        a = math.acos(max(-1.0, min(1.0, y @ yyy)))
+        if z @ zz > 0.0:
+            b = c = 0.0
+            if yyy[0] > 0.0:
+                a = 2 * math.pi - a
+        else:
+            c, a, b = a, 0.0, math.pi
+            if yyy[0] < 0.0:
+                c = 2 * math.pi - c
+        return a, b, c
+    yy = yy / y_norm
+    a = math.acos(max(-1.0, min(1.0, y @ yy)))
+    b = math.acos(max(-1.0, min(1.0, z @ zz)))
+    c = math.acos(max(-1.0, min(1.0, yy @ yyy)))
+    if yy[0] > 0.0:
+        a = 2 * math.pi - a
+    if np.cross(yy, yyy) @ zz < 0.0:
+        c = 2 * math.pi - c
+    return a, b, c
+
+
+def _rloc_spinor_s(ref, op, timeinv):
+    """2x2 spinor Rloc of an s shell under SP+SO (outputqmc.f, l=0): built
+    from the Euler angles of the composed proper rotation
+    krotm(op) krotm(ref), which fixes its SU(2) sign. For a magnetic op the
+    columns carry the time reversal -i sigma_y."""
+    a, b, g = euler((op['krotm'] @ ref['krotm']).T)
+    s = np.zeros((2, 2), dtype=complex)
+    if timeinv:
+        s[1, 0] = np.exp(1j * (a + g) / 2) * math.cos(b / 2)
+        s[0, 1] = -np.conj(s[1, 0])
+        s[1, 1] = -np.exp(-1j * (a - g) / 2) * math.sin(b / 2)
+        s[0, 0] = np.conj(s[1, 1])
+    else:
+        s[0, 0] = np.exp(1j * (a + g) / 2) * math.cos(b / 2)
+        s[1, 1] = np.conj(s[0, 0])
+        s[0, 1] = np.exp(-1j * (a - g) / 2) * math.sin(b / 2)
+        s[1, 0] = -np.conj(s[0, 1])
+    return s
+
+
 def rotloc_rotl_so(l, ref, ops, iatom, iref):
     """The composed 2(2l+1) Rloc rotation rotloc(iatom)%rotl(l) under SP+SO
     (setsym.f:496-528 + set_rotloc.f): the representative spinor rotloc
@@ -233,6 +287,12 @@ def rotloc_rotl_so(l, ref, ops, iatom, iref):
     Returns (rotl, timeinv); callers apply their own basis transform (the full
     transmat for a mixing basis, blkdiag(transmat, transmat) otherwise)."""
     d = 2 * l + 1
+    op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+    det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
+            - op['krotm'][0, 1] * op['krotm'][1, 0])
+    timeinv = det2 < 0.0
+    if l == 0:
+        return _rloc_spinor_s(ref, op, timeinv), timeinv
     Dref = dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
     f = (ref['a'] + ref['g']) / 2.0
     spmt = np.zeros((2, 2), dtype=complex)
@@ -247,10 +307,6 @@ def rotloc_rotl_so(l, ref, ops, iatom, iref):
     rotl[:d, d:] = spmt[0, 1] * Dref
     rotl[d:, :d] = spmt[1, 0] * Dref
 
-    op = next(o for o in ops if o['perm'][iref - 1] == iatom)
-    det2 = (op['krotm'][0, 0] * op['krotm'][1, 1]
-            - op['krotm'][0, 1] * op['krotm'][1, 0])
-    timeinv = det2 < 0.0
     srot_phase = (op['g'] - op['a']) if timeinv else (op['a'] + op['g'])
     rotl_sym = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
     if timeinv:
@@ -263,6 +319,13 @@ def rotloc_rotl_so(l, ref, ops, iatom, iref):
     return rotl, timeinv
 
 
+def rotloc_orbital(l, ref, ops, iatom, iref):
+    """Orbital part of the full local rotation R[isym] * rotloc_ref."""
+    op = next(o for o in ops if o['perm'][iref - 1] == iatom)
+    return dmat(l, op['a'], op['b'], op['g'], float(op['iprop'])) @ \
+        dmat(l, ref['a'], ref['b'], ref['g'], float(ref['iprop']))
+
+
 def mixing_rotrep(op, l, P, ti):
     """The full 2(2l+1) spinor representation D(R)_{new_i} = P spinrot P^dag of
     one symmetry operation in a spin-coupling (mixing) basis, with the spinor
@@ -625,18 +688,6 @@ def fmt_scalar(x):
     return '%.16f' % x if abs(x) < 1e5 else '%.16E' % x
 
 
-def require_identity_struct(case):
-    """Fail explicitly until representative struct rotations are ported."""
-    with open(case + '.struct') as f:
-        lines = f.readlines()
-    for i, line in enumerate(lines):
-        if 'LOCAL ROT MATRIX:' in line:
-            rows = [line.split(':', 1)[1], *lines[i + 1:i + 3]]
-            rotation = np.array([[to_float(x) for x in row.split()] for row in rows])
-            if not np.allclose(rotation, np.eye(3), atol=1e-7, rtol=0):
-                raise NotImplementedError('non-identity struct local rotations are not implemented yet')
-
-
 def validate_spin_windows(spins, info, windows, below=None):
     """Require valid selected windows and the same absolute bands in both spins."""
     if len(spins) != 2:

python/triqs_dftkit/wien2k/ctqmcout.py

@@ -60,9 +60,9 @@
 from scipy.linalg.lapack import zheev
 from scipy.linalg.blas import zgemm
 
-from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_spin_windows, band_index_window, dmat, fmt_scalar, read_almblm,
+from ._dmftproj import (build_orbs as _build_orbs, validate_spin_windows, band_index_window, dmat, fmt_scalar, read_almblm,
                         read_dmftsym, read_fromfile, read_indmftpr, reptrans,
-                        rotloc_rotl_so, select_band_window, select_window,
+                        rotloc_orbital, rotloc_rotl_so, select_band_window, select_window,
                         tmat, write_row)
 
 
@@ -102,34 +102,23 @@ def _build_crorbs(info):
 
 # --- Rloc spinor representation (set_rotloc.f / setsym.f) ---------------------
 
-def _rloc_rotrep(op, l, transmat, ifSO):
-    """rotloc%rotrep(l)%mat, the Rloc rotation in the new basis for a non-mixing
-    shell, with identity struct local rotation (rotloc_ref Euler = 0).
-
-    Under SO it is the 2(2l+1) spinor rotation: rotloc%rotl =
-    blkdiag(ephase*D, conj(ephase)*D), rotrep = S rotl S^H, S =
-    blkdiag(transmat, transmat). Without SO srot%timeinv is always false, there
-    is no spin phase, and it reduces to the bare (2l+1) transmat D transmat^H
+def _rloc_rotrep(l, ref, ops, iatom, iref, transmat, ifSO):
+    """rotloc(iatom)%rotrep(l)%mat for a non-mixing shell: the composed Rloc
+    spinor rotation rotloc_ref then the symmetry op mapping iref onto iatom
+    (set_rotloc.f / setsym.f), put into the new basis with S = blkdiag(transmat,
+    transmat). Shares rotloc_rotl_so with the mixing path so the non-trivial
+    struct reference rotation rotloc_ref (e.g. Te, beta=pi/2) is not dropped.
+    Without SO it reduces to the bare (2l+1) D(op)D(ref) transmat sandwich
     (set_rotloc.f non-SO branch). Returns (rotrep, timeinv)."""
-    a, b, g, iprop = op['a'], op['b'], op['g'], op['iprop']
-    D = dmat(l, a, b, g, float(iprop))
     if not ifSO:
+        D = rotloc_orbital(l, ref, ops, iatom, iref)
         return transmat @ D @ np.conj(transmat.T), False
-    krotm = op['krotm']
-    det2 = krotm[0, 0] * krotm[1, 1] - krotm[0, 1] * krotm[1, 0]
-    timeinv = det2 < 0.0
-    phase = (g - a) if timeinv else (a + g)
-    if timeinv:
-        D = tmat(l) @ np.conj(D)       # setsym.f:320-326 orbital time reversal
-    ephase = np.exp(1j * phase / 2)
+    rotl, timeinv = rotloc_rotl_so(l, ref, ops, iatom, iref)
     d = 2 * l + 1
-    rotl = np.zeros((2 * d, 2 * d), dtype=complex)
-    rotl[:d, :d] = ephase * D
-    rotl[d:, d:] = np.conj(ephase) * D
     S = np.zeros((2 * d, 2 * d), dtype=complex)
     S[:d, :d] = transmat
     S[d:, d:] = transmat
-    rotrep = S @ rotl @ np.conj(S.T)
+    rotrep = (S @ rotl @ S.T) if timeinv else (S @ rotl @ np.conj(S.T))
     return rotrep, timeinv
 
 
@@ -148,7 +137,6 @@ def write_ctqmcout(case):
     """Read <case>.almblm{up,dn}, <case>.indmftpr, <case>.struct,
     <case>.dmftsym and write <case>.ctqmcout in the dmftproj format."""
     info = read_indmftpr(case + '.indmftpr')
-    require_identity_struct(case)
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
 
     up, dn = case + '.almblmup', case + '.almblmdn'
@@ -230,12 +218,6 @@ def write_ctqmcout(case):
     # rot_projectmat applies dmat(l, a, b, g, iprop) in the |lm> basis
     # (rot_projectmat.f:59-65). For atom 2 this is a non-trivial in-plane
     # rotation (a=270) that mixes the m-rows; omitting it flips the m=1,2 rows.
-    rotloc_op = {}
-    for icr, cr in enumerate(crorbs):
-        iref = sum(info['mult'][:cr['sort'] - 1]) + 1
-        rotloc_op[icr] = next(o for o in ops
-                              if o['perm'][iref - 1] == cr['atom'])
-
     # ---- raw correlated projector mat_rep, then Loewdin orthonormalize ----
     # Non-mixing: mat_rep[(icr, ik, is)] -> (2l+1, nbsel) per spin block.
     # Mixing: mat_rep[(icr, ik)] -> (2*(2l+1), nbsel), the stacked spinor block.
@@ -247,8 +229,8 @@ def write_ctqmcout(case):
         d = 2 * l + 1
         transmat = transmats[(l, sort)]
         ismix = mixing[(l, sort)]
-        op = rotloc_op[icr]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[sort - 1]
+        rot = rotloc_orbital(l, ref, ops, atom, sum(info['mult'][:sort - 1]) + 1)
         for ik in range(nk):
             incl, nb_bot, nb_top = windows[ik]
             if not incl:
@@ -328,8 +310,8 @@ def write_ctqmcout(case):
             rloc_blocks.append(_rloc_rotrep_mixing(
                 l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat))
         else:
-            op = next(o for o in ops if o['perm'][iref - 1] == cr['atom'])
-            rloc_blocks.append(_rloc_rotrep(op, l, transmat, ifSO))
+            rloc_blocks.append(_rloc_rotrep(
+                l, rotloc_ref[cr['sort'] - 1], ops, cr['atom'], iref, transmat, ifSO))
 
     # ---- write ----
     with open(case + '.ctqmcout', 'w') as f:

python/triqs_dftkit/wien2k/outband.py

@@ -51,8 +51,8 @@
 import numpy as np
 from scipy.linalg.blas import zgemm
 
-from ._dmftproj import (require_identity_struct, validate_spin_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
-                        read_indmftpr, reptrans, select_window, to_float,
+from ._dmftproj import (validate_spin_windows, dmat, read_almblm, read_dmftsym, read_fromfile,
+                        read_indmftpr, reptrans, rotloc_orbital, select_window, to_float,
                         write_row)
 from .ctqmcout import _build_crorbs, _sqrt_inv
 from .parproj import (_build_matn_rep, _build_matn_rep_mixing, _build_orbs)
@@ -108,7 +108,6 @@ def write_outband(case):
     info, eferm = _read_indmftpr_band(case + '.indmftpr')
     if not info['so']:
         raise NotImplementedError('non-spin-orbit projectors are not implemented yet')
-    require_identity_struct(case)
     if info['proj_mode'] != 0:
         raise NotImplementedError('this output currently supports projection mode 0')
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
@@ -158,12 +157,6 @@ def write_outband(case):
 
     # rot_projectmat local rotation: the op mapping the sort representative onto
     # this atom (set_projections.f via rot_projectmat). Same as ctqmcout/parproj.
-    rotloc_op = {}
-    for o in orbs:
-        iref = sum(info['mult'][:o['sort'] - 1]) + 1
-        rotloc_op[o['atom']] = next(op for op in ops
-                                    if op['perm'][iref - 1] == o['atom'])
-
     # ---- correlated projector mat_rep, then Loewdin orthonormalize ----
     mat_rep = {}
     for icr, cr in enumerate(crorbs):
@@ -173,8 +166,8 @@ def write_outband(case):
         d = 2 * l + 1
         transmat = transmats[(l, sort)]
         ismix = mixing[(l, sort)]
-        op = rotloc_op[atom]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[sort - 1]
+        rot = rotloc_orbital(l, ref, ops, atom, sum(info['mult'][:sort - 1]) + 1)
         for ik in range(nk):
             incl, nb_bot, nb_top = windows[ik]
             if not incl:
@@ -236,8 +229,9 @@ def write_outband(case):
     for o in orbs:
         l = o['l']
         transmat = transmats[(l, o['sort'])]
-        op = rotloc_op[o['atom']]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[o['sort'] - 1]
+        iref = sum(info['mult'][:o['sort'] - 1]) + 1
+        rot = rotloc_orbital(l, ref, ops, o['atom'], iref)
         if mixing[(l, o['sort'])]:
             matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
                 o, info, spins, windows, transmat, rot)

python/triqs_dftkit/wien2k/parproj.py

@@ -57,9 +57,9 @@
 import os
 import numpy as np
 
-from ._dmftproj import (build_orbs as _build_orbs, require_identity_struct, validate_spin_windows, band_index_window, dmat, mixing_rotrep, read_almblm, read_dmftsym,
-                        read_fromfile, read_indmftpr, reptrans,
-                        rotloc_rotl_so, select_band_window, select_window, tmat, write_row)
+from ._dmftproj import (build_orbs as _build_orbs, validate_spin_windows, band_index_window, dmat, mixing_rotrep, paramag_time_op,
+                        read_almblm, read_dmftsym, read_fromfile, read_indmftpr, reptrans,
+                        rotloc_orbital, rotloc_rotl_so, select_band_window, select_window, tmat, write_row)
 
 
 def _op_timeinv(op, ifSO=True):
@@ -117,15 +117,14 @@ def _srot_rotrep_nonmixing(op, l, transmat, ifSO=True):
 
 # --- rotloc rotrep (set_rotloc.f) under SP+SO, non-mixing --------------------
 
-def _rotloc_rotrep_nonso(orb, ops, info, transmat):
+def _rotloc_rotrep_nonso(orb, ops, info, ref, transmat):
     """rotloc(iatom)%rotrep(l)%mat without SO: the bare (2l+1) Rloc rotation in
     the new basis, transmat D(Rloc)_{lm} transmat^H with D(Rloc) the symmetry op
     mapping the sort representative onto this atom (set_rotloc.f non-SO branch).
     srot%timeinv is always false under non-SO, so timeinv=False."""
     l = orb['l']
     iref = sum(info['mult'][:orb['sort'] - 1]) + 1
-    op = next(o for o in ops if o['perm'][iref - 1] == orb['atom'])
-    D = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+    D = rotloc_orbital(l, ref, ops, orb['atom'], iref)
     return transmat @ D @ np.conj(transmat.T), False
 
 
@@ -389,6 +388,14 @@ def _symmetrize_densmat_nonso(orbs, info, ops, nsym, srot_rotreps, dens_raw):
     return out
 
 
+def _add_timeinv_nonso(blocks, l, transmat):
+    """add_timeinv for a paramagnetic non-SO shell (timeinv.f, nsp=1): the
+    single spin block D becomes (D + U conj(D) U^H)/2, U the time-reversal
+    operator in the shell basis."""
+    U = paramag_time_op(dict(l=l, P=transmat))
+    return [(D + U @ np.conj(D) @ np.conj(U.T)) / 2 for D in blocks]
+
+
 def _rotdens_densmat_nonso(blocks, rotrep_loc):
     """rotdens_mat for the non-SO non-mixing path (rot_dens.f:197-228):
     each (2l+1) spin block is rotated as conj(Rloc^T) D Rloc with the bare
@@ -422,7 +429,6 @@ def write_parproj(case):
     """Read <case>.almblm{up,dn}, <case>.indmftpr, <case>.struct,
     <case>.dmftsym and write <case>.parproj in the dmftproj format."""
     info = read_indmftpr(case + '.indmftpr')
-    require_identity_struct(case)
     nsym, ops, rotloc_ref = read_dmftsym(case + '.dmftsym', rotloc=True)
 
     spin_files = [case + '.almblmup', case + '.almblmdn']
@@ -476,21 +482,16 @@ def write_parproj(case):
 
     # rot_projectmat local rotation: the op mapping the representative atom of
     # the sort onto this atom (set_projections.f via rot_projectmat).
-    rotloc_op = {}
-    for o in orbs:
-        iref = sum(info['mult'][:o['sort'] - 1]) + 1
-        rotloc_op[o['atom']] = next(op for op in ops
-                                    if op['perm'][iref - 1] == o['atom'])
-
     # ---- projectors matn_rep per orbital (energy window for block A, full
     # band range for the density matrix) ----
     matn_reps = {}
     matn_reps_full = {}
     for o in orbs:
         l = o['l']
         transmat = transmats[(l, o['sort'])]
-        op = rotloc_op[o['atom']]
-        rot = dmat(l, op['a'], op['b'], op['g'], float(op['iprop']))
+        ref = rotloc_ref[o['sort'] - 1]
+        iref = sum(info['mult'][:o['sort'] - 1]) + 1
+        rot = rotloc_orbital(l, ref, ops, o['atom'], iref)
         if mixing[(l, o['sort'])]:
             matn_reps[(o['atom'], o['l'])] = _build_matn_rep_mixing(
                 o, info, spins, windows, transmat, rot)
@@ -506,11 +507,11 @@ def write_parproj(case):
     rotloc_rotrep = {}
     for o in orbs:
         transmat = transmats[(o['l'], o['sort'])]
+        ref = rotloc_ref[o['sort'] - 1]
         if not ifSO:
             rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_nonso(
-                o, ops, info, transmat)
+                o, ops, info, ref, transmat)
         else:
-            ref = rotloc_ref[o['sort'] - 1]
             rotloc_rotrep[(o['atom'], o['l'])] = _rotloc_rotrep_so(
                 o, ops, info, ref, transmat, mixing[(o['l'], o['sort'])])
 
@@ -544,9 +545,12 @@ def write_parproj(case):
                                                  srot_rotreps, dens_raw)
             for o in nonmix_orbs:
                 key = (o['atom'], o['l'])
+                blocks = dens_sym[key]
+                if not ifSP:
+                    blocks = _add_timeinv_nonso(
+                        blocks, o['l'], transmats[(o['l'], o['sort'])])
                 rotrep_loc, _ = rotloc_rotrep[key]
-                densprint[key] = _rotdens_densmat_nonso(
-                    dens_sym[key], rotrep_loc)
+                densprint[key] = _rotdens_densmat_nonso(blocks, rotrep_loc)
 
     mixing_groups = {}
     for o in orbs: