From 390e8564b7bcdd8eb3685143062871a10b3ad563 Mon Sep 17 00:00:00 2001 From: Priyanka Seth Date: Mon, 9 May 2016 10:19:56 +0200 Subject: [PATCH] Minor clean up, pep-ified to allow doc compilation to run smoothly --- python/__init__.py | 7 +- python/clear_h5_output.py | 13 +- python/converters/__init__.py | 8 +- python/converters/converter_tools.py | 56 +- python/converters/hk_converter.py | 212 +++--- python/converters/wannier90_converter.py | 104 ++- python/converters/wien2k_converter.py | 588 ++++++++------- python/sumk_dft.py | 865 +++++++++++++---------- python/sumk_dft_tools.py | 676 +++++++++++------- python/symmetry.py | 86 ++- python/trans_basis.py | 137 ++-- python/update_archive.py | 95 +-- 12 files changed, 1663 insertions(+), 1184 deletions(-) diff --git a/python/__init__.py b/python/__init__.py index d61d896e..b6c3ad2d 100644 --- a/python/__init__.py +++ b/python/__init__.py @@ -1,5 +1,5 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,11 +18,12 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## from sumk_dft import SumkDFT from symmetry import Symmetry from sumk_dft_tools import SumkDFTTools from converters import * -__all__=['SumkDFT','Symmetry','SumkDFTTools','Wien2kConverter','HkConverter'] +__all__ = ['SumkDFT', 'Symmetry', 'SumkDFTTools', + 'Wien2kConverter', 'HkConverter'] diff --git a/python/clear_h5_output.py b/python/clear_h5_output.py index 9a110cf0..a9135771 100644 --- a/python/clear_h5_output.py +++ b/python/clear_h5_output.py @@ -3,8 +3,8 @@ import sys import subprocess if len(sys.argv) < 2: - print "Usage: python clear_h5_output.py archive" - sys.exit() + print "Usage: python clear_h5_output.py archive" + sys.exit() print """ This script is to remove any SumkDFT generated output from the h5 archive @@ -13,13 +13,14 @@ and to restore it to the original post-converter state. filename = sys.argv[1] A = h5py.File(filename) -for group in ['dmft_output','user_data']: - if group in A: del(A[group]) +for group in ['dmft_output', 'user_data']: + if group in A: + del(A[group]) A.close() # Repack to reclaim disk space -retcode = subprocess.call(["h5repack","-i%s"%filename, "-otemphgfrt.h5"]) +retcode = subprocess.call(["h5repack", "-i%s" % filename, "-otemphgfrt.h5"]) if retcode != 0: print "h5repack failed!" else: - subprocess.call(["mv","-f","temphgfrt.h5","%s"%filename]) + subprocess.call(["mv", "-f", "temphgfrt.h5", "%s" % filename]) diff --git a/python/converters/__init__.py b/python/converters/__init__.py index dcbb6e08..f753f455 100644 --- a/python/converters/__init__.py +++ b/python/converters/__init__.py @@ -1,5 +1,5 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,12 +18,10 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## from wien2k_converter import Wien2kConverter from hk_converter import HkConverter from wannier90_converter import Wannier90Converter -__all__ =['Wien2kConverter','HkConverter','Wannier90Converter'] - - +__all__ = ['Wien2kConverter', 'HkConverter', 'Wannier90Converter'] diff --git a/python/converters/converter_tools.py b/python/converters/converter_tools.py index 84051b13..126ffd57 100644 --- a/python/converters/converter_tools.py +++ b/python/converters/converter_tools.py @@ -1,5 +1,5 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,16 +18,17 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## from pytriqs.cmake_info import hdf5_command_path import pytriqs.utility.mpi as mpi + class ConverterTools: def __init__(self): pass - def read_fortran_file(self,filename,to_replace): + def read_fortran_file(self, filename, to_replace): """ Returns a generator that yields all numbers in the Fortran file as float, with possible replacements. @@ -37,7 +38,7 @@ class ConverterTools: Name of Fortran-produced file. to_replace : dict of str:str Dictionary defining old_char:new_char. - + Yields ------ string @@ -46,11 +47,13 @@ class ConverterTools: """ import os.path import string - if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename - for line in open(filename,'r') : - for old,new in to_replace.iteritems(): line = line.replace(old,new) - for x in line.split(): yield string.atof(x) - + if not(os.path.exists(filename)): + raise IOError, "File %s does not exist." % filename + for line in open(filename, 'r'): + for old, new in to_replace.iteritems(): + line = line.replace(old, new) + for x in line.split(): + yield string.atof(x) def repack(self): """ @@ -65,17 +68,18 @@ class ConverterTools: import subprocess - if not (mpi.is_master_node()): return - mpi.report("Repacking the file %s"%self.hdf_file) + if not (mpi.is_master_node()): + return + mpi.report("Repacking the file %s" % self.hdf_file) - retcode = subprocess.call([hdf5_command_path+"/h5repack","-i%s"%self.hdf_file,"-otemphgfrt.h5"]) + retcode = subprocess.call( + [hdf5_command_path + "/h5repack", "-i%s" % self.hdf_file, "-otemphgfrt.h5"]) if retcode != 0: mpi.report("h5repack failed!") else: - subprocess.call(["mv","-f","temphgfrt.h5","%s"%self.hdf_file]) - + subprocess.call(["mv", "-f", "temphgfrt.h5", "%s" % self.hdf_file]) - def det_shell_equivalence(self,corr_shells): + def det_shell_equivalence(self, corr_shells): """ Determine the equivalence of correlated shells. @@ -83,7 +87,7 @@ class ConverterTools: ---------- corr_shells : list of dicts See documentation of necessary hdf5 elements. - + Returns ------- n_inequiv_shells : integer @@ -105,19 +109,19 @@ class ConverterTools: n_inequiv_shells = 1 if len(corr_shells) > 1: - inequiv_sort = [ corr_shells[0]['sort'] ] - inequiv_l = [ corr_shells[0]['l'] ] - for i in range(len(corr_shells)-1): + inequiv_sort = [corr_shells[0]['sort']] + inequiv_l = [corr_shells[0]['l']] + for i in range(len(corr_shells) - 1): is_equiv = False for j in range(n_inequiv_shells): - if (inequiv_sort[j]==corr_shells[i+1]['sort']) and (inequiv_l[j]==corr_shells[i+1]['l']): + if (inequiv_sort[j] == corr_shells[i + 1]['sort']) and (inequiv_l[j] == corr_shells[i + 1]['l']): is_equiv = True - corr_to_inequiv[i+1] = j - if is_equiv==False: - corr_to_inequiv[i+1] = n_inequiv_shells + corr_to_inequiv[i + 1] = j + if is_equiv == False: + corr_to_inequiv[i + 1] = n_inequiv_shells n_inequiv_shells += 1 - inequiv_sort.append( corr_shells[i+1]['sort'] ) - inequiv_l.append( corr_shells[i+1]['l'] ) - inequiv_to_corr.append( i+1 ) + inequiv_sort.append(corr_shells[i + 1]['sort']) + inequiv_l.append(corr_shells[i + 1]['l']) + inequiv_to_corr.append(i + 1) return n_inequiv_shells, corr_to_inequiv, inequiv_to_corr diff --git a/python/converters/hk_converter.py b/python/converters/hk_converter.py index dc372c65..26fe1656 100644 --- a/python/converters/hk_converter.py +++ b/python/converters/hk_converter.py @@ -1,5 +1,5 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,7 +18,7 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## from types import * import numpy @@ -27,12 +27,13 @@ import pytriqs.utility.mpi as mpi from math import sqrt from converter_tools import * + class HkConverter(ConverterTools): """ Conversion from general H(k) file to an hdf5 file that can be used as input for the SumKDFT class. """ - def __init__(self, filename, hdf_filename = None, dft_subgrp = 'dft_input', symmcorr_subgrp = 'dft_symmcorr_input', repacking = False): + def __init__(self, filename, hdf_filename=None, dft_subgrp='dft_input', symmcorr_subgrp='dft_symmcorr_input', repacking=False): """ Initialise the class. @@ -49,24 +50,25 @@ class HkConverter(ConverterTools): The group is actually empty; it is just included for compatibility. repacking : boolean, optional Does the hdf5 archive need to be repacked to save space? - + """ - assert type(filename)==StringType,"HkConverter: filename must be a filename." - if hdf_filename is None: hdf_filename = filename+'.h5' + assert type( + filename) == StringType, "HkConverter: filename must be a filename." + if hdf_filename is None: + hdf_filename = filename + '.h5' self.hdf_file = hdf_filename self.dft_file = filename self.dft_subgrp = dft_subgrp self.symmcorr_subgrp = symmcorr_subgrp - self.fortran_to_replace = {'D':'E', '(':' ', ')':' ', ',':' '} + self.fortran_to_replace = {'D': 'E', '(': ' ', ')': ' ', ',': ' '} # Checks if h5 file is there and repacks it if wanted: import os.path if (os.path.exists(self.hdf_file) and repacking): ConverterTools.repack(self) - - def convert_dft_input(self, first_real_part_matrix = True, only_upper_triangle = False, weights_in_file = False): + def convert_dft_input(self, first_real_part_matrix=True, only_upper_triangle=False, weights_in_file=False): """ Reads the appropriate files and stores the data for the dft_subgrp in the hdf5 archive. @@ -80,71 +82,97 @@ class HkConverter(ConverterTools): Are the k-point weights to be read in? """ - - # Read and write only on the master node - if not (mpi.is_master_node()): return - mpi.report("Reading input from %s..."%self.dft_file) - # R is a generator : each R.Next() will return the next number in the file - R = ConverterTools.read_fortran_file(self,self.dft_file,self.fortran_to_replace) + # Read and write only on the master node + if not (mpi.is_master_node()): + return + mpi.report("Reading input from %s..." % self.dft_file) + + # R is a generator : each R.Next() will return the next number in the + # file + R = ConverterTools.read_fortran_file( + self, self.dft_file, self.fortran_to_replace) try: - energy_unit = 1.0 # the energy conversion factor is 1.0, we assume eV in files - n_k = int(R.next()) # read the number of k points - k_dep_projection = 0 + # the energy conversion factor is 1.0, we assume eV in files + energy_unit = 1.0 + # read the number of k points + n_k = int(R.next()) + k_dep_projection = 0 SP = 0 # no spin-polarision - SO = 0 # no spin-orbit - charge_below = 0.0 # total charge below energy window is set to 0 - density_required = R.next() # density required, for setting the chemical potential + SO = 0 # no spin-orbit + # total charge below energy window is set to 0 + charge_below = 0.0 + # density required, for setting the chemical potential + density_required = R.next() symm_op = 0 # No symmetry groups for the k-sum - # the information on the non-correlated shells is needed for defining dimension of matrices: - n_shells = int(R.next()) # number of shells considered in the Wanniers - # corresponds to index R in formulas + # the information on the non-correlated shells is needed for + # defining dimension of matrices: + # number of shells considered in the Wanniers + n_shells = int(R.next()) + # corresponds to index R in formulas # now read the information about the shells (atom, sort, l, dim): shell_entries = ['atom', 'sort', 'l', 'dim'] - shells = [ {name: int(val) for name, val in zip(shell_entries, R)} for ish in range(n_shells) ] + shells = [{name: int(val) for name, val in zip( + shell_entries, R)} for ish in range(n_shells)] - n_corr_shells = int(R.next()) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, - # corresponds to index R in formulas - # now read the information about the shells (atom, sort, l, dim, SO flag, irep): + # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, + n_corr_shells = int(R.next()) + # corresponds to index R in formulas + # now read the information about the shells (atom, sort, l, dim, SO + # flag, irep): corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep'] - corr_shells = [ {name: int(val) for name, val in zip(corr_shell_entries, R)} for icrsh in range(n_corr_shells) ] + corr_shells = [{name: int(val) for name, val in zip( + corr_shell_entries, R)} for icrsh in range(n_corr_shells)] - # determine the number of inequivalent correlated shells and maps, needed for further reading - [n_inequiv_shells, corr_to_inequiv, inequiv_to_corr] = ConverterTools.det_shell_equivalence(self,corr_shells) + # determine the number of inequivalent correlated shells and maps, + # needed for further reading + [n_inequiv_shells, corr_to_inequiv, + inequiv_to_corr] = ConverterTools.det_shell_equivalence(self, corr_shells) use_rotations = 0 - rot_mat = [numpy.identity(corr_shells[icrsh]['dim'],numpy.complex_) for icrsh in range(n_corr_shells)] + rot_mat = [numpy.identity( + corr_shells[icrsh]['dim'], numpy.complex_) for icrsh in range(n_corr_shells)] rot_mat_time_inv = [0 for i in range(n_corr_shells)] - + # Representative representations are read from file n_reps = [1 for i in range(n_inequiv_shells)] dim_reps = [0 for i in range(n_inequiv_shells)] T = [] for ish in range(n_inequiv_shells): - n_reps[ish] = int(R.next()) # number of representatives ("subsets"), e.g. t2g and eg - dim_reps[ish] = [int(R.next()) for i in range(n_reps[ish])] # dimensions of the subsets - + # number of representatives ("subsets"), e.g. t2g and eg + n_reps[ish] = int(R.next()) + dim_reps[ish] = [int(R.next()) for i in range( + n_reps[ish])] # dimensions of the subsets + # The transformation matrix: - # is of dimension 2l+1, it is taken to be standard d (as in Wien2k) - ll = 2*corr_shells[inequiv_to_corr[ish]]['l']+1 + # is of dimension 2l+1, it is taken to be standard d (as in + # Wien2k) + ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1 lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1) - T.append(numpy.zeros([lmax,lmax],numpy.complex_)) - + T.append(numpy.zeros([lmax, lmax], numpy.complex_)) + T[ish] = numpy.array([[0.0, 0.0, 1.0, 0.0, 0.0], - [1.0/sqrt(2.0), 0.0, 0.0, 0.0, 1.0/sqrt(2.0)], - [-1.0/sqrt(2.0), 0.0, 0.0, 0.0, 1.0/sqrt(2.0)], - [0.0, 1.0/sqrt(2.0), 0.0, -1.0/sqrt(2.0), 0.0], - [0.0, 1.0/sqrt(2.0), 0.0, 1.0/sqrt(2.0), 0.0]]) + [1.0 / sqrt(2.0), 0.0, 0.0, + 0.0, 1.0 / sqrt(2.0)], + [-1.0 / sqrt(2.0), 0.0, 0.0, + 0.0, 1.0 / sqrt(2.0)], + [0.0, 1.0 / + sqrt(2.0), 0.0, -1.0 / sqrt(2.0), 0.0], + [0.0, 1.0 / sqrt(2.0), 0.0, 1.0 / sqrt(2.0), 0.0]]) # Spin blocks to be read: - n_spin_blocs = SP + 1 - SO # number of spins to read for Norbs and Ham, NOT Projectors - - # define the number of n_orbitals for all k points: it is the number of total bands and independent of k! - n_orbitals = numpy.ones([n_k,n_spin_blocs],numpy.int) * sum([ sh['dim'] for sh in shells ]) + # number of spins to read for Norbs and Ham, NOT Projectors + n_spin_blocs = SP + 1 - SO + + # define the number of n_orbitals for all k points: it is the + # number of total bands and independent of k! + n_orbitals = numpy.ones( + [n_k, n_spin_blocs], numpy.int) * sum([sh['dim'] for sh in shells]) # Initialise the projectors: - proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max([crsh['dim'] for crsh in corr_shells]),max(n_orbitals)],numpy.complex_) + proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max( + [crsh['dim'] for crsh in corr_shells]), max(n_orbitals)], numpy.complex_) # Read the projectors from the file: for ik in range(n_k): @@ -155,76 +183,90 @@ class HkConverter(ConverterTools): offset = 0 n_orb = 0 for ish in range(n_shells): - if (n_orb==0): - if (shells[ish]['atom']==corr_shells[icrsh]['atom']) and (shells[ish]['sort']==corr_shells[icrsh]['sort']): + if (n_orb == 0): + if (shells[ish]['atom'] == corr_shells[icrsh]['atom']) and (shells[ish]['sort'] == corr_shells[icrsh]['sort']): n_orb = corr_shells[icrsh]['dim'] else: offset += shells[ish]['dim'] - proj_mat[ik,isp,icrsh,0:n_orb,offset:offset+n_orb] = numpy.identity(n_orb) - + proj_mat[ik, isp, icrsh, 0:n_orb, + offset:offset + n_orb] = numpy.identity(n_orb) + # now define the arrays for weights and hopping ... - bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation - hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_) + # w(k_index), default normalisation + bz_weights = numpy.ones([n_k], numpy.float_) / float(n_k) + hopping = numpy.zeros([n_k, n_spin_blocs, max( + n_orbitals), max(n_orbitals)], numpy.complex_) if (weights_in_file): # weights in the file - for ik in range(n_k) : bz_weights[ik] = R.next() - + for ik in range(n_k): + bz_weights[ik] = R.next() + # if the sum over spins is in the weights, take it out again!! sm = sum(bz_weights) - bz_weights[:] /= sm + bz_weights[:] /= sm # Grab the H for isp in range(n_spin_blocs): - for ik in range(n_k) : - n_orb = n_orbitals[ik,isp] + for ik in range(n_k): + n_orb = n_orbitals[ik, isp] + + # first read all real components for given k, then read + # imaginary parts + if (first_real_part_matrix): - if (first_real_part_matrix): # first read all real components for given k, then read imaginary parts - for i in range(n_orb): if (only_upper_triangle): istart = i else: istart = 0 - for j in range(istart,n_orb): - hopping[ik,isp,i,j] = R.next() - + for j in range(istart, n_orb): + hopping[ik, isp, i, j] = R.next() + for i in range(n_orb): if (only_upper_triangle): istart = i else: istart = 0 - for j in range(istart,n_orb): - hopping[ik,isp,i,j] += R.next() * 1j - if ((only_upper_triangle)and(i!=j)): hopping[ik,isp,j,i] = hopping[ik,isp,i,j].conjugate() - - else: # read (real,im) tuple - + for j in range(istart, n_orb): + hopping[ik, isp, i, j] += R.next() * 1j + if ((only_upper_triangle)and(i != j)): + hopping[ik, isp, j, i] = hopping[ + ik, isp, i, j].conjugate() + + else: # read (real,im) tuple + for i in range(n_orb): if (only_upper_triangle): istart = i else: istart = 0 - for j in range(istart,n_orb): - hopping[ik,isp,i,j] = R.next() - hopping[ik,isp,i,j] += R.next() * 1j - - if ((only_upper_triangle)and(i!=j)): hopping[ik,isp,j,i] = hopping[ik,isp,i,j].conjugate() + for j in range(istart, n_orb): + hopping[ik, isp, i, j] = R.next() + hopping[ik, isp, i, j] += R.next() * 1j + + if ((only_upper_triangle)and(i != j)): + hopping[ik, isp, j, i] = hopping[ + ik, isp, i, j].conjugate() # keep some things that we need for reading parproj: - things_to_set = ['n_shells','shells','n_corr_shells','corr_shells','n_spin_blocs','n_orbitals','n_k','SO','SP','energy_unit'] - for it in things_to_set: setattr(self,it,locals()[it]) - except StopIteration : # a more explicit error if the file is corrupted. + things_to_set = ['n_shells', 'shells', 'n_corr_shells', 'corr_shells', + 'n_spin_blocs', 'n_orbitals', 'n_k', 'SO', 'SP', 'energy_unit'] + for it in things_to_set: + setattr(self, it, locals()[it]) + except StopIteration: # a more explicit error if the file is corrupted. raise "HK Converter : reading file dft_file failed!" R.close() # Save to the HDF5: - ar = HDFArchive(self.hdf_file,'a') - if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp) - things_to_save = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required', - 'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat', - 'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping', + ar = HDFArchive(self.hdf_file, 'a') + if not (self.dft_subgrp in ar): + ar.create_group(self.dft_subgrp) + things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required', + 'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat', + 'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping', 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr'] - for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it] + for it in things_to_save: + ar[self.dft_subgrp][it] = locals()[it] del ar diff --git a/python/converters/wannier90_converter.py b/python/converters/wannier90_converter.py index ae748172..14f2f71b 100644 --- a/python/converters/wannier90_converter.py +++ b/python/converters/wannier90_converter.py @@ -91,7 +91,8 @@ class Wannier90Converter(ConverterTools): self.dft_subgrp = dft_subgrp self.symmcorr_subgrp = symmcorr_subgrp self.fortran_to_replace = {'D': 'E'} - # threshold below which matrix elements from wannier90 should be considered equal + # threshold below which matrix elements from wannier90 should be + # considered equal self._w90zero = 2.e-6 # Checks if h5 file is there and repacks it if wanted: @@ -114,12 +115,14 @@ class Wannier90Converter(ConverterTools): return mpi.report("Reading input from %s..." % self.inp_file) - # R is a generator : each R.Next() will return the next number in the file + # R is a generator : each R.Next() will return the next number in the + # file R = ConverterTools.read_fortran_file( self, self.inp_file, self.fortran_to_replace) shell_entries = ['atom', 'sort', 'l', 'dim'] corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep'] - # First, let's read the input file with the parameters needed for the conversion + # First, let's read the input file with the parameters needed for the + # conversion try: # read k - point mesh generation option kmesh_mode = int(R.next()) @@ -135,7 +138,8 @@ class Wannier90Converter(ConverterTools): # and the data will be copied from corr_shells into shells (see below) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, n_corr_shells = int(R.next()) - # now read the information about the correlated shells (atom, sort, l, dim, SO flag, irep): + # now read the information about the correlated shells (atom, sort, + # l, dim, SO flag, irep): corr_shells = [{name: int(val) for name, val in zip( corr_shell_entries, R)} for icrsh in range(n_corr_shells)] except StopIteration: # a more explicit error if the file is corrupted. @@ -147,7 +151,7 @@ class Wannier90Converter(ConverterTools): # Set or derive some quantities # Wannier90 does not use symmetries to reduce the k-points # the following might change in future versions - symm_op = 0 + symm_op = 0 # copy corr_shells into shells (see above) n_shells = n_corr_shells shells = [] @@ -166,7 +170,8 @@ class Wannier90Converter(ConverterTools): mpi.report( "Total number of WFs expected in the correlated shells: %d" % dim_corr_shells) - # determine the number of inequivalent correlated shells and maps, needed for further processing + # determine the number of inequivalent correlated shells and maps, + # needed for further processing n_inequiv_shells, corr_to_inequiv, inequiv_to_corr = ConverterTools.det_shell_equivalence( self, corr_shells) mpi.report("Number of inequivalent shells: %d" % n_inequiv_shells) @@ -176,7 +181,8 @@ class Wannier90Converter(ConverterTools): mpi.report("Mapping: " + format(shells_map)) # build the k-point mesh, if its size was given on input (kmesh_mode >= 0), - # otherwise it is built according to the data in the hr file (see below) + # otherwise it is built according to the data in the hr file (see + # below) if kmesh_mode >= 0: n_k, k_mesh, bz_weights = self.kmesh_build(nki, kmesh_mode) self.n_k = n_k @@ -197,7 +203,8 @@ class Wannier90Converter(ConverterTools): # TODO: generalise to SP=1 (only partially done) rot_mat_time_inv = [0 for i in range(n_corr_shells)] - # Second, let's read the file containing the Hamiltonian in WF basis produced by Wannier90 + # Second, let's read the file containing the Hamiltonian in WF basis + # produced by Wannier90 for isp in range(n_spin): # begin loop on isp @@ -212,20 +219,24 @@ class Wannier90Converter(ConverterTools): mpi.report( "The Hamiltonian in MLWF basis is extracted from %s ..." % hr_file) nr, rvec, rdeg, nw, hamr = self.read_wannier90hr(hr_file) - # number of R vectors, their indices, their degeneracy, number of WFs, H(R) + # number of R vectors, their indices, their degeneracy, number of + # WFs, H(R) mpi.report("... done: %d R vectors, %d WFs found" % (nr, nw)) if isp == 0: - # set or check some quantities that must be the same for both spins + # set or check some quantities that must be the same for both + # spins self.nrpt = nr # k-point grid: (if not defined before) if kmesh_mode == -1: - # the size of the k-point mesh is determined from the largest R vector + # the size of the k-point mesh is determined from the + # largest R vector nki = [2 * rvec[:, idir].max() + 1 for idir in range(3)] # it will be the same as in the win only when nki is odd, because of the # wannier90 convention: if we have nki k-points along the i-th direction, - # then we should get 2*(nki/2)+nki%2 R points along that direction + # then we should get 2*(nki/2)+nki%2 R points along that + # direction n_k, k_mesh, bz_weights = self.kmesh_build(nki) self.n_k = n_k self.k_mesh = k_mesh @@ -237,33 +248,41 @@ class Wannier90Converter(ConverterTools): self.nwfs = nw # check that the total number of WFs makes sense if self.nwfs < dim_corr_shells: - mpi.report("ERROR: number of WFs in the file smaller than number of correlated orbitals!") + mpi.report( + "ERROR: number of WFs in the file smaller than number of correlated orbitals!") elif self.nwfs > dim_corr_shells: - # NOTE: correlated shells must appear before uncorrelated ones inside the file + # NOTE: correlated shells must appear before uncorrelated + # ones inside the file mpi.report("Number of WFs larger than correlated orbitals:\n" + "WFs from %d to %d treated as uncorrelated" % (dim_corr_shells + 1, self.nwfs)) else: - mpi.report("Number of WFs equal to number of correlated orbitals") + mpi.report( + "Number of WFs equal to number of correlated orbitals") - # we assume spin up and spin down always have same total number of WFs + # we assume spin up and spin down always have same total number + # of WFs n_orbitals = numpy.ones( [self.n_k, n_spin], numpy.int) * self.nwfs else: # consistency check between the _up and _down file contents if nr != self.nrpt: - mpi.report("Different number of R vectors for spin-up/spin-down!") + mpi.report( + "Different number of R vectors for spin-up/spin-down!") if nw != self.nwfs: - mpi.report("Different number of WFs for spin-up/spin-down!") + mpi.report( + "Different number of WFs for spin-up/spin-down!") hamr_full.append(hamr) # FIXME: when do we actually need deepcopy()? # hamr_full.append(deepcopy(hamr)) for ir in range(nr): - # checks if the Hamiltonian is real (it should, if wannierisation worked fine) + # checks if the Hamiltonian is real (it should, if + # wannierisation worked fine) if numpy.abs((hamr[ir].imag.max()).max()) > self._w90zero: - mpi.report("H(R) has large complex components at R %d" % ir) + mpi.report( + "H(R) has large complex components at R %d" % ir) # copy the R=0 block corresponding to the correlated shells # into another variable (needed later for finding rot_mat) if rvec[ir, 0] == 0 and rvec[ir, 1] == 0 and rvec[ir, 2] == 0: @@ -273,17 +292,22 @@ class Wannier90Converter(ConverterTools): if not numpy.allclose(ham_corr0.transpose().conjugate(), ham_corr0, atol=self._w90zero, rtol=1.e-9): raise ValueError("H(R=0) matrix is not Hermitian!") - # find rot_mat symmetries by diagonalising the on-site Hamiltonian of the first spin + # find rot_mat symmetries by diagonalising the on-site Hamiltonian + # of the first spin if isp == 0: - use_rotations, rot_mat = self.find_rot_mat(n_corr_shells, corr_shells, shells_map, ham_corr0) + use_rotations, rot_mat = self.find_rot_mat( + n_corr_shells, corr_shells, shells_map, ham_corr0) else: # consistency check - use_rotations_, rot_mat_ = self.find_rot_mat(n_corr_shells, corr_shells, shells_map, ham_corr0) + use_rotations_, rot_mat_ = self.find_rot_mat( + n_corr_shells, corr_shells, shells_map, ham_corr0) if (use_rotations and not use_rotations_): - mpi.report("Rotations cannot be used for spin component n. %d" % isp) + mpi.report( + "Rotations cannot be used for spin component n. %d" % isp) for icrsh in range(n_corr_shells): if not numpy.allclose(rot_mat_[icrsh], rot_mat[icrsh], atol=self._w90zero, rtol=1.e-15): - mpi.report("Rotations for spin component n. %d do not match!" % isp) + mpi.report( + "Rotations for spin component n. %d do not match!" % isp) # end loop on isp mpi.report("The k-point grid has dimensions: %d, %d, %d" % tuple(nki)) @@ -292,11 +316,14 @@ class Wannier90Converter(ConverterTools): bz_weights = 0.5 * bz_weights # Third, compute the hoppings in reciprocal space - hopping = numpy.zeros([self.n_k, n_spin, numpy.max(n_orbitals), numpy.max(n_orbitals)], numpy.complex_) + hopping = numpy.zeros([self.n_k, n_spin, numpy.max( + n_orbitals), numpy.max(n_orbitals)], numpy.complex_) for isp in range(n_spin): - # make Fourier transform H(R) -> H(k) : it can be done one spin at a time + # make Fourier transform H(R) -> H(k) : it can be done one spin at + # a time hamk = self.fourier_ham(self.nwfs, hamr_full[isp]) - # copy the H(k) in the right place of hoppings... is there a better way to do this?? + # copy the H(k) in the right place of hoppings... is there a better + # way to do this?? for ik in range(self.n_k): #hopping[ik,isp,:,:] = deepcopy(hamk[ik][:,:])*energy_unit hopping[ik, isp, :, :] = hamk[ik][:, :] * energy_unit @@ -309,7 +336,8 @@ class Wannier90Converter(ConverterTools): # Projectors simply consist in identity matrix blocks selecting those MLWFs that # correspond to the specific correlated shell indexed by icrsh. # NOTE: we assume that the correlated orbitals appear at the beginning of the H(R) - # file and that the ordering of MLWFs matches the corr_shell info from the input. + # file and that the ordering of MLWFs matches the corr_shell info from + # the input. for icrsh in range(n_corr_shells): norb = corr_shells[icrsh]['dim'] proj_mat[:, :, icrsh, 0:norb, iorb:iorb + @@ -320,7 +348,8 @@ class Wannier90Converter(ConverterTools): ar = HDFArchive(self.hdf_file, 'a') if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp) - # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten! + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten! things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required', 'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat', 'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping', @@ -373,7 +402,8 @@ class Wannier90Converter(ConverterTools): except ValueError: mpi.report("Could not read number of WFs or R vectors") - # allocate arrays to save the R vector indexes and degeneracies and the Hamiltonian + # allocate arrays to save the R vector indexes and degeneracies and the + # Hamiltonian rvec_idx = numpy.zeros((nrpt, 3), dtype=int) rvec_deg = numpy.zeros(nrpt, dtype=int) h_of_r = [numpy.zeros((num_wf, num_wf), dtype=numpy.complex_) @@ -383,7 +413,8 @@ class Wannier90Converter(ConverterTools): currpos = 2 try: ir = 0 - # read the degeneracy of the R vectors (needed for the Fourier transform) + # read the degeneracy of the R vectors (needed for the Fourier + # transform) while ir < nrpt: currpos += 1 for x in hr_data[currpos].split(): @@ -540,7 +571,8 @@ class Wannier90Converter(ConverterTools): kmesh = numpy.zeros((nkpt, 3), dtype=float) ii = 0 for ix, iy, iz in product(range(msize[0]), range(msize[1]), range(msize[2])): - kmesh[ii, :] = [float(ix) / msize[0], float(iy) / msize[1], float(iz) / msize[2]] + kmesh[ii, :] = [float(ix) / msize[0], float(iy) / + msize[1], float(iz) / msize[2]] ii += 1 # weight is equal for all k-points because wannier90 uses uniform grid on whole BZ # (normalization is always 1 and takes into account spin degeneracy) @@ -568,11 +600,13 @@ class Wannier90Converter(ConverterTools): """ twopi = 2 * numpy.pi - h_of_k = [numpy.zeros((norb, norb), dtype=numpy.complex_) for ik in range(self.n_k)] + h_of_k = [numpy.zeros((norb, norb), dtype=numpy.complex_) + for ik in range(self.n_k)] ridx = numpy.array(range(self.nrpt)) for ik, ir in product(range(self.n_k), ridx): rdotk = twopi * numpy.dot(self.k_mesh[ik], self.rvec[ir]) - factor = (math.cos(rdotk) + 1j * math.sin(rdotk)) / float(self.rdeg[ir]) + factor = (math.cos(rdotk) + 1j * math.sin(rdotk)) / \ + float(self.rdeg[ir]) h_of_k[ik][:, :] += factor * h_of_r[ir][:, :] return h_of_k diff --git a/python/converters/wien2k_converter.py b/python/converters/wien2k_converter.py index 8be718c1..211dfc1c 100644 --- a/python/converters/wien2k_converter.py +++ b/python/converters/wien2k_converter.py @@ -1,5 +1,5 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,7 +18,7 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## from types import * import numpy @@ -26,16 +26,17 @@ from pytriqs.archive import * from converter_tools import * import os.path + class Wien2kConverter(ConverterTools): """ Conversion from Wien2k output to an hdf5 file that can be used as input for the SumkDFT class. """ - def __init__(self, filename, hdf_filename = None, - dft_subgrp = 'dft_input', symmcorr_subgrp = 'dft_symmcorr_input', - parproj_subgrp='dft_parproj_input', symmpar_subgrp='dft_symmpar_input', - bands_subgrp = 'dft_bands_input', misc_subgrp = 'dft_misc_input', - transp_subgrp = 'dft_transp_input', repacking = False): + def __init__(self, filename, hdf_filename=None, + dft_subgrp='dft_input', symmcorr_subgrp='dft_symmcorr_input', + parproj_subgrp='dft_parproj_input', symmpar_subgrp='dft_symmpar_input', + bands_subgrp='dft_bands_input', misc_subgrp='dft_misc_input', + transp_subgrp='dft_transp_input', repacking=False): """ Initialise the class. @@ -61,21 +62,23 @@ class Wien2kConverter(ConverterTools): Name of subgroup storing transport data. repacking : boolean, optional Does the hdf5 archive need to be repacked to save space? - + """ - assert type(filename)==StringType, "Wien2kConverter: Please provide the DFT files' base name as a string." - if hdf_filename is None: hdf_filename = filename+'.h5' + assert type( + filename) == StringType, "Wien2kConverter: Please provide the DFT files' base name as a string." + if hdf_filename is None: + hdf_filename = filename + '.h5' self.hdf_file = hdf_filename - self.dft_file = filename+'.ctqmcout' - self.symmcorr_file = filename+'.symqmc' - self.parproj_file = filename+'.parproj' - self.symmpar_file = filename+'.sympar' - self.band_file = filename+'.outband' - self.bandwin_file = filename+'.oubwin' - self.struct_file = filename+'.struct' - self.outputs_file = filename+'.outputs' - self.pmat_file = filename+'.pmat' + self.dft_file = filename + '.ctqmcout' + self.symmcorr_file = filename + '.symqmc' + self.parproj_file = filename + '.parproj' + self.symmpar_file = filename + '.sympar' + self.band_file = filename + '.outband' + self.bandwin_file = filename + '.oubwin' + self.struct_file = filename + '.struct' + self.outputs_file = filename + '.outputs' + self.pmat_file = filename + '.pmat' self.dft_subgrp = dft_subgrp self.symmcorr_subgrp = symmcorr_subgrp self.parproj_subgrp = parproj_subgrp @@ -83,13 +86,12 @@ class Wien2kConverter(ConverterTools): self.bands_subgrp = bands_subgrp self.misc_subgrp = misc_subgrp self.transp_subgrp = transp_subgrp - self.fortran_to_replace = {'D':'E'} + self.fortran_to_replace = {'D': 'E'} # Checks if h5 file is there and repacks it if wanted: if (os.path.exists(self.hdf_file) and repacking): ConverterTools.repack(self) - def convert_dft_input(self): """ Reads the appropriate files and stores the data for the @@ -101,149 +103,180 @@ class Wien2kConverter(ConverterTools): in the hdf5 archive. """ - - # Read and write only on the master node - if not (mpi.is_master_node()): return - mpi.report("Reading input from %s..."%self.dft_file) - # R is a generator : each R.Next() will return the next number in the file - R = ConverterTools.read_fortran_file(self,self.dft_file,self.fortran_to_replace) + # Read and write only on the master node + if not (mpi.is_master_node()): + return + mpi.report("Reading input from %s..." % self.dft_file) + + # R is a generator : each R.Next() will return the next number in the + # file + R = ConverterTools.read_fortran_file( + self, self.dft_file, self.fortran_to_replace) try: energy_unit = R.next() # read the energy convertion factor - n_k = int(R.next()) # read the number of k points - k_dep_projection = 1 - SP = int(R.next()) # flag for spin-polarised calculation - SO = int(R.next()) # flag for spin-orbit calculation + # read the number of k points + n_k = int(R.next()) + k_dep_projection = 1 + # flag for spin-polarised calculation + SP = int(R.next()) + # flag for spin-orbit calculation + SO = int(R.next()) charge_below = R.next() # total charge below energy window - density_required = R.next() # total density required, for setting the chemical potential + # total density required, for setting the chemical potential + density_required = R.next() symm_op = 1 # Use symmetry groups for the k-sum - # the information on the non-correlated shells is not important here, maybe skip: - n_shells = int(R.next()) # number of shells (e.g. Fe d, As p, O p) in the unit cell, - # corresponds to index R in formulas + # the information on the non-correlated shells is not important + # here, maybe skip: + # number of shells (e.g. Fe d, As p, O p) in the unit cell, + n_shells = int(R.next()) + # corresponds to index R in formulas # now read the information about the shells (atom, sort, l, dim): shell_entries = ['atom', 'sort', 'l', 'dim'] - shells = [ {name: int(val) for name, val in zip(shell_entries, R)} for ish in range(n_shells) ] + shells = [{name: int(val) for name, val in zip( + shell_entries, R)} for ish in range(n_shells)] - n_corr_shells = int(R.next()) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, - # corresponds to index R in formulas - # now read the information about the shells (atom, sort, l, dim, SO flag, irep): + # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, + n_corr_shells = int(R.next()) + # corresponds to index R in formulas + # now read the information about the shells (atom, sort, l, dim, SO + # flag, irep): corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep'] - corr_shells = [ {name: int(val) for name, val in zip(corr_shell_entries, R)} for icrsh in range(n_corr_shells) ] + corr_shells = [{name: int(val) for name, val in zip( + corr_shell_entries, R)} for icrsh in range(n_corr_shells)] - # determine the number of inequivalent correlated shells and maps, needed for further reading - n_inequiv_shells, corr_to_inequiv, inequiv_to_corr = ConverterTools.det_shell_equivalence(self,corr_shells) + # determine the number of inequivalent correlated shells and maps, + # needed for further reading + n_inequiv_shells, corr_to_inequiv, inequiv_to_corr = ConverterTools.det_shell_equivalence( + self, corr_shells) use_rotations = 1 - rot_mat = [numpy.identity(corr_shells[icrsh]['dim'],numpy.complex_) for icrsh in range(n_corr_shells)] - + rot_mat = [numpy.identity( + corr_shells[icrsh]['dim'], numpy.complex_) for icrsh in range(n_corr_shells)] + # read the matrices rot_mat_time_inv = [0 for i in range(n_corr_shells)] for icrsh in range(n_corr_shells): for i in range(corr_shells[icrsh]['dim']): # read real part: for j in range(corr_shells[icrsh]['dim']): - rot_mat[icrsh][i,j] = R.next() - for i in range(corr_shells[icrsh]['dim']): # read imaginary part: + rot_mat[icrsh][i, j] = R.next() + # read imaginary part: + for i in range(corr_shells[icrsh]['dim']): for j in range(corr_shells[icrsh]['dim']): - rot_mat[icrsh][i,j] += 1j * R.next() + rot_mat[icrsh][i, j] += 1j * R.next() - if (SP==1): # read time inversion flag: + if (SP == 1): # read time inversion flag: rot_mat_time_inv[icrsh] = int(R.next()) - + # Read here the info for the transformation of the basis: n_reps = [1 for i in range(n_inequiv_shells)] dim_reps = [0 for i in range(n_inequiv_shells)] T = [] for ish in range(n_inequiv_shells): - n_reps[ish] = int(R.next()) # number of representatives ("subsets"), e.g. t2g and eg - dim_reps[ish] = [int(R.next()) for i in range(n_reps[ish])] # dimensions of the subsets - + # number of representatives ("subsets"), e.g. t2g and eg + n_reps[ish] = int(R.next()) + dim_reps[ish] = [int(R.next()) for i in range( + n_reps[ish])] # dimensions of the subsets + # The transformation matrix: # is of dimension 2l+1 without SO, and 2*(2l+1) with SO! - ll = 2*corr_shells[inequiv_to_corr[ish]]['l']+1 + ll = 2 * corr_shells[inequiv_to_corr[ish]]['l'] + 1 lmax = ll * (corr_shells[inequiv_to_corr[ish]]['SO'] + 1) - T.append(numpy.zeros([lmax,lmax],numpy.complex_)) - + T.append(numpy.zeros([lmax, lmax], numpy.complex_)) + # now read it from file: for i in range(lmax): for j in range(lmax): - T[ish][i,j] = R.next() + T[ish][i, j] = R.next() for i in range(lmax): for j in range(lmax): - T[ish][i,j] += 1j * R.next() - + T[ish][i, j] += 1j * R.next() + # Spin blocks to be read: - n_spin_blocs = SP + 1 - SO - + n_spin_blocs = SP + 1 - SO + # read the list of n_orbitals for all k points - n_orbitals = numpy.zeros([n_k,n_spin_blocs],numpy.int) + n_orbitals = numpy.zeros([n_k, n_spin_blocs], numpy.int) for isp in range(n_spin_blocs): for ik in range(n_k): - n_orbitals[ik,isp] = int(R.next()) - + n_orbitals[ik, isp] = int(R.next()) + # Initialise the projectors: - proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max([crsh['dim'] for crsh in corr_shells]),numpy.max(n_orbitals)],numpy.complex_) + proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, max( + [crsh['dim'] for crsh in corr_shells]), numpy.max(n_orbitals)], numpy.complex_) # Read the projectors from the file: for ik in range(n_k): for icrsh in range(n_corr_shells): n_orb = corr_shells[icrsh]['dim'] - # first Real part for BOTH spins, due to conventions in dmftproj: + # first Real part for BOTH spins, due to conventions in + # dmftproj: for isp in range(n_spin_blocs): for i in range(n_orb): for j in range(n_orbitals[ik][isp]): - proj_mat[ik,isp,icrsh,i,j] = R.next() + proj_mat[ik, isp, icrsh, i, j] = R.next() # now Imag part: for isp in range(n_spin_blocs): for i in range(n_orb): for j in range(n_orbitals[ik][isp]): - proj_mat[ik,isp,icrsh,i,j] += 1j * R.next() - + proj_mat[ik, isp, icrsh, i, j] += 1j * R.next() + # now define the arrays for weights and hopping ... - bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation - hopping = numpy.zeros([n_k,n_spin_blocs,numpy.max(n_orbitals),numpy.max(n_orbitals)],numpy.complex_) + # w(k_index), default normalisation + bz_weights = numpy.ones([n_k], numpy.float_) / float(n_k) + hopping = numpy.zeros([n_k, n_spin_blocs, numpy.max( + n_orbitals), numpy.max(n_orbitals)], numpy.complex_) # weights in the file - for ik in range(n_k) : bz_weights[ik] = R.next() - + for ik in range(n_k): + bz_weights[ik] = R.next() + # if the sum over spins is in the weights, take it out again!! sm = sum(bz_weights) - bz_weights[:] /= sm + bz_weights[:] /= sm # Grab the H - # we use now the convention of a DIAGONAL Hamiltonian -- convention for Wien2K. + # we use now the convention of a DIAGONAL Hamiltonian -- convention + # for Wien2K. for isp in range(n_spin_blocs): - for ik in range(n_k) : - n_orb = n_orbitals[ik,isp] + for ik in range(n_k): + n_orb = n_orbitals[ik, isp] for i in range(n_orb): - hopping[ik,isp,i,i] = R.next() * energy_unit - + hopping[ik, isp, i, i] = R.next() * energy_unit + # keep some things that we need for reading parproj: - things_to_set = ['n_shells','shells','n_corr_shells','corr_shells','n_spin_blocs','n_orbitals','n_k','SO','SP','energy_unit'] - for it in things_to_set: setattr(self,it,locals()[it]) - except StopIteration : # a more explicit error if the file is corrupted. - raise "Wien2k_converter : reading file %s failed!"%self.dft_file + things_to_set = ['n_shells', 'shells', 'n_corr_shells', 'corr_shells', + 'n_spin_blocs', 'n_orbitals', 'n_k', 'SO', 'SP', 'energy_unit'] + for it in things_to_set: + setattr(self, it, locals()[it]) + except StopIteration: # a more explicit error if the file is corrupted. + raise "Wien2k_converter : reading file %s failed!" % self.dft_file R.close() # Reading done! - + # Save it to the HDF: - ar = HDFArchive(self.hdf_file,'a') - if not (self.dft_subgrp in ar): ar.create_group(self.dft_subgrp) - # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten! - things_to_save = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required', - 'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat', - 'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping', + ar = HDFArchive(self.hdf_file, 'a') + if not (self.dft_subgrp in ar): + ar.create_group(self.dft_subgrp) + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten! + things_to_save = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required', + 'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat', + 'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping', 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr'] - for it in things_to_save: ar[self.dft_subgrp][it] = locals()[it] + for it in things_to_save: + ar[self.dft_subgrp][it] = locals()[it] del ar - # Symmetries are used, so now convert symmetry information for *correlated* orbitals: - self.convert_symmetry_input(orbits=self.corr_shells,symm_file=self.symmcorr_file,symm_subgrp=self.symmcorr_subgrp,SO=self.SO,SP=self.SP) + # Symmetries are used, so now convert symmetry information for + # *correlated* orbitals: + self.convert_symmetry_input(orbits=self.corr_shells, symm_file=self.symmcorr_file, + symm_subgrp=self.symmcorr_subgrp, SO=self.SO, SP=self.SP) self.convert_misc_input() - def convert_parproj_input(self): """ Reads the appropriate files and stores the data for the @@ -255,31 +288,37 @@ class Wien2kConverter(ConverterTools): """ - if not (mpi.is_master_node()): return + if not (mpi.is_master_node()): + return # get needed data from hdf file - ar = HDFArchive(self.hdf_file,'a') - things_to_read = ['SP','SO','n_shells','n_k','n_orbitals','shells'] + ar = HDFArchive(self.hdf_file, 'a') + things_to_read = ['SP', 'SO', 'n_shells', + 'n_k', 'n_orbitals', 'shells'] for it in things_to_read: - if not hasattr(self,it): setattr(self,it,ar[self.dft_subgrp][it]) + if not hasattr(self, it): + setattr(self, it, ar[self.dft_subgrp][it]) self.n_spin_blocs = self.SP + 1 - self.SO del ar - mpi.report("Reading input from %s..."%self.parproj_file) + mpi.report("Reading input from %s..." % self.parproj_file) - dens_mat_below = [ [numpy.zeros([self.shells[ish]['dim'],self.shells[ish]['dim']],numpy.complex_) for ish in range(self.n_shells)] - for isp in range(self.n_spin_blocs) ] + dens_mat_below = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], numpy.complex_) for ish in range(self.n_shells)] + for isp in range(self.n_spin_blocs)] - R = ConverterTools.read_fortran_file(self,self.parproj_file,self.fortran_to_replace) + R = ConverterTools.read_fortran_file( + self, self.parproj_file, self.fortran_to_replace) n_parproj = [int(R.next()) for i in range(self.n_shells)] n_parproj = numpy.array(n_parproj) - + # Initialise P, here a double list of matrices: - proj_mat_all = numpy.zeros([self.n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max([sh['dim'] for sh in self.shells]),max(self.n_orbitals)],numpy.complex_) - - rot_mat_all = [numpy.identity(self.shells[ish]['dim'],numpy.complex_) for ish in range(self.n_shells)] + proj_mat_all = numpy.zeros([self.n_k, self.n_spin_blocs, self.n_shells, max( + n_parproj), max([sh['dim'] for sh in self.shells]), max(self.n_orbitals)], numpy.complex_) + + rot_mat_all = [numpy.identity( + self.shells[ish]['dim'], numpy.complex_) for ish in range(self.n_shells)] rot_mat_all_time_inv = [0 for i in range(self.n_shells)] for ish in range(self.n_shells): @@ -288,35 +327,40 @@ class Wien2kConverter(ConverterTools): for ir in range(n_parproj[ish]): for isp in range(self.n_spin_blocs): - for i in range(self.shells[ish]['dim']): # read real part: + # read real part: + for i in range(self.shells[ish]['dim']): for j in range(self.n_orbitals[ik][isp]): - proj_mat_all[ik,isp,ish,ir,i,j] = R.next() - + proj_mat_all[ik, isp, ish, ir, i, j] = R.next() + for isp in range(self.n_spin_blocs): - for i in range(self.shells[ish]['dim']): # read imaginary part: + # read imaginary part: + for i in range(self.shells[ish]['dim']): for j in range(self.n_orbitals[ik][isp]): - proj_mat_all[ik,isp,ish,ir,i,j] += 1j * R.next() - - - # now read the Density Matrix for this orbital below the energy window: + proj_mat_all[ik, isp, ish, + ir, i, j] += 1j * R.next() + + # now read the Density Matrix for this orbital below the energy + # window: for isp in range(self.n_spin_blocs): for i in range(self.shells[ish]['dim']): # read real part: for j in range(self.shells[ish]['dim']): - dens_mat_below[isp][ish][i,j] = R.next() + dens_mat_below[isp][ish][i, j] = R.next() for isp in range(self.n_spin_blocs): - for i in range(self.shells[ish]['dim']): # read imaginary part: + # read imaginary part: + for i in range(self.shells[ish]['dim']): for j in range(self.shells[ish]['dim']): - dens_mat_below[isp][ish][i,j] += 1j * R.next() - if (self.SP==0): dens_mat_below[isp][ish] /= 2.0 + dens_mat_below[isp][ish][i, j] += 1j * R.next() + if (self.SP == 0): + dens_mat_below[isp][ish] /= 2.0 # Global -> local rotation matrix for this shell: for i in range(self.shells[ish]['dim']): # read real part: for j in range(self.shells[ish]['dim']): - rot_mat_all[ish][i,j] = R.next() + rot_mat_all[ish][i, j] = R.next() for i in range(self.shells[ish]['dim']): # read imaginary part: for j in range(self.shells[ish]['dim']): - rot_mat_all[ish][i,j] += 1j * R.next() - + rot_mat_all[ish][i, j] += 1j * R.next() + if (self.SP): rot_mat_all_time_inv[ish] = int(R.next()) @@ -324,16 +368,21 @@ class Wien2kConverter(ConverterTools): # Reading done! # Save it to the HDF: - ar = HDFArchive(self.hdf_file,'a') - if not (self.parproj_subgrp in ar): ar.create_group(self.parproj_subgrp) - # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten! - things_to_save = ['dens_mat_below','n_parproj','proj_mat_all','rot_mat_all','rot_mat_all_time_inv'] - for it in things_to_save: ar[self.parproj_subgrp][it] = locals()[it] + ar = HDFArchive(self.hdf_file, 'a') + if not (self.parproj_subgrp in ar): + ar.create_group(self.parproj_subgrp) + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten! + things_to_save = ['dens_mat_below', 'n_parproj', + 'proj_mat_all', 'rot_mat_all', 'rot_mat_all_time_inv'] + for it in things_to_save: + ar[self.parproj_subgrp][it] = locals()[it] del ar - # Symmetries are used, so now convert symmetry information for *all* orbitals: - self.convert_symmetry_input(orbits=self.shells,symm_file=self.symmpar_file,symm_subgrp=self.symmpar_subgrp,SO=self.SO,SP=self.SP) - + # Symmetries are used, so now convert symmetry information for *all* + # orbitals: + self.convert_symmetry_input(orbits=self.shells, symm_file=self.symmpar_file, + symm_subgrp=self.symmpar_subgrp, SO=self.SO, SP=self.SP) def convert_bands_input(self): """ @@ -341,117 +390,134 @@ class Wien2kConverter(ConverterTools): """ - if not (mpi.is_master_node()): return + if not (mpi.is_master_node()): + return try: # get needed data from hdf file - ar = HDFArchive(self.hdf_file,'a') - things_to_read = ['SP','SO','n_corr_shells','n_shells','corr_shells','shells','energy_unit'] + ar = HDFArchive(self.hdf_file, 'a') + things_to_read = ['SP', 'SO', 'n_corr_shells', + 'n_shells', 'corr_shells', 'shells', 'energy_unit'] for it in things_to_read: - if not hasattr(self,it): setattr(self,it,ar[self.dft_subgrp][it]) + if not hasattr(self, it): + setattr(self, it, ar[self.dft_subgrp][it]) self.n_spin_blocs = self.SP + 1 - self.SO del ar - mpi.report("Reading input from %s..."%self.band_file) - R = ConverterTools.read_fortran_file(self,self.band_file,self.fortran_to_replace) + mpi.report("Reading input from %s..." % self.band_file) + R = ConverterTools.read_fortran_file( + self, self.band_file, self.fortran_to_replace) n_k = int(R.next()) # read the list of n_orbitals for all k points - n_orbitals = numpy.zeros([n_k,self.n_spin_blocs],numpy.int) + n_orbitals = numpy.zeros([n_k, self.n_spin_blocs], numpy.int) for isp in range(self.n_spin_blocs): for ik in range(n_k): - n_orbitals[ik,isp] = int(R.next()) + n_orbitals[ik, isp] = int(R.next()) # Initialise the projectors: - proj_mat = numpy.zeros([n_k,self.n_spin_blocs,self.n_corr_shells,max([crsh['dim'] for crsh in self.corr_shells]),numpy.max(n_orbitals)],numpy.complex_) + proj_mat = numpy.zeros([n_k, self.n_spin_blocs, self.n_corr_shells, max( + [crsh['dim'] for crsh in self.corr_shells]), numpy.max(n_orbitals)], numpy.complex_) # Read the projectors from the file: for ik in range(n_k): for icrsh in range(self.n_corr_shells): n_orb = self.corr_shells[icrsh]['dim'] - # first Real part for BOTH spins, due to conventions in dmftproj: + # first Real part for BOTH spins, due to conventions in + # dmftproj: for isp in range(self.n_spin_blocs): for i in range(n_orb): - for j in range(n_orbitals[ik,isp]): - proj_mat[ik,isp,icrsh,i,j] = R.next() + for j in range(n_orbitals[ik, isp]): + proj_mat[ik, isp, icrsh, i, j] = R.next() # now Imag part: for isp in range(self.n_spin_blocs): for i in range(n_orb): - for j in range(n_orbitals[ik,isp]): - proj_mat[ik,isp,icrsh,i,j] += 1j * R.next() + for j in range(n_orbitals[ik, isp]): + proj_mat[ik, isp, icrsh, i, j] += 1j * R.next() + + hopping = numpy.zeros([n_k, self.n_spin_blocs, numpy.max( + n_orbitals), numpy.max(n_orbitals)], numpy.complex_) - hopping = numpy.zeros([n_k,self.n_spin_blocs,numpy.max(n_orbitals),numpy.max(n_orbitals)],numpy.complex_) - # Grab the H # we use now the convention of a DIAGONAL Hamiltonian!!!! for isp in range(self.n_spin_blocs): - for ik in range(n_k) : - n_orb = n_orbitals[ik,isp] + for ik in range(n_k): + n_orb = n_orbitals[ik, isp] for i in range(n_orb): - hopping[ik,isp,i,i] = R.next() * self.energy_unit + hopping[ik, isp, i, i] = R.next() * self.energy_unit # now read the partial projectors: n_parproj = [int(R.next()) for i in range(self.n_shells)] n_parproj = numpy.array(n_parproj) - + # Initialise P, here a double list of matrices: - proj_mat_all = numpy.zeros([n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max([sh['dim'] for sh in self.shells]),numpy.max(n_orbitals)],numpy.complex_) + proj_mat_all = numpy.zeros([n_k, self.n_spin_blocs, self.n_shells, max(n_parproj), max( + [sh['dim'] for sh in self.shells]), numpy.max(n_orbitals)], numpy.complex_) for ish in range(self.n_shells): for ik in range(n_k): for ir in range(n_parproj[ish]): for isp in range(self.n_spin_blocs): - - for i in range(self.shells[ish]['dim']): # read real part: - for j in range(n_orbitals[ik,isp]): - proj_mat_all[ik,isp,ish,ir,i,j] = R.next() - - for i in range(self.shells[ish]['dim']): # read imaginary part: - for j in range(n_orbitals[ik,isp]): - proj_mat_all[ik,isp,ish,ir,i,j] += 1j * R.next() + + # read real part: + for i in range(self.shells[ish]['dim']): + for j in range(n_orbitals[ik, isp]): + proj_mat_all[ik, isp, ish, + ir, i, j] = R.next() + + # read imaginary part: + for i in range(self.shells[ish]['dim']): + for j in range(n_orbitals[ik, isp]): + proj_mat_all[ik, isp, ish, + ir, i, j] += 1j * R.next() R.close() except KeyError: raise "convert_bands_input : Needed data not found in hdf file. Consider calling convert_dft_input first!" - except StopIteration : # a more explicit error if the file is corrupted. + except StopIteration: # a more explicit error if the file is corrupted. raise "Wien2k_converter : reading file band_file failed!" # Reading done! # Save it to the HDF: - ar = HDFArchive(self.hdf_file,'a') - if not (self.bands_subgrp in ar): ar.create_group(self.bands_subgrp) - # The subgroup containing the data. If it does not exist, it is created. If it exists, the data is overwritten! - things_to_save = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_all'] - for it in things_to_save: ar[self.bands_subgrp][it] = locals()[it] + ar = HDFArchive(self.hdf_file, 'a') + if not (self.bands_subgrp in ar): + ar.create_group(self.bands_subgrp) + # The subgroup containing the data. If it does not exist, it is + # created. If it exists, the data is overwritten! + things_to_save = ['n_k', 'n_orbitals', 'proj_mat', + 'hopping', 'n_parproj', 'proj_mat_all'] + for it in things_to_save: + ar[self.bands_subgrp][it] = locals()[it] del ar - def convert_misc_input(self): """ Reads additional information on: - + - the band window from :file:`case.oubwin`, - lattice parameters from :file:`case.struct`, - symmetries from :file:`case.outputs`, - + if those Wien2k files are present and stores the data in the hdf5 archive. This function is automatically called by :meth:`convert_dft_input `. """ - if not (mpi.is_master_node()): return - + if not (mpi.is_master_node()): + return + # Check if SP, SO and n_k are already in h5 ar = HDFArchive(self.hdf_file, 'r') - if not (self.dft_subgrp in ar): raise IOError, "convert_misc_input: No %s subgroup in hdf file found! Call convert_dft_input first." %self.dft_subgrp + if not (self.dft_subgrp in ar): + raise IOError, "convert_misc_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp SP = ar[self.dft_subgrp]['SP'] SO = ar[self.dft_subgrp]['SO'] n_k = ar[self.dft_subgrp]['n_k'] del ar - + things_to_save = [] # Read relevant data from .oubwin/up/dn files @@ -459,32 +525,35 @@ class Wien2kConverter(ConverterTools): # band_window: Contains the index of the lowest and highest band within the # projected subspace (used by dmftproj) for each k-point. - if (SP == 0 or SO == 1): + if (SP == 0 or SO == 1): files = [self.bandwin_file] elif SP == 1: - files = [self.bandwin_file+'up', self.bandwin_file+'dn'] - else: # SO and SP can't both be 1 + files = [self.bandwin_file + 'up', self.bandwin_file + 'dn'] + else: # SO and SP can't both be 1 assert 0, "convert_misc_input: Reading oubwin error! Check SP and SO!" band_window = [None for isp in range(SP + 1 - SO)] for isp, f in enumerate(files): if os.path.exists(f): - mpi.report("Reading input from %s..."%f) - R = ConverterTools.read_fortran_file(self, f, self.fortran_to_replace) + mpi.report("Reading input from %s..." % f) + R = ConverterTools.read_fortran_file( + self, f, self.fortran_to_replace) n_k_oubwin = int(R.next()) if (n_k_oubwin != n_k): - mpi.report("convert_misc_input : WARNING : n_k in case.oubwin is different from n_k in case.klist") - assert int(R.next()) == SO, "convert_misc_input: SO is inconsistent in oubwin file!" + mpi.report( + "convert_misc_input : WARNING : n_k in case.oubwin is different from n_k in case.klist") + assert int( + R.next()) == SO, "convert_misc_input: SO is inconsistent in oubwin file!" - band_window[isp] = numpy.zeros((n_k_oubwin, 2), dtype=int) + band_window[isp] = numpy.zeros((n_k_oubwin, 2), dtype=int) for ik in xrange(n_k_oubwin): R.next() - band_window[isp][ik,0] = R.next() # lowest band - band_window[isp][ik,1] = R.next() # highest band + band_window[isp][ik, 0] = R.next() # lowest band + band_window[isp][ik, 1] = R.next() # highest band R.next() things_to_save.append('band_window') - R.close() # Reading done! + R.close() # Reading done! # Read relevant data from .struct file ###################################### @@ -493,39 +562,44 @@ class Wien2kConverter(ConverterTools): # lattice_angles: unit cell angles in rad if (os.path.exists(self.struct_file)): - mpi.report("Reading input from %s..."%self.struct_file) - + mpi.report("Reading input from %s..." % self.struct_file) + with open(self.struct_file) as R: try: R.readline() lattice_type = R.readline().split()[0] R.readline() temp = R.readline() - lattice_constants = numpy.array([float(temp[0+10*i:10+10*i].strip()) for i in range(3)]) - lattice_angles = numpy.array([float(temp[30+10*i:40+10*i].strip()) for i in range(3)]) * numpy.pi / 180.0 - things_to_save.extend(['lattice_type', 'lattice_constants', 'lattice_angles']) + lattice_constants = numpy.array( + [float(temp[0 + 10 * i:10 + 10 * i].strip()) for i in range(3)]) + lattice_angles = numpy.array( + [float(temp[30 + 10 * i:40 + 10 * i].strip()) for i in range(3)]) * numpy.pi / 180.0 + things_to_save.extend( + ['lattice_type', 'lattice_constants', 'lattice_angles']) except IOError: - raise "convert_misc_input: reading file %s failed" %self.struct_file + raise "convert_misc_input: reading file %s failed" % self.struct_file # Read relevant data from .outputs file ####################################### - # rot_symmetries: matrix representation of all (space group) symmetry operations - + # rot_symmetries: matrix representation of all (space group) symmetry + # operations + if (os.path.exists(self.outputs_file)): - mpi.report("Reading input from %s..."%self.outputs_file) - + mpi.report("Reading input from %s..." % self.outputs_file) + rot_symmetries = [] with open(self.outputs_file) as R: try: while 1: temp = R.readline().strip(' ').split() - if (temp[0] =='PGBSYM:'): + if (temp[0] == 'PGBSYM:'): n_symmetries = int(temp[-1]) break for i in range(n_symmetries): while 1: - if (R.readline().strip().split()[0] == 'Symmetry'): break - sym_i = numpy.zeros((3, 3), dtype = float) + if (R.readline().strip().split()[0] == 'Symmetry'): + break + sym_i = numpy.zeros((3, 3), dtype=float) for ir in range(3): temp = R.readline().strip().split() for ic in range(3): @@ -535,30 +609,33 @@ class Wien2kConverter(ConverterTools): things_to_save.extend(['n_symmetries', 'rot_symmetries']) things_to_save.append('rot_symmetries') except IOError: - raise "convert_misc_input: reading file %s failed" %self.outputs_file + raise "convert_misc_input: reading file %s failed" % self.outputs_file # Save it to the HDF: - ar = HDFArchive(self.hdf_file,'a') - if not (self.misc_subgrp in ar): ar.create_group(self.misc_subgrp) - for it in things_to_save: ar[self.misc_subgrp][it] = locals()[it] + ar = HDFArchive(self.hdf_file, 'a') + if not (self.misc_subgrp in ar): + ar.create_group(self.misc_subgrp) + for it in things_to_save: + ar[self.misc_subgrp][it] = locals()[it] del ar - def convert_transport_input(self): """ Reads the necessary information for transport calculations on: - + - the optical band window and the velocity matrix elements from :file:`case.pmat` and stores the data in the hdf5 archive. - + """ - - if not (mpi.is_master_node()): return - + + if not (mpi.is_master_node()): + return + # Check if SP, SO and n_k are already in h5 ar = HDFArchive(self.hdf_file, 'r') - if not (self.dft_subgrp in ar): raise IOError, "convert_transport_input: No %s subgroup in hdf file found! Call convert_dft_input first." %self.dft_subgrp + if not (self.dft_subgrp in ar): + raise IOError, "convert_transport_input: No %s subgroup in hdf file found! Call convert_dft_input first." % self.dft_subgrp SP = ar[self.dft_subgrp]['SP'] SO = ar[self.dft_subgrp]['SO'] n_k = ar[self.dft_subgrp]['n_k'] @@ -571,20 +648,22 @@ class Wien2kConverter(ConverterTools): # velocities_k: velocity (momentum) matrix elements between all bands in band_window_optics # and each k-point. - if (SP == 0 or SO == 1): + if (SP == 0 or SO == 1): files = [self.pmat_file] elif SP == 1: - files = [self.pmat_file+'up', self.pmat_file+'dn'] - else: # SO and SP can't both be 1 + files = [self.pmat_file + 'up', self.pmat_file + 'dn'] + else: # SO and SP can't both be 1 assert 0, "convert_transport_input: Reading velocity file error! Check SP and SO!" velocities_k = [[] for f in files] band_window_optics = [] for isp, f in enumerate(files): - if not os.path.exists(f) : raise IOError, "convert_transport_input: File %s does not exist" %f - mpi.report("Reading input from %s..."%f) + if not os.path.exists(f): + raise IOError, "convert_transport_input: File %s does not exist" % f + mpi.report("Reading input from %s..." % f) - R = ConverterTools.read_fortran_file(self, f, {'D':'E','(':'',')':'',',':' '}) + R = ConverterTools.read_fortran_file( + self, f, {'D': 'E', '(': '', ')': '', ',': ' '}) band_window_optics_isp = [] for ik in xrange(n_k): R.next() @@ -592,26 +671,34 @@ class Wien2kConverter(ConverterTools): nu2 = int(R.next()) band_window_optics_isp.append((nu1, nu2)) n_bands = nu2 - nu1 + 1 - for _ in range(4): R.next() + for _ in range(4): + R.next() if n_bands <= 0: - velocity_xyz = numpy.zeros((1, 1, 3), dtype = complex) + velocity_xyz = numpy.zeros((1, 1, 3), dtype=complex) else: - velocity_xyz = numpy.zeros((n_bands, n_bands, 3), dtype = complex) + velocity_xyz = numpy.zeros( + (n_bands, n_bands, 3), dtype=complex) for nu_i in range(n_bands): for nu_j in range(nu_i, n_bands): for i in range(3): - velocity_xyz[nu_i][nu_j][i] = R.next() + R.next()*1j - if (nu_i != nu_j): velocity_xyz[nu_j][nu_i][i] = velocity_xyz[nu_i][nu_j][i].conjugate() + velocity_xyz[nu_i][nu_j][ + i] = R.next() + R.next() * 1j + if (nu_i != nu_j): + velocity_xyz[nu_j][nu_i][i] = velocity_xyz[ + nu_i][nu_j][i].conjugate() velocities_k[isp].append(velocity_xyz) band_window_optics.append(numpy.array(band_window_optics_isp)) - R.close() # Reading done! + R.close() # Reading done! # Put data to HDF5 file 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!!! + 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 = ['band_window_optics', 'velocities_k'] - for it in things_to_save: ar[self.transp_subgrp][it] = locals()[it] + for it in things_to_save: + ar[self.transp_subgrp][it] = locals()[it] del ar def convert_symmetry_input(self, orbits, symm_file, symm_subgrp, SO, SP): @@ -635,59 +722,70 @@ class Wien2kConverter(ConverterTools): """ - if not (mpi.is_master_node()): return - mpi.report("Reading input from %s..."%symm_file) + if not (mpi.is_master_node()): + return + mpi.report("Reading input from %s..." % symm_file) n_orbits = len(orbits) - R = ConverterTools.read_fortran_file(self,symm_file,self.fortran_to_replace) + R = ConverterTools.read_fortran_file( + self, symm_file, self.fortran_to_replace) try: n_symm = int(R.next()) # Number of symmetry operations n_atoms = int(R.next()) # number of atoms involved - perm = [ [int(R.next()) for i in range(n_atoms)] for j in range(n_symm) ] # list of permutations of the atoms - if SP: - time_inv = [ int(R.next()) for j in range(n_symm) ] # time inversion for SO coupling + perm = [[int(R.next()) for i in range(n_atoms)] + for j in range(n_symm)] # list of permutations of the atoms + if SP: + # time inversion for SO coupling + time_inv = [int(R.next()) for j in range(n_symm)] else: - time_inv = [ 0 for j in range(n_symm) ] + time_inv = [0 for j in range(n_symm)] # Now read matrices: - mat = [] + mat = [] for i_symm in range(n_symm): - - mat.append( [ numpy.zeros([orbits[orb]['dim'], orbits[orb]['dim']],numpy.complex_) for orb in range(n_orbits) ] ) + + mat.append([numpy.zeros([orbits[orb]['dim'], orbits[orb][ + 'dim']], numpy.complex_) for orb in range(n_orbits)]) for orb in range(n_orbits): for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): - mat[i_symm][orb][i,j] = R.next() # real part + # real part + mat[i_symm][orb][i, j] = R.next() for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): - mat[i_symm][orb][i,j] += 1j * R.next() # imaginary part + mat[i_symm][orb][i, j] += 1j * \ + R.next() # imaginary part - mat_tinv = [numpy.identity(orbits[orb]['dim'],numpy.complex_) + mat_tinv = [numpy.identity(orbits[orb]['dim'], numpy.complex_) for orb in range(n_orbits)] - if ((SO==0) and (SP==0)): - # here we need an additional time inversion operation, so read it: + if ((SO == 0) and (SP == 0)): + # here we need an additional time inversion operation, so read + # it: for orb in range(n_orbits): for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): - mat_tinv[orb][i,j] = R.next() # real part + # real part + mat_tinv[orb][i, j] = R.next() for i in range(orbits[orb]['dim']): for j in range(orbits[orb]['dim']): - mat_tinv[orb][i,j] += 1j * R.next() # imaginary part - + mat_tinv[orb][i, j] += 1j * \ + R.next() # imaginary part - - except StopIteration : # a more explicit error if the file is corrupted. + except StopIteration: # a more explicit error if the file is corrupted. raise "Wien2k_converter : reading file symm_file failed!" - + R.close() # Reading done! # Save it to the HDF: - ar = HDFArchive(self.hdf_file,'a') - if not (symm_subgrp in ar): ar.create_group(symm_subgrp) - things_to_save = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv'] - for it in things_to_save: ar[symm_subgrp][it] = locals()[it] + ar = HDFArchive(self.hdf_file, 'a') + if not (symm_subgrp in ar): + ar.create_group(symm_subgrp) + things_to_save = ['n_symm', 'n_atoms', 'perm', + 'orbits', 'SO', 'SP', 'time_inv', 'mat', 'mat_tinv'] + for it in things_to_save: + ar[symm_subgrp][it] = locals()[it] del ar diff --git a/python/sumk_dft.py b/python/sumk_dft.py index 629b58df..a500b478 100644 --- a/python/sumk_dft.py +++ b/python/sumk_dft.py @@ -1,5 +1,5 @@ - -################################################################################ + +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,7 +18,7 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## from types import * import numpy @@ -30,14 +30,14 @@ from symmetry import * from sets import Set from itertools import product + class SumkDFT: """This class provides a general SumK method for combining ab-initio code and pytriqs.""" - - def __init__(self, hdf_file, h_field = 0.0, use_dft_blocks = False, - dft_data = 'dft_input', symmcorr_data = 'dft_symmcorr_input', parproj_data = 'dft_parproj_input', - symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input', transp_data = 'dft_transp_input', - misc_data = 'dft_misc_input'): + def __init__(self, hdf_file, h_field=0.0, use_dft_blocks=False, + dft_data='dft_input', symmcorr_data='dft_symmcorr_input', parproj_data='dft_parproj_input', + symmpar_data='dft_symmpar_input', bands_data='dft_bands_input', transp_data='dft_transp_input', + misc_data='dft_misc_input'): r""" Initialises the class from data previously stored into an hdf5 archive. @@ -87,50 +87,63 @@ class SumkDFT: self.h_field = h_field # Read input from HDF: - things_to_read = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required', - 'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat', - 'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping', + things_to_read = ['energy_unit', 'n_k', 'k_dep_projection', 'SP', 'SO', 'charge_below', 'density_required', + 'symm_op', 'n_shells', 'shells', 'n_corr_shells', 'corr_shells', 'use_rotations', 'rot_mat', + 'rot_mat_time_inv', 'n_reps', 'dim_reps', 'T', 'n_orbitals', 'proj_mat', 'bz_weights', 'hopping', 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr'] - self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.dft_data, things_to_read = things_to_read) - if self.symm_op: self.symmcorr = Symmetry(hdf_file,subgroup=self.symmcorr_data) + self.subgroup_present, self.value_read = self.read_input_from_hdf( + subgrp=self.dft_data, things_to_read=things_to_read) + if self.symm_op: + self.symmcorr = Symmetry(hdf_file, subgroup=self.symmcorr_data) if self.SO and (abs(self.h_field) > 0.000001): self.h_field = 0.0 - mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!") + mpi.report( + "For SO, the external magnetic field is not implemented, setting it to 0!") - self.spin_block_names = [ ['up','down'], ['ud'] ] - self.n_spin_blocks = [2,1] - # Convert spin_block_names to indices -- if spin polarized, differentiate up and down blocks + self.spin_block_names = [['up', 'down'], ['ud']] + self.n_spin_blocks = [2, 1] + # Convert spin_block_names to indices -- if spin polarized, + # differentiate up and down blocks self.spin_names_to_ind = [{}, {}] - for iso in range(2): # SO = 0 or 1 + for iso in range(2): # SO = 0 or 1 for isp in range(self.n_spin_blocks[iso]): - self.spin_names_to_ind[iso][self.spin_block_names[iso][isp]] = isp * self.SP + self.spin_names_to_ind[iso][ + self.spin_block_names[iso][isp]] = isp * self.SP # GF structure used for the local things in the k sums - # Most general form allowing for all hybridisation, i.e. largest blocks possible - self.gf_struct_sumk = [ [ (sp, range( self.corr_shells[icrsh]['dim'])) for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']] ] - for icrsh in range(self.n_corr_shells) ] + # Most general form allowing for all hybridisation, i.e. largest + # blocks possible + self.gf_struct_sumk = [[(sp, range(self.corr_shells[icrsh]['dim'])) for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]] + for icrsh in range(self.n_corr_shells)] # First set a standard gf_struct solver: - self.gf_struct_solver = [ dict([ (sp, range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']) ) - for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']] ]) - for ish in range(self.n_inequiv_shells) ] - # Set standard (identity) maps from gf_struct_sumk <-> gf_struct_solver - self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ] - self.solver_to_sumk = [ {} for ish in range(self.n_inequiv_shells) ] - self.solver_to_sumk_block = [ {} for ish in range(self.n_inequiv_shells) ] + self.gf_struct_solver = [dict([(sp, range(self.corr_shells[self.inequiv_to_corr[ish]]['dim'])) + for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]]) + for ish in range(self.n_inequiv_shells)] + # Set standard (identity) maps from gf_struct_sumk <-> + # gf_struct_solver + self.sumk_to_solver = [{} for ish in range(self.n_inequiv_shells)] + self.solver_to_sumk = [{} for ish in range(self.n_inequiv_shells)] + self.solver_to_sumk_block = [{} + for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells): - for block,inner_list in self.gf_struct_sumk[self.inequiv_to_corr[ish]]: + for block, inner_list in self.gf_struct_sumk[self.inequiv_to_corr[ish]]: self.solver_to_sumk_block[ish][block] = block for inner in inner_list: - self.sumk_to_solver[ish][(block,inner)] = (block,inner) - self.solver_to_sumk[ish][(block,inner)] = (block,inner) - self.deg_shells = [ [] for ish in range(self.n_inequiv_shells) ] # assume no shells are degenerate + self.sumk_to_solver[ish][ + (block, inner)] = (block, inner) + self.solver_to_sumk[ish][ + (block, inner)] = (block, inner) + # assume no shells are degenerate + self.deg_shells = [[] for ish in range(self.n_inequiv_shells)] - self.chemical_potential = 0.0 # initialise mu - self.init_dc() # initialise the double counting + self.chemical_potential = 0.0 # initialise mu + self.init_dc() # initialise the double counting - # Analyse the block structure and determine the smallest gf_struct blocks and maps, if desired - if use_dft_blocks: self.analyse_block_structure() + # Analyse the block structure and determine the smallest gf_struct + # blocks and maps, if desired + if use_dft_blocks: + self.analyse_block_structure() ################ # hdf5 FUNCTIONS @@ -157,36 +170,39 @@ class SumkDFT: """ value_read = True - # initialise variables on all nodes to ensure mpi broadcast works at the end - for it in things_to_read: setattr(self,it,0) + # initialise variables on all nodes to ensure mpi broadcast works at + # the end + for it in things_to_read: + setattr(self, it, 0) subgroup_present = 0 if mpi.is_master_node(): - ar = HDFArchive(self.hdf_file,'r') + ar = HDFArchive(self.hdf_file, 'r') if subgrp in ar: subgroup_present = True # first read the necessary things: for it in things_to_read: if it in ar[subgrp]: - setattr(self,it,ar[subgrp][it]) + setattr(self, it, ar[subgrp][it]) else: - mpi.report("Loading %s failed!"%it) + mpi.report("Loading %s failed!" % it) value_read = False else: - if (len(things_to_read) != 0): mpi.report("Loading failed: No %s subgroup in hdf5!"%subgrp) + if (len(things_to_read) != 0): + mpi.report( + "Loading failed: No %s subgroup in hdf5!" % subgrp) subgroup_present = False value_read = False del ar # now do the broadcasting: - for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it))) + for it in things_to_read: + setattr(self, it, mpi.bcast(getattr(self, it))) subgroup_present = mpi.bcast(subgroup_present) value_read = mpi.bcast(value_read) return subgroup_present, value_read - def save(self, things_to_save, subgrp='user_data'): - r""" Saves data from a list into the HDF file. Prints a warning if a requested data is not found in SumkDFT object. @@ -198,17 +214,18 @@ class SumkDFT: Name of hdf5 file subgroup in which the data are to be stored. """ - if not (mpi.is_master_node()): return # do nothing on nodes - ar = HDFArchive(self.hdf_file,'a') - if not subgrp in ar: ar.create_group(subgrp) - for it in things_to_save: + if not (mpi.is_master_node()): + return # do nothing on nodes + ar = HDFArchive(self.hdf_file, 'a') + if not subgrp in ar: + ar.create_group(subgrp) + for it in things_to_save: try: - ar[subgrp][it] = getattr(self,it) + ar[subgrp][it] = getattr(self, it) except: - mpi.report("%s not found, and so not saved."%it) + mpi.report("%s not found, and so not saved." % it) del ar - def load(self, things_to_load, subgrp='user_data'): r""" Loads user data from the HDF file. Raises an exeption if a requested dataset is not found. @@ -226,15 +243,17 @@ class SumkDFT: A list containing data read from hdf5. """ - if not (mpi.is_master_node()): return # do nothing on nodes - ar = HDFArchive(self.hdf_file,'r') - if not subgrp in ar: mpi.report("Loading %s failed!"%subgrp) + if not (mpi.is_master_node()): + return # do nothing on nodes + ar = HDFArchive(self.hdf_file, 'r') + if not subgrp in ar: + mpi.report("Loading %s failed!" % subgrp) list_to_return = [] - for it in things_to_load: + for it in things_to_load: try: list_to_return.append(ar[subgrp][it]) except: - raise ValueError, "load: %s not found, and so not loaded."%it + raise ValueError, "load: %s not found, and so not loaded." % it del ar return list_to_return @@ -242,7 +261,7 @@ class SumkDFT: # CORE FUNCTIONS ################ - def downfold(self,ik,ish,bname,gf_to_downfold,gf_inp,shells='corr',ir=None): + def downfold(self, ik, ish, bname, gf_to_downfold, gf_inp, shells='corr', ir=None): r""" Downfolds a block of the Green's function for a given shell and k-point using the corresponding projector matrices. @@ -269,30 +288,32 @@ class SumkDFT: ir : integer, optional Index of equivalent site in the non-correlated shell 'ish', only used if shells='all'. - + Returns ------- gf_downfolded : Gf Downfolded block of the lattice Green's function. """ - + gf_downfolded = gf_inp.copy() - isp = self.spin_names_to_ind[self.SO][bname] # get spin index for proj. matrices - n_orb = self.n_orbitals[ik,isp] + # get spin index for proj. matrices + isp = self.spin_names_to_ind[self.SO][bname] + n_orb = self.n_orbitals[ik, isp] if shells == 'corr': dim = self.corr_shells[ish]['dim'] - projmat = self.proj_mat[ik,isp,ish,0:dim,0:n_orb] + projmat = self.proj_mat[ik, isp, ish, 0:dim, 0:n_orb] elif shells == 'all': - if ir is None: raise ValueError, "downfold: provide ir if treating all shells." + if ir is None: + raise ValueError, "downfold: provide ir if treating all shells." dim = self.shells[ish]['dim'] - projmat = self.proj_mat_all[ik,isp,ish,ir,0:dim,0:n_orb] - - gf_downfolded.from_L_G_R(projmat,gf_to_downfold,projmat.conjugate().transpose()) - + projmat = self.proj_mat_all[ik, isp, ish, ir, 0:dim, 0:n_orb] + + gf_downfolded.from_L_G_R( + projmat, gf_to_downfold, projmat.conjugate().transpose()) + return gf_downfolded - - def upfold(self,ik,ish,bname,gf_to_upfold,gf_inp,shells='corr',ir=None): + def upfold(self, ik, ish, bname, gf_to_upfold, gf_inp, shells='corr', ir=None): r""" Upfolds a block of the Green's function for a given shell and k-point using the corresponding projector matrices. @@ -319,30 +340,32 @@ class SumkDFT: ir : integer, optional Index of equivalent site in the non-correlated shell 'ish', only used if shells='all'. - + Returns ------- gf_upfolded : Gf Upfolded block of the lattice Green's function. """ - + gf_upfolded = gf_inp.copy() - isp = self.spin_names_to_ind[self.SO][bname] # get spin index for proj. matrices - n_orb = self.n_orbitals[ik,isp] + # get spin index for proj. matrices + isp = self.spin_names_to_ind[self.SO][bname] + n_orb = self.n_orbitals[ik, isp] if shells == 'corr': dim = self.corr_shells[ish]['dim'] - projmat = self.proj_mat[ik,isp,ish,0:dim,0:n_orb] + projmat = self.proj_mat[ik, isp, ish, 0:dim, 0:n_orb] elif shells == 'all': - if ir is None: raise ValueError, "upfold: provide ir if treating all shells." + if ir is None: + raise ValueError, "upfold: provide ir if treating all shells." dim = self.shells[ish]['dim'] - projmat = self.proj_mat_all[ik,isp,ish,ir,0:dim,0:n_orb] - - gf_upfolded.from_L_G_R(projmat.conjugate().transpose(),gf_to_upfold,projmat) - + projmat = self.proj_mat_all[ik, isp, ish, ir, 0:dim, 0:n_orb] + + gf_upfolded.from_L_G_R( + projmat.conjugate().transpose(), gf_to_upfold, projmat) + return gf_upfolded - - def rotloc(self,ish,gf_to_rotate,direction,shells='corr'): + def rotloc(self, ish, gf_to_rotate, direction, shells='corr'): r""" Rotates a block of the local Green's function from the local frame to the global frame and vice versa. @@ -373,7 +396,8 @@ class SumkDFT: Rotated block of the local Green's function. """ - assert ((direction == 'toLocal') or (direction == 'toGlobal')),"rotloc: Give direction 'toLocal' or 'toGlobal'." + assert ((direction == 'toLocal') or (direction == 'toGlobal') + ), "rotloc: Give direction 'toLocal' or 'toGlobal'." gf_rotated = gf_to_rotate.copy() if shells == 'corr': rot_mat_time_inv = self.rot_mat_time_inv @@ -386,21 +410,24 @@ class SumkDFT: if (rot_mat_time_inv[ish] == 1) and self.SO: gf_rotated << gf_rotated.transpose() - gf_rotated.from_L_G_R(rot_mat[ish].conjugate(),gf_rotated,rot_mat[ish].transpose()) + gf_rotated.from_L_G_R(rot_mat[ish].conjugate( + ), gf_rotated, rot_mat[ish].transpose()) else: - gf_rotated.from_L_G_R(rot_mat[ish],gf_rotated,rot_mat[ish].conjugate().transpose()) + gf_rotated.from_L_G_R(rot_mat[ish], gf_rotated, rot_mat[ + ish].conjugate().transpose()) elif direction == 'toLocal': if (rot_mat_time_inv[ish] == 1) and self.SO: gf_rotated << gf_rotated.transpose() - gf_rotated.from_L_G_R(rot_mat[ish].transpose(),gf_rotated,rot_mat[ish].conjugate()) + gf_rotated.from_L_G_R(rot_mat[ish].transpose( + ), gf_rotated, rot_mat[ish].conjugate()) else: - gf_rotated.from_L_G_R(rot_mat[ish].conjugate().transpose(),gf_rotated,rot_mat[ish]) + gf_rotated.from_L_G_R(rot_mat[ish].conjugate( + ).transpose(), gf_rotated, rot_mat[ish]) return gf_rotated - def lattice_gf(self, ik, mu=None, iw_or_w="iw", beta=40, broadening=None, mesh=None, with_Sigma=True, with_dc=True): r""" Calculates the lattice Green function for a given k-point from the DFT Hamiltonian and the self energy. @@ -431,88 +458,109 @@ class SumkDFT: If with_Sigma=True but self.Sigmaimp_(w/iw) is not present, with_Sigma is reset to False. with_dc : boolean, optional if True and with_Sigma=True, the dc correction is substracted from the self-energy before it is included into GF. - + Returns ------- G_latt : BlockGf Lattice Green's function. - + """ - if mu is None: mu = self.chemical_potential + if mu is None: + mu = self.chemical_potential ntoi = self.spin_names_to_ind[self.SO] spn = self.spin_block_names[self.SO] - if (iw_or_w != "iw") and (iw_or_w != "w"): raise ValueError, "lattice_gf: Implemented only for Re/Im frequency functions." - if not hasattr(self,"Sigma_imp_"+iw_or_w): with_Sigma = False + if (iw_or_w != "iw") and (iw_or_w != "w"): + raise ValueError, "lattice_gf: Implemented only for Re/Im frequency functions." + if not hasattr(self, "Sigma_imp_" + iw_or_w): + with_Sigma = False if broadening is None: if mesh is None: broadening = 0.01 - else: # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points - broadening = 2.0 * ( (mesh[1]-mesh[0])/(mesh[2]-1) ) + else: # broadening = 2 * \Delta omega, where \Delta omega is the spacing of omega points + broadening = 2.0 * ((mesh[1] - mesh[0]) / (mesh[2] - 1)) # Are we including Sigma? if with_Sigma: - Sigma_imp = getattr(self,"Sigma_imp_"+iw_or_w) + Sigma_imp = getattr(self, "Sigma_imp_" + iw_or_w) sigma_minus_dc = [s.copy() for s in Sigma_imp] - if with_dc: sigma_minus_dc = self.add_dc(iw_or_w) + if with_dc: + sigma_minus_dc = self.add_dc(iw_or_w) if iw_or_w == "iw": - beta = Sigma_imp[0].mesh.beta # override beta if Sigma_iw is present + # override beta if Sigma_iw is present + beta = Sigma_imp[0].mesh.beta mesh = Sigma_imp[0].mesh elif iw_or_w == "w": mesh = Sigma_imp[0].mesh else: if iw_or_w == "iw": - if beta is None: raise ValueError, "lattice_gf: Give the beta for the lattice GfReFreq." - mesh = MeshImFreq(beta=beta, S='Fermion', n_max=1025) # Default number of Matsubara frequencies + if beta is None: + raise ValueError, "lattice_gf: Give the beta for the lattice GfReFreq." + # Default number of Matsubara frequencies + mesh = MeshImFreq(beta=beta, S='Fermion', n_max=1025) elif iw_or_w == "w": - if mesh is None: raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq." - mesh = MeshReFreq(mesh[0],mesh[1],mesh[2]) + if mesh is None: + raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq." + mesh = MeshReFreq(mesh[0], mesh[1], mesh[2]) # Check if G_latt is present set_up_G_latt = False # Assume not - if not hasattr(self,"G_latt_"+iw_or_w): - set_up_G_latt = True # Need to create G_latt_(i)w + if not hasattr(self, "G_latt_" + iw_or_w): + # Need to create G_latt_(i)w + set_up_G_latt = True else: # Check that existing GF is consistent - G_latt = getattr(self,"G_latt_"+iw_or_w) - GFsize = [ gf.N1 for bname,gf in G_latt] - unchangedsize = all( [ self.n_orbitals[ik,ntoi[spn[isp]]]==GFsize[isp] for isp in range(self.n_spin_blocks[self.SO]) ] ) - if not unchangedsize: set_up_G_latt = True - if (iw_or_w == "iw") and (self.G_latt_iw.mesh.beta != beta): set_up_G_latt = True # additional check for ImFreq + G_latt = getattr(self, "G_latt_" + iw_or_w) + GFsize = [gf.N1 for bname, gf in G_latt] + unchangedsize = all([self.n_orbitals[ik, ntoi[spn[isp]]] == GFsize[ + isp] for isp in range(self.n_spin_blocks[self.SO])]) + if not unchangedsize: + set_up_G_latt = True + if (iw_or_w == "iw") and (self.G_latt_iw.mesh.beta != beta): + set_up_G_latt = True # additional check for ImFreq # Set up G_latt if set_up_G_latt: - block_structure = [ range(self.n_orbitals[ik,ntoi[sp]]) for sp in spn ] - gf_struct = [ (spn[isp], block_structure[isp]) for isp in range(self.n_spin_blocks[self.SO]) ] - block_ind_list = [block for block,inner in gf_struct] + block_structure = [ + range(self.n_orbitals[ik, ntoi[sp]]) for sp in spn] + gf_struct = [(spn[isp], block_structure[isp]) + for isp in range(self.n_spin_blocks[self.SO])] + block_ind_list = [block for block, inner in gf_struct] if iw_or_w == "iw": - glist = lambda : [ GfImFreq(indices=inner,mesh=mesh) for block,inner in gf_struct ] + glist = lambda: [GfImFreq(indices=inner, mesh=mesh) + for block, inner in gf_struct] elif iw_or_w == "w": - glist = lambda : [ GfReFreq(indices=inner,mesh=mesh) for block,inner in gf_struct ] - G_latt = BlockGf(name_list = block_ind_list, block_list = glist(), make_copies = False) + glist = lambda: [GfReFreq(indices=inner, mesh=mesh) + for block, inner in gf_struct] + G_latt = BlockGf(name_list=block_ind_list, + block_list=glist(), make_copies=False) G_latt.zero() if iw_or_w == "iw": G_latt << iOmega_n elif iw_or_w == "w": - G_latt << Omega + 1j*broadening + G_latt << Omega + 1j * broadening - idmat = [numpy.identity(self.n_orbitals[ik,ntoi[sp]],numpy.complex_) for sp in spn] + idmat = [numpy.identity( + self.n_orbitals[ik, ntoi[sp]], numpy.complex_) for sp in spn] M = copy.deepcopy(idmat) for ibl in range(self.n_spin_blocks[self.SO]): ind = ntoi[spn[ibl]] - n_orb = self.n_orbitals[ik,ind] - M[ibl] = self.hopping[ik,ind,0:n_orb,0:n_orb] - (idmat[ibl]*mu) - (idmat[ibl] * self.h_field * (1-2*ibl)) + n_orb = self.n_orbitals[ik, ind] + M[ibl] = self.hopping[ik, ind, 0:n_orb, 0:n_orb] - \ + (idmat[ibl] * mu) - (idmat[ibl] * self.h_field * (1 - 2 * ibl)) G_latt -= M if with_Sigma: for icrsh in range(self.n_corr_shells): - for bname,gf in G_latt: gf -= self.upfold(ik,icrsh,bname,sigma_minus_dc[icrsh][bname],gf) + for bname, gf in G_latt: + gf -= self.upfold(ik, icrsh, bname, + sigma_minus_dc[icrsh][bname], gf) G_latt.invert() - setattr(self,"G_latt_"+iw_or_w,G_latt) + setattr(self, "G_latt_" + iw_or_w, G_latt) return G_latt - def set_Sigma(self,Sigma_imp): + def set_Sigma(self, Sigma_imp): self.put_Sigma(Sigma_imp) def put_Sigma(self, Sigma_imp): @@ -527,39 +575,46 @@ class SumkDFT: The self-energies can be of the real or imaginary-frequency type. """ - assert isinstance(Sigma_imp,list), "put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!" - assert len(Sigma_imp) == self.n_inequiv_shells, "put_Sigma: give exactly one Sigma for each inequivalent corr. shell!" + assert isinstance( + Sigma_imp, list), "put_Sigma: Sigma_imp has to be a list of Sigmas for the correlated shells, even if it is of length 1!" + assert len( + Sigma_imp) == self.n_inequiv_shells, "put_Sigma: give exactly one Sigma for each inequivalent corr. shell!" # init self.Sigma_imp_(i)w: - if all(type(gf) == GfImFreq for bname,gf in Sigma_imp[0]): + if all(type(gf) == GfImFreq for bname, gf in Sigma_imp[0]): # Imaginary frequency Sigma: - self.Sigma_imp_iw = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = Sigma_imp[0].mesh)) - for block,inner in self.gf_struct_sumk[icrsh] ], make_copies = False) - for icrsh in range(self.n_corr_shells) ] + self.Sigma_imp_iw = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, mesh=Sigma_imp[0].mesh)) + for block, inner in self.gf_struct_sumk[icrsh]], make_copies=False) + for icrsh in range(self.n_corr_shells)] SK_Sigma_imp = self.Sigma_imp_iw - elif all(type(gf) == GfReFreq for bname,gf in Sigma_imp[0]): + elif all(type(gf) == GfReFreq for bname, gf in Sigma_imp[0]): # Real frequency Sigma: - self.Sigma_imp_w = [ BlockGf( name_block_generator = [ (block,GfReFreq(indices = inner, mesh = Sigma_imp[0].mesh)) - for block,inner in self.gf_struct_sumk[icrsh] ], make_copies = False) - for icrsh in range(self.n_corr_shells) ] + self.Sigma_imp_w = [BlockGf(name_block_generator=[(block, GfReFreq(indices=inner, mesh=Sigma_imp[0].mesh)) + for block, inner in self.gf_struct_sumk[icrsh]], make_copies=False) + for icrsh in range(self.n_corr_shells)] SK_Sigma_imp = self.Sigma_imp_w else: raise ValueError, "put_Sigma: This type of Sigma is not handled." # transform the CTQMC blocks to the full matrix: for icrsh in range(self.n_corr_shells): - ish = self.corr_to_inequiv[icrsh] # ish is the index of the inequivalent shell corresponding to icrsh - for block,inner in self.gf_struct_solver[ish].iteritems(): + # ish is the index of the inequivalent shell corresponding to icrsh + ish = self.corr_to_inequiv[icrsh] + for block, inner in self.gf_struct_solver[ish].iteritems(): for ind1 in inner: for ind2 in inner: - block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] - block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)] - SK_Sigma_imp[icrsh][block_sumk][ind1_sumk,ind2_sumk] << Sigma_imp[ish][block][ind1,ind2] + block_sumk, ind1_sumk = self.solver_to_sumk[ + ish][(block, ind1)] + block_sumk, ind2_sumk = self.solver_to_sumk[ + ish][(block, ind2)] + SK_Sigma_imp[icrsh][block_sumk][ + ind1_sumk, ind2_sumk] << Sigma_imp[ish][block][ind1, ind2] # rotation from local to global coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bname,gf in SK_Sigma_imp[icrsh]: gf << self.rotloc(icrsh,gf,direction='toGlobal') + for bname, gf in SK_Sigma_imp[icrsh]: + gf << self.rotloc(icrsh, gf, direction='toGlobal') def extract_G_loc(self, mu=None, iw_or_w='iw', with_Sigma=True, with_dc=True, broadening=None): r""" @@ -583,65 +638,80 @@ class SumkDFT: G_loc_inequiv : list of BlockGf (Green's function) objects List of the local Green's functions for all inequivalent correlated shells, rotated into the corresponding local frames. - + """ - if mu is None: mu = self.chemical_potential - - if iw_or_w == "iw": - G_loc = [ self.Sigma_imp_iw[icrsh].copy() for icrsh in range(self.n_corr_shells) ] # this list will be returned + if mu is None: + mu = self.chemical_potential + + if iw_or_w == "iw": + G_loc = [self.Sigma_imp_iw[icrsh].copy() for icrsh in range( + self.n_corr_shells)] # this list will be returned beta = G_loc[0].mesh.beta - G_loc_inequiv = [ BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = G_loc[0].mesh)) for block,inner in self.gf_struct_solver[ish].iteritems() ], - make_copies = False) for ish in range(self.n_inequiv_shells) ] - elif iw_or_w == "w": - G_loc = [ self.Sigma_imp_w[icrsh].copy() for icrsh in range(self.n_corr_shells) ] # this list will be returned + G_loc_inequiv = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, mesh=G_loc[0].mesh)) for block, inner in self.gf_struct_solver[ish].iteritems()], + make_copies=False) for ish in range(self.n_inequiv_shells)] + elif iw_or_w == "w": + G_loc = [self.Sigma_imp_w[icrsh].copy() for icrsh in range( + self.n_corr_shells)] # this list will be returned mesh = G_loc[0].mesh - G_loc_inequiv = [ BlockGf( name_block_generator = [ (block,GfReFreq(indices = inner, mesh = mesh)) for block,inner in self.gf_struct_solver[ish].iteritems() ], - make_copies = False) for ish in range(self.n_inequiv_shells) ] - - for icrsh in range(self.n_corr_shells): G_loc[icrsh].zero() # initialize to zero + G_loc_inequiv = [BlockGf(name_block_generator=[(block, GfReFreq(indices=inner, mesh=mesh)) for block, inner in self.gf_struct_solver[ish].iteritems()], + make_copies=False) for ish in range(self.n_inequiv_shells)] + + for icrsh in range(self.n_corr_shells): + G_loc[icrsh].zero() # initialize to zero ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): if iw_or_w == 'iw': - G_latt = self.lattice_gf(ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, beta=beta) + G_latt = self.lattice_gf( + ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, beta=beta) elif iw_or_w == 'w': - G_latt = self.lattice_gf(ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening, mesh=mesh) + G_latt = self.lattice_gf( + ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening, mesh=mesh) G_latt *= self.bz_weights[ik] for icrsh in range(self.n_corr_shells): - tmp = G_loc[icrsh].copy() # init temporary storage - for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,G_latt[bname],gf) + # init temporary storage + tmp = G_loc[icrsh].copy() + for bname, gf in tmp: + tmp[bname] << self.downfold( + ik, icrsh, bname, G_latt[bname], gf) G_loc[icrsh] += tmp # Collect data from mpi - for icrsh in range(self.n_corr_shells): - G_loc[icrsh] << mpi.all_reduce(mpi.world, G_loc[icrsh], lambda x,y : x+y) + for icrsh in range(self.n_corr_shells): + G_loc[icrsh] << mpi.all_reduce( + mpi.world, G_loc[icrsh], lambda x, y: x + y) mpi.barrier() # G_loc[:] is now the sum over k projected to the local orbitals. # here comes the symmetrisation, if needed: - if self.symm_op != 0: G_loc = self.symmcorr.symmetrize(G_loc) + if self.symm_op != 0: + G_loc = self.symmcorr.symmetrize(G_loc) # G_loc is rotated to the local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bname,gf in G_loc[icrsh]: G_loc[icrsh][bname] << self.rotloc(icrsh,gf,direction='toLocal') + for bname, gf in G_loc[icrsh]: + G_loc[icrsh][bname] << self.rotloc( + icrsh, gf, direction='toLocal') # transform to CTQMC blocks: for ish in range(self.n_inequiv_shells): - for block,inner in self.gf_struct_solver[ish].iteritems(): + for block, inner in self.gf_struct_solver[ish].iteritems(): for ind1 in inner: for ind2 in inner: - block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] - block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)] - G_loc_inequiv[ish][block][ind1,ind2] << G_loc[self.inequiv_to_corr[ish]][block_sumk][ind1_sumk,ind2_sumk] + block_sumk, ind1_sumk = self.solver_to_sumk[ + ish][(block, ind1)] + block_sumk, ind2_sumk = self.solver_to_sumk[ + ish][(block, ind2)] + G_loc_inequiv[ish][block][ind1, ind2] << G_loc[ + self.inequiv_to_corr[ish]][block_sumk][ind1_sumk, ind2_sumk] # return only the inequivalent shells: return G_loc_inequiv - - def analyse_block_structure(self, threshold = 0.00001, include_shells = None, dm = None): + def analyse_block_structure(self, threshold=0.00001, include_shells=None, dm=None): r""" Determines the block structure of local Green's functions by analysing the structure of the corresponding density matrices. The resulting block structures for correlated shells @@ -661,86 +731,102 @@ class SumkDFT: If not provided, dm will be calculated from the DFT Hamiltonian by a simple-point BZ integration. """ - self.gf_struct_solver = [ {} for ish in range(self.n_inequiv_shells) ] - self.sumk_to_solver = [ {} for ish in range(self.n_inequiv_shells) ] - self.solver_to_sumk = [ {} for ish in range(self.n_inequiv_shells) ] - self.solver_to_sumk_block = [ {} for ish in range(self.n_inequiv_shells) ] + self.gf_struct_solver = [{} for ish in range(self.n_inequiv_shells)] + self.sumk_to_solver = [{} for ish in range(self.n_inequiv_shells)] + self.solver_to_sumk = [{} for ish in range(self.n_inequiv_shells)] + self.solver_to_sumk_block = [{} + for ish in range(self.n_inequiv_shells)] - if dm is None: dm = self.density_matrix(method = 'using_point_integration') - dens_mat = [ dm[self.inequiv_to_corr[ish]] for ish in range(self.n_inequiv_shells) ] + if dm is None: + dm = self.density_matrix(method='using_point_integration') + dens_mat = [dm[self.inequiv_to_corr[ish]] + for ish in range(self.n_inequiv_shells)] - if include_shells is None: include_shells = range(self.n_inequiv_shells) + if include_shells is None: + include_shells = range(self.n_inequiv_shells) for ish in include_shells: - for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: n_orb = self.corr_shells[self.inequiv_to_corr[ish]]['dim'] - dmbool = (abs(dens_mat[ish][sp]) > threshold) # gives an index list of entries larger that threshold + # gives an index list of entries larger that threshold + dmbool = (abs(dens_mat[ish][sp]) > threshold) - # Determine off-diagonal entries in upper triangular part of density matrix + # Determine off-diagonal entries in upper triangular part of + # density matrix offdiag = Set([]) for i in range(n_orb): - for j in range(i+1,n_orb): - if dmbool[i,j]: offdiag.add((i,j)) + for j in range(i + 1, n_orb): + if dmbool[i, j]: + offdiag.add((i, j)) # Determine the number of non-hybridising blocks in the gf - blocs = [ [i] for i in range(n_orb) ] + blocs = [[i] for i in range(n_orb)] while len(offdiag) != 0: pair = offdiag.pop() - for b1,b2 in product(blocs,blocs): - if (pair[0] in b1) and (pair[1] in b2): - if blocs.index(b1) != blocs.index(b2): # In separate blocks? - b1.extend(blocs.pop(blocs.index(b2))) # Merge two blocks - break # Move on to next pair in offdiag + for b1, b2 in product(blocs, blocs): + if (pair[0] in b1) and (pair[1] in b2): + if blocs.index(b1) != blocs.index(b2): # In separate blocks? + # Merge two blocks + b1.extend(blocs.pop(blocs.index(b2))) + break # Move on to next pair in offdiag # Set the gf_struct for the solver accordingly num_blocs = len(blocs) for i in range(num_blocs): blocs[i].sort() - self.gf_struct_solver[ish].update( [('%s_%s'%(sp,i),range(len(blocs[i])))] ) + self.gf_struct_solver[ish].update( + [('%s_%s' % (sp, i), range(len(blocs[i])))]) # Construct sumk_to_solver taking (sumk_block, sumk_index) --> (solver_block, solver_inner) - # and solver_to_sumk taking (solver_block, solver_inner) --> (sumk_block, sumk_index) + # and solver_to_sumk taking (solver_block, solver_inner) --> + # (sumk_block, sumk_index) for i in range(num_blocs): for j in range(len(blocs[i])): block_sumk = sp inner_sumk = blocs[i][j] - block_solv = '%s_%s'%(sp,i) + block_solv = '%s_%s' % (sp, i) inner_solv = j - self.sumk_to_solver[ish][(block_sumk,inner_sumk)] = (block_solv,inner_solv) - self.solver_to_sumk[ish][(block_solv,inner_solv)] = (block_sumk,inner_sumk) + self.sumk_to_solver[ish][(block_sumk, inner_sumk)] = ( + block_solv, inner_solv) + self.solver_to_sumk[ish][(block_solv, inner_solv)] = ( + block_sumk, inner_sumk) self.solver_to_sumk_block[ish][block_solv] = block_sumk # Now calculate degeneracies of orbitals dm = {} - for block,inner in self.gf_struct_solver[ish].iteritems(): + for block, inner in self.gf_struct_solver[ish].iteritems(): # get dm for the blocks: - dm[block] = numpy.zeros([len(inner),len(inner)],numpy.complex_) + dm[block] = numpy.zeros( + [len(inner), len(inner)], numpy.complex_) for ind1 in inner: for ind2 in inner: - block_sumk,ind1_sumk = self.solver_to_sumk[ish][(block,ind1)] - block_sumk,ind2_sumk = self.solver_to_sumk[ish][(block,ind2)] - dm[block][ind1,ind2] = dens_mat[ish][block_sumk][ind1_sumk,ind2_sumk] + block_sumk, ind1_sumk = self.solver_to_sumk[ + ish][(block, ind1)] + block_sumk, ind2_sumk = self.solver_to_sumk[ + ish][(block, ind2)] + dm[block][ind1, ind2] = dens_mat[ish][ + block_sumk][ind1_sumk, ind2_sumk] for block1 in self.gf_struct_solver[ish].iterkeys(): - for block2 in self.gf_struct_solver[ish].iterkeys(): + for block2 in self.gf_struct_solver[ish].iterkeys(): if dm[block1].shape == dm[block2].shape: - if ( (abs(dm[block1] - dm[block2]) < threshold).all() ) and (block1 != block2): + if ((abs(dm[block1] - dm[block2]) < threshold).all()) and (block1 != block2): ind1 = -1 ind2 = -2 # check if it was already there: - for n,ind in enumerate(self.deg_shells[ish]): - if block1 in ind: ind1 = n - if block2 in ind: ind2 = n + for n, ind in enumerate(self.deg_shells[ish]): + if block1 in ind: + ind1 = n + if block2 in ind: + ind2 = n if (ind1 < 0) and (ind2 >= 0): self.deg_shells[ish][ind2].append(block1) elif (ind1 >= 0) and (ind2 < 0): self.deg_shells[ish][ind1].append(block2) elif (ind1 < 0) and (ind2 < 0): - self.deg_shells[ish].append([block1,block2]) + self.deg_shells[ish].append([block1, block2]) - - def density_matrix(self, method = 'using_gf', beta = 40.0): + def density_matrix(self, method='using_gf', beta=40.0): """Calculate density matrices in one of two ways. Parameters @@ -760,17 +846,19 @@ class SumkDFT: dens_mat : list of dicts Density matrix for each spin in each correlated shell. """ - dens_mat = [ {} for icrsh in range(self.n_corr_shells)] + dens_mat = [{} for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells): for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]: - dens_mat[icrsh][sp] = numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']], numpy.complex_) + dens_mat[icrsh][sp] = numpy.zeros( + [self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_) ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): if method == "using_gf": - G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw", beta = beta) + G_latt_iw = self.lattice_gf( + ik=ik, mu=self.chemical_potential, iw_or_w="iw", beta=beta) G_latt_iw *= self.bz_weights[ik] dm = G_latt_iw.density() MMat = [dm[sp] for sp in self.spin_block_names[self.SO]] @@ -779,55 +867,62 @@ class SumkDFT: ntoi = self.spin_names_to_ind[self.SO] spn = self.spin_block_names[self.SO] - unchangedsize = all( [self.n_orbitals[ik,ntoi[sp]] == self.n_orbitals[0,ntoi[sp]] for sp in spn] ) + unchangedsize = all( + [self.n_orbitals[ik, ntoi[sp]] == self.n_orbitals[0, ntoi[sp]] for sp in spn]) if unchangedsize: - dim = self.n_orbitals[0,ntoi[sp]] + dim = self.n_orbitals[0, ntoi[sp]] else: - dim = self.n_orbitals[ik,ntoi[sp]] - MMat = [numpy.zeros( [dim,dim], numpy.complex_) for sp in spn] + dim = self.n_orbitals[ik, ntoi[sp]] + MMat = [numpy.zeros([dim, dim], numpy.complex_) for sp in spn] for isp, sp in enumerate(spn): ind = ntoi[sp] - for inu in range(self.n_orbitals[ik,ind]): - if (self.hopping[ik,ind,inu,inu] - self.h_field*(1-2*isp)) < 0.0: # only works for diagonal hopping matrix (true in wien2k) - MMat[isp][inu,inu] = 1.0 + for inu in range(self.n_orbitals[ik, ind]): + # only works for diagonal hopping matrix (true in + # wien2k) + if (self.hopping[ik, ind, inu, inu] - self.h_field * (1 - 2 * isp)) < 0.0: + MMat[isp][inu, inu] = 1.0 else: - MMat[isp][inu,inu] = 0.0 + MMat[isp][inu, inu] = 0.0 - else: raise ValueError, "density_matrix: the method '%s' is not supported."%method + else: + raise ValueError, "density_matrix: the method '%s' is not supported." % method for icrsh in range(self.n_corr_shells): for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]): - ind = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp] + ind = self.spin_names_to_ind[ + self.corr_shells[icrsh]['SO']][sp] dim = self.corr_shells[icrsh]['dim'] - n_orb = self.n_orbitals[ik,ind] - projmat = self.proj_mat[ik,ind,icrsh,0:dim,0:n_orb] + n_orb = self.n_orbitals[ik, ind] + projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb] if method == "using_gf": - dens_mat[icrsh][sp] += numpy.dot( numpy.dot(projmat,MMat[isp]), - projmat.transpose().conjugate() ) + dens_mat[icrsh][sp] += numpy.dot(numpy.dot(projmat, MMat[isp]), + projmat.transpose().conjugate()) elif method == "using_point_integration": - dens_mat[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat[isp]) , - projmat.transpose().conjugate() ) + dens_mat[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.dot(projmat, MMat[isp]), + projmat.transpose().conjugate()) # get data from nodes: for icrsh in range(self.n_corr_shells): for sp in dens_mat[icrsh]: - dens_mat[icrsh][sp] = mpi.all_reduce(mpi.world, dens_mat[icrsh][sp], lambda x,y : x+y) + dens_mat[icrsh][sp] = mpi.all_reduce( + mpi.world, dens_mat[icrsh][sp], lambda x, y: x + y) mpi.barrier() - if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat) + if self.symm_op != 0: + dens_mat = self.symmcorr.symmetrize(dens_mat) # Rotate to local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): for sp in dens_mat[icrsh]: - if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate() - dens_mat[icrsh][sp] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh][sp]), - self.rot_mat[icrsh] ) + if self.rot_mat_time_inv[icrsh] == 1: + dens_mat[icrsh][sp] = dens_mat[icrsh][sp].conjugate() + dens_mat[icrsh][sp] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh][sp]), + self.rot_mat[icrsh]) return dens_mat - # For simple dft input, get crystal field splittings. def eff_atomic_levels(self): r""" @@ -854,50 +949,58 @@ class SumkDFT: """ # define matrices for inequivalent shells: - eff_atlevels = [ {} for ish in range(self.n_inequiv_shells) ] + eff_atlevels = [{} for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells): for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: - eff_atlevels[ish][sp] = numpy.identity(self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_) + eff_atlevels[ish][sp] = numpy.identity( + self.corr_shells[self.inequiv_to_corr[ish]]['dim'], numpy.complex_) eff_atlevels[ish][sp] *= -self.chemical_potential - eff_atlevels[ish][sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp] + eff_atlevels[ish][ + sp] -= self.dc_imp[self.inequiv_to_corr[ish]][sp] # sum over k: - if not hasattr(self,"Hsumk"): - # calculate the sum over k. Does not depend on mu, so do it only once: - self.Hsumk = [ {} for icrsh in range(self.n_corr_shells) ] + if not hasattr(self, "Hsumk"): + # calculate the sum over k. Does not depend on mu, so do it only + # once: + self.Hsumk = [{} for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh]['dim'] for sp in self.spin_block_names[self.corr_shells[icrsh]['SO']]: - self.Hsumk[icrsh][sp] = numpy.zeros([dim,dim],numpy.complex_) + self.Hsumk[icrsh][sp] = numpy.zeros( + [dim, dim], numpy.complex_) for isp, sp in enumerate(self.spin_block_names[self.corr_shells[icrsh]['SO']]): - ind = self.spin_names_to_ind[self.corr_shells[icrsh]['SO']][sp] + ind = self.spin_names_to_ind[ + self.corr_shells[icrsh]['SO']][sp] for ik in range(self.n_k): - n_orb = self.n_orbitals[ik,ind] + n_orb = self.n_orbitals[ik, ind] MMat = numpy.identity(n_orb, numpy.complex_) - MMat = self.hopping[ik,ind,0:n_orb,0:n_orb] - (1-2*isp) * self.h_field * MMat - projmat = self.proj_mat[ik,ind,icrsh,0:dim,0:n_orb] - self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot( numpy.dot(projmat,MMat), - projmat.conjugate().transpose() ) + MMat = self.hopping[ + ik, ind, 0:n_orb, 0:n_orb] - (1 - 2 * isp) * self.h_field * MMat + projmat = self.proj_mat[ik, ind, icrsh, 0:dim, 0:n_orb] + self.Hsumk[icrsh][sp] += self.bz_weights[ik] * numpy.dot(numpy.dot(projmat, MMat), + projmat.conjugate().transpose()) # symmetrisation: - if self.symm_op != 0: self.Hsumk = self.symmcorr.symmetrize(self.Hsumk) + if self.symm_op != 0: + self.Hsumk = self.symmcorr.symmetrize(self.Hsumk) # Rotate to local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): for sp in self.Hsumk[icrsh]: - if self.rot_mat_time_inv[icrsh] == 1: self.Hsumk[icrsh][sp] = self.Hsumk[icrsh][sp].conjugate() - self.Hsumk[icrsh][sp] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),self.Hsumk[icrsh][sp]) , - self.rot_mat[icrsh] ) + if self.rot_mat_time_inv[icrsh] == 1: + self.Hsumk[icrsh][sp] = self.Hsumk[ + icrsh][sp].conjugate() + self.Hsumk[icrsh][sp] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), self.Hsumk[icrsh][sp]), + self.rot_mat[icrsh]) # add to matrix: for ish in range(self.n_inequiv_shells): for sp in eff_atlevels[ish]: - eff_atlevels[ish][sp] += self.Hsumk[self.inequiv_to_corr[ish]][sp] - + eff_atlevels[ish][ + sp] += self.Hsumk[self.inequiv_to_corr[ish]][sp] return eff_atlevels - def init_dc(self): r""" Initializes the double counting terms. @@ -907,32 +1010,31 @@ class SumkDFT: None """ - self.dc_imp = [ {} for icrsh in range(self.n_corr_shells)] + self.dc_imp = [{} for icrsh in range(self.n_corr_shells)] for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh]['dim'] spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] - for sp in spn: self.dc_imp[icrsh][sp] = numpy.zeros([dim,dim],numpy.float_) + for sp in spn: + self.dc_imp[icrsh][sp] = numpy.zeros([dim, dim], numpy.float_) self.dc_energ = [0.0 for icrsh in range(self.n_corr_shells)] - - def set_dc(self,dc_imp,dc_energ): + def set_dc(self, dc_imp, dc_energ): r""" Sets double counting corrections to given values. - + Parameters ---------- dc_imp : gf_struct_sumk like Double-counting self-energy term. dc_energ : list of floats Double-counting energy corrections for each correlated shell. - + """ self.dc_imp = dc_imp self.dc_energ = dc_energ - - def calc_dc(self,dens_mat,orb=0,U_interact=None,J_hund=None,use_dc_formula=0,use_dc_value=None): + def calc_dc(self, dens_mat, orb=0, U_interact=None, J_hund=None, use_dc_formula=0, use_dc_value=None): r""" Calculates and sets the double counting corrections. @@ -974,68 +1076,87 @@ class SumkDFT: for icrsh in range(self.n_corr_shells): - ish = self.corr_to_inequiv[icrsh] # ish is the index of the inequivalent shell corresponding to icrsh - if ish != orb: continue # ignore this orbital - dim = self.corr_shells[icrsh]['dim'] #*(1+self.corr_shells[icrsh]['SO']) + # ish is the index of the inequivalent shell corresponding to icrsh + ish = self.corr_to_inequiv[icrsh] + if ish != orb: + continue # ignore this orbital + # *(1+self.corr_shells[icrsh]['SO']) + dim = self.corr_shells[icrsh]['dim'] spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] - Ncr = { sp: 0.0 for sp in spn } - for block,inner in self.gf_struct_solver[ish].iteritems(): + Ncr = {sp: 0.0 for sp in spn} + for block, inner in self.gf_struct_solver[ish].iteritems(): bl = self.solver_to_sumk_block[ish][block] Ncr[bl] += dens_mat[block].real.trace() Ncrtot = sum(Ncr.itervalues()) - for sp in spn: - self.dc_imp[icrsh][sp] = numpy.identity(dim,numpy.float_) - if self.SP == 0: # average the densities if there is no SP: + for sp in spn: + self.dc_imp[icrsh][sp] = numpy.identity(dim, numpy.float_) + if self.SP == 0: # average the densities if there is no SP: Ncr[sp] = Ncrtot / len(spn) - elif self.SP == 1 and self.SO == 1: # correction for SO: we have only one block in this case, but in DC we need N/2 + # correction for SO: we have only one block in this case, but + # in DC we need N/2 + elif self.SP == 1 and self.SO == 1: Ncr[sp] = Ncrtot / 2.0 if use_dc_value is None: - if U_interact is None and J_hund is None: raise ValueError, "set_dc: either provide U_interact and J_hund or set use_dc_value to dc value." + if U_interact is None and J_hund is None: + raise ValueError, "set_dc: either provide U_interact and J_hund or set use_dc_value to dc value." - if use_dc_formula == 0: # FLL + if use_dc_formula == 0: # FLL - self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0) + self.dc_energ[icrsh] = U_interact / \ + 2.0 * Ncrtot * (Ncrtot - 1.0) for sp in spn: - Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[sp] - 0.5) + Uav = U_interact * (Ncrtot - 0.5) - \ + J_hund * (Ncr[sp] - 0.5) self.dc_imp[icrsh][sp] *= Uav - self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[sp]) * (Ncr[sp]-1.0) - mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals()) + self.dc_energ[icrsh] -= J_hund / \ + 2.0 * (Ncr[sp]) * (Ncr[sp] - 1.0) + mpi.report( + "DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals()) - elif use_dc_formula == 1: # Held's formula, with U_interact the interorbital onsite interaction + elif use_dc_formula == 1: # Held's formula, with U_interact the interorbital onsite interaction - self.dc_energ[icrsh] = (U_interact + (dim-1)*(U_interact-2.0*J_hund) + (dim-1)*(U_interact-3.0*J_hund))/(2*dim-1) / 2.0 * Ncrtot * (Ncrtot-1.0) + self.dc_energ[icrsh] = (U_interact + (dim - 1) * (U_interact - 2.0 * J_hund) + ( + dim - 1) * (U_interact - 3.0 * J_hund)) / (2 * dim - 1) / 2.0 * Ncrtot * (Ncrtot - 1.0) for sp in spn: - Uav =(U_interact + (dim-1)*(U_interact-2.0*J_hund) + (dim-1)*(U_interact-3.0*J_hund))/(2*dim-1) * (Ncrtot-0.5) + Uav = (U_interact + (dim - 1) * (U_interact - 2.0 * J_hund) + (dim - 1) + * (U_interact - 3.0 * J_hund)) / (2 * dim - 1) * (Ncrtot - 0.5) self.dc_imp[icrsh][sp] *= Uav - mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals()) + mpi.report( + "DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals()) - elif use_dc_formula == 2: # AMF + elif use_dc_formula == 2: # AMF self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot for sp in spn: - Uav = U_interact*(Ncrtot - Ncr[sp]/dim) - J_hund * (Ncr[sp] - Ncr[sp]/dim) + Uav = U_interact * \ + (Ncrtot - Ncr[sp] / dim) - \ + J_hund * (Ncr[sp] - Ncr[sp] / dim) self.dc_imp[icrsh][sp] *= Uav - self.dc_energ[icrsh] -= (U_interact + (dim-1)*J_hund)/dim * 0.5 * Ncr[sp] * Ncr[sp] - mpi.report("DC for shell %(icrsh)i and block %(sp)s = %(Uav)f"%locals()) + self.dc_energ[ + icrsh] -= (U_interact + (dim - 1) * J_hund) / dim * 0.5 * Ncr[sp] * Ncr[sp] + mpi.report( + "DC for shell %(icrsh)i and block %(sp)s = %(Uav)f" % locals()) - mpi.report("DC energy for shell %s = %s"%(icrsh,self.dc_energ[icrsh])) + mpi.report("DC energy for shell %s = %s" % + (icrsh, self.dc_energ[icrsh])) - else: # use value provided for user to determine dc_energ and dc_imp + else: # use value provided for user to determine dc_energ and dc_imp self.dc_energ[icrsh] = use_dc_value * Ncrtot - for sp in spn: self.dc_imp[icrsh][sp] *= use_dc_value + for sp in spn: + self.dc_imp[icrsh][sp] *= use_dc_value - mpi.report("DC for shell %(icrsh)i = %(use_dc_value)f"%locals()) - mpi.report("DC energy = %s"%self.dc_energ[icrsh]) + mpi.report( + "DC for shell %(icrsh)i = %(use_dc_value)f" % locals()) + mpi.report("DC energy = %s" % self.dc_energ[icrsh]) - - def add_dc(self,iw_or_w="iw"): + def add_dc(self, iw_or_w="iw"): r""" Subtracts the double counting term from the impurity self energy. - + Parameters ---------- iw_or_w : string, optional @@ -1051,17 +1172,18 @@ class SumkDFT: """ # Be careful: Sigma_imp is already in the global coordinate system!! - sigma_minus_dc = [s.copy() for s in getattr(self,"Sigma_imp_"+iw_or_w)] + sigma_minus_dc = [s.copy() + for s in getattr(self, "Sigma_imp_" + iw_or_w)] for icrsh in range(self.n_corr_shells): - for bname,gf in sigma_minus_dc[icrsh]: + for bname, gf in sigma_minus_dc[icrsh]: # Transform dc_imp to global coordinate system - dccont = numpy.dot(self.rot_mat[icrsh],numpy.dot(self.dc_imp[icrsh][bname],self.rot_mat[icrsh].conjugate().transpose())) + dccont = numpy.dot(self.rot_mat[icrsh], numpy.dot(self.dc_imp[icrsh][ + bname], self.rot_mat[icrsh].conjugate().transpose())) sigma_minus_dc[icrsh][bname] -= dccont return sigma_minus_dc - - def symm_deg_gf(self,gf_to_symm,orb): + def symm_deg_gf(self, gf_to_symm, orb): r""" Averages a GF over degenerate shells. @@ -1082,9 +1204,10 @@ class SumkDFT: ss = gf_to_symm[degsh[0]].copy() ss.zero() n_deg = len(degsh) - for bl in degsh: ss += gf_to_symm[bl] / (1.0*n_deg) - for bl in degsh: gf_to_symm[bl] << ss - + for bl in degsh: + ss += gf_to_symm[bl] / (1.0 * n_deg) + for bl in degsh: + gf_to_symm[bl] << ss def total_density(self, mu=None, iw_or_w="iw", with_Sigma=True, with_dc=True, broadening=None): r""" @@ -1129,20 +1252,21 @@ class SumkDFT: Total charge :math:`n_{tot}`. """ - if mu is None: mu = self.chemical_potential + if mu is None: + mu = self.chemical_potential dens = 0.0 ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - G_latt = self.lattice_gf(ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening) + G_latt = self.lattice_gf( + ik=ik, mu=mu, iw_or_w=iw_or_w, with_Sigma=with_Sigma, with_dc=with_dc, broadening=broadening) dens += self.bz_weights[ik] * G_latt.total_density() # collect data from mpi: - dens = mpi.all_reduce(mpi.world, dens, lambda x,y : x+y) + dens = mpi.all_reduce(mpi.world, dens, lambda x, y: x + y) mpi.barrier() return dens - - def set_mu(self,mu): + def set_mu(self, mu): r""" Sets a new chemical potential. @@ -1154,8 +1278,7 @@ class SumkDFT: """ self.chemical_potential = mu - - def calc_mu(self, precision=0.01, iw_or_w='iw',broadening=None): + def calc_mu(self, precision=0.01, iw_or_w='iw', broadening=None): r""" Searches for the chemical potential that gives the DFT total charge. A simple bisection method is used. @@ -1179,28 +1302,28 @@ class SumkDFT: within specified precision. """ - F = lambda mu : self.total_density(mu=mu,iw_or_w=iw_or_w,broadening=broadening) + F = lambda mu: self.total_density( + mu=mu, iw_or_w=iw_or_w, broadening=broadening) density = self.density_required - self.charge_below - self.chemical_potential = dichotomy.dichotomy(function = F, - x_init = self.chemical_potential, y_value = density, - precision_on_y = precision, delta_x = 0.5, max_loops = 100, - x_name = "Chemical Potential", y_name = "Total Density", - verbosity = 3)[0] + self.chemical_potential = dichotomy.dichotomy(function=F, + x_init=self.chemical_potential, y_value=density, + precision_on_y=precision, delta_x=0.5, max_loops=100, + x_name="Chemical Potential", y_name="Total Density", + verbosity=3)[0] return self.chemical_potential - def calc_density_correction(self, filename='dens_mat.dat'): r""" Calculates the charge density correction and stores it into a file. - + The charge density correction is needed for charge-self-consistent DFT+DMFT calculations. It represents a density matrix of the interacting system defined in Bloch basis and it is calculated from the sum over Matsubara frequecies of the full GF, ..math:: N_{\nu\nu'}(k) = \sum_{i\omega_{n}} G_{\nu\nu'}(k, i\omega_{n}) - + The density matrix for every `k`-point is stored into a file. Parameters @@ -1216,7 +1339,8 @@ class SumkDFT: """ - assert type(filename) == StringType, "calc_density_correction: filename has to be a string!" + assert type( + filename) == StringType, "calc_density_correction: filename has to be a string!" ntoi = self.spin_names_to_ind[self.SO] spn = self.spin_block_names[self.SO] @@ -1225,61 +1349,74 @@ class SumkDFT: # Set up deltaN: deltaN = {} for sp in spn: - deltaN[sp] = [numpy.zeros([self.n_orbitals[ik,ntoi[sp]],self.n_orbitals[ik,ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)] + deltaN[sp] = [numpy.zeros([self.n_orbitals[ik, ntoi[sp]], self.n_orbitals[ + ik, ntoi[sp]]], numpy.complex_) for ik in range(self.n_k)] ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - G_latt_iw = self.lattice_gf(ik = ik, mu = self.chemical_potential, iw_or_w = "iw") - for bname,gf in G_latt_iw: + G_latt_iw = self.lattice_gf( + ik=ik, mu=self.chemical_potential, iw_or_w="iw") + for bname, gf in G_latt_iw: deltaN[bname][ik] = G_latt_iw[bname].density() - dens[bname] += self.bz_weights[ik] * G_latt_iw[bname].total_density() + dens[bname] += self.bz_weights[ik] * \ + G_latt_iw[bname].total_density() # mpi reduce: for bname in deltaN: for ik in range(self.n_k): - deltaN[bname][ik] = mpi.all_reduce(mpi.world, deltaN[bname][ik], lambda x,y : x+y) - dens[bname] = mpi.all_reduce(mpi.world, dens[bname], lambda x,y : x+y) + deltaN[bname][ik] = mpi.all_reduce( + mpi.world, deltaN[bname][ik], lambda x, y: x + y) + dens[bname] = mpi.all_reduce( + mpi.world, dens[bname], lambda x, y: x + y) mpi.barrier() # now save to file: if mpi.is_master_node(): if self.SP == 0: - f = open(filename,'w') + f = open(filename, 'w') else: - f = open(filename+'up','w') - f1 = open(filename+'dn','w') + f = open(filename + 'up', 'w') + f1 = open(filename + 'dn', 'w') # write chemical potential (in Rydberg): - f.write("%.14f\n"%(self.chemical_potential/self.energy_unit)) - if self.SP != 0: f1.write("%.14f\n"%(self.chemical_potential/self.energy_unit)) + f.write("%.14f\n" % (self.chemical_potential / self.energy_unit)) + if self.SP != 0: + f1.write("%.14f\n" % + (self.chemical_potential / self.energy_unit)) # write beta in rydberg-1 - f.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit)) - if self.SP != 0: f1.write("%.14f\n"%(G_latt_iw.mesh.beta*self.energy_unit)) + f.write("%.14f\n" % (G_latt_iw.mesh.beta * self.energy_unit)) + if self.SP != 0: + f1.write("%.14f\n" % (G_latt_iw.mesh.beta * self.energy_unit)) - if self.SP == 0: # no spin-polarization + if self.SP == 0: # no spin-polarization for ik in range(self.n_k): - f.write("%s\n"%self.n_orbitals[ik,0]) - for inu in range(self.n_orbitals[ik,0]): - for imu in range(self.n_orbitals[ik,0]): - valre = (deltaN['up'][ik][inu,imu].real + deltaN['down'][ik][inu,imu].real) / 2.0 - valim = (deltaN['up'][ik][inu,imu].imag + deltaN['down'][ik][inu,imu].imag) / 2.0 - f.write("%.14f %.14f "%(valre,valim)) + f.write("%s\n" % self.n_orbitals[ik, 0]) + for inu in range(self.n_orbitals[ik, 0]): + for imu in range(self.n_orbitals[ik, 0]): + valre = (deltaN['up'][ik][ + inu, imu].real + deltaN['down'][ik][inu, imu].real) / 2.0 + valim = (deltaN['up'][ik][ + inu, imu].imag + deltaN['down'][ik][inu, imu].imag) / 2.0 + f.write("%.14f %.14f " % (valre, valim)) f.write("\n") f.write("\n") f.close() - elif self.SP == 1: # with spin-polarization + elif self.SP == 1: # with spin-polarization # dict of filename: (spin index, block_name) - if self.SO == 0: to_write = {f: (0, 'up'), f1: (1, 'down')} - if self.SO == 1: to_write = {f: (0, 'ud'), f1: (0, 'ud')} + if self.SO == 0: + to_write = {f: (0, 'up'), f1: (1, 'down')} + if self.SO == 1: + to_write = {f: (0, 'ud'), f1: (0, 'ud')} for fout in to_write.iterkeys(): isp, sp = to_write[fout] for ik in range(self.n_k): - fout.write("%s\n"%self.n_orbitals[ik,isp]) - for inu in range(self.n_orbitals[ik,isp]): - for imu in range(self.n_orbitals[ik,isp]): - fout.write("%.14f %.14f "%(deltaN[sp][ik][inu,imu].real,deltaN[sp][ik][inu,imu].imag)) + fout.write("%s\n" % self.n_orbitals[ik, isp]) + for inu in range(self.n_orbitals[ik, isp]): + for imu in range(self.n_orbitals[ik, isp]): + fout.write("%.14f %.14f " % (deltaN[sp][ik][ + inu, imu].real, deltaN[sp][ik][inu, imu].imag)) fout.write("\n") fout.write("\n") fout.close() @@ -1290,67 +1427,69 @@ class SumkDFT: # FIXME LEAVE UNDOCUMENTED ################ - def calc_dc_for_density(self,orb,dc_init,dens_mat,density=None,precision=0.01): + def calc_dc_for_density(self, orb, dc_init, dens_mat, density=None, precision=0.01): """Searches for DC in order to fulfill charge neutrality. If density is given, then DC is set such that the LOCAL charge of orbital orb coincides with the given density.""" def F(dc): - self.calc_dc(dens_mat = dens_mat, U_interact = 0, J_hund = 0, orb = orb, use_dc_value = dc) + self.calc_dc(dens_mat=dens_mat, U_interact=0, + J_hund=0, orb=orb, use_dc_value=dc) if dens_req is None: - return self.total_density(mu = mu) + return self.total_density(mu=mu) else: return self.extract_G_loc()[orb].total_density() - if density is None: density = self.density_required - self.charge_below + if density is None: + density = self.density_required - self.charge_below - dc = dichotomy.dichotomy(function = F, - x_init = dc_init, y_value = density, - precision_on_y = precision, delta_x = 0.5, - max_loops = 100, x_name = "Double Counting", y_name= "Total Density", - verbosity = 3)[0] + dc = dichotomy.dichotomy(function=F, + x_init=dc_init, y_value=density, + precision_on_y=precision, delta_x=0.5, + max_loops=100, x_name="Double Counting", y_name="Total Density", + verbosity=3)[0] return dc - def check_projectors(self): """Calculated the density matrix from projectors (DM = P Pdagger) to check that it is correct and specifically that it matches DFT.""" - dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'],self.corr_shells[icrsh]['dim']],numpy.complex_) - for icrsh in range(self.n_corr_shells)] + dens_mat = [numpy.zeros([self.corr_shells[icrsh]['dim'], self.corr_shells[icrsh]['dim']], numpy.complex_) + for icrsh in range(self.n_corr_shells)] for ik in range(self.n_k): for icrsh in range(self.n_corr_shells): dim = self.corr_shells[icrsh]['dim'] - n_orb = self.n_orbitals[ik,0] - projmat = self.proj_mat[ik,0,icrsh,0:dim,0:n_orb] - dens_mat[icrsh][:,:] += numpy.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik] + n_orb = self.n_orbitals[ik, 0] + projmat = self.proj_mat[ik, 0, icrsh, 0:dim, 0:n_orb] + dens_mat[icrsh][ + :, :] += numpy.dot(projmat, projmat.transpose().conjugate()) * self.bz_weights[ik] - if self.symm_op != 0: dens_mat = self.symmcorr.symmetrize(dens_mat) + if self.symm_op != 0: + dens_mat = self.symmcorr.symmetrize(dens_mat) # Rotate to local coordinate system: if self.use_rotations: for icrsh in range(self.n_corr_shells): - if self.rot_mat_time_inv[icrsh] == 1: dens_mat[icrsh] = dens_mat[icrsh].conjugate() - dens_mat[icrsh] = numpy.dot( numpy.dot(self.rot_mat[icrsh].conjugate().transpose(),dens_mat[icrsh]), - self.rot_mat[icrsh] ) + if self.rot_mat_time_inv[icrsh] == 1: + dens_mat[icrsh] = dens_mat[icrsh].conjugate() + dens_mat[icrsh] = numpy.dot(numpy.dot(self.rot_mat[icrsh].conjugate().transpose(), dens_mat[icrsh]), + self.rot_mat[icrsh]) return dens_mat - - def sorts_of_atoms(self,shells): + def sorts_of_atoms(self, shells): """ Determine the number of inequivalent sorts. """ - sortlst = [ shells[i]['sort'] for i in range(len(shells)) ] + sortlst = [shells[i]['sort'] for i in range(len(shells))] n_sorts = len(set(sortlst)) return n_sorts - - def number_of_atoms(self,shells): + def number_of_atoms(self, shells): """ Determine the number of inequivalent atoms. """ - atomlst = [ shells[i]['atom'] for i in range(len(shells)) ] + atomlst = [shells[i]['atom'] for i in range(len(shells))] n_atoms = len(set(atomlst)) return n_atoms diff --git a/python/sumk_dft_tools.py b/python/sumk_dft_tools.py index 9d369b2a..046735ab 100644 --- a/python/sumk_dft_tools.py +++ b/python/sumk_dft_tools.py @@ -1,4 +1,4 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -17,7 +17,7 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## import sys from types import * import numpy @@ -28,30 +28,29 @@ from sumk_dft import SumkDFT from scipy.integrate import * from scipy.interpolate import * + class SumkDFTTools(SumkDFT): """ Extends the SumkDFT class with some tools for analysing the data. """ - - def __init__(self, hdf_file, h_field = 0.0, use_dft_blocks = False, dft_data = 'dft_input', symmcorr_data = 'dft_symmcorr_input', - parproj_data = 'dft_parproj_input', symmpar_data = 'dft_symmpar_input', bands_data = 'dft_bands_input', - transp_data = 'dft_transp_input', misc_data = 'dft_misc_input'): + def __init__(self, hdf_file, h_field=0.0, use_dft_blocks=False, dft_data='dft_input', symmcorr_data='dft_symmcorr_input', + parproj_data='dft_parproj_input', symmpar_data='dft_symmpar_input', bands_data='dft_bands_input', + transp_data='dft_transp_input', misc_data='dft_misc_input'): """ Initialisation of the class. Parameters are exactly as for SumKDFT. """ - - SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks, - dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data, - symmpar_data=symmpar_data, bands_data=bands_data, transp_data=transp_data, - misc_data=misc_data) + SumkDFT.__init__(self, hdf_file=hdf_file, h_field=h_field, use_dft_blocks=use_dft_blocks, + dft_data=dft_data, symmcorr_data=symmcorr_data, parproj_data=parproj_data, + symmpar_data=symmpar_data, bands_data=bands_data, transp_data=transp_data, + misc_data=misc_data) # Uses .data of only GfReFreq objects. def dos_wannier_basis(self, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True, save_to_file=True): """ Calculates the density of states in the basis of the Wannier functions. - + Parameters ---------- mu : double, optional @@ -66,7 +65,7 @@ class SumkDFTTools(SumkDFT): If True the double counting correction is used. save_to_file : boolean, optional If True, text files with the calculated data will be created. - + Returns ------- DOS : Dict of numpy arrays @@ -76,7 +75,7 @@ class SumkDFTTools(SumkDFT): DOSproj_orb : Dict of numpy arrays DOS projected to atoms and resolved into orbital contributions. """ - + if (mesh is None) and (not with_Sigma): raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq." if mesh is None: @@ -84,93 +83,115 @@ class SumkDFTTools(SumkDFT): om_min = om_mesh[0] om_max = om_mesh[-1] n_om = len(om_mesh) - mesh = (om_min,om_max,n_om) - else: - om_min,om_max,n_om = mesh + mesh = (om_min, om_max, n_om) + else: + om_min, om_max, n_om = mesh om_mesh = numpy.linspace(om_min, om_max, n_om) - + G_loc = [] for icrsh in range(self.n_corr_shells): spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] - glist = [ GfReFreq(indices = inner, window = (om_min,om_max), n_points = n_om) for block,inner in self.gf_struct_sumk[icrsh]] - G_loc.append(BlockGf(name_list = spn, block_list = glist, make_copies=False)) - for icrsh in range(self.n_corr_shells): G_loc[icrsh].zero() + glist = [GfReFreq(indices=inner, window=(om_min, om_max), n_points=n_om) + for block, inner in self.gf_struct_sumk[icrsh]] + G_loc.append( + BlockGf(name_list=spn, block_list=glist, make_copies=False)) + for icrsh in range(self.n_corr_shells): + G_loc[icrsh].zero() - DOS = { sp: numpy.zeros([n_om],numpy.float_) for sp in self.spin_block_names[self.SO] } - DOSproj = [ {} for ish in range(self.n_inequiv_shells) ] - DOSproj_orb = [ {} for ish in range(self.n_inequiv_shells) ] + DOS = {sp: numpy.zeros([n_om], numpy.float_) + for sp in self.spin_block_names[self.SO]} + DOSproj = [{} for ish in range(self.n_inequiv_shells)] + DOSproj_orb = [{} for ish in range(self.n_inequiv_shells)] for ish in range(self.n_inequiv_shells): for sp in self.spin_block_names[self.corr_shells[self.inequiv_to_corr[ish]]['SO']]: dim = self.corr_shells[self.inequiv_to_corr[ish]]['dim'] - DOSproj[ish][sp] = numpy.zeros([n_om],numpy.float_) - DOSproj_orb[ish][sp] = numpy.zeros([n_om,dim,dim],numpy.float_) + DOSproj[ish][sp] = numpy.zeros([n_om], numpy.float_) + DOSproj_orb[ish][sp] = numpy.zeros( + [n_om, dim, dim], numpy.float_) ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening,mesh=mesh,with_Sigma=with_Sigma,with_dc=with_dc) + G_latt_w = self.lattice_gf( + ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc) G_latt_w *= self.bz_weights[ik] # Non-projected DOS for iom in range(n_om): - for bname,gf in G_latt_w: - DOS[bname][iom] -= gf.data[iom,:,:].imag.trace()/numpy.pi + for bname, gf in G_latt_w: + DOS[bname][iom] -= gf.data[iom, :, :].imag.trace() / \ + numpy.pi # Projected DOS: for icrsh in range(self.n_corr_shells): tmp = G_loc[icrsh].copy() - for bname,gf in tmp: tmp[bname] << self.downfold(ik,icrsh,bname,G_latt_w[bname],gf) # downfolding G + for bname, gf in tmp: + tmp[bname] << self.downfold(ik, icrsh, bname, G_latt_w[ + bname], gf) # downfolding G G_loc[icrsh] += tmp # Collect data from mpi: for bname in DOS: - DOS[bname] = mpi.all_reduce(mpi.world, DOS[bname], lambda x,y : x+y) + DOS[bname] = mpi.all_reduce( + mpi.world, DOS[bname], lambda x, y: x + y) for icrsh in range(self.n_corr_shells): - G_loc[icrsh] << mpi.all_reduce(mpi.world, G_loc[icrsh], lambda x,y : x+y) + G_loc[icrsh] << mpi.all_reduce( + mpi.world, G_loc[icrsh], lambda x, y: x + y) mpi.barrier() # Symmetrize and rotate to local coord. system if needed: - if self.symm_op != 0: G_loc = self.symmcorr.symmetrize(G_loc) + if self.symm_op != 0: + G_loc = self.symmcorr.symmetrize(G_loc) if self.use_rotations: for icrsh in range(self.n_corr_shells): - for bname,gf in G_loc[icrsh]: G_loc[icrsh][bname] << self.rotloc(icrsh,gf,direction='toLocal') + for bname, gf in G_loc[icrsh]: + G_loc[icrsh][bname] << self.rotloc( + icrsh, gf, direction='toLocal') # G_loc can now also be used to look at orbitally-resolved quantities for ish in range(self.n_inequiv_shells): - for bname,gf in G_loc[self.inequiv_to_corr[ish]]: # loop over spins - for iom in range(n_om): DOSproj[ish][bname][iom] -= gf.data[iom,:,:].imag.trace()/numpy.pi - DOSproj_orb[ish][bname][:,:,:] -= gf.data[:,:,:].imag/numpy.pi + for bname, gf in G_loc[self.inequiv_to_corr[ish]]: # loop over spins + for iom in range(n_om): + DOSproj[ish][bname][iom] -= gf.data[iom, + :, :].imag.trace() / numpy.pi + DOSproj_orb[ish][bname][ + :, :, :] -= gf.data[:, :, :].imag / numpy.pi # Write to files if save_to_file and mpi.is_master_node(): for sp in self.spin_block_names[self.SO]: - f = open('DOS_wann_%s.dat'%sp, 'w') - for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOS[sp][iom])) + f = open('DOS_wann_%s.dat' % sp, 'w') + for iom in range(n_om): + f.write("%s %s\n" % (om_mesh[iom], DOS[sp][iom])) f.close() # Partial for ish in range(self.n_inequiv_shells): - f = open('DOS_wann_%s_proj%s.dat'%(sp,ish),'w') - for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj[ish][sp][iom])) + f = open('DOS_wann_%s_proj%s.dat' % (sp, ish), 'w') + for iom in range(n_om): + f.write("%s %s\n" % + (om_mesh[iom], DOSproj[ish][sp][iom])) f.close() # Orbitally-resolved for i in range(self.corr_shells[self.inequiv_to_corr[ish]]['dim']): - for j in range(i,self.corr_shells[self.inequiv_to_corr[ish]]['dim']): - f = open('DOS_wann_'+sp+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat','w') - for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][sp][iom,i,j])) + for j in range(i, self.corr_shells[self.inequiv_to_corr[ish]]['dim']): + f = open('DOS_wann_' + sp + '_proj' + str(ish) + + '_' + str(i) + '_' + str(j) + '.dat', 'w') + for iom in range(n_om): + f.write("%s %s\n" % ( + om_mesh[iom], DOSproj_orb[ish][sp][iom, i, j])) f.close() return DOS, DOSproj, DOSproj_orb - # Uses .data of only GfReFreq objects. def dos_parproj_basis(self, mu=None, broadening=None, mesh=None, with_Sigma=True, with_dc=True, save_to_file=True): """ Calculates the orbitally-resolved DOS. Different to dos_Wannier_basis is that here we calculate projections also to non-Wannier projectors, in the flavour of Wien2k QTL calculatuions. - + Parameters ---------- mu : double, optional @@ -185,7 +206,7 @@ class SumkDFTTools(SumkDFT): If True the double counting correction is used. save_to_file : boolean, optional If True, text files with the calculated data will be created. - + Returns ------- DOS : Dict of numpy arrays @@ -196,107 +217,132 @@ class SumkDFTTools(SumkDFT): DOS projected to atoms and resolved into orbital contributions. """ - - things_to_read = ['n_parproj','proj_mat_all','rot_mat_all','rot_mat_all_time_inv'] - value_read = self.read_input_from_hdf(subgrp=self.parproj_data,things_to_read = things_to_read) - if not value_read: return value_read - if self.symm_op: self.symmpar = Symmetry(self.hdf_file,subgroup=self.symmpar_data) + things_to_read = ['n_parproj', 'proj_mat_all', + 'rot_mat_all', 'rot_mat_all_time_inv'] + value_read = self.read_input_from_hdf( + subgrp=self.parproj_data, things_to_read=things_to_read) + if not value_read: + return value_read + if self.symm_op: + self.symmpar = Symmetry(self.hdf_file, subgroup=self.symmpar_data) if (mesh is None) and (not with_Sigma): raise ValueError, "lattice_gf: Give the mesh=(om_min,om_max,n_points) for the lattice GfReFreq." - if mesh is None: + if mesh is None: om_mesh = [x.real for x in self.Sigma_imp_w[0].mesh] om_min = om_mesh[0] om_max = om_mesh[-1] n_om = len(om_mesh) - mesh = (om_min,om_max,n_om) - else: - om_min,om_max,n_om = mesh + mesh = (om_min, om_max, n_om) + else: + om_min, om_max, n_om = mesh om_mesh = numpy.linspace(om_min, om_max, n_om) - + G_loc = [] spn = self.spin_block_names[self.SO] - gf_struct_parproj = [ [ (sp, range(self.shells[ish]['dim'])) for sp in spn ] - for ish in range(self.n_shells) ] - for ish in range(self.n_shells): - glist = [ GfReFreq(indices = inner, window = (om_min,om_max), n_points = n_om) for block,inner in gf_struct_parproj[ish] ] - G_loc.append(BlockGf(name_list = spn, block_list = glist, make_copies=False)) - for ish in range(self.n_shells): G_loc[ish].zero() + gf_struct_parproj = [[(sp, range(self.shells[ish]['dim'])) for sp in spn] + for ish in range(self.n_shells)] + for ish in range(self.n_shells): + glist = [GfReFreq(indices=inner, window=(om_min, om_max), n_points=n_om) + for block, inner in gf_struct_parproj[ish]] + G_loc.append( + BlockGf(name_list=spn, block_list=glist, make_copies=False)) + for ish in range(self.n_shells): + G_loc[ish].zero() - DOS = { sp: numpy.zeros([n_om],numpy.float_) for sp in self.spin_block_names[self.SO] } - DOSproj = [ {} for ish in range(self.n_shells) ] - DOSproj_orb = [ {} for ish in range(self.n_shells) ] + DOS = {sp: numpy.zeros([n_om], numpy.float_) + for sp in self.spin_block_names[self.SO]} + DOSproj = [{} for ish in range(self.n_shells)] + DOSproj_orb = [{} for ish in range(self.n_shells)] for ish in range(self.n_shells): for sp in self.spin_block_names[self.SO]: dim = self.shells[ish]['dim'] - DOSproj[ish][sp] = numpy.zeros([n_om],numpy.float_) - DOSproj_orb[ish][sp] = numpy.zeros([n_om,dim,dim],numpy.float_) + DOSproj[ish][sp] = numpy.zeros([n_om], numpy.float_) + DOSproj_orb[ish][sp] = numpy.zeros( + [n_om, dim, dim], numpy.float_) ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening,mesh=mesh,with_Sigma=with_Sigma,with_dc=with_dc) + G_latt_w = self.lattice_gf( + ik=ik, mu=mu, iw_or_w="w", broadening=broadening, mesh=mesh, with_Sigma=with_Sigma, with_dc=with_dc) G_latt_w *= self.bz_weights[ik] # Non-projected DOS for iom in range(n_om): - for bname,gf in G_latt_w: - DOS[bname][iom] -= gf.data[iom,:,:].imag.trace()/numpy.pi + for bname, gf in G_latt_w: + DOS[bname][iom] -= gf.data[iom, :, :].imag.trace() / \ + numpy.pi # Projected DOS: for ish in range(self.n_shells): tmp = G_loc[ish].copy() for ir in range(self.n_parproj[ish]): - for bname,gf in tmp: tmp[bname] << self.downfold(ik,ish,bname,G_latt_w[bname],gf,shells='all',ir=ir) + for bname, gf in tmp: + tmp[bname] << self.downfold(ik, ish, bname, G_latt_w[ + bname], gf, shells='all', ir=ir) G_loc[ish] += tmp # Collect data from mpi: for bname in DOS: - DOS[bname] = mpi.all_reduce(mpi.world, DOS[bname], lambda x,y : x+y) + DOS[bname] = mpi.all_reduce( + mpi.world, DOS[bname], lambda x, y: x + y) for ish in range(self.n_shells): - G_loc[ish] << mpi.all_reduce(mpi.world, G_loc[ish], lambda x,y : x+y) + G_loc[ish] << mpi.all_reduce( + mpi.world, G_loc[ish], lambda x, y: x + y) mpi.barrier() # Symmetrize and rotate to local coord. system if needed: - if self.symm_op != 0: G_loc = self.symmpar.symmetrize(G_loc) + if self.symm_op != 0: + G_loc = self.symmpar.symmetrize(G_loc) if self.use_rotations: for ish in range(self.n_shells): - for bname,gf in G_loc[ish]: G_loc[ish][bname] << self.rotloc(ish,gf,direction='toLocal',shells='all') + for bname, gf in G_loc[ish]: + G_loc[ish][bname] << self.rotloc( + ish, gf, direction='toLocal', shells='all') # G_loc can now also be used to look at orbitally-resolved quantities for ish in range(self.n_shells): - for bname,gf in G_loc[ish]: - for iom in range(n_om): DOSproj[ish][bname][iom] -= gf.data[iom,:,:].imag.trace()/numpy.pi - DOSproj_orb[ish][bname][:,:,:] -= gf.data[:,:,:].imag/numpy.pi + for bname, gf in G_loc[ish]: + for iom in range(n_om): + DOSproj[ish][bname][iom] -= gf.data[iom, + :, :].imag.trace() / numpy.pi + DOSproj_orb[ish][bname][ + :, :, :] -= gf.data[:, :, :].imag / numpy.pi # Write to files if save_to_file and mpi.is_master_node(): for sp in self.spin_block_names[self.SO]: - f = open('DOS_parproj_%s.dat'%sp, 'w') - for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOS[sp][iom])) + f = open('DOS_parproj_%s.dat' % sp, 'w') + for iom in range(n_om): + f.write("%s %s\n" % (om_mesh[iom], DOS[sp][iom])) f.close() # Partial for ish in range(self.n_shells): - f = open('DOS_parproj_%s_proj%s.dat'%(sp,ish),'w') - for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj[ish][sp][iom])) + f = open('DOS_parproj_%s_proj%s.dat' % (sp, ish), 'w') + for iom in range(n_om): + f.write("%s %s\n" % + (om_mesh[iom], DOSproj[ish][sp][iom])) f.close() # Orbitally-resolved for i in range(self.shells[ish]['dim']): - for j in range(i,self.shells[ish]['dim']): - f = open('DOS_parproj_'+sp+'_proj'+str(ish)+'_'+str(i)+'_'+str(j)+'.dat','w') - for iom in range(n_om): f.write("%s %s\n"%(om_mesh[iom],DOSproj_orb[ish][sp][iom,i,j])) + for j in range(i, self.shells[ish]['dim']): + f = open('DOS_parproj_' + sp + '_proj' + str(ish) + + '_' + str(i) + '_' + str(j) + '.dat', 'w') + for iom in range(n_om): + f.write("%s %s\n" % ( + om_mesh[iom], DOSproj_orb[ish][sp][iom, i, j])) f.close() return DOS, DOSproj, DOSproj_orb - # Uses .data of only GfReFreq objects. - def spaghettis(self,broadening=None,plot_shift=0.0,plot_range=None,ishell=None,mu=None,save_to_file='Akw_'): + def spaghettis(self, broadening=None, plot_shift=0.0, plot_range=None, ishell=None, mu=None, save_to_file='Akw_'): """ Calculates the correlated band structure using a real-frequency self energy. - + Parameters ---------- mu : double, optional @@ -311,119 +357,145 @@ class SumkDFTTools(SumkDFT): Contains the index of the shell on which the spectral function is projected. If ishell=None, the total spectrum without projection is calculated. save_to_file : string, optional Filename where the spectra are stored. - + Returns ------- Akw : Dict of numpy arrays Data as it is also written to the files. """ + assert hasattr( + self, "Sigma_imp_w"), "spaghettis: Set Sigma_imp_w first." + things_to_read = ['n_k', 'n_orbitals', 'proj_mat', + 'hopping', 'n_parproj', 'proj_mat_all'] + value_read = self.read_input_from_hdf( + subgrp=self.bands_data, things_to_read=things_to_read) + if not value_read: + return value_read + things_to_read = ['rot_mat_all', 'rot_mat_all_time_inv'] + value_read = self.read_input_from_hdf( + subgrp=self.parproj_data, things_to_read=things_to_read) + if not value_read: + return value_read - assert hasattr(self,"Sigma_imp_w"), "spaghettis: Set Sigma_imp_w first." - things_to_read = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_all'] - value_read = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=things_to_read) - if not value_read: return value_read - things_to_read = ['rot_mat_all','rot_mat_all_time_inv'] - value_read = self.read_input_from_hdf(subgrp=self.parproj_data,things_to_read = things_to_read) - if not value_read: return value_read - - if mu is None: mu = self.chemical_potential + if mu is None: + mu = self.chemical_potential spn = self.spin_block_names[self.SO] mesh = [x.real for x in self.Sigma_imp_w[0].mesh] n_om = len(mesh) if plot_range is None: om_minplot = mesh[0] - 0.001 - om_maxplot = mesh[n_om-1] + 0.001 + om_maxplot = mesh[n_om - 1] + 0.001 else: om_minplot = plot_range[0] om_maxplot = plot_range[1] if ishell is None: - Akw = { sp: numpy.zeros([self.n_k,n_om],numpy.float_) for sp in spn } + Akw = {sp: numpy.zeros([self.n_k, n_om], numpy.float_) + for sp in spn} else: - Akw = { sp: numpy.zeros([self.shells[ishell]['dim'],self.n_k,n_om],numpy.float_) for sp in spn } + Akw = {sp: numpy.zeros( + [self.shells[ishell]['dim'], self.n_k, n_om], numpy.float_) for sp in spn} if not ishell is None: - gf_struct_parproj = [ (sp, range(self.shells[ishell]['dim'])) for sp in spn ] - G_loc = BlockGf(name_block_generator = [ (block,GfReFreq(indices = inner, mesh = self.Sigma_imp_w[0].mesh)) - for block,inner in gf_struct_parproj ], make_copies = False) + gf_struct_parproj = [ + (sp, range(self.shells[ishell]['dim'])) for sp in spn] + G_loc = BlockGf(name_block_generator=[(block, GfReFreq(indices=inner, mesh=self.Sigma_imp_w[0].mesh)) + for block, inner in gf_struct_parproj], make_copies=False) G_loc.zero() ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - G_latt_w = self.lattice_gf(ik=ik,mu=mu,iw_or_w="w",broadening=broadening) + G_latt_w = self.lattice_gf( + ik=ik, mu=mu, iw_or_w="w", broadening=broadening) if ishell is None: # Non-projected A(k,w) for iom in range(n_om): if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): - for bname,gf in G_latt_w: Akw[bname][ik,iom] += gf.data[iom,:,:].imag.trace()/(-1.0*numpy.pi) - Akw[bname][ik,iom] += ik*plot_shift # shift Akw for plotting stacked k-resolved eps(k) curves + for bname, gf in G_latt_w: + Akw[bname][ik, iom] += gf.data[iom, :, + :].imag.trace() / (-1.0 * numpy.pi) + # shift Akw for plotting stacked k-resolved eps(k) + # curves + Akw[bname][ik, iom] += ik * plot_shift - else: # ishell not None + else: # ishell not None # Projected A(k,w): G_loc.zero() tmp = G_loc.copy() for ir in range(self.n_parproj[ishell]): - for bname,gf in tmp: tmp[bname] << self.downfold(ik,ishell,bname,G_latt_w[bname],gf,shells='all',ir=ir) + for bname, gf in tmp: + tmp[bname] << self.downfold(ik, ishell, bname, G_latt_w[ + bname], gf, shells='all', ir=ir) G_loc += tmp # Rotate to local frame if self.use_rotations: - for bname,gf in G_loc: G_loc[bname] << self.rotloc(ishell,gf,direction='toLocal',shells='all') + for bname, gf in G_loc: + G_loc[bname] << self.rotloc( + ishell, gf, direction='toLocal', shells='all') for iom in range(n_om): if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): for ish in range(self.shells[ishell]['dim']): for sp in spn: - Akw[sp][ish,ik,iom] = G_loc[sp].data[iom,ish,ish].imag/(-1.0*numpy.pi) + Akw[sp][ish, ik, iom] = G_loc[sp].data[ + iom, ish, ish].imag / (-1.0 * numpy.pi) # Collect data from mpi for sp in spn: - Akw[sp] = mpi.all_reduce(mpi.world, Akw[sp], lambda x,y : x+y) + Akw[sp] = mpi.all_reduce(mpi.world, Akw[sp], lambda x, y: x + y) mpi.barrier() if save_to_file and mpi.is_master_node(): if ishell is None: - for sp in spn: # loop over GF blocs: - f = open(save_to_file+sp+'.dat','w') # Open file for storage: + for sp in spn: # loop over GF blocs: + # Open file for storage: + f = open(save_to_file + sp + '.dat', 'w') for ik in range(self.n_k): for iom in range(n_om): if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): if plot_shift > 0.0001: - f.write('%s %s\n'%(mesh[iom],Akw[sp][ik,iom])) + f.write('%s %s\n' % + (mesh[iom], Akw[sp][ik, iom])) else: - f.write('%s %s %s\n'%(ik,mesh[iom],Akw[sp][ik,iom])) + f.write('%s %s %s\n' % + (ik, mesh[iom], Akw[sp][ik, iom])) f.write('\n') f.close() - else: # ishell is not None + else: # ishell is not None for sp in spn: for ish in range(self.shells[ishell]['dim']): - f = open(save_to_file+str(ishell)+'_'+sp+'_proj'+str(ish)+'.dat','w') # Open file for storage: + # Open file for storage: + f = open(save_to_file + str(ishell) + '_' + + sp + '_proj' + str(ish) + '.dat', 'w') for ik in range(self.n_k): for iom in range(n_om): if (mesh[iom] > om_minplot) and (mesh[iom] < om_maxplot): if plot_shift > 0.0001: - f.write('%s %s\n'%(mesh[iom],Akw[sp][ish,ik,iom])) + f.write('%s %s\n' % ( + mesh[iom], Akw[sp][ish, ik, iom])) else: - f.write('%s %s %s\n'%(ik,mesh[iom],Akw[sp][ish,ik,iom])) + f.write('%s %s %s\n' % ( + ik, mesh[iom], Akw[sp][ish, ik, iom])) f.write('\n') f.close() return Akw - def partial_charges(self,beta=40,mu=None,with_Sigma=True,with_dc=True): + def partial_charges(self, beta=40, mu=None, with_Sigma=True, with_dc=True): """ Calculates the orbitally-resolved density matrix for all the orbitals considered in the input, consistent with the definition of Wien2k. Hence, (possibly non-orthonormal) projectors have to be provided in the partial projectors subgroup of the hdf5 archive. - + Parameters ---------- - + with_Sigma : boolean, optional If True, the self energy is used for the calculation. If false, partial charges are calculated without self-energy correction. beta : double, optional @@ -432,96 +504,110 @@ class SumkDFTTools(SumkDFT): Chemical potential, overrides the one stored in the hdf5 archive. with_dc : boolean, optional If True the double counting correction is used. - + Returns ------- dens_mat : list of numpy array A list of density matrices projected to all shells provided in the input. """ - - things_to_read = ['dens_mat_below','n_parproj','proj_mat_all','rot_mat_all','rot_mat_all_time_inv'] - value_read = self.read_input_from_hdf(subgrp=self.parproj_data,things_to_read = things_to_read) - if not value_read: return value_read - if self.symm_op: self.symmpar = Symmetry(self.hdf_file,subgroup=self.symmpar_data) + + things_to_read = ['dens_mat_below', 'n_parproj', + 'proj_mat_all', 'rot_mat_all', 'rot_mat_all_time_inv'] + value_read = self.read_input_from_hdf( + subgrp=self.parproj_data, things_to_read=things_to_read) + if not value_read: + return value_read + if self.symm_op: + self.symmpar = Symmetry(self.hdf_file, subgroup=self.symmpar_data) spn = self.spin_block_names[self.SO] ntoi = self.spin_names_to_ind[self.SO] # Density matrix in the window - self.dens_mat_window = [ [ numpy.zeros([self.shells[ish]['dim'],self.shells[ish]['dim']],numpy.complex_) - for ish in range(self.n_shells) ] - for isp in range(len(spn)) ] + self.dens_mat_window = [[numpy.zeros([self.shells[ish]['dim'], self.shells[ish]['dim']], numpy.complex_) + for ish in range(self.n_shells)] + for isp in range(len(spn))] # Set up G_loc - gf_struct_parproj = [ [ (sp, range(self.shells[ish]['dim'])) for sp in spn ] - for ish in range(self.n_shells) ] + gf_struct_parproj = [[(sp, range(self.shells[ish]['dim'])) for sp in spn] + for ish in range(self.n_shells)] if with_Sigma: - G_loc = [ BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, mesh = self.Sigma_imp_iw[0].mesh)) - for block,inner in gf_struct_parproj[ish] ], make_copies = False) + G_loc = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, mesh=self.Sigma_imp_iw[0].mesh)) + for block, inner in gf_struct_parproj[ish]], make_copies=False) for ish in range(self.n_shells)] beta = self.Sigma_imp_iw[0].mesh.beta else: - G_loc = [ BlockGf(name_block_generator = [ (block,GfImFreq(indices = inner, beta = beta)) - for block,inner in gf_struct_parproj[ish] ], make_copies = False) + G_loc = [BlockGf(name_block_generator=[(block, GfImFreq(indices=inner, beta=beta)) + for block, inner in gf_struct_parproj[ish]], make_copies=False) for ish in range(self.n_shells)] - for ish in range(self.n_shells): G_loc[ish].zero() + for ish in range(self.n_shells): + G_loc[ish].zero() ikarray = numpy.array(range(self.n_k)) for ik in mpi.slice_array(ikarray): - G_latt_iw = self.lattice_gf(ik=ik,mu=mu,iw_or_w="iw",beta=beta,with_Sigma=with_Sigma,with_dc=with_dc) + G_latt_iw = self.lattice_gf( + ik=ik, mu=mu, iw_or_w="iw", beta=beta, with_Sigma=with_Sigma, with_dc=with_dc) G_latt_iw *= self.bz_weights[ik] for ish in range(self.n_shells): tmp = G_loc[ish].copy() for ir in range(self.n_parproj[ish]): - for bname,gf in tmp: tmp[bname] << self.downfold(ik,ish,bname,G_latt_iw[bname],gf,shells='all',ir=ir) + for bname, gf in tmp: + tmp[bname] << self.downfold(ik, ish, bname, G_latt_iw[ + bname], gf, shells='all', ir=ir) G_loc[ish] += tmp # Collect data from mpi: for ish in range(self.n_shells): - G_loc[ish] << mpi.all_reduce(mpi.world, G_loc[ish], lambda x,y : x+y) + G_loc[ish] << mpi.all_reduce( + mpi.world, G_loc[ish], lambda x, y: x + y) mpi.barrier() # Symmetrize and rotate to local coord. system if needed: - if self.symm_op != 0: G_loc = self.symmpar.symmetrize(G_loc) + if self.symm_op != 0: + G_loc = self.symmpar.symmetrize(G_loc) if self.use_rotations: for ish in range(self.n_shells): - for bname,gf in G_loc[ish]: G_loc[ish][bname] << self.rotloc(ish,gf,direction='toLocal',shells='all') + for bname, gf in G_loc[ish]: + G_loc[ish][bname] << self.rotloc( + ish, gf, direction='toLocal', shells='all') for ish in range(self.n_shells): isp = 0 - for bname,gf in G_loc[ish]: + for bname, gf in G_loc[ish]: self.dens_mat_window[isp][ish] = G_loc[ish].density()[bname] isp += 1 # Add density matrices to get the total: - dens_mat = [ [ self.dens_mat_below[ntoi[spn[isp]]][ish] + self.dens_mat_window[isp][ish] - for ish in range(self.n_shells) ] - for isp in range(len(spn)) ] + dens_mat = [[self.dens_mat_below[ntoi[spn[isp]]][ish] + self.dens_mat_window[isp][ish] + for ish in range(self.n_shells)] + for isp in range(len(spn))] return dens_mat - def print_hamiltonian(self): """ Prints the Kohn-Sham Hamiltonian to the text files hamup.dat and hamdn.dat (no spin orbit-coupling), or to ham.dat (with spin-orbit coupling). """ if self.SP == 1 and self.SO == 0: - f1 = open('hamup.dat','w') - f2 = open('hamdn.dat','w') + f1 = open('hamup.dat', 'w') + f2 = open('hamdn.dat', 'w') for ik in range(self.n_k): - for i in range(self.n_orbitals[ik,0]): - f1.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real)) - for i in range(self.n_orbitals[ik,1]): - f2.write('%s %s\n'%(ik,self.hopping[ik,1,i,i].real)) + for i in range(self.n_orbitals[ik, 0]): + f1.write('%s %s\n' % + (ik, self.hopping[ik, 0, i, i].real)) + for i in range(self.n_orbitals[ik, 1]): + f2.write('%s %s\n' % + (ik, self.hopping[ik, 1, i, i].real)) f1.write('\n') f2.write('\n') f1.close() f2.close() else: - f = open('ham.dat','w') + f = open('ham.dat', 'w') for ik in range(self.n_k): - for i in range(self.n_orbitals[ik,0]): - f.write('%s %s\n'%(ik,self.hopping[ik,0,i,i].real)) + for i in range(self.n_orbitals[ik, 0]): + f.write('%s %s\n' % + (ik, self.hopping[ik, 0, i, i].real)) f.write('\n') f.close() @@ -532,16 +618,18 @@ class SumkDFTTools(SumkDFT): r""" Reads the data for transport calculations from the hdf5 archive. """ - thingstoread = ['band_window_optics','velocities_k'] - self.read_input_from_hdf(subgrp=self.transp_data,things_to_read = thingstoread) - thingstoread = ['band_window','lattice_angles','lattice_constants','lattice_type','n_symmetries','rot_symmetries'] - self.read_input_from_hdf(subgrp=self.misc_data,things_to_read = thingstoread) - - + thingstoread = ['band_window_optics', 'velocities_k'] + self.read_input_from_hdf( + subgrp=self.transp_data, things_to_read=thingstoread) + thingstoread = ['band_window', 'lattice_angles', 'lattice_constants', + 'lattice_type', '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 @@ -550,7 +638,7 @@ class SumkDFTTools(SumkDFT): Lattice constants (a, b, c). lattice angles : list of double Lattice angles (:math:`\alpha, \beta, \gamma`). - + Returns ------- vol_c : double @@ -558,20 +646,22 @@ class SumkDFTTools(SumkDFT): 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 + 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): @@ -582,10 +672,10 @@ class SumkDFTTools(SumkDFT): \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 @@ -606,43 +696,52 @@ class SumkDFTTools(SumkDFT): 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. """ - - # Check if wien converter was called and read transport subgroup form hdf file + + # 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 + 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()." - + 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() - + if mpi.is_master_node(): # k-dependent-projections. assert self.k_dep_projection == 1, "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!" + 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") + 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") - n_inequiv_spin_blocks = self.SP + 1 - self.SO # up and down are equivalent if SP = 0 + # up and down are equivalent if SP = 0 + n_inequiv_spin_blocks = self.SP + 1 - self.SO self.directions = directions - dir_to_int = {'x':0, 'y':1, 'z':2} - + 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.Sigma_imp_w[0].mesh]) + self.omega = numpy.array([round(x.real, 12) + for x in self.Sigma_imp_w[0].mesh]) mesh = None mu = self.chemical_potential n_om = len(self.omega) @@ -650,65 +749,81 @@ class SumkDFTTools(SumkDFT): if energy_window is not None: # 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))] + 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: + # 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_w[icrsh].copy() spn = self.spin_block_names[self.corr_shells[icrsh]['SO']] - glist = lambda : [ GfReFreq(indices = inner, window=(self.omega[0], self.omega[-1]),n_points=n_om) for block, inner in self.gf_struct_sumk[icrsh]] - self.Sigma_imp_w[icrsh] = BlockGf(name_list = spn, block_list = glist(),make_copies=False) - for i,g in self.Sigma_imp_w[icrsh]: + glist = lambda: [GfReFreq(indices=inner, window=(self.omega[ + 0], self.omega[-1]), n_points=n_om) for block, inner in self.gf_struct_sumk[icrsh]] + self.Sigma_imp_w[icrsh] = BlockGf( + name_list=spn, block_list=glist(), make_copies=False) + for i, g in self.Sigma_imp_w[icrsh]: for iL in g.indices: for iR in g.indices: for iom in xrange(n_om): - g.data[iom,iL,iR] = Sigma_save[i].data[ioffset+iom,iL,iR] + g.data[iom, iL, iR] = Sigma_save[ + i].data[ioffset + iom, iL, 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 = [energy_window[0]-max(Om_mesh), energy_window[1]+max(Om_mesh), n_om] + self.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] 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]) + 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]), + print "Using n_om = %s points in the energy_window [%s,%s]" % (n_om, self.omega[0], self.omega[-1]), 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=numpy.float_) for direction in self.directions} - + + self.Gamma_w = {direction: numpy.zeros( + (len(self.Om_mesh), n_om), dtype=numpy.float_) for direction in self.directions} + # Sum over all k-points ikarray = numpy.array(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, iw_or_w="w", beta=beta, 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=numpy.complex_) - for isp in range(n_inequiv_spin_blocks)] - + G_w = self.lattice_gf(ik, mu, iw_or_w="w", beta=beta, + 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=numpy.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)); + # 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 xrange(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]))) - - 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) + A_kw[isp][:, :, iw] = -1.0 / (2.0 * numpy.pi * 1j) * ( + A_kw[isp][:, :, iw] - numpy.conjugate(numpy.transpose(A_kw[isp][:, :, iw]))) + + 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: @@ -716,28 +831,30 @@ class SumkDFTTools(SumkDFT): vel_R = copy.deepcopy(self.velocities_k[isp][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 + 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 xrange(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, 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]) - - for direction in self.directions: - 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) + 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, 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]) + + for direction in self.directions: + 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, 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 @@ -751,60 +868,64 @@ class SumkDFTTools(SumkDFT): 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): + 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 + 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) - + 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_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) + A = simps(A_int, self.omega) elif method == 'trapz': # trapezoidal rule for w-grid - A = numpy.trapz(A_int,self.omega) + A = numpy.trapz(A_int, self.omega) else: - # rectangular integration for w-grid (orignal implementation) + # 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) + 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 @@ -814,35 +935,44 @@ class SumkDFTTools(SumkDFT): Seebeck coefficient in each direction. If zero is not present in Om_mesh the Seebeck coefficient 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!" + 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} + + 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} self.seebeck = {direction: numpy.nan for direction in self.directions} - self.optic_cond = {direction: numpy.full((n_q,),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 xrange(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) + 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]): - # Seebeck is overwritten if there is more than one Omega = 0 in Om_mesh - self.seebeck[direction] = - A1[direction][iq] / A0[direction][iq] * 86.17 - self.optic_cond[direction] = beta * A0[direction] * 10700.0 / numpy.pi + # Seebeck is overwritten if there is more than one Omega = + # 0 in Om_mesh + self.seebeck[direction] = - \ + A1[direction][iq] / A0[direction][iq] * 86.17 + self.optic_cond[direction] = beta * \ + A0[direction] * 10700.0 / numpy.pi for iq in xrange(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 "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]) - + return self.optic_cond, self.seebeck - - def fermi_dis(self,w,beta): + def fermi_dis(self, w, beta): r""" Fermi distribution. @@ -855,9 +985,9 @@ class SumkDFTTools(SumkDFT): frequency beta : double inverse temperature - + Returns ------- f : double """ - return 1.0/(numpy.exp(w*beta)+1) + return 1.0 / (numpy.exp(w * beta) + 1) diff --git a/python/symmetry.py b/python/symmetry.py index 2ddec363..4b27b205 100644 --- a/python/symmetry.py +++ b/python/symmetry.py @@ -1,5 +1,5 @@ -################################################################################ +########################################################################## # # TRIQS: a Toolbox for Research in Interacting Quantum Systems # @@ -18,14 +18,16 @@ # You should have received a copy of the GNU General Public License along with # TRIQS. If not, see . # -################################################################################ +########################################################################## -import copy,numpy +import copy +import numpy from types import * from pytriqs.gf.local import * from pytriqs.archive import * import pytriqs.utility.mpi as mpi + class Symmetry: """ This class provides the routines for applying symmetry operations for the k sums. @@ -33,10 +35,10 @@ class Symmetry: rotational matrices for each symmetry operation. """ - def __init__(self, hdf_file, subgroup = None): + def __init__(self, hdf_file, subgroup=None): """ Initialises the class. - + Parameters ---------- hdf_file : string @@ -46,69 +48,80 @@ class Symmetry: the data is stored at the root of the hdf5 archive. """ - assert type(hdf_file) == StringType, "Symmetry: hdf_file must be a filename." + assert type( + hdf_file) == StringType, "Symmetry: hdf_file must be a filename." self.hdf_file = hdf_file - things_to_read = ['n_symm','n_atoms','perm','orbits','SO','SP','time_inv','mat','mat_tinv'] - for it in things_to_read: setattr(self,it,0) + things_to_read = ['n_symm', 'n_atoms', 'perm', + 'orbits', 'SO', 'SP', 'time_inv', 'mat', 'mat_tinv'] + for it in things_to_read: + setattr(self, it, 0) if mpi.is_master_node(): - #Read the stuff on master: - ar = HDFArchive(hdf_file,'r') + # Read the stuff on master: + ar = HDFArchive(hdf_file, 'r') if subgroup is None: ar2 = ar else: ar2 = ar[subgroup] - for it in things_to_read: setattr(self,it,ar2[it]) + for it in things_to_read: + setattr(self, it, ar2[it]) del ar2 del ar # Broadcasting - for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it))) + for it in things_to_read: + setattr(self, it, mpi.bcast(getattr(self, it))) # now define the mapping of orbitals: # self.orb_map[iorb] = jorb gives the permutation of the orbitals as given in the list, when the # permutation of the atoms is done: self.n_orbits = len(self.orbits) - self.orb_map = [ [0 for iorb in range(self.n_orbits)] for i_symm in range(self.n_symm) ] + self.orb_map = [[0 for iorb in range( + self.n_orbits)] for i_symm in range(self.n_symm)] for i_symm in range(self.n_symm): for iorb in range(self.n_orbits): srch = copy.deepcopy(self.orbits[iorb]) - srch['atom'] = self.perm[i_symm][self.orbits[iorb]['atom']-1] + srch['atom'] = self.perm[i_symm][self.orbits[iorb]['atom'] - 1] self.orb_map[i_symm][iorb] = self.orbits.index(srch) - - def symmetrize(self,obj): + def symmetrize(self, obj): """ Symmetrizes a given object. - + Parameters ---------- obj : list object to symmetrize. It has to be given as list, where its length is determined by the number of equivalent members of the object. Two types of objects are supported: - + - BlockGf : list of Green's functions, - Matrices : The format is taken from density matrices as obtained from Green's functions (DictType). - + Returns ------- symm_obj : list Symmetrized object, of the same type as input object. """ - assert isinstance(obj,list), "symmetrize: obj has to be a list of objects." - assert len(obj) == self.n_orbits, "symmetrize: obj has to be a list of the same length as defined in the init." + assert isinstance( + obj, list), "symmetrize: obj has to be a list of objects." + assert len( + obj) == self.n_orbits, "symmetrize: obj has to be a list of the same length as defined in the init." - if isinstance(obj[0],BlockGf): - symm_obj = [ obj[i].copy() for i in range(len(obj)) ] # here the result is stored, it is a BlockGf! - for iorb in range(self.n_orbits): symm_obj[iorb].zero() # set to zero + if isinstance(obj[0], BlockGf): + # here the result is stored, it is a BlockGf! + symm_obj = [obj[i].copy() for i in range(len(obj))] + for iorb in range(self.n_orbits): + symm_obj[iorb].zero() # set to zero else: - # if not a BlockGf, we assume it is a matrix (density matrix), has to be complex since self.mat is complex! - symm_obj = [ copy.deepcopy(obj[i]) for i in range(len(obj)) ] + # if not a BlockGf, we assume it is a matrix (density matrix), has + # to be complex since self.mat is complex! + symm_obj = [copy.deepcopy(obj[i]) for i in range(len(obj))] for iorb in range(self.n_orbits): if type(symm_obj[iorb]) == DictType: - for ii in symm_obj[iorb]: symm_obj[iorb][ii] *= 0.0 + for ii in symm_obj[iorb]: + symm_obj[iorb][ii] *= 0.0 else: symm_obj[iorb] *= 0.0 @@ -118,12 +131,15 @@ class Symmetry: dim = self.orbits[iorb]['dim'] jorb = self.orb_map[i_symm][iorb] - if isinstance(obj[0],BlockGf): + if isinstance(obj[0], BlockGf): tmp = obj[iorb].copy() - if self.time_inv[i_symm]: tmp << tmp.transpose() - for bname,gf in tmp: tmp[bname].from_L_G_R(self.mat[i_symm][iorb],tmp[bname],self.mat[i_symm][iorb].conjugate().transpose()) - tmp *= 1.0/self.n_symm + if self.time_inv[i_symm]: + tmp << tmp.transpose() + for bname, gf in tmp: + tmp[bname].from_L_G_R(self.mat[i_symm][iorb], tmp[bname], self.mat[ + i_symm][iorb].conjugate().transpose()) + tmp *= 1.0 / self.n_symm symm_obj[jorb] += tmp else: @@ -131,17 +147,17 @@ class Symmetry: if type(obj[iorb]) == DictType: for ii in obj[iorb]: if self.time_inv[i_symm] == 0: - symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb][ii]), + symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[i_symm][iorb], obj[iorb][ii]), self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm else: - symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb][ii].conjugate()), + symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[i_symm][iorb], obj[iorb][ii].conjugate()), self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm else: if self.time_inv[i_symm] == 0: - symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb]), + symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[i_symm][iorb], obj[iorb]), self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm else: - symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[i_symm][iorb],obj[iorb].conjugate()), + symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[i_symm][iorb], obj[iorb].conjugate()), self.mat[i_symm][iorb].conjugate().transpose()) / self.n_symm # Markus: This does not what it is supposed to do, check how this should work (keep for now) diff --git a/python/trans_basis.py b/python/trans_basis.py index 1ef8eba2..d21ab66d 100644 --- a/python/trans_basis.py +++ b/python/trans_basis.py @@ -6,6 +6,7 @@ import pytriqs.utility.mpi as mpi import numpy import copy + class TransBasis: """ Computates rotations into a new basis, using the condition that a given property is diagonal in the new basis. @@ -14,19 +15,19 @@ class TransBasis: def __init__(self, SK=None, hdf_datafile=None): """ Initialization of the class. There are two ways to do so: - + - existing SumkLDA class : when you have an existing SumkLDA instance - from hdf5 archive : when you want to use data from hdf5 archive - + Giving the class instance overrides giving the string for the hdf5 archive. - + Parameters ---------- SK : class SumkLDA, optional Existing instance of SumkLDA class. hdf5_datafile : string, optional Name of hdf5 archive to be used. - + """ if SK is None: @@ -35,68 +36,70 @@ class TransBasis: mpi.report("trans_basis: give SK instance or HDF filename!") return 0 - Converter = Wien2kConverter(filename=hdf_datafile,repacking=False) + Converter = Wien2kConverter(filename=hdf_datafile, repacking=False) Converter.convert_dft_input() del Converter - self.SK = SumkDFT(hdf_file=hdf_datafile+'.h5',use_dft_blocks=False) + self.SK = SumkDFT(hdf_file=hdf_datafile + + '.h5', use_dft_blocks=False) else: self.SK = SK self.T = copy.deepcopy(self.SK.T[0]) self.w = numpy.identity(SK.corr_shells[0]['dim']) - - def calculate_diagonalisation_matrix(self, prop_to_be_diagonal = 'eal'): + def calculate_diagonalisation_matrix(self, prop_to_be_diagonal='eal'): """ Calculates the diagonalisation matrix w, and stores it as member of the class. - + Parameters ---------- prop_to_be_diagonal : string, optional Defines the property to be diagonalized. - + - 'eal' : local hamiltonian (i.e. crystal field) - 'dm' : local density matrix - + Returns ------- wsqr : double Measure for the degree of rotation done by the diagonalisation. wsqr=1 means no rotation. - + """ if prop_to_be_diagonal == 'eal': prop = self.SK.eff_atomic_levels()[0] elif prop_to_be_diagonal == 'dm': - prop = self.SK.density_matrix(method = 'using_point_integration')[0] + prop = self.SK.density_matrix(method='using_point_integration')[0] else: - mpi.report("trans_basis: not a valid quantitiy to be diagonal. Choices are 'eal' or 'dm'.") + mpi.report( + "trans_basis: not a valid quantitiy to be diagonal. Choices are 'eal' or 'dm'.") return 0 if self.SK.SO == 0: - self.eig,self.w = numpy.linalg.eigh(prop['up']) + self.eig, self.w = numpy.linalg.eigh(prop['up']) # calculate new Transformation matrix - self.T = numpy.dot(self.T.transpose().conjugate(),self.w).conjugate().transpose() + self.T = numpy.dot(self.T.transpose().conjugate(), + self.w).conjugate().transpose() else: - self.eig,self.w = numpy.linalg.eigh(prop['ud']) + self.eig, self.w = numpy.linalg.eigh(prop['ud']) # calculate new Transformation matrix - self.T = numpy.dot(self.T.transpose().conjugate(),self.w).conjugate().transpose() + self.T = numpy.dot(self.T.transpose().conjugate(), + self.w).conjugate().transpose() # measure for the 'unity' of the transformation: - wsqr = sum(abs(self.w.diagonal())**2)/self.w.diagonal().size + wsqr = sum(abs(self.w.diagonal())**2) / self.w.diagonal().size return wsqr - - def rotate_gf(self,gf_to_rot): + def rotate_gf(self, gf_to_rot): """ Uses the diagonalisation matrix w to rotate a given GF into the new basis. - + Parameters ---------- gf_to_rot : BlockGf Green's function block to rotate. - + Returns ------- gfreturn : BlockGf @@ -104,86 +107,90 @@ class TransBasis: """ # build a full GF - gfrotated = BlockGf( name_block_generator = [ (block,GfImFreq(indices = inner, mesh = gf_to_rot.mesh)) for block,inner in self.SK.gf_struct_sumk[0] ], make_copies = False) + gfrotated = BlockGf(name_block_generator=[(block, GfImFreq( + indices=inner, mesh=gf_to_rot.mesh)) for block, inner in self.SK.gf_struct_sumk[0]], make_copies=False) # transform the CTQMC blocks to the full matrix: - ish = self.SK.corr_to_inequiv[0] # ish is the index of the inequivalent shell corresponding to icrsh + # ish is the index of the inequivalent shell corresponding to icrsh + ish = self.SK.corr_to_inequiv[0] for block, inner in self.gf_struct_solver[ish].iteritems(): for ind1 in inner: for ind2 in inner: - gfrotated[self.SK.solver_to_sumk_block[ish][block]][ind1,ind2] << gf_to_rot[block][ind1,ind2] + gfrotated[self.SK.solver_to_sumk_block[ish][block]][ + ind1, ind2] << gf_to_rot[block][ind1, ind2] # Rotate using the matrix w - for bname,gf in gfrotated: - gfrotated[bname].from_L_G_R(self.w.transpose().conjugate(),gfrotated[bname],self.w) + for bname, gf in gfrotated: + gfrotated[bname].from_L_G_R( + self.w.transpose().conjugate(), gfrotated[bname], self.w) gfreturn = gf_to_rot.copy() # Put back into CTQMC basis: for block, inner in self.gf_struct_solver[ish].iteritems(): for ind1 in inner: for ind2 in inner: - gfreturn[block][ind1,ind2] << gfrotated[self.SK.solver_to_sumk_block[0][block]][ind1,ind2] + gfreturn[block][ind1, ind2] << gfrotated[ + self.SK.solver_to_sumk_block[0][block]][ind1, ind2] return gfreturn - def write_trans_file(self, filename): """ Writes the new transformation T into a file readable by dmftproj. By that, the requested quantity is diagonal already at input. - + Parameters ---------- filename : string Name of the file where the transformation is stored. """ - f = open(filename,'w') + f = open(filename, 'w') Tnew = self.T.conjugate() dim = self.SK.corr_shells[0]['dim'] if self.SK.SO == 0: - for i in range(dim): - st = '' - for k in range(dim): - st += " %9.6f"%(Tnew[i,k].real) - st += " %9.6f"%(Tnew[i,k].imag) - for k in range(2*dim): - st += " 0.0" + for i in range(dim): + st = '' + for k in range(dim): + st += " %9.6f" % (Tnew[i, k].real) + st += " %9.6f" % (Tnew[i, k].imag) + for k in range(2 * dim): + st += " 0.0" - if i < (dim-1): - f.write("%s\n"%(st)) - else: - st1 = st.replace(' ','*',1) - f.write("%s\n"%(st1)) + if i < (dim - 1): + f.write("%s\n" % (st)) + else: + st1 = st.replace(' ', '*', 1) + f.write("%s\n" % (st1)) - for i in range(dim): - st = '' - for k in range(2*dim): - st += " 0.0" - for k in range(dim): - st += " %9.6f"%(Tnew[i,k].real) - st += " %9.6f"%(Tnew[i,k].imag) + for i in range(dim): + st = '' + for k in range(2 * dim): + st += " 0.0" + for k in range(dim): + st += " %9.6f" % (Tnew[i, k].real) + st += " %9.6f" % (Tnew[i, k].imag) - if i < (dim-1): - f.write("%s\n"%(st)) - else: - st1 = st.replace(' ','*',1) - f.write("%s\n"%(st1)) + if i < (dim - 1): + f.write("%s\n" % (st)) + else: + st1 = st.replace(' ', '*', 1) + f.write("%s\n" % (st1)) else: for i in range(dim): - st = '' - for k in range(dim): - st += " %9.6f"%(Tnew[i,k].real) - st += " %9.6f"%(Tnew[i,k].imag) + st = '' + for k in range(dim): + st += " %9.6f" % (Tnew[i, k].real) + st += " %9.6f" % (Tnew[i, k].imag) - if i < (dim-1): - f.write("%s\n"%(st)) - else: - st1 = st.replace(' ','*',1) - f.write("%s\n"%(st1)) + if i < (dim - 1): + f.write("%s\n" % (st)) + else: + st1 = st.replace(' ', '*', 1) + f.write("%s\n" % (st1)) f.close() diff --git a/python/update_archive.py b/python/update_archive.py index d961d0d8..c2af8c69 100644 --- a/python/update_archive.py +++ b/python/update_archive.py @@ -5,8 +5,8 @@ import numpy import subprocess if len(sys.argv) < 2: - print "Usage: python update_archive.py old_archive [v1.0|v1.2]" - sys.exit() + print "Usage: python update_archive.py old_archive [v1.0|v1.2]" + sys.exit() print """ This script is an attempt to update your archive to TRIQS 1.2. @@ -15,13 +15,16 @@ Please keep a copy of your old archive as this script is If you encounter any problem please report it on github! """ + def convert_shells(shells): shell_entries = ['atom', 'sort', 'l', 'dim'] - return [ {name: int(val) for name, val in zip(shell_entries, shells[ish])} for ish in range(len(shells)) ] + return [{name: int(val) for name, val in zip(shell_entries, shells[ish])} for ish in range(len(shells))] + def convert_corr_shells(corr_shells): corr_shell_entries = ['atom', 'sort', 'l', 'dim', 'SO', 'irep'] - return [ {name: int(val) for name, val in zip(corr_shell_entries, corr_shells[icrsh])} for icrsh in range(len(corr_shells)) ] + return [{name: int(val) for name, val in zip(corr_shell_entries, corr_shells[icrsh])} for icrsh in range(len(corr_shells))] + def det_shell_equivalence(corr_shells): corr_to_inequiv = [0 for i in range(len(corr_shells))] @@ -29,20 +32,20 @@ def det_shell_equivalence(corr_shells): n_inequiv_shells = 1 if len(corr_shells) > 1: - inequiv_sort = [ corr_shells[0]['sort'] ] - inequiv_l = [ corr_shells[0]['l'] ] - for i in range(len(corr_shells)-1): + inequiv_sort = [corr_shells[0]['sort']] + inequiv_l = [corr_shells[0]['l']] + for i in range(len(corr_shells) - 1): is_equiv = False for j in range(n_inequiv_shells): - if (inequiv_sort[j]==corr_shells[i+1]['sort']) and (inequiv_l[j]==corr_shells[i+1]['l']): + if (inequiv_sort[j] == corr_shells[i + 1]['sort']) and (inequiv_l[j] == corr_shells[i + 1]['l']): is_equiv = True - corr_to_inequiv[i+1] = j - if is_equiv==False: - corr_to_inequiv[i+1] = n_inequiv_shells + corr_to_inequiv[i + 1] = j + if is_equiv == False: + corr_to_inequiv[i + 1] = n_inequiv_shells n_inequiv_shells += 1 - inequiv_sort.append( corr_shells[i+1]['sort'] ) - inequiv_l.append( corr_shells[i+1]['l'] ) - inequiv_to_corr.append( i+1 ) + inequiv_sort.append(corr_shells[i + 1]['sort']) + inequiv_l.append(corr_shells[i + 1]['l']) + inequiv_to_corr.append(i + 1) return n_inequiv_shells, corr_to_inequiv, inequiv_to_corr @@ -50,48 +53,50 @@ def det_shell_equivalence(corr_shells): ### Main ### filename = sys.argv[1] -if len(sys.argv) > 2: +if len(sys.argv) > 2: from_v = sys.argv[2] -else: # Assume updating an old v1.0 script +else: # Assume updating an old v1.0 script from_v = 'v1.0' A = h5py.File(filename) # Rename groups -old_to_new = {'SumK_LDA':'dft_input', 'SumK_LDA_ParProj':'dft_parproj_input', - 'SymmCorr':'dft_symmcorr_input', 'SymmPar':'dft_symmpar_input', 'SumK_LDA_Bands':'dft_bands_input'} +old_to_new = {'SumK_LDA': 'dft_input', 'SumK_LDA_ParProj': 'dft_parproj_input', + 'SymmCorr': 'dft_symmcorr_input', 'SymmPar': 'dft_symmpar_input', 'SumK_LDA_Bands': 'dft_bands_input'} for old, new in old_to_new.iteritems(): - if old not in A.keys(): continue - print "Changing %s to %s ..."%(old, new) - A.copy(old,new) + if old not in A.keys(): + continue + print "Changing %s to %s ..." % (old, new) + A.copy(old, new) del(A[old]) # Move output items from dft_input to user_data -move_to_output = ['chemical_potential','dc_imp','dc_energ'] +move_to_output = ['chemical_potential', 'dc_imp', 'dc_energ'] for obj in move_to_output: if obj in A['dft_input'].keys(): - if 'user_data' not in A: A.create_group('user_data') - print "Moving %s to user_data ..."%obj - A.copy('dft_input/'+obj,'user_data/'+obj) - del(A['dft_input'][obj]) + if 'user_data' not in A: + A.create_group('user_data') + print "Moving %s to user_data ..." % obj + A.copy('dft_input/' + obj, 'user_data/' + obj) + del(A['dft_input'][obj]) # Delete obsolete quantities -to_delete = ['gf_struct_solver','map_inv','map','deg_shells','h_field'] +to_delete = ['gf_struct_solver', 'map_inv', 'map', 'deg_shells', 'h_field'] for obj in to_delete: if obj in A['dft_input'].keys(): - del(A['dft_input'][obj]) + del(A['dft_input'][obj]) if from_v == 'v1.0': # Update shells and corr_shells to list of dicts - shells_old = HDFArchive(filename,'r')['dft_input']['shells'] - corr_shells_old = HDFArchive(filename,'r')['dft_input']['corr_shells'] + shells_old = HDFArchive(filename, 'r')['dft_input']['shells'] + corr_shells_old = HDFArchive(filename, 'r')['dft_input']['corr_shells'] shells = convert_shells(shells_old) corr_shells = convert_corr_shells(corr_shells_old) del(A['dft_input']['shells']) del(A['dft_input']['corr_shells']) A.close() # Need to use HDFArchive for the following - HDFArchive(filename,'a')['dft_input']['shells'] = shells - HDFArchive(filename,'a')['dft_input']['corr_shells'] = corr_shells + HDFArchive(filename, 'a')['dft_input']['shells'] = shells + HDFArchive(filename, 'a')['dft_input']['corr_shells'] = corr_shells A = h5py.File(filename) # Add shell equivalency quantities @@ -102,32 +107,36 @@ if 'n_inequiv_shells' not in A['dft_input']: A['dft_input']['inequiv_to_corr'] = equiv_shell_info[2] # Rename variables -groups = ['dft_symmcorr_input','dft_symmpar_input'] +groups = ['dft_symmcorr_input', 'dft_symmpar_input'] for group in groups: - if group not in A.keys(): continue - if 'n_s' not in A[group]: continue + if group not in A.keys(): + continue + if 'n_s' not in A[group]: + continue print "Changing n_s to n_symm ..." - A[group].move('n_s','n_symm') + A[group].move('n_s', 'n_symm') # Convert orbits to list of dicts - orbits_old = HDFArchive(filename,'r')[group]['orbits'] + orbits_old = HDFArchive(filename, 'r')[group]['orbits'] orbits = convert_corr_shells(orbits_old) del(A[group]['orbits']) A.close() - HDFArchive(filename,'a')[group]['orbits'] = orbits + HDFArchive(filename, 'a')[group]['orbits'] = orbits A = h5py.File(filename) groups = ['dft_parproj_input'] for group in groups: - if group not in A.keys(): continue - if 'proj_mat_pc' not in A[group]: continue + if group not in A.keys(): + continue + if 'proj_mat_pc' not in A[group]: + continue print "Changing proj_mat_pc to proj_mat_all ..." - A[group].move('proj_mat_pc','proj_mat_all') + A[group].move('proj_mat_pc', 'proj_mat_all') A.close() # Repack to reclaim disk space -retcode = subprocess.call(["h5repack","-i%s"%filename, "-otemphgfrt.h5"]) +retcode = subprocess.call(["h5repack", "-i%s" % filename, "-otemphgfrt.h5"]) if retcode != 0: print "h5repack failed!" else: - subprocess.call(["mv","-f","temphgfrt.h5","%s"%filename]) + subprocess.call(["mv", "-f", "temphgfrt.h5", "%s" % filename])