diff --git a/python/vasp/plotools.py b/python/vasp/plotools.py index 321173ab..c1c087f1 100644 --- a/python/vasp/plotools.py +++ b/python/vasp/plotools.py @@ -1,7 +1,6 @@ import itertools as it import numpy as np -import vasp.atm.c_atm_dos as c_atm_dos np.set_printoptions(suppress=True) @@ -50,39 +49,6 @@ class Projector: return (overlap, over_eig) -################################################################################ -# -# orthogonalize_projector_matrix() -# -################################################################################ -def orthogonalize_projector_matrix(p_matrix): - """ - Orthogonalizes a projector defined by a rectangular matrix `p_matrix`. - - Parameters - ---------- - - p_matrix (numpy.array[complex]) : matrix `Nm x Nb`, where `Nm` is - the number of orbitals, `Nb` number of bands - - Returns - ------- - - Orthogonalized projector matrix, initial overlap matrix and its eigenvalues. - """ -# Overlap matrix O_{m m'} = \sum_{v} P_{m v} P^{*}_{v m'} - overlap = np.dot(p_matrix, p_matrix.conj().T) -# Calculate [O^{-1/2}]_{m m'} - eig, eigv = np.linalg.eigh(overlap) - assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:" - "projectors are ill-defined") - sqrt_eig = np.diag(1.0 / np.sqrt(eig)) - shalf = np.dot(eigv, np.dot(sqrt_eig, eigv.conj().T)) -# Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v} - p_ortho = np.dot(shalf, p_matrix) - - return (p_ortho, overlap, eig) - ################################################################################ # check_data_consistency() ################################################################################ @@ -108,481 +74,6 @@ def check_data_consistency(pars, el_struct): raise Exception(errmsg) -################################################################################ -# select_bands() -################################################################################ -def select_bands(eigvals, emin, emax): - """ - Select a subset of bands lying within a given energy window. - The band energies are assumed to be sorted in an ascending order. - - Parameters - ---------- - - eigvals (numpy.array) : all eigenvalues - emin, emax (float) : energy window - - Returns - ------- - - ib_win, nb_min, nb_max : - """ -# Sanity check - if emin > eigvals.max() or emax < eigvals.min(): - raise Exception("Energy window does not overlap with the band structure") - - nk, nband, ns_band = eigvals.shape - ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32) - - ib_min = 10000000 - ib_max = 0 - for isp in xrange(ns_band): - for ik in xrange(nk): - for ib in xrange(nband): - en = eigvals[ik, ib, isp] - if en >= emin: - break - ib1 = ib - for ib in xrange(ib1, nband): - en = eigvals[ik, ib, isp] - if en > emax: - break - else: -# If we reached the last band add 1 to get the correct bound - ib += 1 - ib2 = ib - 1 - - assert ib1 <= ib2, "No bands inside the window for ik = %s"%(ik) - - ib_win[ik, isp, 0] = ib1 - ib_win[ik, isp, 1] = ib2 - - ib_min = min(ib_min, ib1) - ib_max = max(ib_max, ib2) - - return ib_win, ib_min, ib_max - -################################################################################ -################################################################################ -# -# class ProjectorGroup -# -################################################################################ -################################################################################ -class ProjectorGroup: - """ - Container of projectors defined within a certain energy window. - - The constructor selects a subset of projectors according to - the parameters from the config-file (passed in `pars`). - - Parameters: - - - gr_pars (dict) : group parameters from the config-file - - shells ([ProjectorShell]) : array of ProjectorShell objects - - eigvals (numpy.array) : array of KS eigenvalues - - """ - def __init__(self, gr_pars, shells, eigvals, ferw): - """ - Constructor - """ - self.emin, self.emax = gr_pars['ewindow'] - self.ishells = gr_pars['shells'] - self.ortho = gr_pars['normalize'] - self.normion = gr_pars['normion'] - - self.shells = shells - -# Determine the minimum and maximum band numbers - ib_win, ib_min, ib_max = select_bands(eigvals, self.emin, self.emax) - self.ib_win = ib_win - self.ib_min = ib_min - self.ib_max = ib_max - self.nb_max = ib_max - ib_min + 1 - -# Select projectors within the energy window - for ish in self.ishells: - shell = self.shells[ish] - shell.select_projectors(ib_win, ib_min, ib_max) - - - -################################################################################ -# -# nelect_window -# -################################################################################ - def nelect_window(self, el_struct): - """ - Determines the total number of electrons within the window. - """ - self.nelect = 0 - nk, ns_band, _ = self.ib_win.shape - rspin = 2.0 if ns_band == 1 else 1.0 - for isp in xrange(ns_band): - for ik in xrange(nk): - ib1 = self.ib_win[ik, isp, 0] - ib2 = self.ib_win[ik, isp, 1] - occ = el_struct.ferw[isp, ik, ib1:ib2] - kwght = el_struct.kmesh['kweights'][ik] - self.nelect += occ.sum() * kwght * rspin - - return self.nelect - -################################################################################ -# -# orthogonalize -# -################################################################################ - def orthogonalize(self): - """ - Orthogonalize a group of projectors. - """ -# Quick exit if no normalization is requested - if not self.ortho: - return - -# TODO: add the case of 'normion = True' - assert not self.normion, "'NORMION = True' is not yet implemented" - -# Determine the dimension of the projector matrix -# and map the blocks to the big matrix - i1_bl = 0 - bl_map = [{} for ish in self.ishells] - for ish in self.ishells: - _shell = self.shells[ish] - nion, ns, nk, nlm, nb_max = _shell.proj_win.shape - bmat_bl = [] # indices corresponding to a big block matrix - for ion in xrange(nion): - i2_bl = i1_bl + nlm - bmat_bl.append((i1_bl, i2_bl)) - i1_bl = i2_bl - bl_map[ish]['bmat_blocks'] = bmat_bl - - ndim = i2_bl - p_mat = np.zeros((ndim, nb_max), dtype=np.complex128) - for isp in xrange(ns): - for ik in xrange(nk): - nb = self.ib_win[ik, isp, 1] - self.ib_win[ik, isp, 0] + 1 -# Combine all projectors of the group to one block projector - for ish in self.ishells: - shell = self.shells[ish] - blocks = bl_map[ish]['bmat_blocks'] - for ion in xrange(nion): - i1, i2 = blocks[ion] - p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb] -# Now orthogonalize the obtained block projector - p_orth, overl, eig = orthogonalize_projector_matrix(p_mat) -# print "ik = ", ik -# print overl.real -# Distribute back projectors in the same order - for ish in self.ishells: - shell = self.shells[ish] - blocks = bl_map[ish]['bmat_blocks'] - for ion in xrange(nion): - i1, i2 = blocks[ion] - shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb] - - -################################################################################ -################################################################################ -# -# class ProjectorShell -# -################################################################################ -################################################################################ -class ProjectorShell: - """ - Container of projectors related to a specific shell. - - The constructor pre-selects a subset of projectors according to - the shell parameters passed from the config-file. - - Parameters: - - - sh_pars (dict) : shell parameters from the config-file - - proj_raw (numpy.array) : array of raw projectors - - """ - def __init__(self, sh_pars, proj_raw, proj_params, nc_flag): - self.lorb = sh_pars['lshell'] - self.ion_list = sh_pars['ion_list'] - self.user_index = sh_pars['user_index'] - self.nc_flag = nc_flag -# try: -# self.tmatrix = sh_pars['tmatrix'] -# except KeyError: -# self.tmatrix = None - - self.lm1 = self.lorb**2 - self.lm2 = (self.lorb+1)**2 - - self.ndim = self.extract_tmatrices(sh_pars) -# if self.tmatrix is None: -# self.ndim = self.lm2 - self.lm1 -# else: -## TODO: generalize this to a tmatrix for every ion -# self.ndim = self.tmatrix.shape[0] - -# Pre-select a subset of projectors (this should be an array view => no memory is wasted) -# !!! This sucks but I have to change the order of 'ib' and 'ilm' indices here -# This should perhaps be done right after the projector array is read from PLOCAR -# self.proj_arr = proj_raw[self.ion_list, :, :, :, self.lm1:self.lm2].transpose((0, 1, 2, 4, 3)) -# We want to select projectors from 'proj_raw' and form an array -# self.proj_arr[nion, ns, nk, nlm, nb] -# TODO: think of a smart way of copying the selected projectors -# perhaps, by redesigning slightly the structure of 'proj_arr' and 'proj_win' -# or by storing only a mapping between site/orbitals and indices of 'proj_raw' -# iproj_l = [] - nion = len(self.ion_list) - nlm = self.lm2 - self.lm1 - _, ns, nk, nb = proj_raw.shape - self.proj_arr = np.zeros((nion, ns, nk, nlm, nb), dtype=np.complex128) - for io, ion in enumerate(self.ion_list): - for m in xrange(nlm): -# Here we search for the index of the projector with the given isite/l/m indices - for ip, par in enumerate(proj_params): - if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m: -# iproj_l.append(ip) - self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :] - break - -# self.proj_arr = proj_raw[iproj_l, :, :, :].transpose((1, 2, 0, 3)) - -################################################################################ -# -# extract_tmatrices -# -################################################################################ - def extract_tmatrices(self, sh_pars): - """ - Extracts and interprets transformation matrices provided by the - config-parser. - There are two relevant options in 'sh_pars': - - 'tmatrix' : a transformation matrix applied to all ions in the shell - 'tmatrices': interpreted as a set of transformation matrices for each ion. - - If both of the options are present a warning is issued and 'tmatrices' - supersedes 'tmatrix'. - """ - nion = len(self.ion_list) - nm = self.lm2 - self.lm1 - - if 'tmatrices' in sh_pars: - if 'tmatrix' in sh_pars: - mess = "Both TRANSFORM and TRANSFILE are specified, TRANSFORM will be ignored." - issue_warning(mess) - - raw_matrices = sh_pars['tmatrices'] - nrow, ncol = raw_matrices.shape - - assert nrow%nion == 0, "Number of rows in TRANSFILE must be divisible by the number of ions" - assert ncol%nm == 0, "Number of columns in TRANSFILE must be divisible by the number of orbitals 2*l + 1" - - nr = nrow / nion - nsize = ncol / nm - assert nsize in (1, 2, 4), "Number of columns in TRANSFILE must be divisible by either 1, 2, or 4" -# -# Determine the spin-dimension and whether the matrices are real or complex -# -# if nsize == 1 or nsize == 2: -# Matrices (either real or complex) are spin-independent -# nls_dim = nm -# if msize == 2: -# is_complex = True -# else: -# is_complex = False -# elif nsize = 4: -# Matrices are complex and spin-dependent -# nls_dim = 2 * nm -# is_complex = True -# - is_complex = nsize > 1 - ns_dim = max(1, nsize / 2) - -# Dimension of the orbital subspace - assert nr%ns_dim == 0, "Number of rows in TRANSFILE is not compatible with the spin dimension" - ndim = nr / ns_dim - - self.tmatrices = np.zeros((nion, nr, nm * ns_dim), dtype=np.complex128) - - if is_complex: - raw_matrices = raw_matrices[:, ::2] + raw_matrices[:, 1::2] * 1j - - for io in xrange(nion): - i1 = io * nr - i2 = (io + 1) * nr - self.tmatrices[io, :, :] = raw_matrices[i1:i2, :] - - return ndim - - if 'tmatrix' in sh_pars: - raw_matrix = sh_pars['tmatrix'] - nrow, ncol = raw_matrix.shape - - assert ncol%nm == 0, "Number of columns in TRANSFORM must be divisible by the number of orbitals 2*l + 1" - -# Only spin-independent matrices are expected here - nsize = ncol / nm - assert nsize in (1, 2), "Number of columns in TRANSFORM must be divisible by either 1 or 2" - - is_complex = nsize > 1 - if is_complex: - matrix = raw_matrix[:, ::2] + raw_matrix[:, 1::2] * 1j - else: - matrix = raw_matrix - - ndim = nrow - - self.tmatrices = np.zeros((nion, nrow, nm), dtype=np.complex128) - for io in xrange(nion): - self.tmatrices[io, :, :] = raw_matrix - - return ndim - -# If no transformation matrices are provided define a default one - ns_dim = 2 if self.nc_flag else 1 - ndim = nm * ns_dim - - self.tmatrices = np.zeros((nion, ndim, ndim), dtype=np.complex128) - for io in xrange(nion): - self.tmatrices[io, :, :] = np.identity(ndim, dtype=np.complex128) - - return ndim - -################################################################################ -# -# select_projectors -# -################################################################################ - def select_projectors(self, ib_win, ib_min, ib_max): - """ - Selects a subset of projectors corresponding to a given energy window. - """ - self.ib_win = ib_win - self.ib_min = ib_min - self.ib_max = ib_max - nb_max = ib_max - ib_min + 1 - -# Set the dimensions of the array - nion, ns, nk, nlm, nbtot = self.proj_arr.shape -# !!! Note that the order of the two last indices is different !!! - self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128) - -# Select projectors for a given energy window - ns_band = self.ib_win.shape[1] - for isp in xrange(ns): - for ik in xrange(nk): -# TODO: for non-collinear case something else should be done here - is_b = min(isp, ns_band) - ib1 = self.ib_win[ik, is_b, 0] - ib2 = self.ib_win[ik, is_b, 1] + 1 - ib_win = ib2 - ib1 - self.proj_win[:, isp, ik, :, :ib_win] = self.proj_arr[:, isp, ik, :, ib1:ib2] - -################################################################################ -# -# density_matrix -# -################################################################################ - 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" - - occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) - overlaps = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) - -# self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128) - kweights = el_struct.kmesh['kweights'] - occnums = el_struct.ferw - ib1 = self.ib_min - ib2 = self.ib_max + 1 - for isp in xrange(ns): - for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]): - for io in xrange(nion): - proj_k = self.proj_win[io, isp, ik, ...] - occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], - proj_k.conj().T).real * weight - overlaps[isp, io, :, :] += np.dot(proj_k, - proj_k.conj().T).real * weight - -# if not symops is None: -# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map) - - return occ_mats, overlaps - -################################################################################ -# -# density_of_states -# -################################################################################ - def density_of_states(self, el_struct, emesh): - """ - Returns projected DOS for the shell. - """ - nion, ns, nk, nlm, nbtot = self.proj_win.shape - -# There is a problem with data storage structure of projectors that will -# make life more complicated. The problem is that band-indices of projectors -# for different k-points do not match because we store 'nb_max' values starting -# from 0. - nb_max = self.ib_max - self.ib_min + 1 - ns_band = self.ib_win.shape[1] - - ne = len(emesh) - dos = np.zeros((ne, ns, nion, nlm)) - w_k = np.zeros((nk, nb_max, ns, nion, nlm), dtype=np.complex128) - for isp in xrange(ns): - for ik in xrange(nk): - is_b = min(isp, ns_band) - ib1 = self.ib_win[ik, is_b, 0] - ib2 = self.ib_win[ik, is_b, 1] + 1 - for ib_g in xrange(ib1, ib2): - for io in xrange(nion): -# Note the difference between 'ib' and 'ibn': -# 'ib' counts from 0 to 'nb_k - 1' -# 'ibn' counts from 'ib1 - ib_min' to 'ib2 - ib_min' - ib = ib_g - ib1 - ibn = ib_g - self.ib_min - proj_k = self.proj_win[io, isp, ik, :, ib] - w_k[ik, ib, isp, io, :] = proj_k * proj_k.conj() - -# eigv_ef = el_struct.eigvals[ik, ib, isp] - el_struct.efermi - itt = el_struct.kmesh['itet'].T -# k-indices are starting from 0 in Python - itt[1:, :] -= 1 - for isp in xrange(ns): - for ib, eigk in enumerate(el_struct.eigvals[:, self.ib_min:self.ib_max+1, isp].T): - for ie, e in enumerate(emesh): - eigk_ef = eigk - el_struct.efermi - cti = c_atm_dos.dos_weights_3d(eigk_ef, e, itt) - for im in xrange(nlm): - for io in xrange(nion): - dos[ie, isp, io, im] += np.sum((cti * w_k[itt[1:, :], ib, isp, io, im].real).sum(0) * itt[0, :]) - - dos *= 2 * el_struct.kmesh['volt'] -# for isp in xrange(ns): -# for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]): -# for io in xrange(nion): -# proj_k = self.proj_win[isp, io, ik, ...] -# occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], -# proj_k.conj().T).real * weight -# overlaps[isp, io, :, :] += np.dot(proj_k, -# proj_k.conj().T).real * weight - -# if not symops is None: -# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map) - - return dos - - ################################################################################ # diff --git a/python/vasp/proj_group.py b/python/vasp/proj_group.py new file mode 100644 index 00000000..53f98698 --- /dev/null +++ b/python/vasp/proj_group.py @@ -0,0 +1,269 @@ + +import numpy as np + +np.set_printoptions(suppress=True) + +################################################################################ +# +# orthogonalize_projector_matrix() +# +################################################################################ +def orthogonalize_projector_matrix(p_matrix): + """ + Orthogonalizes a projector defined by a rectangular matrix `p_matrix`. + + Parameters + ---------- + + p_matrix (numpy.array[complex]) : matrix `Nm x Nb`, where `Nm` is + the number of orbitals, `Nb` number of bands + + Returns + ------- + + Orthogonalized projector matrix, initial overlap matrix and its eigenvalues. + """ +# Overlap matrix O_{m m'} = \sum_{v} P_{m v} P^{*}_{v m'} + overlap = np.dot(p_matrix, p_matrix.conj().T) +# Calculate [O^{-1/2}]_{m m'} + eig, eigv = np.linalg.eigh(overlap) + assert np.all(eig > 0.0), ("Negative eigenvalues of the overlap matrix:" + "projectors are ill-defined") + sqrt_eig = np.diag(1.0 / np.sqrt(eig)) + shalf = np.dot(eigv, np.dot(sqrt_eig, eigv.conj().T)) +# Apply \tilde{P}_{m v} = \sum_{m'} [O^{-1/2}]_{m m'} P_{m' v} + p_ortho = np.dot(shalf, p_matrix) + + return (p_ortho, overlap, eig) + +################################################################################ +# +# select_bands() +# +################################################################################ +def select_bands(eigvals, emin, emax): + """ + Select a subset of bands lying within a given energy window. + The band energies are assumed to be sorted in an ascending order. + + Parameters + ---------- + + eigvals (numpy.array) : all eigenvalues + emin, emax (float) : energy window + + Returns + ------- + + ib_win, nb_min, nb_max : + """ +# Sanity check + if emin > eigvals.max() or emax < eigvals.min(): + raise Exception("Energy window does not overlap with the band structure") + + nk, nband, ns_band = eigvals.shape + ib_win = np.zeros((nk, ns_band, 2), dtype=np.int32) + + ib_min = 10000000 + ib_max = 0 + for isp in xrange(ns_band): + for ik in xrange(nk): + for ib in xrange(nband): + en = eigvals[ik, ib, isp] + if en >= emin: + break + ib1 = ib + for ib in xrange(ib1, nband): + en = eigvals[ik, ib, isp] + if en > emax: + break + else: +# If we reached the last band add 1 to get the correct bound + ib += 1 + ib2 = ib - 1 + + assert ib1 <= ib2, "No bands inside the window for ik = %s"%(ik) + + ib_win[ik, isp, 0] = ib1 + ib_win[ik, isp, 1] = ib2 + + ib_min = min(ib_min, ib1) + ib_max = max(ib_max, ib2) + + return ib_win, ib_min, ib_max + +################################################################################ +################################################################################ +# +# class ProjectorGroup +# +################################################################################ +################################################################################ +class ProjectorGroup: + """ + Container of projectors defined within a certain energy window. + + The constructor selects a subset of projectors according to + the parameters from the config-file (passed in `pars`). + + Parameters: + + - gr_pars (dict) : group parameters from the config-file + - shells ([ProjectorShell]) : array of ProjectorShell objects + - eigvals (numpy.array) : array of KS eigenvalues + + """ + def __init__(self, gr_pars, shells, eigvals, ferw): + """ + Constructor + """ + self.emin, self.emax = gr_pars['ewindow'] + self.ishells = gr_pars['shells'] + self.ortho = gr_pars['normalize'] + self.normion = gr_pars['normion'] + + self.shells = shells + +# Determine the minimum and maximum band numbers + ib_win, ib_min, ib_max = select_bands(eigvals, self.emin, self.emax) + self.ib_win = ib_win + self.ib_min = ib_min + self.ib_max = ib_max + self.nb_max = ib_max - ib_min + 1 + +# Select projectors within the energy window + for ish in self.ishells: + shell = self.shells[ish] + shell.select_projectors(ib_win, ib_min, ib_max) + + + +################################################################################ +# +# nelect_window +# +################################################################################ + def nelect_window(self, el_struct): + """ + Determines the total number of electrons within the window. + """ + self.nelect = 0 + nk, ns_band, _ = self.ib_win.shape + rspin = 2.0 if ns_band == 1 else 1.0 + for isp in xrange(ns_band): + for ik in xrange(nk): + ib1 = self.ib_win[ik, isp, 0] + ib2 = self.ib_win[ik, isp, 1] + occ = el_struct.ferw[isp, ik, ib1:ib2] + kwght = el_struct.kmesh['kweights'][ik] + self.nelect += occ.sum() * kwght * rspin + + return self.nelect + +################################################################################ +# +# orthogonalize +# +################################################################################ + def orthogonalize(self): + """ + Orthogonalize a group of projectors. + + There are two options for orthogonalizing projectors: + 1. one ensures orthogonality on each site (NORMION = True); + 2. one ensures orthogonality for subsets of sites (NORMION = False), + as, e.g., in cluster calculations. + + In order to handle various cases the strategy is first to build a + mapping that selects appropriate blocks of raw projectors, forms a + matrix consisting of these blocks, orthogonalize the matrix, and use + the mapping again to write the orthogonalized projectors back to the + projector arrays. Note that the blocks can comprise several projector arrays + contained in different projector shells. + + Mapping is defined as a list of 'block_maps' corresponding to subsets + of projectors to be orthogonalized. Each subset corresponds to a subset of sites + and spans all orbital indices. defined by 'bl_map' as + + bl_map = [((i1_start, i1_end), (i1_shell, i1_ion)), + ((i2_start, i2_end), (i2_shell, i2_ion)), + ...], + + where `iX_start`, `iX_end` is the range of indices of the block matrix + (in Python convention `iX_end = iX_last + 1`, with `iX_last` being the last index + of the range), + `iX_shell` and `iX_ion` the shell and site indices. The length of the range + should be consistent with 'nlm' dimensions of a corresponding shell, i.e., + `iX_end - iX_start = nlm[iX_shell]`. + + Consider particular cases: + 1. Orthogonality is ensured on each site (NORMION = True). + For each site 'ion' we have the following mapping: + + block_maps = [bl_map[ion] for ion in xrange(shell.nion) + for shell in shells] + + bl_map = [((i1_start, i1_end), (i1_shell, ion)), + ((i2_start, i2_end), (i2_shell, ion)), + ...], + + 2. Orthogonality is ensured on all sites within the group (NORMION = True). + The mapping: + + block_maps = [bl_map] + + bl_map = [((i1_start, i1_end), (i1_shell, i1_shell.ion1)), + ((i1_start, i1_end), (i1_shell, i1_shell.ion2)), + ... + ((i2_start, i2_end), (i2_shell, i2_shell.ion1)), + ((i2_start, i2_end), (i2_shell, i2_shell.ion2)), + ...], + + """ +# Quick exit if no normalization is requested + if not self.ortho: + return + +# TODO: add the case of 'normion = True' + assert not self.normion, "'NORMION = True' is not yet implemented" + +# Determine the dimension of the projector matrix +# and map the blocks to the big matrix + i1_bl = 0 + bl_map = [{} for ish in self.ishells] + for ish in self.ishells: + _shell = self.shells[ish] + nion, ns, nk, nlm, nb_max = _shell.proj_win.shape + bmat_bl = [] # indices corresponding to a big block matrix + for ion in xrange(nion): + i2_bl = i1_bl + nlm + bmat_bl.append((i1_bl, i2_bl)) + i1_bl = i2_bl + bl_map[ish]['bmat_blocks'] = bmat_bl + + ndim = i2_bl + p_mat = np.zeros((ndim, nb_max), dtype=np.complex128) + for isp in xrange(ns): + for ik in xrange(nk): + nb = self.ib_win[ik, isp, 1] - self.ib_win[ik, isp, 0] + 1 +# Combine all projectors of the group to one block projector + for ish in self.ishells: + shell = self.shells[ish] + blocks = bl_map[ish]['bmat_blocks'] + for ion in xrange(nion): + i1, i2 = blocks[ion] + p_mat[i1:i2, :nb] = shell.proj_win[ion, isp, ik, :nlm, :nb] +# Now orthogonalize the obtained block projector + p_orth, overl, eig = orthogonalize_projector_matrix(p_mat) +# print "ik = ", ik +# print overl.real +# Distribute back projectors in the same order + for ish in self.ishells: + shell = self.shells[ish] + blocks = bl_map[ish]['bmat_blocks'] + for ion in xrange(nion): + i1, i2 = blocks[ion] + shell.proj_win[ion, isp, ik, :nlm, :nb] = p_orth[i1:i2, :nb] + + + diff --git a/python/vasp/proj_shell.py b/python/vasp/proj_shell.py new file mode 100644 index 00000000..228f1107 --- /dev/null +++ b/python/vasp/proj_shell.py @@ -0,0 +1,314 @@ + +import itertools as it +import numpy as np +import vasp.atm.c_atm_dos as c_atm_dos + +np.set_printoptions(suppress=True) + +def issue_warning(message): + """ + Issues a warning. + """ + print + print " !!! WARNING !!!: " + message + print + +################################################################################ +################################################################################ +# +# class ProjectorShell +# +################################################################################ +################################################################################ +class ProjectorShell: + """ + Container of projectors related to a specific shell. + + The constructor pre-selects a subset of projectors according to + the shell parameters passed from the config-file. + + Parameters: + + - sh_pars (dict) : shell parameters from the config-file + - proj_raw (numpy.array) : array of raw projectors + + """ + def __init__(self, sh_pars, proj_raw, proj_params, nc_flag): + self.lorb = sh_pars['lshell'] + self.ion_list = sh_pars['ion_list'] + self.user_index = sh_pars['user_index'] + self.nc_flag = nc_flag +# try: +# self.tmatrix = sh_pars['tmatrix'] +# except KeyError: +# self.tmatrix = None + + self.lm1 = self.lorb**2 + self.lm2 = (self.lorb+1)**2 + + self.ndim = self.extract_tmatrices(sh_pars) +# if self.tmatrix is None: +# self.ndim = self.lm2 - self.lm1 +# else: +## TODO: generalize this to a tmatrix for every ion +# self.ndim = self.tmatrix.shape[0] + +# Pre-select a subset of projectors (this should be an array view => no memory is wasted) +# !!! This sucks but I have to change the order of 'ib' and 'ilm' indices here +# This should perhaps be done right after the projector array is read from PLOCAR +# self.proj_arr = proj_raw[self.ion_list, :, :, :, self.lm1:self.lm2].transpose((0, 1, 2, 4, 3)) +# We want to select projectors from 'proj_raw' and form an array +# self.proj_arr[nion, ns, nk, nlm, nb] +# TODO: think of a smart way of copying the selected projectors +# perhaps, by redesigning slightly the structure of 'proj_arr' and 'proj_win' +# or by storing only a mapping between site/orbitals and indices of 'proj_raw' +# iproj_l = [] + nion = len(self.ion_list) + nlm = self.lm2 - self.lm1 + _, ns, nk, nb = proj_raw.shape + self.proj_arr = np.zeros((nion, ns, nk, nlm, nb), dtype=np.complex128) + for io, ion in enumerate(self.ion_list): + for m in xrange(nlm): +# Here we search for the index of the projector with the given isite/l/m indices + for ip, par in enumerate(proj_params): + if par['isite'] - 1 == ion and par['l'] == self.lorb and par['m'] == m: +# iproj_l.append(ip) + self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :] + break + +# self.proj_arr = proj_raw[iproj_l, :, :, :].transpose((1, 2, 0, 3)) + +################################################################################ +# +# extract_tmatrices +# +################################################################################ + def extract_tmatrices(self, sh_pars): + """ + Extracts and interprets transformation matrices provided by the + config-parser. + There are two relevant options in 'sh_pars': + + 'tmatrix' : a transformation matrix applied to all ions in the shell + 'tmatrices': interpreted as a set of transformation matrices for each ion. + + If both of the options are present a warning is issued and 'tmatrices' + supersedes 'tmatrix'. + """ + nion = len(self.ion_list) + nm = self.lm2 - self.lm1 + + if 'tmatrices' in sh_pars: + if 'tmatrix' in sh_pars: + mess = "Both TRANSFORM and TRANSFILE are specified, TRANSFORM will be ignored." + issue_warning(mess) + + raw_matrices = sh_pars['tmatrices'] + nrow, ncol = raw_matrices.shape + + assert nrow%nion == 0, "Number of rows in TRANSFILE must be divisible by the number of ions" + assert ncol%nm == 0, "Number of columns in TRANSFILE must be divisible by the number of orbitals 2*l + 1" + + nr = nrow / nion + nsize = ncol / nm + assert nsize in (1, 2, 4), "Number of columns in TRANSFILE must be divisible by either 1, 2, or 4" +# +# Determine the spin-dimension and whether the matrices are real or complex +# +# if nsize == 1 or nsize == 2: +# Matrices (either real or complex) are spin-independent +# nls_dim = nm +# if msize == 2: +# is_complex = True +# else: +# is_complex = False +# elif nsize = 4: +# Matrices are complex and spin-dependent +# nls_dim = 2 * nm +# is_complex = True +# + is_complex = nsize > 1 + ns_dim = max(1, nsize / 2) + +# Dimension of the orbital subspace + assert nr%ns_dim == 0, "Number of rows in TRANSFILE is not compatible with the spin dimension" + ndim = nr / ns_dim + + self.tmatrices = np.zeros((nion, nr, nm * ns_dim), dtype=np.complex128) + + if is_complex: + raw_matrices = raw_matrices[:, ::2] + raw_matrices[:, 1::2] * 1j + + for io in xrange(nion): + i1 = io * nr + i2 = (io + 1) * nr + self.tmatrices[io, :, :] = raw_matrices[i1:i2, :] + + return ndim + + if 'tmatrix' in sh_pars: + raw_matrix = sh_pars['tmatrix'] + nrow, ncol = raw_matrix.shape + + assert ncol%nm == 0, "Number of columns in TRANSFORM must be divisible by the number of orbitals 2*l + 1" + +# Only spin-independent matrices are expected here + nsize = ncol / nm + assert nsize in (1, 2), "Number of columns in TRANSFORM must be divisible by either 1 or 2" + + is_complex = nsize > 1 + if is_complex: + matrix = raw_matrix[:, ::2] + raw_matrix[:, 1::2] * 1j + else: + matrix = raw_matrix + + ndim = nrow + + self.tmatrices = np.zeros((nion, nrow, nm), dtype=np.complex128) + for io in xrange(nion): + self.tmatrices[io, :, :] = raw_matrix + + return ndim + +# If no transformation matrices are provided define a default one + ns_dim = 2 if self.nc_flag else 1 + ndim = nm * ns_dim + + self.tmatrices = np.zeros((nion, ndim, ndim), dtype=np.complex128) + for io in xrange(nion): + self.tmatrices[io, :, :] = np.identity(ndim, dtype=np.complex128) + + return ndim + +################################################################################ +# +# select_projectors +# +################################################################################ + def select_projectors(self, ib_win, ib_min, ib_max): + """ + Selects a subset of projectors corresponding to a given energy window. + """ + self.ib_win = ib_win + self.ib_min = ib_min + self.ib_max = ib_max + nb_max = ib_max - ib_min + 1 + +# Set the dimensions of the array + nion, ns, nk, nlm, nbtot = self.proj_arr.shape +# !!! Note that the order of the two last indices is different !!! + self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128) + +# Select projectors for a given energy window + ns_band = self.ib_win.shape[1] + for isp in xrange(ns): + for ik in xrange(nk): +# TODO: for non-collinear case something else should be done here + is_b = min(isp, ns_band) + ib1 = self.ib_win[ik, is_b, 0] + ib2 = self.ib_win[ik, is_b, 1] + 1 + ib_win = ib2 - ib1 + self.proj_win[:, isp, ik, :, :ib_win] = self.proj_arr[:, isp, ik, :, ib1:ib2] + +################################################################################ +# +# density_matrix +# +################################################################################ + 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" + + occ_mats = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) + overlaps = np.zeros((ns, nion, nlm, nlm), dtype=np.float64) + +# self.proj_win = np.zeros((nion, ns, nk, nlm, nb_max), dtype=np.complex128) + kweights = el_struct.kmesh['kweights'] + occnums = el_struct.ferw + ib1 = self.ib_min + ib2 = self.ib_max + 1 + for isp in xrange(ns): + for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]): + for io in xrange(nion): + proj_k = self.proj_win[io, isp, ik, ...] + occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], + proj_k.conj().T).real * weight + overlaps[isp, io, :, :] += np.dot(proj_k, + proj_k.conj().T).real * weight + +# if not symops is None: +# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map) + + return occ_mats, overlaps + +################################################################################ +# +# density_of_states +# +################################################################################ + def density_of_states(self, el_struct, emesh): + """ + Returns projected DOS for the shell. + """ + nion, ns, nk, nlm, nbtot = self.proj_win.shape + +# There is a problem with data storage structure of projectors that will +# make life more complicated. The problem is that band-indices of projectors +# for different k-points do not match because we store 'nb_max' values starting +# from 0. + nb_max = self.ib_max - self.ib_min + 1 + ns_band = self.ib_win.shape[1] + + ne = len(emesh) + dos = np.zeros((ne, ns, nion, nlm)) + w_k = np.zeros((nk, nb_max, ns, nion, nlm), dtype=np.complex128) + for isp in xrange(ns): + for ik in xrange(nk): + is_b = min(isp, ns_band) + ib1 = self.ib_win[ik, is_b, 0] + ib2 = self.ib_win[ik, is_b, 1] + 1 + for ib_g in xrange(ib1, ib2): + for io in xrange(nion): +# Note the difference between 'ib' and 'ibn': +# 'ib' counts from 0 to 'nb_k - 1' +# 'ibn' counts from 'ib1 - ib_min' to 'ib2 - ib_min' + ib = ib_g - ib1 + ibn = ib_g - self.ib_min + proj_k = self.proj_win[io, isp, ik, :, ib] + w_k[ik, ib, isp, io, :] = proj_k * proj_k.conj() + +# eigv_ef = el_struct.eigvals[ik, ib, isp] - el_struct.efermi + itt = el_struct.kmesh['itet'].T +# k-indices are starting from 0 in Python + itt[1:, :] -= 1 + for isp in xrange(ns): + for ib, eigk in enumerate(el_struct.eigvals[:, self.ib_min:self.ib_max+1, isp].T): + for ie, e in enumerate(emesh): + eigk_ef = eigk - el_struct.efermi + cti = c_atm_dos.dos_weights_3d(eigk_ef, e, itt) + for im in xrange(nlm): + for io in xrange(nion): + dos[ie, isp, io, im] += np.sum((cti * w_k[itt[1:, :], ib, isp, io, im].real).sum(0) * itt[0, :]) + + dos *= 2 * el_struct.kmesh['volt'] +# for isp in xrange(ns): +# for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]): +# for io in xrange(nion): +# proj_k = self.proj_win[isp, io, ik, ...] +# occ_mats[isp, io, :, :] += np.dot(proj_k * occ[ib1:ib2], +# proj_k.conj().T).real * weight +# overlaps[isp, io, :, :] += np.dot(proj_k, +# proj_k.conj().T).real * weight + +# if not symops is None: +# occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map) + + return dos + + + diff --git a/python/vasp/test/_inpconf/test_special_parsers.py b/python/vasp/test/_inpconf/test_special_parsers.py index efd83aab..b4ffc697 100644 --- a/python/vasp/test/_inpconf/test_special_parsers.py +++ b/python/vasp/test/_inpconf/test_special_parsers.py @@ -236,13 +236,13 @@ class TestParseFileTmatrix(arraytest.ArrayTestCase): [ -3.80200000e-04, 0.00000000e+00, 6.04452000e-02, -1.00000000e-07, -9.98171400e-01], [ -5.14500000e-04, -0.00000000e+00, -9.98171400e-01, 0.00000000e+00, -6.04450000e-02]]) - res = self.cpars.parse_file_tmatrix('tmatrix_file.dat') + res = self.cpars.parse_file_tmatrix(_rpath + 'tmatrix_file.dat') self.assertEqual(res, expected) # Scenario 2 def test_wrong_file(self): with self.assertRaises(ValueError): - self.cpars.parse_file_tmatrix('test1.cfg') + self.cpars.parse_file_tmatrix(_rpath + 'test1.cfg') ################################################################################ #