LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-07-17 08:30:35 +02:00
Code Issues Releases Wiki Activity

SumK generated output in new group 'lda_output', cleaning up.

Browse Source 
* updated update_archive.py script to move output elements out
* added clear_lda_output.py script to restore state of converter generated h5
This commit is contained in:
Priyanka Seth 2014-11-07 00:55:40 +01:00
parent 8bbbe81c7d
commit 8d44e5ce91
5 changed files with 112 additions and 92 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
python/TODOFIX
View File

@ -8,6 +8,7 @@ Changed the following:
* retval -> read_value * retval -> read_value
* Gupf -> G_upfold * Gupf -> G_upfold
* read_symmetry_input -> convert_symmetry_input * read_symmetry_input -> convert_symmetry_input
* Symm_corr -> symmcorr
********** **********
* changed default h5 subgroup names * changed default h5 subgroup names

24
python/clear_lda_output.py Normal file
View File

@ -0,0 +1,24 @@
import h5py
import sys
import subprocess
if len(sys.argv) < 2:
print "Usage: python clear_lda_output.py archive"
sys.exit()
print """
This script is to remove any SumkLDA generated output from the h5 archive
and to restore it to the original post-converter state.
"""
filename = sys.argv[1]
A = h5py.File(filename)
del(A["lda_output"])
A.close()
# Repack to reclaim disk space
retcode = subprocess.call(["h5repack","-i%s"%filename, "-otemphgfrt.h5"])
if retcode != 0:
print "h5repack failed!"
else:
subprocess.call(["mv","-f","temphgfrt.h5","%s"%filename])

147
python/sumk_lda.py
View File

@ -32,13 +32,14 @@ class SumkLDA:
"""This class provides a general SumK method for combining ab-initio code and pytriqs.""" """This class provides a general SumK method for combining ab-initio code and pytriqs."""
def __init__(self, hdf_file, mu = 0.0, h_field = 0.0, use_lda_blocks = False, lda_data = 'lda_input', symmcorr_data = 'lda_symmcorr_input', def __init__(self, hdf_file, mu = 0.0, h_field = 0.0, use_lda_blocks = False,
parproj_data = 'lda_parproj_input', symmpar_data = 'lda_symmpar_input', bands_data = 'lda_bands_input'): lda_data = 'lda_input', symmcorr_data = 'lda_symmcorr_input', parproj_data = 'lda_parproj_input',
symmpar_data = 'lda_symmpar_input', bands_data = 'lda_bands_input', lda_output = 'lda_output'):
""" """
Initialises the class from data previously stored into an HDF5 Initialises the class from data previously stored into an HDF5
""" """
if not (type(hdf_file)==StringType): if not (type(hdf_file)==StringType):
mpi.report("Give a string for the HDF5 filename to read the input!") mpi.report("Give a string for the HDF5 filename to read the input!")
else: else:
self.hdf_file = hdf_file self.hdf_file = hdf_file
@ -47,6 +48,7 @@ class SumkLDA:
self.parproj_data = parproj_data self.parproj_data = parproj_data
self.symmpar_data = symmpar_data self.symmpar_data = symmpar_data
self.bands_data = bands_data self.bands_data = bands_data
self.lda_output = lda_output
self.block_names = [ ['up','down'], ['ud'] ] self.block_names = [ ['up','down'], ['ud'] ]
self.n_spin_blocks_gf = [2,1] self.n_spin_blocks_gf = [2,1]
self.G_upfold = None self.G_upfold = None
@ -56,15 +58,12 @@ class SumkLDA:
things_to_read = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required', things_to_read = ['energy_unit','n_k','k_dep_projection','SP','SO','charge_below','density_required',
'symm_op','n_shells','shells','n_corr_shells','corr_shells','use_rotations','rot_mat', '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'] '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'] self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.lda_data, things_to_read = things_to_read)
self.read_value = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read,optional_things=optional_things)
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!")
# determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells) # determine the number of inequivalent correlated shells (self.n_inequiv_corr_shells)
# and related maps (self.shellmap, self.invshellmap) # and related maps (self.shellmap, self.invshellmap)
self.inequiv_shells(self.corr_shells) self.inequiv_shells(self.corr_shells)
@ -73,14 +72,19 @@ class SumkLDA:
self.names_to_ind = [{}, {}] self.names_to_ind = [{}, {}]
for ibl in range(2): for ibl in range(2):
for inm in range(self.n_spin_blocks_gf[ibl]): for inm in range(self.n_spin_blocks_gf[ibl]):
self.names_to_ind[ibl][self.block_names[ibl][inm]] = inm * self.SP #(self.Nspinblocs-1) self.names_to_ind[ibl][self.block_names[ibl][inm]] = inm * self.SP
# GF structure used for the local things in the k sums # GF structure used for the local things in the k sums
# Most general form allowing for all hybridisation, i.e. largest blocks possible # Most general form allowing for all hybridisation, i.e. largest blocks possible
self.gf_struct_corr = [ [ (al, range( self.corr_shells[i][3])) for al in self.block_names[self.corr_shells[i][4]] ] self.gf_struct_corr = [ [ (al, range( self.corr_shells[i][3])) for al in self.block_names[self.corr_shells[i][4]] ]
for i in xrange(self.n_corr_shells) ] for i in xrange(self.n_corr_shells) ]
if not (self.read_value['gf_struct_solver']): #-----
# If these quantities are not in HDF, set them up and save into HDF
optional_things = ['gf_struct_solver','map_inv','map','chemical_potential','dc_imp','dc_energ','deg_shells']
self.subgroup_present, self.value_read = self.read_input_from_hdf(subgrp = self.lda_output, things_to_read = [],
optional_things = optional_things)
if (not self.subgroup_present) or (not self.value_read['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 = [ dict([ (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]] ])
@ -93,100 +97,91 @@ class SumkLDA:
self.map[i][al] = [al for j in range( self.corr_shells[self.invshellmap[i]][3] ) ] self.map[i][al] = [al for j in range( self.corr_shells[self.invshellmap[i]][3] ) ]
self.map_inv[i][al] = al self.map_inv[i][al] = al
if not (self.read_value['dc_imp']): if (not self.subgroup_present) or (not self.value_read['dc_imp']):
# init the double counting: # init the double counting:
self.__init_dc() self.__init_dc()
if not (self.read_value['chemical_potential']): if (not self.subgroup_present) or (not self.value_read['chemical_potential']):
self.chemical_potential = mu self.chemical_potential = mu
if not (self.read_value['deg_shells']): if (not self.subgroup_present) or (not self.value_read['deg_shells']):
self.deg_shells = [ [] for i in range(self.n_inequiv_corr_shells)] self.deg_shells = [ [] for i in range(self.n_inequiv_corr_shells)]
#-----
if self.symm_op: if self.symm_op:
#mpi.report("Do the init for symm:") #mpi.report("Do the init for symm:")
self.Symm_corr = Symmetry(hdf_file,subgroup=self.symmcorr_data) self.symmcorr = Symmetry(hdf_file,subgroup=self.symmcorr_data)
# Analyse the block structure and determine the smallest blocs, if desired # 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()
# Now save things again to HDF5:
# now save things again to HDF5: # FIXME WHAT HAPPENS TO h_field? INPUT TO __INIT__? ADD TO OPTIONAL_THINGS?
if (mpi.is_master_node()): things_to_save=['chemical_potential','dc_imp','dc_energ','h_field']
ar=HDFArchive(self.hdf_file,'a') self.save(things_to_save)
ar[self.lda_data]['h_field'] = self.h_field
del ar
self.save()
def read_input_from_hdf(self, subgrp, things_to_read=[], optional_things=[]):
def read_input_from_hdf(self, subgrp, things_to_read, optional_things=[]):
""" """
Reads data from the HDF file Reads data from the HDF file
""" """
read_value = True value_read = True
# init variables on all nodes to ensure mpi broadcast works at the end # init variables on all nodes to ensure mpi broadcast works at the end
for it in things_to_read: setattr(self,it,0) for it in things_to_read: setattr(self,it,0)
for it in optional_things: setattr(self,it,0) for it in optional_things: setattr(self,it,0)
if (mpi.is_master_node()): if mpi.is_master_node():
ar=HDFArchive(self.hdf_file,'a') ar=HDFArchive(self.hdf_file,'a')
if (subgrp in ar): if subgrp in ar:
subgroup_present = True
# first read the necessary things: # first read the necessary things:
for it in things_to_read: for it in things_to_read:
if (it in ar[subgrp]): if it in ar[subgrp]:
setattr(self,it,ar[subgrp][it]) setattr(self,it,ar[subgrp][it])
else: else:
mpi.report("Loading %s failed!"%it) mpi.report("Loading %s failed!"%it)
read_value = False value_read = False
if ((read_value) and (len(optional_things)>0)): if (value_read and (len(optional_things)>0)):
# if necessary things worked, now read optional things: # if necessary things worked, now read optional things:
read_value = {} value_read = {}
for it in optional_things: for it in optional_things:
if (it in ar[subgrp]): if it in ar[subgrp]:
setattr(self,it,ar[subgrp][it]) setattr(self,it,ar[subgrp][it])
read_value['%s'%it] = True value_read['%s'%it] = True
else: else:
read_value['%s'%it] = False value_read['%s'%it] = False
else: else:
mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp) mpi.report("Loading failed: No %s subgroup in HDF5!"%subgrp)
read_value = False subgroup_present = False
value_read = False
del ar del ar
# now do the broadcasting: # now do the broadcasting:
for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it))) for it in things_to_read: setattr(self,it,mpi.bcast(getattr(self,it)))
for it in optional_things: setattr(self,it,mpi.bcast(getattr(self,it))) for it in optional_things: setattr(self,it,mpi.bcast(getattr(self,it)))
subgroup_present = mpi.bcast(subgroup_present)
value_read = mpi.bcast(value_read)
read_value = mpi.bcast(read_value) return (subgroup_present, value_read)
return read_value
def save(self,things_to_save):
def save(self):
"""Saves some quantities into an HDF5 arxiv""" """Saves some quantities into an HDF5 arxiv"""
if not (mpi.is_master_node()): return # do nothing on nodes if not (mpi.is_master_node()): return # do nothing on nodes
ar=HDFArchive(self.hdf_file,'a') ar = HDFArchive(self.hdf_file,'a')
things_to_save=['chemical_potential','dc_imp','dc_energ'] if not (self.lda_output in ar): ar.create_group(self.lda_output)
for it in things_to_save: ar[self.lda_data][it] = getattr(self,it) for it in things_to_save:
try:
ar[self.lda_output][it] = getattr(self,it)
except:
mpi.report("%s not found, and so not stored."%it)
del ar del ar
def load(self):
"""Loads some quantities from an HDF5 arxiv"""
things_to_read=['chemical_potential','dc_imp','dc_energ']
read_value = self.read_input_from_hdf(subgrp=self.lda_data,things_to_read=things_to_read)
return read_value
def downfold(self,ik,icrsh,sig,gf_to_downfold,gf_inp): def downfold(self,ik,icrsh,sig,gf_to_downfold,gf_inp):
"""Downfolding a block of the Greens function""" """Downfolding a block of the Greens function"""
@ -344,7 +339,7 @@ class SumkLDA:
mpi.barrier() mpi.barrier()
if (self.symm_op!=0): dens_mat = self.Symm_corr.symmetrize(dens_mat) if (self.symm_op!=0): dens_mat = self.symmcorr.symmetrize(dens_mat)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if (self.use_rotations):
@ -391,7 +386,7 @@ class SumkLDA:
mpi.barrier() mpi.barrier()
if (self.symm_op!=0): dens_mat = self.Symm_corr.symmetrize(dens_mat) if (self.symm_op!=0): dens_mat = self.symmcorr.symmetrize(dens_mat)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if (self.use_rotations):
@ -486,16 +481,8 @@ class SumkLDA:
elif ((ind1<0)and(ind2<0)): elif ((ind1<0)and(ind2<0)):
self.deg_shells[ish].append([bl[0],bl2[0]]) self.deg_shells[ish].append([bl[0],bl2[0]])
if (mpi.is_master_node()): things_to_save=['gf_struct_solver','map','map_inv','deg_shells']
ar=HDFArchive(self.hdf_file,'a') self.save(things_to_save)
ar[self.lda_data]['gf_struct_solver'] = self.gf_struct_solver
ar[self.lda_data]['map'] = self.map
ar[self.lda_data]['map_inv'] = self.map_inv
try:
ar[self.lda_data]['deg_shells'] = self.deg_shells
except:
mpi.report("deg_shells not stored, degeneracies not found")
del ar
return dens_mat return dens_mat
@ -552,7 +539,7 @@ class SumkLDA:
self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose() ) self.proj_mat[ik,isp,icrsh,0:dim,0:n_orb].conjugate().transpose() )
# symmetrisation: # symmetrisation:
if (self.symm_op!=0): self.Hsumk = self.Symm_corr.symmetrize(self.Hsumk) if (self.symm_op!=0): self.Hsumk = self.symmcorr.symmetrize(self.Hsumk)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if (self.use_rotations):
@ -743,9 +730,8 @@ class SumkLDA:
def set_mu(self,mu): def set_mu(self,mu):
"""Sets a new chemical potential""" """Sets a new chemical potential"""
self.chemical_potential = mu
#print "Chemical potential in SumK set to ",mu
self.chemical_potential = mu
def inequiv_shells(self,lst): def inequiv_shells(self,lst):
@ -810,13 +796,13 @@ class SumkLDA:
F = lambda mu : self.total_density(mu = mu) F = lambda mu : self.total_density(mu = mu)
Dens_rel = self.density_required - self.charge_below density = self.density_required - self.charge_below
self.chemical_potential = dichotomy.dichotomy(function = F, self.chemical_potential = dichotomy.dichotomy(function = F,
x_init = self.chemical_potential, y_value = Dens_rel, x_init = self.chemical_potential, y_value = density,
precision_on_y = precision, delta_x=0.5, precision_on_y = precision, delta_x = 0.5, max_loops = 100,
max_loops = 100, x_name="chemical_potential", y_name= "Total Density", x_name = "Chemical Potential", y_name = "Total Density",
verbosity = 3)[0] verbosity = 3)[0]
return self.chemical_potential return self.chemical_potential
@ -856,7 +842,7 @@ class SumkLDA:
# Gloc[:] is now the sum over k projected to the local orbitals. # Gloc[:] is now the sum over k projected to the local orbitals.
# here comes the symmetrisation, if needed: # here comes the symmetrisation, if needed:
if (self.symm_op!=0): Gloc = self.Symm_corr.symmetrize(Gloc) if (self.symm_op!=0): Gloc = self.symmcorr.symmetrize(Gloc)
# Gloc is rotated to the local coordinate system: # Gloc is rotated to the local coordinate system:
if (self.use_rotations): if (self.use_rotations):
@ -1003,7 +989,6 @@ class SumkLDA:
def find_mu_nonint(self, dens_req, orb = None, beta = 40, precision = 0.01): def find_mu_nonint(self, dens_req, orb = None, beta = 40, precision = 0.01):
def F(mu): def F(mu):
#gnonint = self.nonint_G(beta=beta,mu=mu)
gnonint = self.extract_G_loc(mu=mu,with_Sigma=False) gnonint = self.extract_G_loc(mu=mu,with_Sigma=False)
if (orb is None): if (orb is None):
@ -1017,10 +1002,10 @@ class SumkLDA:
self.chemical_potential = dichotomy.dichotomy(function = F, self.chemical_potential = dichotomy.dichotomy(function = F,
x_init = self.chemical_potential, y_value = dens_req, x_init = self.chemical_potential, y_value = dens_req,
precision_on_y = precision, delta_x=0.5, precision_on_y = precision, delta_x = 0.5, max_loops = 100,
max_loops = 100, x_name="chemical_potential", y_name= "Local Density", x_name = "Chemical Potential", y_name = "Total Density",
verbosity = 3)[0] verbosity = 3)[0]
return self.chemical_potential return self.chemical_potential
@ -1042,12 +1027,12 @@ class SumkLDA:
if (dens_req is None): if (dens_req is None):
Dens_rel = self.density_required - self.charge_below density = self.density_required - self.charge_below
else: else:
Dens_rel = dens_req density = dens_req
dcnew = dichotomy.dichotomy(function = F, dcnew = dichotomy.dichotomy(function = F,
x_init = guess, y_value = Dens_rel, x_init = guess, y_value = density,
precision_on_y = precision, delta_x=0.5, precision_on_y = precision, delta_x=0.5,
max_loops = 100, x_name="Double-Counting", y_name= "Total Density", max_loops = 100, x_name="Double-Counting", y_name= "Total Density",
verbosity = 3)[0] verbosity = 3)[0]
@ -1068,7 +1053,7 @@ class SumkLDA:
n_orb = self.n_orbitals[ik,0] n_orb = self.n_orbitals[ik,0]
dens_mat[ish][:,:] += numpy.dot(self.proj_mat[ik,0,ish,0:dim,0:n_orb],self.proj_mat[ik,0,ish,0:dim,0:n_orb].transpose().conjugate()) * self.bz_weights[ik] dens_mat[ish][:,:] += numpy.dot(self.proj_mat[ik,0,ish,0:dim,0:n_orb],self.proj_mat[ik,0,ish,0:dim,0:n_orb].transpose().conjugate()) * self.bz_weights[ik]
if (self.symm_op!=0): dens_mat = self.Symm_corr.symmetrize(dens_mat) if (self.symm_op!=0): dens_mat = self.symmcorr.symmetrize(dens_mat)
# Rotate to local coordinate system: # Rotate to local coordinate system:
if (self.use_rotations): if (self.use_rotations):

22
python/sumk_lda_tools.py
View File

@ -202,7 +202,7 @@ class SumkLDATools(SumkLDA):
if (self.symm_op!=0): Gloc = self.Symm_corr.symmetrize(Gloc) if (self.symm_op!=0): Gloc = self.symmcorr.symmetrize(Gloc)
if (self.use_rotations): if (self.use_rotations):
for icrsh in xrange(self.n_corr_shells): for icrsh in xrange(self.n_corr_shells):
@ -243,8 +243,8 @@ class SumkLDATools(SumkLDA):
""" """
things_to_read = ['dens_mat_below','n_parproj','proj_mat_pc','rot_mat_all','rot_mat_all_time_inv'] things_to_read = ['dens_mat_below','n_parproj','proj_mat_pc','rot_mat_all','rot_mat_all_time_inv']
read_value = self.read_input_from_hdf(subgrp=self.parproj_data,things_to_read = things_to_read) value_read = self.read_input_from_hdf(subgrp=self.parproj_data,things_to_read = things_to_read)
return read_value return value_read
@ -254,9 +254,9 @@ class SumkLDATools(SumkLDA):
assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" assert hasattr(self,"Sigma_imp"), "Set Sigma First!!"
#things_to_read = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all'] #things_to_read = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all']
#read_value = self.read_input_from_HDF(SubGrp=self.parproj_data, things_to_read=things_to_read) #value_read = self.read_input_from_HDF(SubGrp=self.parproj_data, things_to_read=things_to_read)
read_value = self.read_parproj_input_from_hdf() value_read = self.read_parproj_input_from_hdf()
if not read_value: return read_value if not value_read: return value_read
if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symmpar_data) if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symmpar_data)
mu = self.chemical_potential mu = self.chemical_potential
@ -348,8 +348,8 @@ class SumkLDATools(SumkLDA):
assert hasattr(self,"Sigma_imp"), "Set Sigma First!!" assert hasattr(self,"Sigma_imp"), "Set Sigma First!!"
things_to_read = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc'] things_to_read = ['n_k','n_orbitals','proj_mat','hopping','n_parproj','proj_mat_pc']
read_value = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=things_to_read) value_read = self.read_input_from_hdf(subgrp=self.bands_data,things_to_read=things_to_read)
if not read_value: return read_value if not value_read: return value_read
if fermi_surface: ishell=None if fermi_surface: ishell=None
@ -515,9 +515,9 @@ class SumkLDATools(SumkLDA):
#things_to_read = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all'] #things_to_read = ['Dens_Mat_below','N_parproj','Proj_Mat_pc','rotmat_all']
#read_value = self.read_input_from_HDF(SubGrp=self.parproj_data,things_to_read=things_to_read) #value_read = self.read_input_from_HDF(SubGrp=self.parproj_data,things_to_read=things_to_read)
read_value = self.read_parproj_input_from_hdf() value_read = self.read_parproj_input_from_hdf()
if not read_value: return read_value if not value_read: return value_read
if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symmpar_data) if self.symm_op: self.Symm_par = Symmetry(self.hdf_file,subgroup=self.symmpar_data)
# Density matrix in the window # Density matrix in the window

10
python/update_archive.py
View File

@ -25,6 +25,16 @@ for old, new in old_to_new.iteritems():
print "Changing %s to %s ..."%(old, new) print "Changing %s to %s ..."%(old, new)
A.copy(old,new) A.copy(old,new)
del(A[old]) del(A[old])
move_to_output = ['gf_struct_solver','map_inv','map',
'chemical_potential','dc_imp','dc_energ','deg_shells',
'h_field']
for obj in move_to_output:
if (obj in A['lda_input'].keys()) and (not 'lda_output' in A): A.create_group('lda_output')
print "Moving %s to lda_output ..."%obj
A.copy('lda_input/'+obj,'lda_output/'+obj)
del(A['lda_input'][obj])
A.close() A.close()
# Repack to reclaim disk space # Repack to reclaim disk space