What can one learn about material structure given a single first-principles calculation?

TL;DR

This paper demonstrates that a variable derived from first-principles electron orbitals in high-symmetry structures can predict the true ground state structure of complex oxides, aiding in material design.

Contribution

It introduces a method to extract a hidden variable from electron orbitals that predicts material ground states without structural distortions.

Findings

Variable X correlates with ground state structures.

Wannier functions can extract X for structure prediction.

Method applicable to various complex oxides.

Abstract

We extract a variable $X$ from electron orbitals $\Psi_{n\bf{k}}$ and energies $E_{n\bf{k}}$ in the parent high-symmetry structure of a wide range of complex oxides: perovskites, rutiles, pyrochlores, and cristobalites. Even though calculation was done only in the parent structure, with no distortions, we show that $X$ dictates material's true ground state structure. We propose using Wannier functions to extract concealed variables such as $X$ both for material structure prediction and for high-throughput approaches.

Figures

Tables

Table 1. Table 1: Variable X 𝑋 X calculated in a high-symmetry parent state correlates well with the preferred ground state structure (third subcolumn). The meanings of the crystal structure abbreviations are given at the bottom of the table.

Perovskite			Rutile			Pyrochlore			Cristobalite
ABO3	X	Struc.	AO2	X	Struc.	A2B2O7	X	Struc.	AO2	X	Struc.
	(eV)			(eV)			(eV)		ABO4	(eV)
CaZrO3	$0.47$	np	ReO2	$-0.27$	m	Mg2P2O7	$-1.47$	th	InSbO4	$3.62$	a
SrZrO3	$0.57$	np	WO2	$-0.26$	m	Mg2As2O7	$-0.96$	th	SnO2	$3.68$	a
CaHfO3	$0.66$	np	MoO2	$-0.18$	m	Y2Si2O7	$-0.73$	th	PbO2	$3.80$	a
PbZrO3	$0.70$	np	TcO2	$-0.09$	m	In2Si2O7	$-0.60$	th	AlSbO4	$3.95$	a
BaZrO3	$0.76$	np	NbO2	$-0.06$	n	Cd2V2O7	$-0.59$	th	GaSbO4	$3.98$	a
SrHfO3	$0.76$	np	VO2	$0.11$	m	Sc2Si2O7	$-0.54$	th	InPO4	$4.00$	c
PbHfO3	$0.90$	np	TiO2	$0.17$	r	In2Ge2O7	$-0.31$	th	InAsO4	$4.15$	c
CdTiO3	$0.90$	p	CrO2	$0.22$	r	La2Sn2O7	$-0.27$	py	GaPO4	$4.43$	c
CaTiO3	$0.91$	np	MnO2	$0.27$	r	Sc2Ge2O7	$-0.22$	th	AlPO4	$4.47$	c
BaHfO3	$0.96$	np	RuO2	$0.27$	r	Y2V2O7	$-0.18$	py	SiO2	$4.52$	c
SrTiO3	$1.03$	p	OsO2	$0.39$	r	Bi2Ti2O7	$-0.14$	py	GaAsO4	$4.54$	c
PbTiO3	$1.20$	p	IrO2	$0.45$	r	Y2Mo2O7	$-0.10$	py	AlAsO4	$4.56$	c
BaTiO3	$1.25$	p	SnO2	$0.77$	r	Y2Ti2O7	$-0.04$	py	GeO2	$4.67$	c
NaNbO3	$1.31$	p	PbO2	$0.82$	r	La2Hf2O7	$-0.04$	py
KNbO3	$1.45$	p	SiO2	$1.00$	r	La2Pb2O7	$-0.03$	py
LaAlO3	$1.54$	np	GeO2	$1.06$	r	Y2Sn2O7	$-0.02$	py
KTaO3	$1.66$	p				La2Zr2O7	$0.01$	py
						Bi2Hf2O7	$0.02$	py
np = non-polar hettotype			r = rutile hettotype			py = pyrochlore			c = cristobalite
p = polar hettotype			m = manganite hettotype			th = thortveitite			a = anatase
			n = NbO2 hettotype

Equations

$$W_n(\mathbf{r})=\frac{1}{N_{\mathbf{k}}}\sum _k{U}_{nm\mathbf{k}}{\Psi }_{m\mathbf{k}}(\mathbf{r})$$

$$H_{nm\mathbf{R}}=\text{matrixelement}{W}_{n}H{W}_{m\mathbf{R}}.$$

$$H_{nn0}-H_{mm0}$$

$$X=\text{matrixelement}{W}_{\perp }H{W}_{\perp }-\text{matrixelement}{W}_{\parallel }H{W}_{\parallel }.$$

Full text

I Introduction and motivation

There is a growing need for machines to learn about matter from results of a large set of first-principles computer calculations. However, a single such calculation produces $\sim 10^{6}$–$10^{9}$ bits of information and it is unclear how machine can make use of them.Curtarolo et al. (2013); Ghiringhelli et al. (2015) Therefore, so far, machines set aside most of these bits — such as those describing electron orbitals $\Psi_{n\bf{k}}$ — and try to learn from the simplest calculated quantities such as the total energy $E^{\rm tot}$ or similar.

One such example where machines are learning from $E^{\rm tot}$ are methods used for predicting positions of atoms in a yet unsynthesized material. These methods look for a set of atom positions $\xi$ that minimize the total energy $E^{\rm tot}(\xi)$ of the material. Total energy of a periodic solid can be computed from first-principles using various approximations to the density functional theory.Hohenberg and Kohn (1964); Kohn and Sham (1965) However, minimization of $E^{\rm tot}(\xi)$ still remains an unsolved problem as $\xi$ is a vector in a very highly dimensional space (it has $\approx 3N$ dimensions and $N$ is a number of atoms in the material). Nevertheless, this problem can be addressed heuristically using machine learning techniques such as evolutionary algorithmsLyakhov et al. (2013) or particle swarm optimization.Wang et al. (2012) Broadly speaking these methods first use $E^{\rm tot}$ calculated for a wide range of different candidate structures $\xi_{1},\xi_{2},...\xi_{m}$, but the same chemical composition, to learn about the underlying interactions in that material. Next, given this knowledge, the machine makes an informed guess for the next structure $\xi_{m+1}$ and the process repeats until an optimal structure is found. Clearly, if one used in this process information contained in electron orbitals $\Psi_{n\bf{k}}$, and not only $E^{\rm tot}$, one could make a more informed guess of structure $\xi_{m+1}$.

The need for a machine to learn in the context of materials science is also relevant for the so-called high-throughput approaches such as the materials project,Jain et al. (2011) the aflow,Curtarolo et al. (2012) the oqmd,Saal et al. (2013), aiida,Pizzi et al. (2016) or the nomadnom materials databases. While in the structure prediction problem one considers a single chemical composition at the time, in the high-throughput approach one wishes to learn about materials with a wide range of chemical compositions. Since these databases contain $\sim 10^{5}$–$10^{6}$ materials their total information content, if one were to store electron orbitals, is about $10^{14}$–$10^{15}$ bits.

In this paper we do not focus on what or how machines can learn from total energy $E^{\rm tot}$. Instead, the goal of this paper is to construct a descriptor of electron orbitals $\Psi_{n\bf{k}}$ and eigenenergies that can be used for learning in structure prediction and high-throughput approaches. As a proof of principle, we construct here a descriptor — denoted as $X$ — that is strongly correlated with the preferred crystal structure of a material, as described later. Since in our proof of principle work $X$ turns out to be a single number, one can establish a correlation between $X$ and structure just by inspection, without using machine learning. However, for materials with lower symmetry or in the cases of other properties of interest (i.e., not crystal structure) similarly constructed descriptors of electron orbitals will correspond to more than one number and one would therefore have to use machine learning.

Crystal structures of materials can be divided somewhat loosely into families of structuresMuller (2011) based on polyhedral units present in the structure and their connectivity. Each structure family is derived from a simple high-symmetry structure also called parent, aristotype, or prototype structure. Remaining structures in the family are then derived from the parent by either displacing or substituting atoms.Megaw (1973); Barnighausen (1980) These derivative structures are also called hettotypes. When structures are related to each other by distortion that does not preserve polyhedral units or their connectivity, we refer to them only as polymorphs.

This paper is structured as follows. In Sec. II we describe our approach and in Sec. III we present and discuss our results. A comparison of our results with atomic descriptors (Pettifor maps) is done in Sec. III.1.

II Approach

In what follows we present a general approach to extract a variable $X$ — given electron orbitals $\Psi_{n\bf{k}}$ and eigenenergies of a parent structure alone — that dictates which structure is a true ground state structure of that material. Since such a variable knows about material’s structure (i) it has the potential to be used in machine learning approaches described earlier. To that end, a useful variable should also satisfy two additional criteria: (ii) $X$ should be represented with as few numbers as possible, so that learning is as efficient as possible; and (iii) $X$ should be calculable automatically for a given material without any input from the user, so that it can be used by a machine.

To satisfy the criterion (i) we should ideally extract $X$ in some way from the fundamental variable of the problem. In the case of quantum mechanical descriptions of matter, the fundamental variable is the electron many-body wave function $\widetilde{\Psi}$. We note that the many-body wave function $\widetilde{\Psi}$ is formally not accessible in the commonly used Kohn-ShamKohn and Sham (1965) scheme for the density functional theory. Therefore we instead use here the so-called Kohn-Sham orbitals $\Psi_{n\bm{k}}$ (and corresponding eigenenergies $E_{n\bm{k}}$) as they are typically reasonably good approximationsHybertsen and Louie (1986) of the quasiparticle wave function and energy (here $n$ is a band index while $\bm{k}$ is a crystal momentum). Figure 1 shows some of the Kohn-Sham orbitals for one of the materials we studied in this work (SrTiO3 cubic perovskite).

Orbitals $\Psi_{n\bm{k}}$ themselves are not a good choice of $X$ as they do not satisfy criterion (ii). Namely, even in the simplest calculations (say, cubic silicon), a machine needs $\sim 10^{6}$ bits to describe $\Psi_{n\bm{k}}$ reasonably accurately. As discussed later, we circumvent this difficulty by defining $X$ in terms of the Wannier functions,Marzari et al. (2012) as they provide a compact and faithful representation of $\Psi_{n\bm{k}}$ (see Fig. 2). Wannier functions can be computed without user inputMustafa et al. (2015) so our approach satisfies condition (iii) as well.

To demonstrate the generality of our approach we consider here a diversified set of four structure families: ABO3 perovskites, AO2 rutiles, A2B2O7 pyrochlores, and AO2 cristobalites for a range of A and B anions. We chose these four families because they display different bonding environments containing both octahedral and tetrahedral oxygen units, as well as different connectivity of these units (see Fig. 3). The main result of this paper is that we find that our general approach yields an energy scale $X$ that is very well correlated with preferred ground state structure. We studied 64 representative compounds: 17 perovskites, 16 rutiles, 18 pyrochlores, and 13 cristobalites.

The variable $X$ is constructed in five steps. First, we relax the structure of each compound while preserving the symmetry of the parent structure (aristotype).111We use generalized gradient approximationPerdew et al. (1996) to the density functional theory within the quantum-espresso packageGiannozzi et al. (2009); Hamann (2013); Garrity et al. (2014) Wannier functions were constructed using Wannier90 package.Mostofi et al. (2008)

Second, we identify a set of electron bands of interest given the band structure from the first step. In our case we simply take as bands of interest a complex of bands with a dominant oxygen $p$-like character, as they are the most dispersive and thus contain the most information about the inter-atomic interactions in the material. Oxygen bands are also fully occupied and isolated in these materials, which simplifies our Wannier function based analysis. Figure 2 shows bands of interest in the representative case of each of the four families we studied.

Third, we convert electron orbitals $\Psi_{n\bm{k}}$ of bands selected in step two into a basis of maximally localized Wannier functions $W_{n}$.Wannier (1937); Marzari and Vanderbilt (1997); Marzari et al. (2012) This procedure significantly simplifies our analysis as Wannier functions are typically strongly localized on a single atomMarzari and Vanderbilt (1997); Brouder et al. (2007) and are therefore more convenient descriptors of the chemical environment of a material. Given a set of extended periodic orbitals $\Psi_{n\bm{k}}$ Wannier functions $W_{n}$ are defined as,

[TABLE]

where $N_{\bm{k}}$ is a total number of $k$-points. While unitary matrices $U_{nm\bm{k}}$ are in principle arbitrary, in the case of maximally localized Wannier functions they are chosen so that the resulting function $W_{n}(\bm{r})$ is as localized in real space as possible.Marzari and Vanderbilt (1997) We stress that basis change in Eq. (1) is exact in the sense that the vector space spanned by $W_{n}$ (and its periodic images) is the same as the space spanned by $\Psi_{m\bm{k}}$. However, $W_{n}$ is a more convenient object to study than $\Psi_{n\bm{k}}$ since a single function $W_{n}$ contains the same information as an entire band of Bloch functions $\Psi_{n\bm{k}}$ (there is one function for each $\bm{k}$-point). Therefore any information contained in bands is also contained in the corresponding Wannier functions. One of the Wannier functions for one of the parent structures (perovskites) is denoted in Fig. 2 (a) as $W_{\parallel}$. This Wannier function corresponds to the oxygen $p$-like function oriented along the B–O–B line (here B is the anion in the ABO3 perovskite that is in the center of the oxygen octahedron). The remaining two functions ($W_{\perp}$) centered on the same oxygen atom are perpendicular to the B–O–B line [see also Fig. 2 (a)]. In the remaining three families (rutile, pyrochlore, cristobalite) we also find that Wannier function on each oxygen atom follow the same logic: there is a singlet orbital (denoted again as $W_{\parallel}$) and a degenerate pair of orbitals ($W_{\perp}$). These Wannier functions are shown in Figs. 2 (b), (c), and (d).

Fourth, we construct the representation of the Hamiltonian operator $H$ in the Wannier basis from the previous step,

[TABLE]

Here $W_{m\bm{R}}$ is defined as $W_{m}$ translated by a lattice vector $\bm{R}$, $W_{m\bm{R}}(\bm{r})=W_{m}(\bm{r}-\bm{R}).$ It is straightforward to show that for structural properties we need not consider $H_{nm\bm{R}}$ when either $n\neq m$ or $\bm{R}\neq\bm{0}$. These terms are often referred to as hopping integrals. Hopping integrals do not contribute to the integrated band energy $\int_{\rm BZ}E_{n\bm{k}}d\bm{k}$ of a fully occupied band and therefore they don’t have the learning potential for determining ground state structure.222Total energy $E^{\rm tot}$ in the Kohn-Sham framework does not depend only on $\int_{\rm BZ}E_{n\bm{k}}d\bm{k}$ but also on additional terms that are explicit functions of electron density $n(\bm{r})$. However, density can be computed from the Wannier function as $n(\bm{r})=\sum_{n\bm{R}}|W_{n}(\bm{r})|^{2}$. Therefore, in some cases it is possible that one would have to include a descriptor for the shape of the Wannier function in order to learn about the material structure. That is not the case here as we consider materials in their parent structure, so oxygen Wannier function have nearly constant shape regardless of anion A or B. Therefore, we are left with $H_{nm\bm{R}}$ when both $n=m$ and $\bm{R}=\bm{0}$. The matrix element $H_{nn\bm{0}}$ is usually referred to as the onsite energy of the $n$-th Wannier function. The absolute value of onsite energy is ill-defined for a periodic solid as one can change its value by adding an arbitrary constant $C$ to the Hamiltonian of the periodic solid: $H\rightarrow H+C$. However, changing $H$ to $H+C$ changes onsite energies of all Wannier functions by the same amount. Therefore, even though absolute values are ill-defined, the differences of onsite energies,

[TABLE]

are well defined.

Fifth, we use symmetry to find all distinct values of $H_{nn\bm{0}}-H_{mm\bm{0}}$ from the previous step. Since the symmetry of all four parent structures is high, we find that in all four families there remains only one distinct numerical value in which $n$ corresponds to $W_{\perp}$ and $m$ to $W_{\parallel}$. This difference we denote as $X$,

[TABLE]

As an example, Fig. 2 (a) shows a complex of oxygen $p$-like electron bands (black) in SrTiO3 along with calculated $\matrixelement{W_{\perp}}{H}{W_{\perp}}$ and $\matrixelement{W_{\parallel}}{H}{W_{\parallel}}$ (dotted/dashed red line). The black arrow indicates their difference ($X$). The numerical value of $X$ for all materials we studied is provided in Tab. 1.

We stress here that the energy scale $X$ can’t be inferred from the electron band-energies $E_{n\bm{k}}$ alone. Instead, one also needs to use wave functions $\Psi_{n\bm{k}}$ in the construction. In addition, $X$ is gauge-dependent as a different choice of relative phases of wave functions will lead to a different Wannier function and thus different $X$. For example, one can show that a simple unitary rotation that rotates $W_{\parallel}$ and $W_{\perp}$ will produce two Wannier function with the difference in onsite energy having any value between $+X$ and $-X$, including zero. However, this procedure will increase the spread of the Wannier functions and will reduce their symmetry. Therefore, in this work we constructed Wannier functions for each family structure in a consistent way that respects the symmetry of the parent phase.

III Results

Table 1 contains calculated values of $X$ for the 64 compounds we studied. As can be seen from the table, $X$ is well correlated with the preferred ground state structure. We discuss now all four structure families in more detail.

Perovskites are one of the most studied complex oxides. The numerical values of $X$ reported in Tab. 1 for perovskites was calculated in the parent cubic phase with space group Pm$\bar{3}$m. As can be seen from the table, $X$ is positive for all perovskites and it ranges from $0.47$ eV to $1.66$ eV. There is a large number of structures (hettotypes) that derive from the cubic perovskite parent structure. Some of these structures are polar (denoted as p in Tab.1 and Fig. 3) while others are non-polar and have rotated oxygen octahedra (np).Glazer (1972); King-Smith and Vanderbilt (1994); Woodward (1997); Benedek and Fennie (2013) Our analysis shows that perovskites with $X$ less than $\sim 1$ eV tend to condense into a non-polar state (np) while those with $X$ larger than $\sim 1$ eV condense into a polar state (p). The exception is CdTiO3 which is polar but has $X=0.90$ eV and LaAlO3 which is non-polar but has $X=1.54$ eV. While the value of $X$ for CdTiO3 is near the polar–non-polar boundary, LaAlO3 is not. The anomalous value of $X$ for LaAlO3 likely occurs because its perovskite phase is degenerate with another structure.Persson

We note here that the importance of local interactions in perovskites has been discussed earlier in Refs. Cohen (1992); Bersuker (1995); Woodward (1997); Palmer (1997); Koleżyński and Tkacz-Śmiech (2005). However, we are unaware of any other work in which ground state density functional calculation of a high-symmetry perovskite alone can be used to infer its low-symmetry ground state structure.

Rutiles are another common structure of complex oxides.Baur (2007) Unlike perovskites, rutiles have a somewhat more complicated structure as their octahedral units are both corner and edge shared. The numerical values of $X$ reported in Tab. 1 for rutiles was calculated in the parent rutile structure (r) with tetragonal space group P42/mnm. As can be seen from the table compounds with $X$ larger than $0.17$ eV remain in the rutile phase, while those with smaller $X$ distort into a manganite (m) or NbO2-type (n) structure.Baur (2007) Both of these structures contain off-centered A anions and deformed oxygen octahedra. The more common, manganite structure, is monoclinic with space group P21/c while the structure unique to NbO2 is tetragonal I41/a. In rutiles calculated value of $X$ ranges between $-0.27$ eV and $1.06$ eV.

Pyrochlores have the most complex structure among the materials we studied. Variable $X$ was computed in the cubic pyrochlore (py) state with space group Fd$\bar{3}$m. We find that compounds with $X$ above $-0.2$ eV remain in the pyrochlore state at their ground state while others form into thortveitite (th)Brisse et al. (1972) polymorphfam with monoclinic space group C2/m. The only outlier we found is La2Sn2O7 as it has $X=-0.27$ eV but, as far as we are aware, it remains in a pyrochlore structure at a ground state.

Cristobalite is a well known structure in which anions are surrounded with oxygen tetrahedra, unlike the other three families in which oxygen atoms form octahedral units. The numerical values of $X$ reported in Tab. 1 for cristobalites were calculated in the parent cubic idealized cristobalite structure (space group F$\bar{4}$3m). We again find a strong correlation between $X$ and the preferred ground state structure. Compounds with $X$ larger than $\sim 4$ eV remain cristobalite (with small deviations known as $\alpha$ or $\beta$-cristobalite, space groups P41212 and I$\bar{4}$2d). Those with $X$ less than $\sim 4$ eV are unstable in the cristobalite structure and distort instead into anatase phase (a) with tetragonal space group I41/amd.fam In the cristobalites we studied the calculated values of $X$ range from $3.62$ to $4.67$ eV.

III.1 Comparison with Pettifor maps

While descriptor $X$ in this work was extracted from a first-principles calculation of a periodic solid, there are several descriptors in use that can be extracted just based on the chemical composition of the solid. One such example is Pettifor mapsPettifor (1986) which assign a phenomenological chemical scale $\chi$ to each chemical elementPettifor (1984) in the periodic table. Figure 4 shows such maps for all of the materials studied in this work. Horizontal and vertical axes in these plots correspond to the ordering of $\chi$ for metal atoms appearing in each compound (this ordering is also called Mendeleev number $m$).

As can be seen from the Fig. 4, Pettifor maps for these materials correlate well with the ground state structure only in the case of cristobalites. Namely, compounds with smaller $m$ ($\chi$) tend to deform into anatase phase (blue cross symbols) while those with larger $m$ tend to remain in the cristobalite phase (red square symbols). In the other three classes of materials correlations with Pettifor maps are not as good. For example, for rutiles we find that those with very large $m$ tend to remain in a rutile phase, but those with lower $m$ can be either rutile or manganite. In fact, the one with the lowest $m$ (TiO2) is in the rutile phase. Similarly, no obvious correlation is found for perovskites and pyrochlores.

IV Conclusion

Our work shows that it is possible to extract information from a first-principles calculation on a parent phase and correlate it with the preference for a ground state structure. While few of these correlations can be found using empirical parameters such as Pettifor maps, or ionic radii, our work shows that one can also extract relevant parameters from the fundamental variable in the problem ($\Psi_{n{\bm{k}}}$).

We expect that one could use this approach not only to characterize Hamiltonian operator $H$ but also any other quantum-mechanical operator, such as electron-light interaction,Ong et al. (2017) electron-phonon interaction,Giustino et al. (2007) or other interactions.