LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-09-27 20:11:00 +02:00
Code Issues Releases Wiki Activity

Tidy up of trans_basis

Browse Source
This commit is contained in:
Priyanka Seth 2014-11-15 15:30:10 +01:00
parent c23f6fd720
commit 628f774234
1 changed files with 38 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

96
python/trans_basis.py
View File

@ -2,22 +2,21 @@ from pytriqs.applications.dft.sumk_lda import *
from pytriqs.applications.dft.converters import Wien2kConverter from pytriqs.applications.dft.converters import Wien2kConverter
from pytriqs.gf.local.block_gf import BlockGf from pytriqs.gf.local.block_gf import BlockGf
from pytriqs.gf.local.gf_imfreq import GfImFreq from pytriqs.gf.local.gf_imfreq import GfImFreq
import numpy
from pytriqs.archive import * from pytriqs.archive import *
import copy
import pytriqs.utility.mpi as mpi import pytriqs.utility.mpi as mpi
import numpy
import copy
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):
'''Inits the class by reading the input.''' '''Inits the class by reading the input.'''
if SK is None: if SK is None:
# build our own SK instance # build our own SK instance
if hdf_datafile is None: if hdf_datafile is None:
mpi.report("Give SK instance or HDF filename!") mpi.report("trans_basis: give SK instance or HDF filename!")
return 0 return 0
Converter = Wien2kConverter(filename=hdf_datafile,repacking=False) Converter = Wien2kConverter(filename=hdf_datafile,repacking=False)
@ -32,38 +31,26 @@ class TransBasis:
self.w = numpy.identity(SK.corr_shells[0][3]) self.w = numpy.identity(SK.corr_shells[0][3])
def __call__(self, prop_to_be_diagonal = 'eal'): def __call__(self, prop_to_be_diagonal = 'eal'):
'''Calculates the diagonalisation.''' '''Calculates the diagonalisation.'''
if (prop_to_be_diagonal=='eal'): if prop_to_be_diagonal == 'eal':
eal = self.SK.eff_atomic_levels()[0] prop = self.SK.eff_atomic_levels()[0]
elif (prop_to_be_diagonal=='dm'): elif prop_to_be_diagonal == 'dm':
eal = self.SK.simple_point_dens_mat()[0] prop = self.SK.simple_point_dens_mat()[0]
else: else:
mpi.report("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 return 0
if (self.SK.SO==0): if self.SK.SO == 0:
self.eig,self.w = numpy.linalg.eigh(eal['up']) self.eig,self.w = numpy.linalg.eigh(prop['up'])
# calculate new Transformation matrix
# now 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()
#return numpy.dot(self.w.transpose().conjugate(),numpy.dot(eal['up'],self.w))
else: else:
self.eig,self.w = numpy.linalg.eigh(prop['ud'])
self.eig,self.w = numpy.linalg.eigh(eal['ud']) # calculate new Transformation matrix
# now 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()
#MPI.report("SO not implemented yet!")
#return 0
# measure for the 'unity' of the transformation: # 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 return wsqr
@ -76,13 +63,11 @@ class TransBasis:
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: # transform the CTQMC blocks to the full matrix:
s = self.SK.shellmap[0] # s is the index of the inequivalent shell corresponding to icrsh ish = self.SK.shellmap[0] # ish is the index of the inequivalent shell corresponding to icrsh
for block, inner in self.gf_struct_solver[s].iteritems(): for block, inner in self.gf_struct_solver[ish].iteritems():
for i in range(len(inner)): for ind1 in inner:
for j in range(len(inner)): for ind2 in inner:
ind1 = inner[i] gfrotated[self.SK.solver_to_sumk_block[ish][block]][ind1,ind2] << gf_to_rot[block][ind1,ind2]
ind2 = inner[j]
gfrotated[self.SK.solver_to_sumk_block[s][block]][ind1,ind2] << gf_to_rot[block][ind1,ind2]
# Rotate using the matrix w # Rotate using the matrix w
for bname,gf in gfrotated: for bname,gf in gfrotated:
@ -90,11 +75,9 @@ class TransBasis:
gfreturn = gf_to_rot.copy() gfreturn = gf_to_rot.copy()
# Put back into CTQMC basis: # Put back into CTQMC basis:
for block, inner in self.gf_struct_solver[s].iteritems(): for block, inner in self.gf_struct_solver[ish].iteritems():
for i in range(len(inner)): for ind1 in inner:
for j in range(len(inner)): for ind2 in inner:
ind1 = inner[i]
ind2 = inner[j]
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 return gfreturn
@ -103,55 +86,52 @@ class TransBasis:
def write_trans_file(self, filename): def write_trans_file(self, filename):
'''Writes the new transformation into a file readable by dmftproj.''' '''Writes the new transformation into a file readable by dmftproj.'''
f=open(filename,'w') f = open(filename,'w')
Tnew = self.T.conjugate() Tnew = self.T.conjugate()
N = self.SK.corr_shells[0][3] dim = self.SK.corr_shells[0][3]
if (self.SK.SO==0): if self.SK.SO == 0:
for i in range(N): for i in range(dim):
st = '' st = ''
for k in range(N): for k in range(dim):
st += " %9.6f"%(Tnew[i,k].real) st += " %9.6f"%(Tnew[i,k].real)
st += " %9.6f"%(Tnew[i,k].imag) st += " %9.6f"%(Tnew[i,k].imag)
for k in range(2*N): for k in range(2*dim):
st += " 0.0" st += " 0.0"
if (i<(N-1)): if i < (dim-1):
f.write("%s\n"%(st)) f.write("%s\n"%(st))
else: else:
st1=st.replace(' ','*',1) st1 = st.replace(' ','*',1)
f.write("%s\n"%(st1)) f.write("%s\n"%(st1))
for i in range(dim):
for i in range(N):
st = '' st = ''
for k in range(2*N): for k in range(2*dim):
st += " 0.0" st += " 0.0"
for k in range(N): for k in range(dim):
st += " %9.6f"%(Tnew[i,k].real) st += " %9.6f"%(Tnew[i,k].real)
st += " %9.6f"%(Tnew[i,k].imag) st += " %9.6f"%(Tnew[i,k].imag)
if (i<(N-1)): if i < (dim-1):
f.write("%s\n"%(st)) f.write("%s\n"%(st))
else: else:
st1=st.replace(' ','*',1) st1 = st.replace(' ','*',1)
f.write("%s\n"%(st1)) f.write("%s\n"%(st1))
else: else:
for i in range(N): for i in range(dim):
st = '' st = ''
for k in range(N): for k in range(dim):
st += " %9.6f"%(Tnew[i,k].real) st += " %9.6f"%(Tnew[i,k].real)
st += " %9.6f"%(Tnew[i,k].imag) st += " %9.6f"%(Tnew[i,k].imag)
if (i<(N-1)): if i < (dim-1):
f.write("%s\n"%(st)) f.write("%s\n"%(st))
else: else:
st1=st.replace(' ','*',1) st1 = st.replace(' ','*',1)
f.write("%s\n"%(st1)) f.write("%s\n"%(st1))
#MPI.report("SO not implemented!")
f.close() f.close()