From b6e33ecc2316e5d4f06ab7935aba7a3fae9aa325 Mon Sep 17 00:00:00 2001 From: Manuel Zingl Date: Tue, 28 Jul 2015 17:21:13 +0200 Subject: [PATCH] Add more integrators for the transport integral It is now possible to use trapz, simps and quadl (with cubic spline) to perform the omega integration needed in the transport code. --- python/sumk_dft_tools.py | 76 ++++++++++++++++++++++++++-------------- 1 file changed, 49 insertions(+), 27 deletions(-) diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py index 93f014e5..41870c30 100644 --- a/python/sumk_dft_tools.py +++ b/python/sumk_dft_tools.py @@ -18,13 +18,15 @@ # TRIQS. If not, see . # ################################################################################ - +import sys from types import * import numpy from pytriqs.gf.local import * import pytriqs.utility.mpi as mpi from symmetry import * from sumk_dft import SumkDFT +from scipy.integrate import * +from scipy.interpolate import * class SumkDFTTools(SumkDFT): """ @@ -613,8 +615,8 @@ class SumkDFTTools(SumkDFT): 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): + 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") @@ -718,8 +720,8 @@ class SumkDFTTools(SumkDFT): self.Gamma_w[direction] = (mpi.all_reduce(mpi.world, self.Gamma_w[direction], lambda x, y : x + y) / self.cellvolume(self.lattice_type, self.lattice_constants, self.lattice_angles)[1] / self.n_symmetries) - - def transport_coefficient(self, direction, iq, n, beta): + + 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. @@ -734,6 +736,9 @@ class SumkDFTTools(SumkDFT): 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 ------- @@ -744,24 +749,39 @@ class SumkDFTTools(SumkDFT): if not (mpi.is_master_node()): return assert hasattr(self,'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!" - A = 0.0 - omegaT = self.omega * beta - d_omega = self.omega[1] - self.omega[0] - if (self.Om_mesh[iq] == 0.0): - for iw in xrange(self.Gamma_w[direction].shape[1]): - A += self.Gamma_w[direction][iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw])**n * d_omega - elif (n == 0.0): - for iw in xrange(self.Gamma_w[direction].shape[1]): - A += (self.Gamma_w[direction][iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_mesh[iq] * beta)) - / (self.Om_mesh[iq] * beta) * d_omega) - else: - A = numpy.nan - A = A * numpy.pi * (2.0-self.SP) + 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 xrange(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): + + def conductivity_and_seebeck(self, beta, method=None): r""" Calculates the Seebeck coefficient and the optical conductivity by calling :meth:`transport_coefficient `. @@ -794,8 +814,8 @@ class SumkDFTTools(SumkDFT): for direction in self.directions: for iq in xrange(n_q): - A0[direction][iq] = self.transport_coefficient(direction, iq=iq, n=0, beta=beta) - A1[direction][iq] = self.transport_coefficient(direction, iq=iq, n=1, beta=beta) + 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) 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]) if ~numpy.isnan(A1[direction][iq]): @@ -810,7 +830,7 @@ class SumkDFTTools(SumkDFT): return self.optic_cond, self.seebeck - def fermi_dis(self, x): + def fermi_dis(self,w,beta): r""" Fermi distribution. @@ -819,11 +839,13 @@ class SumkDFTTools(SumkDFT): Parameters ---------- - x : double - Inverse temperature times frequency :math:`\beta\omega`. - + w : double + frequency + beta : double + inverse temperature + Returns ------- f : double """ - return 1.0/(numpy.exp(x)+1) + return 1.0/(numpy.exp(w*beta)+1)