dft_tools/python/converters/plovasp/proj_shell.py

 
################################################################################
#
# TRIQS: a Toolbox for Research in Interacting Quantum Systems
#
# Copyright (C) 2011 by M. Ferrero, O. Parcollet
#
# DFT tools: Copyright (C) 2011 by M. Aichhorn, L. Pourovskii, V. Vildosola
#
# PLOVasp: Copyright (C) 2015 by O. E. Peil
#
# TRIQS is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# TRIQS. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
r"""
    vasp.proj_shell
    ===============

    Storage and manipulation on projector shells.
"""
def issue_warning(message):
    """
    Issues a warning.
    """
    print
    print "  !!! WARNING !!!: " + message
    print

import itertools as it
import numpy as np
try:
    import atm
    atmlib_present = True
except ImportError:
    issue_warning("Error importing ATM libray, DOS calculation will fail!")
    atmlib_present = False

np.set_printoptions(suppress=True)

################################################################################
################################################################################
#
# 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, kmesh, structure, 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)
        self.nion = len(self.ion_list)

        self.extract_projectors(proj_raw, proj_params, kmesh, structure)

################################################################################
#
# 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'.

        Flag 'self.do_transform' is introduced for the optimization purposes
        to avoid superfluous matrix multiplications.
        """
        nion = len(self.ion_list)
        nm = self.lm2 - self.lm1

        if 'tmatrices' in sh_pars:
            self.do_transform = True

            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:
            self.do_transform = True

            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
        self.do_transform = False

        ns_dim = 2 if self.nc_flag else 1
        ndim = nm * ns_dim

# We still need the matrices for the output
        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

################################################################################
#
# extract_projectors
#
################################################################################
    def extract_projectors(self, proj_raw, proj_params, kmesh, structure):
        """
        Extracts projectors for the given shell.
        
        Projectors are selected from the raw-projector array 'proj_raw'
        according to the shell parameters.
        If necessary the projectors are transformed usin 'self.tmatrices'.
        """
        assert self.nc_flag == False, "Non-collinear case is not implemented"

        nion = len(self.ion_list)
        nlm = self.lm2 - self.lm1
        _, ns, nk, nb = proj_raw.shape

        if self.do_transform:
# TODO: implement a non-collinear case
#       for a non-collinear case 'ndim' is 'ns * nm'
            ndim = self.tmatrices.shape[1]
            self.proj_arr = np.zeros((nion, ns, nk, ndim, nb), dtype=np.complex128)
            for ik in xrange(nk):
                kp = kmesh['kpoints'][ik]
                for io, ion in enumerate(self.ion_list):
                    proj_k = np.zeros((ns, nlm, nb), dtype=np.complex128)
                    qcoord = structure['qcoords'][ion]
#                    kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
#                    kphase = 1.0
                    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:
                                proj_k[:, m, :] = proj_raw[ip, :, ik, :]  #* kphase
                                break
                    for isp in xrange(ns):
                        self.proj_arr[io, isp, ik, :, :] = np.dot(self.tmatrices[io, :, :], proj_k[isp, :, :])

        else:
# No transformation: just copy the projectors as they are
            self.proj_arr = np.zeros((nion, ns, nk, nlm, nb), dtype=np.complex128)
            for io, ion in enumerate(self.ion_list):
                qcoord = structure['qcoords'][ion]
                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:
                            self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :]
#                            for ik in xrange(nk):
#                                kp = kmesh['kpoints'][ik]
##                                kphase = np.exp(-2.0j * np.pi * np.dot(kp, qcoord))
#                                kphase = 1.0
#                                self.proj_arr[io, :, :, m, :] = proj_raw[ip, :, :, :] # * kphase
                            break


################################################################################
#
# 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"

        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)

#        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
        if site_diag:
            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
        else:
            proj_k = np.zeros((ndim, nbtot), dtype=np.complex128)
            for isp in xrange(ns):
                for ik, weight, occ in it.izip(it.count(), kweights, occnums[isp, :, :]):
                    for io in xrange(nion):
                        i1 = io * nlm
                        i2 = (io + 1) * nlm
                        proj_k[i1:i2, :] = self.proj_win[io, isp, ik, ...]
                    occ_mats[isp, 0, :, :] += np.dot(proj_k * occ[ib1:ib2],
                                                 proj_k.conj().T).real * weight
                    overlaps[isp, 0, :, :] += 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

################################################################################
#
# local_hamiltonian
#
################################################################################
    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)

#        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, eigk in it.izip(it.count(), kweights, occnums[isp, :, :],
                                          el_struct.eigvals[:, ib1:ib2, isp]):
                for io in xrange(nion):
                    proj_k = self.proj_win[io, isp, ik, ...]
                    loc_ham[isp, io, :, :] += np.dot(proj_k * (eigk - el_struct.efermi),
                                                 proj_k.conj().T).real * weight

#        if not symops is None:
#            occ_mats = symmetrize_matrix_set(occ_mats, symops, ions, perm_map)

        return loc_ham

################################################################################
#
# 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

        assert atmlib_present, "ATM library was not imported; cannot calculate DOS"
# 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 = atm.dos_tetra_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