refactor & fix: merge the two read functions and small fixes

fix: add asserts for spin calculations and disent

python/triqs_dft_tools/sumk_dft_transport.py

@@ -4,6 +4,7 @@
 # TRIQS: a Toolbox for Research in Interacting Quantum Systems
 #
 # Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
+# Copyright (c) 2022-2023 Simons Foundation
 #
 # TRIQS is free software: you can redistribute it and/or modify it under the
 # terms of the GNU General Public License as published by the Free Software
@@ -18,6 +19,8 @@
 # You should have received a copy of the GNU General Public License along with
 # TRIQS. If not, see <http://www.gnu.org/licenses/>.
 #
+# Authors: M. Aichhorn, S. Beck, A. Hampel, L. Pourovskii, V. Vildosola
+
 ##########################################################################
 import sys
 import numpy
@@ -25,7 +28,6 @@ from warnings import warn
 from triqs.gf import *
 import triqs.utility.mpi as mpi
 from .symmetry import *
-from .sumk_dft import SumkDFT
 import scipy.constants as cst
 import os.path
 
@@ -38,12 +40,10 @@ __all__ = ['transport_distribution', 'conductivity_and_seebeck', 'write_output_t
 def read_transport_input_from_hdf(sum_k):
     r"""
     Reads the data for transport calculations from the hdf5 archive.
-
     Parameters
     ----------
     sum_k : sum_k object
             triqs SumkDFT object
-
     Returns
     -------
     sum_k : sum_k object
@@ -214,7 +214,7 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname='
 
     # set-up k mesh depending on input shape
     # read in transport input and some checks
-    read_transport_input_from_hdf_wannier90(sum_k)
+    read_transport_input_from_hdf(sum_k)
 
     # first check for right formatting of sum_k.nk_optics
     assert len(nk_optics) in [1, 3], '"nk_optics" must be given as three integers or one float'
@@ -230,6 +230,9 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname='
     else:
         nk_x, nk_y, nk_z = nk_optics
 
+    # check for spin calculation (not supported)
+    assert sum_k.SP == 0, 'spin dependent transport calculations are not supported.'
+
     n_orb = numpy.max([sum_k.n_orbitals[ik][0] for ik in range(sum_k.n_k)])
 
     # temporarily recompute the following quantities on a different mesh
@@ -268,6 +271,8 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname='
         grid = wb.Grid(wberri, NKdiv=1, NKFFT=[nk_x, nk_y, nk_z])
         dataK = wb.data_K.Data_K(wberri, dK=shift_gamma, grid=grid, fftlib='numpy')
 
+        assert dataK.HH_K.shape == hopping[:, 0, :, :].shape, 'wberri / wannier Hamiltonian has different number of orbitals than SumK object. Disentanglement is not supported as of now.'
+
         # read in hoppings and proj_mat
         if oc_basis == 'h':
             hopping[:, 0, range(hopping.shape[2]), range(hopping.shape[3])] = dataK.E_K
@@ -769,6 +774,9 @@ def transport_coefficient(Gamma_w, omega, Om_mesh, spin_polarization, direction,
         Transport coefficient.
     """
 
+    from scipy.interpolate import interp1d
+    from scipy.integrate import simpson, quad
+
     if not (mpi.is_master_node()):
         return
 
@@ -791,7 +799,7 @@ def transport_coefficient(Gamma_w, omega, Om_mesh, spin_polarization, direction,
             A = A[0]
         elif method == 'simps':
             # simpson rule for w-grid
-            A = simps(A_int, omega)
+            A = simpson(A_int, omega)
         elif method == 'trapz':
             # trapezoidal rule for w-grid
             A = numpy.trapz(A_int, omega)

test/python/w90_optics_Sr2RuO4.py

@@ -4,8 +4,7 @@ import numpy as np
 import triqs.utility.mpi as mpi
 from h5 import HDFArchive
 from triqs.gf import MeshReFreq
-from triqs.utility.h5diff import h5diff
+from triqs.utility import h5diff
-
 from triqs_dft_tools.sumk_dft import SumkDFT
 from triqs_dft_tools.sumk_dft_transport import transport_distribution, init_spectroscopy, conductivity_and_seebeck, write_output_to_hdf
 from triqs_dft_tools.converters.wannier90 import Wannier90Converter
@@ -74,4 +73,8 @@ if mpi.is_master_node():
                    'seebeck': seebeck, 'optic_cond': optic_cond, 'kappa': kappa}
     write_output_to_hdf(sum_k, output_dict, 'transp_output')
 
-    h5diff(h5_archive, "sr2ruo4_transp.ref.h5")
+    # comparison of the output transport data
+    # velocities can differ depending on LAPACK diagonalization
+    out = HDFArchive(h5_archive,'r')
+    ref = HDFArchive('sr2ruo4_transp.ref.h5', 'r')
+    h5diff.compare('', out['transp_output'], ref['transp_output'], 0, 1e-8)