From 0aeb54296c84a2c6d5ad9c7facf594ed4a5bbc48 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Wed, 4 Oct 2023 13:55:56 +0200 Subject: [PATCH 1/8] Update sumk_dft_transport.py inverse_mass as second derivative of Wannier Hamiltonian --- python/triqs_dft_tools/sumk_dft_transport.py | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index 8c5d9e3b..7cd3556c 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -336,13 +336,18 @@ def _commutator(A, B): velocities_k = (Hw_alpha + 1j * c_Hw_Aw_alpha) / HARTREETOEV / BOHRTOANG if calc_inverse_mass: - V_dot_D = numpy.einsum('kmnab, knoab -> kmoab', dataK.Xbar('Ham', 1) - [:, :, :, :, None], dataK.D_H[:, :, :, None, :]) - V_dot_D_dagger = V_dot_D.conj().transpose(0, 2, 1, 3, 4) - V_curly = numpy.einsum('knnab -> knab', V_dot_D + V_dot_D_dagger) - del2E_H_diag = numpy.einsum('knnab->knab', dataK.Xbar('Ham', 2)).real - inverse_mass = del2E_H_diag + V_curly - + # in the band basis + # ToDo: change units of inverse_mass consistently + if oc_basis == 'h': + inverse_mass = dataK.Xbar('Ham', 2) + # in the orbital basis + elif oc_basis == 'w': + Hw_alphabeta_R = dataK.Ham_R.copy() + for i in range(2): + shape_cR = numpy.shape(dataK.cRvec_wcc) + Hw_alphabeta_R = 1j * Hw_alphabeta_R.reshape((Hw_alphabeta_R.shape) + (1, )) * dataK.cRvec_wcc.reshape( + (shape_cR[0], shape_cR[1], dataK.system.nRvec) + (1, ) * len(Hw_alphabeta_R.shape[3:]) + (3, )) + inverse_mass = dataK.fft_R_to_k(Hw_alphabeta_R, hermitean=False)[dataK.select_K] # read in rest from dataK cell_volume = dataK.cell_volume / BOHRTOANG ** 3 kpts = dataK.kpoints_all From dc45b601c10f25b0bca2e729e2875ce74c793983 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Wed, 4 Oct 2023 14:12:51 +0200 Subject: [PATCH 2/8] Update sumk_dft_transport.py Add raman_vertex function --- python/triqs_dft_tools/sumk_dft_transport.py | 64 +++++++++++++++++++- 1 file changed, 63 insertions(+), 1 deletion(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index 7cd3556c..e1e7ffc2 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -32,7 +32,7 @@ import os.path __all__ = ['transport_distribution', 'conductivity_and_seebeck', 'write_output_to_hdf', - 'init_spectroscopy', 'transport_function'] + 'init_spectroscopy', 'transport_function', 'raman_vertex'] # ----------------- helper functions ----------------------- @@ -401,6 +401,68 @@ def _commutator(A, B): # ----------------- transport ----------------------- +def raman_vertex(sumk,ik,direction,code,isp, options=None): + r""" + Reads all necessary quantities for transport calculations from transport subgroup of the hdf5 archive. + Performs checks on input. Uses interpolation if code=wannier90. + + Parameters + ---------- + sum_k : sum_k object + triqs SumkDFT object + ik : integer + the ik index in which will be calculate the raman vertex + direction : string + direction to calculate the Raman vertex. + code : string + DFT code from which velocities are being read. Options: 'wien2k', 'wannier90' + isp : integer + spin index + options : dictionary, optional + additional keywords necessary in case a new direction wants to be implemented + + Returns + ------- + raman_vertex : 2D array + numpy array size [n_orb, n_orb] of complex type + """ + if code in ('wien2k'): + assert 0, 'Raman for wien2k not yet implemented' #ToDo + dir_names=['xx','yy','zz','B2g','B1g','A1g','Eg'] + dir_array=[ [[1.,0.,0.],[0.,0.,0.],[0.,0.,0.]], + [[0.,0.,0.],[0.,1.,0.],[0.,0.,0.]], + [[0.,0.,0.],[0.,0.,0.],[0.,0.,1.]], + [[0.,1.,0.],[0.,0.,0.],[0.,0.,0.]], + [[0.5,0.,0.],[0.,-0.5,0.],[0.,0.,0.]], + [[0.,0.,0.],[0.,0.,0.],[0.,0.,1.]], + [[0.,0.,1.],[0.,0.,0.],[0.,0.,0.]] ] + # Load custom directions + if "custom_dir" in options: + assert isinstance(options["custom_dir"],dict), "raman_vertex: in options, custom_dir must be a dictionary" + for dire in options["custom_dir"]: + assert numpy.shape(numpy.array(options["custom_dir"][dire]))==(3,3), "raman_vertex: custom_dir must have shape 3x3 (a numpy array or a nested list)" + if dire in dir_names: + if direction==dire and ik==0: #reports only when ik==0 + mpi.report("Warning: the direction %s was already loaded and will be replace with the one provided in custom_dir"%dire) + idir = dir_names.index(dire) + dir_array[idir]=options["custom_dir"][dire] + else: + dir_names.append(dire) + dir_array.append(options["custom_dir"][dire]) + + dir_array=[numpy.array(el,dtype=numpy.float_) for el in dir_array] # convert list to numpy array + + if direction in dir_names: + idir = dir_names.index(direction) + else: + assert 0,'raman_vertex: direction %s is not supported by default. Try to add it by raman_vertex options custom_dir.'%direction + n_bands = sumk.n_orbitals[ik][isp] + ram_vert = numpy.zeros( (n_bands, n_bands), dtype=numpy.complex_) + for i in range(n_bands): + for j in range(n_bands): + ram_vert[i,j]=numpy.dot(dir_array[idir],sumk.inverse_mass[ik,i,j,:,:]).trace() + + return ram_vert def init_spectroscopy(sum_k, code='wien2k', w90_params={}): r""" From a2f5d986cf04d69b18ce38e1a5e5b0cced1bb0a0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Wed, 4 Oct 2023 14:36:21 +0200 Subject: [PATCH 3/8] Update sumk_dft_transport.py Include Raman in transport_distribution --- python/triqs_dft_tools/sumk_dft_transport.py | 69 +++++++++++++------- 1 file changed, 46 insertions(+), 23 deletions(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index e1e7ffc2..8cd91a55 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -556,7 +556,7 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}): # Uses .data of only GfReFreq objects. -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'): +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', mode='optics', raman_options={}): r""" Calculates the transport distribution @@ -589,6 +589,10 @@ def transport_distribution(sum_k, beta, directions=['xx'], energy_window=None, O 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' + mode : string + Choose between optical ('optics') or Raman ('raman') transport distribution. + raman_options : dictionary + additional keywords necessary in case mode == 'raman'. Depending on the situation, the allow keys could be 'custom_dir'. Returns ------- @@ -713,28 +717,47 @@ def glist(): return [GfReFreq(target_shape=(block_dim, block_dim), window=(omega v_i = slice(b_min - sum_k.band_window_optics[isp][ ik, 0], b_max - sum_k.band_window_optics[isp][ik, 0] + 1) - # loop over all symmetries - for R in sum_k.rot_symmetries: - # get transformed velocity under symmetry R - if code in ('wien2k'): - vel_R = copy.deepcopy(sum_k.velocities_k[isp][ik]) - elif code in ('wannier90'): - vel_R = copy.deepcopy(sum_k.velocities_k[ik]) - for nu1 in range(sum_k.band_window_optics[isp][ik, 1] - sum_k.band_window_optics[isp][ik, 0] + 1): - for nu2 in range(sum_k.band_window_optics[isp][ik, 1] - sum_k.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 directions: - for iw in range(n_om): - for iq in range(len(temp_Om_mesh)): - if (iw + iOm_mesh[iq] >= n_om or omega[iw] < -temp_Om_mesh[iq] + energy_window[0] or omega[iw] > temp_Om_mesh[iq] + energy_window[1]): - continue - - 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 * sum_k.bz_weights[ik]) + if mode in ('optics'): + # loop over all symmetries + for R in sum_k.rot_symmetries: + # get transformed velocity under symmetry R + if code in ('wien2k'): + vel_R = copy.deepcopy(sum_k.velocities_k[isp][ik]) + elif code in ('wannier90'): + vel_R = copy.deepcopy(sum_k.velocities_k[ik]) + for nu1 in range(sum_k.band_window_optics[isp][ik, 1] - sum_k.band_window_optics[isp][ik, 0] + 1): + for nu2 in range(sum_k.band_window_optics[isp][ik, 1] - sum_k.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 directions: + for iw in range(n_om): + for iq in range(len(temp_Om_mesh)): + if (iw + iOm_mesh[iq] >= n_om or omega[iw] < -temp_Om_mesh[iq] + energy_window[0] or omega[iw] > temp_Om_mesh[iq] + energy_window[1]): + continue + + 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 * sum_k.bz_weights[ik]) + elif mode in ('raman'): + if code in ('wannier90'): + assert hasattr(sum_k,"inverse_mass"), 'inverse_mass not available in sum_k. Set calc_inverse_mass=True in w90_params.' + elif code in ('wien2k'): + assert 0, 'Raman for wien2k not yet implemented' #ToDo + # loop over all symmetries + for R in sum_k.rot_symmetries: + for direction in directions: + # calculate the raman vertex for each direction + vert = raman_vertex(sum_k, ik, direction, code, isp, raman_options) + + for iw in range(n_om): + for iq in range(len(Om_mesh)): + if (iw + iOm_mesh[iq] >= n_om or omega[iw] < -Om_mesh[iq] + energy_window[0] or omega[iw] > Om_mesh[iq] + energy_window[1]): + continue + + Gamma_w[direction][iq, iw] += (numpy.dot(numpy.dot(numpy.dot(vert[v_i, v_i], A_kw[isp][A_i, A_i, int(iw + iOm_mesh[iq])]), + vert[v_i, v_i]), A_kw[isp][A_i, A_i, iw]).trace().real * sum_k.bz_weights[ik]) for direction in directions: Gamma_w[direction] = (mpi.all_reduce(Gamma_w[direction]) / sum_k.cell_vol / sum_k.n_symmetries) From 297e4ecd9891d8264b250005903a5a41471d5019 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Wed, 4 Oct 2023 14:50:04 +0200 Subject: [PATCH 4/8] Update sumk_dft_transport.py Implement Raman in conductivity_and_seebeck function. --- python/triqs_dft_tools/sumk_dft_transport.py | 102 +++++++++++-------- 1 file changed, 59 insertions(+), 43 deletions(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index 8cd91a55..eae3783e 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -899,7 +899,7 @@ def transport_coefficient(Gamma_w, omega, Om_mesh, spin_polarization, direction, return A -def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, method=None): +def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, method=None, mode='optics'): r""" Calculates the Seebeck coefficient and the optical conductivity by calling :meth:`transport_coefficient `. @@ -923,7 +923,9 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth method : string Integration method: cubic spline and scipy.integrate.quad ('quad'), simpson rule ('simpson'), trapezoidal rule ('trapz'), rectangular integration (otherwise) Note that the sampling points of the the self-energy are used! - + mode : string + Choose between optical conductivity/seebeck/Kappa ('optics') or Raman conductivity ('raman') + Returns ------- optic_cond : dictionary of double vectors @@ -945,45 +947,59 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth # initialization A0 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} - A1 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} - A2 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} - optic_cond = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} - seebeck = {direction: numpy.nan for direction in directions} - kappa = {direction: numpy.nan for direction in directions} + if mode in ('optics'): + A1 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} + A2 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} + optic_cond = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} + seebeck = {direction: numpy.nan for direction in directions} + kappa = {direction: numpy.nan for direction in directions} - for direction in directions: - for iq in range(n_q): - A0[direction][iq] = transport_coefficient( - Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=0, beta=beta, method=method) - A1[direction][iq] = transport_coefficient( - Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=1, beta=beta, method=method) - A2[direction][iq] = transport_coefficient( - Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=2, beta=beta, method=method) - print("A_0 in direction %s for Omega = %.2f %e a.u." % - (direction, Om_mesh[iq], A0[direction][iq])) - print("A_1 in direction %s for Omega = %.2f %e a.u." % - (direction, Om_mesh[iq], A1[direction][iq])) - print("A_2 in direction %s for Omega = %.2f %e a.u." % - (direction, 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 - seebeck[direction] = - A1[direction][iq] / A0[direction][iq] * 86.17 - kappa[direction] = A2[direction][iq] - \ - A1[direction][iq]*A1[direction][iq]/A0[direction][iq] - kappa[direction] *= 293178.0 - - # factor for optical conductivity: hbar * velocity_Hartree_to_SI * volume_Hartree_to_SI * m_to_cm * 10^-4 final unit - convert_to_SI = cst.hbar * (cst.c * cst.fine_structure) ** 2 * \ - (1/cst.physical_constants['Bohr radius'][0]) ** 3 * 1e-6 - optic_cond[direction] = beta * convert_to_SI * A0[direction] - for iq in range(n_q): - print("Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % - (direction, Om_mesh[iq], 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, seebeck[direction])) - print("kappa in direction %s for Omega = 0.00 %f W/(m * K)" % - (direction, kappa[direction])) - - return optic_cond, seebeck, kappa + for direction in directions: + for iq in range(n_q): + A0[direction][iq] = transport_coefficient( + Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=0, beta=beta, method=method) + A1[direction][iq] = transport_coefficient( + Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=1, beta=beta, method=method) + A2[direction][iq] = transport_coefficient( + Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=2, beta=beta, method=method) + print("A_0 in direction %s for Omega = %.2f %e a.u." % + (direction, Om_mesh[iq], A0[direction][iq])) + print("A_1 in direction %s for Omega = %.2f %e a.u." % + (direction, Om_mesh[iq], A1[direction][iq])) + print("A_2 in direction %s for Omega = %.2f %e a.u." % + (direction, 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 + seebeck[direction] = - A1[direction][iq] / A0[direction][iq] * 86.17 + kappa[direction] = A2[direction][iq] - \ + A1[direction][iq]*A1[direction][iq]/A0[direction][iq] + kappa[direction] *= 293178.0 + + # factor for optical conductivity: hbar * velocity_Hartree_to_SI * volume_Hartree_to_SI * m_to_cm * 10^-4 final unit + convert_to_SI = cst.hbar * (cst.c * cst.fine_structure) ** 2 * \ + (1/cst.physical_constants['Bohr radius'][0]) ** 3 * 1e-6 + optic_cond[direction] = beta * convert_to_SI * A0[direction] + for iq in range(n_q): + print("Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % + (direction, Om_mesh[iq], 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, seebeck[direction])) + print("kappa in direction %s for Omega = 0.00 %f W/(m * K)" % + (direction, kappa[direction])) + + return optic_cond, seebeck, kappa + elif mode in ('raman'): + # ToDo: correct units + raman_cond = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} + + for direction in directions: + for iq in range(n_q): + A0[direction][iq] = transport_coefficient(Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=0, beta=beta, method=method) + print("A_0 in direction %s for Omega = %.2f %e a.u." % (direction, Om_mesh[iq], A0[direction][iq])) + raman_cond[direction] = beta * A0[direction] * 10700.0 / numpy.pi + for iq in range(n_q): + print("Raman conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % (direction, Om_mesh[iq], raman_cond[direction][iq])) + + return raman_cond From 6dde56d1c5d19a258d54f6a0c562af0c2a694f3e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Thu, 5 Oct 2023 11:15:51 +0200 Subject: [PATCH 5/8] Update sumk_dft_transport.py Add optic_kappa to conductivity_and_seebeck --- python/triqs_dft_tools/sumk_dft_transport.py | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index eae3783e..8e75d462 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -899,7 +899,7 @@ def transport_coefficient(Gamma_w, omega, Om_mesh, spin_polarization, direction, return A -def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, method=None, mode='optics'): +def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, method=None, mode='optics', optic_kappa=False): r""" Calculates the Seebeck coefficient and the optical conductivity by calling :meth:`transport_coefficient `. @@ -925,6 +925,8 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth Note that the sampling points of the the self-energy are used! mode : string Choose between optical conductivity/seebeck/Kappa ('optics') or Raman conductivity ('raman') + optic_kappa : bool + If calculates $\kappa(\omega)$ or not. Only if mode is 'optics'. Returns ------- @@ -935,7 +937,8 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth 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 + thermal conductivity in each direction. If zero is not present in Om_mesh the thermal conductivity is set to NaN. + If optic_kappa is True, then returns the kappa for all frequencies. """ mpi.report('Computing optical conductivity and kinetic coefficients...') @@ -953,6 +956,7 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth optic_cond = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} seebeck = {direction: numpy.nan for direction in directions} kappa = {direction: numpy.nan for direction in directions} + if optic_kappa: optic_kappa = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} for direction in directions: for iq in range(n_q): @@ -980,16 +984,22 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth convert_to_SI = cst.hbar * (cst.c * cst.fine_structure) ** 2 * \ (1/cst.physical_constants['Bohr radius'][0]) ** 3 * 1e-6 optic_cond[direction] = beta * convert_to_SI * A0[direction] + if optic_kappa: optic_kappa[direction] = (A2[direction] - A1[direction]*A1[direction]/A0[direction])*293178.0 for iq in range(n_q): print("Conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % (direction, Om_mesh[iq], optic_cond[direction][iq])) + if optic_kappa: + print("kappa(Omega) in direction %s for Omega = %.2f %f W/(m * K)" % + (direction, Om_mesh[iq], optic_kappa[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, seebeck[direction])) print("kappa in direction %s for Omega = 0.00 %f W/(m * K)" % (direction, kappa[direction])) - - return optic_cond, seebeck, kappa + if optic_kappa: + return optic_cond, seebeck, optic_kappa + else: + return optic_cond, seebeck, kappa elif mode in ('raman'): # ToDo: correct units raman_cond = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} From 4711e1d4b00afd5dcb8f822aee67569a1bc737f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Thu, 5 Oct 2023 14:35:19 +0200 Subject: [PATCH 6/8] Update sumk_dft_transport.py Raman conductivity in atomic units. --- python/triqs_dft_tools/sumk_dft_transport.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index 8e75d462..48cffdc7 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -337,7 +337,6 @@ def _commutator(A, B): if calc_inverse_mass: # in the band basis - # ToDo: change units of inverse_mass consistently if oc_basis == 'h': inverse_mass = dataK.Xbar('Ham', 2) # in the orbital basis @@ -348,6 +347,7 @@ def _commutator(A, B): Hw_alphabeta_R = 1j * Hw_alphabeta_R.reshape((Hw_alphabeta_R.shape) + (1, )) * dataK.cRvec_wcc.reshape( (shape_cR[0], shape_cR[1], dataK.system.nRvec) + (1, ) * len(Hw_alphabeta_R.shape[3:]) + (3, )) inverse_mass = dataK.fft_R_to_k(Hw_alphabeta_R, hermitean=False)[dataK.select_K] + inverse_mass = inverse_mass / HARTREETOEV / BOHRTOANG**2 # read in rest from dataK cell_volume = dataK.cell_volume / BOHRTOANG ** 3 kpts = dataK.kpoints_all @@ -427,7 +427,7 @@ def raman_vertex(sumk,ik,direction,code,isp, options=None): numpy array size [n_orb, n_orb] of complex type """ if code in ('wien2k'): - assert 0, 'Raman for wien2k not yet implemented' #ToDo + assert 0, 'Raman for wien2k not yet implemented' dir_names=['xx','yy','zz','B2g','B1g','A1g','Eg'] dir_array=[ [[1.,0.,0.],[0.,0.,0.],[0.,0.,0.]], [[0.,0.,0.],[0.,1.,0.],[0.,0.,0.]], @@ -744,7 +744,7 @@ def glist(): return [GfReFreq(target_shape=(block_dim, block_dim), window=(omega if code in ('wannier90'): assert hasattr(sum_k,"inverse_mass"), 'inverse_mass not available in sum_k. Set calc_inverse_mass=True in w90_params.' elif code in ('wien2k'): - assert 0, 'Raman for wien2k not yet implemented' #ToDo + assert 0, 'Raman for wien2k not yet implemented' # loop over all symmetries for R in sum_k.rot_symmetries: for direction in directions: @@ -1001,15 +1001,15 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth else: return optic_cond, seebeck, kappa elif mode in ('raman'): - # ToDo: correct units + # TODO: express in SI units raman_cond = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} for direction in directions: for iq in range(n_q): A0[direction][iq] = transport_coefficient(Gamma_w, omega, Om_mesh, SP, direction, iq=iq, n=0, beta=beta, method=method) - print("A_0 in direction %s for Omega = %.2f %e a.u." % (direction, Om_mesh[iq], A0[direction][iq])) - raman_cond[direction] = beta * A0[direction] * 10700.0 / numpy.pi + print("A_0 in direction %s for Omega = %.2f %e atomic units" % (direction, Om_mesh[iq], A0[direction][iq])) + raman_cond[direction] = beta * A0[direction] for iq in range(n_q): - print("Raman conductivity in direction %s for Omega = %.2f %f x 10^4 Ohm^-1 cm^-1" % (direction, Om_mesh[iq], raman_cond[direction][iq])) + print("Raman conductivity in direction %s for Omega = %.4f %f atomic units" % (direction, Om_mesh[iq], raman_cond[direction][iq])) return raman_cond From dd7001a89a775b0b8daeaaa26b237b965b6ade61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Germ=C3=A1n=20Blesio?= <48755986+gblesio@users.noreply.github.com> Date: Mon, 9 Oct 2023 17:11:13 +0200 Subject: [PATCH 7/8] Update sumk_dft_transport.py Allow calc_inverse_mass for both basis --- python/triqs_dft_tools/sumk_dft_transport.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index 48cffdc7..16523a7c 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -536,7 +536,6 @@ def init_spectroscopy(sum_k, code='wien2k', w90_params={}): oc_select = 'both' # further checks for calc_inverse_mass if calc_inverse_mass: - assert oc_basis == 'h', '"calc_inverse_mass" only implemented for "oc_basis" == "h"' assert oc_select == 'both', '"oc_select" not implemented for "calc_inverse_mass"' # print some information mpi.report(f'{"Basis choice [h (Hamiltonian), w (Wannier)]:":<60s} {oc_basis}') From 5c87dec2e4c9708da8015cd8ee65715a4f6b6b7c Mon Sep 17 00:00:00 2001 From: phibeck Date: Fri, 16 Aug 2024 10:40:35 -0400 Subject: [PATCH 8/8] adding comments regarding PR #246 --- python/triqs_dft_tools/sumk_dft_transport.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/python/triqs_dft_tools/sumk_dft_transport.py b/python/triqs_dft_tools/sumk_dft_transport.py index 16523a7c..3fc96a08 100644 --- a/python/triqs_dft_tools/sumk_dft_transport.py +++ b/python/triqs_dft_tools/sumk_dft_transport.py @@ -336,17 +336,15 @@ def _commutator(A, B): velocities_k = (Hw_alpha + 1j * c_Hw_Aw_alpha) / HARTREETOEV / BOHRTOANG if calc_inverse_mass: + # approximating this as the first term in eq (28) in PRB 75, 195121 (2007) + # application published in PHYSICAL REVIEW RESEARCH 6, 023124 (2024) + # Procedure is the same as for calc_velocity, but using second derivative of Hamiltonian. # in the band basis if oc_basis == 'h': inverse_mass = dataK.Xbar('Ham', 2) # in the orbital basis elif oc_basis == 'w': - Hw_alphabeta_R = dataK.Ham_R.copy() - for i in range(2): - shape_cR = numpy.shape(dataK.cRvec_wcc) - Hw_alphabeta_R = 1j * Hw_alphabeta_R.reshape((Hw_alphabeta_R.shape) + (1, )) * dataK.cRvec_wcc.reshape( - (shape_cR[0], shape_cR[1], dataK.system.nRvec) + (1, ) * len(Hw_alphabeta_R.shape[3:]) + (3, )) - inverse_mass = dataK.fft_R_to_k(Hw_alphabeta_R, hermitean=False)[dataK.select_K] + inverse_mass = wb.data_K.Data_K_R._R_to_k_H(dataK, dataK.Ham_R, der=2, hermitian=True) inverse_mass = inverse_mass / HARTREETOEV / BOHRTOANG**2 # read in rest from dataK cell_volume = dataK.cell_volume / BOHRTOANG ** 3 @@ -949,6 +947,7 @@ def conductivity_and_seebeck(Gamma_w, omega, Om_mesh, SP, directions, beta, meth # initialization A0 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} + # The conductivity is proportional to A0, A1 and A2 are only needed for the Seebeck coefficient and thermal conductivity if mode in ('optics'): A1 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions} A2 = {direction: numpy.full((n_q,), numpy.nan) for direction in directions}