################################################################################ # # 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 . # ################################################################################ r""" vasp.vaspio =========== Input of required VASP data. Six VASP files are required: - PROJCAR - LOCPROJ - POSCAR - IBZKPT - EIGENVAL - DOSCAR """ import numpy as np import re #import plocar_io.c_plocar_io as c_plocar_io def read_lines(filename): r""" Generator of lines for a file Parameters ---------- filename (str) : name of the file """ with open(filename, 'r') as f: for line in f: yield line ################################################################################ ################################################################################ # # class VaspData # ################################################################################ ################################################################################ class VaspData: """ Container class for all VASP data. """ def __init__(self, vasp_dir, read_all=True, efermi_required=True): self.vasp_dir = vasp_dir self.plocar = Plocar() self.poscar = Poscar() self.kpoints = Kpoints() self.eigenval = Eigenval() self.doscar = Doscar() if read_all: self.plocar.from_file(vasp_dir) self.poscar.from_file(vasp_dir) self.kpoints.from_file(vasp_dir) try: self.eigenval.from_file(vasp_dir) except (IOError, StopIteration): self.eigenval.eigs = None self.eigenval.ferw = None print "!!! WARNING !!!: Error reading from EIGENVAL, trying LOCPROJ" try: self.doscar.from_file(vasp_dir) except (IOError, StopIteration): if efermi_required: # raise Exception("Efermi cannot be read from DOSCAR") pass else: # TODO: This a hack. Find out a way to determine ncdij without DOSCAR print "!!! WARNING !!!: Error reading from DOSCAR, taking Efermi from config" self.doscar.ncdij = self.plocar.nspin ################################################################################ ################################################################################ # # class Plocar # ################################################################################ ################################################################################ class Plocar: r""" Class containing raw PLO data from VASP. Properties ---------- - *plo* (numpy.array((nion, ns, nk, nb, nlmmax))) : raw projectors - *params* (dict) : parameters read from PLOCAR - *ferw* (array(nion, ns, nk, nb)) : Fermi weights from VASP """ def from_file(self, vasp_dir='./', plocar_filename='PLOCAR'): r""" Reads non-normalized projectors from a binary file (`PLOCAR' by default) generated by VASP PLO interface. Parameters ---------- vasp_dir (str) : path to the VASP working directory [default = `./'] plocar_filename (str) : filename [default = `PLOCAR'] """ # Add a slash to the path name if necessary if vasp_dir[-1] != '/': vasp_dir += '/' # self.params, self.plo, self.ferw = c_plocar_io.read_plocar(vasp_dir + plocar_filename) # self.proj_params, self.plo = self.temp_parser(projcar_filename=vasp_dir + "PROJCAR", locproj_filename=vasp_dir + "LOCPROJ") self.proj_params, self.plo = self.locproj_parser(locproj_filename=vasp_dir + "LOCPROJ") def temp_parser(self, projcar_filename='PROJCAR', locproj_filename='LOCPROJ'): r""" Parses PROJCAR (and partially LOCPROJ) to get VASP projectors. This is a prototype parser that should eventually be written in C for better performance on large files. Returns projector parameters (site/orbital indices etc.) and an array with projectors. """ orb_labels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2", "fz3", "fxz2", "fyz2", "fz(x2-y2)", "fxyz", "fx(x2-3y2)", "fy(3x2-y2)"] def lm_to_l_m(lm): l = int(np.sqrt(lm)) m = lm - l*l return l, m # Read the first line of LOCPROJ to get the dimensions with open(locproj_filename, 'rt') as f: line = f.readline() nproj, nspin, nk, nband = map(int, line.split()) plo = np.zeros((nproj, nspin, nk, nband), dtype=np.complex128) proj_params = [{} for i in xrange(nproj)] iproj_site = 0 is_first_read = True with open(projcar_filename, 'rt') as f: line = self.search_for(f, "^ *ISITE") while line: isite = int(line.split()[1]) if not is_first_read: for il in xrange(norb): ip_new = iproj_site * norb + il ip_prev = (iproj_site - 1) * norb + il proj_params[ip_new]['label'] = proj_params[ip_prev]['label'] proj_params[ip_new]['isite'] = isite proj_params[ip_new]['l'] = proj_params[ip_prev]['l'] proj_params[ip_new]['m'] = proj_params[ip_prev]['m'] for ispin in xrange(nspin): for ik in xrange(nk): # Parse the orbital labels and convert them to l,m-indices line = self.search_for(f, "^ *band") if is_first_read: cpatt = re.compile("lm= *([^\s]+)") labels = re.findall(cpatt, line) norb = len(labels) for il, label in enumerate(labels): lm = orb_labels.index(label) l, m = lm_to_l_m(lm) # For the first read 'iproj_site = 0' and only orbital index 'il' is used proj_params[il]['label'] = label proj_params[il]['isite'] = isite proj_params[il]['l'] = l proj_params[il]['m'] = m is_first_read = False # Read the block of nk * ns * nband complex numbers for ib in xrange(nband): line = f.readline() rtmp = map(float, line.split()[1:]) for il in xrange(norb): ctmp = complex(rtmp[2 * il], rtmp[2 * il + 1]) plo[iproj_site * norb + il, ispin, ik, ib] = ctmp # End of site-block iproj_site += 1 line = self.search_for(f, "^ *ISITE") print "Read parameters:" for il, par in enumerate(proj_params): print il, " -> ", par return proj_params, plo def locproj_parser(self, locproj_filename='LOCPROJ'): r""" Parses LOCPROJ (for VASP >= 5.4.2) to get VASP projectors. This is a prototype parser that should eventually be written in C for better performance on large files. Returns projector parameters (site/orbital indices etc.) and an array with projectors. """ orb_labels = ["s", "py", "pz", "px", "dxy", "dyz", "dz2", "dxz", "dx2-y2", "fy(3x2-y2)", "fxyz", "fyz2", "fz3", "fxz2", "fz(x2-y2)", "fx(x2-3y2)"] def lm_to_l_m(lm): l = int(np.sqrt(lm)) m = lm - l*l return l, m # Read the first line of LOCPROJ to get the dimensions with open(locproj_filename, 'rt') as f: line = f.readline() line = line.split("#")[0] sline = line.split() self.ncdij, nk, self.nband, nproj = map(int, sline[:4]) self.nspin = 1 if self.ncdij == 1 else 2 self.nspin_band = 2 if self.ncdij == 2 else 1 self.efermi = float(sline[4]) plo = np.zeros((nproj, self.nspin, nk, self.nband), dtype=np.complex128) proj_params = [{} for i in xrange(nproj)] iproj_site = 0 is_first_read = True # First read the header block with orbital labels line = self.search_for(f, "^ *ISITE") ip = 0 while line: sline = line.split(':') isite = int(sline[1].split()[0]) label = sline[-1].strip() lm = orb_labels.index(label) l, m = lm_to_l_m(lm) # ip_new = iproj_site * norb + il # ip_prev = (iproj_site - 1) * norb + il proj_params[ip]['label'] = label proj_params[ip]['isite'] = isite proj_params[ip]['l'] = l proj_params[ip]['m'] = m ip += 1 line = f.readline().strip() assert ip == nproj, "Number of projectors in the header is wrong in LOCPROJ" self.eigs = np.zeros((nk, self.nband, self.nspin_band)) self.ferw = np.zeros((nk, self.nband, self.nspin_band)) patt = re.compile("^orbital") # FIXME: fix spin indices for NCDIJ = 4 (non-collinear) assert self.ncdij < 4, "Non-collinear case is not implemented" for ispin in xrange(self.nspin): for ik in xrange(nk): for ib in xrange(self.nband): line = "" while not line: line = f.readline().strip() sline = line.split() isp_, ik_, ib_ = map(int, sline[1:4]) assert isp_ == ispin + 1 and ik_ == ik + 1 and ib_ == ib + 1, "Inconsistency in reading LOCPROJ" self.eigs[ik, ib, ispin] = float(sline[4]) self.ferw[ik, ib, ispin] = float(sline[5]) for ip in xrange(nproj): line = f.readline() sline = line.split() ctmp = complex(float(sline[1]), float(sline[2])) plo[ip, ispin, ik, ib] = ctmp print "Read parameters:" for il, par in enumerate(proj_params): print il, " -> ", par return proj_params, plo def search_for(self, f, patt): r""" Reads file 'f' until pattern 'patt' is encountered and returns the corresponding line. """ cpatt = re.compile(patt) line = "x" while not re.match(cpatt, line) and line: line = f.readline() return line ################################################################################ ################################################################################ # # class Poscar # ################################################################################ ################################################################################ class Poscar: """ Class containing POSCAR data from VASP. Properties ---------- nq (int) : total number of ions ntypes ([int]) : number of ion types nions (int) : a list of number of ions of each type a_brav (numpy.array((3, 3), dtype=float)) : lattice vectors q_types ([numpy.array((nions, 3), dtype=float)]) : a list of arrays each containing fractional coordinates of ions of a given type """ def __init__(self): self.q_cart = None self.b_rec = None def from_file(self, vasp_dir='./', poscar_filename='POSCAR'): """ Reads POSCAR and returns a dictionary. Parameters ---------- vasp_dir (str) : path to the VASP working directory [default = `./'] plocar_filename (str) : filename [default = `PLOCAR'] """ # Convenince local function def readline_remove_comments(): return f.next().split('!')[0].strip() # Add a slash to the path name if necessary if vasp_dir[-1] != '/': vasp_dir += '/' f = read_lines(vasp_dir + poscar_filename) # Comment line comment = f.next().rstrip() print " Found POSCAR, title line: %s"%(comment) # Read scale sline = readline_remove_comments() ascale = float(sline) # Read lattice vectors self.a_brav = np.zeros((3, 3)) for ia in xrange(3): sline = readline_remove_comments() self.a_brav[ia, :] = map(float, sline.split()) # Negative scale means that it is a volume scale if ascale < 0: vscale = -ascale vol = np.linalg.det(self.a_brav) ascale = (vscale / vol)**(1.0/3) self.a_brav *= ascale # Depending on the version of VASP there could be # an extra line with element names sline = readline_remove_comments() try: # Old v4.6 format: no element names self.nions = map(int, sline.split()) self.el_names = ['El%i'%(i) for i in xrange(len(self.nions))] except ValueError: # New v5.x format: read element names first self.el_names = sline.split() sline = readline_remove_comments() self.nions = map(int, sline.split()) # Set the number of atom sorts (types) and the total # number of atoms in the unit cell self.ntypes = len(self.nions) self.nq = sum(self.nions) # Check for the line 'Selective dynamics' (and ignore it) sline = readline_remove_comments() if sline[0].lower() == 's': sline = readline_remove_comments() # Check whether coordinates are cartesian or fractional cartesian = (sline[0].lower() in 'ck') if cartesian: brec = np.linalg.inv(self.a_brav.T) # Read atomic positions self.q_types = [] self.type_of_ion = [] for it in xrange(self.ntypes): # Array mapping ion index to type self.type_of_ion += self.nions[it] * [it] q_at_it = np.zeros((self.nions[it], 3)) for iq in xrange(self.nions[it]): sline = readline_remove_comments() qcoord = map(float, sline.split()[:3]) if cartesian: qcoord = np.dot(brec, qcoord) q_at_it[iq, :] = qcoord self.q_types.append(q_at_it) print " Total number of ions:", self.nq print " Number of types:", self.ntypes print " Number of ions for each type:", self.nions # print # print " Coords:" # for it in xrange(ntypes): # print " Element:", el_names[it] # print q_at[it] ################################################################################ ################################################################################ # # class Kpoints # ################################################################################ ################################################################################ class Kpoints: """ Class describing k-points and optionally tetrahedra. Properties ---------- - nktot (int) : total number of k-points in the IBZ - kpts (numpy.array((nktot, 3), dtype=float)) : k-point vectors (fractional coordinates) - ntet (int) : total number of k-point tetrahedra - itet (numpy.array((ntet, 5), dtype=float) : array of tetrahedra - volt (float) : volume of a tetrahedron (the k-grid is assumed to be uniform) """ # # Reads IBZKPT file # def from_file(self, vasp_dir='./', ibz_filename='IBZKPT'): """ Reads from IBZKPT: k-points and optionally tetrahedra topology (if present). Parameters ---------- vasp_dir (str) : path to the VASP working directory [default = `./'] plocar_filename (str) : filename [default = `PLOCAR'] """ # Add a slash to the path name if necessary if vasp_dir[-1] != '/': vasp_dir += '/' ibz_file = read_lines(vasp_dir + ibz_filename) # Skip comment line line = ibz_file.next() # Number of k-points line = ibz_file.next() self.nktot = int(line.strip().split()[0]) print print " {0:>26} {1:d}".format("Total number of k-points:", self.nktot) self.kpts = np.zeros((self.nktot, 3)) self.kwghts = np.zeros((self.nktot)) # Skip comment line line = ibz_file.next() for ik in xrange(self.nktot): line = ibz_file.next() sline = line.strip().split() self.kpts[ik, :] = map(float, sline[:3]) self.kwghts[ik] = float(sline[3]) self.kwghts /= self.nktot # Attempt to read tetrahedra # Skip comment line ("Tetrahedra") try: line = ibz_file.next() # Number of tetrahedra and volume = 1/(6*nkx*nky*nkz) line = ibz_file.next() sline = line.split() self.ntet = int(sline[0]) self.volt = float(sline[1]) print " {0:>26} {1:d}".format("Total number of tetrahedra:", self.ntet) # Traditionally, itet[it, 0] contains multiplicity self.itet = np.zeros((self.ntet, 5), dtype=int) for it in xrange(self.ntet): line = ibz_file.next() self.itet[it, :] = map(int, line.split()[:5]) except StopIteration, ValueError: print " No tetrahedron data found in %s. Skipping..."%(ibz_filename) self.ntet = 0 # data = { 'nktot': nktot, # 'kpts': kpts, # 'ntet': ntet, # 'itet': itet, # 'volt': volt } # # return data ################################################################################ ################################################################################ # # class Eigenval # ################################################################################ ################################################################################ class Eigenval: """ Class containing Kohn-Sham-eigenvalues data from VASP (EIGENVAL file). """ def __init__(self): self.eigs = None self.ferw = None def from_file(self, vasp_dir='./', eig_filename='EIGENVAL'): """ Reads eigenvalues from EIGENVAL. Note that the file also contains k-points with weights. They are also stored and then used to check the consistency of files read. """ # Add a slash to the path name if necessary if vasp_dir[-1] != '/': vasp_dir += '/' f = read_lines(vasp_dir + eig_filename) # First line: only the first and the last number out of four # are used; these are 'nions' and 'ispin' sline = f.next().split() self.nq = int(sline[0]) self.ispin = int(sline[3]) # Second line: cell volume and lengths of lattice vectors (skip) sline = f.next() # Third line: temperature (skip) sline = f.next() # Fourth and fifth line: useless sline = f.next() sline = f.next() # Sixth line: NELECT, NKTOT, NBTOT sline = f.next().split() self.nelect = int(sline[0]) self.nktot = int(sline[1]) self.nband = int(sline[2]) # Set of eigenvalues and k-points self.kpts = np.zeros((self.nktot, 3)) self.kwghts = np.zeros((self.nktot,)) self.eigs = np.zeros((self.nktot, self.nband, self.ispin)) self.ferw = np.zeros((self.nktot, self.nband, self.ispin)) for ik in xrange(self.nktot): sline = f.next() # Empty line sline = f.next() # k-point info tmp = map(float, sline.split()) self.kpts[ik, :] = tmp[:3] self.kwghts[ik] = tmp[3] for ib in xrange(self.nband): sline = f.next().split() tmp = map(float, sline) assert len(tmp) == 2 * self.ispin + 1, "EIGENVAL file is incorrect (probably from old versions of VASP)" self.eigs[ik, ib, :] = tmp[1:self.ispin+1] self.ferw[ik, ib, :] = tmp[self.ispin+1:] ################################################################################ ################################################################################ # # class Doscar # ################################################################################ ################################################################################ class Doscar: """ Class containing some data from DOSCAR """ def from_file(self, vasp_dir='./', dos_filename='DOSCAR'): """ Reads only E_Fermi from DOSCAR. """ # Add a slash to the path name if necessary if vasp_dir[-1] != '/': vasp_dir += '/' f = read_lines(vasp_dir + dos_filename) # First line: NION, NION, JOBPAR, NCDIJ sline = f.next().split() self.ncdij = int(sline[3]) # Skip next 4 lines for _ in xrange(4): sline = f.next() # Sixth line: EMAX, EMIN, NEDOS, EFERMI, 1.0 sline = f.next().split() self.efermi = float(sline[3]) # TODO: implement output of SYMMCAR in VASP and read it here ################################################################ # # Reads SYMMCAR # ################################################################ def read_symmcar(vasp_dir, symm_filename='SYMMCAR'): """ Reads SYMMCAR. """ # Shorthand for simple parsing def extract_int_par(parname): return int(re.findall(parname + '\s*=\s*(\d+)', line)[-1]) # Add a slash to the path name if necessary if vasp_dir[-1] != '/': vasp_dir += '/' symmcar_exist = False sym_file = read_lines(vasp_dir + symm_filename) line = sym_file.next() nrot = extract_int_par('NROT') line = sym_file.next() ntrans = extract_int_par('NPCELL') # Lmax line = sym_file.next() lmax = extract_int_par('LMAX') mmax = 2 * lmax + 1 # Nion line = sym_file.next() nion = extract_int_par('NION') print " {0:>26} {1:d}".format("Number of rotations:", nrot) print " {0:>26} {1:d}".format("Number of translations:", ntrans) print " {0:>26} {1:d}".format("Number of ions:", nion) print " {0:>26} {1:d}".format("L_max:", lmax) rot_mats = np.zeros((nrot, lmax+1, mmax, mmax)) rot_map = np.zeros((nrot, ntrans, nion), dtype=np.int32) for irot in xrange(nrot): # Empty line line = sym_file.next() # IROT index (skip it) line = sym_file.next() # ISYMOP matrix (can be also skipped) line = sym_file.next() line = sym_file.next() line = sym_file.next() # Skip comment " Permutation map..." line = sym_file.next() # Permutations (in chunks of 20 indices per line) for it in xrange(ntrans): for ibl in xrange((nion - 1) / 20 + 1): i1 = ibl * 20 i2 = (ibl + 1) * 20 line = sym_file.next() rot_map[irot, it, i1:i2] = map(int, line.split()) for l in xrange(lmax + 1): mmax = 2 * l + 1 # Comment: "L = ..." line = sym_file.next() for m in xrange(mmax): line = sym_file.next() rot_mats[irot, l, m, :mmax] = map(float, line.split()[:mmax]) data.update({ 'nrot': nrot, 'ntrans': ntrans, 'lmax': lmax, 'nion': nion, 'sym_rots': rot_mats, 'perm_map': rot_map })