[wien2k] rotate the symmetry matrices into the local solver frame (dft_tools#148)

python/triqs_dftkit/wien2k/_symmetry.py

@@ -0,0 +1,37 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+"""Small, NumPy-only helpers for Wien2k symmetry conversion."""
+
+import numpy
+
+
+def rotate_to_local_frames(mat, perm, time_inv, orbits, rot_mat):
+    """Rotate global-frame symmetry matrices into each shell's solver frame.
+
+    A symmetry maps source shell ``i`` to ``j``.  In local frames its matrix is
+    ``R_j^dag M_ji R_i``.  For an antiunitary operation Symmetry.symmetrize
+    conjugates the source quantity, so the source rotation is conjugated too.
+    """
+    lookup = {}
+    for iorb, orbit in enumerate(orbits):
+        key = (orbit['atom'], orbit['sort'], orbit['l'], orbit['dim'])
+        if key in lookup:
+            raise ValueError('duplicate shell identity while mapping symmetry frames')
+        lookup[key] = iorb
+
+    for isym, atom_perm in enumerate(perm):
+        targets = []
+        for orbit in orbits:
+            key = (atom_perm[orbit['atom'] - 1], orbit['sort'],
+                   orbit['l'], orbit['dim'])
+            if key not in lookup:
+                raise ValueError('symmetry operation maps to an unknown shell')
+            targets.append(lookup[key])
+        for iorb, jorb in enumerate(targets):
+            source = rot_mat[iorb].conjugate() if time_inv[isym] else rot_mat[iorb]
+            mat[isym][iorb] = (rot_mat[jorb].conjugate().transpose()
+                               @ mat[isym][iorb] @ source)
+    return mat

python/triqs_dftkit/wien2k/converter.py

@@ -25,6 +25,8 @@
 
 from types import *
 import numpy
+
+from ._symmetry import rotate_to_local_frames
 from h5 import *
 from ..converter_tools import *
 import os.path
@@ -364,7 +366,8 @@ def convert_dft_input(self):
         # Symmetries are used, so now convert symmetry information for
         # *correlated* orbitals:
         self.convert_symmetry_input(orbits=self.corr_shells, symm_file=self.symmcorr_file,
-                                    symm_subgrp=self.symmcorr_subgrp, SO=self.SO, SP=self.SP)
+                                    symm_subgrp=self.symmcorr_subgrp, SO=self.SO, SP=self.SP,
+                                    rot_mat=rot_mat)
         self.convert_misc_input()
 
     def convert_parproj_input(self):
@@ -470,7 +473,8 @@ def convert_parproj_input(self):
         # Symmetries are used, so now convert symmetry information for *all*
         # orbitals:
         self.convert_symmetry_input(orbits=self.shells, symm_file=self.symmpar_file,
-                                    symm_subgrp=self.symmpar_subgrp, SO=self.SO, SP=self.SP)
+                                    symm_subgrp=self.symmpar_subgrp, SO=self.SO, SP=self.SP,
+                                    rot_mat=rot_mat_all)
 
     def convert_bands_input(self):
         """
@@ -810,7 +814,7 @@ def convert_transport_input(self):
             for it in things_to_save:
                 ar[self.transp_subgrp][it] = locals()[it]
 
-    def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP):
+    def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP, rot_mat=None):
         """
         Reads and stores symmetrisation data from symm_file, which can be is case.sympar or case.symqmc.
 
@@ -828,6 +832,10 @@ def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP):
              Is spin-orbit coupling considered?
         SP : integer
              Is the system spin-polarised?
+        rot_mat : list of ndarray, optional
+                  The local-frame rotation per orbit (rot_mat / rot_mat_all). When
+                  given, the symmetry matrices are rotated into that local solver
+                  frame.
 
         """
 
@@ -889,6 +897,11 @@ def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP):
         R.close()
         # Reading done!
 
+        # Symmetrization acts in each shell's local solver frame; antiunitary
+        # operations also conjugate the source-side rotation.
+        if rot_mat is not None:
+            rotate_to_local_frames(mat, perm, time_inv, orbits, rot_mat)
+
         # Save it to the HDF:
         with HDFArchive(self.hdf_file, 'a') as ar:
             if not (symm_subgrp in ar):

test/python/wien2k/CMakeLists.txt

@@ -91,3 +91,9 @@ add_test(NAME Py_wien2k_outband_python
          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})
+
+add_test(NAME Py_wien2k_symmetry_frames
+         COMMAND ${TRIQS_PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/wien2k_symmetry_frames.py
+         WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+set_property(TEST Py_wien2k_symmetry_frames APPEND PROPERTY ENVIRONMENT
+             PYTHONPATH=${PROJECT_BINARY_DIR}/python:$ENV{PYTHONPATH} ${SANITIZER_RT_PRELOAD})

test/python/wien2k/wien2k_symmetry_frames.py

@@ -0,0 +1,53 @@
+################################################################################
+# TRIQS: a Toolbox for Research in Interacting Quantum Systems
+# (GPL-3.0-or-later)
+################################################################################
+
+import numpy as np
+
+from triqs_dftkit.wien2k._symmetry import rotate_to_local_frames
+
+
+# Two angular momenta on each atom make an atom/sort-only lookup ambiguous.
+orbits = [
+    dict(atom=1, sort=1, l=0, dim=1),
+    dict(atom=1, sort=1, l=1, dim=2),
+    dict(atom=2, sort=1, l=0, dim=1),
+    dict(atom=2, sort=1, l=1, dim=2),
+]
+r0 = np.array([[1j]])
+r1 = np.array([[0, 1], [1j, 0]], dtype=complex)
+r2 = np.array([[1.]])
+r3 = np.array([[1, 1j], [1j, 1]], dtype=complex) / np.sqrt(2)
+rotations = [r0, r1, r2, r3]
+permutations = [[2, 1], [1, 2]]
+rng = np.random.default_rng(1729)
+
+
+def random_matrix(dim):
+    return rng.normal(size=(dim, dim)) + 1j * rng.normal(size=(dim, dim))
+
+
+def action(matrix, density, antiunitary):
+    source = density.conj() if antiunitary else density
+    return matrix @ source @ matrix.conj().T
+
+for antiunitary in (False, True):
+    matrices = [[np.linalg.qr(random_matrix(orbit['dim']))[0] for orbit in orbits]
+                for _ in permutations]
+    global_matrices = [[matrix.copy() for matrix in row] for row in matrices]
+    rotate_to_local_frames(matrices, permutations,
+                           [antiunitary, antiunitary], orbits, rotations)
+    # Independent covariance oracle: acting in the global frame and then
+    # changing frame must equal changing frame first and acting locally.
+    for isym, targets in enumerate(((2, 3, 0, 1), (0, 1, 2, 3))):
+        for iorb, jorb in enumerate(targets):
+            seed = random_matrix(orbits[iorb]['dim'])
+            density = seed + seed.conj().T
+            local_density = rotations[iorb].conj().T @ density @ rotations[iorb]
+            global_result = action(global_matrices[isym][iorb], density, antiunitary)
+            expected = rotations[jorb].conj().T @ global_result @ rotations[jorb]
+            actual = action(matrices[isym][iorb], local_density, antiunitary)
+            np.testing.assert_allclose(actual, expected, atol=1e-13, rtol=0)
+
+print('wien2k_symmetry_frames: ok')