Investigating the Energetic Ordering of Stable and Metastable \ce{TiO_2} Polymorphs Using DFT+$U$ and Hybrid Functionals
\date{\today} \maketitle
\tableofcontents \newpage
The Supporting Information document for the article, “Investigating the Energetic Ordering of Stable and Metastable \ce{TiO_2} Polymorphs Using DFT+$U$ and Hybrid Functionals”, provides all of the information needed to reproduce the results presented in this article. Additionally, this document will provide explanations and additional data for analyses completed in the article to substantiate several of the insights and claims made in the article itself. Both this Supporting Information document and the article explained by it were prepared using the org-mode mark-up language, which facilitates the formatting, drafting, organization, and quick PDF conversion of manuscripts prepared in \LaTeX cite:Dominik2010. The original source of this document can be found here: attachfile:supporting-information.org.
Sections of this Supporting Information document provide additional information regarding how the results contained within the article were achieved, illustrating how calculations yielding the results depicted were completed, detailing how results were aggregated and visualized as plotted data, and verifying the claims made by sections alluding to a supporting information document. Minimally, this document will contain information concerning the input files used to perform calculations, such that the replication of each different type of calculation performed in the Results and Discussion section of the article is possible. In addition, sections will also contain Python cite:Python scripts capable of generating all plots within the article, accompanied by the numerical values of the formation energies depicted on those plots.
All of the direct and indirect allusions to a Supporting Information document in the article are addressed in this document. These allusions can be made for multiple reasons, including the inability to efficiently portray all combinations of data within the article itself or the need to more broadly contextualize results calculated for several TiO2 polymorphs with respect to other BO2 polymorphs (B = metal cation) covered in past research cite:Mehta2014. For example, the number of combinations of different pseudopotentials, functionals, and software packages used to calculate formation energy orderings over ranges of
All of the results of calculations and input parameters responsible for these calculations are stored in a database consisting of relations that are largely consistent with either just Boyce-Codd Normal Form (BCNF) cite:Codd1972 or both BCNF and 4th Normal Form (4NF) cite:Fagin1977. The information needed to generate and sort all input files and results of calculations will be stored in several files of the .csv file format, which can be imported into a wide variety of relational database management systems (RDBMS) such as MySQL cite:MySQL or SQLite cite:sqlite to appropriately process information. In the following section, the relational schema categorizing the data will be detailed, illustrating how the information can be recalled. To facilitate the recollection of information needed for generating plots contained within the article, relational algebraic expressions and matching SQL queries suitable for generating all of the plots contained in the article via the database are provided.
The relational hierarchy of the database detailed below can be illustrated using the class diagram generation feature of the PlantUML Java tool, cite:PlantUML representing each relation (R) as an object in the diagram, each attribute as an element listed within an object or relation, and each key as a morphism or link between two objects or relations.
In the following sections, each feature of the data structure illustrated above is described in detail. Assembly of the relational schema dictating the structure of the database presented here proceeds from a decomposition directly corresponding to that shown in previous research cite:Curnan2014, thus only the fully decomposed relations are characterized below with their constituent attributes explained. This schema, though inspired by the design principles of BCNF and 4NF, also has several attributes organized so as to ease the querying of necessary, polymorph-specific energetic and volumetric information needed to reproduce and verify presented results. Excluding the relations within this schema devoted to defining structures, this schema can also be considered a set of hierarchical depths of a tree diagram, in which each depth contains organized sets of parameters pertaining to a different step in the process of calculating results to be represented in plots within the article. Prior to reviewing these hiearchical depths, the relations within this schema relating to the definition of structures (located on the leftmost branch of the relation Metadata in Figure ref:fig-HierarchyFormE) are characterized in the following sections.
The Symmetry relation contains the attributes “Sym” and “Struct”, the former of which is a key attribute and the latter of which relates Symmetry to Structure. The attribute “Sym”, which features TEXT formatted instances, refers to the crystallographic point group of a particular periodic structure. The attribute “Struct”, which also features TEXT formatted instances, refers to the Bravais lattice of a particular periodic, crystallographic structure. Values for “Sym” relating to calculations detailed in this article or its Supporting Information document sections include “P4|2/mnm”, “I4|1/amd”, “Pbcn”, “Pbca”, “Pnma”, “Pa-3”, “Fm-3m”, and “P2|1/c”. Corresponding values for “Struct” include “tetragonal”, “orthorhombic”, “cubic”, and “monoclinic”. Note that vertical bar characters (|) are used instead of underscore characters (_) to specify subscripts in symmetry groups recorded as data instances in this database, for underscore characters serve as string literals (i.e.: wildcards) in SQL-based string comparison commands (involving data instances, not attributes or relations) cite:MySQL.
The Structure relation contains the attributes “Struct”, Morph”, and “SID”, with the key attribute “SID” relating Structure to Composition1. In this article, the attribute “Morph”, which features TEXT formatted instances, refers to the structural polymorph represented by a particular system. The attribute “SID”, which features INTEGER formatted instances, is a unique integer assigned to each lattice structure, serving as an identification of the lattice itself independent of the compositions of particular lattice sites (i.e.: a structural ID, or sID). Values for “Morph” relating to calculations detailed in article or its Supporting Information document sections include “Rutile”, “Anatase”, “Columbite”, “Brookite”, “Cotunnite”, “Pyrite”, “Fluorite”, and “Baddeleyite”.
The relation Composition1 contains the attributes “MID”, “At1”, “Stoich1”, and “SID”. The attribute “MID”, which features INTEGER formatted instances, is a unique integer assigned to each set of atomic compositions ordered in a periodic lattice (also known as a motif identifier, or mID). Combination of the lattice and motif identifications allows identification of any periodic solid, thus in combination, “SID” and “MID” serve as an attribute key for the relation Composition1. These attributes are linked to the relation Metadata, which is the first relation in queries performed later that contains data explicitly involving calculations. Composition1 contains information concerning the portion of the motif represented by the first atomic type in a calculation, serving as the only information on motif directly linked to other relations in the schema. The composition of the first atomic type, namely “At1”, features TEXT formatted instances. The attribute “Stoich1”, which features INTEGER formatted instances, indicates the number of atoms of the first atomic type that are present in a system. Values for “At1” relating to calculations detailed in article or its Supporting Information document sections include “Ti”, “V”, “Ru”, and “Ir”.
Each relation containing information concerning each subsequent atomic type is directly linked to the relation directly preceding it via the key attribute “MID”, which is shared for all atomic types represented by a single motif within a schema. Given that each system studied in this article contains two atomic types (i.e.: a metal cation, generally Ti, and O), Composition2 is the terminal relation in the structural branch of the schema, which represents information concerning the portion of the motif represented by the second atomic type in a calculation. Generally, the metal cation (usually Ti) in each calculation is represented in Composition1, while the O atom of each calculation is usually represented in Composition2. The attributes “At2” and “Stoich2” represent the second atomic type and the number of atoms of that type in a calculation, respectively. They possess the same formatting as corresponding instances of “At1” and “Stoich1”.
The hierarchical depths within the schema not relating to the leftmost branch depicting structural information in Figure ref:fig-HierarchyFormE are stated as prefixes in the relations constituting them, namely “Parameters”, “CalcResult”, “FormE”, and the head node (relative to nodes containing results) “Metadata”. Relations within the depth “Parameters” describe the relaxed or final structural information of a system (PosFinal), necessary input file values needed to achieve those structures and related energetic properties (Input), or information needed to calculate response matrices (Response). These relations divide data based on the software platform employed to calculate them, using VASP or Quantum Espresso (abbreviated as QE) in this study. Relations of this form are linked to “CalcResult” relations via the key attribute “Energy” (contains FLOAT formatted data instances), the DFT total energy of a single calculation that serves as a unique identifier of each calculation. Unless otherwise noted, all instances contained within relations prefixed by “Parameters” are TEXT formatted instances.
The relations ParametersPosFinalVASP and ParametersPosFinalQE store information concerning the final, relaxed structural information of a system of interest, given its software platform. For VASP calculations, this includes all information in a CONTCAR file, including primitive lattice vectors (“PrimVec1”, “PrimVec2”, “PrimVec3”) and the magnitude scale for those vectors (“Scale”). Additionally, the atomic coordinates of each atom are instances within the attributes “Coord\textit{N}”, where \textit{N} represents a number between one and the total number of atoms featured in the largest system studied in this article. For QE calculations, this only includes the “Scale” and “Coord\textit{N}” attributes due to differences in the formats of the input and output files presented in VASP and QE. The “Scale” features FLOAT formatted instances, while “PrimVec” and “Coord\textit{N}” attributes feature TEXT formatted instances. TEXT formatted instances are employed in this case, as atomic coordinates and primitive lattice vectors feature numbers separated by semicolons, in order to differentiate between different coordinates or indices, respectively. For both software platforms, atomic coordinates are ordered in accordance with atomic type instances, which always start with the metal cation species (usually Ti) and end with the O species. In linear response calculations, the first atomic coordinate present per system always represents the perturbed metal cation.
The relations ParametersInputVASP and ParametersInputQE store information concerning the other input file data necessary for reproducing data in the article beyond relaxed structural information. For VASP calculations, this includes the tags associated with the use of particular pseudopotentials (“Pseudopotential1” and “Pseudopotential2”) in the order that they are presented in a POTCAR file, the
The relation ParametersResponseVASP stores information concerning the outputted data received from calculations involving the linear response approach. Each atom, which is classified by its atomic index (INTEGER formatting, attribute “AtomIndex”), is mapped to its
Relations within the depth “CalcResult” represent individual calculation results achieved from the input information described in “Parameters”. These relations contain data based on the types of calculation results placed within them, namely formation energy calculations (“Energetics” suffix) primarily applied to plots depicting incremental change of
The relation CalcResultEnergetics contains the attributes “EOpt”, “Energy”, “Volume”, and “NumAtoms”. Attributes serving as keys to link one relation to another, namely “Energy” and “EOpt”, represent the total, final DFT energy associated with a calculation and the fitted energy associated with a set of calculations of a shared system performed over a structural displacement, respectively. The attribute “Volume” represents the total relaxed system volume associated with a single calculation, whereas the attribute “NumAtoms” represents the number of atoms within that calculation. Note that all of these attributes feature FLOAT formatted instances.
The relation CalcResultLRCalc contains the attributes “CID”, “U3dIn”, “U3dOut”, “Energy”, “PturbMin”, “PturbMax”, Chi0\textit{N150-P150} (or explicitly, “Chi0N150”, “Chi0N100”, “Chi0N050”, “Chi0P000”, “Chi0P050”, “Chi0P100”, “Chi0P150”), and Chi\textit{N150-P150} (or explicitly, “ChiN150”, “ChiN100”, “ChiN050”, “ChiP000”, “ChiP050”, “ChiP100”, “ChiP150”). Attributes serving as keys to link one relation to another, namely “Energy” and “CID”, represent the total, final DFT energy associated with a calculation and the unique identifiers of particular plotted data points in Metadata (to be explained later), respectively. The attribute “U3dOut” represents the first-principles resolved value of
Relations within the depth “FormE” represent sets of calculation results taken over displacement of a structural parameter, generally cell volume or “Volume”. Data within this set of relations can be directly manipulated by orders of operation to reproduce data visualized in plots contained within the article or its Supporting Information document. These relations are divided based on the type of data stored within them, whether this data is used to produce plots featuring incremental change of a
The relation FormEURange contains the attributes “CID”, “UValue”, “VOpt”, and “EOpt”. With the exception of “CID”, all instances within these attributes are FLOAT formatted. Attributes “CID” and “EOpt” have already been respectively defined as attributes uniquely identifying sets of calculations within particular articles and uniquely identifying sets of calculations modeling a particular system over a structural displacement. The attribute “VOpt” is the complement to the Birch-Murnaghan fitted energy “EOpt” of a system, defining the fitted volume of the system. The attribute “UValue” defines the constant value of \textit{U}\textit{3d} at which a set of calculations plotted in a figure incrementally changing
The relation FormEHybrid contains the attributes “CID”, “FractionHF”, “VOpt”, and “EOpt”. With the exception of “CID”, all instances within these attributes are FLOAT formatted. All attributes besides “FractionHF” have been formerly defined in either Section FormEURange or another section. Attribute “FractionHF” represents the fraction of exact exchange employed in a hybrid functional calculation using either the HSE06 or PBE0 functional. In VASP calculations, this quantity is available in the input file data presented in ParametersInputVASP as each TEXT formatted instance of attribute “AEXX”. However, for the purposes of easily retrieving common data, the “FractionHF” attribute has been reproduced in this relation.
The relation FormEFigure contains the attributes “CID” and “Figure”. The attribute “Figure” links each cID to the figure (figure name minus file extension) in which it is represented in the article, containing TEXT formatted instances that give the name of each figure.
The relation Metadata serves as the head node in the relational schema containing all of the data within this article and its supporting information section. Despite slight redundancy, it contains attributes spanning all of the unique identifiers (cID, mID, sID) linking structures to calculated results, in addition to several attributes that conveniently delineate between calculations that are featured in different Figures and underscore different arguments within the article. These additional attributes of convenience contain all TEXT formatted attributes. With regard to these attributes, “PseudoB” details the type of metal cation (B, usually Ti) pseudopotential employed in a particular calculation, “PseudoO” details the type of O pseudopotential employed in the same calculation, “Functional” specifies the type of functional used in that calculation, “Software” indicates whether QE or VASP was used to complete the calculation, and “Method” indicates the type of calculation mapped to a particular cID value.
Values for “PseudoB” relating to calculations detailed in the article or its Supporting Information document sections include “PAW-B” (standard VASP PAW pseudopotential, where B = Ti, V, Ru, or Ir), “PAW-Ti-pv” (
The .csv files accompanying this article and its supporting information file can be imported into any RDBMS platform to reproduce results or perform other operations on aggregated data. However, in order to facilitate the easy integration of .csv files into org-mode and enable full reproducibility of all plots from initial data within this contained document, this article will assemble all .csv files into a single .sqlite file, which can be read into Python scripts within the “Plot Generation” subsections of this document to generate all plots contained within the article. Plots and tables not associated with plots in the article can also be reproduced using this file. Upon being equipped with proper extensions, Mozilla Firefox cite:Mozilla can also import the .sqlite file, perform queries on its contained data, and provide a graphic-user interface suitable for browsing and examining the data.
The created database can be downloaded here. ITEOO_data.sqlite: attachfile:ITEOO_data.sqlite
In order to generate a single .sqlite file containing all .csv file data needed to reproduce all plots listed in the following section, execute the Python scripts listed below in the order that they are presented. Note that the .csv files and this document should be located in the current working directory prior to executing these scripts, which upload tables that correspond to aforementioned RDBMS subrelations one at a time:
Data corresponding to this relation is stored in Structure.csv: attachfile:Structure.csv
import sqlite3
import os
# We start from scratch. Delete the database if it already exists.
if os.path.exists('ITEOO_data.sqlite'):
os.remove('ITEOO_data.sqlite')
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table Structure(Morph TEXT, Struct TEXT, SID INTEGER PRIMARY KEY)''')
with open('Structure.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into Structure(Morph,Struct,SID)
VALUES(?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in Symmetry.csv: attachfile:Symmetry.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table Symmetry(Sym TEXT, Struct TEXT)''')
with open('Symmetry.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into Symmetry(Sym, Struct)
VALUES(?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in Composition1.csv: attachfile:Composition1.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table Composition1(At1 TEXT, Stoich1 INTEGER, SID INTEGER, MID INTEGER)''')
with open('Composition1.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into Composition1(At1,Stoich1,SID,MID)
VALUES(?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in Composition2.csv: attachfile:Composition2.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table Composition2(At2 TEXT, Stoich2 INTEGER, SID INTEGER, MID INTEGER)''')
with open('Composition2.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into Composition2(At2,Stoich2,MID,SID)
VALUES(?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in Metadata.csv: attachfile:Metadata.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table Metadata(PseudoB TEXT, PseudoO TEXT,
Functional TEXT, Software TEXT, Method TEXT, CID TEXT,
MID INTEGER, SID INTEGER)''')
with open('Metadata.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into Metadata(PseudoB, PseudoO,
Functional, Software, Method, CID, MID, SID)
VALUES(?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in FormEURange.csv: attachfile:FormEURange.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table FormEURange(UValue FLOAT, VOpt FLOAT, CID TEXT, EOpt FLOAT)''')
with open('FormEURange.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into FormEURange(UValue, VOpt, CID, EOpt)
VALUES(?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in FormEHFRange.csv: attachfile:FormEHFRange.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table FormEHFRange(FractionHF FLOAT, VOpt FLOAT, CID TEXT, EOpt FLOAT)''')
with open('FormEHFRange.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into FormEHFRange(FractionHF,VOpt,CID,EOpt)
VALUES(?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in CalcResultLRCalc.csv: attachfile:CalcResultLRCalc.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table CalcResultLRCalc(CID TEXT, U3dIn FLOAT,
U3dOut FLOAT, Energy FLOAT, Chi0N150 FLOAT, Chi0N100 FLOAT,
Chi0N050 FLOAT, Chi0P000 FLOAT, Chi0P050 FLOAT, Chi0P100 FLOAT,
Chi0P150 FLOAT, ChiN150 FLOAT, ChiN100 FLOAT, ChiN050 FLOAT,
ChiP000 FLOAT, ChiP050 FLOAT, ChiP100 FLOAT, ChiP150 FLOAT,
PturbMax FLOAT, PturbMin FLOAT)''')
with open('CalcResultLRCalc.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into CalcResultLRCalc(CID, U3dIn, U3dOut,
Energy, Chi0N150, Chi0N100, Chi0N050, Chi0P000, Chi0P050, Chi0P100,
Chi0P150, ChiN150, ChiN100, ChiN050, ChiP000, ChiP050, ChiP100,
ChiP150, PturbMax, PturbMin)
VALUES(?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in CalcResultEnergetics.csv: attachfile:CalcResultEnergetics.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table CalcResultEnergetics(Volume FLOAT, NumAtoms FLOAT, EOpt FLOAT, Energy FLOAT)''')
with open('CalcResultEnergetics.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into CalcResultEnergetics(Volume, NumAtoms, EOpt, Energy)
VALUES(?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in ParametersResponseVASP.csv: attachfile:ParametersResponseVASP.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table ParametersResponseVASP(Energy FLOAT,
CalcType TEXT, AtomIndex INTEGER, DChargeRHigh FLOAT, DChargeRLow FLOAT,
DChargeRCenter FLOAT, PChargeRHigh FLOAT, PChargeRLow FLOAT,
PChargeRCenter FLOAT, DMagP FLOAT, PMagP FLOAT)''')
with open('ParametersResponseVASP.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into ParametersResponseVASP(Energy,
CalcType, AtomIndex, DChargeRHigh, DChargeRLow, DChargeRCenter,
PChargeRHigh, PChargeRLow, PChargeRCenter, DMagP, PMagP)
VALUES(?,?,?,?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in ParametersInputVASP.csv: attachfile:ParametersInputVASP.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table ParametersInputVASP(Energy FLOAT, NELMDL TEXT,
SYMPREC TEXT, ISYM TEXT, IBRION TEXT, SIGMA TEXT, NELMIN TEXT,
KPointSampling TEXT, KPoints TEXT, LDAUL TEXT, LDAUJ TEXT, ENCUT TEXT,
ISIF TEXT, ICHARG TEXT, GGA TEXT, LDAUPRINT TEXT, LDAUU TEXT,
Pseudopotential1 TEXT, Pseudopotential2 TEXT, NELM TEXT, NSW TEXT,
LASPH TEXT, MAXMIX TEXT, EDIFF TEXT, ISMEAR TEXT, ISTART TEXT,
LDAU TEXT, LMAXMIX TEXT, EDIFFG TEXT, ISPIN TEXT, LDAUTYPE TEXT,
LORBIT TEXT, AEXX TEXT, NKRED TEXT, PRECFOCK TEXT, NBANDS TEXT,
LMAXFOCK TEXT, ALGO TEXT, LHFCALC TEXT, TIME TEXT, HFSCREEN TEXT)''')
with open('ParametersInputVASP.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into ParametersInputVASP(Energy, NELMDL,
SYMPREC, ISYM, IBRION, SIGMA, NELMIN, KPointSampling, KPoints, LDAUL, LDAUJ,
ENCUT, ISIF, ICHARG, GGA, LDAUPRINT, LDAUU, Pseudopotential1, Pseudopotential2,
NELM, NSW, LASPH, MAXMIX, EDIFF, ISMEAR, ISTART, LDAU, LMAXMIX, EDIFFG, ISPIN,
LDAUTYPE, LORBIT, AEXX, NKRED, PRECFOCK, NBANDS, LMAXFOCK, ALGO, LHFCALC, TIME, HFSCREEN)
VALUES(?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,
?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in ParametersInputQE.csv: attachfile:ParametersInputQE.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table ParametersInputQE(Energy FLOAT,
calculation TEXT, verbosity TEXT, restartmode TEXT, pseudodir TEXT,
outdir TEXT, prefix TEXT, nstep TEXT, wfcollect TEXT,
forcconvthr TEXT, etotconvthr TEXT, ibrav TEXT, nat TEXT,
ntyp TEXT, ecutwfc TEXT, ecutrho TEXT, nspin TEXT,
startingmagnetization1 TEXT, occupations TEXT, smearing TEXT,
degauss TEXT, ldaplusu TEXT, ldaplusukind TEXT, celldm1 TEXT,
celldm2 TEXT, celldm3 TEXT, HubbardU1 TEXT, HubbardU2 TEXT,
electronmaxstep TEXT, convthr TEXT, diagonalization TEXT,
diagothrinit TEXT, diagofullacc TEXT, startingwfc TEXT,
mixingmode TEXT, mixingbeta TEXT, iondynamics TEXT,
upscale TEXT, Kpoints TEXT, Pseudopotential1 TEXT,
Pseudopotential2 TEXT)''')
with open('ParametersInputQE.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into ParametersInputQE(Energy, calculation,
verbosity, restartmode, pseudodir, outdir, prefix, nstep, wfcollect,
forcconvthr, etotconvthr, ibrav, nat, ntyp, ecutwfc, ecutrho, nspin,
startingmagnetization1, occupations, smearing, degauss, ldaplusu,
ldaplusukind, celldm1, celldm2, celldm3, HubbardU1,
HubbardU2, electronmaxstep, convthr, diagonalization,
diagothrinit, diagofullacc, startingwfc, mixingmode, mixingbeta,
iondynamics, upscale, KPoints, Pseudopotential1, Pseudopotential2)
VALUES(?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,
?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in ParametersPosFinalVASP.csv: attachfile:ParametersPosFinalVASP.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table ParametersPosFinalVASP(Energy FLOAT, Scale FLOAT,
PrimVec1 TEXT, PrimVec2 TEXT, PrimVec3 TEXT, Coord1 TEXT, Coord2 TEXT,
Coord3 TEXT, Coord4 TEXT, Coord5 TEXT, Coord6 TEXT, Coord7 TEXT,
Coord8 TEXT, Coord9 TEXT, Coord10 TEXT, Coord11 TEXT, Coord12 TEXT,
Coord13 TEXT, Coord14 TEXT, Coord15 TEXT, Coord16 TEXT, Coord17 TEXT,
Coord18 TEXT, Coord19 TEXT, Coord20 TEXT, Coord21 TEXT, Coord22 TEXT,
Coord23 TEXT, Coord24 TEXT, Coord25 TEXT, Coord26 TEXT, Coord27 TEXT,
Coord28 TEXT, Coord29 TEXT, Coord30 TEXT, Coord31 TEXT, Coord32 TEXT,
Coord33 TEXT, Coord34 TEXT, Coord35 TEXT, Coord36 TEXT, Coord37 TEXT,
Coord38 TEXT, Coord39 TEXT, Coord40 TEXT, Coord41 TEXT, Coord42 TEXT,
Coord43 TEXT, Coord44 TEXT, Coord45 TEXT, Coord46 TEXT, Coord47 TEXT,
Coord48 TEXT, Coord49 TEXT, Coord50 TEXT, Coord51 TEXT, Coord52 TEXT,
Coord53 TEXT, Coord54 TEXT, Coord55 TEXT, Coord56 TEXT, Coord57 TEXT,
Coord58 TEXT, Coord59 TEXT, Coord60 TEXT, Coord61 TEXT, Coord62 TEXT,
Coord63 TEXT, Coord64 TEXT, Coord65 TEXT, Coord66 TEXT, Coord67 TEXT,
Coord68 TEXT, Coord69 TEXT, Coord70 TEXT, Coord71 TEXT, Coord72 TEXT,
Coord73 TEXT, Coord74 TEXT, Coord75 TEXT, Coord76 TEXT, Coord77 TEXT,
Coord78 TEXT, Coord79 TEXT, Coord80 TEXT, Coord81 TEXT, Coord82 TEXT,
Coord83 TEXT, Coord84 TEXT, Coord85 TEXT, Coord86 TEXT, Coord87 TEXT,
Coord88 TEXT, Coord89 TEXT, Coord90 TEXT, Coord91 TEXT, Coord92 TEXT,
Coord93 TEXT, Coord94 TEXT, Coord95 TEXT, Coord96 TEXT)''')
with open('ParametersPosFinalVASP.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into ParametersPosFinalVASP(Energy, Scale,
PrimVec1, PrimVec2, PrimVec3, Coord1, Coord2,
Coord3, Coord4, Coord5, Coord6, Coord7,
Coord8, Coord9, Coord10, Coord11, Coord12,
Coord13, Coord14, Coord15, Coord16, Coord17,
Coord18, Coord19, Coord20, Coord21, Coord22,
Coord23, Coord24, Coord25, Coord26, Coord27,
Coord28, Coord29, Coord30, Coord31, Coord32,
Coord33, Coord34, Coord35, Coord36, Coord37,
Coord38, Coord39, Coord40, Coord41, Coord42,
Coord43, Coord44, Coord45, Coord46, Coord47,
Coord48, Coord49, Coord50, Coord51, Coord52,
Coord53, Coord54, Coord55, Coord56, Coord57,
Coord58, Coord59, Coord60, Coord61, Coord62,
Coord63, Coord64, Coord65, Coord66, Coord67,
Coord68, Coord69, Coord70, Coord71, Coord72,
Coord73, Coord74, Coord75, Coord76, Coord77,
Coord78, Coord79, Coord80, Coord81, Coord82,
Coord83, Coord84, Coord85, Coord86, Coord87,
Coord88, Coord89, Coord90, Coord91, Coord92,
Coord93, Coord94, Coord95, Coord96)
VALUES(?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,
?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,
?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in ParametersPosFinalQE.csv: attachfile:ParametersPosFinalQE.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table ParametersPosFinalQE(Energy FLOAT, Attyps TEXT,
Coord1 TEXT, Coord2 TEXT, Coord3 TEXT, Coord4 TEXT, Coord5 TEXT,
Coord6 TEXT, Coord7 TEXT, Coord8 TEXT, Coord9 TEXT, Coord10 TEXT,
Coord11 TEXT, Coord12 TEXT)''')
with open('ParametersPosFinalQE.csv') as f:
lines = f.readlines() #[0].split('\r')
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into ParametersPosFinalQE(Energy, Attyps,
Coord1, Coord2, Coord3, Coord4, Coord5, Coord6, Coord7, Coord8,
Coord9, Coord10, Coord11, Coord12)
VALUES(?,?,?,?,?,?,?,?,?,?,?,?,?,?)''', fields)
db.commit()
db.close()Data corresponding to this relation is stored in FormEFigure.csv: attachfile:FormEFigure.csv
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
db.execute('''create table FormEFigure(CID TEXT, Figure TEXT)''')
with open('FormEFigure.csv') as f:
lines = f.readlines()
for line in lines[1:]:
fields = [x.strip() for x in line.split(',')]
db.execute('''insert into FormEFigure(CID, Figure)
VALUES(?,?)''', fields)
db.commit()
db.close()In the Supporting Information document of past research cite:Curnan2014, example codes pursuant to the reproduction of all input files needed to complete energetic calculations were created, querying data from the “.sqlite” database file. The input needed to reproduce any calculation, which includes converged atomic position data, as well as input parameters derived from “.in” files in QE and INCAR, KPOINTS, and POTCAR tags in VASP, was reproduced in these examples for VASP calculations. Improvements made to the storage conventions for data and implementation of a different relational schema necessitate that these examples be reproduced in this paper as well. Thus, corresponding example codes are transcribed below for first retrieving the INCAR, KPOINTS, and POTCAR associated with an example calculation or cID, then retrieving its CONTCAR associated information. The example below depicts the calculations required to calculate a Birch-Murnaghan EOS fitted volume and energy, namely calculations that employ a VASP calculator, the PBE functional,
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
NULLchar = '-'
datapts_dict = {}
for row in db.execute('''
select distinct input.ISTART, input.ICHARG, input.ENCUT, input.ISMEAR,
input.SIGMA, input.ISYM, input.IBRION, input.EDIFF, input.EDIFFG, input.MAXMIX,
input.NELMIN, input.NELM, input.NSW, input.ISPIN, input.ISIF, input.GGA,
input.LDAU, input.LDAUU, input.LDAUJ, input.LDAUL, input.LDAUPRINT, input.LASPH,
input.LMAXMIX, input.LORBIT, input.SYMPREC, input.NELMDL, input.LHFCALC,
input.ALGO, input.TIME, input.PRECFOCK, input.LMAXFOCK, input.AEXX, input.NKRED,
input.NBANDS, input.KPoints, input.KPointSampling, input.Pseudopotential1,
input.Pseudopotential2, input.Energy from Structure as s
inner join Symmetry as sy on sy.Struct=s.Struct
inner join Composition1 as c1 on c1.SID=s.SID
inner join Composition2 as c2 on c1.MID=c2.MID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feu.EOpt
inner join ParametersInputVASP as input on input.Energy=cre.Energy
where s.Morph='Rutile'
and mt.PseudoB='PAW-Ti-pv'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and feu.UValue='3.00'
;'''):
datapts_list = []
( a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,
u,v,w,x,y,aa,ab,ac,ad,ae,af,ag,ah,ai,aj,ak,al,am,an ) = row
series = (a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,
aa,ab,ac,ad,ae,af,ag,ah,ai,aj,ak,al,am,an)
for z in series:
datapts_list.append(z)
datapts_dict[an] = datapts_list
del datapts_dict[an][-1]
Input_check = {}
Input_match = {}
att_list = [ 'ISTART', 'ICHARG', 'ENCUT', 'ISMEAR', 'SIGMA', 'ISYM', 'IBRION',
'EDIFF', 'EDIFFG', 'MAXMIX', 'NELMIN', 'NELM', 'NSW', 'ISPIN', 'ISIF', 'GGA',
'LDAU', 'LDAUU', 'LDAUJ', 'LDAUL', 'LDAUPRINT', 'LASPH', 'LMAXMIX', 'LORBIT',
'SYMPREC', 'NELMDL', 'LHFCALC', 'ALGO', 'TIME', 'PRECFOCK', 'LMAXFOCK', 'AEXX',
'NKRED', 'NBANDS', 'KPoints', 'KPointSampling', 'Pseudopotential1', 'Pseudopotential2' ]
for i,j in enumerate( datapts_dict.keys() ):
if datapts_dict[j] not in Input_check.values():
list_match = []
list_check = []
for k,l in enumerate( att_list ):
list_check.append( datapts_dict[j][k] )
if datapts_dict[j][k] != NULLchar:
list_match.append( (att_list[k], datapts_dict[j][k]) )
Input_check[i] = list_check
Input_match[i] = list_match
SortCAR = { "ISTART": 0, "ICHARG": 1, "ENCUT": 2, "ISMEAR": 3,
"SIGMA": 4, "ISYM": 5, "IBRION": 6, "EDIFF": 7,
"EDIFFG": 8, "MAXMIX": 9, "NELMIN": 10, "NELM": 11,
"NSW": 12, "ISPIN": 13, "ISIF": 14, "GGA": 10,
"LDAU": 11, "LDAUU": 12, "LDAUJ": 13, "LDAUL": 14,
"LDAUPRINT": 15, "LASPH": 16, "LMAXMIX": 17,
"LORBIT": 17, "SYMPREC": 18, "NELMDL": 19,
"KPointSampling": 20, "KPoints": 21,
"Pseudopotential1": 22, "Pseudopotential2": 23 }
Input_INCKPTPOT = sorted(Input_match[0], key=lambda tag: SortCAR[tag[0]])
num_KPT = 2
num_POT = 2
num_INCAR = len(Input_INCKPTPOT) - num_KPT - num_POT
for i in range( num_INCAR ):
print Input_INCKPTPOT[i][0],'=',Input_INCKPTPOT[i][1]
print NULLchar
print "0"
print Input_INCKPTPOT[num_INCAR+1][0]
print str(Input_INCKPTPOT[num_INCAR+1][1]).replace(';',' ')
print "0 0 0"
print NULLchar
for i in range( num_POT ):
print Input_INCKPTPOT[num_INCAR+num_KPT+i][0],Input_INCKPTPOT[num_INCAR+num_KPT+i][1]
The script below generates the ordered CONTCAR files used to generate the fitted energy of the system described above, which is derived from the Birch-Murnaghan equation of state fitting procedure and is used to generate plotted data associated with cID values. When available in the database, these CONTCAR input coordinates form structurally relaxed systems that can be outputted to POSCAR files.
import sqlite3
import re
Title = 'Ti O'
Attyps = '2 4'
db = sqlite3.connect('ITEOO_data.sqlite')
NULLchar = '-'
datapts_dict = {}
for row in db.execute('''
select pos.Scale, pos.PrimVec1, pos.PrimVec2, pos.PrimVec3,
pos.Coord1, pos.Coord2, pos.Coord3, pos.Coord4, pos.Coord5,
pos.Coord6, pos.Coord7, pos.Coord8, pos.Coord9, pos.Coord10,
pos.Coord11, pos.Coord12, pos.Coord13, pos.Coord14, pos.Coord15,
pos.Coord16, pos.Coord17, pos.Coord18, pos.Coord19, pos.Coord20,
pos.Coord21, pos.Coord22, pos.Coord23, pos.Coord24, pos.Coord25,
pos.Coord26, pos.Coord27, pos.Coord28, pos.Coord29, pos.Coord30,
pos.Coord31, pos.Coord32, pos.Coord33, pos.Coord34, pos.Coord35,
pos.Coord36, pos.Coord37, pos.Coord38, pos.Coord39, pos.Coord40,
pos.Coord41, pos.Coord42, pos.Coord43, pos.Coord44, pos.Coord45,
pos.Coord46, pos.Coord47, pos.Coord48, pos.Coord49, pos.Coord50,
pos.Coord51, pos.Coord52, pos.Coord53, pos.Coord54, pos.Coord55,
pos.Coord56, pos.Coord57, pos.Coord58, pos.Coord59, pos.Coord60,
pos.Coord61, pos.Coord62, pos.Coord63, pos.Coord64, pos.Coord65,
pos.Coord66, pos.Coord67, pos.Coord68, pos.Coord69, pos.Coord70,
pos.Coord71, pos.Coord72, pos.Coord73, pos.Coord74, pos.Coord75,
pos.Coord76, pos.Coord77, pos.Coord78, pos.Coord79, pos.Coord80,
pos.Coord81, pos.Coord82, pos.Coord83, pos.Coord84, pos.Coord85,
pos.Coord86, pos.Coord87, pos.Coord88, pos.Coord89, pos.Coord90,
pos.Coord91, pos.Coord92, pos.Coord93, pos.Coord94, pos.Coord95,
pos.Coord96, pos.Energy from Structure as s
inner join Symmetry as sy on sy.Struct=s.Struct
inner join Composition1 as c1 on c1.SID=s.SID
inner join Composition2 as c2 on c1.MID=c2.MID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feu.EOpt
inner join ParametersPosFinalVASP as pos on pos.Energy=cre.Energy
where s.Morph='Rutile'
and mt.PseudoB='PAW-Ti-pv'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and feu.UValue='3.00'
;'''):
datapts_list = []
(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,
aa,ab,ac,ad,ae,af,ag,ah,ai,aj,ak,al,am,an,ao,ap,aq,
ar,at,au,av,aw,ax,ay,
ba,bb,bc,bd,be,bf,bg,bh,bi,bj,bk,bl,bm,bn,bo,bp,bq,
br,bs,bt,bu,bv,bw,bx,by,
ca,cb,cc,cd,ce,cf,cg,ch,ci,cj,ck,cl,cm,cn,co,cp,cq,
cr,ct,cu,cv,cw,cx,cy,
da,db,dc) = row
series = ( a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,
aa,ab,ac,ad,ae,af,ag,ah,ai,aj,ak,al,am,an,ao,ap,aq,
ar,at,au,av,aw,ax,ay,
ba,bb,bc,bd,be,bf,bg,bh,bi,bj,bk,bl,bm,bn,bo,bp,bq,
br,bs,bt,bu,bv,bw,bx,by,
ca,cb,cc,cd,ce,cf,cg,ch,ci,cj,ck,cl,cm,cn,co,cp,cq,
cr,ct,cu,cv,cw,cx,cy,
da,db,dc)
for z in series:
datapts_list.append(z)
datapts_dict[a] = datapts_list
del datapts_dict[a][-1]
Input_check = {}
Input_match = {}
att_list = [ 'Scale', 'PrimVec1', 'PrimVec2', 'PrimVec3',
'Coord1', 'Coord2', 'Coord3', 'Coord4', 'Coord5',
'Coord6', 'Coord7', 'Coord8', 'Coord9', 'Coord10',
'Coord11', 'Coord12', 'Coord13', 'Coord14', 'Coord15',
'Coord16', 'Coord17', 'Coord18', 'Coord19', 'Coord20',
'Coord21', 'Coord22', 'Coord23', 'Coord24', 'Coord25',
'Coord26', 'Coord27', 'Coord28', 'Coord29', 'Coord30',
'Coord31', 'Coord32', 'Coord33', 'Coord34', 'Coord35',
'Coord36', 'Coord37', 'Coord38', 'Coord39', 'Coord40',
'Coord41', 'Coord42', 'Coord43', 'Coord44', 'Coord45',
'Coord46', 'Coord47', 'Coord48', 'Coord49', 'Coord50',
'Coord51', 'Coord52', 'Coord53', 'Coord54', 'Coord55',
'Coord56', 'Coord57', 'Coord58', 'Coord59', 'Coord60',
'Coord61', 'Coord62', 'Coord63', 'Coord64', 'Coord65',
'Coord66', 'Coord67', 'Coord68', 'Coord69', 'Coord70',
'Coord71', 'Coord72', 'Coord73', 'Coord74', 'Coord75',
'Coord76', 'Coord77', 'Coord78', 'Coord79', 'Coord80',
'Coord81', 'Coord82', 'Coord83', 'Coord84', 'Coord85',
'Coord86', 'Coord87', 'Coord88', 'Coord89', 'Coord90',
'Coord91', 'Coord92', 'Coord93', 'Coord94', 'Coord95',
'Coord96' ]
SortCAR = { "Scale": 0, "PrimVec1": 1, "PrimVec2": 2, "PrimVec3": 3,
"Coord1": 4, "Coord2": 5, "Coord3": 6, "Coord4": 7,
"Coord5": 8, "Coord6": 9 }
for i,j in enumerate( datapts_dict.keys() ):
list_match = []
for k,l in enumerate( att_list ):
if datapts_dict[j][k] != NULLchar:
list_match.append( (att_list[k], str(datapts_dict[j][k]) ) )
Input_match[datapts_dict[j][0]] = list_match
Input_match[datapts_dict[j][0]].sort(key=lambda val: SortCAR[val[0]])
Input_POSCAR = sorted(Input_match.iteritems(), key=lambda (k,v): k)
Scales_POSCAR = len( Input_match.keys() )
Lines_POSCAR = len( SortCAR )
for i in range( Scales_POSCAR ):
print Title
for j in range( Lines_POSCAR ):
linePOSCAR = re.sub(';', ' ', Input_POSCAR[i][1][j][1].rstrip())
if Input_POSCAR[i][1][j][0] == 'Coord1':
print Attyps
print linePOSCAR
else:
print linePOSCAR
print NULLchar
Supporting information pursuant to explaining the results achieved in this article can largely be divided into two sections, namely sections that primarily explain information pertinent to Hubbard
In calculations featuring the incrementation of Hubbard
In VASP, only input file information corresponding to (final) fixed cell volume, variable cell shape, and atomic coordinate structural relaxation calculations (e.g.: ISIF = 4) are stored in the database, as this is the information that is necessary to reproduce results demonstrated in this article. However, the steps implemented to perform structural relaxations in this study are detailed below, with sequentially ordered input files listed to illustrate the process through which structural relaxations were performed and pertinent energetics were resolved. For an example calculation of Rutile completed using standard PAW pseudopotentials and the PBE functional, calculations employing ISIF = 4 were completed over ranges of volumes encompassing an equilibrium volume determined from experiment cite:Muscat2002 at
Energetic convergence is evaluated iteratively between distinct sets of freedom in each structure, while convergence is determined using a volumetric criterion. Thus, for Rutile, ISIF = 4 relaxations are first performed over the
import math
from decimal import *
V_it = []
dof_num = 2
for i in (1 1 N):
Energies = []
Volumes = []
for v in (95 S 105):
Energies.append( get.Energy(v) )
Volumes.append( get.Volume(v) )
E_opt, V_opt = BMurn( Energies, Volumes )
V_it.append( V_opt )
if i =< dof_num:
pass
else:
S_float = math.log10( float(S) ) / 100.
if Decimal( V_it[i] ).quantize( Decimal(S_float) ) == Decimal( V_it[i - dof_num] ).quantize(
Decimal(S_float) ):
E_final = E_opt
V_final = V_opt
break
print E_final, V_finalThe iterative cycle above, which illustrates the aggregation of data that has already been calculated within a fixed range of volumes, requires that several variables be defined or explained. Firstly, “N” represents the total number of iterations accomplished by the structural relaxation, while “i” represents the current iteration being performed by the loop. However, if the pseudocode above was adapted to execute calculations as well, the “for” loop above would be substituted with a “while” loop, “i” would become a counter for the loop itself, and “N” would not be defined. The volume of a current system, or at least the measure of a structural parameter monotonically related to volume over which ISIF = 4 calculations are performed (e.g.:
Using an additive scale as opposed to a multiplicative scale (shown above), as well as a
ISTART = 0 ; ICHARG = 2
ENCUT = 600
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 1
IBRION = 1
EDIFF = 5E-06; EDIFFG = -0.01
MAXMIX = -100
NELMIN = 5
NELMDL = -10
NELM = 200
NSW = 100
ISPIN = 2
ISIF = 4
8x8x8
0
Monkhorst
8 8 8
0 0 0
The first structural iteration (it1) takes its
4.66
4.66
1 0 0
0 1 0
0 0 0.643993896
2 4
direct
0.000000000 0.000000000 0.000000000
0.500000000 0.500000000 0.500000000
0.304880731 0.304880731 0.000000000
-0.304880731 -0.304880731 0.000000000
0.804875587 0.195124413 0.500000000
-0.804875587 -0.195124413 0.500000000
The second iteration (it2) varies the
4.65922056943
4.65922056943
1 0 0
0 1 0
0 0 0.64
2 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.3048134906850396 0.3048134906850396 0.0000000000000000
0.6951865093149605 0.6951865093149605 0.0000000000000000
0.8048075306488891 0.1951924693511108 0.5000000000000000
0.1951924693511108 0.8048075306488891 0.5000000000000000
The final iteration takes the finalized atomic coordinates and fitted
4.66
4.66
1 0 0
0 1 0
0 0 0.638740876486
2 4
direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.3048082437656310 0.3048082437656310 0.0000000000000000
0.6951917562343689 0.6951917562343689 0.0000000000000000
0.8048024643031844 0.1951975356968156 0.5000000000000000
0.1951975356968156 0.8048024643031844 0.5000000000000000
A calculation employing ISIF = 3 that uses the atomic coordinates of the system tested in the previous step closest to the energetic minimum is completed. This calculation uses the following INCAR file and input file information that is otherwise comparable to the first step:
ISTART = 0 ; ICHARG = 2
ENCUT = 600
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 1
IBRION = 1
EDIFF = 5E-06; EDIFFG = -0.01
MAXMIX = -100
NELMIN = 5
NELMDL = -10
NELM = 200
NSW = 100
ISPIN = 2
ISIF = 3
4.66222984355
4.66222984355
1 0 0
0 1 0
0 0 0.6375189084585891
2 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.3047756242591685 0.3047756242591685 0.0000000000000000
0.6952243757408315 0.6952243757408315 0.0000000000000000
0.8047805609392189 0.1952194390607812 0.5000000000000000
0.1952194390607812 0.8047805609392189 0.5000000000000000
The structures resolved from these ISIF = 3 calculations are applied to Hubbard
ISTART = 0 ; ICHARG = 2
ENCUT = 600
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 1; SYMPREC = 1E-06
IBRION = 1
EDIFF = 5E-06; EDIFFG = -0.01
MAXMIX = -100
NELMIN = 5
NELMDL = -10
NELM = 200
NSW = 100
ISPIN = 2
ISIF = 4
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1
LDAUPRINT = 1
LASPH = .TRUE.
LMAXMIX = 4
LDAUJ = 0.00 0.00
LDAUU = 1.00 0.00
&CONTROL
calculation = "relax",
verbosity = "low",
restart_mode = "from_scratch",
pseudo_dir = "../../",
outdir = "./",
prefix = "BO2_relax",
nstep = 100,
wf_collect = .false.
forc_conv_thr = 1.0E-3
etot_conv_thr = 1.0E-6
/
&SYSTEM
ibrav = 6
nat = 6
ntyp = 2
ecutwfc = 50.0
ecutrho = 600.0
nspin = 2
starting_magnetization(1) = 0.0
occupations = "smearing",
smearing = "gaussian",
degauss = 0.01
lda_plus_u = .true.
lda_plus_u_kind = 0
celldm(1) = 8.85
celldm(3) = 0.637562623236
Hubbard_U(1) = 1.00
/
&ELECTRONS
electron_maxstep = 100
conv_thr = 1D-9
diagonalization = "david",
diago_thr_init = 1D-2
diago_full_acc = .false.
startingwfc = "atomic+random",
mixing_mode = "plain",
mixing_beta = 0.3
/
&IONS
ion_dynamics = "bfgs",
upscale = 1000
/
ATOMIC_SPECIES
Ti 1.0 Ti.pbe-sp-van_ak.UPF
O 1.0 O.pbe-rrkjus.UPF
ATOMIC_POSITIONS (crystal)
Ti 0.000000000 0.000000000 0.000000000
Ti 0.500000000 0.500000000 0.500000000
O 0.305356333 0.305356333 0.000000000
O 0.694643667 0.694643667 0.000000000
O 0.805354772 0.194645228 0.500000000
O 0.194645228 0.805354772 0.500000000
K_POINTS (automatic)
8 8 8 0 0 0
In the article, there are four plots featuring Hubbard
\begin{flalign*}
&ΠPseudoO, Functional, Morph, UValue, EOpt, Stoich1
&σ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) \
&σ ∨ (Morph=’Columbite’) ∨ (Morph=’Brookite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’) ∨ (PseudoO=’PAW-O-s’)] \
&σ ∧ [(Functional=’PBE’) ∨ (Functional=’LDA’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query PBE LDA Ti O O_s
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Pseudo_O, Functional, Polymorph, Uvalue, Eopt, atomnum = [], [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select distinct mt.PseudoO, mt.Functional, s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or s.Morph='Brookite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and (mt.PseudoO='PAW-O' or mt.PseudoO='PAW-O-s')
and (mt.Functional='PBE' or mt.Functional='LDA')
and mt.Software='VASP'
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
datapts_list = []
a,b,c,d,e,f = row
datapts_list.append( (a,b,c,d,e,f) )
Pseudo_O += [a]
Functional += [b]
Polymorph += [c]
datapts_dict[row] = datapts_list
Pseudo_list = list( set( Pseudo_O ) )
Functional_list = list( set( Functional ) )
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Pseudo_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Functional_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = {}
for k in Morph_list:
SortElist[i][j][k] = {}
EBsiteMatch[i][j][k] = []
del_list = []
for l in datapts_dict:
if l[0] == i and l[1] == j and l[2] == k:
SortElist[ l[0] ][ l[1] ][ l[2] ][ l[3] ] = float( l[4] ) / float( l[5] )
del_list.append( datapts_dict[l] )
else:
pass
for l in del_list:
del l
for i in Pseudo_list:
for j in Functional_list:
for k in Morph_list:
U_list = SortElist[i][j][k].keys()
U_list.sort()
EmapU_list = []
for l in U_list:
EmapU_value = SortElist[i][j][k][l] - SortElist[i][j][WRT_Morph][l]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j][k].append( EmapU_list )
EBsiteMatch[i][j][k].append( U_list )
ax = plt.gca()
print "PBE Functional, O pseudopotential:"
print "Rutile: " + str(EBsiteMatch['PAW-O']['PBE']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O']['PBE']['Rutile'][0][2:4])
print str(EBsiteMatch['PAW-O']['PBE']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['PAW-O']['PBE']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O']['PBE']['Anatase'][0][2:4])
print str(EBsiteMatch['PAW-O']['PBE']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['PAW-O']['PBE']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O']['PBE']['Columbite'][0][2:4])
print str(EBsiteMatch['PAW-O']['PBE']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['PAW-O']['PBE']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O']['PBE']['Brookite'][0][2:4])
print str(EBsiteMatch['PAW-O']['PBE']['Brookite'][0][4:]) + "\n"
print "LDA Functional, O pseudopotential:"
print "Rutile: " + str(EBsiteMatch['PAW-O']['LDA']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O']['LDA']['Rutile'][0][2:4])
print str(EBsiteMatch['PAW-O']['LDA']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['PAW-O']['LDA']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O']['LDA']['Anatase'][0][2:4])
print str(EBsiteMatch['PAW-O']['LDA']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['PAW-O']['LDA']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O']['LDA']['Columbite'][0][2:4])
print str(EBsiteMatch['PAW-O']['LDA']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['PAW-O']['LDA']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O']['LDA']['Brookite'][0][2:4])
print str(EBsiteMatch['PAW-O']['LDA']['Brookite'][0][4:]) + "\n"
print "PBE Functional, O_s pseudopotential:"
print "Rutile: " + str(EBsiteMatch['PAW-O-s']['PBE']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['PBE']['Rutile'][0][2:4])
print str(EBsiteMatch['PAW-O-s']['PBE']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['PAW-O-s']['PBE']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['PBE']['Anatase'][0][2:4])
print str(EBsiteMatch['PAW-O-s']['PBE']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['PAW-O-s']['PBE']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['PBE']['Columbite'][0][2:4])
print str(EBsiteMatch['PAW-O-s']['PBE']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['PAW-O-s']['PBE']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['PBE']['Brookite'][0][2:4])
print str(EBsiteMatch['PAW-O-s']['PBE']['Brookite'][0][4:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['PAW-O']['PBE']['Rutile'][1], EBsiteMatch['PAW-O']['PBE']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['PAW-O']['PBE']['Anatase'][1], EBsiteMatch['PAW-O']['PBE']['Anatase'][0],
'bo-', label = r'$\Delta$E$_{R-A}$,PBE')
plt.plot(EBsiteMatch['PAW-O']['PBE']['Columbite'][1], EBsiteMatch['PAW-O']['PBE']['Columbite'][0],
'ro-', label = r'$\Delta$E$_{R-C}$,PBE')
plt.plot(EBsiteMatch['PAW-O']['PBE']['Brookite'][1], EBsiteMatch['PAW-O']['PBE']['Brookite'][0],
'go-', label = r'$\Delta$E$_{R-B}$,PBE')
plt.plot(EBsiteMatch['PAW-O-s']['PBE']['Anatase'][1], EBsiteMatch['PAW-O-s']['PBE']['Anatase'][0],
'bv-', label = r'$\Delta$E$_{R-A}$,PBEs')
plt.plot(EBsiteMatch['PAW-O-s']['PBE']['Columbite'][1], EBsiteMatch['PAW-O-s']['PBE']['Columbite'][0],
'rv-', label = r'$\Delta$E$_{R-C}$,PBEs')
plt.plot(EBsiteMatch['PAW-O-s']['PBE']['Brookite'][1], EBsiteMatch['PAW-O-s']['PBE']['Brookite'][0],
'gv-', label = r'$\Delta$E$_{R-B}$,PBEs')
plt.plot(EBsiteMatch['PAW-O']['LDA']['Anatase'][1], EBsiteMatch['PAW-O']['LDA']['Anatase'][0],
'bs-', label = r'$\Delta$E$_{R-A}$,LDA')
plt.plot(EBsiteMatch['PAW-O']['LDA']['Columbite'][1], EBsiteMatch['PAW-O']['LDA']['Columbite'][0],
'rs-', label = r'$\Delta$E$_{R-C}$,LDA')
plt.plot(EBsiteMatch['PAW-O']['LDA']['Brookite'][1], EBsiteMatch['PAW-O']['LDA']['Brookite'][0],
'gs-', label = r'$\Delta$E$_{R-B}$,LDA')
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 0.3) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper center', prop={'size':6.5}, ncol = 2)
plt.axvspan(2.79, 4.3, facecolor='m', alpha=0.5)
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-RACB-PBELDAPBEs' + '.' + ext, dpi=300)
plt.clf()
The second plot features solely the PBE functional, albeit pseudopotential sets in this plot feature
\begin{flalign*}
&ΠPseudoO, Morph, UValue, EOpt, Stoich1
&σ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) \
&σ ∨ (Morph=’Columbite’) ∨ (Morph=’Brookite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti-pv’)] \
&σ ∧ [(PseudoO=’PAW-O’) ∨ (PseudoO=’PAW-O-s’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’) ∨ (UValue=’7’) \
&σ ∨ (UValue=’8’) ∨ (UValue=’9’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query PBE Ti_pv O O_s
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Pseudo_O, Polymorph, Uvalue, Eopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select distinct mt.PseudoO, s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or s.Morph='Brookite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti-pv'
and (mt.PseudoO='PAW-O' or mt.PseudoO='PAW-O-s')
and mt.Functional='PBE'
and mt.Software='VASP'
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6' or feu.UValue='7'
or feu.UValue='8' or feu.UValue='9')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
Pseudo_O += [a]
Polymorph += [b]
datapts_dict[row] = datapts_list
Pseudo_list = list( set( Pseudo_O ) )
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Pseudo_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Morph_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in Pseudo_list:
for j in Morph_list:
U_list = SortElist[i][j].keys()
U_list.sort()
EmapU_list = []
for k in U_list:
EmapU_value = SortElist[i][j][k] - SortElist[i][WRT_Morph][k]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j].append( EmapU_list )
EBsiteMatch[i][j].append( U_list )
ax = plt.gca()
print "PBE Functional, Ti_pv O pseudopotentials:"
print "Rutile: " + str(EBsiteMatch['PAW-O']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O']['Rutile'][0][2:5])
print str(EBsiteMatch['PAW-O']['Rutile'][0][5:8])
print str(EBsiteMatch['PAW-O']['Rutile'][0][8:])
print "Anatase: " + str(EBsiteMatch['PAW-O']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O']['Anatase'][0][2:5])
print str(EBsiteMatch['PAW-O']['Anatase'][0][5:8])
print str(EBsiteMatch['PAW-O']['Anatase'][0][8:])
print "Columbite: " + str(EBsiteMatch['PAW-O']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O']['Columbite'][0][2:5])
print str(EBsiteMatch['PAW-O']['Columbite'][0][5:8])
print str(EBsiteMatch['PAW-O']['Columbite'][0][8:])
print "Brookite: " + str(EBsiteMatch['PAW-O']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O']['Brookite'][0][2:5])
print str(EBsiteMatch['PAW-O']['Brookite'][0][5:8])
print str(EBsiteMatch['PAW-O']['Brookite'][0][8:]) + "\n"
print "PBE Functional, Ti_pv O_s pseudopotentials:"
print "Rutile: " + str(EBsiteMatch['PAW-O-s']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Rutile'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Rutile'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Rutile'][0][8:])
print "Anatase: " + str(EBsiteMatch['PAW-O-s']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Anatase'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Anatase'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Anatase'][0][8:])
print "Columbite: " + str(EBsiteMatch['PAW-O-s']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Columbite'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Columbite'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Columbite'][0][8:])
print "Brookite: " + str(EBsiteMatch['PAW-O-s']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Brookite'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Brookite'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Brookite'][0][8:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['PAW-O']['Rutile'][1], EBsiteMatch['PAW-O']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['PAW-O']['Anatase'][1], EBsiteMatch['PAW-O']['Anatase'][0],
'bo-', label = r'$\Delta$E$_{R-A}$,PBEpv')
plt.plot(EBsiteMatch['PAW-O']['Columbite'][1], EBsiteMatch['PAW-O']['Columbite'][0],
'ro-', label = r'$\Delta$E$_{R-C}$,PBEpv')
plt.plot(EBsiteMatch['PAW-O']['Brookite'][1], EBsiteMatch['PAW-O']['Brookite'][0],
'go-', label = r'$\Delta$E$_{R-B}$,PBEpv')
plt.plot(EBsiteMatch['PAW-O-s']['Anatase'][1], EBsiteMatch['PAW-O-s']['Anatase'][0],
'bs-', label = r'$\Delta$E$_{R-A}$,PBEpvs')
plt.plot(EBsiteMatch['PAW-O-s']['Columbite'][1], EBsiteMatch['PAW-O-s']['Columbite'][0],
'rs-', label = r'$\Delta$E$_{R-C}$,PBEpvs')
plt.plot(EBsiteMatch['PAW-O-s']['Brookite'][1], EBsiteMatch['PAW-O-s']['Brookite'][0],
'gs-', label = r'$\Delta$E$_{R-B}$,PBEpvs')
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 0.15) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper center', prop={'size':6}, ncol = 2)
plt.axvspan(4.74, 7.0, facecolor='m', alpha=0.5)
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-RACB-pvpvs' + '.' + ext, dpi=300)
plt.clf()
The third plot features solely the PBE functional, albeit pseudopotential sets in this plot feature
\begin{flalign*}
&ΠPseudoO, Morph, UValue, EOpt, Stoich1
&σ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) \
&σ ∨ (Morph=’Columbite’) ∨ (Morph=’Brookite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti-sv’)] \
&σ ∧ [(PseudoO=’PAW-O’) ∨ (PseudoO=’PAW-O-s’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’) ∨ (UValue=’7’) \
&σ ∨ (UValue=’8’) ∨ (UValue=’9’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query PBE Ti_sv O O_s
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Pseudo_O, Polymorph, Uvalue, Eopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select distinct mt.PseudoO, s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or s.Morph='Brookite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti-sv'
and (mt.PseudoO='PAW-O' or mt.PseudoO='PAW-O-s')
and mt.Functional='PBE'
and mt.Software='VASP'
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6' or feu.UValue='7'
or feu.UValue='8' or feu.UValue='9')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
Pseudo_O += [a]
Polymorph += [b]
datapts_dict[row] = datapts_list
Pseudo_list = list( set( Pseudo_O ) )
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Pseudo_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Morph_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in Pseudo_list:
for j in Morph_list:
U_list = SortElist[i][j].keys()
U_list.sort()
EmapU_list = []
for k in U_list:
EmapU_value = SortElist[i][j][k] - SortElist[i][WRT_Morph][k]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j].append( EmapU_list )
EBsiteMatch[i][j].append( U_list )
ax = plt.gca()
print "PBE Functional, Ti_sv O pseudopotentials:"
print "Rutile: " + str(EBsiteMatch['PAW-O']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O']['Rutile'][0][2:5])
print str(EBsiteMatch['PAW-O']['Rutile'][0][5:8])
print str(EBsiteMatch['PAW-O']['Rutile'][0][8:])
print "Anatase: " + str(EBsiteMatch['PAW-O']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O']['Anatase'][0][2:5])
print str(EBsiteMatch['PAW-O']['Anatase'][0][5:8])
print str(EBsiteMatch['PAW-O']['Anatase'][0][8:])
print "Columbite: " + str(EBsiteMatch['PAW-O']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O']['Columbite'][0][2:5])
print str(EBsiteMatch['PAW-O']['Columbite'][0][5:8])
print str(EBsiteMatch['PAW-O']['Columbite'][0][8:])
print "Brookite: " + str(EBsiteMatch['PAW-O']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O']['Brookite'][0][2:5])
print str(EBsiteMatch['PAW-O']['Brookite'][0][5:8])
print str(EBsiteMatch['PAW-O']['Brookite'][0][8:]) + "\n"
print "PBE Functional, Ti_sv O_s pseudopotentials:"
print "Rutile: " + str(EBsiteMatch['PAW-O-s']['Rutile'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Rutile'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Rutile'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Rutile'][0][8:])
print "Anatase: " + str(EBsiteMatch['PAW-O-s']['Anatase'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Anatase'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Anatase'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Anatase'][0][8:])
print "Columbite: " + str(EBsiteMatch['PAW-O-s']['Columbite'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Columbite'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Columbite'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Columbite'][0][8:])
print "Brookite: " + str(EBsiteMatch['PAW-O-s']['Brookite'][0][0:2])
print str(EBsiteMatch['PAW-O-s']['Brookite'][0][2:5])
print str(EBsiteMatch['PAW-O-s']['Brookite'][0][5:8])
print str(EBsiteMatch['PAW-O-s']['Brookite'][0][8:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['PAW-O']['Rutile'][1], EBsiteMatch['PAW-O']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['PAW-O']['Anatase'][1], EBsiteMatch['PAW-O']['Anatase'][0],
'bo-', label = r'$\Delta$E$_{R-A}$,PBEsv')
plt.plot(EBsiteMatch['PAW-O']['Columbite'][1], EBsiteMatch['PAW-O']['Columbite'][0],
'ro-', label = r'$\Delta$E$_{R-C}$,PBEsv')
plt.plot(EBsiteMatch['PAW-O']['Brookite'][1], EBsiteMatch['PAW-O']['Brookite'][0],
'go-', label = r'$\Delta$E$_{R-B}$,PBEsv')
plt.plot(EBsiteMatch['PAW-O-s']['Anatase'][1], EBsiteMatch['PAW-O-s']['Anatase'][0],
'bs-', label = r'$\Delta$E$_{R-A}$,PBEsvs')
plt.plot(EBsiteMatch['PAW-O-s']['Columbite'][1], EBsiteMatch['PAW-O-s']['Columbite'][0],
'rs-', label = r'$\Delta$E$_{R-C}$,PBEsvs')
plt.plot(EBsiteMatch['PAW-O-s']['Brookite'][1], EBsiteMatch['PAW-O-s']['Brookite'][0],
'gs-', label = r'$\Delta$E$_{R-B}$,PBEsvs')
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 0.1) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper center', prop={'size':6}, ncol = 2)
plt.axvspan(5.8, 8.2, facecolor='m', alpha=0.5)
plt.gcf().subplots_adjust(left=0.27)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-RACB-svsvs' + '.' + ext, dpi=300)
plt.clf()
The fourth plot features solely standard pseudpotentials, albeit calculation in this plot use the PBEsol, PW91, and AM05 functionals. The relational algebraic expression and complementary MySQL query used to generate these plots are transcribed below:
\begin{flalign*}
&ΠFunctional, Morph, UValue, EOpt, Stoich1
&σ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) \
&σ ∨ (Morph=’Columbite’) ∨ (Morph=’Brookite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Functional=’PS’) ∨ (Functional=’PW91’) ∨ (Functional=’AM05’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query PBEsol PW91 AM05 Ti O
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Functional, Polymorph, Uvalue, Eopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select mt.Functional, s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or s.Morph='Brookite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and (mt.Functional='PS' or mt.Functional='PW91' or mt.Functional='AM05')
and mt.Software='VASP'
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
Functional += [a]
Polymorph += [b]
datapts_dict[row] = datapts_list
Functional_list = list( set( Functional ) )
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Functional_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Morph_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in Functional_list:
for j in Morph_list:
U_list = SortElist[i][j].keys()
U_list.sort()
EmapU_list = []
for k in U_list:
EmapU_value = SortElist[i][j][k] - SortElist[i][WRT_Morph][k]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j].append( EmapU_list )
EBsiteMatch[i][j].append( U_list )
ax = plt.gca()
print "PBEsol Functional:"
print "Rutile: " + str(EBsiteMatch['PS']['Rutile'][0][0:2])
print str(EBsiteMatch['PS']['Rutile'][0][2:4])
print str(EBsiteMatch['PS']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['PS']['Anatase'][0][0:2])
print str(EBsiteMatch['PS']['Anatase'][0][2:4])
print str(EBsiteMatch['PS']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['PS']['Columbite'][0][0:2])
print str(EBsiteMatch['PS']['Columbite'][0][2:4])
print str(EBsiteMatch['PS']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['PS']['Brookite'][0][0:2])
print str(EBsiteMatch['PS']['Brookite'][0][2:4])
print str(EBsiteMatch['PS']['Brookite'][0][4:]) + "\n"
print "PW91 Functional:"
print "Rutile: " + str(EBsiteMatch['PW91']['Rutile'][0][0:2])
print str(EBsiteMatch['PW91']['Rutile'][0][2:4])
print str(EBsiteMatch['PW91']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['PW91']['Anatase'][0][0:2])
print str(EBsiteMatch['PW91']['Anatase'][0][2:4])
print str(EBsiteMatch['PW91']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['PW91']['Columbite'][0][0:2])
print str(EBsiteMatch['PW91']['Columbite'][0][2:4])
print str(EBsiteMatch['PW91']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['PW91']['Brookite'][0][0:2])
print str(EBsiteMatch['PW91']['Brookite'][0][2:4])
print str(EBsiteMatch['PW91']['Brookite'][0][4:]) + "\n"
print "AM05 Functional:"
print "Rutile: " + str(EBsiteMatch['AM05']['Rutile'][0][0:2])
print str(EBsiteMatch['AM05']['Rutile'][0][2:4])
print str(EBsiteMatch['AM05']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['AM05']['Anatase'][0][0:2])
print str(EBsiteMatch['AM05']['Anatase'][0][2:4])
print str(EBsiteMatch['AM05']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['AM05']['Columbite'][0][0:2])
print str(EBsiteMatch['AM05']['Columbite'][0][2:4])
print str(EBsiteMatch['AM05']['Columbite'][0][4:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['PS']['Rutile'][1], EBsiteMatch['PS']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['PS']['Anatase'][1], EBsiteMatch['PS']['Anatase'][0],
'bo-', label = r'$\Delta$E$_{R-A}$,PBEsol')
plt.plot(EBsiteMatch['PS']['Columbite'][1], EBsiteMatch['PS']['Columbite'][0],
'ro-', label = r'$\Delta$E$_{R-C}$,PBEsol')
plt.plot(EBsiteMatch['PS']['Brookite'][1], EBsiteMatch['PS']['Brookite'][0],
'go-', label = r'$\Delta$E$_{R-B}$,PBEsol')
plt.plot(EBsiteMatch['PW91']['Anatase'][1], EBsiteMatch['PW91']['Anatase'][0],
'bs-', label = r'$\Delta$E$_{R-A}$,PW91')
plt.plot(EBsiteMatch['PW91']['Columbite'][1], EBsiteMatch['PW91']['Columbite'][0],
'rs-', label = r'$\Delta$E$_{R-C}$,PW91')
plt.plot(EBsiteMatch['PW91']['Brookite'][1], EBsiteMatch['PW91']['Brookite'][0],
'gs-', label = r'$\Delta$E$_{R-B}$,PW91')
plt.plot(EBsiteMatch['AM05']['Anatase'][1], EBsiteMatch['AM05']['Anatase'][0],
'bv-', label = r'$\Delta$E$_{R-A}$,AM05')
plt.plot(EBsiteMatch['AM05']['Columbite'][1], EBsiteMatch['AM05']['Columbite'][0],
'rv-', label = r'$\Delta$E$_{R-C}$,AM05')
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 0.2) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper center', prop={'size':6}, ncol = 2)
plt.axvspan(2.79, 4.0, facecolor='m', alpha=0.5)
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-RACB-PSPW91AM05' + '.' + ext, dpi=300)
plt.clf()
Upon modification of the above queries in the org-mode source file introduced in Section Introduction of this Supporting Information document, the four example queries shown above can serve as templates through which any combination of BO2 polymorphic trend data featuring
The Columbite polymorph can alternatively be classified via the α-PbO2 structure (
\begin{flalign*}
&ΠMorph, UValue, EOpt, NumAtoms
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Columbite’) \
&σ ∨ (Morph=’Baddeleyite’) ∨ (Morph=’Cotunnite’) \
&σ ∨ (Morph=’Fluorite’) ∨ (Morph=’Pyrite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query high pressure polymorph
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, Uvalue, Eopt, atomnum = [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Columbite' or s.Morph='Baddeleyite' or
s.Morph='Cotunnite' or s.Morph='Fluorite' or s.Morph='Pyrite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Functional='PBE'
and mt.Software='VASP'
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
Polymorph += [a]
datapts_dict[row] = datapts_list
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Morph_list:
SortElist[i] = {}
EBsiteMatch[i] = []
del_list = []
for j in datapts_dict:
if j[0] == i:
SortElist[ j[0] ][ j[1] ] = float( j[2] ) / float( j[3] )
del_list.append( datapts_dict[j] )
else:
pass
for l in del_list:
del l
for i in Morph_list:
U_list = SortElist[i].keys()
U_list.sort()
EmapU_list = []
for j in U_list:
EmapU_value = SortElist[i][j] - SortElist[WRT_Morph][j]
EmapU_list.append( EmapU_value )
EBsiteMatch[i].append( EmapU_list )
EBsiteMatch[i].append( U_list )
ax = plt.gca()
print "PBE Functional, Ti O pseudopotentials:"
print "Rutile: " + str(EBsiteMatch['Rutile'][0][0:2])
print str(EBsiteMatch['Rutile'][0][2:4])
print str(EBsiteMatch['Rutile'][0][4:])
print "Columbite: " + str(EBsiteMatch['Columbite'][0][0:2])
print str(EBsiteMatch['Columbite'][0][2:4])
print str(EBsiteMatch['Columbite'][0][4:])
print "Baddeleyite: " + str(EBsiteMatch['Baddeleyite'][0][0:2])
print str(EBsiteMatch['Baddeleyite'][0][2:4])
print str(EBsiteMatch['Baddeleyite'][0][4:])
print "Cotunnite: " + str(EBsiteMatch['Cotunnite'][0][0:2])
print str(EBsiteMatch['Cotunnite'][0][2:4])
print str(EBsiteMatch['Cotunnite'][0][4:])
print "Fluorite: " + str(EBsiteMatch['Fluorite'][0][0:2])
print str(EBsiteMatch['Fluorite'][0][2:4])
print str(EBsiteMatch['Fluorite'][0][4:])
print "Pyrite: " + str(EBsiteMatch['Pyrite'][0][0:2])
print str(EBsiteMatch['Pyrite'][0][2:4])
print str(EBsiteMatch['Pyrite'][0][4:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Rutile'][1], EBsiteMatch['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['Columbite'][1], EBsiteMatch['Columbite'][0],
'bo-', label = r'$\Delta$E$_{R-C}$')
plt.plot(EBsiteMatch['Baddeleyite'][1], EBsiteMatch['Baddeleyite'][0],
'ro-', label = r'$\Delta$E$_{R-Ba}$')
plt.plot(EBsiteMatch['Cotunnite'][1], EBsiteMatch['Cotunnite'][0],
'go-', label = r'$\Delta$E$_{R-Co}$')
plt.plot(EBsiteMatch['Fluorite'][1], EBsiteMatch['Fluorite'][0],
'co-', label = r'$\Delta$E$_{R-F}$')
plt.plot(EBsiteMatch['Pyrite'][1], EBsiteMatch['Pyrite'][0],
'mo-', label = r'$\Delta$E$_{R-P}$')
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 1.2) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper center', prop={'size':9}, ncol = 2)
plt.gcf().subplots_adjust(left=0.20)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-RCBaCoFP-PBE' + '.' + ext, dpi=300)
plt.clf()
In the plot shown above, note that calculated results for the Baddeleyite polymorph are only physically representative of experimental structures up to a
In TiO2, the Baddeleyite polymorph observes experimentally resolved lattice vectors of
\begin{flalign*}
&ΠEnergy, LDAUU, Eopt, PrimVec1, PrimVec2
&σ ∧ [(Morph=’Baddeleyite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange \
& \Join CalcResultEnergetics \Join ParametersInputVASP \Join ParametersPosFinalVASP]
\end{flalign*}
#+caption: MySQL query high pressure Baddeleyite criterion
import matplotlib.pyplot as plt
import sqlite3
import numpy as np
db = sqlite3.connect('ITEOO_data.sqlite')
U_vals = []
EtoEOptU_dict = {}
NULL_CHAR = '-'
SColon_CHAR = ';'
WRT_U = 0.
Morph_list = ['Baddeleyite']
for row in db.execute('''
select distinct cre.Energy, cre.EOpt from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Composition2 as c2 on c1.MID=c2.MID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feu.EOpt
where (s.Morph='Baddeleyite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and (mt.Functional='PBE')
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
a,b = row
try:
EtoEOptU_dict[b].append(a)
except KeyError:
EtoEOptU_dict[b] = []
EtoEOptU_dict[b].append(a)
Emin_list = []
for i in EtoEOptU_dict.keys():
Ediff_list = []
for j in EtoEOptU_dict[i]:
Ediff_value = abs( float(i) - float(j) )
Ediff_list.append( Ediff_value )
Emindex = Ediff_list.index( min( Ediff_list ) )
Emin_list.append( EtoEOptU_dict[i][ Emindex ] )
datapts_dict = {}
for row in db.execute('''
select distinct input.Energy, input.LDAUU, pos.PrimVec1, pos.PrimVec2 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Composition2 as c2 on c1.MID=c2.MID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feu.EOpt
inner join ParametersPosFinalVASP as pos on pos.Energy=cre.Energy
inner join ParametersInputVASP as input on input.Energy=pos.Energy
where (s.Morph='Baddeleyite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and (mt.Functional='PBE')
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
U_vals += [b]
if a in Emin_list:
datapts_dict[row] = datapts_list
U_vals_list = list( set( U_vals ) )
SortElist = {}
EBsiteMatch = {}
Ulabel_list = []
for i in Morph_list:
SortElist[i] = {}
EBsiteMatch[i] = []
for j in U_vals_list:
Utype_val = str(j).split(SColon_CHAR)[0]
Ulabel_list.append( Utype_val )
SortElist[i][ Utype_val ] = {}
del_list = []
for k in datapts_dict:
if k[1] == j:
a_string = str( k[2] ).split(SColon_CHAR)
c_string = str( k[3] ).split(SColon_CHAR)
a_vec = np.zeros( len(a_string) )
c_vec = np.zeros( len(c_string) )
for l in range( len(a_string) ):
a_vec[l] = float( a_string[l] )
c_vec[l] = float( c_string[l] )
a_dot = np.dot( a_vec, a_vec )
c_dot = np.dot( c_vec, c_vec )
ca_ratio = ( c_dot / a_dot )**(0.5)
SortElist[ i ][ Utype_val ] = ca_ratio
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in Morph_list:
U_list = SortElist[i].keys()
camapU_list = []
U_list.sort()
for k in U_list:
camapU_value = SortElist[i][k]
camapU_list.append( camapU_value )
EBsiteMatch[i].append( camapU_list )
EBsiteMatch[i].append( U_list )
print "U values: " + str(EBsiteMatch['Baddeleyite'][1])
print "b/a ratios: " + str(EBsiteMatch['Baddeleyite'][0][0:3])
print str(EBsiteMatch['Baddeleyite'][0][3:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Baddeleyite'][1], EBsiteMatch['Baddeleyite'][0], 'bo-', label = r'$\Delta$b/a(Baddeleyite)')
plt.xlabel( 'U value (eV)' )
plt.ylabel('b/a ratio')
plt.legend(loc = 'lower left', prop={'size':9})
plt.gcf().subplots_adjust(left=0.25)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-baratio-Ba-PBE' + '.' + ext, dpi=300)
plt.clf()
Additionally, the Quantum Espresso software package is used to calculate first-principles
\begin{flalign*}
&ΠMorph, Software, UValue, EOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Columbite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’) ∨ (PseudoB=’US-Ti-pv-sv’)] \
&σ ∧ [(PseudoO=’PAW-O’) ∨ (PseudoO=’US-O’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’) ∨ (Software=’QE’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query VASP vs. QE Rutile Columbite
import matplotlib.pyplot as plt
import sqlite3
convRytoeV = 13.605698066
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, Calculator, Uvalue, Eopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select s.Morph, mt.Software, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on ( mt.MID=c1.MID and mt.SID=c1.SID )
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Columbite')
and c1.At1='Ti'
and (mt.PseudoB='PAW-Ti' or mt.PseudoB='US-Ti-pv-sv')
and (mt.PseudoO='PAW-O' or mt.PseudoO='US-O')
and mt.Functional='PBE'
and (mt.Software='VASP' or mt.Software='QE')
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
Polymorph += [a]
Calculator += [b]
datapts_dict[row] = datapts_list
Morph_list = list( set( Polymorph ) )
Calculator_list = list( set( Calculator ) )
SortElist = {}
EBsiteMatch = {}
for i in Morph_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Calculator_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in Morph_list:
for j in Calculator_list:
U_list = SortElist[i][j].keys()
U_list.sort()
EmapU_list = []
for k in U_list:
if j == 'QE':
EmapU_value = ( SortElist[i][j][k] - SortElist[WRT_Morph][j][k] ) * convRytoeV
EmapU_list.append( EmapU_value )
else:
EmapU_value = SortElist[i][j][k] - SortElist[WRT_Morph][j][k]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j].append( EmapU_list )
EBsiteMatch[i][j].append( U_list )
ax = plt.gca()
print "PBE Functional, VASP Calculator:"
print "Rutile: " + str(EBsiteMatch['Rutile']['VASP'][0])
print "Columbite: " + str(EBsiteMatch['Columbite']['VASP'][0][0:2])
print str(EBsiteMatch['Columbite']['VASP'][0][2:4])
print str(EBsiteMatch['Columbite']['VASP'][0][4:]) + "\n"
print "PBE Functional, QE Calculator:"
print "Rutile: " + str(EBsiteMatch['Rutile']['QE'][0])
print "Columbite: " + str(EBsiteMatch['Columbite']['QE'][0][0:2])
print str(EBsiteMatch['Columbite']['QE'][0][2:4])
print str(EBsiteMatch['Columbite']['QE'][0][4:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Rutile']['VASP'][1], EBsiteMatch['Rutile']['VASP'][0], 'k-')
plt.plot(EBsiteMatch['Columbite']['VASP'][1], EBsiteMatch['Columbite']['VASP'][0],
'bo-', label = r'$\Delta$E$_{R-C,VASP}$')
plt.plot(EBsiteMatch['Columbite']['QE'][1], EBsiteMatch['Columbite']['QE'][0],
'ro-', label = r'$\Delta$E$_{R-C,QE}$')
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.02, 0.08) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper left', prop={'size':10})
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-VASPQE-PBE' + '.' + ext, dpi=300)
plt.clf()
Given that the Rutile-Columbite formation energetics evaluated in VASP and QE using the PBE functional largely overlap, the first-principles Hubbard
The first-principles method of linear response theory was applied to several TiO2 polymorphic systems, in order to assess the extent to which formation energy ordering relating to these polymorphs could be predicted, given known experimental outcomes (e.g.: Rutile is the most stable bulk TiO2 polymorph, Columbite is the least stable of the four studied, etc.). The types of systems evaluated in this study reflect the results of Hubbard
Regardless of the magnetic state observed or induced in a system, linear response calculations performed on a particular periodic system are accomplished by perturbing the orbital occupations of
\begin{equation} Uout = χ0-1 - χ-1 \label{eq-DFTULinResp} \end{equation}
Using VASP, this calculation can be completed via a four step procedure that employs the “LDAUTYPE = 3” setting revealed on the official VASP forum cite:VASPforum,Zhou2004. Firstly, take the relevant TiO2 structures achieved after variable cell volume relaxation using the “ISIF = 3” setting in Subsection Sample Input Files: Hubbard U, produce an appropriate supercell representation of that structure, and then complete an equivalent relaxation on that supercell representation. An example of this calculation step using the 2
ISTART = 0 ; ICHARG = 2
ENCUT = 600
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 1
IBRION = 1
EDIFF = 5E-05; EDIFFG = -0.02
MAXMIX = -100
NELMIN = 5
NELM = 200
NSW = 100
ISPIN = 2
ISIF = 3
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1
LDAUU = 0.00 0.00
LDAUJ = 0.00 0.00
LDAUPRINT= 2
LASPH = .TRUE.
LMAXMIX = 4
LORBIT = 11
4x4x4
0
Monkhorst-Pack
4 4 4
0 0 0
Ti O
4.66
1 0 0
0 1 0
0 0 0.643993896
2 4
direct
0.000000000 0.000000000 0.000000000
0.500000000 0.500000000 0.500000000
0.304880731 0.304880731 0.000000000
-0.304880731 -0.304880731 0.000000000
0.804875587 0.195124413 0.500000000
-0.804875587 -0.195124413 0.500000000
The calculation above was completed for the PM case, as is indicated by the presence of “ISPIN = 2” in the INCAR file. Completion of the NM case (with “ISPIN” = 1 in the INCAR file) reveals equivalent results, as will be shown later for both Rutile and Anatase polymorphs. The results achieved using the input files shown above are not directly related to any final calculated values in the paper, thus the input and output files associated with the step above are not present in the database. Creation of the supercells needed to complete this first step can be accomplished via the ase.io.vasp module and associated supercell command found in the Atomic Simulation Environment (ASE) software package. The command, which is reproduced in the code below, is performed on a POSCAR file in the directory containing the code itself: cite:ASE
import ase.io.vasp
cell = ase.io.vasp.read_vasp("POSCAR")
ase.io.vasp.write_vasp("POSCAR_2",cell*(2,2,2), label='Ti O', direct=True, sort=True)In the second step of this four step procedure, the CONTCAR of the structurally relaxed system from the first step is imported into a second calculation, namely one that performs solely electronic relaxation at a higher precision specified by the “EDIFF” tag. This step in the procedure is used to converge a CHGCAR file used to store orbital occupation information, which is subsequently imported into the third step of this calculation procedure: cite:VASPforum,Zhou2004
ENCUT = 600
ISMEAR = 0
SIGMA = 0.05
ISYM = 1
EDIFF = 1E-06
ISPIN = 2
LDAU = .TRUE.
LDAUTYPE = 2
LDAUL = 2 -1 -1
LDAUU = 0.00 0.00 0.00
LDAUJ = 0.00 0.00 0.00
LDAUPRINT= 2
LASPH = .TRUE.
LMAXMIX = 4
LORBIT = 11
4x4x4
0
Monkhorst-Pack
4 4 4
0 0 0
Ti Ti O
1.00000000000000
9.3217564470847751 0.0000000000000000 0.0000000000000000
0.0000000000000000 9.3217564470847751 0.0000000000000000
0.0000000000000000 0.0000000000000000 5.9372793380803959
1 15 32
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.7500000000000000 0.7500000000000000 0.7500000000000000
0.2500000000000000 0.2500000000000000 0.2500000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
0.7500000000000000 0.7500000000000000 0.2500000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.2500000000000000 0.2500000000000000 0.7500000000000000
0.2500000000000000 0.7500000000000000 0.7500000000000000
0.7500000000000000 0.2500000000000000 0.7500000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.2500000000000000 0.7500000000000000 0.2500000000000000
0.7500000000000000 0.2500000000000000 0.2500000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.0976626291255371 0.9023373708744629 0.7500000000000000
0.4023373708744629 0.5976626291255371 0.7500000000000000
0.9023373708744629 0.5976626291255371 0.7500000000000000
0.6523349025344345 0.1523349025344345 0.0000000000000000
0.8476650974655655 0.3476650974655655 0.0000000000000000
0.9023373708744629 0.0976626291255371 0.2500000000000000
0.5976626291255371 0.4023373708744629 0.2500000000000000
0.6523349025344345 0.1523349025344345 0.5000000000000000
0.3476650974655655 0.8476650974655655 0.5000000000000000
0.8476650974655655 0.3476650974655655 0.5000000000000000
0.5976626291255371 0.4023373708744629 0.7500000000000000
0.6523349025344345 0.6523349025344345 0.0000000000000000
0.8476650974655655 0.8476650974655655 0.0000000000000000
0.9023373708744629 0.5976626291255371 0.2500000000000000
0.5976626291255371 0.9023373708744629 0.2500000000000000
0.6523349025344345 0.6523349025344345 0.5000000000000000
0.8476650974655655 0.8476650974655655 0.5000000000000000
0.9023373708744629 0.0976626291255371 0.7500000000000000
0.1523349025344345 0.6523349025344345 0.5000000000000000
0.5976626291255371 0.9023373708744629 0.7500000000000000
0.3476650974655655 0.3476650974655655 0.5000000000000000
0.0976626291255371 0.9023373708744629 0.2500000000000000
0.4023373708744629 0.5976626291255371 0.2500000000000000
0.3476650974655655 0.8476650974655655 0.0000000000000000
0.1523349025344345 0.1523349025344345 0.5000000000000000
0.1523349025344345 0.1523349025344345 0.0000000000000000
0.1523349025344345 0.6523349025344345 0.0000000000000000
0.0976626291255371 0.4023373708744629 0.2500000000000000
0.4023373708744629 0.0976626291255371 0.7500000000000000
0.0976626291255371 0.4023373708744629 0.7500000000000000
0.3476650974655655 0.3476650974655655 0.0000000000000000
0.4023373708744629 0.0976626291255371 0.2500000000000000
Note that, in all steps of this procedure directly relating to the calculation of perturbations, the perturbed cation must be distinguished as a separate atomic species of the same composition as matching non-perturbed species in both atomic type (atom indices of the POSCAR file) and POTCAR designations. Thus, for TiO2 species, atomic type based information must be represented through three indices (e.g.: Ti Ti O) as oppposed to two indices (e.g.: Ti O). In the third step of this procedure, CHGCAR information created in the second step is imported into a non-selfconsistent calculation completed using “ICHARG = 11”. The final orbital occupations of this non-selfconsistent calculation are then used to produce the initial response (chi0) needed to calculate a first-principles
ICHARG = 11
ENCUT = 600
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 0
ISPIN = 2
EDIFF = 1E-07
LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 -1 -1
LDAUU = 0.10 0.00 0.00
LDAUJ = 0.10 0.00 0.00
LDAUPRINT = 2
LASPH = .TRUE.
LMAXMIX = 4
LORBIT = 11
Lastly, a self-consistent calculation is completed using calculation parameters matching that of the third step, except without implementing the converged “CHGCAR” file from the second step. An example of this final response INCAR file (chi), completed at a perturbation potential (α) of 0.10 eV performed in the fourth and final step of this procedure, is transcribed below (the KPOINTS and POSCAR files match that shown in the “Convergence” step above). Note that the perturbation range and incrementation (α = -0.150 – 0.150 eV, 0.05 eV increments) seen in the third step is also employed in this step: cite:VASPforum,Zhou2004
ENCUT = 600
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 0
ISPIN = 2
EDIFF = 1E-07
LDAU = .TRUE.
LDAUTYPE = 3
LDAUL = 2 -1 -1
LDAUU = 0.10 0.00 0.00
LDAUJ = 0.10 0.00 0.00
LDAUPRINT = 2
LASPH = .TRUE.
LMAXMIX = 4
LORBIT = 11
The input files above, which correspond to the case of the PM TiO2 Rutile polymorph, yielded results that can be compared to their NM analogues using a PBE functional and standard pseudopotentials. The effects of incorporating or neglecting spin polarization and magnetism during different steps of the procedure above were evaluated to assess the impact of considering PM and NM states in TiO2 systems. In the case of Rutile, calculations were completed while neglecting spin polarization and magnetism (NM case) in both the second (“convergence”) and third/fourth steps (“perturbation”) above (“NS-NS”), incorporating spin polarization and magnetism (PM case) in the “convergence” step but ignoring it in “perturbation” steps (“S-NS”), performing a NM “convergence” step and PM “perturbation” step (“NS-S”), and completing PM calculations across both steps (“S-S”). The results of these calculations, which are queried, displayed, and effectively compared below, reveal that spin polarization cannot be ignored in a convergence calculation and then incorporated into perturbation calculations due to the unphysical magnetic moments achieved when performing this combination (“NS-S”). However, in the case of TiO2 polymorphs, these combinations of calculations achieved otherwise equivalent perturbed atom occupation and magnetism results, thus NM and PM modeled TiO2 systems would possess the same value of
\begin{flalign*}
&ΠCalcType, AtomIndex, DChargeRHigh, DMagP
&σ [(Morph=’Rutile’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(CalcType=’Chi0NS-NS’) ∨ (CalcType=’Chi0NS-S’) \
&σ ∨ (CalcType=’Chi0S-NS’) ∨ (CalcType=’Chi0S-S’) \
&σ ∨ (CalcType=’ChiNS-NS’) ∨ (CalcType=’ChiNS-S’) \
&σ ∨ (CalcType=’ChiS-NS’) ∨ (CalcType=’ChiS-S’)] \
&σ ∧ [(AtomIndex=’1’)] \
&[Structure \Join Composition1 \Join Metadata \Join CalcResultLRCalc \
& \Join ParametersResponseVASP]
\end{flalign*}
#+caption: MySQL query Linear Response NS-NS NS-S S-NS S-S
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Calc_type, Atom_index, d_charge, d_mag = [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select prv.CalcType, prv.AtomIndex, prv.DChargeRHigh, prv.DMagP from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join CalcResultLRCalc as crlrc on crlrc.CID=mt.CID
inner join ParametersResponseVASP as prv on crlrc.Energy = prv.Energy
where s.Morph='Rutile'
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Functional='PBE'
and mt.Software='VASP'
and crlrc.U3dIn='0'
and (prv.CalcType='Chi0NS-NS' or prv.CalcType='Chi0NS-S' or prv.CalcType='Chi0S-NS'
or prv.CalcType='Chi0S-S' or prv.CalcType='ChiNS-NS'
or prv.CalcType='ChiNS-S' or prv.CalcType='ChiS-NS' or prv.CalcType='ChiS-S')
and prv.AtomIndex='1'
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
Calc_type += [a]
Atom_index += [b]
datapts_dict[row] = datapts_list
CalcType_list = list( set( Calc_type ) )
AtomIndex_list = list( set( Atom_index ) )
SortElist = {}
EBsiteMatch = {}
for i in CalcType_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in AtomIndex_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ] = [ float( k[2] ), float( k[3] ) ]
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
print "Linear response orbital occupancy perturbations:"
print "No spin (convergence step), no spin (perturbation step) [NS-NS]: " + str(SortElist['Chi0NS-NS'][1][0])
print "No spin (convergence step), spin (perturbation step) [NS-S]: " + str(SortElist['Chi0NS-S'][1][0])
print "Spin (convergence step), no spin (perturbation step) [S-NS]: " + str(SortElist['Chi0S-NS'][1][0])
print "Spin (convergence step), spin (perturbation step) [S-S]: " + str(SortElist['Chi0S-S'][1][0]) + "\n"
print "Linear response magnetic moments (perturbed atom, index 1):"
print "No spin (convergence step), no spin (perturbation step) [NS-NS]: " + str(SortElist['Chi0NS-NS'][1][1])
print "No spin (convergence step), spin (perturbation step) [NS-S]: " + str(SortElist['Chi0NS-S'][1][1])
print "Spin (convergence step), no spin (perturbation step) [S-NS]: " + str(SortElist['Chi0S-NS'][1][1])
print "Spin (convergence step), spin (perturbation step) [S-S]: " + str(SortElist['Chi0S-S'][1][1]) + "\n"
Particular linear response calculations were completed using a previous response matrix calculation methodology adapted into a Python code shown below. cite:Cococcioni2005 Each matrix is an
\begin{equation} Uout = (χ0-1)i,i - (χ-1)i,i ← \frac{n0(α = 0.15)i,i - n0(α = -0.15)i,i}{0.15 - (-0.15)} - \frac{n(α = 0.15)i,i - n(α = -0.15)i,i}{0.15 - (-0.15)} \label{eq-ActualLinResp} \end{equation}
Upon populating the first column with terms of the above form, equivalencies between any two pairs of atomic sites are assessed in terms of several measures. For periodic crystal structures, atomic site pairs are equated via interatomic distance, spin state (magnetism), atomic composition, and Wyckoff position cite:Cococcioni2005,Aroyo2006 Given that two atomic pairs share equivalent information relative to all of these measures, the first-order difference created in one equated response matrix entry (i.e.: interaction) is repeated in the other equated entry. Inversion of the matrices corresponding to the initial (chi0) and final (chi) responses, upon population of each entire matrix in accordance with the aforementioned atomic pair equivalence criteria mentioned above, reveals the quantities that are to be subtracted to calculate a first-principles
In order to calculate the uncertainty imposed on a single matrix entry by inversion of an entire matrix exactly, the implementation of Cayley-Hamilton theorem or a similar principle would be required to allow each matrix term contributing to uncertainty to be subjectable to common error propagation techniques. cite:Lefebvre2000,Ghaoui2002 In the case of matrices wherein single diagonal elements contribute predominately to matrix inversion operations, the uncertainty for those operations can be approximated as the uncertainty of inverting that single entry, namely by applying multiplication (division or inversion) in quadrature. In the linear response calculations performed in this paper, the condition of predominately large single diagonal elements can always be considered satisfied, given that the diagonal matrix terms associated with Ti cation perturbation are always at least one order of magnitude larger than implemented and observed off-diagonal terms. Thus, the uncertainty associated with matrix inversion is approximated as that of the single perturbed cation matrix entry of interest in each case. After matrix inversion, the subtraction of the perturbed cation inverted matrix entries renders a final contribution to uncertainty, namely via application of addition (subtraction) in quadrature. cite:Lefebvre2000,Ghaoui2002 Thus, the overall expression for the uncertainty of a
\begin{equation} σprec = \frac{\sqrt{ (ε)2 + (ε)2 }}{Δ(α)} \end{equation}
\begin{equation} fdiag,chi_{0-1/chi-1} = \frac{n(α = 0.15) - n(α = -0.15)}{Δ(α)} \end{equation}
\begin{equation} σchi_{0-1/chi-1} = \sqrt{ (fdiag,chi_{0-1/chi-1})-2 × \frac{(σprec)2}{(fdiag,chi_{0-1/chi-1})2} } = \frac{σprec}{(fdiag,chi_{0-1/chi-1})2} \end{equation}
\begin{equation} σ = \sqrt{ (σchi_{0-1})2 + (σchi^{-1})2 } \end{equation}
In the first expression above, the uncertainty associated with the measurement precision of any response matrix entry (σprec) is calculated. In the second expression above, the denominators of the relative errors needed for the error propagation calculations performed in the third expression is completed. These denominator values, fdiag,chi_{0-1/chi-1}, correspond to the initial (chi0-1) and final (chi-1) response contributions to the linear response calculated value of
Calculations of all response matrices featuring these uncertainty propagation techniques, the calculation of
\begin{flalign*}
&ΠPseudoB, Functional, Morph, AtomIndex, U3dOut, CalcType, Chi0[N150-P150], Chi[N150-P150]
&ΠDChargeRHigh, DChargeRLow, DChargeRCenter, PChargeRHigh, PChargeRLow, PChargeRCenter \
&ΠCoord[1-96], PturbMax, PturbMin \
&σ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) \
&σ ∨ (Morph=’Columbite’) ∨ (Morph=’Brookite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’) ∨ (PseudoB=’PAW-Ti-pv’) ∨ (PseudoB=’PAW-Ti-sv’)] \
&σ ∧ [(Functional=’PBE’) ∨ (Functional=’PS’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’LR’)] \
&σ ∧ [(U3dIn=’0’)] \
&[Structure \Join Composition1 \Join Metadata \Join CalcResultLRCalc \
& \Join ParametersResponseVASP \Join ParametersPosFinalVASP]
\end{flalign*}
#+caption: MySQL query Linear Response Calculation
import matplotlib.pyplot as plt
import sqlite3
import numpy as np
from scipy import stats
NULL_CHAR = '-' # NULL value in database
SColon_CHAR = ';'
comp_back = -1E-3 # computed background placed on on-diagonal unperturbed
# atomic indices distinguished from the perturbed atom
pturb_index = 0 # index of perturbed atom
epsilon_error = 0.001 # measurement error
WyckoffSym = '1' # to be improved in future work
MagVal = '0' # to be improved in future work
db = sqlite3.connect('ITEOO_data.sqlite')
Pseudo_B, Functional, Polymorph, CalcType, AtomIndices = [], [], [], [], []
coord_dict = {}
for row in db.execute('''
select distinct ppfv.Energy,
ppfv.Coord1, ppfv.Coord2, ppfv.Coord3, ppfv.Coord4, ppfv.Coord5,
ppfv.Coord6, ppfv.Coord7, ppfv.Coord8, ppfv.Coord9, ppfv.Coord10,
ppfv.Coord11, ppfv.Coord12, ppfv.Coord13, ppfv.Coord14, ppfv.Coord15,
ppfv.Coord16, ppfv.Coord17, ppfv.Coord18, ppfv.Coord19, ppfv.Coord20,
ppfv.Coord21, ppfv.Coord22, ppfv.Coord23, ppfv.Coord24, ppfv.Coord25,
ppfv.Coord26, ppfv.Coord27, ppfv.Coord28, ppfv.Coord29, ppfv.Coord30,
ppfv.Coord31, ppfv.Coord32, ppfv.Coord33, ppfv.Coord34, ppfv.Coord35,
ppfv.Coord36, ppfv.Coord37, ppfv.Coord38, ppfv.Coord39, ppfv.Coord40,
ppfv.Coord41, ppfv.Coord42, ppfv.Coord43, ppfv.Coord44, ppfv.Coord45,
ppfv.Coord46, ppfv.Coord47, ppfv.Coord48, ppfv.Coord49, ppfv.Coord50,
ppfv.Coord51, ppfv.Coord52, ppfv.Coord53, ppfv.Coord54, ppfv.Coord55,
ppfv.Coord56, ppfv.Coord57, ppfv.Coord58, ppfv.Coord59, ppfv.Coord60,
ppfv.Coord61, ppfv.Coord62, ppfv.Coord63, ppfv.Coord64, ppfv.Coord65,
ppfv.Coord66, ppfv.Coord67, ppfv.Coord68, ppfv.Coord69, ppfv.Coord70,
ppfv.Coord71, ppfv.Coord72, ppfv.Coord73, ppfv.Coord74, ppfv.Coord75,
ppfv.Coord76, ppfv.Coord77, ppfv.Coord78, ppfv.Coord79, ppfv.Coord80,
ppfv.Coord81, ppfv.Coord82, ppfv.Coord83, ppfv.Coord84, ppfv.Coord85,
ppfv.Coord86, ppfv.Coord87, ppfv.Coord88, ppfv.Coord89, ppfv.Coord90,
ppfv.Coord91, ppfv.Coord92, ppfv.Coord93, ppfv.Coord94, ppfv.Coord95,
ppfv.Coord96 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join CalcResultLRCalc as crlrc on crlrc.CID=mt.CID
inner join ParametersResponseVASP as prv on prv.Energy=crlrc.Energy
inner join ParametersPosFinalVASP as ppfv on ppfv.Energy=prv.Energy
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or s.Morph='Brookite')
and c1.At1='Ti'
and (mt.PseudoB='PAW-Ti' or mt.PseudoB='PAW-Ti-pv' or mt.PseudoB='PAW-Ti-sv')
and (mt.Functional='PBE' or mt.Functional='PS')
and mt.Software='VASP'
and mt.Method='LR'
and crlrc.U3dIn='0'
and (prv.CalcType='Chi' or prv.CalcType='Chi0' or prv.CalcType='ChiS-S'
or prv.CalcType='Chi0S-S')
;'''):
coord_list = []
for i in range(1, len(row), 1):
if row[i] != NULL_CHAR:
coord_list.append( row[i] )
coord_dict[ row[0] ] = coord_list
charge_dict = {}
CalcType, AtomIndices = [], []
for row in db.execute('''
select prv.Energy, prv.CalcType, prv.AtomIndex, prv.DChargeRHigh, prv.DChargeRLow, prv.DChargeRCenter,
prv.PChargeRHigh, prv.PChargeRLow, prv.PChargeRCenter from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join CalcResultLRCalc as crlrc on crlrc.CID=mt.CID
inner join ParametersResponseVASP as prv on prv.Energy=crlrc.Energy
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or
s.Morph='Brookite')
and c1.At1='Ti'
and (mt.PseudoB='PAW-Ti' or mt.PseudoB='PAW-Ti-pv' or mt.PseudoB='PAW-Ti-sv')
and (mt.Functional='PBE' or mt.Functional='PS')
and mt.Software='VASP'
and mt.Method='LR'
and crlrc.U3dIn='0'
and (prv.CalcType='Chi' or prv.CalcType='Chi0' or prv.CalcType='ChiS-S'
or prv.CalcType='Chi0S-S')
;'''):
CalcType.append( row[1] )
AtomIndices.append( row[2] )
charge_list = []
for i in range(1, len(row), 1):
charge_list.append( row[i] )
if row[0] not in charge_dict:
charge_dict[ row[0] ] = {}
try:
charge_dict[ row[0] ][ row[1] ][ row[2] ] = charge_list
except KeyError:
charge_dict[ row[0] ][ row[1] ] = {}
charge_dict[ row[0] ][ row[1] ][ row[2] ] = charge_list
Uturb_dict = {}
Pseudo, Functional, Polymorph = [], [], []
for row in db.execute('''
select crlrc.Energy, mt.PseudoB, mt.Functional, s.Morph, crlrc.U3dOut,
crlrc.Chi0N150, crlrc.Chi0N100, crlrc.Chi0N050, crlrc.Chi0P000, crlrc.Chi0P050,
crlrc.Chi0P100, crlrc.Chi0P150, crlrc.ChiN150, crlrc.ChiN100, crlrc.ChiN050, crlrc.ChiP000, crlrc.ChiP050, crlrc.ChiP100,
crlrc.ChiP150, crlrc.PturbMax, crlrc.PturbMin from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join CalcResultLRCalc as crlrc on crlrc.CID=mt.CID
inner join ParametersResponseVASP as prv on prv.Energy=crlrc.Energy
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite' or s.Morph='Brookite')
and c1.At1='Ti'
and (mt.PseudoB='PAW-Ti' or mt.PseudoB='PAW-Ti-pv' or mt.PseudoB='PAW-Ti-sv')
and (mt.Functional='PBE' or mt.Functional='PS')
and mt.Software='VASP'
and mt.Method='LR'
and crlrc.U3dIn='0'
and (prv.CalcType='Chi' or prv.CalcType='Chi0' or prv.CalcType='ChiS-S'
or prv.CalcType='Chi0S-S')
;'''):
Pseudo.append( row[1] )
Functional.append( row[2] )
Polymorph.append( row[3] )
Uturb_list = []
for i in range(1, len(row), 1):
Uturb_list.append( row[i] )
Uturb_dict[ row[0] ] = Uturb_list
combine_dict = {}
for i in Uturb_dict.keys():
combine_dict[i] = [ coord_dict[i], charge_dict[i], Uturb_dict[i] ]
Pseudo_list = list( set (Pseudo) )
Functional_list = list( set( Functional ) )
Morph_list = list( set( Polymorph ) )
CalcType_list = list( set( CalcType ) )
AtomIndices_list = list( set( AtomIndices ) )
SortElist = {}
EBsiteMatch = {}
for i in Pseudo_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Functional_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = {}
for k in Morph_list:
EBsiteMatch[i][j][k] = {}
EBsiteMatch[i][j][k]['Chi'] = {}
EBsiteMatch[i][j][k]['Chi0'] = {}
del_list = []
for l in combine_dict.keys():
if combine_dict[l][2][0] == i and combine_dict[l][2][1] == j and combine_dict[l][2][2] == k:
SortElist[i][j][k] = combine_dict[l]
del_list.append( combine_dict[l] )
for k in del_list:
del k
for i in Pseudo_list:
for j in Functional_list:
for k in Morph_list:
try:
coord_list = []
for l,m in enumerate( SortElist[i][j][k][0] ):
if SortElist[i][j][k][0][l] != NULL_CHAR:
coord_value = str(m).split(SColon_CHAR)
coord_list.append( [ float( coord_value[0] ), float(
coord_value[1] ), float( coord_value[2] ) ] )
EBsiteMatch[i][j][k]['coord'] = coord_list
EBsiteMatch[i][j][k]['U_diag'] = SortElist[i][j][k][2][3]
EBsiteMatch[i][j][k]['Chi0_Pturb'] = [ SortElist[i][j][k][2][4],
SortElist[i][j][k][2][5],
SortElist[i][j][k][2][6],
SortElist[i][j][k][2][7],
SortElist[i][j][k][2][8],
SortElist[i][j][k][2][9],
SortElist[i][j][k][2][10]]
EBsiteMatch[i][j][k]['Chi_Pturb'] = [ SortElist[i][j][k][2][11],
SortElist[i][j][k][2][12],
SortElist[i][j][k][2][13],
SortElist[i][j][k][2][14],
SortElist[i][j][k][2][15],
SortElist[i][j][k][2][16],
SortElist[i][j][k][2][17] ]
EBsiteMatch[i][j][k]['PturbMax'] = SortElist[i][j][k][2][18]
EBsiteMatch[i][j][k]['PturbMin'] = SortElist[i][j][k][2][19]
Pturb_divisor = float( len( EBsiteMatch[i][j][k]['Chi0_Pturb'] )
- 1 )
Pturbcept_index = int( Pturb_divisor / 2. )
EBsiteMatch[i][j][k]['List_Pturb'] = np.arange(
EBsiteMatch[i][j][k]['PturbMin'],
EBsiteMatch[i][j][k]['PturbMax'] + 0.0001,
( EBsiteMatch[i][j][k]['PturbMax'] -
EBsiteMatch[i][j][k]['PturbMin'] ) / Pturb_divisor )
EBsiteMatch[i][j][k]['Intercept'] = abs(
EBsiteMatch[i][j][k]['Chi0_Pturb'][Pturbcept_index] -
EBsiteMatch[i][j][k]['Chi_Pturb'][Pturbcept_index] )
slope, intercept, r_value, p_value, std_err = stats.linregress(
EBsiteMatch[i][j][k]['List_Pturb'],
EBsiteMatch[i][j][k]['Chi0_Pturb'] )
EBsiteMatch[i][j][k]['Chi0_R2'] = r_value**2
slope, intercept, r_value, p_value, std_err = stats.linregress(
EBsiteMatch[i][j][k]['List_Pturb'],
EBsiteMatch[i][j][k]['Chi_Pturb'] )
EBsiteMatch[i][j][k]['Chi_R2'] = r_value**2
for l in range( len(EBsiteMatch[i][j][k]['coord']) ):
try:
if SortElist[i][j][k][1]['ChiS-S'][l+1][2] == 0.0:
EBsiteMatch[i][j][k]['Chi'][l+1] = (
float(SortElist[i][j][k][1]['ChiS-S'][l+1][5]) -
float(SortElist[i][j][k][1]['ChiS-S'][l+1][6]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
EBsiteMatch[i][j][k]['Chi0'][l+1] = (
float(SortElist[i][j][k][1]['Chi0S-S'][l+1][5]) -
float(SortElist[i][j][k][1]['Chi0S-S'][l+1][6]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
else:
EBsiteMatch[i][j][k]['Chi'][l+1] = (
float(SortElist[i][j][k][1]['ChiS-S'][l+1][2]) -
float(SortElist[i][j][k][1]['ChiS-S'][l+1][3]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
EBsiteMatch[i][j][k]['Chi0'][l+1] = (
float(SortElist[i][j][k][1]['Chi0S-S'][l+1][2]) -
float(SortElist[i][j][k][1]['Chi0S-S'][l+1][3]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
except KeyError:
if SortElist[i][j][k][1]['Chi'][l+1][2] == 0.0:
EBsiteMatch[i][j][k]['Chi'][l+1] = (
float(SortElist[i][j][k][1]['Chi'][l+1][5]) -
float(SortElist[i][j][k][1]['Chi'][l+1][6]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
EBsiteMatch[i][j][k]['Chi0'][l+1] = (
float(SortElist[i][j][k][1]['Chi0'][l+1][5]) -
float(SortElist[i][j][k][1]['Chi0'][l+1][6]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
else:
EBsiteMatch[i][j][k]['Chi'][l+1] = (
float(SortElist[i][j][k][1]['Chi'][l+1][2]) -
float(SortElist[i][j][k][1]['Chi'][l+1][3]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
EBsiteMatch[i][j][k]['Chi0'][l+1] = (
float(SortElist[i][j][k][1]['Chi0'][l+1][2]) -
float(SortElist[i][j][k][1]['Chi0'][l+1][3]) ) / (
float(EBsiteMatch[i][j][k]['PturbMax'])
- float(EBsiteMatch[i][j][k]['PturbMin']) )
except KeyError:
pass
for i in Pseudo_list:
for j in Functional_list:
for k in Morph_list:
try:
sigma_error_1 = ( ( epsilon_error**(2) +
epsilon_error**(2) )**(0.5)
/ ( EBsiteMatch[i][j][k]['PturbMax'] -
EBsiteMatch[i][j][k]['PturbMin'] ) )
sigma_error_2_chi = ( EBsiteMatch[i][j][k]['Chi'][1] )**(-2.0)
sigma_error_2_chi0 = ( EBsiteMatch[i][j][k]['Chi0'][1] )**(-2.0)
EBsiteMatch[i][j][k]['Sigma'] = (
(sigma_error_1 * sigma_error_2_chi)**(2.0) +
(sigma_error_1 * sigma_error_2_chi0)**(2.0) )**(0.5)
except KeyError:
pass
for i in Pseudo_list:
for j in Functional_list:
for k in Morph_list:
try:
respmat_dim = len( EBsiteMatch[i][j][k]['coord'] )
respmat_chi0 = np.zeros( (respmat_dim, respmat_dim) )
respmat_chi = np.zeros( (respmat_dim, respmat_dim) )
for l in range( respmat_dim ):
respmat_chi0[0][l] = EBsiteMatch[i][j][k]['Chi0'][l+1]
respmat_chi[0][l] = EBsiteMatch[i][j][k]['Chi'][l+1]
for l in range( respmat_dim ):
respmat_chi0[l][0] = respmat_chi0[0][l]
respmat_chi[l][0] = respmat_chi[0][l]
respmat_intxn = {}
for l in range( 0, respmat_dim, 1 ):
respmat_intxn[l] = {}
for l in range( 0, respmat_dim, 1 ):
respmat_intxn[l][l] = [ abs( np.subtract(
EBsiteMatch[i][j][k]['coord'][l],
EBsiteMatch[i][j][k]['coord'][l]) ),
[ WyckoffSym, WyckoffSym ], [ MagVal, MagVal ] ]
for m in range( 1, respmat_dim, 1 ):
respmat_intxn[l][m] = [ abs( np.subtract(
EBsiteMatch[i][j][k]['coord'][l],
EBsiteMatch[i][j][k]['coord'][m]) ),
[ WyckoffSym, WyckoffSym ], [ MagVal, MagVal ] ]
respmat_intxn[m][l] = respmat_intxn[l][m]
for l in range( 1, respmat_dim, 1 ):
for m in range( 2, respmat_dim, 1 ):
for n in range( 0, respmat_dim, 1 ):
counter_sym = 0
for o, p in zip( range( len(respmat_intxn[l][m][1])
), range( len(respmat_intxn[0][n][1]) ) ):
if ( respmat_intxn[l][m][1][o] !=
respmat_intxn[l][m][1][p] ):
counter_sym += 1
if ( respmat_intxn[l][m][2][o] !=
respmat_intxn[l][m][2][p] ):
counter_sym += 1
if counter_sym == 0:
respmat_chi0[l][m] = respmat_chi0[0][n]
respmat_chi[l][m] = respmat_chi[0][n]
for l in range( 0, respmat_dim, 1 ):
for m in range( 1, respmat_dim, 1 ):
respmat_chi0[m][l] = respmat_chi0[l][m]
respmat_chi[m][l] = respmat_chi[l][m]
for l in range(1, respmat_dim, 1 ):
if l < int( len( EBsiteMatch[i][j][k]['coord'] ) / 3. ):
respmat_chi0[l][l] = respmat_chi0[0][0]
respmat_chi[l][l] = respmat_chi[0][0]
else:
respmat_chi0[l][l] = comp_back
respmat_chi[l][l] = comp_back
np.savetxt('./figures/respmat_chi0_' + str(i) + '_' + str(j)
+ '_' + str(k) + '.txt', respmat_chi0)
np.savetxt('./figures/respmat_chi_' + str(i) + '_' + str(j)
+ '_' + str(k) + '.txt', respmat_chi)
try:
invmat_chi0 = np.linalg.inv(respmat_chi0)
invmat_chi = np.linalg.inv(respmat_chi)
except (np.linalg.linalg.LinAlgError):
print "Error for system: " + str(i) + "," + str(j) + "," + str(k)
print respmat_chi0
print respmat_chi
import sys; sys.exit()
EBsiteMatch[i][j][k]['U_val'] = (invmat_chi[ pturb_index + 1 ][ pturb_index + 1 ]
- invmat_chi0[ pturb_index + 1 ][ pturb_index + 1 ])
except KeyError:
pass
print "PBEsol Functional, Ti pseudopotential, Rutile: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti']['PS']['Rutile']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti']['PS']['Rutile']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti']['PS']['Rutile']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti']['PS']['Rutile']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti']['PS']['Rutile']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti']['PS']['Rutile']['Sigma']) + "\n"
print "PBEsol Functional, Ti pseudopotential, Anatase: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti']['PS']['Anatase']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti']['PS']['Anatase']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti']['PS']['Anatase']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti']['PS']['Anatase']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti']['PS']['Anatase']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti']['PS']['Anatase']['Sigma']) + "\n"
print "PBE Functional, Ti-pv pseudopotential, Rutile: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Rutile']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Rutile']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti-pv']['PBE']['Rutile']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Rutile']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Rutile']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Rutile']['Sigma']) + "\n"
print "PBE Functional, Ti-pv pseudopotential, Anatase: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Anatase']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Anatase']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti-pv']['PBE']['Anatase']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Anatase']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Anatase']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti-pv']['PBE']['Anatase']['Sigma']) + "\n"
print "PBE Functional, Ti-sv pseudopotential, Rutile: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Rutile']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Rutile']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti-sv']['PBE']['Rutile']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Rutile']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Rutile']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Rutile']['Sigma']) + "\n"
print "PBE Functional, Ti-sv pseudopotential, Anatase: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Anatase']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Anatase']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti-sv']['PBE']['Anatase']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Anatase']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Anatase']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti-sv']['PBE']['Anatase']['Sigma']) + "\n"
print "PBE Functional, Ti pseudopotential, Rutile: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['Sigma']) + "\n"
print "PBE Functional, Ti pseudopotential, Anatase: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Anatase']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Anatase']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti']['PBE']['Anatase']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Anatase']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Anatase']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti']['PBE']['Anatase']['Sigma']) + "\n"
print "PBE Functional, Ti pseudopotential, Columbite: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Columbite']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Columbite']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti']['PBE']['Columbite']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Columbite']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Columbite']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti']['PBE']['Rutile']['Sigma']) + "\n"
print "PBE Functional, Ti pseudopotential, Brookite: "
print r'Coefficient of correlation (R^2) of Chi = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Brookite']['Chi_R2'])
print r'Coefficient of correlation (R^2) of Chi0 = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Brookite']['Chi0_R2'])
print r'Intercept (|Chi0(alpha = 0 eV) - Chi(alpha = 0 eV)|) = '
print str(EBsiteMatch['PAW-Ti']['PBE']['Brookite']['Intercept'])
print r'U(diag) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Brookite']['U_diag'])
print r'U(3d) = ' + str(EBsiteMatch['PAW-Ti']['PBE']['Brookite']['U_val'])
print ' +/- ' + str(EBsiteMatch['PAW-Ti']['PBE']['Brookite']['Sigma']) + "\n"
As shown in the output above and mentioned in the article, the linear response resolved values of
Considering that standard linear response calculations are initialized using electronic structure information provided from GGA rather than GGA+\(U\) ground states, changes in the electronic and crystal structure of the material induced by addition of the \(U\) parameter are handled approximately. In general, GGA-based functionals incorrectly represent energetic ground states of strongly correlated systems by significantly overestimating electron-electron interaction in pertinent orbitals (frequently
For a set of inputted values of
\begin{flalign*}
&ΠU3dIn, U3dOut
&σ [(Morph=’Rutile’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’LR’)] \
&σ ∧ [(U3dIn=’0’) ∨ (U3dIn=’1’) ∨ (U3dIn=’2’) ∨ (U3dIn=’2.5’) \
&σ ∨ (U3dIn=’3’) ∨ (U3dIn=’3.5’) ∨ (U3dIn=’4’)] \
&[Structure \Join Composition1 \Join Metadata \Join CalcResultLRCalc]
\end{flalign*}
#+caption: MySQL query Self-Consistent Linear Response Theory
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
U3dIn, U3dOut = [], []
NULL_CHAR = '-'
for row in db.execute('''
select distinct crlrc.U3dIn, crlrc.U3dOut from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join CalcResultLRCalc as crlrc on crlrc.CID=mt.CID
where s.Morph='Rutile'
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.Functional='PBE'
and mt.Software='VASP'
and mt.Method='LR'
and (crlrc.U3dIn='0' or crlrc.U3dIn='1' or crlrc.U3dIn='2' or crlrc.U3dIn='2.5'
or crlrc.U3dIn='3' or crlrc.U3dIn='3.5' or crlrc.U3dIn='4')
;'''):
a,b = row
if b != NULL_CHAR:
U3dIn.append( a )
U3dOut.append( b )
ax = plt.gca()
print "PBE Functional, Ti O pseudopotentials, Rutile: \n"
print "U (input): " + str(U3dIn) + "\n"
print "U (output): " + str(U3dOut[0:3])
print str(U3dOut[3:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(U3dIn, U3dOut, 'ro-', label = r'U$_{out}$, Rutile')
plt.xticks( [0, 1, 2, 3, 4] )
plt.xlabel( 'Input U (eV)' )
plt.ylim( (3.0, 4.0) )
plt.ylabel( 'Output U (eV)' )
plt.legend(loc = 'upper left', prop={'size':10})
plt.gcf().subplots_adjust(left=0.20)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-LRSC-R' + '.' + ext, dpi=300)
plt.clf()
The query and plot above illustrates an attempt to achieve a self-consistent value of
In this study, hybrid functional calculations are completed primarily for the Rutile and Anatase TiO2 polymorphs, applying the PBE functional for non-exact exchange, the PBE0 cite:Perdew1996_PRL and HSE06 cite:Heyd2003,Heyd2006 methodologies for mixing PBE and exact HF exchange, and the standard Ti and O pseudopotentials to complete calculations. In a multi-step procedure similar to that completed in Hubbard
A single example set of input files characterizing the predominate implementation of the procedure detailed above is described as follows. The first step of the procedure developed in Section Sample Input Files: Hubbard U, namely that involving the iterative use of fixed cell shape, variable cell volume and atomic coordinate (ISIF = 4 in VASP) calculations, is performed for all PBE calculations of interest in hybrid functional calculations. Subsequently, the CONTCAR and WAVECAR produced from that calculation is applied to a hybrid functional structural relaxation calculation (also ISIF = 4) with a particular value of ENCUT and NKRED. In the event that electronic convergence was difficult to reach for a particular system, hybrid functional calculations employing lower values of ENCUT (e.g.: 400 eV) and/or higher values of NKRED (e.g.: 3, rather than 2) initially received the WAVECAR and CONTCAR inputs from PBE calculations. Subsequently, the output of those low ENCUT and/or high NKRED calculations served as inputs for calculations used within this article, the input files of which are reproduced below. Note that, in cases where initial calculations with lower ENCUT or higher NKRED values were used as intermediates for more expensive calculations, system volume was always conserved over PBE, initial hybrid, and implemented hybrid calculations. Also, note that the INCAR and KPOINTS input file information used in these intermediate steps is, excluding differences in ENCUT and NKRED information, equivalent to that of analogous implemented hybrid functional calculations. The example below corresponds to an HSE06 calculation involving TiO2 Rutile with 25% exact HF exchange, 75% PBE exchange at a fixed cell volume of 62.35 A3.
ISTART = 1 ; ICHARG = 0
ENCUT = 550
ISMEAR = 0 ; SIGMA = 0.05
ISYM = 1; SYMPREC = 1E-06
IBRION = 1
EDIFF = 5E-05; EDIFFG = -0.02
MAXMIX = -75
NELMIN = 5
NELM = 250
NSW = 100
ISPIN = 2
ISIF = 4
LHFCALC = .TRUE.
HFSCREEN = 0.2
ALGO = Damped
TIME = 0.4
PRECFOCK = Normal
LMAXFOCK = 4
LASPH = .TRUE.
AEXX = 0.25
NKRED = 2
NBANDS = 25
6x6x6
0
Gamma
6 6 6
0 0 0
Ti O
4.61000000000000
0.9980395576792923 0.0000037613141482 0.0000000000000000
0.0000037613141482 0.9980395576792923 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.6388758158322786
Ti O
2 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.3047876991736018 0.3047876991736018 0.0000000000000000
0.6952123008263982 0.6952123008263982 0.0000000000000000
0.8047916748910851 0.1952083251089149 0.5000000000000000
0.1952083251089149 0.8047916748910851 0.5000000000000000
Results from the calculations of the following form presented above, namely WAVECAR and CONTCAR files, can be integrated into calculations employing different fractions of exact exchange (AEXX). This can be performed to complete calculations at a higher fraction of exact exchange with reduced computational expense. All calculations completed using this technique used input parameters identical to those presented above, except for appropriately substituted values of AEXX.
For all fractions of exact exchange evaluated for Rutile-Anatase hybrid functional formation energetics (
\begin{flalign*}
&ΠMorph, Functional, FractionHF, EOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) ∨ (Morph=’Columbite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(NKRED=’2’)] \
&σ ∧ [(Functional=’HSE06’) ∨ (Functional=’PBE0’)] \
&σ ∧ [(FractionHF=’0.25’) ∨ (FractionHF=’0.5’) ∨ (FractionHF=’0.75’) \
&σ ∨ (FractionHF=’0.825’) ∨ (FractionHF=’0.875’) ∨ (FractionHF=’0.95’) \
&σ ∨ (FractionHF=’1’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEHFRange \
& \Join CalcResultEnergetics \Join ParametersInputVASP]
\end{flalign*}
\begin{flalign*}
&ΠMorph, UValue, EOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(UValue=’0’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query PBE0 vs. HSE06 Rutile Anatase
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, Functional, FractionHF, Eopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select distinct s.Morph, mt.Functional, feh.FractionHF, feh.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEHFRange as feh on feh.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feh.EOpt
inner join ParametersInputVASP as input on input.Energy=cre.Energy
where (s.Morph='Rutile' or s.Morph='Anatase' or s.Morph='Columbite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and input.NKRED='2'
and (mt.Functional='HSE06' or mt.Functional='PBE0')
and (feh.FractionHF='0.25' or feh.FractionHF='0.5' or feh.FractionHF='0.75'
or feh.FractionHF='0.825' or feh.FractionHF='0.875'
or feh.FractionHF='0.95' or feh.FractionHF='1')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
Polymorph += [a]
Functional += [b]
datapts_dict[row] = datapts_list
Morph_list = list( set( Polymorph ) )
Functional_list = list( set( Functional ) )
SortElist = {}
EBsiteMatch = {}
for i in Morph_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Functional_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
if not SortElist[i][j]:
del SortElist[i][j]
Polymorph, Uvalue, Eopt, atomnum = [], [], [], []
datapts_dict = {}
for row in db.execute('''
select distinct s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on mt.CID=feu.CID
where (s.Morph='Rutile' or s.Morph='Anatase')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and mt.Functional='PBE'
and feu.UValue='0'
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
Polymorph += [a]
datapts_dict[row] = datapts_list
Morph_list_0 = list( set( Polymorph ) )
for i in Morph_list_0:
SortElist[i]['PBE0'][float(0.)] = []
SortElist[i]['HSE06'][float(0.)] = []
del_list = []
for j in datapts_dict:
if j[0] == i:
SortElist[ j[0] ]['PBE0'][float(0.)] = float( j[2] ) / float( j[3] )
SortElist[ j[0] ]['HSE06'][float(0.)] = float( j[2] ) / float( j[3] )
del_list.append( datapts_dict[j] )
else:
pass
for i in Morph_list:
for j in Functional_list:
try:
HF_list = SortElist[i][j].keys()
HF_list.sort()
EmapHF_list = []
for k in HF_list:
EmapHF_value = SortElist[i][j][k] - SortElist[WRT_Morph][j][k]
EmapHF_list.append( EmapHF_value )
EBsiteMatch[i][j].append( EmapHF_list )
EBsiteMatch[i][j].append( HF_list )
except KeyError:
pass
ax = plt.gca()
print "PBE0 Functional:"
print "Rutile: " + str(EBsiteMatch['Rutile']['PBE0'][0][0:2])
print str(EBsiteMatch['Rutile']['PBE0'][0][2:5])
print str(EBsiteMatch['Rutile']['PBE0'][0][5:])
print "Anatase: " + str(EBsiteMatch['Anatase']['PBE0'][0][0:2])
print str(EBsiteMatch['Anatase']['PBE0'][0][2:5])
print str(EBsiteMatch['Anatase']['PBE0'][0][5:]) + "\n"
print "HSE06 Functional:"
print "Rutile: " + str(EBsiteMatch['Rutile']['HSE06'][0][0:2])
print str(EBsiteMatch['Rutile']['HSE06'][0][2:5])
print str(EBsiteMatch['Rutile']['HSE06'][0][5:])
print "Anatase: " + str(EBsiteMatch['Anatase']['HSE06'][0][0:2])
print str(EBsiteMatch['Anatase']['HSE06'][0][2:5])
print str(EBsiteMatch['Anatase']['HSE06'][0][5:])
print "Columbite (a = " + str(EBsiteMatch['Columbite']['HSE06'][1]) + "): "
print str(EBsiteMatch['Columbite']['HSE06'][0]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Rutile']['PBE0'][1], EBsiteMatch['Rutile']['HSE06'][0], 'k-')
plt.plot(EBsiteMatch['Anatase']['PBE0'][1], EBsiteMatch['Anatase']['PBE0'][0],
'ro-', label = r'$\Delta$E$_{R-A,PBE0}$')
plt.plot(EBsiteMatch['Anatase']['HSE06'][1], EBsiteMatch['Anatase']['HSE06'][0],
'bo-', label = r'$\Delta$E$_{R-A,HSE06}$')
plt.plot(EBsiteMatch['Columbite']['HSE06'][1],
EBsiteMatch['Columbite']['HSE06'][0], 'go-', label = r'$\Delta$E$_{R-C,HSE06}$')
plt.xlabel( 'Exact Exchange Fraction (%)' )
plt.ylim( (-0.1, 0.05) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'lower right', prop={'size':9.5})
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-stability-RAC-HSE06PBE0' + '.' + ext, dpi=300)
plt.clf()
As was shown in Section Standard vs. Self-consistent Linear Response, corresponding PBE functional results for the Rutile, Anatase, Columbite, and Brookite TiO2 polymorphs were all equal to 3.0 eV within uncertainty and generally closely encompassed this value. Consistent with these results, those shown in the PBE, standard PAW pseudopotential subsection of Section Plot Generation: Hubbard U, and those resolved in previous work, cite:Finazzi2008,Morgan2010 a Rutile-Anatase formation energy of approximately +0.005 eV/f.u. TiO2 can be predictively estimated. This result is in very strong agreement with that derived by Rao et al. cite:Rao1961,Dompablo2011 When this Rutile-Anatase formation energy is transferred from predictive Hubbard
Validation of this result proceeds from performing several variations on the calculations presented above, such as uniformly changing their
\begin{flalign*}
&ΠMorph, NKRED, EOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(NKRED=’1’) or (NKRED=’2’)] \
&σ ∧ [(Functional=’HSE06’)] \
&σ ∧ [(FractionHF=’0.250’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEHFRange \
& \Join CalcResultEnergetics \Join ParametersInputVASP]
\end{flalign*}
#+caption: MySQL query HSE06 NKRED Rutile Anatase
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, NKRED, Eopt, atomnum = [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
EtoNKRED_dict = {}
for row in db.execute('''
select distinct s.Morph, input.NKRED, cre.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEHFRange as feh on feh.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feh.EOpt
inner join ParametersInputVASP as input on input.Energy=cre.Energy
where (s.Morph='Rutile' or s.Morph='Anatase')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and (input.NKRED='1' or input.NKRED='2')
and mt.Functional='HSE06'
and feh.FractionHF='0.25'
;'''):
a,b,c,d = row
try:
EtoNKRED_dict[ (a,b) ].append( [c,d] )
except KeyError:
EtoNKRED_dict[ (a,b) ] = []
EtoNKRED_dict[ (a,b) ].append( [c,d] )
E_NKRED2 = ( EtoNKRED_dict[ ('Anatase','2') ][0][0] / EtoNKRED_dict[ ('Anatase','2') ][0][1]
) - ( EtoNKRED_dict[ ('Rutile','2') ][0][0] / EtoNKRED_dict[ ('Anatase','2') ][0][1] )
E_NKRED1 = ( EtoNKRED_dict[ ('Anatase','1') ][0][0] / EtoNKRED_dict[
('Anatase','1') ][0][1]
) - ( EtoNKRED_dict[ ('Rutile','1') ][0][0] / EtoNKRED_dict[ ('Anatase','1') ][0][1] )
print "HSE06 Functional (a = 0.25):"
print "E(R-A)(NKRED = 2) = " + str(E_NKRED2)
print "E(R-A)(NKRED = 1) = " + str(E_NKRED1) + "\n"Further validation of this result proceeds from investigating the discontinuity of the Rutile-Anatase formation energy in the limit of exact HF exchange (
\begin{flalign*}
&ΠMorph, FractionHF, EOpt, Stoich1
&σ ∧ [(Morph=’Rutile’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(NKRED=’2’)] \
&σ ∧ [(Functional=’HSE06’)] \
&σ ∧ [(FractionHF=’0.250’) ∨ (FractionHF=’0.500’) ∨ (FractionHF=’0.750’) \
&σ ∨ (FractionHF=’0.825’) ∨ (FractionHF=’0.875’) ∨ (FractionHF=’0.950’) \
&σ ∨ (FractionHF=’0.975’) ∨ (FractionHF=’0.990’) ∨ (FractionHF=’1.000’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEHFRange \
& \Join CalcResultEnergetics \Join ParametersInputVASP]
\end{flalign*}
#+caption: MySQL query HSE06 AEXX limit Rutile
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, FractionHF, Eopt, atomnum = [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select distinct s.Morph, feh.FractionHF, feh.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEHFRange as feh on feh.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feh.EOpt
inner join ParametersInputVASP as input on input.Energy=cre.Energy
where s.Morph='Rutile'
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and input.NKRED='2'
and mt.Functional='HSE06'
and (feh.FractionHF='0.250' or feh.FractionHF='0.500' or feh.FractionHF='0.750'
or feh.FractionHF='0.825' or feh.FractionHF='0.875' or feh.FractionHF='0.950'
or feh.FractionHF='0.975' or feh.FractionHF='0.99' or feh.FractionHF='1.000')
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
Polymorph += [a]
FractionHF += [b]
datapts_dict[row] = datapts_list
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Morph_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
del_list = []
for k in datapts_dict:
if k[0] == i:
SortElist[ k[0] ][ k[1] ] = float( k[2] ) / float( k[3] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in Morph_list:
HF_list = SortElist[i].keys()
HF_list.sort()
print "HSE06 Functional: \n"
for k in HF_list:
print "E(R, a = " + str(k) + ") = " + str(SortElist[i][k])
In the previous section, Section Sample Input Files: Hybrid Functionals, the value of the fraction of exact exchange (
\begin{flalign*}
&ΠMorph, UValue, VOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Columbite’) \
&σ ∧ (Morph=’Anatase’) ∨ (Morph=’Brookite’) \
&σ ∨ (Morph=’Baddeleyite’) ∨ (Morph=’Cotunnite’) \
&σ ∨ (Morph=’Fluorite’) ∨ (Morph=’Pyrite’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query Volumetric Expansion Hubbard U
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, Uvalue, Vopt, atomnum = [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select s.Morph, feu.UValue, feu.VOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Columbite' or s.Morph='Anatase'
or s.Morph='Brookite' or s.Morph='Baddeleyite' or
s.Morph='Cotunnite' or s.Morph='Fluorite' or s.Morph='Pyrite')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Functional='PBE'
and mt.Software='VASP'
and mt.Method='E'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
Polymorph += [a]
datapts_dict[row] = datapts_list
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in Morph_list:
SortElist[i] = {}
EBsiteMatch[i] = []
del_list = []
for j in datapts_dict:
if j[0] == i:
SortElist[ j[0] ][ j[1] ] = float( j[2] )
del_list.append( datapts_dict[j] )
else:
pass
for l in del_list:
del l
for i in Morph_list:
U_list = SortElist[i].keys()
U_list.sort()
VmapU_list = []
for j in U_list:
VmapU_value = SortElist[i][j] / SortElist[i][0]
VmapU_list.append( VmapU_value )
EBsiteMatch[i].append( VmapU_list )
EBsiteMatch[i].append( U_list )
ax = plt.gca()
print "PBE Functional, Ti O pseudopotentials:"
print "Rutile: " + str(EBsiteMatch['Rutile'][0][0:2])
print str(EBsiteMatch['Rutile'][0][2:4])
print str(EBsiteMatch['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['Anatase'][0][0:2])
print str(EBsiteMatch['Anatase'][0][2:4])
print str(EBsiteMatch['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['Columbite'][0][0:2])
print str(EBsiteMatch['Columbite'][0][2:4])
print str(EBsiteMatch['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['Brookite'][0][0:2])
print str(EBsiteMatch['Brookite'][0][2:4])
print str(EBsiteMatch['Brookite'][0][4:])
print "Baddeleyite: " + str(EBsiteMatch['Baddeleyite'][0][0:2])
print str(EBsiteMatch['Baddeleyite'][0][2:4])
print str(EBsiteMatch['Baddeleyite'][0][4:])
print "Cotunnite: " + str(EBsiteMatch['Cotunnite'][0][0:2])
print str(EBsiteMatch['Cotunnite'][0][2:4])
print str(EBsiteMatch['Cotunnite'][0][4:])
print "Fluorite: " + str(EBsiteMatch['Fluorite'][0][0:2])
print str(EBsiteMatch['Fluorite'][0][2:4])
print str(EBsiteMatch['Fluorite'][0][4:])
print "Pyrite: " + str(EBsiteMatch['Pyrite'][0][0:2])
print str(EBsiteMatch['Pyrite'][0][2:4])
print str(EBsiteMatch['Pyrite'][0][4:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Rutile'][1], EBsiteMatch['Rutile'][0], 'ko-')
plt.plot(EBsiteMatch['Anatase'][1], EBsiteMatch['Anatase'][0],
'bo-', label = r'$\Delta$V$_{R-A}$')
plt.plot(EBsiteMatch['Columbite'][1], EBsiteMatch['Columbite'][0],
'ro-', label = r'$\Delta$V$_{R-C}$')
plt.plot(EBsiteMatch['Brookite'][1], EBsiteMatch['Brookite'][0],
'go-', label = r'$\Delta$V$_{R-B}$')
plt.plot(EBsiteMatch['Baddeleyite'][1], EBsiteMatch['Baddeleyite'][0],
'mo-', label = r'$\Delta$V$_{R-Ba}$')
plt.plot(EBsiteMatch['Cotunnite'][1], EBsiteMatch['Cotunnite'][0],
'co--', label = r'$\Delta$V$_{R-Co}$')
plt.plot(EBsiteMatch['Fluorite'][1], EBsiteMatch['Fluorite'][0],
'yo-', label = r'$\Delta$V$_{R-F}$')
plt.plot(EBsiteMatch['Pyrite'][1], EBsiteMatch['Pyrite'][0],
'co-', label = r'$\Delta$V$_{R-P}$')
plt.xlabel( 'U value (eV)' )
plt.ylim( (1, 1.16) )
plt.ylabel('Volume Expansion Ratio (V$_{0}$[$U$]/V$_{0}$[$U$ = 0])')
plt.legend(loc = 'upper center', prop={'size':8.75}, ncol = 2)
plt.gcf().subplots_adjust(left=0.25)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-volume-RACBBaCoFP-PBE' + '.' + ext, dpi=300)
plt.clf()
In the hybrid functional calculations performed in this article, a noteworthy discontinuity in the trend depicting the relationship between Rutile-Anatase formation energy vs.
\begin{flalign*}
&ΠMorph, Functional, FractionHF, VOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(NKRED=’2’)] \
&σ ∧ [(Functional=’HSE06’) ∨ (Functional=’PBE0’)] \
&σ ∧ [(FractionHF=’0.250’) ∨ (FractionHF=’0.500’) ∨ (FractionHF=’0.750’) \
&σ ∨ (FractionHF=’0.825’) ∨ (FractionHF=’0.875’) ∨ (FractionHF=’0.950’) \
&σ ∨ (FractionHF=’1.000’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEHFRange \
& \Join CalcResultEnergetics \Join ParametersInputVASP]
\end{flalign*}
\begin{flalign*}
&ΠMorph, UValue, VOpt, Stoich1
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(UValue=’0’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query Rutile vs. Anatase Hybrid Functional Volumetric Expansion
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
Polymorph, Functional, FractionHF, Vopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select distinct s.Morph, mt.Functional, feh.FractionHF, feh.VOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEHFRange as feh on feh.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feh.EOpt
inner join ParametersInputVASP as input on input.Energy=cre.Energy
where (s.Morph='Rutile' or s.Morph='Anatase')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and input.NKRED='2'
and (mt.Functional='HSE06' or mt.Functional='PBE0')
and (feh.FractionHF='0.250' or feh.FractionHF='0.500' or feh.FractionHF='0.750'
or feh.FractionHF='0.825' or feh.FractionHF='0.875'
or feh.FractionHF='0.950' or feh.FractionHF='1.000')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
Polymorph += [a]
Functional += [b]
datapts_dict[row] = datapts_list
Morph_list = list( set( Polymorph ) )
Functional_list = list( set( Functional ) )
SortElist = {}
for i in Morph_list:
SortElist[i] = {}
for j in Functional_list:
SortElist[i][j] = {}
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
if not SortElist[i][j]:
del SortElist[i][j]
Polymorph, Uvalue, Eopt, atomnum = [], [], [], []
datapts_dict = {}
for row in db.execute('''
select s.Morph, feu.UValue, feu.VOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where (s.Morph='Rutile' or s.Morph='Anatase')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and mt.Functional='PBE'
and feu.UValue='0'
;'''):
datapts_list = []
a,b,c,d = row
datapts_list.append( (a,b,c,d) )
Polymorph += [a]
datapts_dict[row] = datapts_list
Morph_list_0 = list( set( Polymorph ) )
for i in Morph_list_0:
SortElist[i]['PBE0'][float(0.)] = []
SortElist[i]['HSE06'][float(0.)] = []
del_list = []
for j in datapts_dict:
if j[0] == i:
SortElist[ j[0] ]['PBE0'][float(0.)] = float( j[2] )
SortElist[ j[0] ]['HSE06'][float(0.)] = float( j[2] )
del_list.append( datapts_dict[j] )
else:
pass
for l in del_list:
del l
EBsiteMatch = {}
for i in Morph_list:
EBsiteMatch[i] = {}
for j in Functional_list:
EBsiteMatch[i][j] = []
HF_list = SortElist[i][j].keys()
HF_list.sort()
VmapHF_list = []
for k in HF_list:
VmapHF_value = SortElist[i][j][k] / SortElist[i][j][float(0.0)]
VmapHF_list.append( VmapHF_value )
EBsiteMatch[i][j].append( VmapHF_list )
EBsiteMatch[i][j].append( HF_list )
ax = plt.gca()
print "PBE0 Functional:"
print "Rutile: " + str(EBsiteMatch['Rutile']['PBE0'][0][0:2])
print str(EBsiteMatch['Rutile']['PBE0'][0][2:5])
print str(EBsiteMatch['Rutile']['PBE0'][0][5:])
print "Anatase: " + str(EBsiteMatch['Anatase']['PBE0'][0][0:2])
print str(EBsiteMatch['Anatase']['PBE0'][0][2:5])
print str(EBsiteMatch['Anatase']['PBE0'][0][5:]) + "\n"
print "HSE06 Functional:"
print "Rutile: " + str(EBsiteMatch['Rutile']['HSE06'][0][0:2])
print str(EBsiteMatch['Rutile']['HSE06'][0][2:5])
print str(EBsiteMatch['Rutile']['HSE06'][0][5:])
print "Anatase: " + str(EBsiteMatch['Anatase']['HSE06'][0][0:2])
print str(EBsiteMatch['Anatase']['HSE06'][0][2:5])
print str(EBsiteMatch['Anatase']['HSE06'][0][5:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Rutile']['PBE0'][1], EBsiteMatch['Rutile']['PBE0'][0],
'ko-', label = r'$\Delta$V$_{R,PBE0}$')
plt.plot(EBsiteMatch['Rutile']['HSE06'][1], EBsiteMatch['Rutile']['HSE06'][0],
'bo-', label = r'$\Delta$V$_{R,HSE06}$')
plt.plot(EBsiteMatch['Anatase']['PBE0'][1], EBsiteMatch['Anatase']['PBE0'][0],
'ro-', label = r'$\Delta$V$_{A,PBE0}$')
plt.plot(EBsiteMatch['Anatase']['HSE06'][1], EBsiteMatch['Anatase']['HSE06'][0],
'go-', label = r'$\Delta$V$_{A,HSE06}$')
plt.xlabel( 'Exact Exchange Fraction (%)' )
plt.ylim( (0.9, 1.0) )
plt.ylabel('Volume Expansion Ratio (V$_{0}$[%]/V$_{0}$[% = 0])')
plt.legend(loc = 'lower left', prop={'size':9})
plt.gcf().subplots_adjust(left=0.25)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-volume-RA-PBE0HSE06' + '.' + ext, dpi=300)
plt.clf()
In previous work, cite:Dompablo2011 the relationship between several structural factors and incrementation of the
\begin{flalign*}
&ΠEnergy, EOpt, HFSCREEN, AEXX, PrimVec1, PrimVec3
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’)] \
&σ ∧ [(At1=’Ti’)] \
&σ ∧ [(PseudoB=’PAW-Ti’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(NKRED=’2’)] \
&σ ∧ [(Functional=’HSE06’) ∨ (Functional=’PBE0’)] \
&σ ∧ [(FractionHF=’0.250’) ∨ (FractionHF=’0.500’) ∨ (FractionHF=’0.750’) \
&σ ∨ (FractionHF=’0.825’) ∨ (FractionHF=’0.875’) ∨ (FractionHF=’0.950’) \
&σ ∨ (FractionHF=’1.000’)] \
&[Structure \Join Composition1 \Join Metadata \Join FormEHFRange \
& \Join CalcResultEnergetics \Join ParametersInputVASP \Join ParametersPosFinal]
\end{flalign*}
#+caption: MySQL query Rutile vs. Anatase Hybrid Functional c/a Expansion
import matplotlib.pyplot as plt
import sqlite3
import numpy as np
db = sqlite3.connect('ITEOO_data.sqlite')
HFtype, FractionHF = [], []
EtoEOptHF_dict = {}
NULL_CHAR = '-'
SColon_CHAR = ';'
Morph_list = [ 'Rutile', 'Anatase' ]
WRT_Frac = 0.25
for row in db.execute('''
select distinct cre.Energy, cre.EOpt from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Composition2 as c2 on c1.MID=c2.MID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEHFRange as feh on feh.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feh.EOpt
inner join ParametersInputVASP as input on input.Energy=cre.Energy
where (s.Morph='Rutile' or s.Morph='Anatase')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and input.NKRED='2'
and (mt.Functional='HSE06' or mt.Functional='PBE0')
and (feh.FractionHF='0.250' or feh.FractionHF='0.500' or feh.FractionHF='0.750'
or feh.FractionHF='0.825' or feh.FractionHF='0.875'
or feh.FractionHF='0.950' or feh.FractionHF='1.000')
;'''):
a,b = row
try:
EtoEOptHF_dict[b].append(a)
except KeyError:
EtoEOptHF_dict[b] = []
EtoEOptHF_dict[b].append(a)
Emin_list = []
for i in EtoEOptHF_dict.keys():
Ediff_list = []
for j in EtoEOptHF_dict[i]:
Ediff_value = abs( float(i) - float(j) )
Ediff_list.append( Ediff_value )
Emindex = Ediff_list.index( min( Ediff_list ) )
Emin_list.append( EtoEOptHF_dict[i][ Emindex ] )
datapts_dict = {}
for row in db.execute('''
select distinct input.Energy, input.HFSCREEN, input.AEXX, pos.PrimVec1, pos.PrimVec3 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Composition2 as c2 on c1.MID=c2.MID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEHFRange as feh on feh.CID=mt.CID
inner join CalcResultEnergetics as cre on cre.EOpt=feh.EOpt
inner join ParametersPosFinalVASP as pos on pos.Energy=cre.Energy
inner join ParametersInputVASP as input on input.Energy=pos.Energy
where (s.Morph='Rutile' or s.Morph='Anatase')
and c1.At1='Ti'
and mt.PseudoB='PAW-Ti'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and mt.Method='E'
and input.NKRED='2'
and (mt.Functional='HSE06' or mt.Functional='PBE0')
and (feh.FractionHF='0.250' or feh.FractionHF='0.500' or feh.FractionHF='0.750'
or feh.FractionHF='0.825' or feh.FractionHF='0.875'
or feh.FractionHF='0.950' or feh.FractionHF='1.000')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
HFtype += [b]
FractionHF += [c]
if a in Emin_list:
datapts_dict[row] = datapts_list
HFtype_list = list( set( HFtype ) )
FractionHF_list = list( set( FractionHF ) )
SortElist = {}
EBsiteMatch = {}
HFlabel_list = []
for i in HFtype_list:
if i == NULL_CHAR:
HFtype_val = 'PBE0'
else:
HFtype_val = 'HSE06'
HFlabel_list.append( HFtype_val )
SortElist[ HFtype_val ] = {}
EBsiteMatch[ HFtype_val ] = {}
for j in Morph_list:
SortElist[ HFtype_val ][j] = {}
EBsiteMatch[ HFtype_val ][j] = []
del_list = []
for k in datapts_dict:
if k[1] == i:
a_string = str( k[3] ).split(SColon_CHAR)
c_string = str( k[4] ).split(SColon_CHAR)
a_vec = np.zeros( len(a_string) )
c_vec = np.zeros( len(c_string) )
for l in range( len(a_string) ):
a_vec[l] = float( a_string[l] )
c_vec[l] = float( c_string[l] )
a_dot = np.dot( a_vec, a_vec )
c_dot = np.dot( c_vec, c_vec )
ca_ratio = ( c_dot / a_dot )**(0.5)
if ca_ratio < 1.:
SortElist[ HFtype_val ]['Rutile'][ float(k[2]) ] = ca_ratio
else:
SortElist[ HFtype_val ]['Anatase'][ float(k[2]) ] = (2. * c_vec[2] ) / (
a_dot )**(0.5)
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in HFlabel_list:
for j in Morph_list:
HF_list = SortElist[i][j].keys()
camapHF_list = []
HF_list.sort()
for k in HF_list:
camapHF_value = SortElist[i][j][k] / SortElist[i][j][WRT_Frac]
camapHF_list.append( camapHF_value )
EBsiteMatch[i][j].append( camapHF_list )
EBsiteMatch[i][j].append( HF_list )
ax = plt.gca()
print "PBE0 Functional:"
print "Rutile: " + str(EBsiteMatch['PBE0']['Rutile'][0][0:2])
print str(EBsiteMatch['PBE0']['Rutile'][0][2:4])
print str(EBsiteMatch['PBE0']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['PBE0']['Anatase'][0][0:2])
print str(EBsiteMatch['PBE0']['Anatase'][0][2:4])
print str(EBsiteMatch['PBE0']['Anatase'][0][4:])+ "\n"
print "HSE06 Functional:"
print "Rutile: " + str(EBsiteMatch['HSE06']['Rutile'][0][0:2])
print str(EBsiteMatch['HSE06']['Rutile'][0][2:4])
print str(EBsiteMatch['HSE06']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['HSE06']['Anatase'][0][0:2])
print str(EBsiteMatch['HSE06']['Anatase'][0][2:4])
print str(EBsiteMatch['HSE06']['Anatase'][0][4:]) + "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['PBE0']['Rutile'][1], EBsiteMatch['PBE0']['Rutile'][0],
'ko-', label = r'$\Delta$c/a$_{R,PBE0}$')
plt.plot(EBsiteMatch['HSE06']['Rutile'][1], EBsiteMatch['HSE06']['Rutile'][0],
'bo-', label = r'$\Delta$c/a$_{R,HSE06}$')
plt.plot(EBsiteMatch['PBE0']['Anatase'][1], EBsiteMatch['PBE0']['Anatase'][0],
'ro-', label = r'$\Delta$c/a$_{A,PBE0}$')
plt.plot(EBsiteMatch['HSE06']['Anatase'][1], EBsiteMatch['HSE06']['Anatase'][0],
'go-', label = r'$\Delta$c/a$_{A,HSE06}$')
plt.xlabel( 'Exact Exchange Fraction (%)' )
labels = [ WRT_Frac, 0.50, 0.75, 1.00 ]
plt.xticks(labels)
plt.ylim( (0.985, 1.020) )
plt.ylabel('c/a Ratio (c/a[%]/c/a[% = ' + str(WRT_Frac) + '])')
plt.legend(loc = 'upper left', prop={'size':8.5})
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-caratio-RA-PBE0HSE06' + '.' + ext, dpi=300)
plt.clf()
As shown in the article and Section Plot Generation: Hybrid Functionals, the PBE0 and HSE06 functionals approach respective Rutile-Anatase formation energy maxima of approximately +0.045 and +0.030 eV/f.u. TiO2 as
In extending this error analysis to other BO2 (B = Ti, V, Ru, Ir) systems, pertinent formation energetic results cite:Mehta2014 for the TiO2, VO2, RuO2, and IrO2 systems are queried, plotted, and analyzed below. Though these plots are assembled as a function of
\begin{flalign*}
&ΠAt1, Morph, UValue, EOpt, Stoich1
&σ ∧ [(Functional=’PBE’)] \
&σ ∧ [(PseudoO=’PAW-O’)] \
&σ ∧ [(Software=’VASP’)] \
&σ ∧ [(Method=’E’)] \
&σ ∧ [(UValue=’0’) ∨ (UValue=’1’) ∨ (UValue=’2’) ∨ (UValue=’3’) \
&σ ∨ (UValue=’4’) ∨ (UValue=’5’) ∨ (UValue=’6’)] \
&σ ∧ [[(PseudoB=’PAW-Ti’) ∨ (PseudoB=’PAW-V’)] \
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Anatase’) \
&σ ∨ (Morph=’Columbite’) ∨ (Morph=’Brookite’)]] \
&σ ∧ [[(PseudoB=’PAW-Ru’) ∨ (PseudoB=’PAW-Ir’)] \
&σ ∧ [(Morph=’Rutile’) ∨ (Morph=’Pyrite’) ∨ (Morph=’Columbite’)]] \
&[Structure \Join Composition1 \Join Metadata \Join FormEURange]
\end{flalign*}
#+caption: MySQL query Epitaxial Comparison Ti V Ru Ir polymorph
import matplotlib.pyplot as plt
import sqlite3
db = sqlite3.connect('ITEOO_data.sqlite')
atom1, Polymorph, Uvalue, Eopt, atomnum = [], [], [], [], []
datapts_dict = {}
WRT_Morph = 'Rutile'
for row in db.execute('''
select c1.At1, s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where mt.Functional='PBE'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
and (mt.PseudoB='PAW-Ti' or mt.PseudoB='PAW-V')
and (s.Morph='Rutile' or s.Morph='Anatase'
or s.Morph='Columbite' or s.Morph='Brookite')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
atom1 += [a]
Polymorph += [b]
datapts_dict[row] = datapts_list
At1_list = list( set( atom1 ) )
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in At1_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Morph_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in At1_list:
for j in Morph_list:
U_list = SortElist[i][j].keys()
U_list.sort()
EmapU_list = []
for k in U_list:
EmapU_value = SortElist[i][j][k] - SortElist[i][WRT_Morph][k]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j].append( EmapU_list )
EBsiteMatch[i][j].append( U_list )
for row in db.execute('''
select c1.At1, s.Morph, feu.UValue, feu.EOpt, c1.Stoich1 from Structure as s
inner join Composition1 as c1 on c1.SID=s.SID
inner join Metadata as mt on (mt.MID=c1.MID and mt.SID=c1.SID)
inner join FormEURange as feu on feu.CID=mt.CID
where mt.Functional='PBE'
and mt.PseudoO='PAW-O'
and mt.Software='VASP'
and (feu.UValue='0' or feu.UValue='1' or feu.UValue='2' or feu.UValue='3'
or feu.UValue='4' or feu.UValue='5' or feu.UValue='6')
and mt.PseudoB='PAW-Ru' or mt.PseudoB='PAW-Ir'
and (s.Morph='Rutile' or s.Morph='Pyrite'
or s.Morph='Columbite')
;'''):
datapts_list = []
a,b,c,d,e = row
datapts_list.append( (a,b,c,d,e) )
atom1 += [a]
Polymorph += [b]
datapts_dict[row] = datapts_list
At1_list = list( set( atom1 ) )
Morph_list = list( set( Polymorph ) )
SortElist = {}
EBsiteMatch = {}
for i in At1_list:
SortElist[i] = {}
EBsiteMatch[i] = {}
for j in Morph_list:
SortElist[i][j] = {}
EBsiteMatch[i][j] = []
del_list = []
for k in datapts_dict:
if k[0] == i and k[1] == j:
SortElist[ k[0] ][ k[1] ][ k[2] ] = float( k[3] ) / float( k[4] )
del_list.append( datapts_dict[k] )
else:
pass
for l in del_list:
del l
for i in At1_list:
for j in Morph_list:
U_list = SortElist[i][j].keys()
U_list.sort()
EmapU_list = []
for k in U_list:
EmapU_value = SortElist[i][j][k] - SortElist[i][WRT_Morph][k]
EmapU_list.append( EmapU_value )
EBsiteMatch[i][j].append( EmapU_list )
EBsiteMatch[i][j].append( U_list )
ax = plt.gca()
print r'PBE Functional, TiO2:'
print "Rutile: " + str(EBsiteMatch['Ti']['Rutile'][0][0:2])
print str(EBsiteMatch['Ti']['Rutile'][0][2:4])
print str(EBsiteMatch['Ti']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['Ti']['Anatase'][0][0:2])
print str(EBsiteMatch['Ti']['Anatase'][0][2:4])
print str(EBsiteMatch['Ti']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['Ti']['Columbite'][0][0:2])
print str(EBsiteMatch['Ti']['Columbite'][0][2:4])
print str(EBsiteMatch['Ti']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['Ti']['Brookite'][0][0:2])
print str(EBsiteMatch['Ti']['Brookite'][0][2:4])
print str(EBsiteMatch['Ti']['Brookite'][0][4:]) + "\n"
print r'PBE Functional, VO2:'
print "Rutile: " + str(EBsiteMatch['V']['Rutile'][0][0:2])
print str(EBsiteMatch['V']['Rutile'][0][2:4])
print str(EBsiteMatch['V']['Rutile'][0][4:])
print "Anatase: " + str(EBsiteMatch['V']['Anatase'][0][0:2])
print str(EBsiteMatch['V']['Anatase'][0][2:4])
print str(EBsiteMatch['V']['Anatase'][0][4:])
print "Columbite: " + str(EBsiteMatch['V']['Columbite'][0][0:2])
print str(EBsiteMatch['V']['Columbite'][0][2:4])
print str(EBsiteMatch['V']['Columbite'][0][4:])
print "Brookite: " + str(EBsiteMatch['V']['Brookite'][0][0:2])
print str(EBsiteMatch['V']['Brookite'][0][2:4])
print str(EBsiteMatch['V']['Brookite'][0][4:])+ "\n"
print r'PBE Functional, RuO2:'
print "Rutile: " + str(EBsiteMatch['Ru']['Rutile'][0][0:2])
print str(EBsiteMatch['Ru']['Rutile'][0][2:4])
print str(EBsiteMatch['Ru']['Rutile'][0][4:])
print "Pyrite: " + str(EBsiteMatch['Ru']['Pyrite'][0][0:2])
print str(EBsiteMatch['Ru']['Pyrite'][0][2:4])
print str(EBsiteMatch['Ru']['Pyrite'][0][4:])
print "Columbite: " + str(EBsiteMatch['Ru']['Columbite'][0][0:2])
print str(EBsiteMatch['Ru']['Columbite'][0][2:4])
print str(EBsiteMatch['Ru']['Columbite'][0][4:]) + "\n"
print r'PBE Functional, IrO2:'
print "Rutile: " + str(EBsiteMatch['Ir']['Rutile'][0][0:2])
print str(EBsiteMatch['Ir']['Rutile'][0][2:4])
print str(EBsiteMatch['Ir']['Rutile'][0][4:])
print "Pyrite: " + str(EBsiteMatch['Ir']['Pyrite'][0][0:2])
print str(EBsiteMatch['Ir']['Pyrite'][0][2:4])
print str(EBsiteMatch['Ir']['Pyrite'][0][4:])
print "Columbite: " + str(EBsiteMatch['Ir']['Columbite'][0][0:2])
print str(EBsiteMatch['Ir']['Columbite'][0][2:4])
print str(EBsiteMatch['Ir']['Columbite'][0][4:])+ "\n"
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Ti']['Rutile'][1], EBsiteMatch['Ti']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['Ti']['Anatase'][1], EBsiteMatch['Ti']['Anatase'][0],
'bo-', label = r'$\Delta$E$_{R-A}$')
plt.plot(EBsiteMatch['Ti']['Columbite'][1], EBsiteMatch['Ti']['Columbite'][0],
'ro-', label = r'$\Delta$E$_{R-C}$')
plt.plot(EBsiteMatch['Ti']['Brookite'][1], EBsiteMatch['Ti']['Brookite'][0],
'go-', label = r'$\Delta$E$_{R-B}$')
plt.text(0.595, 2.75, 'B = Ti',
style = 'italic', color = 'white', fontsize = 24,
transform = ax.transAxes,
bbox={'facecolor':'b', 'alpha':0.5, 'pad':5})
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 0.1) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper left', prop={'size':10})
plt.gcf().subplots_adjust(left=0.27)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/TiO2-epitax-RACB-PBE' + '.' + ext, dpi=300)
plt.clf()
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['V']['Rutile'][1], EBsiteMatch['V']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['V']['Anatase'][1], EBsiteMatch['V']['Anatase'][0],
'bo-', label = r'$\Delta$E$_{R-A}$')
plt.plot(EBsiteMatch['V']['Columbite'][1], EBsiteMatch['V']['Columbite'][0],
'ro-', label = r'$\Delta$E$_{R-C}$')
plt.plot(EBsiteMatch['V']['Brookite'][1], EBsiteMatch['V']['Brookite'][0],
'go-', label = r'$\Delta$E$_{R-B}$')
plt.text(0.565, 2.75, 'B = V',
style = 'italic', color = 'white', fontsize = 24,
transform = ax.transAxes,
bbox={'facecolor':'b', 'alpha':0.5, 'pad':5})
plt.xlabel( 'U value (eV)' )
plt.ylim( (0.0, 0.8) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper left', prop={'size':10})
plt.gcf().subplots_adjust(left=0.23)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/VO2-epitax-RACB-PBE' + '.' + ext, dpi=300)
plt.clf()
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Ru']['Rutile'][1], EBsiteMatch['Ru']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['Ru']['Pyrite'][1], EBsiteMatch['Ru']['Pyrite'][0],
'ro-', label = r'$\Delta$E$_{R-P}$')
plt.plot(EBsiteMatch['Ru']['Columbite'][1], EBsiteMatch['Ru']['Columbite'][0],
'bo-', label = r'$\Delta$E$_{R-C}$')
plt.text(0.54, 2.75, 'B = Ru',
style = 'italic', color = 'white', fontsize = 24,
transform = ax.transAxes,
bbox={'facecolor':'b', 'alpha':0.5, 'pad':5})
plt.xlabel( 'U value (eV)' )
plt.ylim( (-0.1, 0.5) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper left', prop={'size':10})
plt.gcf().subplots_adjust(left=0.25)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/RuO2-epitax-RPC-PBE' + '.' + ext, dpi=300)
plt.clf()
plt.figure(figsize=(3,4))
plt.plot(EBsiteMatch['Ir']['Rutile'][1], EBsiteMatch['Ir']['Rutile'][0], 'k-')
plt.plot(EBsiteMatch['Ir']['Pyrite'][1], EBsiteMatch['Ir']['Pyrite'][0],
'ro-', label = r'$\Delta$E$_{R-P}$')
plt.plot(EBsiteMatch['Ir']['Columbite'][1], EBsiteMatch['Ir']['Columbite'][0],
'bo-', label = r'$\Delta$E$_{R-C}$')
plt.text(0.580, 2.75, 'B = Ir',
style = 'italic', color = 'white', fontsize = 24,
transform = ax.transAxes,
bbox={'facecolor':'b', 'alpha':0.5, 'pad':5})
plt.xlabel( 'U value (eV)' )
plt.ylim( (0.2, 0.45) )
plt.ylabel('Energy Difference (eV/f.u.)')
plt.legend(loc = 'upper left', prop={'size':10})
plt.gcf().subplots_adjust(left=0.25)
plt.gcf().subplots_adjust(bottom=0.11)
for ext in ['png', 'pdf', 'eps']:
plt.savefig('./figures/IrO2-epitax-RPC-PBE' + '.' + ext, dpi=300)
plt.clf()
For the TiO2 system depicted above, Rutile, Anatase, Brookite, and Columbite polymorphs are featured in formation energy calculations set with respect to Rutile. Given an epitaxial stabilization window of 0.1-0.2 eV, stabilization of Anatase, Columbite, and Brookite polymorphs from a Rutile precursor is highly possible in the TiO2 system. The errors represented above, which range from 0.068-0.127 eV, are within an order of magnitude of the epitaxial stabilization range. When using the
For the VO2 system depicted above, Rutile, Anatase, Brookite, and Columbite polymorphs are featured in formation energy calculations set with respect to Rutile. Given an epitaxial stabilization window of 0.1-0.2 eV and the errors mentioned above, stabilization of Anatase, Columbite, and Brookite polymorphs from a Rutile precursor is possible in the VO2 system and error in calculations involving them can affect their predicted energetic stabilities. For VO2, an appropriate selection of a
For the RuO2 system depicted above, Rutile, Pyrite, and Columbite polymorphs are featured in formation energy calculations set with respect to Rutile. Given an epitaxial stabilization window of 0.1-0.2 eV and the errors mentioned previously, stabilization of the Columbite polymorph from a Rutile precursor is arguably possible in the RuO2 system. For RuO2, an appropriate selection of a
For the IrO2 system depicted above, Rutile, Pyrite, and Columbite polymorphs are featured in formation energy calculations set with respect to Rutile. Given a value of
bibliography:shorttitles.bib,references.bib