Source code for pdos_overlap.vasp_dos

# -*- coding: utf-8 -*-
"""
Created on Tue May 19 22:28:09 2020

@author: lansf
"""

from __future__ import absolute_import, division, print_function
import os
import pkg_resources
import numpy as np
from scipy.ndimage import gaussian_filter1d

[docs]def get_data_path(): """ Get default paths to package data. Returns ------- data_path : str path to package data """ data_path = pkg_resources.resource_filename(__name__, 'data/') return data_path
[docs]def get_example_data(): """ Get default path to vasp_dos data Returns ------- data_path : str path to example VASP data """ data_path = pkg_resources.resource_filename(__name__, 'data/pdos') return data_path
[docs]def get_all_VASP_files(directory): """ Get all DOSCAR and CONTCAR file paths Parameters ---------- directory : str root directory to look for DOSCAR files Returns ------- DOSCAR_files : list list of all DOSCAR files CONTCAR_files : list list of all CONTCAR files """ DOSCAR_directories = [os.path.join(r,subdirectory) for r,d,f in os.walk(directory) \ for subdirectory in d \ if 'DOSCAR' in os.listdir(os.path.join(r,subdirectory))] DOSCAR_files = [os.path.join(d,'DOSCAR') for d in DOSCAR_directories] CONTCAR_files = [os.path.join(d,'CONTCAR') for d in DOSCAR_directories] return DOSCAR_files, CONTCAR_files
[docs]def get_band_center(energies, densities, max_energy=None, min_energy=None, axis=-1): """ Get band center given energies and densities Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities max_energy : float cutoff energy (often the fermi energy) min_energy : float cutoff energy (often the fermi energy) axis : int axis of densities on which to integrate Returns ------- band_center : numpy.float64 or numpy.ndarray center of the band(s) up to max_energy Notes ----- trapezoidal rule is better for narrow gaussian peaks and for "rough" functions https://doi.org/10.1016/j.chemolab.2018.06.001 http://emis.icm.edu.pl/journals/JIPAM/v3n4/031_02.html """ densities = np.atleast_2d(densities) if max_energy is None: idx_stop = len(energies) else: idx_stop = (np.abs(energies-max_energy)).argmin()+1 if min_energy is None: idx_start = 0 else: idx_start = (np.abs(energies-min_energy)).argmin() Integrated_Energy = np.trapz(densities[:,idx_start:idx_stop]\ * energies[idx_start:idx_stop]\ , energies[idx_start:idx_stop], axis=axis) Integrated_Filling = np.trapz(densities[:,idx_start:idx_stop]\ , energies[idx_start:idx_stop], axis=axis) band_center = Integrated_Energy/Integrated_Filling return band_center.squeeze()
[docs]def get_bond_energy(energies, densities, e_fermi): """ Get bond energy from density Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities e_fermi : float highest occupied energy level Returns ------- bond_energy : numpy.float64 or numpy.ndarray total energy in the bonds """ densities = np.atleast_2d(densities) band_center = get_band_center(energies, densities) if len(np.array(band_center).shape) == 1: band_center = band_center.reshape(-1,1) integrand = (energies - band_center) * densities idx_stop = (np.abs(energies - e_fermi)).argmin()+1 bond_energy = 2 * np.trapz(integrand[:,0:idx_stop], energies[0:idx_stop]) return bond_energy.squeeze()
[docs]def get_band_width(energies, densities, fraction=0): """ Get band width given energies and densities Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities fraction : float maximum projected density considered part of a band Returns ------- band_width : numpy.float64 or numpy.ndarray width of the band(s) """ densities = np.atleast_2d(densities) index_values = np.arange(densities.shape[-1]) min_index = np.zeros(densities.shape[0]) max_index = np.zeros(densities.shape[0]) for count, density in enumerate(densities): min_index[count] = np.min(index_values[density > fraction * density.max()]) max_index[count] = np.max(index_values[density > fraction * density.max()]) band_width = energies[max_index.astype(int)] - energies[min_index.astype(int)] return band_width.squeeze()
def get_center_width(energies, densities, energy): """ Get width between to band centers given a division Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities energy : int energy that divides band centers Returns ------- center_width : numpy.float64 or numpy.ndarray width of two band centers separated by some ennergy """ center_lower = get_band_center(energies, densities, max_energy=energy) center_upper = get_band_center(energies, densities, min_energy=energy) center_width = center_upper - center_lower return center_width.squeeze()
[docs]def get_orbital_proximity(energies, densities, energy, moment=1): """ Get width between to band centers given a division Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities energy : int some energy with which the proximity of states is calculated moment : int moment with which to calculate the orbital proximity Returns ------- orbital_proximity : numpy.float64 or numpy.ndarray proximity of orbitals to some energy """ if moment == 1: integrand = np.abs(energies - energy) * densities else: integrand = (energies - energy)**moment * densities orbital_proximity = np.trapz(integrand, energies) / np.trapz(densities, energies) return orbital_proximity.squeeze()
[docs]def get_second_moment(energies, densities): """ Get the second moment Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities Returns ------- second_moment : numpy.float64 or numpy.ndarray second moment of the densities Notes ----- Utilizeds get_orbital_proximity """ band_center = get_band_center(energies, densities) if len(np.array(band_center).shape) == 1: band_center = band_center.reshape(-1,1) second_moment = get_orbital_proximity(energies, densities, band_center, moment=2) return second_moment.squeeze()
[docs]def get_filling(energies, densities, max_energy=None, min_energy=None, axis=-1): """ Get band center given energies and densities Parameters ---------- energies : numpy.ndarray discretized orbital energies densities : numpy.ndarray projected state densities max_energy : float cutoff energy (often the fermi energy) min_energy : float cutoff energy (often the fermi energy) axis : int axis of densities on which to integrate Returns ------- Integrated_Filling : float or numpy.ndarray Total filling between min and max energy Notes ----- trapezoidal rule is better for narrow gaussian peaks and for "rough" functions https://doi.org/10.1016/j.chemolab.2018.06.001 http://emis.icm.edu.pl/journals/JIPAM/v3n4/031_02.html """ densities = np.atleast_2d(densities) if max_energy is None: idx_stop = len(energies) else: idx_stop = (np.abs(energies-max_energy)).argmin()+1 if min_energy is None: idx_start = 0 else: idx_start = (np.abs(energies-min_energy)).argmin() Integrated_Filling = np.trapz(densities[:,idx_start:idx_stop]\ , energies[idx_start:idx_stop], axis=axis) return Integrated_Filling.squeeze()
[docs]class VASP_DOS: """Class for extracting projected density of states from VASP""" def __init__(self, file_name="DOSCAR",no_negatives=True, add_p2s=False): """ Parameters ---------- file_name : str full DOSCAR file location no_negatives : bool Indicates wheather negative occupances will be converted to zero. Negative occupances can occur if Methfessel-Paxton is used. Attributes ---------- no_negatives : bool Indicates wheather negative occupances will be converted to zero. Negative occupances can occur if Methfessel-Paxton is used. natoms : int number of atoms in the system emax : float maximum energy level emin : float minimum energy level ndos : int number of descritized energy levels e_fermi : float highest occupied energy level is_spin : bool indicates if projected density is spin resolved m_projected : bool indicates if projected density is orbital resolved orbital_dictionary: dict dictionary that maps resolved sublevels/orbitals to indices add_p2s : bool if true adds the p density to the s density (useful for metals) """ if file_name[-14:] == 'DOSCAR.lobster': natoms, emax, emin, ndos, e_fermi, is_spin, m_projected\ , orbital_dictionary = self._read_lobster(file_name=file_name, no_negatives=no_negatives) elif file_name[-6:] == 'DOSCAR': natoms, emax, emin, ndos, e_fermi, is_spin, m_projected\ , orbital_dictionary = self._read_doscar(file_name=file_name, no_negatives=no_negatives) self.no_negatives = no_negatives self.natoms = natoms self.emax = emax self.emin = emin self.ndos = ndos self.e_fermi = e_fermi self.is_spin = is_spin self.m_projected = m_projected self.add_p2s = add_p2s self.orbital_dictionary = orbital_dictionary if add_p2s == True and is_spin == False: for i in range(natoms): self._site_dos[i,orbital_dictionary['s'],:]\ += self.get_site_dos(i, ['p'], sum_density=True)[1][0] elif add_p2s == True and is_spin == True: for i in range(natoms): self._site_dos[i,orbital_dictionary['s+'],:]\ += self.get_site_dos(i, ['p+'], sum_density=True)[1][0] self._site_dos[i,orbital_dictionary['s-'],:]\ += self.get_site_dos(i, ['p-'], sum_density=True)[1][0]
[docs] def apply_gaussian_filter(self, sigma): """Applies Gaussian filter to self._total_dos and self._site_dos Parameters ---------- sigma : float standard deviation of Gaussian Kernel Attributes ---------- _total_dos_original : numpy.ndarray self._total_dos array without filter _site_dos_original : numpy.ndarray self._site_dos array without filter """ total_dos_exists = False try: self._total_dos total_dos_exists = True try: self._total_dos_original except: self._total_dos_original = self._total_dos.copy() except: pass site_dos_exists = False try: self._site_dos site_dos_exists = True try: self._site_dos_original except: self._site_dos_original = self._site_dos.copy() except: pass if total_dos_exists == True: self._total_dos[1:,:] = gaussian_filter1d(self._total_dos_original[1:,:].copy(), sigma) if site_dos_exists == True: self._site_dos[:,1:,:] = gaussian_filter1d(self._site_dos_original[:,1:,:].copy(), sigma)
[docs] def get_band_center(self, atom_indices, orbital_list, sum_density=False, sum_spin=True\ , max_energy=None, min_energy=None, axis=-1): """ Get band center for a given atom and list of orbitals Parameters ---------- atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. max_energy : float cutoff energy (often the fermi energy) min_energy : float cutoff energy (often the fermi energy) axis : int axis of densities on which to integrate Returns ------- band_center : float or numpy.ndarray center of the band(s) up to max_energy """ energies = self.get_energies() get_site_dos = self.get_site_dos orbitals, density = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) band_center = get_band_center(energies, density, max_energy=max_energy\ , min_energy = min_energy, axis=-1) return band_center
[docs] def get_filling(self, atom_indices, orbital_list, sum_density=False, sum_spin=True\ , max_energy=None, min_energy=None, axis=-1): """ Get band center for a given atom and list of orbitals Parameters ---------- atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. max_energy : float cutoff energy (often the fermi energy) min_energy : float cutoff energy (often the fermi energy) axis : int axis of densities on which to integrate Returns ------- band_center : float or numpy.ndarray center of the band(s) up to max_energy """ energies = self.get_energies() get_site_dos = self.get_site_dos orbitals, density = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) filling = get_filling(energies, density, max_energy=max_energy\ , min_energy = min_energy, axis=-1) return filling
[docs] def get_band_width(self, atom_indices, orbital_list, sum_density=False\ , sum_spin=True, fraction=0): """ Get band width given energies and densities Parameters ---------- atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. fraction : float maximum projected density considered part of a band Returns ------- band_width : float or numpy.ndarray width of the band(s) """ energies = self.get_energies() get_site_dos = self.get_site_dos orbitals, densities = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) band_width = get_band_width(energies, densities, fraction) return band_width
[docs] def get_center_width(self, energy, atom_indices, orbital_list, sum_density=False\ , sum_spin=True): """ Get width between to band centers given a division Parameters ---------- energy : int energy that divides band centers atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. Returns ------- center_width : float or numpy.ndarray width of two band centers separated by some ennergy """ energies = self.get_energies() get_site_dos = self.get_site_dos orbitals, densities = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) center_width = get_center_width(energies, densities, energy) return center_width
[docs] def get_orbital_proximity(self, energy, atom_indices, orbital_list, moment=1 , sum_density=False, sum_spin=True): """ Get width between to band centers given a division Parameters ---------- energy : int some energy with which the proximity of states is calculated moment : int moment with which to calculate the orbital proximity atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. Returns ------- orbital_proximity : float or numpy.ndarray proximity of orbitals to some energy """ energies = self.get_energies() get_site_dos = self.get_site_dos orbitals, densities = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) orbital_proximity = get_orbital_proximity(energies, densities, energy, moment) return orbital_proximity
[docs] def get_second_moment(self, atom_indices, orbital_list, sum_density=False\ , sum_spin=True): """ Get second moment of the density projected onto atomic orbitals Parameters ---------- atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. Returns ------- second_moment : float or numpy.ndarray width of two band centers separated by some ennergy """ energies = self.get_energies() get_site_dos = self.get_site_dos orbitals, densities = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) second_moment = get_second_moment(energies, densities) return second_moment
[docs] def get_bond_energy(self, atom_indices, orbital_list, sum_density=False\ , sum_spin=True): """ Get second moment of the density projected onto atomic orbitals Parameters ---------- atom_indices : list[int] list of atom indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual sub-level densities should be summed sum_spin : bool different spin densities are summed. Returns ------- bond_energy : float or numpy.ndarray bond enegy of the orbitals on atom_indices """ energies = self.get_energies() get_site_dos = self.get_site_dos e_fermi = self.e_fermi orbitals, densities = get_site_dos(atom_indices,orbital_list\ , sum_density=sum_density\ , sum_spin=sum_spin) second_moment = get_bond_energy(energies, densities, e_fermi) return second_moment
[docs] def get_energies(self): """ method for obtaining energy levels Returns ------- energies : numpy.ndarray 1-D array of energies """ energies = self._total_dos[0,:].copy() return energies
[docs] def get_total_dos(self, sum_spin=True): """ method for obtaining total density of states Returns ------- total_dos : numpy.ndarray 1-D or 2-D array of state densities """ _total_dos = self._total_dos if _total_dos.shape[0] == 3: total_dos = _total_dos[1, :] elif _total_dos.shape[0] == 5: if sum_spin == True: total_dos = _total_dos[1:3, :].sum(axis=0) else: total_dos = _total_dos[1:3, :] return total_dos
[docs] def get_integrated_dos(self, sum_spin=True): """ method for obtaining total integrated density of states Returns ------- integrated_dos : numpy.ndarray 1-D or 2-D array of state integrated densities """ _total_dos = self._total_dos if _total_dos.shape[0] == 3: integrated_dos = _total_dos[2, :] elif _total_dos.shape[0] == 5: if sum_spin == True: integrated_dos = _total_dos[3:5, :].sum(axis=0) else: integrated_dos = _total_dos[3:5, :] return integrated_dos
[docs] def get_site_dos(self, atom_indices, orbital_list=[], sum_density=False\ , sum_spin=True): """Return an NDOSxM array with dos for the chosen atom and orbital(s). Parameters ---------- atom_list : list[int] List of atom index/indices orbital_list : list[str] Which orbitals to return sum_density : bool if a sub-level is provided instead of an orbital, sum_density indicates if the individual orbital densities should be summed sum_spin : bool different spin densities are summed. Returns ------- new_orbital_list : list[str] If sum_density is True, new_orbital_list is the same as orbital_list. If sum_density is False, new_orbital_list is resolved by both orbital and spin (if available) projected_density : np.array Array of shape (len(new_orbital_list), ndos) """ # Integer indexing for orbitals starts from 1 in the _site_dos array # since the 0th column contains the energies orbital_dictionary = self.orbital_dictionary is_spin = self.is_spin m_projected = self.m_projected _site_dos = self._site_dos ndos = self.ndos try: len(atom_indices) except: atom_indices = [atom_indices] if len(orbital_list) == 0: orbital_list = list(orbital_dictionary.keys()) def get_orbitals(orbital): #case where spin polarization is false and m level is resloved if is_spin == False and m_projected == True: if orbital == 'p': orbitals = ['py', 'pz', 'px'] elif orbital == 'd': orbitals = ['dxy', 'dyz', 'dz2', 'dxz', 'dx2-y2'] elif orbital == 'f': orbitals = ['fy(3x2-y2)', 'fxyz', 'fyz2', 'fz3', 'fxz2', 'fz(x2-y2)', 'fx(x2-3y2)'] else: orbitals = [orbital] #case where spin polarization is true and m level is resloved elif is_spin == True and m_projected == True: #case where neither spin or m-level are provided if orbital == 's': orbitals = ['s+','s-'] elif orbital == 'p': orbitals = ['py+', 'py-', 'pz+', 'pz-', 'px+', 'px-'] elif orbital == 'd': orbitals = ['dxy+', 'dxy-', 'dyz+', 'dyz-', 'dz2+', 'dz2-', 'dxz+', 'dxz-', 'dx2-y2+', 'dx2-y2-'] elif orbital == 'f': orbitals = ['fy(3x2-y2)+', 'fy(3x2-y2)-','fxyz+', 'fxyz-', 'fyz2+', 'fyz2-', 'fz3+', 'fz3-', 'fxz2+', 'fxz2-', 'fz(x2-y2)+', 'fz(x2-y2)-', 'fx(x2-3y2)+', 'fx(x2-3y2)-'] #case where spin is provided but m-level is not #up spin elif orbital == 's+': orbitals = ['s+'] elif orbital == 's-': orbitals = ['s-'] elif orbital == 'p+': orbitals = ['py+', 'pz+', 'px+'] elif orbital == 'd+': orbitals = ['dxy+', 'dyz+', 'dz2+', 'dxz+', 'dx2-y2+'] elif orbital == 'f+': orbitals = ['fy(3x2-y2)+', 'fxyz+', 'fyz2+', 'fz3+', 'fxz2+', 'fz(x2-y2)+', 'fx(x2-3y2)+'] # down spin elif orbital == 'p-': orbitals = ['py-', 'pz-', 'px-'] elif orbital == 'd-': orbitals = ['dxy-', 'dyz-', 'dz2-', 'dxz-', 'dx2-y2-'] elif orbital == 'f-': orbitals = ['fy(3x2-y2)-', 'fxyz-', 'fyz2-', 'fz3-', 'fxz2-', 'fz(x2-y2)-', 'fx(x2-3y2)-'] elif '+' not in orbital and '-' not in orbital: orbitals = [orbital + '+', orbital + '-'] else: orbitals = [orbital] else: orbitals = [orbital] return orbitals #force list of atom_list is an int if sum_density == True: projected_density = np.zeros((len(orbital_list), ndos)) for atom in atom_indices: for count, orbital in enumerate(orbital_list): new_orbital_list = get_orbitals(orbital) indices = [orbital_dictionary[key] for key in new_orbital_list] projected_density[count]+= _site_dos[atom, indices,:].sum(axis=0) #return the list of orbitals provided new_orbital_list = orbital_list else: #get complete set of orbitals new_orbital_list = [] for orbital in orbital_list: new_orbital_list += get_orbitals(orbital) projected_density = np.zeros((len(new_orbital_list), ndos)) indices = [orbital_dictionary[key] for key in new_orbital_list] for atom in atom_indices: projected_density += _site_dos[atom, indices, :] if is_spin == True and sum_spin == True and sum_density == False: projected_density = projected_density[0::2,:] + projected_density[1::2,:] new_orbital_list = [key.rstrip('+') for key in new_orbital_list[::2]] return new_orbital_list, projected_density
[docs] def _read_doscar(self, file_name="DOSCAR", no_negatives=True): """Read VASP DOSCAR and extract projected densities Parameters ---------- file_name : str file location of the DOSCAR file no_negatives : bool Indicates wheather negative occupances will be converted to zero. Negative occupances can occur if Methfessel-Paxton is used. Attributes ---------- _total_dos : numpy.ndarray numpy array that contains the energy of the orbitals and the total projected and integrated density _site_dos : numpy.ndarray numpy array that contains the energy of the orbitals and the site and orbital projected density of states. Only available if a site projected calculation was performed. Returns ------- emax : float maximum energy level emin : float minimum energy level ndos : int number of descritized energy levels e_fermi : float highest occupied energy level is_spin : bool indicates if projected density is spin resolved m_projected : bool indicates if projected density is orbital resolved orbital_dictionary : dict dictionary that maps resolved sublevels/orbitals to indices """ #Accepts a file and reads through to get the density of states def get_dos(f, ndos): #get first line of density line = f.readline().split() dos = np.zeros((ndos, len(line))) dos[0] = np.array(line) for nd in range(1, ndos): line = f.readline().split() dos[nd] = np.array([float(x) for x in line]) return dos.T f = open(file_name) natoms = int(f.readline().split()[0]) [f.readline() for lines in range(4)] # Skip next 4 lines. # First we have a block with total and total integrated DOS descriptive_line = f.readline().split() emax = float(descriptive_line[0]) emin = float(descriptive_line[1]) ndos = int(descriptive_line[2]) e_fermi = float(descriptive_line[3]) _total_dos = get_dos(f,ndos) if no_negatives == True: _total_dos[1:][_total_dos[1:][...] < 0] = 0 # if energies have been shifted by fermi level undo that if np.abs(_total_dos[0].min() + e_fermi - emin)\ < np.abs(_total_dos[0].min() - emin): _total_dos[0] += e_fermi if _total_dos.shape[0] == 5: is_spin = True elif _total_dos.shape[0] == 3: is_spin = False # Next we have one block per atom, if DOSCAR contains the stuff # necessary for generating site-projected DOS dos = [] for na in range(natoms): line = f.readline() if line == '': # No site-projected DOS break #ndos = int(line.split()[2]) ndos does not change pdos = get_dos(f,ndos) dos.append(pdos) _site_dos = np.array(dos) if no_negatives == True and len(_site_dos.shape) > 1: _site_dos[:,1:,:][_site_dos[:,1:,:][...] < 0] = 0 # Integer indexing for orbitals starts from 1 in the _site_dos array # since the 0th column contains the energies if len(_site_dos.shape) > 1: norbs = _site_dos.shape[1] - 1 if norbs == 3: m_projected = False orbitals = {'s': 1, 'p': 2, 'd': 3} elif norbs == 4: m_projected = False orbitals = {'s': 1, 'p': 2, 'd': 3, 'f': 4} elif norbs == 6: m_projected = False orbitals = {'s+': 1, 's-': 2, 'p+': 3, 'p-': 4, 'd+': 5, 'd-': 6} elif norbs == 8: m_projected = False orbitals = {'s+': 1, 's-': 2, 'p+': 3, 'p-': 4, 'd+': 5, 'd-': 6, 'f+': 7, 'f-': 8} elif norbs == 9: m_projected = True orbitals = {'s': 1, 'py': 2, 'pz': 3, 'px': 4, 'dxy': 5, 'dyz': 6, 'dz2': 7, 'dxz': 8, 'dx2-y2': 9} elif norbs == 16: m_projected = True orbitals = {'s': 1, 'py': 2, 'pz': 3, 'px': 4, 'dxy': 5, 'dyz': 6, 'dz2': 7, 'dxz': 8, 'dx2': 9, 'fy(3x2-y2)': 10, 'fxyz': 11, 'fyz2': 12, 'fz3': 13, 'fxz2': 14, 'fz(x2-y2)': 15, 'fx(x2-3y2)': 16} elif norbs == 18: m_projected = True orbitals = {'s+': 1, 's-': 2, 'py+': 3, 'py-': 4, 'pz+': 5, 'pz-': 6, 'px+': 7, 'px-': 8, 'dxy+': 9, 'dxy-': 10, 'dyz+': 11, 'dyz-': 12, 'dz2+': 13, 'dz2-': 14, 'dxz+': 15, 'dxz-': 16, 'dx2-y2+': 17, 'dx2-y2-': 18} elif norbs == 32: m_projected = True orbitals = {'s+': 1, 's-': 2, 'py+': 3, 'py-': 4, 'pz+': 5, 'pz-': 6, 'px+': 7, 'px-': 8, 'dxy+': 9, 'dxy-': 10, 'dyz+': 11, 'dyz-': 12, 'dz2+': 13, 'dz2-': 14, 'dxz+': 15, 'dxz-': 16, 'dx2-y2+': 17, 'dx2-y2-': 18, 'fy(3x2-y2)+': 19, 'fy(3x2-y2)-': 20, 'fxyz+': 21, 'fxyz-': 22, 'fyz2+': 23, 'fyz2-': 24, 'fz3+': 25, 'fz3-': 26, 'fxz2+': 27, 'fxz2-': 28, 'fz(x2-y2)+': 29, 'fz(x2-y2)-': 30, 'fx(x2-3y2)+': 31, 'fx(x2-3y2)-': 32} else: norbs = None orbitals = None m_projected = False self._total_dos = _total_dos self._site_dos = _site_dos return natoms, emax, emin, ndos, e_fermi, is_spin, m_projected, orbitals
[docs] def _read_lobster(self, file_name="DOSCAR.lobster", no_negatives=True): """Read VASP DOSCAR and extract projected densities Parameters ---------- file_name : str file location of the DOSCAR file no_negatives : bool Indicates wheather negative occupances will be converted to zero. Negative occupances can occur if Methfessel-Paxton is used. Attributes ---------- _total_dos : numpy.ndarray numpy array that contains the energy of the orbitals and the total projected and integrated density _site_dos : numpy.ndarray numpy array that contains the energy of the orbitals and the site and orbital projected density of states. Only available if a site projected calculation was performed. Returns ------- emax : float maximum energy level emin : float minimum energy level ndos : int number of descritized energy levels e_fermi : float highest occupied energy level is_spin : bool indicates if projected density is spin resolved m_projected : bool indicates if projected density is orbital resolved orbital_dictionary : dict dictionary that maps resolved sublevels/orbitals to indices """ #Accepts a file and reads through to get the density of states def get_dos(f, ndos): #get first line of density line = f.readline().split() dos = np.zeros((ndos, len(line))) dos[0] = np.array(line) for nd in range(1, ndos): line = f.readline().split() dos[nd] = np.array([float(x) for x in line]) return dos.T f = open(file_name) natoms = int(f.readline().split()[0]) [f.readline() for lines in range(4)] # Skip next 4 lines. # First we have a block with total and total integrated DOS descriptive_line = f.readline().split() emax = float(descriptive_line[0]) emin = float(descriptive_line[1]) ndos = int(descriptive_line[2]) e_fermi = float(descriptive_line[3]) _total_dos = get_dos(f,ndos) if no_negatives == True: _total_dos[1:][_total_dos[1:][...] < 0] = 0 # if energies have been shifted by fermi level undo that if np.abs(_total_dos[0].min() + e_fermi - emin)\ < np.abs(_total_dos[0].min() - emin): _total_dos[0] += e_fermi if _total_dos.shape[0] == 5: is_spin = True elif _total_dos.shape[0] == 3: is_spin = False # Next we have one block per atom, if DOSCAR contains the stuff # necessary for generating site-projected DOS dos = [] for na in range(natoms): line = f.readline() if line == '': # No site-projected DOS break #ndos = int(line.split()[2]) ndos does not change pdos = get_dos(f,ndos) dos.append(pdos) # Integer indexing for orbitals starts from 1 in the _site_dos array # since the 0th column contains the energies if is_spin == False: orbitals = {'s': 1, 'py': 2, 'pz': 3, 'px': 4, 'dxy': 5, 'dyz': 6, 'dz2': 7, 'dxz': 8, 'dx2-y2': 9} _site_dos = np.zeros((natoms, 10, ndos)) # not every pdos will have the same orbitals (basis set) for count, pdos in enumerate(dos): if pdos.shape[0] == 2: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s elif pdos.shape[0] == 4: _site_dos[count][0] = pdos[0] # energy _site_dos[count][2] = pdos[1] # py _site_dos[count][3] = pdos[2] # pz _site_dos[count][4] = pdos[3] # px elif pdos.shape[0] == 5: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s _site_dos[count][2] = pdos[2] # py _site_dos[count][3] = pdos[3] # pz _site_dos[count][4] = pdos[4] # px elif pdos.shape[0] == 6: _site_dos[count][0] = pdos[0] # energy _site_dos[count][5] = pdos[1] # dxy _site_dos[count][6] = pdos[2] # dyz _site_dos[count][7] = pdos[3] # dz2 _site_dos[count][8] = pdos[4] # dxz _site_dos[count][9] = pdos[5] # dx2-y2 elif pdos.shape[0] == 7: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s _site_dos[count][5] = pdos[2] # dxy _site_dos[count][6] = pdos[3] # dyz _site_dos[count][7] = pdos[4] # dz2 _site_dos[count][8] = pdos[5] # dxz _site_dos[count][9] = pdos[6] # dx2-y2 elif pdos.shape[0] == 9: _site_dos[count][0] = pdos[0] # energy _site_dos[count][2] = pdos[1] # py _site_dos[count][3] = pdos[2] # pz _site_dos[count][4] = pdos[3] # px _site_dos[count][5] = pdos[4] # dxy _site_dos[count][6] = pdos[5] # dyz _site_dos[count][7] = pdos[6] # dz2 _site_dos[count][8] = pdos[7] # dxz _site_dos[count][9] = pdos[8] # dx2-y2 else: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s _site_dos[count][2] = pdos[2] # py _site_dos[count][3] = pdos[3] # pz _site_dos[count][4] = pdos[4] # px _site_dos[count][5] = pdos[5] # dxy _site_dos[count][6] = pdos[6] # dyz _site_dos[count][7] = pdos[7] # dz2 _site_dos[count][8] = pdos[8] # dxz _site_dos[count][9] = pdos[9] # dx2-y2 else: orbitals = {'s+': 1, 's-': 2, 'py+': 3, 'py-': 4, 'pz+': 5, 'pz-': 6, 'px+': 7, 'px-': 8, 'dxy+': 9, 'dxy-': 10, 'dyz+': 11, 'dyz-': 12, 'dz2+': 13, 'dz2-': 14, 'dxz+': 15, 'dxz-': 16, 'dx2-y2+': 17, 'dx2-y2-': 18} _site_dos = np.zeros((natoms, 19, ndos)) for count, pdos in enumerate(dos): if pdos.shape[0] == 3: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s+ _site_dos[count][2] = pdos[2] # s- elif pdos.shape[0] == 7: _site_dos[count][0] = pdos[0] # energy _site_dos[count][3] = pdos[1] # py+ _site_dos[count][4] = pdos[2] # py- _site_dos[count][5] = pdos[3] # pz+ _site_dos[count][6] = pdos[4] # pz- _site_dos[count][7] = pdos[5] # px+ _site_dos[count][8] = pdos[6] # px- elif pdos.shape[0] == 9: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s+ _site_dos[count][2] = pdos[2] # s- _site_dos[count][3] = pdos[3] # py+ _site_dos[count][4] = pdos[4] # py- _site_dos[count][5] = pdos[5] # pz+ _site_dos[count][6] = pdos[6] # pz- _site_dos[count][7] = pdos[7] # px+ _site_dos[count][8] = pdos[8] # px- elif pdos.shape[0] == 11: _site_dos[count][0] = pdos[0] # energy _site_dos[count][9] = pdos[1] # dxy+ _site_dos[count][10] = pdos[2] # dxy- _site_dos[count][11] = pdos[3] # dyz+ _site_dos[count][12] = pdos[4] # dyz- _site_dos[count][13] = pdos[5] # dz2+ _site_dos[count][14] = pdos[6] # dz2- _site_dos[count][15] = pdos[7] # dxz+ _site_dos[count][16] = pdos[8] # dxz- _site_dos[count][17] = pdos[9] # dx2-y2+ _site_dos[count][18] = pdos[10] # dx2-y2- elif pdos.shape[0] == 13: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s+ _site_dos[count][2] = pdos[2] # s- _site_dos[count][9] = pdos[3] # dxy+ _site_dos[count][10] = pdos[4] # dxy- _site_dos[count][11] = pdos[5] # dyz+ _site_dos[count][12] = pdos[6] # dyz- _site_dos[count][13] = pdos[7] # dz2+ _site_dos[count][14] = pdos[8] # dz2- _site_dos[count][15] = pdos[9] # dxz+ _site_dos[count][16] = pdos[10] # dxz- _site_dos[count][17] = pdos[11] # dx2-y2+ _site_dos[count][18] = pdos[12] # dx2-y2- elif pdos.shape[0] == 17: _site_dos[count][0] = pdos[0] # energy _site_dos[count][3] = pdos[1] # py+ _site_dos[count][4] = pdos[2] # py- _site_dos[count][5] = pdos[3] # pz+ _site_dos[count][6] = pdos[4] # pz- _site_dos[count][7] = pdos[5] # px+ _site_dos[count][8] = pdos[6] # px- _site_dos[count][9] = pdos[7] # dxy+ _site_dos[count][10] = pdos[8] # dxy- _site_dos[count][11] = pdos[9] # dyz+ _site_dos[count][12] = pdos[10] # dyz- _site_dos[count][13] = pdos[11] # dz2+ _site_dos[count][14] = pdos[12] # dz2- _site_dos[count][15] = pdos[13] # dxz+ _site_dos[count][16] = pdos[14] # dxz- _site_dos[count][17] = pdos[15] # dx2-y2+ _site_dos[count][18] = pdos[16] # dx2-y2- else: _site_dos[count][0] = pdos[0] # energy _site_dos[count][1] = pdos[1] # s+ _site_dos[count][2] = pdos[2] # s- _site_dos[count][3] = pdos[3] # py+ _site_dos[count][4] = pdos[4] # py- _site_dos[count][5] = pdos[5] # pz+ _site_dos[count][6] = pdos[6] # pz- _site_dos[count][7] = pdos[7] # px+ _site_dos[count][8] = pdos[8] # px- _site_dos[count][9] = pdos[9] # dxy+ _site_dos[count][10] = pdos[10] # dxy- _site_dos[count][11] = pdos[11] # dyz+ _site_dos[count][12] = pdos[12] # dyz- _site_dos[count][13] = pdos[13] # dz2+ _site_dos[count][14] = pdos[14] # dz2- _site_dos[count][15] = pdos[15] # dxz+ _site_dos[count][16] = pdos[16] # dxz- _site_dos[count][17] = pdos[17] # dx2-y2+ _site_dos[count][18] = pdos[18] # dx2-y2- m_projected = True if no_negatives == True and len(_site_dos.shape) > 1: _site_dos[:,1:,:][_site_dos[:,1:,:][...] < 0] = 0 self._total_dos = _total_dos self._site_dos = _site_dos return natoms, emax, emin, ndos, e_fermi, is_spin, m_projected, orbitals