diff --git a/.repo-hygiene-allow b/.repo-hygiene-allow deleted file mode 100644 index 86f9197..0000000 --- a/.repo-hygiene-allow +++ /dev/null @@ -1 +0,0 @@ -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/python/triqs_dftkit/wien2k/_symmetry.py b/python/triqs_dftkit/wien2k/_symmetry.py new file mode 100644 index 0000000..5c140ef --- /dev/null +++ b/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 diff --git a/python/triqs_dftkit/wien2k/converter.py b/python/triqs_dftkit/wien2k/converter.py index 37e25f3..d0ac4f6 100644 --- a/python/triqs_dftkit/wien2k/converter.py +++ b/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): diff --git a/test/python/wien2k/CMakeLists.txt b/test/python/wien2k/CMakeLists.txt index f94fffa..6c1b428 100644 --- a/test/python/wien2k/CMakeLists.txt +++ b/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}) diff --git a/test/python/wien2k/wien2k_convert.ref.h5 b/test/python/wien2k/wien2k_convert.ref.h5 index f837d17..870b04b 100644 Binary files a/test/python/wien2k/wien2k_convert.ref.h5 and b/test/python/wien2k/wien2k_convert.ref.h5 differ diff --git a/test/python/wien2k/wien2k_soc_convert.ref.h5 b/test/python/wien2k/wien2k_soc_convert.ref.h5 index 84f56d7..705f24b 100644 Binary files a/test/python/wien2k/wien2k_soc_convert.ref.h5 and b/test/python/wien2k/wien2k_soc_convert.ref.h5 differ diff --git a/test/python/wien2k/wien2k_symmetry_frames.py b/test/python/wien2k/wien2k_symmetry_frames.py new file mode 100644 index 0000000..fac6f17 --- /dev/null +++ b/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')