Source code for pyprocar.elkparser.elkparser

"""
Created on Fri March 2 2020
@author: Pedram Tavadze

"""

from re import findall

import numpy as np


[docs]class ElkParser: def __init__(self, elkin="elk.in", kdirect=True): # elk specific inp self.fin = elkin self.file_names = [] self.elkin = None self.nbands = None self.nspin = None self.kticks = None self.kdirect = kdirect self.kpoints = None self.bands = None self.spd = None self.cspd = None # spin polarized parameters self.spinpol = None self.orbitalName = [ "Y00", "Y1-1", "Y10", "Y11", "Y2-2", "Y2-1", "Y20", "Y21", "Y22", "Y3-3", "Y3-2", "Y3-1", "Y30", "Y3-1", "Y3-2", "Y3-3", ] self.orbital_name_short = ["s", "p", "d", "f"] self.norbital = None # number of orbitals # number of spin components (blocks of data), 1: non-magnetic non # polarized, 2: spin polarized collinear, 4: non-collinear # spin. # NOTE: before calling to `self._readOrbital` the case '4' # is marked as '1' self._read_elkin() self._read_bands() return # @property # def nspin(self): # """ # number of spin, default is 1. # """ # nspindict = {1: 1, 2: 2, 4: 2, None: 1} # return nspindict[self.ispin] @property def spd_orb(self): # indices: ikpt, iband, ispin, iion, iorb # remove indices and total from iorb. return self.spd[:, :, :, 1:-1] @property def nhigh_sym(self): """ Returns nu,ber of high symmetry points """ return int(findall("plot1d\n\s*([0-9]*)\s*[0-9]*", self.elkin)[0]) @property def nkpoints(self): """ Returns total number of kpoints """ return int(findall("plot1d\n\s*[0-9]*\s*([0-9]*)", self.elkin)[0]) @property def tasks(self): """ Returns the tasks calculated by elk """ return [int(x) for x in findall("tasks\s*([0-9\s\n]*)", self.elkin)[0].split()] @property def high_symmetry_points(self): """ Returns the corrdinates of high symmtery points provided in elk.in """ if 20 in self.tasks or 21 in self.tasks or 22 in self.tasks: raw_hsp = findall( "plot1d\n\s*[0-9]*\s*[0-9]*.*\n" + self.nhigh_sym * "([0-9\s\.-]*).*\n", self.elkin, )[0] high_symmetry_points = np.zeros(shape=(self.nhigh_sym, 3)) for ihs in range(self.nhigh_sym): high_symmetry_points[ihs, :] = [ float(x) for x in raw_hsp[ihs].split()[0:3] ] return high_symmetry_points @property def nspecies(self): """ Returns number of species """ return int(findall("atoms\n\s*([0-9]*)", self.elkin)[0]) @property def natom(self): """ Returns number of atoms in the structure """ natom = 0 for ispc in findall("'([A-Za-z]*).in'.*\n\s*([0-9]*)", self.elkin): natom += int(ispc[1]) return natom @property def composition(self): """ Returns the composition of the structure """ composition = {} for ispc in findall("'([A-Za-z]*).in'.*\n\s*([0-9]*)", self.elkin): composition[ispc[0]] = int(ispc[1]) return composition @property def knames(self): """ Returns the names of the high symmetry points(x ticks) provided in elk.in """ raw_ticks = findall( "plot1d\n\s*[0-9]*\s*[0-9]*.*\n" + self.nhigh_sym * ".*:(.*)\n", self.elkin )[0] if len(raw_ticks) != self.nhigh_sym: knames = [str(x) for x in range(self.nhigh_sym)] else: knames = [ "$%s$" % (x.replace(",", "").replace("vlvp1d", "").replace(" ", "")) for x in findall( "plot1d\n\s*[0-9]*\s*[0-9]*.*\n" + self.nhigh_sym * ".*:(.*)\n", self.elkin, )[0] ] return knames def _read_elkin(self): """ Reads and parses elk.in """ rf = open(self.fin, "r") self.elkin = rf.read() rf.close() if 20 in self.tasks: self.file_names.append("BANDS.OUT") if 21 in self.tasks or 22 in self.tasks: ispc = 1 for spc in self.composition: for iatom in range(self.composition[spc]): self.file_names.append( "BAND_S{:02d}_A{:04d}.OUT".format(ispc, iatom + 1) ) ispc += 1 # Checking if spinpol = .true. in elk.in try: self.spinpol = findall(r"spinpol\s*([.a-zA-Z]*)", self.elkin)[0] except IndexError: self.spinpol = None if self.spinpol: if self.spinpol == ".true.": print("\nElk colinear spin calculation detected.\n") self.spinpol = True self.nspin = 2 else: print("\nElk non spin calculation detected.\n") self.nspin = 1 else: print( "\nNo spinpol keyword found in elk.in. Assuming non spin calculation.\n" ) self.nspin = 1 def _read_bands(self): """ if the task is any of 20,21,22 it parses the BANDS*.OUT files and prepares the spd, bands and kpoints for bandsplot """ if np.any([x in self.tasks for x in [20, 21, 22]]): rf = open(self.file_names[0], "r") lines = rf.readlines() rf.close() self.nbands = int(len(lines) / (self.nkpoints + 1)) self.bands = np.zeros(shape=(self.nkpoints, self.nbands)) rf = open("BANDLINES.OUT", "r") bandLines = rf.readlines() rf.close() tick_pos = [] # using strings for a better comparision and avoiding rounding by python for iline in range(0, len(bandLines), 3): tick_pos.append(bandLines[iline].split()[0]) x_points = [] for iline in range(self.nkpoints): x_points.append(lines[iline].split()[0]) self.kpoints = np.zeros(shape=(self.nkpoints, 3)) x_points = np.array(x_points) tick_pos = np.array(tick_pos) self.kticks = [] for ihs in range(1, self.nhigh_sym): start = np.where(x_points == tick_pos[ihs - 1])[0][0] end = np.where(x_points == tick_pos[ihs])[0][0] + 1 self.kpoints[start:end][:] = np.linspace( self.high_symmetry_points[ihs - 1], self.high_symmetry_points[ihs], end - start, ) self.kticks.append(start) self.kticks.append(self.nkpoints - 1) if not self.kdirect: self.kpoints = np.dot(self.kpoints, self.reclat) rf = open(self.file_names[0], "r") lines = rf.readlines() rf.close() iline = 0 for iband in range(self.nbands): for ikpoint in range(self.nkpoints): self.bands[ikpoint, iband] = float(lines[iline].split()[1]) iline += 1 if ikpoint == self.nkpoints - 1: iline += 1 self.bands *= 27.21138386 self.norbital = 16 self.spd = np.zeros( shape=( self.nkpoints, self.nbands, self.nspin, self.natom + 1, self.norbital + 2, ) ) idx_bands_out = None for ifile in range(len(self.file_names)): if self.file_names[ifile] == "BANDS.OUT": idx_bands_out = ifile if idx_bands_out != None: del self.file_names[idx_bands_out] for ifile in range(self.natom): rf = open(self.file_names[ifile], "r") lines = rf.readlines() rf.close() iline = 0 for iband in range(self.nbands): for ikpoint in range(self.nkpoints): temp = np.array([float(x) for x in lines[iline].split()]) self.spd[ikpoint, iband, 0, ifile, 0] = ifile + 1 self.spd[ikpoint, iband, 0, ifile, 1:-1] = temp[2:] iline += 1 if ikpoint == self.nkpoints - 1: iline += 1 # self.spd[:,:,:,-1,:] = self.spd.sum(axis=3) self.spd[:, :, :, :, -1] = np.sum(self.spd[:, :, :, :, 1:-1], axis=4) self.spd[:, :, :, -1, :] = self.spd.sum(axis=3) self.spd[:, :, 0, -1, 0] = 0 if self.nspin == 2: # spin up block for spin = 1 self.spd[:, : self.nbands // 2, 1, :, :] = self.spd[ :, : self.nbands // 2, 0, :, : ] # spin down block for spin = 1 self.spd[:, self.nbands // 2 :, 1, :, :] = ( -1 * self.spd[:, self.nbands // 2 :, 0, :, :] ) # 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. # spin up and spin down block for spin = 0 # spd2[:, :, 0, :, :] = self.spd[:, :, 0, :, :] @property def fermi(self): """ Returns the fermi energy read from FERMI.OUT """ rf = open("EFERMI.OUT", "r") fermi = float(rf.readline().split()[0]) rf.close() return fermi @property def reclat(self): """ Returns the reciprocal lattice read from LATTICE.OUT """ rf = open("LATTICE.OUT", "r") data = rf.read() rf.close() lattice_block = findall(r"matrix\s*:([\s0-9.]*)Inverse", data) lattice_array = np.array(lattice_block[1].split(), dtype=float) reclat = lattice_array.reshape((3, 3)) return reclat