LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2025-01-09 12:44:03 +01:00
Code Issues Releases Wiki Activity

plovasp: fixed some bugs for multiple ions per shell

Browse Source
This commit is contained in:
Malte Schüler 2019-07-02 11:00:21 +02:00
parent d2263ae3d5
commit 2dda711aed
3 changed files with 74 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
python/converters/plovasp/plotools.py
View File

@ -81,6 +81,7 @@ def check_data_consistency(pars, el_struct):
else: else:
errmsg = "Projector for isite = %s, l = %s does not match PROJCAR"%(ion + 1, lshell) errmsg = "Projector for isite = %s, l = %s does not match PROJCAR"%(ion + 1, lshell)
raise Exception(errmsg) raise Exception(errmsg)
################################################################################ ################################################################################
# #
@ -106,7 +107,8 @@ def generate_plo(conf_pars, el_struct):
# eigvals(nktot, nband, ispin) are defined with respect to the Fermi level # eigvals(nktot, nband, ispin) are defined with respect to the Fermi level
eigvals = el_struct.eigvals - efermi eigvals = el_struct.eigvals - efermi
# check if at least one shell is correlated
assert np.any([shell['corr'] for shell in conf_pars.shells]), 'at least one shell has be CORR = True'
nshell = len(conf_pars.shells) nshell = len(conf_pars.shells)
print print
print " Generating %i shell%s..."%(nshell, '' if nshell == 1 else 's') print " Generating %i shell%s..."%(nshell, '' if nshell == 1 else 's')
@ -121,6 +123,7 @@ def generate_plo(conf_pars, el_struct):
print " Correlated : %r"%(pshell.corr) print " Correlated : %r"%(pshell.corr)
pshells.append(pshell) pshells.append(pshell)
pgroups = [] pgroups = []
for gr_par in conf_pars.groups: for gr_par in conf_pars.groups:
pgroup = ProjectorGroup(gr_par, pshells, eigvals) pgroup = ProjectorGroup(gr_par, pshells, eigvals)
@ -436,20 +439,22 @@ def hk_output(conf_pars, el_struct, pgroups):
shell = pgroup.shells[ish] shell = pgroup.shells[ish]
sh_dict = {}
sh_dict['shell_index'] = ish
sh_dict['lorb'] = shell.lorb
sh_dict['ndim'] = shell.ndim
# Convert ion indices from the internal representation (starting from 0)
# to conventional VASP representation (starting from 1)
ion_output = [io + 1 for io in shell.ion_list] ion_output = [io + 1 for io in shell.ion_list]
# Derive sorts from equivalence classes
sh_dict['ion_list'] = ion_output for iion in ion_output:
sh_dict['ion_sort'] = shell.ion_sort sh_dict = {}
sh_dict['shell_index'] = ish
sh_dict['lorb'] = shell.lorb
sh_dict['ndim'] = shell.ndim
# Convert ion indices from the internal representation (starting from 0)
# to conventional VASP representation (starting from 1)
# Derive sorts from equivalence classes
sh_dict['ion_list'] = ion_output
sh_dict['ion_sort'] = shell.ion_sort
head_shells.append(sh_dict) head_shells.append(sh_dict)
with open(hk_fname, 'wt') as f: with open(hk_fname, 'wt') as f:
# Eigenvalues within the window # Eigenvalues within the window

12
python/converters/plovasp/proj_group.py
View File

@ -192,12 +192,18 @@ class ProjectorGroup:
H_ij(k) = sum_l P*_il eps_l P_lj H_ij(k) = sum_l P*_il eps_l P_lj
""" """
# here we abuse the get_block_matrix_map(), however, it only works
# if self.normion is false
temp = self.normion
self.normion = False
block_maps, ndim = self.get_block_matrix_map() block_maps, ndim = self.get_block_matrix_map()
self.normion = temp
_, ns, nk, _, _ = self.shells[0].proj_win.shape _, ns, nk, _, _ = self.shells[0].proj_win.shape
#print(p_mat.shape) #print(p_mat.shape)
print block_maps, ndim
self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128) self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128)
# Note that 'ns' and 'nk' are the same for all shells # Note that 'ns' and 'nk' are the same for all shells
for isp in xrange(ns): for isp in xrange(ns):
@ -218,8 +224,8 @@ class ProjectorGroup:
shell = self.shells[ish] shell = self.shells[ish]
p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb] p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb]
self.hk[isp,ik,:,:] = np.dot(p_mat.conjugate()*eigvals[ik,bmin:bmax,isp], self.hk[isp,ik,:,:] = np.dot(p_mat*eigvals[ik,bmin:bmax,isp],
p_mat.transpose()) p_mat.transpose().conjugate())
################################################################################ ################################################################################

92
python/converters/vasp_converter.py
View File

@ -145,7 +145,8 @@ class VaspConverter(ConverterTools):
# TODO: think about multiple shell groups and how to map them on h5 structures # TODO: think about multiple shell groups and how to map them on h5 structures
assert ng == 1, "Only one group is allowed at the moment" assert ng == 1, "Only one group is allowed at the moment"
try: #try:
if True:
for ig in xrange(ng): for ig in xrange(ng):
gr_file = self.basename + '.pg%i'%(ig + 1) gr_file = self.basename + '.pg%i'%(ig + 1)
jheader, rf = self.read_header_and_data(gr_file) jheader, rf = self.read_header_and_data(gr_file)
@ -165,8 +166,9 @@ class VaspConverter(ConverterTools):
shells = [] shells = []
corr_shells = [] corr_shells = []
shion_to_shell = [[] for ish in xrange(len(p_shells))] shion_to_shell = [[] for ish in xrange(len(p_shells))]
cr_shion_to_shell = [[] for ish in xrange(len(p_shells))]
shorbs_to_globalorbs = [[] for ish in xrange(len(p_shells))] shorbs_to_globalorbs = [[] for ish in xrange(len(p_shells))]
shorbs_to_globalorbs[0] = [0,p_shells[0]['ndim'] ] last_dimension = 0
crshorbs_to_globalorbs = [] crshorbs_to_globalorbs = []
icsh = 0 icsh = 0
for ish, sh in enumerate(p_shells): for ish, sh in enumerate(p_shells):
@ -179,23 +181,22 @@ class VaspConverter(ConverterTools):
pars['l'] = sh['lorb'] pars['l'] = sh['lorb']
#pars['corr'] = sh['corr'] #pars['corr'] = sh['corr']
pars['dim'] = sh['ndim'] pars['dim'] = sh['ndim']
pars['ion_list'] = sh['ion_list'] #pars['ion_list'] = sh['ion_list']
pars['SO'] = SO pars['SO'] = SO
# TODO: check what 'irep' entry does (it seems to be very specific to dmftproj) # TODO: check what 'irep' entry does (it seems to be very specific to dmftproj)
pars['irep'] = 0 pars['irep'] = 0
shells.append(pars) shells.append(pars)
shion_to_shell[ish].append(i) shion_to_shell[ish].append(i)
if ish > 0: shorbs_to_globalorbs[ish].append([last_dimension,
shorbs_to_globalorbs[ish] = [shorbs_to_globalorbs[ish-1][1], last_dimension+sh['ndim']])
shorbs_to_globalorbs[ish-1][1]+sh['ndim']] last_dimension = last_dimension+sh['ndim']
if sh['corr']: if sh['corr']:
corr_shells.append(pars) corr_shells.append(pars)
crshorbs_to_globalorbs.append(shorbs_to_globalorbs[ish])
print shorbs_to_globalorbs
# TODO: generalize this to the case of multiple shell groups # TODO: generalize this to the case of multiple shell groups
n_corr_shells = len(corr_shells) n_corr_shells = len(corr_shells)
#n_shells = n_corr_shells # No non-correlated shells at the moment
#shells = corr_shells
n_orbs = sum([sh['dim'] for sh in shells]) n_orbs = sum([sh['dim'] for sh in shells])
# FIXME: atomic sorts in Wien2K are not the same as in VASP. # FIXME: atomic sorts in Wien2K are not the same as in VASP.
# A symmetry analysis from OUTCAR or symmetry file should be used # A symmetry analysis from OUTCAR or symmetry file should be used
@ -260,18 +261,20 @@ class VaspConverter(ConverterTools):
rf_hk = self.read_data(f_hk) rf_hk = self.read_data(f_hk)
for isp in xrange(n_spin_blocs): for isp in xrange(n_spin_blocs):
for ik in xrange(n_k): for ik in xrange(n_k):
print ik
for ib in xrange(n_orbs): for ib in xrange(n_orbs):
for jb in xrange(n_orbs): for jb in xrange(n_orbs):
hopping[ik, isp, ib, jb] = rf_hk.next() hopping[ik, isp, ib, jb] = rf_hk.next()
for ib in xrange(n_orbs): for ib in xrange(n_orbs):
for jb in xrange(n_orbs): for jb in xrange(n_orbs):
hopping[ik, isp, ib, jb] += 1j*rf_hk.next() hopping[ik, isp, ib, jb] += 1j*rf_hk.next()
rf_hk.close()
# Projectors # Projectors
# print n_orbitals # print n_orbitals
# print [crsh['dim'] for crsh in corr_shells] # print [crsh['dim'] for crsh in corr_shells]
proj_mat_all = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, sum([sh['dim'] for sh in shells]), numpy.max(n_orbitals)], numpy.complex_) proj_mat_csc = numpy.zeros([n_k, n_spin_blocs, sum([sh['dim'] for sh in shells]), numpy.max(n_orbitals)], numpy.complex_)
print 'proj_mat_csc', proj_mat_csc.shape
# TODO: implement reading from more than one projector group # TODO: implement reading from more than one projector group
# In 'dmftproj' each ion represents a separate correlated shell. # In 'dmftproj' each ion represents a separate correlated shell.
# In my interface a 'projected shell' includes sets of ions. # In my interface a 'projected shell' includes sets of ions.
@ -284,52 +287,51 @@ class VaspConverter(ConverterTools):
# #
# At the moment I choose i.2 for its simplicity. But one should consider possible # At the moment I choose i.2 for its simplicity. But one should consider possible
# use cases and decide which solution is to be made permanent. # use cases and decide which solution is to be made permanent.
# #
for ish, sh in enumerate(p_shells): for ish, sh in enumerate(p_shells):
for isp in xrange(n_spin_blocs): for isp in xrange(n_spin_blocs):
for ik in xrange(n_k): for ik in xrange(n_k):
for ion in xrange(len(sh['ion_list'])): for ion in xrange(len(sh['ion_list'])):
icsh = shion_to_shell[ish][ion] for ilm in xrange(shorbs_to_globalorbs[ish][ion][0],shorbs_to_globalorbs[ish][ion][1]):
for ilm in xrange(shorbs_to_globalorbs[ish][0],shorbs_to_globalorbs[ish][1]):
for ib in xrange(n_orbitals[ik, isp]): for ib in xrange(n_orbitals[ik, isp]):
# This is to avoid confusion with the order of arguments # This is to avoid confusion with the order of arguments
pr = rf.next() pr = rf.next()
pi = rf.next() pi = rf.next()
proj_mat_all[ik, isp, icsh, ilm, ib] = complex(pr, pi) proj_mat_csc[ik, isp, ilm, ib] = complex(pr, pi)
# now save only projectors with flag 'corr' to proj_mat # now save only projectors with flag 'corr' to proj_mat
proj_mat = numpy.zeros([n_k, n_spin_blocs, n_corr_shells, sum([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_)
print 'proj_mat', proj_mat.shape
if self.proj_or_hk == 'proj': if self.proj_or_hk == 'proj':
addIndex = 0 for ish, sh in enumerate(p_shells):
for ish, corr_shell in enumerate(corr_shells): if sh['corr']:
for isp in xrange(n_spin_blocs): for isp in xrange(n_spin_blocs):
for ik in xrange(n_k): for ik in xrange(n_k):
for ion in xrange(len(corr_shell['ion_list'])): for ion in xrange(len(sh['ion_list'])):
icsh = shion_to_shell[ish][ion] icsh = shion_to_shell[ish][ion]
for iclm,ilm in enumerate(xrange(shorbs_to_globalorbs[ish][ion][0],shorbs_to_globalorbs[ish][ion][1])):
for iclm,ilm in enumerate(xrange(crshorbs_to_globalorbs[ish][0],crshorbs_to_globalorbs[ish][1])): for ib in xrange(n_orbitals[ik, isp]):
for ib in xrange(n_orbitals[ik, isp]): proj_mat[ik,isp,icsh,iclm,ib] = proj_mat_csc[ik,isp,ilm,ib]
proj_mat[ik,isp,icsh,iclm+addIndex,ib] = proj_mat_all[ik,isp,icsh,ilm,ib]
addIndex += corr_shell['dim']
elif self.proj_or_hk == 'hk': elif self.proj_or_hk == 'hk':
addIndex = 0 print crshorbs_to_globalorbs, proj_mat.shape, shorbs_to_globalorbs
for ish, corr_shell in enumerate(corr_shells): for ish, sh in enumerate(p_shells):
for isp in xrange(n_spin_blocs): if sh['corr']:
for ik in xrange(n_k): for ion in xrange(len(sh['ion_list'])):
for ion in xrange(len(corr_shell['ion_list'])): icsh = shion_to_shell[ish][ion]
icsh = shion_to_shell[ish][ion] for isp in xrange(n_spin_blocs):
for ik in xrange(n_k):
for iclm,ilm in enumerate(xrange(crshorbs_to_globalorbs[ish][0],crshorbs_to_globalorbs[ish][1])): for iclm,ilm in enumerate(xrange(shorbs_to_globalorbs[ish][ion][0],shorbs_to_globalorbs[ish][ion][1])):
proj_mat[ik,isp,icsh,iclm+addIndex,ilm] = 1.0 proj_mat[ik,isp,icsh,iclm,ilm] = 1.0
addIndex += corr_shell['dim'] print proj_mat[0,0,0,:,:].real
corr_shell.pop('ion_list') print proj_mat[0,0,1,:,:].real
#corr_shell.pop('ion_list')
things_to_set = ['n_shells','shells','n_corr_shells','corr_shells','n_spin_blocs','n_orbitals','n_k','SO','SP','energy_unit'] 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: for it in things_to_set:
# print "%s:"%(it), locals()[it] # print "%s:"%(it), locals()[it]
setattr(self,it,locals()[it]) setattr(self,it,locals()[it])
except StopIteration: #except StopIteration:
raise "VaspConverter: error reading %s"%self.gr_file # raise "VaspConverter: error reading %s"%self.gr_file
rf.close() rf.close()
@ -340,10 +342,12 @@ class VaspConverter(ConverterTools):
# 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', 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', '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', 'rot_mat_time_inv','n_reps','dim_reps','T','n_orbitals','proj_mat','bz_weights','hopping',
# 'proj_mat_all',
'bz_weights','hopping',
'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr'] 'n_inequiv_shells', 'corr_to_inequiv', 'inequiv_to_corr']
if self.proj_or_hk == 'hk':
proj_or_hk = self.proj_or_hk
things_to_save.append('proj_mat_csc')
things_to_save.append('proj_or_hk')
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]
# Store Fermi weights to 'dft_misc_input' # Store Fermi weights to 'dft_misc_input'