LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-07-11 05:43:48 +02:00
Code Issues Releases Wiki Activity

Changes to old interface files to comply with new gf_struct

Browse Source 
Minor tidy-up too.
This commit is contained in:
Priyanka Seth 2014-09-22 19:24:33 +02:00
parent f803c13285
commit 906398894a
6 changed files with 405 additions and 503 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

62
python/converters/hk_converter.py
View File

@ -28,10 +28,10 @@ import string
from math import sqrt from math import sqrt
def Read_Fortran_File (filename): def read_fortran_file (filename):
""" Returns a generator that yields all numbers in the Fortran file as float, one by one""" """ Returns a generator that yields all numbers in the Fortran file as float, one by one"""
import os.path import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename
for line in open(filename,'r') : for line in open(filename,'r') :
for x in line.replace('D','E').replace('(',' ').replace(')',' ').replace(',',' ').split() : for x in line.replace('D','E').replace('(',' ').replace(')',' ').replace(',',' ').split() :
yield string.atof(x) yield string.atof(x)
@ -40,22 +40,18 @@ def Read_Fortran_File (filename):
class HkConverter: class HkConverter:
""" """
Conversion from general H(k) file to an hdf5 file, that can be used as input for the SumK_LDA class. Conversion from general H(k) file to an hdf5 file that can be used as input for the SumK_LDA class.
""" """
def __init__(self, hk_file, hdf_file, lda_subgrp = 'SumK_LDA', symm_subgrp = 'SymmCorr', repacking = False): def __init__(self, hk_file, hdf_file, lda_subgrp = 'SumK_LDA', symm_subgrp = 'SymmCorr', repacking = False):
""" """
Init of the class. Variable Filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so Init of the class.
on. on.
""" """
assert type(hk_file)==StringType,"hk_file must be a filename" assert type(hk_file)==StringType,"hk_file must be a filename"
self.hdf_file = hdf_file self.hdf_file = hdf_file
self.lda_file = hk_file self.lda_file = hk_file
#self.Symm_file = Filename+'.symqmc'
#self.Parproj_file = Filename+'.parproj'
#self.Symmpar_file = Filename+'.sympar'
#self.Band_file = Filename+'.outband'
self.lda_subgrp = lda_subgrp self.lda_subgrp = lda_subgrp
self.symm_subgrp = symm_subgrp self.symm_subgrp = symm_subgrp
@ -72,12 +68,12 @@ class HkConverter:
""" """
if not (mpi.is_master_node()): return # do it only on master: # Read and write only on the master node
if not (mpi.is_master_node()): return
mpi.report("Reading input from %s..."%self.lda_file) mpi.report("Reading input from %s..."%self.lda_file)
# Read and write only on Master!!!
# 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 = Read_Fortran_File(self.lda_file) R = read_fortran_file(self.lda_file)
try: try:
energy_unit = 1.0 # the energy conversion factor is 1.0, we assume eV in files 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 n_k = int(R.next()) # read the number of k points
@ -101,7 +97,6 @@ class HkConverter:
self.inequiv_shells(corr_shells) # determine the number of inequivalent correlated shells, has to be known for further reading... self.inequiv_shells(corr_shells) # determine the number of inequivalent correlated shells, has to be known for further reading...
use_rotations = 0 use_rotations = 0
rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)] rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)]
rot_mat_time_inv = [0 for i in range(n_corr_shells)] rot_mat_time_inv = [0 for i in range(n_corr_shells)]
@ -109,16 +104,13 @@ class HkConverter:
# Representative representations are read from file # Representative representations are read from file
n_reps = [1 for i in range(self.n_inequiv_corr_shells)] n_reps = [1 for i in range(self.n_inequiv_corr_shells)]
dim_reps = [0 for i in range(self.n_inequiv_corr_shells)] dim_reps = [0 for i in range(self.n_inequiv_corr_shells)]
T = []
for icrsh in range(self.n_inequiv_corr_shells): for icrsh in range(self.n_inequiv_corr_shells):
n_reps[icrsh] = int(R.next()) # number of representatives ("subsets"), e.g. t2g and eg n_reps[icrsh] = int(R.next()) # number of representatives ("subsets"), e.g. t2g and eg
dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])] # dimensions of the subsets dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])] # dimensions of the subsets
# The transformation matrix: # The transformation matrix:
# it is of dimension 2l+1, it is taken to be standard d (as in Wien2k) # is of dimension 2l+1, it is taken to be standard d (as in Wien2k)
T = []
for icrsh in range(self.n_inequiv_corr_shells):
#for ish in xrange(self.N_inequiv_corr_shells):
ll = 2*corr_shells[self.invshellmap[icrsh]][2]+1 ll = 2*corr_shells[self.invshellmap[icrsh]][2]+1
lmax = ll * (corr_shells[self.invshellmap[icrsh]][4] + 1) lmax = ll * (corr_shells[self.invshellmap[icrsh]][4] + 1)
T.append(numpy.zeros([lmax,lmax],numpy.complex_)) T.append(numpy.zeros([lmax,lmax],numpy.complex_))
@ -131,27 +123,21 @@ class HkConverter:
# Spin blocks to be read: # Spin blocks to be read:
n_spin_blocks = SP + 1 - SO # number of spins to read for Norbs and Ham, NOT Projectors n_spin_blocs = SP + 1 - SO # number of spins to read for Norbs and Ham, NOT Projectors
# define the number of N_Orbitals for all k points: it is the number of total bands and independent of k! # define the number of n_orbitals for all k points: it is the number of total bands and independent of k!
n_orb = sum([ shells[ish][3] for ish in range(n_shells)]) n_orb = sum([ shells[ish][3] for ish in range(n_shells) ])
#n_orbitals = [ [n_orb for isp in range(n_spin_blocks)] for ik in xrange(n_k)] n_orbitals = numpy.ones([n_k,n_spin_blocs],numpy.int) * n_orb
n_orbitals = numpy.ones([n_k,n_spin_blocks],numpy.int) * n_orb
#print N_Orbitals
# Initialise the projectors: # Initialise the projectors:
#proj_mat = [ [ [numpy.zeros([corr_shells[icrsh][3], n_orbitals[ik][isp]], numpy.complex_) proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
# for icrsh in range (n_corr_shells)]
# for isp in range(n_spin_blocks)]
# for ik in range(n_k) ]
proj_mat = numpy.zeros([n_k,n_spin_blocks,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
# Read the projectors from the file: # Read the projectors from the file:
for ik in xrange(n_k): for ik in xrange(n_k):
for icrsh in range(n_corr_shells): for icrsh in range(n_corr_shells):
for isp in range(n_spin_blocks): for isp in range(n_spin_blocs):
# calculate the offset: # calculate the offset:
offset = 0 offset = 0
@ -169,9 +155,7 @@ class HkConverter:
# now define the arrays for weights and hopping ... # now define the arrays for weights and hopping ...
bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation
#hopping = [ [numpy.zeros([n_orbitals[ik][isp],n_orbitals[ik][isp]],numpy.complex_) hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_)
# for isp in range(n_spin_blocks)] for ik in xrange(n_k) ]
hopping = numpy.zeros([n_k,n_spin_blocks,max(n_orbitals),max(n_orbitals)],numpy.complex_)
if (weights_in_file): if (weights_in_file):
# weights in the file # weights in the file
@ -181,11 +165,9 @@ class HkConverter:
sm = sum(bz_weights) sm = sum(bz_weights)
bz_weights[:] /= sm bz_weights[:] /= sm
# Grab the H # Grab the H
for isp in range(n_spin_blocs):
for ik in xrange(n_k) : for ik in xrange(n_k) :
for isp in range(n_spin_blocks):
no = n_orbitals[ik,isp] no = n_orbitals[ik,isp]
if (first_real_part_matrix): if (first_real_part_matrix):
@ -220,24 +202,22 @@ class HkConverter:
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: # keep some things that we need for reading parproj:
self.n_shells = n_shells self.n_shells = n_shells
self.shells = shells self.shells = shells
self.n_corr_shells = n_corr_shells self.n_corr_shells = n_corr_shells
self.corr_shells = corr_shells self.corr_shells = corr_shells
self.n_spin_blocks = n_spin_blocks self.n_spin_blocs = n_spin_blocs
self.n_orbitals = n_orbitals self.n_orbitals = n_orbitals
self.n_k = n_k self.n_k = n_k
self.SO = SO self.SO = SO
self.SP = SP self.SP = SP
self.energy_unit = energy_unit self.energy_unit = energy_unit
except StopIteration : # a more explicit error if the file is corrupted. except StopIteration : # a more explicit error if the file is corrupted.
raise "SumK_LDA : reading file HMLT_file failed!" raise "HK Converter : reading file lda_file failed!"
R.close() R.close()
#print Proj_Mat[0]
#----------------------------------------- #-----------------------------------------
# Store the input into HDF5: # Store the input into HDF5:
ar = HDFArchive(self.hdf_file,'a') ar = HDFArchive(self.hdf_file,'a')
@ -276,7 +256,7 @@ class HkConverter:
def __repack(self): def __repack(self):
"""Calls the h5repack routine, in order to reduce the file size of the hdf5 archive. """Calls the h5repack routine, in order to reduce the file size of the hdf5 archive.
Should only be used BEFORE the first invokation of HDF_Archive in the program, otherwise Should only be used BEFORE the first invokation of HDFArchive in the program, otherwise
the hdf5 linking is broken!!!""" the hdf5 linking is broken!!!"""
import subprocess import subprocess

75
python/converters/wien2k_converter.py
View File

@ -30,7 +30,7 @@ import string
def read_fortran_file (filename): def read_fortran_file (filename):
""" Returns a generator that yields all numbers in the Fortran file as float, one by one""" """ Returns a generator that yields all numbers in the Fortran file as float, one by one"""
import os.path import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename
for line in open(filename,'r') : for line in open(filename,'r') :
for x in line.replace('D','E').split() : for x in line.replace('D','E').split() :
yield string.atof(x) yield string.atof(x)
@ -39,13 +39,12 @@ def read_fortran_file (filename):
class Wien2kConverter: class Wien2kConverter:
""" """
Conversion from Wien2k output to an hdf5 file, that can be used as input for the SumkLDA class. Conversion from Wien2k output to an hdf5 file that can be used as input for the SumkLDA class.
""" """
def __init__(self, filename, lda_subgrp = 'SumK_LDA', symm_subgrp = 'SymmCorr', repacking = False): def __init__(self, filename, lda_subgrp = 'SumK_LDA', symm_subgrp = 'SymmCorr', repacking = False):
""" """
Init of the class. Variable filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so Init of the class. Variable filename gives the root of all filenames, e.g. case.ctqmcout, case.h5, and so on.
on.
""" """
assert type(filename)==StringType,"LDA_file must be a filename" assert type(filename)==StringType,"LDA_file must be a filename"
@ -71,10 +70,10 @@ class Wien2kConverter:
""" """
if not (mpi.is_master_node()): return # do it only on master: # Read and write only on the master node
if not (mpi.is_master_node()): return
mpi.report("Reading input from %s..."%self.lda_file) mpi.report("Reading input from %s..."%self.lda_file)
# Read and write only on Master!!!
# 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 = read_fortran_file(self.lda_file) R = read_fortran_file(self.lda_file)
try: try:
@ -91,7 +90,6 @@ class Wien2kConverter:
n_shells = int(R.next()) # number of shells (e.g. Fe d, As p, O p) in the unit cell, n_shells = int(R.next()) # number of shells (e.g. Fe d, As p, O p) in the unit cell,
# corresponds to index R in formulas # corresponds to index R in formulas
shells = [ [ int(R.next()) for i in range(4) ] for icrsh in range(n_shells) ] # reads iatom, sort, l, dim shells = [ [ int(R.next()) for i in range(4) ] for icrsh in range(n_shells) ] # reads iatom, sort, l, dim
#shells = numpy.array(shells)
n_corr_shells = int(R.next()) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell, n_corr_shells = int(R.next()) # number of corr. shells (e.g. Fe d, Ce f) in the unit cell,
# corresponds to index R in formulas # corresponds to index R in formulas
@ -99,7 +97,6 @@ class Wien2kConverter:
corr_shells = [ [ int(R.next()) for i in range(6) ] for icrsh in range(n_corr_shells) ] # reads iatom, sort, l, dim, SO flag, irep corr_shells = [ [ int(R.next()) for i in range(6) ] for icrsh in range(n_corr_shells) ] # reads iatom, sort, l, dim, SO flag, irep
self.inequiv_shells(corr_shells) # determine the number of inequivalent correlated shells, has to be known for further reading... self.inequiv_shells(corr_shells) # determine the number of inequivalent correlated shells, has to be known for further reading...
#corr_shells = numpy.array(corr_shells)
use_rotations = 1 use_rotations = 1
rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)] rot_mat = [numpy.identity(corr_shells[icrsh][3],numpy.complex_) for icrsh in xrange(n_corr_shells)]
@ -120,7 +117,7 @@ class Wien2kConverter:
# Read here the infos for the transformation of the basis: # Read here the info for the transformation of the basis:
n_reps = [1 for i in range(self.n_inequiv_corr_shells)] n_reps = [1 for i in range(self.n_inequiv_corr_shells)]
dim_reps = [0 for i in range(self.n_inequiv_corr_shells)] dim_reps = [0 for i in range(self.n_inequiv_corr_shells)]
T = [] T = []
@ -129,9 +126,7 @@ class Wien2kConverter:
dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])] # dimensions of the subsets dim_reps[icrsh] = [int(R.next()) for i in range(n_reps[icrsh])] # dimensions of the subsets
# The transformation matrix: # The transformation matrix:
# it is of dimension 2l+1, if no SO, and 2*(2l+1) with SO!! # is of dimension 2l+1 without SO, and 2*(2l+1) with SO!
#T = []
#for ish in xrange(self.n_inequiv_corr_shells):
ll = 2*corr_shells[self.invshellmap[icrsh]][2]+1 ll = 2*corr_shells[self.invshellmap[icrsh]][2]+1
lmax = ll * (corr_shells[self.invshellmap[icrsh]][4] + 1) lmax = ll * (corr_shells[self.invshellmap[icrsh]][4] + 1)
T.append(numpy.zeros([lmax,lmax],numpy.complex_)) T.append(numpy.zeros([lmax,lmax],numpy.complex_))
@ -151,19 +146,12 @@ class Wien2kConverter:
# read the list of n_orbitals for all k points # 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)
#n_orbitals = [ [0 for isp in range(n_spin_blocs)] for ik in xrange(n_k)]
for isp in range(n_spin_blocs): for isp in range(n_spin_blocs):
for ik in xrange(n_k): for ik in xrange(n_k):
#n_orbitals[ik][isp] = int(R.next())
n_orbitals[ik,isp] = int(R.next()) n_orbitals[ik,isp] = int(R.next())
#print n_orbitals
# Initialise the projectors: # Initialise the projectors:
#proj_mat = [ [ [numpy.zeros([corr_shells[icrsh][3], n_orbitals[ik][isp]], numpy.complex_)
# for icrsh in range (n_corr_shells)]
# for isp in range(n_spin_blocs)]
# for ik in range(n_k) ]
proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_) proj_mat = numpy.zeros([n_k,n_spin_blocs,n_corr_shells,max(numpy.array(corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
@ -175,20 +163,16 @@ class Wien2kConverter:
for isp in range(n_spin_blocs): for isp in range(n_spin_blocs):
for i in xrange(no): for i in xrange(no):
for j in xrange(n_orbitals[ik][isp]): for j in xrange(n_orbitals[ik][isp]):
#proj_mat[ik][isp][icrsh][i,j] = R.next()
proj_mat[ik,isp,icrsh,i,j] = R.next() proj_mat[ik,isp,icrsh,i,j] = R.next()
# now Imag part: # now Imag part:
for isp in range(n_spin_blocs): for isp in range(n_spin_blocs):
for i in xrange(no): for i in xrange(no):
for j in xrange(n_orbitals[ik][isp]): for j in xrange(n_orbitals[ik][isp]):
#proj_mat[ik][isp][icrsh][i,j] += 1j * R.next()
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 ... # now define the arrays for weights and hopping ...
bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation bz_weights = numpy.ones([n_k],numpy.float_)/ float(n_k) # w(k_index), default normalisation
#hopping = [ [numpy.zeros([n_orbitals[ik][isp],n_orbitals[ik][isp]],numpy.complex_)
# for isp in range(n_spin_blocs)] for ik in xrange(n_k) ]
hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_) hopping = numpy.zeros([n_k,n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_)
# weights in the file # weights in the file
@ -202,12 +186,11 @@ class Wien2kConverter:
# we use now the convention of a DIAGONAL Hamiltonian!!!! # we use now the convention of a DIAGONAL Hamiltonian!!!!
for isp in range(n_spin_blocs): for isp in range(n_spin_blocs):
for ik in xrange(n_k) : for ik in xrange(n_k) :
no = n_orbitals[ik][isp] no = n_orbitals[ik,isp]
for i in xrange(no): for i in xrange(no):
#hopping[ik][isp][i,i] = R.next() * energy_unit
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: # keep some things that we need for reading parproj:
self.n_shells = n_shells self.n_shells = n_shells
self.shells = shells self.shells = shells
self.n_corr_shells = n_corr_shells self.n_corr_shells = n_corr_shells
@ -219,12 +202,10 @@ class Wien2kConverter:
self.SP = SP self.SP = SP
self.energy_unit = energy_unit self.energy_unit = energy_unit
except StopIteration : # a more explicit error if the file is corrupted. except StopIteration : # a more explicit error if the file is corrupted.
raise "SumkLDA : reading file HMLT_file failed!" raise "Wien2k_converter : reading file lda_file failed!"
R.close() R.close()
#print proj_mat[0]
#----------------------------------------- #-----------------------------------------
# Store the input into HDF5: # Store the input into HDF5:
ar = HDFArchive(self.hdf_file,'a') ar = HDFArchive(self.hdf_file,'a')
@ -279,25 +260,17 @@ class Wien2kConverter:
for isp in range(self.n_spin_blocs) ] for isp in range(self.n_spin_blocs) ]
R = read_fortran_file(self.parproj_file) R = read_fortran_file(self.parproj_file)
#try:
n_parproj = [int(R.next()) for i in range(self.n_shells)] n_parproj = [int(R.next()) for i in range(self.n_shells)]
n_parproj = numpy.array(n_parproj) n_parproj = numpy.array(n_parproj)
# Initialise P, here a double list of matrices: # Initialise P, here a double list of matrices:
#proj_mat_pc = [ [ [ [numpy.zeros([self.shells[ish][3], self.n_orbitals[ik][isp]], numpy.complex_)
# for ir in range(n_parproj[ish])]
# for ish in range (self.n_shells) ]
# for isp in range(self.n_spin_blocs) ]
# for ik in range(self.n_k) ]
proj_mat_pc = numpy.zeros([self.n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(self.n_orbitals)],numpy.complex_) proj_mat_pc = numpy.zeros([self.n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(self.n_orbitals)],numpy.complex_)
rot_mat_all = [numpy.identity(self.shells[ish][3],numpy.complex_) for ish in xrange(self.n_shells)] rot_mat_all = [numpy.identity(self.shells[ish][3],numpy.complex_) for ish in xrange(self.n_shells)]
rot_mat_all_time_inv = [0 for i 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): for ish in range(self.n_shells):
#print ish
# read first the projectors for this orbital: # read first the projectors for this orbital:
for ik in xrange(self.n_k): for ik in xrange(self.n_k):
for ir in range(n_parproj[ish]): for ir in range(n_parproj[ish]):
@ -337,8 +310,6 @@ class Wien2kConverter:
if (self.SP): if (self.SP):
rot_mat_all_time_inv[ish] = int(R.next()) rot_mat_all_time_inv[ish] = int(R.next())
#except StopIteration : # a more explicit error if the file is corrupted.
# raise "Wien2kConverter: reading file for Projectors failed!"
R.close() R.close()
#----------------------------------------- #-----------------------------------------
@ -378,10 +349,6 @@ class Wien2kConverter:
n_orbitals[ik,isp] = int(R.next()) n_orbitals[ik,isp] = int(R.next())
# Initialise the projectors: # Initialise the projectors:
#proj_mat = [ [ [numpy.zeros([self.corr_shells[icrsh][3], n_orbitals[ik][isp]], numpy.complex_)
# for icrsh in range (self.n_corr_shells)]
# for isp in range(self.n_spin_blocs)]
# for ik in range(n_k) ]
proj_mat = numpy.zeros([n_k,self.n_spin_blocs,self.n_corr_shells,max(numpy.array(self.corr_shells)[:,3]),max(n_orbitals)],numpy.complex_) proj_mat = numpy.zeros([n_k,self.n_spin_blocs,self.n_corr_shells,max(numpy.array(self.corr_shells)[:,3]),max(n_orbitals)],numpy.complex_)
# Read the projectors from the file: # Read the projectors from the file:
@ -399,8 +366,6 @@ class Wien2kConverter:
for j in xrange(n_orbitals[ik,isp]): for j in xrange(n_orbitals[ik,isp]):
proj_mat[ik,isp,icrsh,i,j] += 1j * R.next() proj_mat[ik,isp,icrsh,i,j] += 1j * R.next()
#hopping = [ [numpy.zeros([n_orbitals[ik][isp],n_orbitals[ik][isp]],numpy.complex_)
# for isp in range(self.n_spin_blocs)] for ik in xrange(n_k) ]
hopping = numpy.zeros([n_k,self.n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_) hopping = numpy.zeros([n_k,self.n_spin_blocs,max(n_orbitals),max(n_orbitals)],numpy.complex_)
# Grab the H # Grab the H
@ -416,11 +381,6 @@ class Wien2kConverter:
n_parproj = numpy.array(n_parproj) n_parproj = numpy.array(n_parproj)
# Initialise P, here a double list of matrices: # Initialise P, here a double list of matrices:
#proj_mat_pc = [ [ [ [numpy.zeros([self.shells[ish][3], n_orbitals[ik][isp]], numpy.complex_)
# for ir in range(n_parproj[ish])]
# for ish in range (self.n_shells) ]
# for isp in range(self.n_spin_blocs) ]
# for ik in range(n_k) ]
proj_mat_pc = numpy.zeros([n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(n_orbitals)],numpy.complex_) proj_mat_pc = numpy.zeros([n_k,self.n_spin_blocs,self.n_shells,max(n_parproj),max(numpy.array(self.shells)[:,3]),max(n_orbitals)],numpy.complex_)
@ -439,7 +399,7 @@ class Wien2kConverter:
proj_mat_pc[ik,isp,ish,ir,i,j] += 1j * R.next() proj_mat_pc[ik,isp,ish,ir,i,j] += 1j * R.next()
except StopIteration : # a more explicit error if the file is corrupted. except StopIteration : # a more explicit error if the file is corrupted.
raise "SumkLDA : reading file HMLT_file failed!" raise "Wien2k_converter : reading file band_file failed!"
R.close() R.close()
# reading done! # reading done!
@ -448,18 +408,11 @@ class Wien2kConverter:
# Store the input into HDF5: # Store the input into HDF5:
ar = HDFArchive(self.hdf_file,'a') ar = HDFArchive(self.hdf_file,'a')
if not (self.bands_subgrp in ar): ar.create_group(self.bands_subgrp) 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. # The subgroup containing the data. If it does not exist, it is created.
# If it exists, the data is overwritten!!! # If it exists, the data is overwritten!!!
thingstowrite = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc'] thingstowrite = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc']
for it in thingstowrite: exec "ar['%s']['%s'] = %s"%(self.bands_subgrp,it,it) for it in thingstowrite: exec "ar['%s']['%s'] = %s"%(self.bands_subgrp,it,it)
#ar[self.bands_subgrp]['n_k'] = n_k
#ar[self.bands_subgrp]['n_orbitals'] = n_orbitals
#ar[self.bands_subgrp]['proj_mat'] = proj_mat
#self.proj_mat = proj_mat
#self.n_orbitals = n_orbitals
#self.n_k = n_k
#self.hopping = hopping
del ar del ar
@ -501,7 +454,7 @@ class Wien2kConverter:
mat[in_s][orb][i,j] += 1j * R.next() # imaginary part mat[in_s][orb][i,j] += 1j * R.next() # imaginary part
# determine the inequivalent shells: # determine the inequivalent shells:
#SHOULD BE FINALLY REMOVED, PUT IT FOR ALL ORBITALS!!!!! #SHOULD BE FINALLY REMOVED, PUT IT FOR ALL ORBITALS!!!!! (PS: FIXME?)
#self.inequiv_shells(orbits) #self.inequiv_shells(orbits)
mat_tinv = [numpy.identity(orbits[orb][3],numpy.complex_) mat_tinv = [numpy.identity(orbits[orb][3],numpy.complex_)
for orb in range(n_orbits)] for orb in range(n_orbits)]
@ -519,7 +472,7 @@ class Wien2kConverter:
except StopIteration : # a more explicit error if the file is corrupted. except StopIteration : # a more explicit error if the file is corrupted.
raise "Symmetry : reading file failed!" raise "Wien2k_converter : reading file symm_file failed!"
R.close() R.close()

176
python/sumk_lda.py
View File

@ -21,18 +21,16 @@
################################################################################ ################################################################################
from types import * from types import *
#from pytriqs.applications.dft.symmetry import *
from symmetry import * from symmetry import *
import numpy import numpy
import pytriqs.utility.dichotomy as dichotomy import pytriqs.utility.dichotomy as dichotomy
from pytriqs.gf.local import * from pytriqs.gf.local import *
#from pytriqs.applications.impurity_solvers.operators import *
from pytriqs.operators import *
from pytriqs.archive import * from pytriqs.archive import *
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
from math import cos,sin from math import cos,sin
# FIXME PS: These two aren't used it seems
from pytriqs.operators.operators2 import *
import string, pickle import string, pickle
class SumkLDA: class SumkLDA:
@ -65,20 +63,16 @@ class SumkLDA:
'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping'] 'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping']
optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells'] optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells']
#ar=HDFArchive(self.hdf_file,'a')
#del ar
self.retval = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read,optional_things=optional_things) self.retval = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read,optional_things=optional_things)
#ar=HDFArchive(self.hdf_file,'a')
#del ar
if (self.SO) and (abs(self.h_field)>0.000001): if (self.SO) and (abs(self.h_field)>0.000001):
self.h_field=0.0 self.h_field=0.0
mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!!") mpi.report("For SO, the external magnetic field is not implemented, setting it to 0!!")
self.inequiv_shells(self.corr_shells) # determine the number of inequivalent correlated shells # determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells)
# and related maps (self.shellmap, self.invshellmap)
self.inequiv_shells(self.corr_shells)
# field to convert block_names to indices # field to convert block_names to indices
self.names_to_ind = [{}, {}] self.names_to_ind = [{}, {}]
@ -92,9 +86,10 @@ class SumkLDA:
if not (self.retval['gf_struct_solver']): if not (self.retval['gf_struct_solver']):
# No gf_struct was stored in HDF, so first set a standard one: # No gf_struct was stored in HDF, so first set a standard one:
self.gf_struct_solver = [ [ (al, range( self.corr_shells[self.invshellmap[i]][3]) ) self.gf_struct_solver = [ dict([ (al, range(self.corr_shells[self.invshellmap[i]][3]) )
for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]] ] for al in self.block_names[self.corr_shells[self.invshellmap[i]][4]] ])
for i in xrange(self.n_inequiv_corr_shells) ] for i in range(self.n_inequiv_corr_shells)
]
self.map = [ {} for i in xrange(self.n_inequiv_corr_shells) ] self.map = [ {} for i in xrange(self.n_inequiv_corr_shells) ]
self.map_inv = [ {} for i in xrange(self.n_inequiv_corr_shells) ] self.map_inv = [ {} for i in xrange(self.n_inequiv_corr_shells) ]
for i in xrange(self.n_inequiv_corr_shells): for i in xrange(self.n_inequiv_corr_shells):
@ -116,7 +111,7 @@ class SumkLDA:
#mpi.report("Do the init for symm:") #mpi.report("Do the init for symm:")
self.Symm_corr = Symmetry(hdf_file,subgroup=self.symm_corr_data) self.Symm_corr = Symmetry(hdf_file,subgroup=self.symm_corr_data)
# determine the smallest blocs, if wanted: # Analyse the block structure and determine the smallest blocs, if desired
if (use_lda_blocks): dm=self.analyse_BS() if (use_lda_blocks): dm=self.analyse_BS()
@ -235,8 +230,7 @@ class SumkLDA:
gf_rotated = gf_to_rotate.copy() gf_rotated = gf_to_rotate.copy()
if (direction=='toGlobal'): if (direction=='toGlobal'):
#if (self.rot_mat_time_inv[icrsh]==1): gf_rotated <<= gf_rotated.transpose()
#gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate())
if ((self.rot_mat_time_inv[icrsh]==1) and (self.SO)): if ((self.rot_mat_time_inv[icrsh]==1) and (self.SO)):
gf_rotated <<= gf_rotated.transpose() gf_rotated <<= gf_rotated.transpose()
gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate(),gf_rotated,self.rot_mat[icrsh].transpose()) gf_rotated.from_L_G_R(self.rot_mat[icrsh].conjugate(),gf_rotated,self.rot_mat[icrsh].transpose())
@ -244,6 +238,7 @@ class SumkLDA:
gf_rotated.from_L_G_R(self.rot_mat[icrsh],gf_rotated,self.rot_mat[icrsh].conjugate().transpose()) gf_rotated.from_L_G_R(self.rot_mat[icrsh],gf_rotated,self.rot_mat[icrsh].conjugate().transpose())
elif (direction=='toLocal'): elif (direction=='toLocal'):
if ((self.rot_mat_time_inv[icrsh]==1)and(self.SO)): if ((self.rot_mat_time_inv[icrsh]==1)and(self.SO)):
gf_rotated <<= gf_rotated.transpose() gf_rotated <<= gf_rotated.transpose()
gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate()) gf_rotated.from_L_G_R(self.rot_mat[icrsh].transpose(),gf_rotated,self.rot_mat[icrsh].conjugate())
@ -266,25 +261,52 @@ class SumkLDA:
stmp = self.add_dc() stmp = self.add_dc()
beta = self.Sigma_imp[0].mesh.beta #override beta if Sigma is present beta = self.Sigma_imp[0].mesh.beta #override beta if Sigma is present
if (self.Gupf is None): #
# first setting up of Gupf # if (self.Gupf is None):
BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] # # first setting up of Gupf
gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] # BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
a_list = [a for a,al in gf_struct] # gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
if (with_Sigma): #take the mesh from Sigma if necessary # a_list = [a for a,al in gf_struct]
glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct] # if (with_Sigma): #take the mesh from Sigma if necessary
# glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct]
# else:
# glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct]
# self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
# self.Gupf.zero()
# self.Gupf_id = self.Gupf.copy()
# self.Gupf_id <<= iOmega_n
#
# GFsize = [ gf.N1 for sig,gf in self.Gupf]
# unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib]
# for ib in range(self.n_spin_blocks_gf[self.SO]) ] )
#
# if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)):
# BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
# gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
# a_list = [a for a,al in gf_struct]
# if (with_Sigma):
# glist = lambda : [ GfImFreq(indices = al, mesh = self.Sigma_imp[0].mesh) for a,al in gf_struct]
# else:
# glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct]
# self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
# self.Gupf.zero()
# self.Gupf_id = self.Gupf.copy()
# self.Gupf_id <<= iOmega_n
#
###
# FIXME PS Remove commented out code above if this works
# Do we need to set up Gupf?
set_up_Gupf = False
if self.Gupf == None:
set_up_Gupf = True
else: else:
glist = lambda : [ GfImFreq(indices = al, beta = beta) for a,al in gf_struct]
self.Gupf = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
self.Gupf.zero()
self.Gupf_id = self.Gupf.copy()
self.Gupf_id <<= iOmega_n
GFsize = [ gf.N1 for sig,gf in self.Gupf] GFsize = [ gf.N1 for sig,gf in self.Gupf]
unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib] unchangedsize = all( [ self.n_orbitals[ik,ntoi[bln[ib]]]==GFsize[ib]
for ib in range(self.n_spin_blocks_gf[self.SO]) ] ) for ib in range(self.n_spin_blocks_gf[self.SO]) ] )
if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)): set_up_Gupf = True
if ((not unchangedsize)or(self.Gupf.mesh.beta!=beta)): # Set up Gupf
if set_up_Gupf:
BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ] BS = [ range(self.n_orbitals[ik,ntoi[ib]]) for ib in bln ]
gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ] gf_struct = [ (bln[ib], BS[ib]) for ib in range(self.n_spin_blocks_gf[self.SO]) ]
a_list = [a for a,al in gf_struct] a_list = [a for a,al in gf_struct]
@ -296,10 +318,9 @@ class SumkLDA:
self.Gupf.zero() self.Gupf.zero()
self.Gupf_id = self.Gupf.copy() self.Gupf_id = self.Gupf.copy()
self.Gupf_id <<= iOmega_n self.Gupf_id <<= iOmega_n
###
idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln] idmat = [numpy.identity(self.n_orbitals[ik,ntoi[bl]],numpy.complex_) for bl in bln]
#for ibl in range(self.n_spin_blocks_gf[self.SO]): mupat[ibl] *= mu
self.Gupf <<= self.Gupf_id self.Gupf <<= self.Gupf_id
M = copy.deepcopy(idmat) M = copy.deepcopy(idmat)
@ -455,7 +476,7 @@ class SumkLDA:
def analyse_BS(self, threshold = 0.00001, include_shells = None, dm = None): def analyse_BS(self, threshold = 0.00001, include_shells = None, dm = None):
""" Determines the Greens function block structure from the simple point integration""" """ Determines the Green function block structure from simple point integration."""
if (dm==None): dm = self.simple_point_dens_mat() if (dm==None): dm = self.simple_point_dens_mat()
@ -464,8 +485,7 @@ class SumkLDA:
if include_shells is None: include_shells=range(self.n_inequiv_corr_shells) if include_shells is None: include_shells=range(self.n_inequiv_corr_shells)
for ish in include_shells: for ish in include_shells:
#self.gf_struct_solver.append([]) self.gf_struct_solver[ish] = {}
self.gf_struct_solver[ish] = []
gf_struct_temp = [] gf_struct_temp = []
a_list = [a for a,al in self.gf_struct_corr[self.invshellmap[ish]] ] a_list = [a for a,al in self.gf_struct_corr[self.invshellmap[ish]] ]
@ -496,16 +516,16 @@ class SumkLDA:
for i in range(NBlocs): for i in range(NBlocs):
blocs[i].sort() blocs[i].sort()
self.gf_struct_solver[ish].append( ('%s%s'%(a,i),range(len(blocs[i]))) ) self.gf_struct_solver[ish].update( [('%s_%s'%(a,i),range(len(blocs[i])))] )
gf_struct_temp.append( ('%s%s'%(a,i),blocs[i]) ) gf_struct_temp.append( ('%s_%s'%(a,i),blocs[i]) )
# map is the mapping of the blocs from the SK blocs to the CTQMC blocs: # map is the mapping of the blocs from the SK blocs to the CTQMC blocs:
self.map[ish][a] = range(len(dmbool)) self.map[ish][a] = range(len(dmbool))
for ibl in range(NBlocs): for ibl in range(NBlocs):
for j in range(len(blocs[ibl])): for j in range(len(blocs[ibl])):
self.map[ish][a][blocs[ibl][j]] = '%s%s'%(a,ibl) self.map[ish][a][blocs[ibl][j]] = '%s_%s'%(a,ibl)
self.map_inv[ish]['%s%s'%(a,ibl)] = a self.map_inv[ish]['%s_%s'%(a,ibl)] = a
# now calculate degeneracies of orbitals: # now calculate degeneracies of orbitals:
@ -638,29 +658,12 @@ class SumkLDA:
def set_lichtenstein_dc(self,Sigma_imp):
"""Sets a double counting term according to Lichtenstein et al. PRL2001"""
assert 0,"Lichtenstein DC not supported any more!"
def set_dc(self,dens_mat,U_interact,J_hund,orb=0,use_dc_formula=0,use_val=None): def set_dc(self,dens_mat,U_interact,J_hund,orb=0,use_dc_formula=0,use_val=None):
"""Sets the double counting term for inequiv orbital orb """Sets the double counting term for inequiv orbital orb:
use_dc_formula=0: LDA+U FLL double counting, use_dc_formula=1: Held's formula. use_dc_formula=0: LDA+U FLL double counting,
use_dc_formula=2: AMF use_dc_formula=1: Held's formula,
Be sure that you use the correct interaction Hamiltonian!""" use_dc_formula=2: AMF.
Be sure that you are using the correct interaction Hamiltonian!"""
#if (not hasattr(self,"dc_imp")): self.__init_dc()
#dm = [ {} for i in xrange(self.n_corr_shells)]
#for i in xrange(self.n_corr_shells):
# l = self.corr_shells[i][3] #*(1+self.corr_shells[i][4])
# for j in xrange(len(self.gf_struct_corr[i])):
# dm[i]['%s'%self.gf_struct_corr[i][j][0]] = numpy.zeros([l,l],numpy.float_)
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
@ -677,13 +680,10 @@ class SumkLDA:
blname = self.gf_struct_corr[icrsh][j][0] blname = self.gf_struct_corr[icrsh][j][0]
Ncr[blname] = 0.0 Ncr[blname] = 0.0
for a,al in self.gf_struct_solver[iorb]: for a,al in self.gf_struct_solver[iorb].iteritems():
#for bl in self.map[iorb][blname]:
bl = self.map_inv[iorb][a] bl = self.map_inv[iorb][a]
#print 'bl, valiue = ',bl,dens_mat[a].real.trace()
Ncr[bl] += dens_mat[a].real.trace() Ncr[bl] += dens_mat[a].real.trace()
#print 'Ncr=',Ncr
M = self.corr_shells[icrsh][3] M = self.corr_shells[icrsh][3]
Ncrtot = 0.0 Ncrtot = 0.0
@ -702,26 +702,22 @@ class SumkLDA:
if (use_val is None): if (use_val is None):
if (use_dc_formula==0): if (use_dc_formula==0): # FLL
self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0) self.dc_energ[icrsh] = U_interact / 2.0 * Ncrtot * (Ncrtot-1.0)
for bl in a_list: for bl in a_list:
Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5) Uav = U_interact*(Ncrtot-0.5) - J_hund*(Ncr[bl] - 0.5)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][bl] *= Uav
self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0) self.dc_energ[icrsh] -= J_hund / 2.0 * (Ncr[bl]) * (Ncr[bl]-1.0)
mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals())
elif (use_dc_formula==1): elif (use_dc_formula==1): # Held's formula, with U_interact the interorbital onsite interaction
self.dc_energ[icrsh] = (U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) / 2.0 * Ncrtot * (Ncrtot-1.0) self.dc_energ[icrsh] = (U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) / 2.0 * Ncrtot * (Ncrtot-1.0)
#self.dc_energ[icrsh] = (U_interact + J_hund * (2.0-(M-1)) / (2*M-1) ) / 2.0 * Ncrtot * (Ncrtot-1.0)
for bl in a_list: for bl in a_list:
# Held's formula, with U_interact the interorbital onsite interaction
Uav =(U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) * (Ncrtot-0.5) Uav =(U_interact + (M-1)*(U_interact-2.0*J_hund) + (M-1)*(U_interact-3.0*J_hund))/(2*M-1) * (Ncrtot-0.5)
#Uav = (U_interact + J_hund * (2.0-(M-1)) / (2*M-1) ) * (Ncrtot-0.5)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][bl] *= Uav
mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals()) mpi.report("DC for shell %(icrsh)i and block %(bl)s = %(Uav)f"%locals())
elif (use_dc_formula==2): elif (use_dc_formula==2): # AMF
self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot self.dc_energ[icrsh] = 0.5 * U_interact * Ncrtot * Ncrtot
for bl in a_list: for bl in a_list:
# AMF
Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M) Uav = U_interact*(Ncrtot - Ncr[bl]/M) - J_hund * (Ncr[bl] - Ncr[bl]/M)
self.dc_imp[icrsh][bl] *= Uav self.dc_imp[icrsh][bl] *= Uav
self.dc_energ[icrsh] -= (U_interact + (M-1)*J_hund)/M * 0.5 * Ncr[bl] * Ncr[bl] self.dc_energ[icrsh] -= (U_interact + (M-1)*J_hund)/M * 0.5 * Ncr[bl] * Ncr[bl]
@ -747,7 +743,7 @@ class SumkLDA:
def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01): def find_dc(self,orb,guess,dens_mat,dens_req=None,precision=0.01):
"""searches for DC in order to fulfill charge neutrality. """Searches for DC in order to fulfill charge neutrality.
If dens_req is given, then DC is set such that the LOCAL charge of orbital If dens_req is given, then DC is set such that the LOCAL charge of orbital
orb coincides with dens_req.""" orb coincides with dens_req."""
@ -806,20 +802,18 @@ class SumkLDA:
for blname in self.map[s]: for blname in self.map[s]:
cnt[blname] = 0 cnt[blname] = 0
for a,al in self.gf_struct_solver[s]: for a,al in self.gf_struct_solver[s].iteritems():
blname = self.map_inv[s][a] blname = self.map_inv[s][a]
map_ind[a] = range(len(al)) map_ind[a] = range(len(al))
for i in al: for i in al:
map_ind[a][i] = cnt[blname] map_ind[a][i] = cnt[blname]
cnt[blname]+=1 cnt[blname]+=1
for bl, orblist in self.gf_struct_solver[s].iteritems():
for ibl in range(len(self.gf_struct_solver[s])): for i in range(len(orblist)):
for i in range(len(self.gf_struct_solver[s][ibl][1])): for j in range(len(orblist)):
for j in range(len(self.gf_struct_solver[s][ibl][1])): ind1 = orblist[i]
bl = self.gf_struct_solver[s][ibl][0] ind2 = orblist[j]
ind1 = self.gf_struct_solver[s][ibl][1][i]
ind2 = self.gf_struct_solver[s][ibl][1][j]
ind1_imp = map_ind[bl][ind1] ind1_imp = map_ind[bl][ind1]
ind2_imp = map_ind[bl][ind2] ind2_imp = map_ind[bl][ind2]
self.Sigma_imp[icrsh][self.map_inv[s][bl]][ind1_imp,ind2_imp] <<= Sigma_imp[s][bl][ind1,ind2] self.Sigma_imp[icrsh][self.map_inv[s][bl]][ind1_imp,ind2_imp] <<= Sigma_imp[s][bl][ind1,ind2]
@ -984,7 +978,7 @@ class SumkLDA:
def extract_G_loc(self, mu=None, with_Sigma = True): def extract_G_loc(self, mu=None, with_Sigma = True):
""" """
extracts the local downfolded Green function at the chemical potential of the class. extracts the local downfolded Green function at the chemical potential of the class.
At the end, the local G is rotated from the gloabl coordinate system to the local system. At the end, the local G is rotated from the global coordinate system to the local system.
if with_Sigma = False: Sigma is not included => non-interacting local GF if with_Sigma = False: Sigma is not included => non-interacting local GF
""" """
@ -1023,7 +1017,7 @@ class SumkLDA:
for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toLocal') for sig,gf in Gloc[icrsh]: Gloc[icrsh][sig] <<= self.rotloc(icrsh,gf,direction='toLocal')
# transform to CTQMC blocks: # transform to CTQMC blocks:
Glocret = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Gloc[0].mesh)) for a,al in self.gf_struct_solver[i] ], Glocret = [ BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = Gloc[0].mesh)) for a,al in self.gf_struct_solver[i].iteritems() ],
make_copies = False) for i in xrange(self.n_inequiv_corr_shells) ] make_copies = False) for i in xrange(self.n_inequiv_corr_shells) ]
for ish in xrange(self.n_inequiv_corr_shells): for ish in xrange(self.n_inequiv_corr_shells):
@ -1033,20 +1027,18 @@ class SumkLDA:
for blname in self.map[ish]: for blname in self.map[ish]:
cnt[blname] = 0 cnt[blname] = 0
for a,al in self.gf_struct_solver[ish]: for a,al in self.gf_struct_solver[ish].iteritems():
blname = self.map_inv[ish][a] blname = self.map_inv[ish][a]
map_ind[a] = range(len(al)) map_ind[a] = range(len(al))
for i in al: for i in al:
map_ind[a][i] = cnt[blname] map_ind[a][i] = cnt[blname]
cnt[blname]+=1 cnt[blname]+=1
for bl, orblist in self.gf_struct_solver[ish].iteritems():
for ibl in range(len(self.gf_struct_solver[ish])): for i in range(len(orblist)):
for i in range(len(self.gf_struct_solver[ish][ibl][1])): for j in range(len(orblist)):
for j in range(len(self.gf_struct_solver[ish][ibl][1])): ind1 = orblist[i]
bl = self.gf_struct_solver[ish][ibl][0] ind2 = orblist[j]
ind1 = self.gf_struct_solver[ish][ibl][1][i]
ind2 = self.gf_struct_solver[ish][ibl][1][j]
ind1_imp = map_ind[bl][ind1] ind1_imp = map_ind[bl][ind1]
ind2_imp = map_ind[bl][ind2] ind2_imp = map_ind[bl][ind2]
Glocret[ish][bl][ind1,ind2] <<= Gloc[self.invshellmap[ish]][self.map_inv[ish][bl]][ind1_imp,ind2_imp] Glocret[ish][bl][ind1,ind2] <<= Gloc[self.invshellmap[ish]][self.map_inv[ish][bl]][ind1_imp,ind2_imp]

36
python/sumk_lda_tools.py
View File

@ -24,23 +24,17 @@ from types import *
import numpy import numpy
import pytriqs.utility.dichotomy as dichotomy import pytriqs.utility.dichotomy as dichotomy
from pytriqs.gf.local import * from pytriqs.gf.local import *
#from pytriqs.applications.impurity_solvers.operators import *
from pytriqs.operators import * from pytriqs.operators import *
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
from datetime import datetime from datetime import datetime
#from pytriqs.applications.dft.symmetry import *
#from pytriqs.applications.dft.sumk_lda import SumkLDA
from symmetry import * from symmetry import *
from sumk_lda import SumkLDA from sumk_lda import SumkLDA
import string import string
def read_fortran_file (filename): def read_fortran_file (filename):
""" Returns a generator that yields all numbers in the Fortran file as float, one by one""" """ Returns a generator that yields all numbers in the Fortran file as float, one by one"""
import os.path import os.path
if not(os.path.exists(filename)) : raise IOError, "File %s does not exists"%filename if not(os.path.exists(filename)) : raise IOError, "File %s does not exist."%filename
for line in open(filename,'r') : for line in open(filename,'r') :
for x in line.replace('D','E').split() : for x in line.replace('D','E').split() :
yield string.atof(x) yield string.atof(x)
@ -516,13 +510,15 @@ class SumkLDATools(SumkLDA):
def constr_Sigma_real_axis(self, filename, hdf=True, hdf_dataset='SigmaReFreq',n_om=0,orb=0, tol_mesh=1e-6): def constr_Sigma_real_axis(self, filename, hdf=True, hdf_dataset='SigmaReFreq',n_om=0,orb=0, tol_mesh=1e-6):
"""Uses Data from files to construct Sigma (or GF) on the real axis.""" """Uses Data from files to construct Sigma (or GF) on the real axis."""
if not hdf: if not hdf: # then read sigma from text files
# read sigma from text files
#first get the mesh out of one of the files: # first get the mesh out of any one of the files:
if (len(self.gf_struct_solver[orb][0][1])==1): bl = self.gf_struct_solver[orb].items()[0][0] # block name
Fname = filename+'_'+self.gf_struct_solver[orb][0][0]+'.dat' ol = self.gf_struct_solver[orb].items()[0][1] # list of orbital indices
if (len(ol)==1): # if blocks are of size one
Fname = filename+'_'+bl+'.dat'
else: else:
Fname = filename+'_'+self.gf_struct_solver[orb][0][0]+'/'+str(self.gf_struct_solver[orb][0][1][0])+'_'+str(self.gf_struct_solver[orb][0][1][0])+'.dat' Fname = filename+'_'+bl+'/'+str(ol[0])+'_'+str(ol[0])+'.dat'
R = read_fortran_file(Fname) R = read_fortran_file(Fname)
mesh = numpy.zeros([n_om],numpy.float_) mesh = numpy.zeros([n_om],numpy.float_)
@ -542,12 +538,11 @@ class SumkLDATools(SumkLDA):
assert abs(i*bin+mesh[0]-mesh[i]) < tol_mesh, 'constr_Sigma_ME: real-axis mesh is non-uniform!' assert abs(i*bin+mesh[0]-mesh[i]) < tol_mesh, 'constr_Sigma_ME: real-axis mesh is non-uniform!'
# construct Sigma # construct Sigma
a_list = [a for a,al in self.gf_struct_solver[orb]] a_list = [a for a,al in self.gf_struct_solver[orb].iteritems()]
glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[n_om-1]),n_points=n_om) for a,al in self.gf_struct_solver[orb]] glist = lambda : [ GfReFreq(indices = al, window=(mesh[0],mesh[n_om-1]),n_points=n_om) for a,al in self.gf_struct_solver[orb].iteritems()]
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
#read Sigma #read Sigma
for i,g in SigmaME: for i,g in SigmaME:
mesh=[w for w in g.mesh] mesh=[w for w in g.mesh]
for iL in g.indices: for iL in g.indices:
@ -568,9 +563,8 @@ class SumkLDATools(SumkLDA):
R.close() R.close()
else: else: # read sigma from hdf
# read sigma from hdf
omega_min=0.0 omega_min=0.0
omega_max=0.0 omega_max=0.0
n_om=0 n_om=0
@ -588,8 +582,8 @@ class SumkLDATools(SumkLDA):
mpi.barrier() mpi.barrier()
# construct Sigma on other nodes # construct Sigma on other nodes
if (not mpi.is_master_node()): if (not mpi.is_master_node()):
a_list = [a for a,al in self.gf_struct_solver[orb]] a_list = [a for a,al in self.gf_struct_solver[orb].iteritems()]
glist = lambda : [ GfReFreq(indices = al, window=(omega_min,omega_max),n_points=n_om) for a,al in self.gf_struct_solver[orb]] glist = lambda : [ GfReFreq(indices = al, window=(omega_min,omega_max),n_points=n_om) for a,al in self.gf_struct_solver[orb].iteritems()]
SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False) SigmaME = BlockGf(name_list = a_list, block_list = glist(),make_copies=False)
# pass SigmaME to other nodes # pass SigmaME to other nodes
SigmaME = mpi.bcast(SigmaME) SigmaME = mpi.bcast(SigmaME)
@ -673,5 +667,3 @@ class SumkLDATools(SumkLDA):
for isp in range(len(bln)) ] for isp in range(len(bln)) ]
return dens_mat return dens_mat

10
python/symmetry.py
View File

@ -107,10 +107,6 @@ class Symmetry:
if (isinstance(obj[0],BlockGf)): if (isinstance(obj[0],BlockGf)):
#if l==0:
# symm_obj[jorb] += obj[iorb]
#else:
tmp = obj[iorb].copy() tmp = obj[iorb].copy()
if (self.time_inv[in_s]): tmp <<= tmp.transpose() if (self.time_inv[in_s]): tmp <<= tmp.transpose()
for sig,gf in tmp: tmp[sig].from_L_G_R(self.mat[in_s][iorb],tmp[sig],self.mat[in_s][iorb].conjugate().transpose()) for sig,gf in tmp: tmp[sig].from_L_G_R(self.mat[in_s][iorb],tmp[sig],self.mat[in_s][iorb].conjugate().transpose())
@ -122,9 +118,6 @@ class Symmetry:
if (type(obj[iorb])==DictType): if (type(obj[iorb])==DictType):
for ii in obj[iorb]: for ii in obj[iorb]:
#if (l==0):
# symm_obj[jorb][ii] += obj[iorb][ii]/self.n_s
#else:
if (self.time_inv[in_s]==0): if (self.time_inv[in_s]==0):
symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb][ii]), symm_obj[jorb][ii] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb][ii]),
self.mat[in_s][iorb].conjugate().transpose()) / self.n_s self.mat[in_s][iorb].conjugate().transpose()) / self.n_s
@ -135,9 +128,6 @@ class Symmetry:
else: else:
#if (l==0):
# symm_obj[jorb] += obj[iorb]/self.n_s
#else:
if (self.time_inv[in_s]==0): if (self.time_inv[in_s]==0):
symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb]),self.mat[in_s][iorb].conjugate().transpose()) / self.n_s symm_obj[jorb] += numpy.dot(numpy.dot(self.mat[in_s][iorb],obj[iorb]),self.mat[in_s][iorb].conjugate().transpose()) / self.n_s
else: else:

31
python/trans_basis.py
View File

@ -8,7 +8,7 @@ import copy
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
class TransBasis: class TransBasis:
'''Computates rotations into a new basis, in order to make certain quantities diagonal.''' '''Computates rotations into a new basis in order to make certain quantities diagonal.'''
def __init__(self, SK=None, hdf_datafile=None): def __init__(self, SK=None, hdf_datafile=None):
@ -70,7 +70,7 @@ class TransBasis:
def rotate_gf(self,gf_to_rot): def rotate_gf(self,gf_to_rot):
'''rotates a given GF into the new basis''' '''Rotates a given GF into the new basis.'''
# build a full GF # build a full GF
gfrotated = BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = gf_to_rot.mesh)) for a,al in self.SK.gf_struct_corr[0] ], make_copies = False) gfrotated = BlockGf( name_block_generator = [ (a,GfImFreq(indices = al, mesh = gf_to_rot.mesh)) for a,al in self.SK.gf_struct_corr[0] ], make_copies = False)
@ -78,12 +78,11 @@ class TransBasis:
# transform the CTQMC blocks to the full matrix: # transform the CTQMC blocks to the full matrix:
s = self.SK.shellmap[0] # s is the index of the inequivalent shell corresponding to icrsh s = self.SK.shellmap[0] # s is the index of the inequivalent shell corresponding to icrsh
for ibl in range(len(self.SK.gf_struct_solver[s])): for bl, orblist in self.gf_struct_solver[s].iteritems():
for i in range(len(self.SK.gf_struct_solver[s][ibl][1])): for i in range(len(orblist)):
for j in range(len(self.SK.gf_struct_solver[s][ibl][1])): for j in range(len(orblist)):
bl = self.SK.gf_struct_solver[s][ibl][0] ind1 = orblist[i]
ind1 = self.SK.gf_struct_solver[s][ibl][1][i] ind2 = orblist[j]
ind2 = self.SK.gf_struct_solver[s][ibl][1][j]
gfrotated[self.SK.map_inv[s][bl]][ind1,ind2] <<= gf_to_rot[bl][ind1,ind2] gfrotated[self.SK.map_inv[s][bl]][ind1,ind2] <<= gf_to_rot[bl][ind1,ind2]
# Rotate using the matrix w # Rotate using the matrix w
@ -92,19 +91,18 @@ class TransBasis:
gfreturn = gf_to_rot.copy() gfreturn = gf_to_rot.copy()
# Put back into CTQMC basis: # Put back into CTQMC basis:
for ibl in range(len(self.SK.gf_struct_solver[0])): for bl, orblist in self.gf_struct_solver[s].iteritems():
for i in range(len(self.SK.gf_struct_solver[0][ibl][1])): for i in range(len(orblist)):
for j in range(len(self.SK.gf_struct_solver[0][ibl][1])): for j in range(len(orblist)):
bl = self.SK.gf_struct_solver[0][ibl][0] ind1 = orblist[i]
ind1 = self.SK.gf_struct_solver[0][ibl][1][i] ind2 = orblist[j]
ind2 = self.SK.gf_struct_solver[0][ibl][1][j]
gfreturn[bl][ind1,ind2] <<= gfrotated[self.SK.map_inv[0][bl]][ind1,ind2] gfreturn[bl][ind1,ind2] <<= gfrotated[self.SK.map_inv[0][bl]][ind1,ind2]
return gfreturn return gfreturn
def write_trans_file(self, filename): def write_trans_file(self, filename):
'''writes the new transformation into a file, readable for dmftproj.''' '''Writes the new transformation into a file readable by dmftproj.'''
f=open(filename,'w') f=open(filename,'w')
@ -158,6 +156,3 @@ class TransBasis:
#MPI.report("SO not implemented!") #MPI.report("SO not implemented!")
f.close() f.close()