Source code for pyprocar.lobsterparser.lobsterparser

# -*- coding: utf-8 -*-
"""
Created on Thu Jun 25 03:12:50 2020
@author: Logan Lang
"""


from re import findall, search, match, DOTALL, MULTILINE, finditer, compile, sub

from numpy import array, dot, linspace, sum, where, zeros, pi, concatenate
import logging


[docs]class LobsterParser: def __init__( self, lobsterin="lobsterin", lobsterout="lobsterout", scfin="scf.in", outfile="scf.out", kdirect=True, lobstercode="qe", ): self.lobsterin = lobsterin self.scfin = scfin self.outfile = outfile self.lobsterout = lobsterout self.dosin = "DOSCAR.lobster" self.file_names = [] self.discontinuities = [] self.kdirect = kdirect ############################ self.high_symmetry_points = None self.nhigh_sym = None self.nspecies = None self.composition = {} self.kticks = None self.knames = None self.speciesList = None self.kpoints = None self.kpointsCount = None self.bands = None self.bandsCount = None self.ionsList = None self.ionsCount = None self.spinCount = None # Used to store atomic states at the top of the kpdos.out file self.states = None self.spd = None """ for l=1: 1 pz (m=0) 2 px (real combination of m=+/-1 with cosine) 3 py (real combination of m=+/-1 with sine) for l=2: 1 dz2 (m=0) 2 dzx (real combination of m=+/-1 with cosine) 3 dzy (real combination of m=+/-1 with sine) 4 dx2-y2 (real combination of m=+/-2 with cosine) 5 dxy (real combination of m=+/-2 with sine)""" # Oribital order. This is the way it goes into the spd Array. index 0 is resered for totals self.orbitals = [ "s", "p_y", "p_z", "p_x", "d_xy", "d_yz", "d_z^2", "d_x^2-y^2", "d_xz", ] # These are not used self.orbitalName_short = ["s", "p", "d", "tot"] self.orbitalCount = None # Opens the needed files rf = open(self.lobsterin, "r") self.lobsterin = rf.read() rf.close() rf = open(self.lobsterout, "r") self.lobsterout = rf.read() rf.close() rf = open(self.scfin, "r") self.scfin = rf.read() rf.close() self.dos = None self._readFileNames() self._readProjFiles() self._readHighKpoint() return ########################################################################################## # Finding high symmetry points ########################################################################################## def _readHighKpoint(self): numK = int(findall("K_POINTS.*\n([0-9]*)", self.scfin)[0]) raw_khigh_sym = findall( "K_POINTS.*\n\s*[0-9]*.*\n" + numK * "(.*)\n", self.scfin, )[0] tickCountIndex = 0 raw_high_symmetry = [] self.knames = [] self.kticks = [] for x in raw_khigh_sym: if len(x.split()) == 5: raw_high_symmetry.append( (float(x.split()[0]), float(x.split()[1]), float(x.split()[2])) ) self.knames.append("$%s$" % x.split()[4].replace("!", "")) self.kticks.append(tickCountIndex) if float(x.split()[3]) == 0: tickCountIndex += 1 self.high_symmetry_points = array(raw_high_symmetry) self.nhigh_sym = len(self.knames) # finds discontinuities for i in range(len(self.kticks)): if(i < len(self.kticks)-1): diff = self.kticks[ i+1 ] - self.kticks[i] if diff == 1 : self.discontinuities.append(self.kticks[i]) discon_name = "$" + self.knames[i].replace("$","") +"|"+ self.knames[i+1].replace("$","") + "$" self.knames.pop(i+1) self.knames[i] = discon_name self.kticks.pop(i+1) ########################################################################################## # Finding file names ########################################################################################## def _readFileNames(self): self.ionsList = findall( "calculating FatBand for Element: (.*) Orbital.*", self.lobsterout ) raw_speciesList = findall( "calculating FatBand for Element: ([a-zA-Z]*)[0-9] Orbital.*", self.lobsterout, ) self.speciesList = [] for x in raw_speciesList: if x not in self.speciesList: self.speciesList.append(x) self.composition.update({str(x): 0}) self.nspecies = len(self.speciesList) for x in raw_speciesList: if x in self.speciesList: self.composition[x] += 1 self.ionsCount = len(self.ionsList) proj_orbitals = findall( "calculating FatBand for Element: (.*) Orbital\(s\):\s*(.*)", self.lobsterout, ) for proj in proj_orbitals: orbitals = proj[1].split() for orbital in orbitals: fileName = "FATBAND_" + proj[0] + "_" + orbital + ".lobster" self.file_names.append(fileName) ########################################################################################## # Reading projection files ########################################################################################## def _readProjFiles(self): rf = open(self.file_names[0], "r") projFile = rf.read() rf.close() ########################################################################################## # kpoints ########################################################################################## raw_kpoints = findall( "# K-Point \d+ :\s*([-\.\\d]*)\s*([-\.\\d]*)\s*([-\.\\d]*)", projFile ) self.kpointsCount = len(raw_kpoints) self.kpoints = zeros(shape=(self.kpointsCount, 3)) for ik in range(len(raw_kpoints)): for coord in range(3): self.kpoints[ik][coord] = raw_kpoints[ik][coord] # raw_kpoints = findall("(# K-Point \d+ :\s*[-\.\\d]*\s*[-\.\\d]*\s*[-\.\\d]*)", projFile) # If kdirect=False, then the kgrid will be in cartesian coordinates. # Requires the reciprocal lattice vectors to be parsed from the output. if not self.kdirect: self.kpoints = dot(self.kpoints, self.reclat) ######################################################################################### # bands ######################################################################################### raw_kpoints = findall( "# K-Point \d+ :\s*([-\.\\d]*\s*[-\.\\d]*\s*[-\.\\d]*)", projFile ) # Checks for spin polarization if len(findall("spillings for spin channel", self.lobsterout)) == 2: self.spinCount = 2 bandsCount = int(findall("NBANDS\s*(\d*)", projFile)[0]) * self.spinCount self.bandsCount = bandsCount // 2 else: self.spinCount = 1 bandsCount = int(findall("NBANDS\s*(\d*)", projFile)[0]) self.bandsCount = bandsCount band_info = [] for ik in range(len(raw_kpoints)): expression = "# K-Point \d+ :\s*" + raw_kpoints[ik] + ".*\n" kpoint_bands = findall(expression + bandsCount * "(.*)\n", projFile)[0] for ikband in kpoint_bands: band_info.append(ikband) # if len(band_info) == self.bandsCount * self.kpointsCount: # print("Number of bands headers match") raw_bands = zeros(shape=(self.kpointsCount, bandsCount)) ik = 0 ib = 0 for i in range(len(band_info)): raw_bands[ik, ib] = float(band_info[i].split()[1]) ib += 1 if int(ib == bandsCount): ik += 1 ib = 0 # Checks for spin polarization if self.spinCount == 2: self.bands = zeros( shape=(self.kpointsCount, self.bandsCount, self.spinCount) ) self.bands[:, 0 : self.bandsCount, 0] = raw_bands[:, 0 : self.bandsCount] self.bands[:, 0 : self.bandsCount, 1] = raw_bands[ :, self.bandsCount : self.bandsCount * 2 ] else: self.bands = zeros( shape=(self.kpointsCount, self.bandsCount, self.spinCount) ) self.bands[:, 0 : self.bandsCount, 0] = raw_bands[:, 0 : self.bandsCount] ######################################################################################### # Forming SPD array ######################################################################################### # Checks for spin polarization if len(findall("spillings for spin channel", self.lobsterout)) == 2: self.spinCount = 2 else: self.spinCount = 1 self.orbitalCount = 10 self.spd = zeros( shape=( self.kpointsCount, self.bandsCount, self.spinCount, self.ionsCount + 1, len(self.orbitals) + 2, ) ) for file in range(len(self.file_names)): rf = open(self.file_names[file], "r") projFile = rf.read() rf.close() raw_kpoints = findall( "# K-Point \d+ :\s*([-\.\\d]*\s*[-\.\\d]*\s*[-\.\\d]*)", projFile ) ionIndex = 0 orbitalIndex = 0 current_orbital = findall("(# FATBAND.*)", projFile)[0].split()[4] for i in range(len(self.ionsList)): if self.ionsList[i] == self.file_names[file].split("_")[1]: ionIndex = i # print(ionIndex) for i in range(len(self.orbitals)): if self.orbitals[i] == sub("[0-9]", "", current_orbital): orbitalIndex = i + 1 # print(orbitalIndex) # print(self.orbitals[orbitalIndex]) band_info = [] for ik in range(len(raw_kpoints)): expression = "# K-Point \d+ :\s*" + raw_kpoints[ik] + ".*\n" bands_wSpin = self.bandsCount * self.spinCount band = findall(expression + bands_wSpin * "(.*)\n", projFile)[0] for iband in range(self.bandsCount): if self.spinCount == 2: self.spd[ik, iband, 0, ionIndex, orbitalIndex] += float( band[iband].split()[2] ) self.spd[ik, iband, 0, ionIndex, 0] = ionIndex + 1 self.spd[ik, iband, 1, ionIndex, orbitalIndex] += float( band[iband + self.bandsCount].split()[2] ) self.spd[ik, iband, 1, ionIndex, 0] = ionIndex + 1 else: self.spd[ik, iband, 0, ionIndex, orbitalIndex] += float( band[iband].split()[2] ) self.spd[ik, iband, 0, ionIndex, 0] = ionIndex + 1 self.spd[:, :, :, :, -1] = sum(self.spd[:, :, :, :, 1:-1], axis=4) self.spd[:, :, :, -1, :] = sum(self.spd[:, :, :, 0:-1, :], axis=3) self.spd[:, :, :, -1, 0] = 0 # colinear spin polarized case # manipulating spd array for spin polarized calculations. # The shape is (nkpoints,2*nbands,2,natoms,norbitals) # The third dimension is for spin. # When this is zero, the bands*2 (all spin up and down bands) have positive projections. # When this is one, the the first half of bands (spin up) will have positive projections # and the second half (spin down) will have negative projections. This is to adhere to # the convention used in PyProcar to obtain spin density and spin magnetization. if self.spinCount == 2: print("\nLobster colinear spin calculation detected.\n") self.spd2 = zeros( shape=( self.kpointsCount, self.bandsCount * 2, self.spinCount, self.ionsCount + 1, len(self.orbitals) + 2, ) ) # spin up block for spin=0 self.spd2[:, : self.bandsCount, 0, :, :] = self.spd[:, :, 0, :, :] # spin down block for spin=0 self.spd2[:, self.bandsCount :, 0, :, :] = self.spd[:, :, 1, :, :] # spin up block for spin=1 self.spd2[:, : self.bandsCount, 1, :, :] = self.spd[:, :, 0, :, :] # spin down block for spin=1 self.spd2[:, self.bandsCount :, 1, :, :] = -1 * self.spd[:, :, 1, :, :] self.spd = self.spd2 # Reshaping bands array to inlcude all bands (spin up and down) # self.bands = concatenate((self.bands, self.bands), axis=1) # Reshaping bands array to inlcude all bands (spin up and down) self.bands = self.bands.reshape( self.kpointsCount, self.bandsCount * 2, order="F" ) else: self.bands = self.bands.reshape( self.kpointsCount, self.bandsCount, order="F" ) @property def fermi(self): """ Returns the fermi energy read from .out """ fi = open(self.outfile, "r") data = fi.read() fi.close() fermi = float(findall(r"the\s*Fermi\s*energy\s*is\s*([\s\d.]*)ev", data)[0]) return fermi @property def reclat(self): """ Returns the reciprocal lattice read from .out """ rf = open(self.outfile, "r") data = rf.read() rf.close() alat = float(findall(r"alat\)\s*=\s*([\d.e+-]*)", data)[0]) b1 = array( findall(r"b\(1\)\s*=\s*\(([\d.\s+-e]*)", data)[0].split(), dtype="float64" ) b2 = array( findall(r"b\(2\)\s*=\s*\(([\d.\s+-e]*)", data)[0].split(), dtype="float64" ) b3 = array( findall(r"b\(3\)\s*=\s*\(([\d.\s+-e]*)", data)[0].split(), dtype="float64" ) reclat = (2 * pi / alat) * (array((b1, b2, b3))) # Transposing to get the correct format reclat = reclat.T return reclat