diff --git a/python/triqs_dft_tools/sumk_dft_tools.py b/python/triqs_dft_tools/sumk_dft_tools.py index 4ec73826..13eb053c 100644 --- a/python/triqs_dft_tools/sumk_dft_tools.py +++ b/python/triqs_dft_tools/sumk_dft_tools.py @@ -31,8 +31,6 @@ from .symmetry import * from .sumk_dft import SumkDFT from scipy.integrate import * from scipy.interpolate import * -from scipy import constants as constants -from itertools import product if not hasattr(numpy, 'full'): # polyfill full for older numpy: @@ -1060,546 +1058,3 @@ class SumkDFTTools(SumkDFT): f.close() -# ----------------- transport ----------------------- - - def read_transport_input_from_hdf(self): - r""" - Reads the data for transport calculations from the hdf5 archive. - """ - assert self.dft_code in ('wien2k','elk'), "Transport has only been implemented for wien2k and elk inputs" - thingstoread = ['band_window_optics', 'velocities_k'] - self.read_input_from_hdf( - subgrp=self.transp_data, things_to_read=thingstoread) - if(self.dft_code=="wien2k"): - thingstoread = ['band_window', 'lattice_angles', 'lattice_constants', - 'lattice_type', 'n_symmetries', 'rot_symmetries'] - elif(self.dft_code=="elk"): - thingstoread = ['band_window', 'n_symmetries', - 'rot_symmetries','cell_vol'] - self.read_input_from_hdf( - subgrp=self.misc_data, things_to_read=thingstoread) - if(self.dft_code=="wien2k"): - self.cell_vol = self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1] - - def read_transport_input_from_hdf_wannier90(self): - r""" - Reads the data for transport calculations from the hdf5 archive. - """ - thingstoread = ['band_window_optics', 'nk_optics'] - self.read_input_from_hdf( - subgrp=self.transp_data, things_to_read=thingstoread) - thingstoread = ['band_window', 'n_symmetries', 'rot_symmetries'] - self.read_input_from_hdf( - subgrp=self.misc_data, things_to_read=thingstoread) - - def cellvolume(self, lattice_type, lattice_constants, latticeangle): - r""" - Determines the conventional und primitive unit cell volumes. - - Parameters - ---------- - lattice_type : string - Lattice type according to the Wien2k convention (P, F, B, R, H, CXY, CYZ, CXZ). - lattice_constants : list of double - Lattice constants (a, b, c). - lattice angles : list of double - Lattice angles (:math:`\alpha, \beta, \gamma`). - - Returns - ------- - vol_c : double - Conventional unit cell volume. - vol_p : double - Primitive unit cell volume. - """ - - a = lattice_constants[0] - b = lattice_constants[1] - c = lattice_constants[2] - c_al = numpy.cos(latticeangle[0]) - c_be = numpy.cos(latticeangle[1]) - c_ga = numpy.cos(latticeangle[2]) - vol_c = a * b * c * \ - numpy.sqrt(1 + 2 * c_al * c_be * c_ga - - c_al ** 2 - c_be ** 2 - c_ga ** 2) - - det = {"P": 1, "F": 4, "B": 2, "R": 3, - "H": 1, "CXY": 2, "CYZ": 2, "CXZ": 2} - vol_p = vol_c / det[lattice_type] - - return vol_c, vol_p - - # Uses .data of only GfReFreq objects. - def transport_distribution(self, beta, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0, code='wien2k'): - r""" - Calculates the transport distribution - - .. math:: - \Gamma_{\alpha\beta}\left(\omega+\Omega/2, \omega-\Omega/2\right) = \frac{1}{V} \sum_k Tr\left(v_{k,\alpha}A_{k}(\omega+\Omega/2)v_{k,\beta}A_{k}\left(\omega-\Omega/2\right)\right) - - in the direction :math:`\alpha\beta`. The velocities :math:`v_{k}` are read from the transport subgroup of the hdf5 archive. - - Parameters - ---------- - - beta : double - Inverse temperature :math:`\beta`. - directions : list of double, optional - :math:`\alpha\beta` e.g.: ['xx','yy','zz','xy','xz','yz']. - energy_window : list of double, optional - Specifies the upper and lower limit of the frequency integration for :math:`\Omega=0.0`. The window is automatically enlarged by the largest :math:`\Omega` value, - hence the integration is performed in the interval [energy_window[0]-max(Om_mesh), energy_window[1]+max(Om_mesh)]. - Om_mesh : list of double, optional - :math:`\Omega` frequency mesh of the optical conductivity. For the conductivity and the Seebeck coefficient :math:`\Omega=0.0` has to be - part of the mesh. In the current version Om_mesh is repined to the mesh provided by the self-energy! The actual mesh is printed on the screen and stored as - member Om_mesh. - with_Sigma : boolean, optional - Determines whether the calculation is performed with or without self energy. If this parameter is set to False the self energy is set to zero (i.e. the DFT band - structure :math:`A(k,\omega)` is used). Note: For with_Sigma=False it is necessary to specify the parameters energy_window, n_om and broadening. - n_om : integer, optional - Number of equidistant frequency points in the interval [energy_window[0]-max(Om_mesh), energy_window[1]+max(Om_mesh)]. This parameters is only used if - with_Sigma = False. - broadening : double, optional - Lorentzian broadening. It is necessary to specify the boradening if with_Sigma = False, otherwise this parameter can be set to 0.0. - code : string - DFT code from which velocities are being read. Options: 'wien2k', 'wannier90' - """ - BOHRTOANG = constants.physical_constants['Bohr radius'][0]/constants.angstrom - HARTREETOEV = constants.physical_constants['Hartree energy'][0]/constants.eV - n_inequiv_spin_blocks = self.SP + 1 - self.SO - - if code in ('wien2k'): - # Check if wien converter was called and read transport subgroup form - # hdf file - if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file, 'r') - if not (self.transp_data in ar): - raise IOError("transport_distribution: No %s subgroup in hdf file found! Call convert_transp_input first." % self.transp_data) - # check if outputs file was converted - if not ('n_symmetries' in ar['dft_misc_input']): - raise IOError("transport_distribution: n_symmetries missing. Check if case.outputs file is present and call convert_misc_input() or convert_dft_input().") - - self.read_transport_input_from_hdf() - cell_volume = self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1] - n_symmetries = self.n_symmetries - - elif code in ('wannier90'): - # check if spin-unpolarized - assert n_inequiv_spin_blocks == 1, "Spin-polarized optical conductivity calculations not implemented with Wannier90" - - # read in transport input - self.read_transport_input_from_hdf_wannier90() - # checks for right formatting of self.nk_optics - assert len(self.nk_optics) in [1,3], '"nk_optics" must be given as three integers or one float' - if len(self.nk_optics) == 1: assert np.array(list(self.nk_optics)).dtype in (int, float), '"nk_optics" single value must be float or integer' - if len(self.nk_optics) == 3: assert np.array(list(self.nk_optics)).dtype == int, '"nk_optics" mesh must be integers' - n_symmetries = 1 - - # calculate velocity - #wberri = wb.System_w90('/mnt/home/sbeck/Dropbox/ccqlin030/sro/I4_mmm_prim/wan_conv_12_v/sro', berry=True) - pathname = './' - seedname = 'sro' - fermi = 0. - if len(self.nk_optics) == 1: - interpolate_factor = self.nk_optics[0] - nk_x, nk_y, nk_z = list(map(lambda i: int(numpy.ceil(interpolate_factor * len(set(self.kpts[:,i])))), range(3))) - else: - nk_x, nk_y, nk_z = self.nk_optics - #nk_x, nk_y, nk_z = 10 - #nk_x, nk_y, nk_z = 12 - #nk_x, nk_y, nk_z = 34 - n_orb = numpy.max([self.n_orbitals[ik][0] for ik in range(self.n_k)]) - shift_gamma = [0.0,0.0,0.0] - #shift_gamma = [0.015,0.015,0.015] - things_to_modify = {'bz_weights': None, 'hopping': None, 'kpt_weights': None, 'kpts': None, - 'n_k': None, 'n_orbitals': None, 'proj_mat': None, 'band_window': None, 'band_window_optics': None} - things_to_store = dict.fromkeys(things_to_modify, None) - - # initialize variables - n_kpts = nk_x * nk_y * nk_z - kpts = numpy.zeros((n_kpts, 3)) - hopping = numpy.zeros((n_kpts, 1, n_orb, n_orb), dtype=complex) - proj_mat = numpy.zeros(numpy.shape(hopping[:,0,0,0]) + numpy.shape(self.proj_mat[0,:]), dtype=complex) - if mpi.is_master_node(): - print(hopping.shape, self.proj_mat.shape, numpy.shape(hopping[:,0,0,0]) + numpy.shape(self.proj_mat[0,:])) - # simple modifications - things_to_modify['n_k'] = n_kpts - things_to_modify['n_orbitals'] = numpy.full((n_kpts, 1), n_orb) - for key in ['bz_weights', 'kpt_weights']: - things_to_modify[key] = numpy.full(n_kpts, 1/n_kpts) - n_inequiv_spin_blocks = self.SP + 1 - self.SO - for key in ['band_window', 'band_window_optics']: - things_to_modify[key] = [numpy.full((n_kpts, 2), self.band_window[isp][0]) for isp in range(n_inequiv_spin_blocks)] - - velocities_k = None - cell_volume = None - kpts = None - - if mpi.is_master_node(): - # try wannierberri import - try: - import wannierberri as wb - except ImportError: - print('ImportError: WannierBerri needs to be installed to run test "Py_w90_optics_Sr2RuO4"') - try: - mpi.MPI.COMM_WORLD.Abort(1) - except: - sys.exit() - # initialize WannierBerri system - wberri = wb.System_w90(pathname + seedname, berry=True) - grid = wb.Grid(wberri, NKdiv=1, NKFFT=[nk_x, nk_y, nk_z]) - dataK = Data_K(wberri, dK=shift_gamma, grid=grid) - - # construct velocities from dataK - V_H_diag = numpy.zeros(numpy.shape(dataK.V_H), dtype=complex) - V_H_diag[:, range(V_H_diag.shape[1]), range(V_H_diag.shape[1]), :] = numpy.diagonal(dataK.V_H[:,:,:,:],axis1=1, axis2=2).transpose(0,2,1).copy() - velocities_k = ( V_H_diag - dataK.A_Hbar * 1j*( dataK.E_K[:,None,:,None] - dataK.E_K[:,:,None,None] ) ) / HARTREETOEV / BOHRTOANG - #velocities_k = V_H_diag / HARTREETOEV / BOHRTOANG - - # read in hoppings and proj_mat - hopping[:,0,range(hopping.shape[2]),range(hopping.shape[3])] = dataK.E_K - for isp in range(n_inequiv_spin_blocks): - iorb = 0 - for icrsh in range(self.n_corr_shells): - dim = self.corr_shells[icrsh]['dim'] - proj_mat[:,isp,icrsh,0:dim,:] = dataK.UU_K[:,iorb:iorb+dim,:] - iorb += dim - - # read in rest from dataK - cell_volume = dataK.cell_volume / BOHRTOANG ** 3 - kpts = dataK.kpoints_all - - # broadcast everything - velocities_k = mpi.bcast(velocities_k) - cell_volume = mpi.bcast(cell_volume) - kpts = mpi.bcast(kpts) - hopping = mpi.bcast(hopping) - proj_mat = mpi.bcast(proj_mat) - - # upgrade sumk quantities for interpolation - things_to_modify['kpts'] = kpts - things_to_modify['hopping'] = hopping - things_to_modify['proj_mat'] = proj_mat - mpi.barrier() - - if mpi.is_master_node(): - print(self.n_k, nk_x, nk_y, nk_z) - for key in things_to_modify: - things_to_store[key] = getattr(self, key) - setattr(self, key, things_to_modify[key]) - #if mpi.is_master_node(): - #print(key, things_to_store[key] ) - #print(getattr(self, key)) - # write velocities to file - if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file, 'a') - ar['dft_transp_input']['velocities_k'] = velocities_k - - if mpi.is_master_node(): - # k-dependent-projections. - # to be checked. But this should be obsolete atm, works for both cases - # k_dep_projection is nowhere used - # assert sum_k.k_dep_projection == 0, "transport_distribution: k dependent projection is not implemented!" - # positive Om_mesh - assert all( - Om >= 0.0 for Om in Om_mesh), "transport_distribution: Om_mesh should not contain negative values!" - - # Check if energy_window is sufficiently large and correct - - if (energy_window[0] >= energy_window[1] or energy_window[0] >= 0 or energy_window[1] <= 0): - assert 0, "transport_distribution: energy_window wrong!" - - if (abs(self.fermi_dis(energy_window[0], beta) * self.fermi_dis(-energy_window[0], beta)) > 1e-5 - or abs(self.fermi_dis(energy_window[1], beta) * self.fermi_dis(-energy_window[1], beta)) > 1e-5): - mpi.report( - "\n####################################################################") - mpi.report( - "transport_distribution: WARNING - energy window might be too narrow!") - mpi.report( - "####################################################################\n") - - # up and down are equivalent if SP = 0 - self.directions = directions - dir_to_int = {'x': 0, 'y': 1, 'z': 2} - - # calculate A(k,w) - ####################################### - - # Define mesh for Green's function and in the specified energy window - if (with_Sigma == True): - self.omega = numpy.array([round(x.real, 12) - for x in self.mesh]) - mesh = None - mu = self.chemical_potential - n_om = len(self.omega) - mpi.report("Using omega mesh provided by Sigma!") - - if energy_window: - # Find according window in Sigma mesh - ioffset = numpy.sum( - self.omega < energy_window[0] - max(Om_mesh)) - self.omega = self.omega[numpy.logical_and(self.omega >= energy_window[ - 0] - max(Om_mesh), self.omega <= energy_window[1] + max(Om_mesh))] - n_om = len(self.omega) - - # Truncate Sigma to given omega window - # In the future there should be an option in gf to manipulate the mesh (e.g. truncate) directly. - # For now we stick with this: - for icrsh in range(self.n_corr_shells): - Sigma_save = self.Sigma_imp[icrsh].copy() - spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] - glist = lambda: [GfReFreq(target_shape=(block_dim, block_dim), window=(self.omega[ - 0], self.omega[-1]), n_points=n_om) for block, block_dim in self.gf_struct_sumk[icrsh]] - self.Sigma_imp[icrsh] = BlockGf( - name_list=spn, block_list=glist(), make_copies=False) - for i, g in self.Sigma_imp[icrsh]: - for iL in g.indices[0]: - for iR in g.indices[0]: - for iom in range(n_om): - g.data[iom, int(iL), int(iR)] = Sigma_save[ - i].data[ioffset + iom, int(iL), int(iR)] - else: - assert n_om is not None, "transport_distribution: Number of omega points (n_om) needed to calculate transport distribution!" - assert energy_window is not None, "transport_distribution: Energy window needed to calculate transport distribution!" - assert broadening != 0.0 and broadening is not None, "transport_distribution: Broadening necessary to calculate transport distribution!" - self.omega = numpy.linspace( - energy_window[0] - max(Om_mesh), energy_window[1] + max(Om_mesh), n_om) - mesh = MeshReFreq(energy_window[0] - - max(Om_mesh), energy_window[1] + max(Om_mesh), n_om) - mu = 0.0 - - # Define mesh for optic conductivity - d_omega = round(numpy.abs(self.omega[0] - self.omega[1]), 12) - iOm_mesh = numpy.array([round((Om / d_omega), 0) for Om in Om_mesh]) - self.Om_mesh = iOm_mesh * d_omega - - if mpi.is_master_node(): - print("Chemical potential: ", mu) - print("Using n_om = %s points in the energy_window [%s,%s]" % (n_om, self.omega[0], self.omega[-1]), end=' ') - print("where the omega vector is:") - print(self.omega) - print("Calculation requested for Omega mesh: ", numpy.array(Om_mesh)) - print("Omega mesh automatically repined to: ", self.Om_mesh) - - self.Gamma_w = {direction: numpy.zeros( - (len(self.Om_mesh), n_om), dtype=float) for direction in self.directions} - max_orb = numpy.max([self.n_orbitals[ik][0] for ik in range(self.n_k)]) - #Akw_write = numpy.zeros((self.n_k, max_orb, max_orb, n_om), dtype=numpy.complex_) - - # Sum over all k-points - ikarray = numpy.array(list(range(self.n_k))) - for ik in mpi.slice_array(ikarray): - # Calculate G_w for ik and initialize A_kw - G_w = self.lattice_gf(ik, mu, broadening=broadening, mesh=mesh, with_Sigma=with_Sigma) - A_kw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=complex) - for isp in range(n_inequiv_spin_blocks)] - - for isp in range(n_inequiv_spin_blocks): - # copy data from G_w (swapaxes is used to have omega in the 3rd - # dimension) - A_kw[isp] = copy.deepcopy(G_w[self.spin_block_names[self.SO][ - isp]].data.swapaxes(0, 1).swapaxes(1, 2)) - # calculate A(k,w) for each frequency - for iw in range(n_om): - A_kw[isp][:, :, iw] = -1.0 / (2.0 * numpy.pi * 1j) * ( - A_kw[isp][:, :, iw] - numpy.conjugate(numpy.transpose(A_kw[isp][:, :, iw]))) - #Akw_write[ik] = A_kw[isp].copy() * self.bz_weights[ik] - - b_min = max(self.band_window[isp][ - ik, 0], self.band_window_optics[isp][ik, 0]) - b_max = min(self.band_window[isp][ - ik, 1], self.band_window_optics[isp][ik, 1]) - A_i = slice( - b_min - self.band_window[isp][ik, 0], b_max - self.band_window[isp][ik, 0] + 1) - v_i = slice(b_min - self.band_window_optics[isp][ - ik, 0], b_max - self.band_window_optics[isp][ik, 0] + 1) - - # loop over all symmetries - for R in self.rot_symmetries: - # get transformed velocity under symmetry R - if code in ('wien2k'): - vel_R = copy.deepcopy(self.velocities_k[isp][ik]) - elif code in ('wannier90'): - vel_R = copy.deepcopy(velocities_k[ik]) - for nu1 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1): - for nu2 in range(self.band_window_optics[isp][ik, 1] - self.band_window_optics[isp][ik, 0] + 1): - vel_R[nu1][nu2][:] = numpy.dot( - R, vel_R[nu1][nu2][:]) - - # calculate Gamma_w for each direction from the velocities - # vel_R and the spectral function A_kw - for direction in self.directions: - for iw in range(n_om): - for iq in range(len(self.Om_mesh)): - if(iw + iOm_mesh[iq] >= n_om or self.omega[iw] < -self.Om_mesh[iq] + energy_window[0] or self.omega[iw] > self.Om_mesh[iq] + energy_window[1]): - continue - - self.Gamma_w[direction][iq, iw] += (numpy.dot(numpy.dot(numpy.dot(vel_R[v_i, v_i, dir_to_int[direction[0]]], - A_kw[isp][A_i, A_i, int(iw + iOm_mesh[iq])]), vel_R[v_i, v_i, dir_to_int[direction[1]]]), - A_kw[isp][A_i, A_i, iw]).trace().real * self.bz_weights[ik]) - - #Akw_write = mpi.all_reduce(mpi.world, Akw_write, lambda x, y: x + y) - #mpi.barrier() - #if mpi.is_master_node(): - # ar = HDFArchive(self.hdf_file, 'a') - # ar.create_group('Akw') - # ar['Akw'] = numpy.sum(numpy.trace(Akw_write, axis1=1, axis2=2), axis=0) - - for direction in self.directions: - self.Gamma_w[direction] = (mpi.all_reduce(mpi.world, self.Gamma_w[direction], lambda x, y: x + y) / self.cell_vol / self.n_symmetries) - - - def transport_coefficient(self, direction, iq, n, beta, method=None): - r""" - Calculates the transport coefficient A_n in a given direction for a given :math:`\Omega`. The required members (Gamma_w, directions, Om_mesh) have to be obtained first - by calling the function :meth:`transport_distribution `. For n>0 A is set to NaN if :math:`\Omega` is not 0.0. - - Parameters - ---------- - direction : string - :math:`\alpha\beta` e.g.: 'xx','yy','zz','xy','xz','yz'. - iq : integer - Index of :math:`\Omega` point in the member Om_mesh. - n : integer - Number of the desired moment of the transport distribution. - beta : double - Inverse temperature :math:`\beta`. - method : string - Integration method: cubic spline and scipy.integrate.quad ('quad'), simpson rule ('simps'), trapezoidal rule ('trapz'), rectangular integration (otherwise) - Note that the sampling points of the the self-energy are used! - - Returns - ------- - A : double - Transport coefficient. - """ - - if not (mpi.is_master_node()): - return - - assert hasattr( - self, 'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!" - - if (self.Om_mesh[iq] == 0.0 or n == 0.0): - A = 0.0 - # setup the integrand - if (self.Om_mesh[iq] == 0.0): - A_int = self.Gamma_w[direction][iq] * (self.fermi_dis( - self.omega, beta) * self.fermi_dis(-self.omega, beta)) * (self.omega * beta)**n - elif (n == 0.0): - A_int = self.Gamma_w[direction][iq] * (self.fermi_dis(self.omega, beta) - self.fermi_dis( - self.omega + self.Om_mesh[iq], beta)) / (self.Om_mesh[iq] * beta) - - # w-integration - if method == 'quad': - # quad on interpolated w-points with cubic spline - A_int_interp = interp1d(self.omega, A_int, kind='cubic') - A = quad(A_int_interp, min(self.omega), max(self.omega), - epsabs=1.0e-12, epsrel=1.0e-12, limit=500) - A = A[0] - elif method == 'simps': - # simpson rule for w-grid - A = simps(A_int, self.omega) - elif method == 'trapz': - # trapezoidal rule for w-grid - A = numpy.trapz(A_int, self.omega) - else: - # rectangular integration for w-grid (orignal implementation) - d_w = self.omega[1] - self.omega[0] - for iw in range(self.Gamma_w[direction].shape[1]): - A += A_int[iw] * d_w - A = A * numpy.pi * (2.0 - self.SP) - else: - A = numpy.nan - return A - - def conductivity_and_seebeck(self, beta, method=None): - r""" - Calculates the Seebeck coefficient and the optical conductivity by calling - :meth:`transport_coefficient `. - The required members (Gamma_w, directions, Om_mesh) have to be obtained first by calling the function - :meth:`transport_distribution `. - - Parameters - ---------- - beta : double - Inverse temperature :math:`\beta`. - - Returns - ------- - optic_cond : dictionary of double vectors - Optical conductivity in each direction and frequency given by Om_mesh. - - seebeck : dictionary of double - Seebeck coefficient in each direction. If zero is not present in Om_mesh the Seebeck coefficient is set to NaN. - - kappa : dictionary of double. - thermal conductivity in each direction. If zero is not present in Om_mesh the thermal conductivity is set to NaN - """ - - if not (mpi.is_master_node()): - return - - assert hasattr( - self, 'Gamma_w'), "conductivity_and_seebeck: Run transport_distribution first or load data from h5!" - n_q = self.Gamma_w[self.directions[0]].shape[0] - - A0 = {direction: numpy.full((n_q,), numpy.nan) - for direction in self.directions} - A1 = {direction: numpy.full((n_q,), numpy.nan) - for direction in self.directions} - A2 = {direction: numpy.full((n_q,), numpy.nan) - for direction in self.directions} - self.seebeck = {direction: numpy.nan for direction in self.directions} - self.kappa = {direction: numpy.nan for direction in self.directions} - self.optic_cond = {direction: numpy.full( - (n_q,), numpy.nan) for direction in self.directions} - - for direction in self.directions: - for iq in range(n_q): - A0[direction][iq] = self.transport_coefficient( - direction, iq=iq, n=0, beta=beta, method=method) - A1[direction][iq] = self.transport_coefficient( - direction, iq=iq, n=1, beta=beta, method=method) - A2[direction][iq] = self.transport_coefficient( - direction, iq=iq, n=2, beta=beta, method=method) - print("A_0 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A0[direction][iq])) - print("A_1 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A1[direction][iq])) - print("A_2 in direction %s for Omega = %.2f %e a.u." % (direction, self.Om_mesh[iq], A2[direction][iq])) - if ~numpy.isnan(A1[direction][iq]): - # Seebeck and kappa are overwritten if there is more than one Omega = - # 0 in Om_mesh - self.seebeck[direction] = - \ - A1[direction][iq] / A0[direction][iq] * 86.17 - self.kappa[direction] = A2[direction][iq] - A1[direction][iq]*A1[direction][iq]/A0[direction][iq] - self.kappa[direction] *= 293178.0 - self.optic_cond[direction] = beta * \ - A0[direction] * 10700.0 / numpy.pi - for iq in range(n_q): - print("Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % (direction, self.Om_mesh[iq], self.optic_cond[direction][iq])) - if not (numpy.isnan(A1[direction][iq])): - print("Seebeck in direction %s for Omega = 0.00 %f x 10^(-6) V/K" % (direction, self.seebeck[direction])) - print("kappa in direction %s for Omega = 0.00 %f W/(m * K)" % (direction, self.kappa[direction])) - - return self.optic_cond, self.seebeck, self.kappa - - - def fermi_dis(self, w, beta): - r""" - Fermi distribution. - - .. math:: - f(x) = 1/(e^x+1). - - Parameters - ---------- - w : double - frequency - beta : double - inverse temperature - - Returns - ------- - f : double - """ - return 1.0 / (numpy.exp(w * beta) + 1) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index d5f718fa..be3a62be 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -47,12 +47,19 @@ def read_transport_input_from_hdf(sum_k): sum_k : sum_k object triqs SumkDFT object """ + assert sum_k.dft_code in ('wien2k','elk'), "read_transport_input_from_hdf() is only implemented for wien2k and elk inputs" thingstoread = ['band_window_optics', 'velocities_k'] sum_k.read_input_from_hdf( subgrp=sum_k.transp_data, things_to_read=thingstoread) - thingstoread = ['band_window', 'lattice_angles', 'lattice_constants', - 'lattice_type', 'n_symmetries', 'rot_symmetries'] + if(sum_k.dft_code=="wien2k"): + thingstoread = ['band_window', 'lattice_angles', 'lattice_constants', + 'lattice_type', 'n_symmetries', 'rot_symmetries'] + elif(sum_k.dft_code=="elk"): + thingstoread = ['band_window', 'n_symmetries', + 'rot_symmetries','cell_vol'] sum_k.read_input_from_hdf(subgrp=sum_k.misc_data, things_to_read=thingstoread) + if(self.dft_code=="wien2k"): + self.cell_vol = self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1] return sum_k @@ -278,7 +285,7 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname=' kpts = dataK.kpoints_all # broadcast everything - cell_volume = mpi.bcast(cell_volume) + sum_k.cell_vol = mpi.bcast(cell_volume) kpts = mpi.bcast(kpts) hopping = mpi.bcast(hopping) proj_mat = mpi.bcast(proj_mat) @@ -319,7 +326,7 @@ def recompute_w90_input_on_different_mesh(sum_k, seedname, nk_optics, pathname=' ar = HDFArchive(sum_k.hdf_file, 'a') ar['dft_transp_input']['inverse_mass'] = inverse_mass - return sum_k, cell_volume, things_to_store + return sum_k, things_to_store # ----------------- transport ----------------------- @@ -341,15 +348,13 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}): ------- sum_k : sum_k object triqs SumkDFT object, interpolated - cell_volume : double - primitive unit cell volume """ n_inequiv_spin_blocks = sum_k.SP + 1 - sum_k.SO # up and down are equivalent if SP = 0 # ----------------- set-up input from DFT ----------------------- - if code in ('wien2k'): + if code in ('wien2k', 'elk'): # Check if wien converter was called and read transport subgroup form # hdf file if mpi.is_master_node(): @@ -361,7 +366,6 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}): raise IOError("transport_distribution: n_symmetries missing. Check if case.outputs file is present and call convert_misc_input() or convert_dft_input().") sum_k = read_transport_input_from_hdf(sum_k) - _, cell_volume = cellvolume(sum_k.lattice_type, sum_k.lattice_constants, sum_k.lattice_angles) elif code in ('wannier90'): required_entries = ['seedname', 'nk_optics'] @@ -380,7 +384,7 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}): assert all(isinstance(name, bool) for name in [calc_velocity, calc_inverse_mass]), f'Parameter {calc_velocity} or {calc_inverse_mass} not bool!' # recompute sum_k instances on denser grid - sum_k, cell_volume, _ = recompute_w90_input_on_different_mesh(sum_k, w90_params['seedname'], nk_optics=w90_params['nk_optics'], pathname=pathname, + sum_k, _ = recompute_w90_input_on_different_mesh(sum_k, w90_params['seedname'], nk_optics=w90_params['nk_optics'], pathname=pathname, calc_velocity=calc_velocity, calc_inverse_mass=calc_inverse_mass) # k-dependent-projections. @@ -388,10 +392,10 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}): # k_dep_projection is nowhere used # assert sum_k.k_dep_projection == 0, "transport_distribution: k dependent projection is not implemented!" - return sum_k, cell_volume + return sum_k # Uses .data of only GfReFreq objects. -def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0, code='wien2k'): +def transport_distribution(sum_k, beta, directions=['xx'], energy_window=None, Om_mesh=[0.0], with_Sigma=False, n_om=None, broadening=0.0, code='wien2k'): r""" Calculates the transport distribution @@ -406,8 +410,6 @@ def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_w triqs SumkDFT object beta : double Inverse temperature :math:`\beta`. - cell_volume : double - primitive unit cell volume directions : list of string, optional :math:`\alpha\beta` e.g.: ['xx','yy','zz','xy','xz','yz']. energy_window : list of double, optional @@ -495,8 +497,8 @@ def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_w assert broadening != 0.0 and broadening is not None, "transport_distribution: Broadening necessary to calculate transport distribution!" omega = numpy.linspace( energy_window[0] - max(Om_mesh), energy_window[1] + max(Om_mesh), n_om) - mesh = [energy_window[0] - - max(Om_mesh), energy_window[1] + max(Om_mesh), n_om] + mesh = MeshReFreq(energy_window[0] - + max(Om_mesh), energy_window[1] + max(Om_mesh), n_om) mu = 0.0 dir_to_int = {'x': 0, 'y': 1, 'z': 2} @@ -569,7 +571,7 @@ def transport_distribution(sum_k, beta, cell_volume, directions=['xx'], energy_w for direction in directions: - Gamma_w[direction] = (mpi.all_reduce(mpi.world, Gamma_w[direction], lambda x, y: x + y) / cell_volume / sum_k.n_symmetries) + Gamma_w[direction] = (mpi.all_reduce(mpi.world, Gamma_w[direction], lambda x, y: x + y) / sum_k.cell_vol / sum_k.n_symmetries) return Gamma_w, omega, temp_Om_mesh