%This script extracts energies and material composition from an xml file 
%containing compounds containing only Ti and O atoms directly downloaded 
%from the open quantum materials database (www.oqmd.org)

%read in file using function parseXML
XMLfile='formationenergy.xml';
TiO_data=parseXML(XMLfile);

%parseXML reads data in as a nested structure
%to extract data from structure need to dive into the matlab structure
all_data=TiO_data.Children(3).Children;

%determine number of compounds in xml file
[d,len_all_data]=size(all_data);

%intialize arrays for storing data (makes code more efficient)
energies=1:len_all_data;
atoms=[1:len_all_data;1:len_all_data]';
composition=strings();

%loop over each material in the xml file and save the name of the structure
%(its chemical composition) and the associated DFT energies in eV/atom.
%This loop also processes the name to determine number of atoms of each 
%type in the material
for material=1:len_all_data
    %energy data extraction from data structure
    energies(material)=eval(all_data(material).Children(11).Children.Data);

    %composition data extraction from data structure
    composition(material)=string(all_data(material).Children(1).Children.Data);
    
    %process name to determine number of atoms
    atoms_tmp=strsplit(composition(material),'Ti'); %split at Ti
    atoms_tmp=strsplit(atoms_tmp(2),'O'); %split second part of string at O
    if atoms_tmp(1)==""
        ti_atoms=1; %if no number must be 1 Ti atom
    else
        ti_atoms=eval(atoms_tmp(1)); %otherwise number of atoms taken from string
    end
    if atoms_tmp(2)==""
        o_atoms=1; %if no number must be 1 O atom
    else
        o_atoms=eval(atoms_tmp(2)); %otherwise number of atoms taken from string
    end
    %save values in array 'atoms' for each material
    atoms(material,1)=ti_atoms; 
    atoms(material,2)=o_atoms;
end

%determine average energy for TiO2 materials
tio2_sum=0; %initialize variables
count=0;
tio2_energies=[];
%loop over each material and find when the number of Ti atoms is 1 and the
%number of oxygen atoms is 2
for material=1:len_all_data
    if atoms(material,1)==1 && atoms(material,2)==2
        tio2_sum=tio2_sum+energies(material); %add to running sum
        count=count+1; %update count
        tio2_energies(count)=energies(material); %what does this store?
    end
end
tio2_avg=tio2_sum/count; %average is sum/count

%print formatted average results
fprintf('The average energy of all Ti and O containing compounds is %f eV\n',mean(energies));
fprintf('The average energy for TiO_2 compounds is %f eV\n',tio2_avg);

%generate plots
plot(energies)
xlabel('Material Index')
ylabel('Energy (eV)')
f2=figure;
plot(atoms(:,2),energies,'o')
xlabel('Number of oxygen atoms per formula unit')
ylabel('Energy (eV)')