Add test for calculate_diagonalization_matrix

python/sumk_dft.py

@@ -1088,7 +1088,7 @@ class SumkDFT(object):
                 full_structure, ish, mesh=G[ish].mesh.copy(), show_warnings=threshold,
                 gf_function=type(G[ish]._first()), space_from='sumk', space_to='solver')
             for ish in range(self.n_inequiv_shells)]
-        #print 'c'
+
         if analyse_deg_shells:
             self.analyse_deg_shells(G_transformed, threshold, include_shells)
         return G_transformed
@@ -1359,6 +1359,11 @@ class SumkDFT(object):
                The transformation matrix for each spin-block in the correlated shell
         """
 
+        if self.block_structure.transformation:
+            mpi.report(
+                    "calculate_diagonalization_matrix: requires block_structure.transformation = None.")
+            return 0
+
         # Use all shells
         if shells is None:
             shells = range(self.n_corr_shells)

test/CMakeLists.txt

@@ -5,7 +5,7 @@ file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/${all_h5_files} DESTINATION ${CMAKE_CURREN
 FILE(COPY SrVO3.pmat SrVO3.struct SrVO3.outputs SrVO3.oubwin SrVO3.ctqmcout SrVO3.symqmc SrVO3.sympar SrVO3.parproj SrIrO3_rot.h5 hk_convert_hamiltonian.hk LaVO3-Pnma_hr.dat LaVO3-Pnma.inp DESTINATION ${CMAKE_CURRENT_BINARY_DIR})
 
 # List all tests
-set(all_tests wien2k_convert hk_convert w90_convert sumkdft_basic srvo3_Gloc srvo3_transp sigma_from_file blockstructure analyse_block_structure_from_gf analyse_block_structure_from_gf2)
+set(all_tests wien2k_convert hk_convert w90_convert sumkdft_basic srvo3_Gloc srvo3_transp sigma_from_file blockstructure analyse_block_structure_from_gf analyse_block_structure_from_gf2 basis_transformation)
 
 set(python_executable python)
 

test/basis_transformation.py

@@ -0,0 +1,70 @@
+from pytriqs.utility.comparison_tests import *
+from triqs_dft_tools.sumk_dft import *
+import numpy as np
+
+def is_diagonal_matrix(M):
+    return abs(np.sum(M-np.diag(np.diagonal(M)))) < 1e-10
+
+def call_diagonalize(SK):
+    SK.block_structure.transformation = None
+    t_sumk_eal = SK.calculate_diagonalization_matrix(prop_to_be_diagonal='eal', calc_in_solver_blocks=False, write_to_blockstructure = True)
+    SK.block_structure.transformation = None
+    t_solver_eal = SK.calculate_diagonalization_matrix(prop_to_be_diagonal='eal', calc_in_solver_blocks=True, write_to_blockstructure = True)
+    SK.block_structure.transformation = None
+    t_solver_dm = SK.calculate_diagonalization_matrix(prop_to_be_diagonal='dm', calc_in_solver_blocks=False, write_to_blockstructure = True)
+    SK.block_structure.transformation = None
+    t_sumk_dm = SK.calculate_diagonalization_matrix(prop_to_be_diagonal='dm', calc_in_solver_blocks=True, write_to_blockstructure = True)
+    SK.block_structure.transformation = None
+    return t_sumk_eal, t_solver_eal, t_sumk_dm, t_solver_dm
+
+SK = SumkDFT(hdf_file = 'SrVO3.h5', use_dft_blocks=True)
+
+# only eal and dm are allowed
+SK.block_structure.transformation = None
+assert not SK.calculate_diagonalization_matrix(prop_to_be_diagonal='test')
+
+# check for shell index
+assert not SK.calculate_diagonalization_matrix(shells = [15])
+
+# calling the function twice leads to block_structure.transformation already being set
+SK.calculate_diagonalization_matrix()
+assert not SK.calculate_diagonalization_matrix()
+SK.block_structure.transformation = None
+
+# Check writing to block_structure
+SK.calculate_diagonalization_matrix(write_to_blockstructure=False)
+assert SK.block_structure.transformation is None
+SK.block_structure.transformation = None
+SK.calculate_diagonalization_matrix(write_to_blockstructure=True)
+assert SK.block_structure.transformation is not None
+SK.block_structure.transformation = None
+
+t_sumk_eal, t_solver_eal, t_sumk_dm, t_solver_dm = call_diagonalize(SK)
+
+# All matrices should be identities
+for orb in range(SK.n_corr_shells):
+    for block in t_solver_eal[orb]:
+        assert_arrays_are_close(t_sumk_eal[orb][block],np.identity(3), precision=1e-6)
+        assert_arrays_are_close(t_sumk_dm[orb][block],np.identity(3), precision=1e-6)
+        assert_arrays_are_close(t_solver_eal[orb][block],np.identity(3), precision=1e-6)
+        assert_arrays_are_close(t_solver_dm[orb][block],np.identity(3), precision=1e-6)
+
+SK = SumkDFT(hdf_file = 'w90_convert.ref.h5', use_dft_blocks=True)
+
+t_sumk_eal, t_solver_eal, t_sumk_dm, t_solver_dm = call_diagonalize(SK)
+
+# In this example solver and sumk should be the same
+for orb in range(SK.n_corr_shells):
+    for block in t_solver_eal[orb]:
+        assert_arrays_are_close(t_sumk_eal[orb][block],t_solver_eal[orb][block], precision=1e-6)
+        assert_arrays_are_close(t_sumk_dm[orb][block],t_solver_dm[orb][block], precision=1e-6)
+
+# Check if transformations make eal and dm really diagonal
+eal = SK.eff_atomic_levels()[0]
+for e in eal:
+    assert is_diagonal_matrix(np.dot(np.dot(t_solver_eal[0][e], eal[e].conj().T),t_solver_eal[0][e].conj().T))
+
+dm = SK.density_matrix(method='using_point_integration')
+for dmi in dm:
+    for e in dmi:
+        assert is_diagonal_matrix(np.dot(np.dot(t_solver_dm[0][e], dmi[e].conj().T),t_solver_dm[0][e].conj().T))