From 96a6891f64c82dd1fab5e0c0880d373071939826 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Malte=20Sch=C3=BCler?= Date: Thu, 27 Jun 2019 15:21:07 +0200 Subject: [PATCH] plovasp: first working version of H(k) and complement --- python/converters/plovasp/plotools.py | 160 ++++++++++++++++++++---- python/converters/plovasp/proj_group.py | 108 +++++++++++++++- python/converters/plovasp/proj_shell.py | 81 ++++++++++++ 3 files changed, 322 insertions(+), 27 deletions(-) diff --git a/python/converters/plovasp/plotools.py b/python/converters/plovasp/plotools.py index e980858d..000c2123 100644 --- a/python/converters/plovasp/plotools.py +++ b/python/converters/plovasp/plotools.py @@ -38,6 +38,7 @@ import itertools as it import numpy as np from proj_group import ProjectorGroup from proj_shell import ProjectorShell +from proj_shell import ComplementShell np.set_printoptions(suppress=True) @@ -124,24 +125,29 @@ def generate_plo(conf_pars, el_struct): for gr_par in conf_pars.groups: pgroup = ProjectorGroup(gr_par, pshells, eigvals) pgroup.orthogonalize() + if gr_par['complement']: + pgroup.complement(eigvals) + if conf_pars.general['hk']: + pgroup.calc_hk(eigvals) # DEBUG output print "Density matrix:" nimp = 0.0 ov_all = [] for ish in pgroup.ishells: - print " Shell %i"%(ish + 1) - dm_all, ov_all_ = pshells[ish].density_matrix(el_struct) - ov_all.append(ov_all_[0]) - spin_fac = 2 if dm_all.shape[0] == 1 else 1 - for io in xrange(dm_all.shape[1]): - print " Site %i"%(io + 1) - dm = spin_fac * dm_all[:, io, : ,:].sum(0) - for row in dm: - print ''.join(map("{0:14.7f}".format, row)) - ndm = dm.trace() - if pshells[ish].corr: - nimp += ndm - print " trace: ", ndm + if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): + print " Shell %i"%(ish + 1) + dm_all, ov_all_ = pshells[ish].density_matrix(el_struct) + ov_all.append(ov_all_[0]) + spin_fac = 2 if dm_all.shape[0] == 1 else 1 + for io in xrange(dm_all.shape[1]): + print " Site %i"%(io + 1) + dm = spin_fac * dm_all[:, io, : ,:].sum(0) + for row in dm: + print ''.join(map("{0:14.7f}".format, row)) + ndm = dm.trace() + if pshells[ish].corr: + nimp += ndm + print " trace: ", ndm print print " Impurity density:", nimp print @@ -152,12 +158,13 @@ def generate_plo(conf_pars, el_struct): print print "Local Hamiltonian:" for ish in pgroup.ishells: - print " Shell %i"%(ish + 1) - loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct) - for io in xrange(loc_ham.shape[1]): - print " Site %i"%(io + 1) - for row in loc_ham[:, io, :, :].sum(0): - print ''.join(map("{0:14.7f}".format, row)) + if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): + print " Shell %i"%(ish + 1) + loc_ham = pshells[pgroup.ishells[ish]].local_hamiltonian(el_struct) + for io in xrange(loc_ham.shape[1]): + print " Site %i"%(io + 1) + for row in loc_ham[:, io, :, :].sum(0): + print ''.join(map("{0:14.7f}".format, row)) # END DEBUG output if 'dosmesh' in conf_pars.general: print @@ -173,13 +180,14 @@ def generate_plo(conf_pars, el_struct): emesh = np.linspace(dos_emin, dos_emax, n_points) for ish in pgroup.ishells: - print " Shell %i"%(ish + 1) - dos = pshells[pgroup.ishells[ish]].density_of_states(el_struct, emesh) - de = emesh[1] - emesh[0] - ntot = (dos[1:,...] + dos[:-1,...]).sum(0) / 2 * de - print " Total number of states:", ntot - for io in xrange(dos.shape[2]): - np.savetxt('pdos_%i_%i.dat'%(ish,io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T) + if not isinstance(pshells[pgroup.ishells[ish]],ComplementShell): + print " Shell %i"%(ish + 1) + dos = pshells[pgroup.ishells[ish]].density_of_states(el_struct, emesh) + de = emesh[1] - emesh[0] + ntot = (dos[1:,...] + dos[:-1,...]).sum(0) / 2 * de + print " Total number of states:", ntot + for io in xrange(dos.shape[2]): + np.savetxt('pdos_%i_%i.dat'%(ish,io), np.vstack((emesh.T, dos[:, 0, io, :].T)).T) pgroups.append(pgroup) @@ -196,6 +204,8 @@ def output_as_text(pars, el_struct, pshells, pgroups): """ ctrl_output(pars, el_struct, len(pgroups)) plo_output(pars, el_struct, pshells, pgroups) + if pars.general['hk']: + hk_output(pars, el_struct, pshells, pgroups) # TODO: k-points with weights should be stored once and for all @@ -388,4 +398,102 @@ def plo_output(conf_pars, el_struct, pshells, pgroups): f.write("{0:16.10f}{1:16.10f}\n".format(p.real, p.imag)) f.write("\n") +################################################################################ +# +# plo_output +# +################################################################################ +def hk_output(conf_pars, el_struct, pshells, pgroups): + """ + Outputs HK into text file. + Filename is defined by that is passed from config-file. + + The Hk for all groups is stored in a '.hk' file. The format is as + defined in the Hk dft_tools format + + # Energy window: emin, emax + ib_min, ib_max + nelect + # Eigenvalues + isp, ik1, kx, ky, kz, kweight + ib1, ib2 + eig1 + eig2 + ... + eigN + ik2, kx, ky, kz, kweight + ... + + # Projected shells + Nshells + # Shells: + # Shell <1> + Shell 1 + ndim + # complex arrays: plo(ns, nion, ndim, nb) + ... + # Shells: + # Shell <2> + Shell 2 + ... + + """ + + print 'writing hk' + for ig, pgroup in enumerate(pgroups): + hk_fname = conf_pars.general['basename'] + '.hk%i'%(ig + 1) + + + head_corr_shells = [] + head_shells = [] + for ish in pgroup.ishells: + + 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] +# Derive sorts from equivalence classes + sh_dict['ion_list'] = ion_output + sh_dict['ion_sort'] = shell.ion_sort + +# TODO: add the output of transformation matrices + + head_shells.append(sh_dict) + if shell.corr: + head_corr_shells.append(sh_dict) + + #head_dict['shells'] = head_shells + + #header = json.dumps(head_dict, indent=4, separators=(',', ': ')) + + with open(hk_fname, 'wt') as f: + +# Eigenvalues within the window + nk, nband, ns_band = el_struct.eigvals.shape + f.write('%i # number of kpoints\n'%nk) + f.write('%f # electron density\n'%1.0) + f.write('%i # number of shells\n'%len(head_shells)) + for head in head_shells: + f.write('%i %i %i %i # atom sort l dim\n'%(head['ion_list'][0],head['ion_sort'][0],head['lorb'],head['ndim'])) + f.write('%i # number of correlated shells\n'%len(head_corr_shells)) + for head in head_corr_shells: + f.write('%i %i %i %i 0 0 # atom sort l dim SO irrep\n'%(head['ion_list'][0],head['ion_sort'][0],head['lorb'],head['ndim'])) + f.write('1 5 # number of ireps, dim of irep\n') + norbs = pgroup.hk.shape[2] + for isp in xrange(ns_band): + #f.write("# is = %i\n"%(isp + 1)) + for ik in xrange(nk): + for io in xrange(norbs): + for iop in xrange(norbs): + f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].real)) + f.write("\n") + for io in xrange(norbs): + for iop in xrange(norbs): + f.write(" {0:14.10f}".format(pgroup.hk[isp,ik,io,iop].imag)) + f.write("\n") diff --git a/python/converters/plovasp/proj_group.py b/python/converters/plovasp/proj_group.py index 1ad99dd6..1fff7a04 100644 --- a/python/converters/plovasp/proj_group.py +++ b/python/converters/plovasp/proj_group.py @@ -30,7 +30,7 @@ r""" Storage and manipulation of projector groups. """ import numpy as np - +from proj_shell import ComplementShell np.set_printoptions(suppress=True) ################################################################################ @@ -165,7 +165,113 @@ class ProjectorGroup: nlm = i2 - i1 + 1 shell = self.shells[ish] shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb] + + +################################################################################ +# +# calc_hk +# +################################################################################ + def calc_hk(self, eigvals): + """ + Calculate H(k) for a group by applying the projectors to the eigenvalues. + + """ + block_maps, ndim = self.get_block_matrix_map() + + _, ns, nk, _, _ = self.shells[0].proj_win.shape + p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128) + #print(p_mat.shape) + self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128) +# Note that 'ns' and 'nk' are the same for all shells + for isp in xrange(ns): + for ik in xrange(nk): + bmin = self.ib_win[ik, isp, 0] + bmax = self.ib_win[ik, isp, 1]+1 + + nb = bmax - bmin + #print(bmin,bmax,nb) +# Combine all projectors of the group to one block projector + for bl_map in block_maps: + p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point and block! + for ibl, block in enumerate(bl_map): + i1, i2 = block['bmat_range'] + ish, ion = block['shell_ion'] + nlm = i2 - i1 + 1 + shell = self.shells[ish] + p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb] + #print(p_mat.shape, eigvals[ik,:,isp].shape) + self.hk[isp,ik,:,:] = np.einsum('ij,jk->ik',p_mat*eigvals[ik,bmin:bmax,isp], + p_mat.transpose().conjugate()) + #print(np.linalg.eigvals(self.hk[isp,ik,...])) + #print(eigvals[ik,bmin:bmax,isp]) + #print(self.hk.shape) + #self.hk[isp,ik,:,:] = np.dot(p_mat.conjugate()*eigvals[ik,bmin:bmax,isp], + # p_mat.transpose()) + + +################################################################################ +# +# complement +# +################################################################################ + def complement(self,eigvals): + """ + Calculate the complement for a group of projectors. + + """ + + print 'claculating complemet' + + block_maps, ndim = self.get_block_matrix_map() + _, ns, nk, _, _ = self.shells[0].proj_win.shape + p_mat = np.zeros((ndim, self.nb_max), dtype=np.complex128) + p_full = np.zeros((1,ns,nk,self.nb_max, self.nb_max), dtype=np.complex128) + + #print(p_mat.shape) + self.hk = np.zeros((ns,nk,ndim,ndim), dtype=np.complex128) +# Note that 'ns' and 'nk' are the same for all shells + orbs_done = 1*ndim + while orbs_done < self.nb_max: + + for isp in xrange(ns): + for ik in xrange(nk): + bmin = self.ib_win[ik, isp, 0] + bmax = self.ib_win[ik, isp, 1]+1 + + nb = bmax - bmin + #print(bmin,bmax,nb) + # Combine all projectors of the group to one block projector + for bl_map in block_maps: + p_mat[:, :] = 0.0j # !!! Clean-up from the last k-point and block! + for ibl, block in enumerate(bl_map): + i1, i2 = block['bmat_range'] + ish, ion = block['shell_ion'] + nlm = i2 - i1 + 1 + shell = self.shells[ish] + p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb] + p_full[0,isp,ik,:ndim,:] = p_mat + proj_work = p_full[0,isp,ik,:,:] + overlap = np.dot(proj_work.transpose().conjugate(),proj_work) + work = np.eye(self.nb_max) - overlap + norm = np.sqrt(np.sum(work*work.transpose().conjugate(),axis=1)) + max_ind = np.argmax(norm) + p_full[0,isp,ik,orbs_done,:] = work[max_ind]/norm[max_ind] + + orbs_done += 1 + sh_pars = {} + sh_pars['lshell'] = -1 + sh_pars['ions'] = {'nion':1,'ion_list':[[1]]} + sh_pars['user_index'] = 'complement' + sh_pars['corr'] = False + self.shells.append(ComplementShell(sh_pars,p_full[:,:,:,ndim:,:],False)) + self.ishells.append(self.ishells[-1]+1) + print(self.shells) + print(self.ishells) + #p_full -= p_mat + + ################################################################################ # # gen_block_matrix_map diff --git a/python/converters/plovasp/proj_shell.py b/python/converters/plovasp/proj_shell.py index 4aa6d8ed..463544a0 100644 --- a/python/converters/plovasp/proj_shell.py +++ b/python/converters/plovasp/proj_shell.py @@ -450,5 +450,86 @@ class ProjectorShell: return dos +################################################################################ +################################################################################ +# +# class ProjectorShell +# +################################################################################ +################################################################################ +class ComplementShell(ProjectorShell): + """ + Container of projectors related to a complement shell. + Parameters: + + - sh_pars (dict) : shell parameters from the config-file + - proj_compl (numpy.array) : array of complement projectors + + """ + def __init__(self, sh_pars, proj_compl, nc_flag): + self.lorb = sh_pars['lshell'] + self.ions = sh_pars['ions'] + self.user_index = sh_pars['user_index'] + self.corr = sh_pars['corr'] + self.nc_flag = nc_flag + + + #self.lm1 = self.lorb**2 + #self.lm2 = (self.lorb+1)**2 + + self.nion = self.ions['nion'] +# Extract ion list and equivalence classes (ion sorts) + self.ion_list = sorted(it.chain(*self.ions['ion_list'])) + self.ion_sort = [] + for ion in self.ion_list: + for icl, eq_cl in enumerate(self.ions['ion_list']): + if ion in eq_cl: + self.ion_sort.append(icl + 1) # Enumerate classes starting from 1 + break + + self.ndim = proj_compl.shape[3] + self.proj_win = proj_compl + + def extract_tmatrices(self, sh_pars): + pass + + def local_hamiltonian(self, el_struct, site_diag=True, spin_diag=True): + """ + Returns occupation matrix/matrices for the shell. + """ + nion, ns, nk, nlm, nbtot = self.proj_win.shape + + assert site_diag, "site_diag = False is not implemented" + assert spin_diag, "spin_diag = False is not implemented" + + loc_ham = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) + + return loc_ham + + def density_matrix(self, el_struct, site_diag=True, spin_diag=True): + """ + Returns occupation matrix/matrices for the shell. + """ + nion, ns, nk, nlm, nbtot = self.proj_win.shape + +# assert site_diag, "site_diag = False is not implemented" + assert spin_diag, "spin_diag = False is not implemented" + + if site_diag: + occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) + overlaps = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) + else: + ndim = nion * nlm + occ_mats = np.zeros((ns, 1, ndim, ndim), dtype=np.float64) + overlaps = np.zeros((ns, 1, ndim, ndim), dtype=np.float64) + + return occ_mats, overlaps + + def density_of_states(self, el_struct, emesh): + nion, ns, nk, nlm, nbtot = self.proj_win.shape + ne = len(emesh) + dos = np.zeros((ne, ns, nion, nlm)) + return dos + \ No newline at end of file