Parametrically Constrained Geometry Relaxations for High-Throughput Materials Science

TL;DR

This paper introduces a flexible parametric relaxation method that maintains symmetry constraints during electronic-structure calculations, significantly accelerating high-throughput materials discovery even with broken symmetries.

Contribution

The authors develop a generalizable parametric relaxation scheme that preserves symmetry constraints at any level, improving relaxation efficiency for distorted structures.

Findings

Successfully relaxes metastable structures

Reduces relaxation steps by an order of magnitude

Enhances high-throughput materials screening

Abstract

Reducing parameter spaces via exploiting symmetries has greatly accelerated and increased the quality of electronic-structure calculations. Unfortunately, many of the traditional methods fail when the global crystal symmetry is broken, even when the distortion is only a slight perturbation (e.g. Jahn-Teller like distortions). Here we introduce a flexible and generalizable parametric relaxation scheme, and implement it in the all-electron code FHI-aims. This approach utilizes parametric constraints to maintain symmetry at any level. After demonstrating the method's ability to relax metastable structures, we highlight its adaptability and performance over a test set of 359 materials, across thirteen lattice prototypes. Finally we show how these constraints can reduce the number of steps needed to relax local lattice distortions by an order of magnitude. The flexibility of these constraints enables a significant acceleration of the high-throughput searches for novel materials for numerous applications.

Tables

Table 1. Table 1: Parametric expressions for each component of the lattice vectors in the twelve atom cubic and tetragonal ZrO 2 supercell

	Cubic			Tetragonal
	$L{1}^x$	$L^y$	$L^z$	$L^x$	$L^y$	$L^z$
$L_1$	a	0	0	a	0	0
$L_2$	0	a	0	0	a	0
$L_3$	0	0	a	0	0	c

Table 2. Table 2: Parametric expressions for each component of the atomic positions in the twelve atom cubic and tetragonal ZrO 2 supercell

	Cubic			Tetragonal
Atom	$L_1$	$L_2$	$L_3$	$L_1$	$L_2$	$L_3$
Zr	0.00	0.00	0.00	0.00	0.00	0.00
Zr	0.50	0.50	0.00	0.50	0.50	0.00
Zr	0.00	0.50	0.50	0.00	0.50	0.50
Zr	0.50	0.00	0.50	0.50	0.00	0.50
O	0.25	0.25	0.25	0.25	0.25	0.25 - z2
O	0.25	0.75	0.25	0.25	0.75	0.25 - z2
O	0.75	0.75	0.75	0.75	0.75	0.75 + z2
O	0.25	0.25	0.75	0.25	0.25	0.75 + z2
O	0.75	0.25	0.25	0.75	0.25	0.25 - z2
O	0.25	0.75	0.75	0.25	0.75	0.75 - z2
O	0.75	0.75	0.75	0.75	0.75	0.75 + z2
O	0.25	0.25	0.25	0.25	0.25	0.25 + z2

Table 3. Table 3: Summary of the materials used in the test dataset

AFLOW Prototype	Space Group	# of Materials	Atoms per Unit cell	Free Parameters	Full d.o.f. / # Free Parameters
AB_oP8_62_c_c	62	8	8	7	4.71
A2B_oP12_62_2c_c	62	35	12	9	5.00
A2BC4_tI14_82_bc_a_g	82	35	7	5	6.00
A2BC4D_tI16_121_d_a_i_b	121	29	8	4	8.25
AB2_hP3_164_a_d	164	25	3	3	6.00
AB_hP4_186_b_b	186	37	4	2	5.25
AB_cF8_216_c_a	216	37	2	1	15.00
ABC_cF12_216_b_c_a	216	54	3	1	18.00
AB2_cF12_225_a_c	225	13	3	1	18.00
AB2C_cF16_225_a_c_b	225	14	4	1	21.00
AB_cF8_225_a_b	225	19	2	1	15.00
A_cF8_227_a	227	3	2	1	15.00
A2BC4_cF56_227_d_a_e	227	50	14	2	25.50

Table 4. Table 4: Summary of the free and constrained relaxation performance by AFLOW prototype.

			PBE			PBEsol
AFLOW Prototype	Space Group	# of Materials	Average Savings	# Preserved Space Group Free	# XTAL Match	Average Savings	# Preserved Space Group Free	# XTAL Match
AB_oP8_62_c_c	62	8	10.23	3	8	24.61	4	8
A2B_oP12_62_2c_c	62	35	10.84	19	29	18.32	20	33
A2BC4_tI14_82_bc_a_g	82	35	40.32	29	32	59.98	28	34
A2BC4D_tI16_121_d_a_i_b	121	29	44.01	23	26	58.13	21	26
AB2_hP3_164_a_d	164	25	7.03	10	24	19.68	5	24
AB_hP4_186_b_b	186	37	30.34	23	36	41.68	19	36
AB_cF8_216_c_a	216	37	29.71	28	36	54.83	31	36
ABC_cF12_216_b_c_a	216	54	33.35	44	54	74.03	36	46
AB2_cF12_225_a_c	225	13	57.63	8	9	54.20	7	9
AB2C_cF16_225_a_c_b	225	14	19.77	12	14	80.02	7	12
AB_cF8_225_a_b	225	19	32.50	11	19	51.72	7	19
A_cF8_227_a	227	3	35.12	3	3	47.62	3	3
A2BC4_cF56_227_d_a_e	227	50	67.10	37	47	73.58	37	47
Full Dataset		359	34.68	69.64%	93.87%	53.80	62.67%	92.76%

Equations

$$\mathbf{F}=-\frac{dE}{d\mathbf{R}},$$

$$\sigma =\frac{1}{V}\frac{dE}{d\mathbf{L}},$$

$$\frac{\mathrm{dE}}{\mathrm{d}\mathcal{L}}={{\mathcal{L}}^T}^{-1}V\cdot \sigma$$

$${\mathcal{F}}_{\mathcal{L}}=-\frac{\mathrm{dE}}{\mathrm{d}\mathcal{L}}-{{\mathcal{R}}_{\mathcal{F}}}^T{\mathcal{F}}_{\mathcal{R}}.$$

$$\mathbf{r}$$

$$\mathbf{l}$$

$${\mathbf{F}}_{\mathbf{r}}$$

$${\mathbf{F}}_{\mathbf{l}}$$

$$\bm{\mathcal{J}_{R}}=\left(\begin{array}[]{cccc}{\bm{\mathcal{L}}^{T}}&{0}&{\dots}&{0}\\ {0}&{\bm{\mathcal{L}}^{T}}&{\dots}&{0}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {0}&{0}&{\dots}&{\bm{\mathcal{L}}^{T}}\end{array}\right)\bm{\mathcal{J}_{Rf}}.$$

$${\mathcal{A}}^{-1,L}={\left({\mathcal{A}}^T\mathcal{A}\right)}^{-1}{\mathcal{A}}^T,$$

$${\mathcal{H}}_{\mathbf{r}}$$

$${\mathcal{H}}_{\mathbf{l}}$$

$$\bm{\mathcal{H}}=\left(\begin{array}[]{cc}{\bm{\mathcal{H}_{r}}}&{0}\\ {0}&{\bm{\mathcal{H}_{l}}}\end{array}\right).$$

$$m=1-(1-\text{dev})(1-\text{disp})(1-\text{fail}),$$

$$m\le 0.1:$$

$$0.1<m\le 0.2:$$

$$m>0.2:$$

$$S=\frac{N_{\text{free}}-N_{\text{constrained}}}{N_{\text{constrained}}}\times 100\%$$

Full text

I Introduction

Symmetry preservation and breaking is one of the most fundamental processes in physics and chemistry. Many properties and applications such as, piezoelectricity Panda2009 ; Wang2012 ; Ok2016 ; Nguyen2013a , pyroelectricity Moure2015a ; Isaenko2015a , ferroelectricity Damjanovic1998 ; Hang2011a ; Sun2014c ; Shi2016a , topological insulators Cava2013 ; Fiete2012a , and non-linear optics Evans2002a ; Wu2018a ; Xu2019a , require certain selection rules to be met, and therefore require certain crystallographic symmetries to be maintained. Furthermore, it is not only global symmetry, the space and point groups of a material, but also local symmetry breaking that matters. For example, defects can cause significant changes in a material’s mechanical Le2014 ; Sayle and optical Pan ; Wang2016d properties as well as in its electronic An2014 ; Anno1999 ; Senanayak2017 and thermal transport Proshchenko2019 ; Toberer2010 ; Xie2013 coefficients. These effects can be particularly important at thin-film interfaces where interactions between different layers can induce systematic distortions in the structure of the film. In turn, these distortions can lead to novel properties in the materials, such as artificial ferroelectricity in layered perovskite supercell structures Wang2016 ; Bousquet2008 . Clearly, such problems cannot be addressed by enforcing a global symmetry constraint, e.g., space group conservation, for the whole system, but require to selectively preserve and break material’s symmetries locally, e.g., around defects or in the individual perovskite layers. This is paramount in computational material science, especially in high-throughput studies, which often aim at calculating and exploring yet unknown properties of already known materials, e.g., band structures Ricci2017 , defect formation energies Broberg2018 , elastic and thermal properties Plata2017a and topological constants Tang2019 ; Maghirang2019 . Similarly, many high-throughput studies aim at discovering potentially stable or meta-stable materials by decorating complex, well-known crystal structures such as Hauslers Roy and perovskites Mazaheri2019 with different species, or by systematically exploring a given alloy system Sutton2018 . To streamline such calculations it is essential to keep both global and local symmetries under control, especially when complex materials or material properties are targeted. In this work, we achieve this goal by proposing and implementing parametric geometric constraints that allow for enforcing or breaking symmetries both globally and locally. Before applying these constraints, an understanding of how the constraints map onto the target systems (e.g., the number of atoms in the unit cell, space group, or the effects of local distortions on the structure) is necessary. To facilitate setting up such constraints, we rely on the AFLOW Library of Crystallographic Prototypes Mehl2017 ; Hicks2019 to generate the initial mapping of real space onto a reduced parameter space that fully describes a system. One can then manually alter the initial mapping to add or lift constraints as needed. This allows for the efficient targeting of specific geometric configurations and avoids revisiting and recalculating already investigated configurations.

Traditionally, crystallographic symmetries are incorporated in first-principles codes already at the electronic-structure level (e.g., by sampling $\bm{k}$-space grids in the irreducible part of the Brillouin zone, as implemented in Vasp Hafner2008 or by sampling real space in symmetry defined “irreducible wedges”, as done in parsec Kronik2006 ) since it leads to significant savings in memory and computational workload for highly symmetric crystals. Also, by this means the obtained forces on the atoms and stresses on the lattice vectors fully reflect the crystallographic symmetries. Since geometry relaxation algorithms such as steepest descent, conjugate gradient, Newton-Raphson, quasi-Newton (e.g. BFGS Fletcher1987 ), and truncated-Newton methods NoceWrig06 rely on the forces and stresses to update the atomic and lattice degrees of freedom, global symmetries are inherently preserved in such approaches. However, this does not allow for the partial, local symmetry breaking discussed in the introduction. To address such cases in first-principles calculations, it is typically necessary to lift all crystallographic symmetry constraints and treat the atomic and lattice degrees of freedom as a set of freely changing parameters. Besides the increased computational cost, such unconstrained structure optimization can lead to long and inefficient relaxation trajectories, resulting in structures far from the ideal geometry. While in some cases this problem can be circumvented by fixing atomic, lattice, or internal Freysoldt2017 ; Panosetti2015 degrees of freedom (as done in Quantum Espresso Giannozzi2009 or Vasp Hafner2008 ), mapping local distortions onto these constraints requires cumbersome manual inspection and analysis, if it is even possible. A more direct approach targeting how the distortion changes the native crystal structure provides an easier and better way of treating these systems.

Here we present a scheme to incorporate parametric constraints in structure optimizations that treats all levels of symmetry equally. The proposed approach employs a mapping of the relevant degrees of freedom onto a lower-dimensional representation of the structure; the respective forces and stresses are then automatically mapped in this reduced representation. With that, the implemented formalism does not require to alter the employed relaxation algorithm, while still allowing the introduction of arbitrary constraints in a user-friendly manner. We first describe how the methodology works and the tools that can be used to quickly generate new structures. We demonstrate that these constraints allow for performing geometry optimizations on dynamically stabilized structures, which are not easily addressable otherwise. By analyzing the constrained and unconstrained relaxations of a test set of 359 materials, we then show that these constraints are also computationally beneficial for the relaxation of stable materials. Finally we illustrate how the parameters can be used to selectively break symmetries and accelerate relaxations in supercells.

II Transformation to reduced space

For a free relaxation, the optimizer acts on the full $3N+9$ dimensional potential-energy surface $E(\bm{R},\bm{L})$ of a material, which is encoded by the atomic, $\bm{R}$, and lattice, $\bm{L}$, degrees of freedom. The lattice degrees of freedom are stored as the three components of the three lattice vectors in the chosen unit cell, and the atomic degrees of freedom are the $3N$ components of the positions of the $N$ atoms in a unit cell, represented by Cartesian or fractional coordinates. The forces, $\bm{F}$, acting on the atomic degrees of freedom are the derivatives of the energy with respect to $\bm{R}$

[TABLE]

while the forces acting on the lattice vectors stem from the stresses, $\sigma$,

[TABLE]

where $V$ is the volume of the unit cell. In ab initio approaches, $E$ is determined by solving the electronic-structure problem, and the respective derivatives are obtained analytically via the Hellmann-Feynman Theorem. However, in practice this requires to account for additional terms, such as the Pulay terms and multipole corrections, as done in FHI-aims Blum2009 ; Knuth2015 .

Because the underlying potential-energy surfaces (PES) are complex, relaxing certain polymorphs of a material on these surfaces can be challenging or even impossible. As an example, zirconia (ZrO2) can exist in its pure form in three different crystal phases: a high temperature ($T>2370^{\,\circ}$C) cubic phase, an intermediate temperature ($1170^{\,\circ}\text{C}\leq T\leq 2370^{\,\circ}$C) tetragonal phase, and a low temperature ($T<1170^{\,\circ}$C) monoclinic phase Bocanegra-Bernal2002 . In particular, the cubic phase is a dynamically stabilized phase representing an average structure that is rarely, if ever actualized in pristine ZrO2 Carbogno2014 . Accordingly, this phase constitutes a saddle point of the PES and its phonon band structure exhibits an imaginary mode at the X point Parlinski1997 , as illustrated in Figure 1a. The associated eigenvector, illustrated in the respective inset, describes a pairwise, antiparallel distortion of the oxygen atoms that goes hand-in-hand with a stretching of the lattice and leads to the tetragonal structure Bocanegra-Bernal2002 . To help illustrate this, in Figure 1b we plot a two dimensional PES for a twelve atom zirconia unit cell over a reduced parameter set describing the cubic lattice parameter $a$ and the motion along the imaginary mode, $z_{2}$. The PES has two wells corresponding to the equivalent tetragonal structures, and a saddle point between them representing a high-symmetry configuration, i.e., the high-temperature cubic phase. On this PES, a free relaxation of the cubic phase would result in the material relaxing towards a local minima; however, by constraining the relaxation to act only on $a$, the cubic phase can be obtained as shown in the outsets in Figure 1b.

To help demonstrate how parametric constraints can facilitate addressing these stable/unstable polymorphs, we provide the respective constraints in Tables 1 and 2. As discussed later, these symmetry blocks correspond to the input files that are required in the proposed formalism for the twelve atom zirconia unit cell. For the cubic polymorph, all atoms are at fixed fractional coordinates and only one parameter, i.e., the lattice constant $a$ can change in a parametrically constrained relaxation. For the relaxation algorithm, this means that the optimizer can now act only on $a$, all other degrees of freedom remain untouched. Clearly, this ensures that the initial space group is retained during relaxation. As mentioned above, such constraints are necessary, since this polymorph corresponds to a saddle point of the PES, so that an unconstrained or free relaxation would effectively allow the system to break the symmetry and to descend in a local minimum. To explicitly explore these local minima, the constraints imposed on the cubic cell can be stepwise lifted, whereby the information contained in the imaginary phonon eigenvectors can be incorporated as parametric constraints. As shown in Tables 1 and 2, the pairwise distortion of the oxygen atoms for the imaginary mode at X, cf. Figure 1a, can be described by introducing one additional parameter for the oxygen distortion, $z_{2}$, and one for the tetragonality of the lattice $c$. These constraints ensure that the geometry optimization occurs along the imaginary phonon mode and leads to the tetragonal minimum. Conversely, a free relaxations can again lead to other local minima of the PES. In this textbook example, the same constraints could have been imposed by relaxing cubic and tetragonal zirconia in their primitive cells with 6 and 3 atoms, respectively. Generally, this is however not the case: some materials, e.g., the bismuth oxide Matsumoto2010 ; Drache2007 discussed later, have high-temperature polymorphs with the same number of sites as their stable structures, or more.

The previous input example illustrates the flexibility of these constraints, but knowledge of which reduced parameters to use and their relation to the full geometry, must be known before generating an input file. For crystals, these are determined by the space group and the Wyckoff positions and can therefore be manually constructed. Luckily, these parameters are already a part of the definitions used in the AFLOW Library of Crystallographic Prototypes, allowing for an easy way to define these constraints for numerous materials via their utilities Mehl2017 ; Hicks2019 . The library sorts materials by their space group, stoichiometry, and occupied Wyckoff sites, as calculated with AFLOW-SYM Hicks2018b , placing all materials that share those features into the same crystal prototype Mehl2017 ; Hicks2019 . A reduced parameter space can then be generated from a prototype definition, and used to describe that class of materials. For example the tetragonal phase of zirconia can be described by only three parameters: length of the lattice vectors in the $\bm{a}$ and $\bm{b}$ directions ($a$), the ratio of the lattice vectors ($\frac{c}{a}$), and the magnitude of the oxygen distortions ($z_{2}$). These parameters represent the same ones we defined earlier from analysis of the phonon bandstructure, with an additional parameter allowing for the relaxation of the lattice upon the atomic distortions.

For all prototypes defined in the library, the automatized generation of input geometries for VASP Hafner2008 , FHI-aims Blum2009 , Quantum Espresso Giannozzi2009 , Abinit Gonze2002 and more codes is supported by AFLOW. The AFLOW library contains 590 unique structure prototypes across all 230 space groups and is thus extremely suitable as a starting point for high-throughput studies. As of version 3.1.204, the option --add_equations can be added to the AFLOW command to generate FHI-aims geometry.in files already containing the additional block required for the constrained relaxation. Because of this, we use the crystal prototypes defined by the AFLOW Library of Crystallographic Prototypes throughout this work. Due to the analytic representation of the parametric expressions it is also straight-forward to add additional parameters to allow for lower symmetric structures or distortions as well as removing parameters to constrain specific components even further. Additionally, because the AFLOW prototypes are partially based on the space group and occupied Wyckoff sites, it is also straight-forward to adapt their techniques to include structure classes not currently in the library.

In practice, the parametric constraints are implemented in the following fashion. Let us assume a $(3\times 3)$-dimensional lattice vector matrix, $\bm{\mathcal{L}}$, and a $(N\times 3)$-dimensional matrix, $\bm{\mathcal{R_{F}}}$, for the fractional atomic positions. Given the atomic forces, $\bm{\mathcal{F_{R}}}$, on the Cartesian atomic positions and the stress tensor $\sigma$, we can calculate the derivatives of the energy with respect to the lattice components doll

[TABLE]

where $V$ is the unit cell volume and obtain the generalized forces on the lattice, $\bm{\mathcal{F_{L}}}$, after cleaning from the atomic contributions

[TABLE]

Each of these matrices denoted by calligraphic letters, $\bm{\mathcal{R_{F}}}$, $\bm{\mathcal{F_{R}}}$, $\bm{\mathcal{L}}$ and $\bm{\mathcal{F_{L}}}$, can be flattened to one-dimensional vectors that we will name $\bm{R_{F}}$, $\bm{F_{R}}$, $\bm{L}$, and $\bm{F_{L}}$ respectively. In the parameter representation these quantities reduce to their small-letter counterparts, the $M_{R}$-dimensional $\bm{r}$, $\bm{F_{r}}$ and $M_{L}$-dimensional $\bm{l}$ and $\bm{F_{l}}$, via

[TABLE]

where $M_{R}$ and $M_{L}$ are the number of free parameters in the atomic and lattice degrees of freedom; $\bm{\mathcal{J}_{R}}$, $\bm{\mathcal{J}_{Rf}}$ and $\bm{\mathcal{J}_{L}}$ are the Jacobian matrix for the transformations and $\bm{t_{Rf}}$ and $\bm{t_{L}}$ are the translation vectors for the respective fractional atomic and lattice degrees of freedom. Here $\bm{\mathcal{J}_{R}}$ represents the transformation of the atomic coordinates from Cartesian space to the reduced space, which is calculated from $\bm{\mathcal{J}_{Rf}}$ by

[TABLE]

The translation vectors are used to include any constant shifts, which are not captured by the Jacobians. Because $\bm{\mathcal{J}_{R}}$, $\bm{\mathcal{J}_{Rf}}$ and $\bm{\mathcal{J_{L}}}$ are not square and therefore not regularly invertible, we use the generalized left inverse rao1972 defined for a matrix $\bm{\mathcal{A}}$ as

[TABLE]

provided $\bm{\mathcal{A}}$ has full column rank. The transformation back to real space can then be performed by inverting Equations 5a-5d. The back-transformation of the forces to the full space is not necessary but can be helpful to obtain symmetrized Cartesian or Fractional forces, to check for the convergence of the relaxation.

To facilitate the construction of the Jacobian matrices, we assume a linear relationship between the full coordinates and the parameters. In principle,$\bm{\mathcal{J}_{Rf}}$ and $\bm{\mathcal{J}_{L}}$ can be constructed at each step by using analytical expressions to describe each real space degree of freedom as a function of the reduce parameter set; however, by assuming a linear relationship between the spaces they can be initialized at the start of the calculation and used at every step. For the atomic positions, this assumption is already fulfilled by using fractional instead of Cartesian coordinates. If we allow angles as unit cell parameters, which is the case for the monoclinic and triclinic lattice systems, the relations become non-linear containing for example expressions like $c\cdot\cos{(\beta)}$. In these cases, the easiest solution is to substitute each non-zero lattice vector component with an independent parameter.

Before the relaxation the Hessian, $\bm{\mathcal{H}}$, is initialized in the full coordinate space, split into atomic and lattice blocks ($\bm{\mathcal{H}_{R}}$ and $\bm{\mathcal{H}_{L}}$ respectively), and individually transformed into reduced coordinate space, $\bm{\mathcal{H}_{r}}$ and $\bm{\mathcal{H}_{l}}$ via separate Jacobians, $\bm{\mathcal{J}_{R}}$ and $\bm{\mathcal{J}_{L}}$

[TABLE]

Here $\bm{\mathcal{J}_{R}}$ is also divided by the average unit vector length, $V^{1/3}$, so $\bm{\mathcal{H}_{r}}$ and $\bm{\mathcal{H}_{l}}$ are on a similar scale. The total Hessian is then recombined resulting in

[TABLE]

Figure 2 illustrates the procedure for relaxing structures with these constraints. During the relaxation, a full SCF cycle is completed to obtain the forces and the stress tensor for the current geometry, at each step. If the convergence criterion is fulfilled, i.e. if the forces are below a given threshold, then the relaxation stops and returns the current geometry. Otherwise the lattice vectors as well as the atomic positions and their respective forces are mapped onto the reduced space using the transformation described in Equations 5a-5d. The atomic coordinates and forces are respectively scaled by $V^{1/3}$ and $V^{-1/3}$, and then passed on to the optimizer. In FHI-aims this is usually a BFGS/TRM optimizer. Once the optimized parameters are obtained, the full geometry is reconstructed from the parameters and a new relaxation step can begin.

III Applications and Benchmarks

III.1 Relaxing Metastable and Unstable Systems

In some cases, constraining a relaxation is necessary to keep the structure in its given polymorph. Similar to what was seen for zirconia, a material can have many phases that are metastable or unstable at zero point conditions, but are stabilized by entropic contributions at higher temperatures or pressures. Here we define a metastable phase to be one that is in a local minimum on its PES. While freely relaxing stable or metastable structures is possible by respectively using an initial geometry near its corresponding global or local minimum on the PES, unstable systems will tend to relax towards lower energy and usually lower-symmetric structures, unless they are somehow constrained.

To demonstrate the ability of these constraints to optimize such structures, we relax the twelve atom cubic zirconia unit cell from Figure 1. While most relaxations will be performed on the primitive cells of structures, we use this system as a simple, demonstrative example. The initial structure for the cubic phase was taken from the AFLOW Library of Crystallographic Prototypes, while the initial structure for the tetragonal phase is the same structure with a minor perturbation along the imaginary mode seen in Figure 1.

All calculations are done using the FHI-aims package, a full-potential, all-electron electronic structure code. FHI-aims utilizes numeric atom-centered orbital basis functions, grouped into different tiers beyond the minimal set needed to describe free atoms. For these calculations we use tier 1 (double numeric plus polarization basis set) with light basis settings which were shown to calculate the lattice parameter and cohesive energy of face-centered cubic gold within 0.001 Å  and 20 meV Blum2009 , respectively. We use the PBEsol as the exchange-correlation functional; SCF convergence criteria of $10^{-6}$ eV/Å and $5\times 10^{-4}$ eV/Å for the density and forces, respectively; and the structures are relaxed until the maximum forces on the degrees of freedom are below 0.005 eV/Å. All other inputs were taken to be the default values in FHI-aims. While a larger basis set and using a hybrid functional would increase the accuracy of the calculations, we do not expect it to affect the performance of the relaxation scheme.

Figure 3a shows that using the constraints both the cubic and tetragonal phase of ZrO2 can be converged in 4 and 10 steps, respectively, while only the tetragonal phase can be obtained in 30 steps with a free relaxation. The free relaxation of the cubic phase proceeds towards the tetragonal phase, but initially stalls at a non-physical simple cubic phase in 37 steps. If the relaxation convergence criteria is further reduced to 0.001 eV/Å the structure reaches the tetragonal phase in 114 steps.

Another example of a material with many metastable phases is bismuth oxide. Bismuth oxide exists in several different polymorphs Matsumoto2010 including the low temperature monoclinic phase, $\alpha$-; the high-temperature, face-centered cubic phase, $\delta$-; the metastable, body-centered cubic phase, $\gamma$-; and the metastable, tetragonal phase $\beta$-Bi2O3 Drache2007 . Upon heating $\alpha$-Bi2O3 transforms into $\delta$-Bi2O3 at around 730∘C, and remains stable until it melts at approximately 825∘C Drache2007 . Depending on the cooling procedure, $\delta$-Bi2O3 transitions to one of the two metastable phases, $\beta$- or $\gamma$-Bi2O3, at approximately 650 ∘C or approximately 640 ∘C, respectively Drache2007 . Upon further cooling $\beta$-Bi2O3 and $\gamma$-Bi2O3 respectively return to the $\alpha$-phase at $\sim$300∘C and at a temperature dependent on the cooling rate Drache2007 . Importantly, unlike ZrO2, both the tetragonal and body-centered cubic phases have as many or more atoms in their primitive cells than the $\alpha$ phase, so freely relaxing both structures in their primitive cells should not prevent them from exploring lower symmetry structures. For these calculations we use the same computational settings as those used for ZrO2, but with the intermediate settings for the basis set. The intermediate settings and basis sets in FHI-aims increases the accuracy of the default numerical settings, but more importantly adds an $f$-orbital to the double numeric plus polarization basis set for oxygen, which is necessary for describing monoclinic structures containing oxygen.

Figure 3b shows the relaxation for the $\alpha$-, $\beta$-, and $\gamma$- bismuth oxide phases. Multiple experimental crystal structure refinements exist Drache2007 for both the $\beta$- and $\gamma$- phases, so we take the relaxed structures from the Materials Project, which were initialized from structures taken from the ICSD. In all cases the constrained relaxation takes less time to converge with the $\gamma$-, $\beta$-, and $\alpha$-phase respectively needing 42, 17, and 20 steps to converge. Freely relaxing both the $\alpha$ and $\beta$-phases simply increases the number of steps needed to converge the system to 25 and 66 steps, respectively, but removing the constraints for $\gamma$-Bi2O3 causes it to relax to an unknown, non-symmetric structure in 160 steps. While it is possible that the divergent relaxation is a result of an incorrect structure refinement, if that is the case, then it suggests that comparing the final energy and forces of the constrained and free relaxation trajectories can be used to indicate if a structure is correct. Since the free relaxation of the $\gamma$-phase departs from the known phases of the material, a more in-depth study of it would be impossible without the constraints as the fully relaxed structures no longer represent the same material.

Both of these cases demonstrate the need for constraining relaxations to their crystal prototype for high-throughput applications. For the high-symmetry phases of both zirconia and bismuth oxide, the free relaxation not only initially converged to a different phase, but also unknown and potentially physically unrealizable ones. While the relaxation of ZrO2 does eventually reach one of its known phases, bismuth oxide remains in incorrect structure. Any further calculations on those structures would be erroneous and could lead to both false positives and skewed property descriptors. While integrity checks could be made for some of the materials, the consistent breaking of the system’s symmetry and the inconsistent degree of that symmetry breaking makes developing standardized checks impractical. This will be particularly useful for testing the predictions from crystal structure discovery methods where exact knowledge of lattice type and decorations is necessary. Combining these constraints with an automatically generated parametric representation, would provide an efficient means to optimize the newly predicted structures.

III.2 Bench-marking the Algorithm

As the above examples show, the implemented constraints not only ensure that symmetry is preserved, but also accelerate the relaxation because the optimization of a reduced representation with less free parameters is by definition a less demanding task. To quantify this, the new relaxation algorithm is tested on a set of 359 materials across multiple lattice systems and AFLOW prototypes, as summarized in Table 3 and listed in the Supplementary Information. The chosen prototypes represent a sample of common materials with varied space groups and parametric relation complexity. To further understand how the method performs within a class, we try to include only prototypes with a significant amount of available structures.

The initial geometry of each material is taken from either the AFLOW Curtarolo2012 or Materials Project Jain2013c database and converted into the right format using the Atomic Simulation Environment HjorthLarsen2017 . All relaxations are done in FHI-aims using the same settings as the zirconia calculations, using both the PBE and PBEsol functionals.

Table 4 and Figure 4 illustrate the largest benefit of using the constraints: the preservation of the initial crystallographic prototypes. Approximately 7% of all materials tested relax to a different structure according to the AFLOW-XTAL-MATCH tool Hicks2019b . This tool measures the similarity between two materials using techniques similar to those of Burzlaff, et al.Burzlaff1997 and produces a misfit value, $m$,

[TABLE]

where dev, disp, and fail are normalized representations of deviations in the lattice vectors, atomic positions, and a failure indicator for when an atomic deviation is more than half the shortest distance in the coordination polyhedra. From $m$ we can determine if a structure is considered a match using the following mapping

[TABLE]

All of the 26 materials with divergent relaxations in the data set follow the same pattern as cubic-ZrO2 and $\gamma$-Bi2O3. Therefore constraining these material’s relaxation is vital because without them all further calculations on these materials would no longer be physically relevant.

Constraining the relaxation has benefits even when the final structures are similar as it significantly reduces the number of steps needed for the trajectories to converge. When the constrained and free relaxations proceed towards the same structures, the constraints reduce the number of relaxation steps by an average of 33.11% and 52.43% for calculations using the PBE and PBEsol functional, respectively. Including the divergent structures respectively increases the saving to 34.68% and 53.80%, but this is expected as the constrained relaxation and free relaxations are acting on qualitatively different PESs. Here we define the savings, $S$, as

[TABLE]

where $N_{\text{free}}$ is the number of steps need to converge the free relaxation and $N_{\text{constrained}}$ is the number of steps need to converge the constrained relaxation. The seemingly better performance of the constraints when using the PBEsol, is likely because of larger differences between the initial and final structures when using this functional. On average comparing the starting geometries with the freely relaxed ones gave an average $m$ of 0.07654 (0.01630 when $m\neq 1$) and 0.1322 (0.02038 when $m\neq 1$) for the PBE and PBEsol calculations, respectively. And since increasing the distance to the final structure also increases the likelihood of the free relaxation deviating from the constrained trajectory, larger savings using the PBEsol functional is expected. This is also supported by the larger differences between the final structures of constrained and free relaxations for the PBEsol functional.

Unfortunately, the savings are not consistent across the various prototypes. Figure 5 shows the total number of steps needed to relax the structures with and without constraints for the PBEsol calculations, sorted by the space group and then the maximum number of relaxation steps needed. Similar to the average savings shown in Table 4, as the number of free parameters approaches the number of degrees of freedom in a material the savings from constraining the relaxation decrease. In some cases using the constraints actually increase the number of steps needed for the relaxation to converge, but in all of these cases for PBEsol the relaxation trajectories take unproductive steps shown in the inset of Figure 5 for PtS2. In these cases the extra degrees of freedom allow the relaxation trajectories to pass the problematic regions faster and therefore converge to the final structures in fewer steps. For the PBE calculations, there are also some cases where the relaxation takes a few extra steps at the end of the trajectory where the total forces are near convergence criteria. These results suggest that for lower symmetry structures, the potential savings from constraining the relaxation can be considerably smaller.

Constraining the relaxation can also decrease the computational time needed to perform further calculations after the relaxations. Despite their similarity to the final geometries of the constrained relaxation, freely relaxed structures always break symmetry to some degree. To measure how well the relaxation preserves symmetry we compare the spglib calculated space groups of the initial and converged structures Togo2018 . spglib calculates the space group of a material by iteratively searching for a given structure’s primitive cell and symmetry operations and using those to generate the space group. The algorithm determines if a structure’s space group has a certain symmetry operation, by checking if the operator transforms all the atoms in the structure to sites occupied by the same type of atom within a given small euclidean distance, $\varepsilon$. By default $\varepsilon$ is set to 10*-5* Å, and none of the tested materials remain in their initial space groups using this setting. To preserve the symmetry for the freely relaxed structures, a larger user-defined tolerance threshold is needed (0.01 Å is used in Table II), consistent with findings reported by Hicks, et al. Hicks2018b . However, when using the constraints, the symmetry is always perfectly preserved for all materials. This is advantageous as exploiting symmetry can greatly reduce the time needed to do further calculations for many applications. In particular phonon calculations using the finite difference approach where the symmetry of the system determines the number of atomic displacements, and therefore the number of force calculations, needed to calculate a complete set of force constants. By default phonopy, a very popular Python package for calculating phonon spectra with finite differences phonopy , uses the default spglib settings to calculate the symmetry of a material, which would misrepresent all of the tested materials using free relaxation. Commensurate phonon calculations can be achieved by performing symmetry constrained relaxations and using more robust symmetry tools, such as AFLOW-SYM Hicks2018b .

III.3 Systems with local symmetries or distortions

One of the key advantages of defining constraints in this way is the ability to locally break symmetry to account for point defects. While the previously discussed benefits of lower relaxation times and relaxing unstable structures could also be achieved by using symmetrized forces, that method would not be able to locally break its symmetry. However many of these defects exhibit some type of short range order, such as Jahn-Teller-type lattice distortions Rodriguez2018 ; Prentice2017 ; Evarestov2012 , that can be parametrically added on top of the standard crystal structure. By including these distortions one can reduce the computational cost of relaxing supercells with defects, which is important when using a supercell approach to study point defects. This approach uses density functional methods to calculate the total energy of a system in a series of supercells of increasing size to study the effects of a defect on the system. If a scaling law is known for the system the results of these calculation can be extrapolated to the experimentally relevant dilute limit. As the supercell size increases the time needed to calculate each step also increases, therefore any reduction in the number of steps needed can save a significant amount of computational time.

To illustrate the ability of the new relaxation scheme to study point defects in materials, we study a polaronic distortion in MgO previously studied in our department at the Fritz Haber Institute Kokott2018 . Polarons are quasiparticles that couple point charges with lattice distortions in a material, reducing its charge carrier mobility.

Typically, the size of the polaron is controlled by the electron-phonon interaction, with a stronger interaction leading to a smaller polaron Kokott2018 ; Nag2014 . Because of the reduced charge carrier mobility, understanding how polarons form and migrate through a material is vital for a number of applications ranging from catalysis Dohnalek2010 ; Henderson2011 ; Rettie2016 to thermoelectricity Nag2014 . For rock salt MgO the lattice distortions for an electron hole polaron are Jahn-Teller like, with the oxygen and magnesium ions being respectively attracted and repelled from the hole. To model this system, we assume the electron hole is located on a fixed oxygen atom placed in the center of the supercell. We then allow the other ions to relax from their initial position along a line going through the center of the cell, with the magnitude and sign of the motion being controlled by a single free parameter for each atom. This constraint choice reduces the number of degrees of freedom to $N-1$, where $N$ is number of atoms in the supercell. While these constraints do not account for periodic images of the electron hole or changes in the positions of the other ions, it does preserve the main polaronic distortion. Moreover, the restrictiveness of the constraints can be decreased further on subsequent calculations until the desired level of accuracy is achieved.

Figure 6 illustrates the uncorrected polaron energy for the free and constrained relaxation trajectories for these supercells Kokott2018 . Because charge localization is necessary when studying polarons, the HSE 2006 functional is used with a screening parameter of 0.11 Bohr*-1* with exact exchange. For the unrestricted spin relaxation the SCF density, forces, total energy, and eigenvalues were all converged to 10*-4* eV/Å, 10*-4* eV/Å, 10*-5* eV, and 10*-2* eV, respectively. Each atom calculated 5 empty states above those used by the Kohn-Sham orbitals and the structures were relaxed until the total forces on the free parameters were below 10*-4* eV/Å. The constrained relaxation for the distorted geometry converges the 64 atom supercell in 11 steps which is only one-eighth of the steps needed by the free relaxation. Significantly better performance is seen for the 216 atom supercell, which converges in 10 constrained relaxation steps, 96% less than in the unconstrained case needing 234 steps. This increase in efficiency is likely a result of the constraints focusing on the primary Coloumbic interaction, reducing the amount of fine-tuning it needs to do to balance out the weaker perturbation in the other inter-ionic interactions.

Using the constraints does cause the relaxation to converge to a structure with slightly higher energy, but this is because the chosen constraints are valid for an isolated polaron and thus do not account for the artificial interactions with periodic images due to the finite supercell size. In practice, this effect is rather negligible, since we find that the obtained polaron energy is 78.4 meV higher in energy than the free relaxation in the 2x2x2 supercell and 69.1 meV higher in the 3x3x3 supercell.

IV Conclusions and Outlook

In this work we presented a new scheme for parametrically relaxing structures in a general, symmetry-reduced space. After explaining the algorithm, we test it on 359 different materials across a broad range of material classes. In all cases the new method was able to strictly preserve the symmetry of the materials, and on average reduced the number of steps needed to converge a material by 50%. We also demonstrated for the example of bismuth oxide how the constraints can be used to relax to metastable phases. Finally, we showcased the relaxation of structures with local symmetry breaking with known distortion patterns for polarons in MgO.

This new method will have a profound impact on computational materials discovery. Not only does the decreased cost of relaxing a material increase the velocity of high-throughput search, but it also allows for those searches to explore metastable and dynamically stabilized structures. The method also has the promise to improve the efficiency of supercell calculations and study only the physically relevant structures. Finally by monitoring the difference between the full forces and symmetrized forces new stable phases can potentially be discovered from metastable or unstable polymorphs. Although we showed that the proposed algorithm is applicable to accelerate and improve standard solid-state physics calculations, its flexibility allows it to be applied to a much wider range of problems, e.g., transition-state searches or interface relaxations. Similarly, it is easily generalizable to other form of coordinates and straightforwardly implementable in any electronic-structure theory code, and has already been implemented as constraints within ASE.

IV.1 Data Availability

The full input and output files for the calculations in this work are available in the NOMAD repository with the identifier [DOI:10.17172/NOMAD/2019.10.19-1] Lenz2019 .

IV.2 Code Availability

The algorithms presented in this paper have been fully implemented in FHI-aims and is available in any version after 20191019. The constraints have also been implemented in ASEhttps://gitlab.com/ase/ase/merge_requests/1298.

IV.3 Author Contributions

MOL and TARP contributed equally on this work by implementing and bench-marking the discussed relaxation scheme. DH worked to incorporate the constraints into AFLOW. SC, MS, and CC directed the project. All authors analysed the data and wrote the manuscript.

IV.4 Acknowledgements

TARP would like to acknowledge Florian Knoop for valuable discussions. The authors also would like to Acknowledge reviewer 3 for pointing out different refinements of the metastable Bi2O3 phase. This project was supported by TEC1p (the European Research Council (ERC) Horizon 2020 research and innovation programme, grant agreement No 740233), BigMax (the Max Planck Society’s Research Network on Big-Data-Driven Materials-Science), and the NOMAD pillar of the FAIR-DI e.V. association. SC and DH acknowledges U.S. DOD-ONR (Grants No. N00014-17-1-2090). D.H. acknowledges support from the U.S. DOD through the National Defense Science and Engineering Graduate (NDSEG) Fellowship Program. We thank the Max Planck Computing and Data Facility for computational resources.