From f24913b8a735f3f3306e3369e077b79aff4ee0a2 Mon Sep 17 00:00:00 2001 From: Priyanka Seth Date: Thu, 11 Dec 2014 12:11:48 +0100 Subject: [PATCH] [transport] more minor changes --- python/converters/wien2k_converter.py | 32 +++++++++------ python/sumk_dft_tools.py | 57 +++++++++++---------------- 2 files changed, 41 insertions(+), 48 deletions(-) diff --git a/python/converters/wien2k_converter.py b/python/converters/wien2k_converter.py index 4184caec..a3c691b5 100644 --- a/python/converters/wien2k_converter.py +++ b/python/converters/wien2k_converter.py @@ -368,7 +368,7 @@ class Wien2kConverter(ConverterTools): del ar - def convert_transport_input(self, spinbl=['']): + def convert_transport_input(self): """ Reads the input files necessary for transport calculations and stores the data in the HDFfile @@ -388,19 +388,26 @@ class Wien2kConverter(ConverterTools): # Read relevant data from .pmat file ############################################ - vk = [[] for sp in spinbl] - kp = [[] for sp in spinbl] + if (SP == 0 or SO == 1): + files = [self.vel_file] + elif SP == 1: + files = [self.vel_file+'up', self.vel_file+'dn'] + else: # SO and SP can't both be 1 + assert 0, "convert_transport_input: reading velocity file error! Check SP and SO!" + + vk = [[] for f in files] + kp = [[] for f in files] bandwin_opt = [] - for isp, sp in enumerate(spinbl): + for isp, f in enumerate(files): bandwin_opt_isp = [] - if not (os.path.exists(self.vel_file + sp)) : raise IOError, "File %s does not exist" %self.vel_file+sp - print "Reading input from %s..."%self.vel_file+sp + if not os.path.exists(f) : raise IOError, "File %s does not exist" %f + print "Reading input from %s..."%f - with open(self.vel_file + sp) as f: + with open(f) as R: while 1: try: - s = f.readline() + s = R.readline() if (s == ''): break except: @@ -410,12 +417,12 @@ class Wien2kConverter(ConverterTools): bandwin_opt_isp.append((nu1,nu2)) dim = nu2 - nu1 +1 v_xyz = numpy.zeros((dim,dim,3), dtype = complex) - temp = f.readline().strip().split() + temp = R.readline().strip().split() kp[isp].append(numpy.array([float(t) for t in temp[0:3]])) for nu_i in xrange(dim): for nu_j in xrange(nu_i, dim): for i in xrange(3): - s = f.readline().strip("\n ()").split(',') + s = R.readline().strip("\n ()").split(',') v_xyz[nu_i][nu_j][i] = float(s[0]) + float(s[1])*1j if (nu_i != nu_j): v_xyz[nu_j][nu_i][i] = v_xyz[nu_i][nu_j][i].conjugate() @@ -438,7 +445,6 @@ class Wien2kConverter(ConverterTools): try: f.readline() #title temp = f.readline() #lattice - #latticetype = temp[0:10].split()[0] latticetype = temp.split()[0] f.readline() temp = f.readline().strip().split() # lattice constants @@ -525,8 +531,8 @@ class Wien2kConverter(ConverterTools): ar = HDFArchive(self.hdf_file, 'a') if not (self.transp_subgrp in ar): ar.create_group(self.transp_subgrp) # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten!!! - things_to_save = ['bandwin_opt', 'kp', 'vk', 'latticetype', 'latticeconstants', 'latticeangles', 'nsymm', 'symmcartesian', - 'taucartesian', 'bandwin'] + things_to_save = ['bandwin', 'bandwin_opt', 'kp', 'vk', 'latticetype', 'latticeconstants', 'latticeangles', 'nsymm', 'symmcartesian', + 'taucartesian'] for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it] del ar diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py index c6a7ab1f..087a9bb9 100644 --- a/python/sumk_dft_tools.py +++ b/python/sumk_dft_tools.py @@ -610,10 +610,10 @@ class SumkDFTTools(SumkDFT): ik = 0 - bln = self.spin_block_names[self.SO] + spn = self.spin_block_names[self.SO] ntoi = self.spin_names_to_ind[self.SO] - S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window=(self.omega[0], self.omega[-1]), n_points = n_om)) + G_w = BlockGf(name_block_generator=[(spn[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window=(self.omega[0], self.omega[-1]), n_points = n_om)) for isp in range(n_inequiv_spin_blocks) ], make_copies=False) mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(n_inequiv_spin_blocks)] # construct mupat Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_) for isp in range(n_inequiv_spin_blocks)] @@ -623,7 +623,7 @@ class SumkDFTTools(SumkDFT): unchangedsize = all([ self.n_orbitals[ik][isp] == mupat[isp].shape[0] for isp in range(n_inequiv_spin_blocks)]) if not unchangedsize: # recontruct green functions. - S = BlockGf(name_block_generator=[(bln[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om)) + G_w = BlockGf(name_block_generator=[(spn[isp], GfReFreq(indices=range(self.n_orbitals[ik][isp]), window = (self.omega[0], self.omega[-1]), n_points = n_om)) for isp in range(n_inequiv_spin_blocks) ], make_copies=False) # construct mupat mupat = [numpy.identity(self.n_orbitals[ik][isp], numpy.complex_) * mu for isp in range(n_inequiv_spin_blocks)] @@ -631,32 +631,30 @@ class SumkDFTTools(SumkDFT): Annkw = [numpy.zeros((self.n_orbitals[ik][isp], self.n_orbitals[ik][isp], n_om), dtype=numpy.complex_) for isp in range(n_inequiv_spin_blocks)] # get lattice green function - S << 1*Omega + 1j*broadening + G_w << 1*Omega + 1j*broadening MS = copy.deepcopy(mupat) for ibl in range(n_inequiv_spin_blocks): - ind = ntoi[bln[ibl]] + ind = ntoi[spn[ibl]] n_orb = self.n_orbitals[ik][ibl] MS[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb].real - mupat[ibl] - S -= MS + G_w -= MS if (with_Sigma == True): - tmp = S.copy() # init temporary storage + tmp = G_w.copy() # init temporary storage # form self energy from impurity self energy and double counting term. sigma_minus_dc = self.add_dc(iw_or_w="w") # substract self energy - for icrsh in xrange(self.n_corr_shells): + for icrsh in range(self.n_corr_shells): for sig, gf in tmp: tmp[sig] << self.upfold(ik, icrsh, sig, sigma_minus_dc[icrsh][sig], gf) - S -= tmp + G_w -= tmp - S.invert() + G_w.invert() for isp in range(n_inequiv_spin_blocks): - Annkw[isp].real = -copy.deepcopy(S[self.spin_block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi + Annkw[isp].real = -copy.deepcopy(G_w[self.spin_block_names[self.SO][isp]].data.swapaxes(0,1).swapaxes(1,2)).imag / numpy.pi for isp in range(n_inequiv_spin_blocks): - if(ik%100==0): - print "ik,isp", ik, isp kvel = velocities[isp][ik] Pwtem = numpy.zeros((len(dir_list), len(Om_mesh_ex), n_om), dtype=numpy.float_) @@ -671,8 +669,8 @@ class SumkDFTTools(SumkDFT): for Rmat in self.symmcartesian: # get new velocity. Rkvel = copy.deepcopy(kvel) - for vnb1 in xrange(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1): - for vnb2 in xrange(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1): + for vnb1 in range(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1): + for vnb2 in range(self.bandwin_opt[isp][ik, 1] - self.bandwin_opt[isp][ik, 0] + 1): Rkvel[vnb1][vnb2][:] = numpy.dot(Rmat, Rkvel[vnb1][vnb2][:]) ipw = 0 for (ir, ic) in dir_list: @@ -722,18 +720,18 @@ class SumkDFTTools(SumkDFT): # if Om_meshr contains multiple entries with w=0, A1 is overwritten! q_0 = True A1 = numpy.zeros((n_direction,), dtype=numpy.float_) - for im in xrange(n_direction): + for idir in range(n_direction): for iw in xrange(n_w): - A0[im, iq] += beta * self.Pw_optic[im, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) - A1[im] += beta * self.Pw_optic[im, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw]) - self.seebeck[im] = -A1[im] / A0[im, iq] - print "A0", A0[im, iq] *d_omega/beta - print "A1", A1[im] *d_omega/beta + A0[idir, iq] += beta * self.Pw_optic[idir, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) + A1[idir] += beta * self.Pw_optic[idir, iq, iw] * self.fermi_dis(omegaT[iw]) * self.fermi_dis(-omegaT[iw]) * numpy.float(omegaT[iw]) + self.seebeck[idir] = -A1[idir] / A0[idir, iq] + print "A0", A0[idir, iq] *d_omega/beta + print "A1", A1[idir] *d_omega/beta # treat q ~= 0, calculate optical conductivity else: - for im in xrange(n_direction): + for idir in range(n_direction): for iw in xrange(n_w): - A0[im, iq] += self.Pw_optic[im, iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_meshr[iq] * beta)) / self.Om_meshr[iq] + A0[idir, iq] += self.Pw_optic[idir, iq, iw] * (self.fermi_dis(omegaT[iw]) - self.fermi_dis(omegaT[iw] + self.Om_meshr[iq] * beta)) / self.Om_meshr[iq] A0 *= d_omega #cond = beta * self.tdintegral(beta, 0)[index] @@ -743,7 +741,7 @@ class SumkDFTTools(SumkDFT): self.seebeck *= 86.17 # print - for im in xrange(n_direction): + for im in range(n_direction): for iq in xrange(n_q): print "Conductivity in direction %s for Om_mesh %d %.4f x 10^4 Ohm^-1 cm^-1" % (self.dir_list[im], iq, self.optic_cond[im, iq]) print "Resistivity in direction %s for Om_mesh %d %.4f x 10^-4 Ohm cm" % (self.dir_list[im], iq, 1.0 / self.optic_cond[im, iq]) @@ -751,17 +749,6 @@ class SumkDFTTools(SumkDFT): print "Seebeck in direction %s for q = 0 %.4f x 10^(-6) V/K" % (self.dir_list[im], self.seebeck[im]) - # ar = HDFArchive(self.hdf_file, 'a') - # if not (res_subgrp in ar): ar.create_group(res_subgrp) - # things_to_save = ['seebeck', 'optic_cond'] - # for it in things_to_save: ar[res_subgrp][it] = locals()[it] - # ar[res_subgrp]['seebeck'] = seebeck - # ar[res_subgrp]['optic_cond'] = optic_cond - # del ar - - # return optic_cond, seebeck - - def fermi_dis(self, x): return 1.0/(numpy.exp(x)+1)