|
| 1 | +################################################################################ |
| 2 | +# |
| 3 | +# TRIQS: a Toolbox for Research in Interacting Quantum Systems |
| 4 | +# |
| 5 | +# Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola |
| 6 | +# |
| 7 | +# TRIQS is free software: you can redistribute it and/or modify it under the |
| 8 | +# terms of the GNU General Public License as published by the Free Software |
| 9 | +# Foundation, either version 3 of the License, or (at your option) any later |
| 10 | +# version. |
| 11 | +# |
| 12 | +# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY |
| 13 | +# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
| 14 | +# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more |
| 15 | +# details. |
| 16 | +# |
| 17 | +# You should have received a copy of the GNU General Public License along with |
| 18 | +# TRIQS. If not, see <http://www.gnu.org/licenses/>. |
| 19 | +# |
| 20 | +################################################################################ |
| 21 | + |
| 22 | +"""Pure-Python generation of the dmftproj correlated-shell symmetry file |
| 23 | +(case.symqmc), replacing the corresponding output of the dmftproj Fortran |
| 24 | +executable. |
| 25 | +
|
| 26 | +Given the dmftproj symmetry input (case.dmftsym), the projector definition |
| 27 | +(case.indmftpr) and the structure (case.struct), this builds the spinor symmetry |
| 28 | +matrices and writes them in the case.symqmc format the converter reads. It |
| 29 | +reproduces the Fortran construction: the Wigner D matrix (dmat), the orbital |
| 30 | +time-reversal operator for the magnetic (SP+SO) operations, the basis transform |
| 31 | +to the chosen angular harmonics, and the spin-1/2 phase blocks. |
| 32 | +
|
| 33 | +Currently covers the non-mixing spin-diagonal bases (complex, cubic) with |
| 34 | +spin-orbit; this is the path used by the Wien2k SOC + spin-polarized workflow. |
| 35 | +""" |
| 36 | + |
| 37 | +import math |
| 38 | +import numpy as np |
| 39 | + |
| 40 | +# --- angular bases (transpose(P) = <new|m>, m = -l..l) ----------------------- |
| 41 | + |
| 42 | +_COMPLEX = {l: np.eye(2 * l + 1, dtype=complex) for l in range(4)} |
| 43 | + |
| 44 | +# standard cubic harmonics, Wien2k convention |
| 45 | +_CUBIC = { |
| 46 | + 1: np.array([[0, 1, 0], [-1j, 0, -1j], [1, 0, -1]], dtype=complex) / 1.0, |
| 47 | + 2: np.array([ |
| 48 | + [0, 0, 1, 0, 0], |
| 49 | + [2 ** -0.5, 0, 0, 0, 2 ** -0.5], |
| 50 | + [-(2 ** -0.5), 0, 0, 0, 2 ** -0.5], |
| 51 | + [0, 2 ** -0.5, 0, -(2 ** -0.5), 0], |
| 52 | + [0, 2 ** -0.5, 0, 2 ** -0.5, 0], |
| 53 | + ], dtype=complex), |
| 54 | +} |
| 55 | + |
| 56 | + |
| 57 | +def _reptrans(basis, l): |
| 58 | + if basis == 'cubic' and l in _CUBIC: |
| 59 | + return _CUBIC[l] |
| 60 | + return _COMPLEX[l] |
| 61 | + |
| 62 | + |
| 63 | +# --- Wigner D matrix (dmftproj convention, setsym.f) ------------------------- |
| 64 | + |
| 65 | +def _small_d(l, m, n, b): |
| 66 | + f1 = (math.factorial(l + m) * math.factorial(l - m)) / \ |
| 67 | + (math.factorial(l + n) * math.factorial(l - n)) |
| 68 | + s = 0.0 |
| 69 | + for t in range(0, 2 * l + 1): |
| 70 | + if (l - m - t) >= 0 and (l - n - t) >= 0 and (t + n + m) >= 0: |
| 71 | + f2 = (math.factorial(l + n) * math.factorial(l - n)) / \ |
| 72 | + (math.factorial(l - m - t) * math.factorial(m + n + t) * |
| 73 | + math.factorial(l - n - t) * math.factorial(t)) |
| 74 | + f3 = 1.0 if (2 * l - m - n - 2 * t) == 0 else math.sin(b / 2) ** (2 * l - m - n - 2 * t) |
| 75 | + f4 = 1.0 if (2 * t + n + m) == 0 else math.cos(b / 2) ** (2 * t + n + m) |
| 76 | + s += (-1) ** (l - m - t) * f2 * f3 * f4 |
| 77 | + return math.sqrt(f1) * s |
| 78 | + |
| 79 | + |
| 80 | +def _dmat(l, a, b, c, det): |
| 81 | + D = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex) |
| 82 | + for m in range(-l, l + 1): |
| 83 | + for n in range(-l, l + 1): |
| 84 | + v = np.exp(1j * n * a) * np.exp(1j * m * c) * _small_d(l, m, n, b) |
| 85 | + if det < -0.5: |
| 86 | + v *= (-1) ** l |
| 87 | + D[m + l, n + l] = v |
| 88 | + return D |
| 89 | + |
| 90 | + |
| 91 | +def _timeinv_orbital(l, mat): |
| 92 | + T = np.zeros((2 * l + 1, 2 * l + 1), dtype=complex) |
| 93 | + for m in range(-l, l + 1): |
| 94 | + T[-m + l, m + l] = (-1) ** m |
| 95 | + return T @ np.conj(mat) |
| 96 | + |
| 97 | + |
| 98 | +# --- input parsing ----------------------------------------------------------- |
| 99 | + |
| 100 | +def _read_dmftsym(path): |
| 101 | + lines = open(path).read().split('\n') |
| 102 | + nsym = int(lines[0].split()[0]) |
| 103 | + perms = [[int(x) for x in lines[1 + i].split()] for i in range(nsym)] |
| 104 | + rest = lines[1 + nsym:] |
| 105 | + starts = [i for i, l in enumerate(rest) if 'Sym. op.' in l] |
| 106 | + ops = [] |
| 107 | + for k, s in enumerate(starts): |
| 108 | + ang = rest[s + 1].split() |
| 109 | + a, b, c = (math.radians(float(x)) for x in ang[:3]) |
| 110 | + krotm = np.array([[float(x) for x in rest[s + 2 + r].split()] for r in range(3)]) |
| 111 | + ops.append(dict(perm=perms[k], a=a, b=b, c=c, krotm=krotm)) |
| 112 | + return nsym, ops |
| 113 | + |
| 114 | + |
| 115 | +def _read_correlated_shells(indmftpr, struct): |
| 116 | + """Return the list of correlated shells (l, basis), one entry per correlated |
| 117 | + atom, plus the SO flag, from case.indmftpr and case.struct multiplicities.""" |
| 118 | + raw = [l.split('!')[0].strip() for l in open(indmftpr)] |
| 119 | + raw = [l for l in raw if l != ''] |
| 120 | + nsort = int(raw[0].split()[0]) |
| 121 | + mult = [int(x) for x in raw[1].split()][:nsort] |
| 122 | + i = 3 |
| 123 | + so = 0 |
| 124 | + shells = [] |
| 125 | + for isort in range(nsort): |
| 126 | + basis = raw[i] |
| 127 | + i += 1 |
| 128 | + l_inc = [int(x) for x in raw[i].split()] |
| 129 | + i += 1 |
| 130 | + ireps = [int(x) for x in raw[i].split()] |
| 131 | + i += 1 |
| 132 | + correlated_ls = [l for l in range(len(l_inc)) if l_inc[l] == 2] |
| 133 | + if any(n > 0 for n in ireps): |
| 134 | + i += 1 # skip the correps line |
| 135 | + if correlated_ls: |
| 136 | + so = int(raw[i].split()[0]) # SO flag follows a correlated sort |
| 137 | + i += 1 |
| 138 | + for l in correlated_ls: |
| 139 | + for _ in range(mult[isort]): |
| 140 | + shells.append(dict(l=l, basis=basis)) |
| 141 | + return shells, so |
| 142 | + |
| 143 | + |
| 144 | +def write_symqmc(case): |
| 145 | + """Write <case>.symqmc from <case>.dmftsym, <case>.indmftpr, <case>.struct.""" |
| 146 | + nsym, ops = _read_dmftsym(case + '.dmftsym') |
| 147 | + shells, so = _read_correlated_shells(case + '.indmftpr', case + '.struct') |
| 148 | + natom = len(ops[0]['perm']) |
| 149 | + |
| 150 | + timeinv = [] |
| 151 | + for op in ops: |
| 152 | + det2 = op['krotm'][0, 0] * op['krotm'][1, 1] - op['krotm'][0, 1] * op['krotm'][1, 0] |
| 153 | + timeinv.append(1 if (so and det2 < 0.0) else 0) |
| 154 | + |
| 155 | + with open(case + '.symqmc', 'w') as f: |
| 156 | + f.write('%6d %6d\n' % (nsym, natom)) |
| 157 | + for op in ops: |
| 158 | + f.write(''.join('%6d ' % p for p in op['perm']) + '\n') |
| 159 | + if so: |
| 160 | + f.write(''.join('%6d ' % t for t in timeinv) + '\n') |
| 161 | + for isym, op in enumerate(ops): |
| 162 | + for sh in shells: |
| 163 | + f.write(_format_matrix(_shell_matrix(op, sh, timeinv[isym]))) |
| 164 | + |
| 165 | + |
| 166 | +def _shell_matrix(op, shell, ti): |
| 167 | + l, basis = shell['l'], shell['basis'] |
| 168 | + a, b, c = op['a'], op['b'], op['c'] |
| 169 | + det = np.linalg.det(op['krotm']) |
| 170 | + rotl = _dmat(l, a, b, c, det) |
| 171 | + if ti: |
| 172 | + rotl = _timeinv_orbital(l, rotl) |
| 173 | + P = _reptrans(basis, l) |
| 174 | + rotrep = P @ rotl @ np.conj(P.T) |
| 175 | + phase = (c - a) if ti else (a + c) |
| 176 | + e = np.exp(1j * phase / 2) |
| 177 | + d = 2 * l + 1 |
| 178 | + mat = np.zeros((2 * d, 2 * d), dtype=complex) |
| 179 | + mat[:d, :d] = e * rotrep |
| 180 | + mat[d:, d:] = np.conj(e) * rotrep |
| 181 | + return mat |
| 182 | + |
| 183 | + |
| 184 | +def _format_matrix(mat): |
| 185 | + out = [] |
| 186 | + for part in (mat.real, mat.imag): |
| 187 | + for row in part: |
| 188 | + out.append(''.join(' %.14E' % x for x in row) + '\n') |
| 189 | + return ''.join(out) |
0 commit comments