LCPQ
/
dft_tools
mirror of https://github.com/triqs/dft_tools
3
0
Fork 0
Code Issues Releases Wiki Activity

Minor clean up, pep-ified to allow doc compilation to run smoothly

This commit is contained in:
Priyanka Seth 2016-05-09 10:19:56 +02:00
parent 841f840df5
commit 390e8564b7
12 changed files with 1663 additions and 1184 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
python/__init__.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
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']

9
python/clear_h5_output.py
View File

@ -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])

8
python/converters/__init__.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
from wien2k_converter import Wien2kConverter
from hk_converter import HkConverter
from wannier90_converter import Wannier90Converter
__all__ =['Wien2kConverter','HkConverter','Wannier90Converter']
__all__ = ['Wien2kConverter', 'HkConverter', 'Wannier90Converter']

52
python/converters/converter_tools.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
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.
@ -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.
@ -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

172
python/converters/hk_converter.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
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.
@ -52,21 +53,22 @@ class HkConverter(ConverterTools):
"""
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.
@ -82,39 +84,55 @@ class HkConverter(ConverterTools):
"""
# Read and write only on the master node
if not (mpi.is_master_node()): return
mpi.report("Reading input from %s..."%self.dft_file)
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)
# 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
# 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
# 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
# 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,
# 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):
# 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
@ -122,29 +140,39 @@ class HkConverter(ConverterTools):
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
# 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 ])
# 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,21 +183,25 @@ 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)
@ -177,27 +209,31 @@ class HkConverter(ConverterTools):
# 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]
if (first_real_part_matrix): # first read all real components for given k, then read imaginary parts
# first read all real components for given k, then read
# imaginary parts
if (first_real_part_matrix):
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()
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
@ -206,25 +242,31 @@ class HkConverter(ConverterTools):
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
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()
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

102
python/converters/wannier90_converter.py
View File

@ -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.
@ -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

486
python/converters/wien2k_converter.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
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',
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):
bands_subgrp='dft_bands_input', misc_subgrp='dft_misc_input',
transp_subgrp='dft_transp_input', repacking=False):
"""
Initialise the class.
@ -64,18 +65,20 @@ class Wien2kConverter(ConverterTools):
"""
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
@ -103,39 +105,55 @@ class Wien2kConverter(ConverterTools):
"""
# Read and write only on the master node
if not (mpi.is_master_node()): return
mpi.report("Reading input from %s..."%self.dft_file)
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)
# 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
# read the number of k points
n_k = int(R.next())
k_dep_projection = 1
SP = int(R.next()) # flag for spin-polarised calculation
SO = int(R.next()) # flag for spin-orbit calculation
# 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,
# 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,
# 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):
# 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)]
@ -143,12 +161,13 @@ class Wien2kConverter(ConverterTools):
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:
@ -156,94 +175,108 @@ class Wien2kConverter(ConverterTools):
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
# 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
# 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_)
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 = [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,34 +327,39 @@ 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()
proj_mat_all[ik, isp, ish,
ir, i, j] += 1j * R.next()
# now read the Density Matrix for this orbital below the energy window:
# 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,94 +390,109 @@ 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()
# 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()
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 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:
@ -442,11 +506,13 @@ class Wien2kConverter(ConverterTools):
"""
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']
@ -462,25 +528,28 @@ class Wien2kConverter(ConverterTools):
if (SP == 0 or SO == 1):
files = [self.bandwin_file]
elif SP == 1:
files = [self.bandwin_file+'up', self.bandwin_file+'dn']
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)
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')
@ -493,7 +562,7 @@ 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:
@ -501,31 +570,36 @@ class Wien2kConverter(ConverterTools):
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,15 +609,16 @@ 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:
@ -554,11 +629,13 @@ class Wien2kConverter(ConverterTools):
"""
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']
@ -574,17 +651,19 @@ class Wien2kConverter(ConverterTools):
if (SP == 0 or SO == 1):
files = [self.pmat_file]
elif SP == 1:
files = [self.pmat_file+'up', self.pmat_file+'dn']
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!
# 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
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
# 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 = []
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

783
python/sumk_dft.py
View File

File diff suppressed because it is too large Load Diff

538
python/sumk_dft_tools.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
import sys
from types import *
import numpy
@ -28,15 +28,15 @@ 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.
"""
@ -46,7 +46,6 @@ class SumkDFTTools(SumkDFT):
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):
"""
@ -84,86 +83,108 @@ 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)
mesh = (om_min, om_max, n_om)
else:
om_min,om_max,n_om = mesh
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):
"""
@ -196,11 +217,14 @@ 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."
@ -209,91 +233,113 @@ 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)
mesh = (om_min, om_max, n_om)
else:
om_min,om_max,n_om = mesh
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) ]
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()
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.
@ -318,104 +364,130 @@ class SumkDFTTools(SumkDFT):
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
# 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:
# 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
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
@ -439,89 +511,103 @@ class SumkDFTTools(SumkDFT):
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,11 +618,13 @@ 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"""
@ -565,14 +653,16 @@ class SumkDFTTools(SumkDFT):
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)
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}
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):
r"""
@ -607,12 +697,15 @@ class SumkDFTTools(SumkDFT):
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()
@ -620,29 +713,35 @@ class SumkDFTTools(SumkDFT):
# 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,8 +749,10 @@ 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
@ -660,55 +761,69 @@ class SumkDFTTools(SumkDFT):
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)
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])))
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)
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,23 +831,25 @@ 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][:])
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
# 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
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])
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.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
@ -758,41 +875,45 @@ class SumkDFTTools(SumkDFT):
Transport coefficient.
"""
if not (mpi.is_master_node()): return
if not (mpi.is_master_node()):
return
assert hasattr(self,'Gamma_w'), "transport_coefficient: Run transport_distribution first or load data from h5!"
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)
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
@ -814,26 +935,36 @@ 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
if not (mpi.is_master_node()):
return
assert hasattr(self,'Gamma_w'), "conductivity_and_seebeck: Run transport_distribution first or load data from h5!"
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])):
@ -841,8 +972,7 @@ class SumkDFTTools(SumkDFT):
return self.optic_cond, self.seebeck
def fermi_dis(self,w,beta):
def fermi_dis(self, w, beta):
r"""
Fermi distribution.
@ -860,4 +990,4 @@ class SumkDFTTools(SumkDFT):
-------
f : double
"""
return 1.0/(numpy.exp(w*beta)+1)
return 1.0 / (numpy.exp(w * beta) + 1)

78
python/symmetry.py
View File

@ -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 <http://www.gnu.org/licenses/>.
#
################################################################################
##########################################################################
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,7 +35,7 @@ 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.
@ -46,39 +48,44 @@ 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.
@ -97,18 +104,24 @@ class Symmetry:
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)

89
python/trans_basis.py
View File

@ -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.
@ -35,19 +36,19 @@ 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.
@ -69,26 +70,28 @@ class TransBasis:
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.
@ -104,29 +107,33 @@ 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
@ -138,7 +145,7 @@ class TransBasis:
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']
@ -147,43 +154,43 @@ class TransBasis:
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 += " %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))
if i < (dim - 1):
f.write("%s\n" % (st))
else:
st1 = st.replace(' ','*',1)
f.write("%s\n"%(st1))
st1 = st.replace(' ', '*', 1)
f.write("%s\n" % (st1))
for i in range(dim):
st = ''
for k in range(2*dim):
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)
st += " %9.6f" % (Tnew[i, k].real)
st += " %9.6f" % (Tnew[i, k].imag)
if i < (dim-1):
f.write("%s\n"%(st))
if i < (dim - 1):
f.write("%s\n" % (st))
else:
st1 = st.replace(' ','*',1)
f.write("%s\n"%(st1))
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 += " %9.6f" % (Tnew[i, k].real)
st += " %9.6f" % (Tnew[i, k].imag)
if i < (dim-1):
f.write("%s\n"%(st))
if i < (dim - 1):
f.write("%s\n" % (st))
else:
st1 = st.replace(' ','*',1)
f.write("%s\n"%(st1))
st1 = st.replace(' ', '*', 1)
f.write("%s\n" % (st1))
f.close()

83
python/update_archive.py
View File

@ -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
@ -57,41 +60,43 @@ else: # Assume updating an old v1.0 script
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)
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])
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])