Source code for pyprocar.vaspxml.vaspxml

# -*- coding: utf-8 -*-
"""
Created on Wed Aug 19 20:49:03 2020

@author: Pedram Tavadze, Logan Lang
"""

from ..core import Structure, DensityOfStates
from numpy import array
import xml.etree.ElementTree as ET
import os
import numpy as np
import collections


[docs]class VaspXML(collections.abc.Mapping): """contains.""" def __init__(self, filename='vasprun.xml', dos_interpolation_factor=None): self.variables = {} self.dos_interpolation_factor = dos_interpolation_factor if not os.path.isfile(filename): raise ValueError('File not found ' + filename) else: self.filename = filename self.spins_dict = {'spin 1': 'Spin-up', 'spin 2': 'Spin-down'} # self.positions = None # self.stress = None self.bands = None # self.array_sizes = {} self.data = self.read()
[docs] def read(self): """ Read and parse vasprun.xml. Returns ------- TYPE DESCRIPTION. """ return parse_vasprun(self.filename)
def _get_dos_total(self): spins = list( self.data['general']['dos']['total']['array']['data'].keys()) energies = np.array(self.data['general']['dos']['total']['array'] ['data'][spins[0]])[:, 0] dos_total = {'energies': energies} for ispin in spins: dos_total[self.spins_dict[ispin]] = np.array( self.data['general']['dos']['total']['array']['data'] [ispin])[:, 1] return dos_total, list(dos_total.keys()) def _get_dos_projected(self, atoms=[]): if len(atoms) == 0: atoms = np.arange(self.initial_structure.natoms) if 'partial' in self.data['general']['dos']: dos_projected = {} ion_list = ["ion %s" % str(x + 1) for x in atoms ] # using this name as vasrun.xml uses ion # for i in range(len(ion_list)): iatom = ion_list[i] name = self.initial_structure.atoms[atoms[i]] + str(atoms[i]) spins = list(self.data['general']['dos']['partial']['array'] ['data'][iatom].keys()) energies = np.array( self.data['general']['dos']['partial']['array']['data'] [iatom][spins[0]][spins[0]])[:, 0] dos_projected[name] = {'energies': energies} for ispin in spins: dos_projected[name][self.spins_dict[ispin]] = np.array( self.data['general']['dos']['partial']['array']['data'] [iatom][ispin][ispin])[:, 1:] return dos_projected, self.data['general']['dos']['partial'][ 'array']['info'] else: print( "This calculation does not include partial density of states") return None, None @property def dos(self): energies = self.dos_total['energies'] total = [] for ispin in self.dos_total: if ispin == 'energies': continue total.append(self.dos_total[ispin]) # total = np.array(total).T return DensityOfStates( energies=energies, total=total, projected=self.dos_projected, interpolation_factor=self.dos_interpolation_factor) @property def dos_to_dict(self): """ Returns the complete density (total,projected) of states as a python dictionary """ return { 'total': self._get_dos_total(), 'projected': self._get_dos_projected() } @property def dos_total(self): """ Returns the total density of states as a pychemia.visual.DensityOfSates object """ dos_total, labels = self._get_dos_total() dos_total['energies'] -= self.fermi return dos_total @property def dos_projected(self): """ Returns the projected DOS as a multi-dimentional array, to be used in the pyprocar.core.dos object """ ret = [] dos_projected, info = self._get_dos_projected() if dos_projected is None: return None norbitals = len(info) - 1 info[0] = info[0].capitalize() labels = [] labels.append(info[0]) ret = [] for iatom in dos_projected: temp_atom = [] for iorbital in range(norbitals): temp_spin = [] for key in dos_projected[iatom]: if key == 'energies': continue temp_spin.append(dos_projected[iatom][key][:, iorbital]) temp_atom.append(temp_spin) ret.append([temp_atom]) return ret @property def kpoints(self): """ Returns the kpoints used in the calculation in form of a pychemia.core.KPoints object """ if self.data['kpoints_info']['mode'] == 'listgenerated': kpoints = dict( mode='path', kvertices=self.data['kpoints_info']['kpoint_vertices']) else: kpoints = dict(mode=self.data['kpoints_info']['mode'].lower(), grid=self.data['kpoints_info']['kgrid'], shifts=self.data['kpoints_info']['user_shift']) return kpoints @property def kpoints_list(self): """ Returns the list of kpoints and weights used in the calculation in form of a pychemia.core.KPoints object """ return dict(mode='reduced', kpoints_list=self.data['kpoints']['kpoints_list'], weights=self.data['kpoints']['k_weights']) @property def incar(self): """ Returns the incar parameters used in the calculation as pychemia.code.vasp.VaspIncar object """ return self.data['incar'] @property def vasp_parameters(self): """ Returns all of the parameters vasp has used in this calculation """ return self.data['vasp_params'] @property def potcar_info(self): """ Returns the information about pseudopotentials(POTCAR) used in this calculation """ return self.data['atom_info']['atom_types'] @property def fermi(self): """ Returns the fermi energy """ return self.data['general']['dos']['efermi'] @property def species(self): """ Returns the species in POSCAR """ return self.initial_structure.species @property def structures(self): """ Returns a list of pychemia.core.Structure representing all the ionic step structures """ symbols = [x.strip() for x in self.data['atom_info']['symbols']] structures = [] for ist in self.data['structures']: st = Structure(atoms=symbols, fractional_coordinates=ist['reduced'], lattice=ist['cell']) structures.append(st) return structures @property def structure(self): """ crystal structure of the last step """ return self.structures[-1] @property def forces(self): """ Returns all the forces in ionic steps """ return self.data['forces'] @property def initial_structure(self): """ Returns the initial Structure as a pychemia structure """ return self.structures[0] @property def final_structure(self): """ Returns the final Structure as a pychemia structure """ return self.structures[-1] @property def iteration_data(self): """ Returns a list of information in each electronic and ionic step of calculation """ return self.data['calculation'] @property def energies(self): """ Returns a list of energies in each electronic and ionic step [ionic step,electronic step, energy] """ scf_step = 0 ion_step = 0 double_counter = 1 energies = [] for calc in self.data['calculation']: if 'ewald' in calc['energy']: if double_counter == 0: double_counter += 1 scf_step += 1 elif double_counter == 1: double_counter = 0 ion_step += 1 scf_step = 1 else: scf_step += 1 energies.append([ion_step, scf_step, calc['energy']['e_0_energy']]) return energies @property def last_energy(self): """ Returns the last calculated energy of the system """ return self.energies[-1][-1] @property def energy(self): """ Returns the last calculated energy of the system """ return self.last_energy @property def convergence_electronic(self): """ Returns a boolian representing if the last electronic self-consistent calculation converged """ ediff = self.vasp_parameters['electronic']['EDIFF'] last_dE = abs(self.energies[-1][-1] - self.energies[-2][-1]) if last_dE < ediff: return True else: return False @property def convergence_ionic(self): """ Returns a boolian representing if the ionic part of the calculation converged """ energies = np.array(self.energies) nsteps = len(np.unique(np.array(self.energies)[:, 0])) if nsteps == 1: print('This calculation does not have ionic steps') return True else: ediffg = self.vasp_parameters['ionic']['EDIFFG'] if ediffg < 0: last_forces_abs = np.abs(np.array(self.forces[-1])) return not (np.any(last_forces_abs > abs(ediffg))) else: last_ionic_energy = energies[( energies[:, 0] == nsteps)][-1][-1] penultimate_ionic_energy = energies[(energies[:, 0] == ( nsteps - 1))][-1][-1] last_dE = abs(last_ionic_energy - penultimate_ionic_energy) if last_dE < ediffg: return True return False @property def convergence(self): """ Returns a boolian representing if the the electronic self-consistent and ionic calculation converged """ return (self.convergence_electronic and self.convergence_ionic) @property def is_finished(self): """ Always returns True, need to fix this according to reading the xml as if the calc is not finished we will have errors in xml parser """ # if vasprun.xml is read the calculation is finished return True def __contains__(self, x): return x in self.variables def __getitem__(self, x): return self.variables.__getitem__(x) def __iter__(self): return self.variables.__iter__() def __len__(self): return self.variables.__len__()
[docs]def text_to_bool(text): """boolians in vaspxml are stores as T or F in str format, this function coverts them to python boolians """ text = text.strip(' ') if text == 'T' or text == '.True.' or text == '.TRUE.': return True else: return False
[docs]def conv(ele, _type): """This function converts the xml text to the type specified in the attrib of xml tree """ if _type == 'string': return ele.strip() elif _type == 'int': return int(ele) elif _type == 'logical': return text_to_bool(ele) elif _type == 'float': if '*' in ele: return np.nan else: return float(ele)
[docs]def get_varray(xml_tree): """Returns an array for each varray tag in vaspxml """ ret = [] for ielement in xml_tree: ret.append([float(x) for x in ielement.text.split()]) return ret
[docs]def get_params(xml_tree, dest): """dest should be a dictionary This function is recurcive #check spelling""" for ielement in xml_tree: if ielement.tag == 'separator': dest[ielement.attrib['name'].strip()] = {} dest[ielement.attrib['name'].strip()] = get_params( ielement, dest[ielement.attrib['name']]) else: if 'type' in ielement.attrib: _type = ielement.attrib['type'] else: _type = 'float' if ielement.text is None: dest[ielement.attrib['name'].strip()] = None elif len(ielement.text.split()) > 1: dest[ielement.attrib['name'].strip()] = [ conv(x, _type) for x in ielement.text.split() ] else: dest[ielement.attrib['name'].strip()] = conv( ielement.text, _type) return dest
[docs]def get_structure(xml_tree): """Returns a dictionary of the structure """ ret = {} for ielement in xml_tree: if ielement.tag == 'crystal': for isub in ielement: if isub.attrib['name'] == 'basis': ret['cell'] = get_varray(isub) elif isub.attrib['name'] == 'volume': ret['volume'] = float(isub.text) elif isub.attrib['name'] == 'rec_basis': ret['rec_cell'] = get_varray(isub) elif ielement.tag == 'varray': if ielement.attrib['name'] == 'positions': ret['reduced'] = get_varray(ielement) return ret
[docs]def get_scstep(xml_tree): """This function extracts the self-consistent step information """ scstep = {'time': {}, 'energy': {}} for isub in xml_tree: if isub.tag == 'time': scstep['time'][isub.attrib['name']] = [ float(x) for x in isub.text.split() ] elif isub.tag == 'energy': for ienergy in isub: scstep['energy'][ienergy.attrib['name']] = float(ienergy.text) return scstep
[docs]def get_set(xml_tree, ret): """ This function will extract any element taged set recurcively""" if xml_tree[0].tag == 'r': ret[xml_tree.attrib['comment']] = get_varray(xml_tree) return ret else: ret[xml_tree.attrib['comment']] = {} for ielement in xml_tree: if ielement.tag == 'set': ret[xml_tree.attrib['comment']][ ielement.attrib['comment']] = {} ret[xml_tree.attrib['comment']][ ielement.attrib['comment']] = get_set( ielement, ret[xml_tree.attrib['comment']][ ielement.attrib['comment']]) return ret
[docs]def get_general(xml_tree, ret): """ This function will parse any element in calculatio other than the structures, scsteps""" if 'dimension' in [x.tag for x in xml_tree]: ret['info'] = [] ret['data'] = {} for ielement in xml_tree: if ielement.tag == 'field': ret['info'].append(ielement.text.strip(' ')) elif ielement.tag == 'set': for iset in ielement: ret['data'] = get_set(iset, ret['data']) return ret else: for ielement in xml_tree: if ielement.tag == 'i': if 'name' in ielement.attrib: if ielement.attrib['name'] == 'efermi': ret['efermi'] = float(ielement.text) continue ret[ielement.tag] = {} ret[ielement.tag] = get_general(ielement, ret[ielement.tag]) return ret
[docs]def parse_vasprun(vasprun): tree = ET.parse(vasprun) root = tree.getroot() calculation = [] structures = [] forces = [] stresses = [] orbital_magnetization = {} run_info = {} incar = {} general = {} kpoints_info = {} vasp_params = {} kpoints_list = [] k_weights = [] atom_info = {} for ichild in root: if ichild.tag == 'generator': for ielement in ichild: run_info[ielement.attrib['name']] = ielement.text elif ichild.tag == 'incar': incar = get_params(ichild, incar) # Skipping 1st structure which is primitive cell elif ichild.tag == 'kpoints': for ielement in ichild: if ielement.items()[0][0] == 'param': kpoints_info['mode'] = ielement.items()[0][1] if kpoints_info['mode'] == 'listgenerated': kpoints_info['kpoint_vertices'] = [] for isub in ielement: if isub.attrib == 'divisions': kpoints_info['ndivision'] = int(isub.text) else: if len(isub.text.split()) != 3: continue kpoints_info['kpoint_vertices'].append( [float(x) for x in isub.text.split()]) else: for isub in ielement: if isub.attrib['name'] == 'divisions': kpoints_info['kgrid'] = [ int(x) for x in isub.text.split() ] elif isub.attrib['name'] == 'usershift': kpoints_info['user_shift'] = [ float(x) for x in isub.text.split() ] elif isub.attrib['name'] == 'genvec1': kpoints_info['genvec1'] = [ float(x) for x in isub.text.split() ] elif isub.attrib['name'] == 'genvec2': kpoints_info['genvec2'] = [ float(x) for x in isub.text.split() ] elif isub.attrib['name'] == 'genvec3': kpoints_info['genvec3'] = [ float(x) for x in isub.text.split() ] elif isub.attrib['name'] == 'shift': kpoints_info['shift'] = [ float(x) for x in isub.text.split() ] elif ielement.items()[0][1] == 'kpointlist': for ik in ielement: kpoints_list.append( [float(x) for x in ik.text.split()]) kpoints_list = array(kpoints_list) elif ielement.items()[0][1] == 'weights': for ik in ielement: k_weights.append(float(ik.text)) k_weights = array(k_weights) # Vasp Parameters elif ichild.tag == 'parameters': vasp_params = get_params(ichild, vasp_params) # Atom info elif ichild.tag == 'atominfo': for ielement in ichild: if ielement.tag == 'atoms': atom_info['natom'] = int(ielement.text) elif ielement.tag == 'types': atom_info['nspecies'] = int(ielement.text) elif ielement.tag == 'array': if ielement.attrib['name'] == 'atoms': for isub in ielement: if isub.tag == 'set': atom_info['symbols'] = [] for isym in isub: atom_info['symbols'].append(isym[0].text) elif ielement.attrib['name'] == 'atomtypes': atom_info['atom_types'] = {} for isub in ielement: if isub.tag == 'set': for iatom in isub: atom_info['atom_types'][iatom[1].text] = {} atom_info['atom_types'][iatom[1].text][ 'natom_per_specie'] = int( iatom[0].text) atom_info['atom_types'][ iatom[1].text]['mass'] = float( iatom[2].text) atom_info['atom_types'][ iatom[1].text]['valance'] = float( iatom[3].text) atom_info['atom_types'][iatom[1].text][ 'pseudopotential'] = iatom[ 4].text.strip() elif ichild.tag == 'structure': if ichild.attrib['name'] == 'initialpos': initial_pos = get_structure(ichild) elif ichild.attrib['name'] == 'finalpos': final_pos = get_structure(ichild) elif ichild.tag == 'calculation': for ielement in ichild: if ielement.tag == 'scstep': calculation.append(get_scstep(ielement)) elif ielement.tag == 'structure': structures.append(get_structure(ielement)) elif ielement.tag == 'varray': if ielement.attrib['name'] == 'forces': forces.append(get_varray(ielement)) elif ielement.attrib['name'] == 'stress': stresses.append(get_varray(ielement)) # elif ielement.tag == 'eigenvalues': # for isub in ielement[0] : # if isub.tag == 'set': # for iset in isub : # eigen_values[iset.attrib['comment']] = {} # for ikpt in iset : # eigen_values[iset.attrib['comment']][ikpt.attrib['comment']] = get_varray(ikpt) elif ielement.tag == 'separator': if ielement.attrib['name'] == "orbital magnetization": for isub in ielement: orbital_magnetization[isub.attrib['name']] = [ float(x) for x in isub.text.split() ] # elif ielement.tag == 'dos': # for isub in ielement : # if 'name' in isub.attrib: # if isub.attrib['name'] == 'efermi' : # dos['efermi'] = float(isub.text) # else : # dos[isub.tag] = {} # dos[isub.tag]['info'] = [] # for iset in isub[0] : # if iset.tag == 'set' : # for isub_set in iset: # dos[isub.tag] = get_set(isub_set,dos[isub.tag]) # elif iset.tag == 'field' : # dos[isub.tag]['info'].append(iset.text.strip(' ')) else: general[ielement.tag] = {} general[ielement.tag] = get_general( ielement, general[ielement.tag]) # NEED TO ADD ORBITAL MAGNETIZATION return { 'calculation': calculation, 'structures': structures, 'forces': forces, 'run_info': run_info, 'incar': incar, 'general': general, 'kpoints_info': kpoints_info, 'vasp_params': vasp_params, 'kpoints': { 'kpoints_list': kpoints_list, 'k_weights': k_weights }, 'atom_info': atom_info }