Source code for pynao.m_siesta_xml

import re
import numpy as np

from pynao.m_color import color as bc
from pynao.m_siesta_units import siesta_conv_coefficients
import xml.etree.ElementTree as ET

pref = "{http://www.xml-cml.org/schema}"

def extract_energy_parameter(parent, field):
    """
    Extract an energy parameter from the siesta.xml file and convert it
    to Hartree
    """

    param = parent.find(field)[0]
    if "units" not in param.keys():
        raise ValueError("units not in param keys")

    unit = param.get("units").split(":")[1]
    if unit == "Ry":
        unit = "ry"
    elif unit == "eV":
        unit = "ev"
    elif unit == "Ha":
        unit = "ha"

    if unit == "ha":
        return float(param.text)
    else:
        key = "{}2ha".format(unit)
        if key not in siesta_conv_coefficients.keys():
            raise ValueError("No conversion defined between {} and ha?".format(unit))

        return float(param.text)*siesta_conv_coefficients[key]

def siesta_xml(fname="siesta.xml"):
    """
    Read the siesta.xml file
    """

    try :
        tree = ET.parse(fname)
    except:
        raise SystemError(fname+" cannot be parsed: calculation did not finish?")

    roo = tree.getroot()
    params = roo.find(pref + "parameterList")
    fin = roo.find(pref+"module[@title='Finalization']")
    mol = fin.find(pref+"molecule")
    coo = mol.find(pref+"atomArray")
    atoms = coo.findall(pref+"atom")

    natoms = len(atoms)
    atom2ref = [atom.attrib["ref"] for atom in atoms]
    atom2sp = [int(re.findall(r'\d+', ref)[0])-1 for ref in atom2ref]
    nspecies = len(set(atom2ref))
    atom2elem = [atom.attrib["elementType"] for atom in atoms]
    telec = extract_energy_parameter(params, pref+"parameter[@dictRef='siesta:etemp']")

    sp2elem = [None]*nspecies
    for a in range(len(atom2elem)):
        sp2elem[atom2sp[a]] = atom2elem[a]

    # read atoms coordinates
    x3 = list(map(float, [atom.attrib["x3"] for atom in atoms]))
    y3 = list(map(float, [atom.attrib["y3"] for atom in atoms]))
    z3 = list(map(float, [atom.attrib["z3"] for atom in atoms]))
    atom2coord = np.empty((natoms,3), dtype='double')
    for a in range(natoms):
        atom2coord[a,0:3] = [x3[a],y3[a],z3[a]]
    atom2coord = atom2coord*siesta_conv_coefficients["ang2bohr"]

    eigvals = fin.find(pref+"propertyList[@title='Eigenvalues']")

    lat = fin.find(pref+"lattice[@dictRef='siesta:ucell']")
    ucell = np.empty((3, 3), dtype='double')
    for i in range(len(lat)):
        ucell[i,0:3] = list(map(float, filter(None, re.split(r'\s+|=', lat[i].text))))
    ucell = ucell * siesta_conv_coefficients["ang2bohr"]

    fermi_energy = extract_energy_parameter(eigvals,
                                            pref+"property[@dictRef='siesta:E_Fermi']")

    nkp = int(eigvals.find(pref+"property[@dictRef='siesta:nkpoints']")[0].text)
    k2xyzw = np.empty((nkp,4), dtype='double')
      
    kpt_band = eigvals.findall(pref+"propertyList[@dictRef='siesta:kpt_band']")
    nspin = len(kpt_band)
    norbs = int(kpt_band[0].find(pref+"property[@dictRef='siesta:eigenenergies']")[0].attrib['size'])
    ksn2e = np.empty((nkp, nspin, norbs), dtype='double')

    siesta_ev2ha = siesta_conv_coefficients["ev2ha"]
    for s in range(nspin):
        eigv = kpt_band[s].findall(pref + "property[@dictRef='siesta:eigenenergies']")
        kpnt = kpt_band[s].findall(pref + "kpoint")
        for k in range(nkp):
            ksn2e[k, s, 0:norbs] = \
                    list(map(lambda x : siesta_ev2ha*float(x), \
                             filter(None, re.split(r'\s+|=', eigv[k][0].text))))
            k2xyzw[k, 0:3] = list(map(float, filter(None, re.split(r'\s+|=', \
                                                    kpnt[k].attrib['coords']))))
            k2xyzw[k, 3] = float(kpnt[k].attrib['weight'])

    return dict({'fermi_energy':fermi_energy,
                 'telec': telec,
                 'atom2coord':atom2coord, 
                 'atom2sp':atom2sp, 
                 'sp2elem':sp2elem, 
                 'k2xyzw':k2xyzw, 
                 'ksn2e':ksn2e, 
                 'ucell':ucell})