Read_Siesta_XML.m

Go to the documentation of this file.

%> @file Read_Siesta_XML.m
%> @brief Parses the SIESTA XML output and looks for relevant data
%> @Available
%> EQuUs v4.9 or later
%>
%> @param 'xmlfile' the xml output of the SIESTA calculation.
%> @param 'fields' what entries to read from the file.
    function xml_params = Read_Siesta_XML(xmlfile, fields, readanyway)

     %SOME CONSTANTS USED IN SIESTA
     RytoeV = 1/13.60580;                      % Rydberg to eV
     Bohr   =  0.529177;                       % Angstrom

     %THAT IS THE XML FILE THAT HOPEFULLY COMES OUT OF SIESTA AS DEFAULT
     if exist(xmlfile, 'file') == 0
          disp(['ERROR: There is no file: ' xmlfile])
        throw(['ERROR: There is no file: ' xmlfile])
    end
    
    xml_params = [];
    
     tmp = read_SiestaXML( xmlfile, fields );
    
    %%convert tmp to make in unfallable
    if isfield(tmp, 'property')
        if isstruct(tmp.property)
            property{1} = tmp.property;
        else
            property = tmp.property;
        end
    else
        property = nan;
    end
    
    if isfield(tmp, 'propertyList')
        if isstruct(tmp.propertyList)
            propertyList{1} = tmp.propertyList;
        else
            propertyList = tmp.propertyList;
        end
    else
        propertyList = nan;
    end
    
    if isfield(tmp, 'module')
        module{1} = tmp.module;
    else
        module = nan;
    end
    
    if iscell(propertyList)
        for i=1:length(propertyList)
            if strcmp(propertyList{i}.dictRef,'siesta:kpoints')
                    %READ MONKHORST PACK K POINTS  FROM "system_label.xml" ('K')
                    Tkpoints=propertyList{i}.subnodes.kpoint;
                    if length(Tkpoints)==1
                        kp(1,1:4)=[0,0,0,1];
                    else
                        for ik=1:length(Tkpoints)
                            kp(ik,1:4) = [str2num(Tkpoints{ik}.coords) str2num(Tkpoints{ik}.weight)];
                        end
                    end
                    xml_params.kpoints = kp;
            end
        end
    end

     %% READ FERMI ENERGY FROM FILE
    if iscell(property)
        for i=1:length(property)
            if strcmp(property{i}.dictRef,'siesta:E_Fermi')
                EF = str2num(property{i}.subnodes.scalar.value);
                %EF = EF - EF_shift;
                xml_params.E_Fermi = EF;
            end
        end
    end
    
    %% READ K-SUPERCELL SIZE FROM FILE
    if iscell(property)
        for i=1:length(property)
            if strcmp(property{i}.dictRef,'siesta:kscell')
                xml_params.kscell = reshape(str2num(property{i}.subnodes.matrix.value), [3,3]);
            end
        end
    end
    
     %COORDS and %LATTICE VECTORS
    if  iscell(module)
        subnodes = module{1}.subnodes;
        coords=[];
        Tatom=subnodes.molecule.subnodes.atomArray.subnodes.atom;
        for i=1:length(Tatom)
          coords = [coords; str2num(Tatom{i}.x3) str2num(Tatom{i}.y3) str2num(Tatom{i}.z3)];
          %elmtxv(i) = get_atomindex(Tatom{i}.elementType);
          spxv(i,:) = Tatom{i}.elementType;
          elmtxv(i) = 71;
        end
        xml_params.coordinates = coords;
        xml_params.elmtxv = elmtxv;
        xml_params.spxv = spxv;
        
        vect(1,:)=str2num(subnodes.lattice.subnodes.latticeVector{1}.value);
        vect(2,:)=str2num(subnodes.lattice.subnodes.latticeVector{2}.value);
        vect(3,:)=str2num(subnodes.lattice.subnodes.latticeVector{3}.value);

        %if module.subnodes.lattice.subnodes.latticeVector{1}.units == 'siestaUnits:angstrom' %THIS IS IN BOHR :|
        if strcmp(subnodes.lattice.subnodes.latticeVector{1}.units,'siestaUnits:Ang')
            coords = coords/Bohr;
            vect=vect/Bohr;
        end
        xml_params.Lattice_vectors = vect;
    end
     %coords = coords*vect;
    

    end