Introduction

Symmetry preservation and breaking is one of the most fundamental processes in physics and chemistry. Many properties and applications, such as piezoelectricity,1,2,3,4 pyroelectricity,5,6 ferroelectricity,7,8,9,10 topological insulators,11,12 and nonlinear optics,13,14,15 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 mechanical16,17 and optical18,19 properties, as well as in its electronic20,21,22 and thermal transport23,24,25 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.26,27 Clearly, such problems cannot be addressed by enforcing a global symmetry constraint, e.g., space group conservation, for the whole system, but require selectively preserving and breaking a 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,28 defect-formation energies,29 elastic and thermal properties,30 and topological constants.31,32 Similarly, many high-throughput studies aim at discovering potentially stable or metastable materials by decorating complex, well-known crystal structures such as Heuslers33 and perovskites34 with different species, or by systematically exploring a given alloy system.35 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 the enforcing or breaking of symmetries, both globally and locally. Before applying these constraints, an understanding of how they 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 Prototypes36,37 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 \({\boldsymbol{k}}\)-space grids in the irreducible part of the Brillouin zone, as implemented in Vasp38 or by sampling real space in symmetry defined “irreducible wedges”, as done in PARSEC39), 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., BFGS40), and truncated-Newton methods41 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 above. 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 internal42,43 degrees of freedom (as done in Quantum Espresso44 or VASP38), 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 altering 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.

Transformation to reduced space

For a free relaxation, the optimizer acts on the full \(3N+9\) dimensional potential-energy surface \(E({\boldsymbol{R}},{\boldsymbol{L}})\) of a material, which is encoded by the atomic, \({\boldsymbol{R}}\), and lattice, \({\boldsymbol{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, \({\boldsymbol{F}}\), acting on the atomic degrees of freedom are the derivatives of the energy with respect to \({\boldsymbol{R}}\)

$${\boldsymbol{F}}=-\frac{{\mathrm{d}}E}{{\mathrm{d}}{\boldsymbol{R}}},$$
(1)

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

$$\sigma =\frac{1}{V}\frac{{\mathrm{d}}E}{{\mathrm{d}}{\boldsymbol{L}}},$$
(2)

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.45,46

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 (ZrO\({}_{2}\)) can exist in its pure form in three different crystal phases: a high-temperature (\(T\,>\,237\,{0} ^{\circ }\)C) cubic phase, an intermediate temperature (\(117\,{0} \, ^{\circ }{\text{C}}\,\le T\le 237\,{0} \, ^{\circ }\)C) tetragonal phase, and a low-temperature (\(T\,<\,117\,{0} \, ^{\circ }\)C) monoclinic phase.47 In particular, the cubic phase is a dynamically stabilized phase representing an average structure that is rarely, if ever actualized in pristine ZrO\({}_{2}\).48 Accordingly, this phase constitutes a saddle point of the PES, and its phonon band structure exhibits an imaginary mode at the X point,49 as illustrated in Fig. 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.47 To help illustrate this, in Fig. 1b we plot a two-dimensional PES for a 12-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 toward 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 Fig. 1b.

Fig. 1: Reduced Parameter Space for 12-Atom Zirconia Supercell.
figure 1

An illustration of the new relaxation scheme. a The phonon band structures of cubic ZrO\({}_{2}\). The inset illustrates the eigenvector of the imaginary mode that drives the system toward the tetragonal structure in a six-atom supercell. b A two-dimensional potential-energy surface for ZrO2. The minimum energy structure is set to 0.0 eV, and the contour lines correspond to a 0.2 eV increase in energy. The maroon dot represents a structure corresponding to a high-temperature cubic phase of the material. The outsets show the one-dimensional potential-energy surface along each mode.

Full size image

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 12-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\), and all other degrees of freedom remain untouched. Clearly, this ensures that the initial space group is retained during the relaxation. As mentioned above, such constraints are necessary, since this polymorph corresponds to a saddle point of the PES. This means that an unconstrained or free relaxation would effectively allow the system to break its symmetry and to descend into a local minimum. To explicitly explore these local minima, the constraints imposed on the cubic cell can be lifted in a stepwise manner, 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 (Fig. 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 relaxation 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 six and three atoms, respectively. However, this is not always the case: some materials, e.g., bismuth oxide,50,51 discussed later, have high-temperature polymorphs with the same number of sites as their stable structures, or more.

Table 1 Example of lattice vector parametric constraints.
Full size table
Table 2 Example of atomic parametric constraints.
Full size table

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.36,37 The library sorts materials by their space group, stoichiometry, and occupied Wyckoff sites, as calculated with AFLOW-SYM,52 placing all materials that share those features into the same crystal prototype.36,37 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: the length of the lattice vectors in the \({\boldsymbol{a}}\) and \({\boldsymbol{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 the analysis of the phonon band structure, 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,38 FHI-aims,45 Quantum Espresso,44 Abinit,53 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 straightforward to add additional parameters to allow for lower-symmetric structures or distortions as well as removing parameters to constrain specific components even further. In addition, because the AFLOW prototypes are partially based on the space group and occupied Wyckoff sites, it is also straightforward to adapt their techniques to include structure classes not currently in the library.

Results

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 of the PES. While freely relaxing stable or metastable structures is possible by 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 12-atom cubic zirconia unit cell from Fig. 1. While most relaxations will be performed on the primitive cells of structures, we use this system as a simple, demonstrative example. The initial geometry for the cubic phase was taken from the AFLOW Library of Crystallographic Prototypes, while the one for the tetragonal phase is the same as the cubic phase with a minor perturbation along the imaginary mode seen in Fig. 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’ (as defined in Blum et al.45) 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,45 respectively. We use PBEsol as the exchange-correlation functional; SCF convergence criteria of \(1{0}^{-6}\) eV/Å and \(5\times 1{0}^{-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 2a shows that using the constraints both the cubic and tetragonal phase of ZrO\({}_{2}\) can be converged in four and ten 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 toward the tetragonal phase, but initially stalls at a non-physical simple cubic structure in 37 steps. If the relaxation convergence criteria is further reduced to 0.001 eV/Å, the structure reaches the tetragonal phase in 114 steps.

Fig. 2: Using parameteric constraints to converge to unstable and metastable phases.
figure 2

Convergence behavior of the free (squares) and constrained (circles) relaxations for a the tetragonal phase (green) and cubic (purple) phases of ZrO\({}_{2}\) and b the \(\alpha\)-phase (blue), \(\beta\)-phase (gray), and \(\gamma\)-phase (red) of Bi\({}_{2}\)O\({}_{3}\). For both the cubic phase of ZrO\({}_{2}\) and the \(\gamma\)-phase of Bi\({}_{2}\)O\({}_{3}\), the free relaxation breaks the symmetry and finds an energetically lower structure which is the tetragonal phase for ZrO\({}_{2}\) and an unknown phase for Bi\({}_{2}\)O\({}_{3}\). The energy scale on the y-axis is set to 1 meV below the minimum energy for each material.

Full size image

Another example of a material with many metastable phases is bismuth oxide. Bismuth oxide exists in several different polymorphs50 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\)-Bi\({}_{2}\)O\({}_{3}\).51 Upon heating \(\alpha\)-Bi\({}_{2}\)O\({}_{3}\) transforms into \(\delta\)-Bi\({}_{2}\)O\({}_{3}\) at around 730 °C, and remains stable until it melts at ~825 °C.51 Depending on the cooling procedure, \(\delta\)-Bi\({}_{2}\)O\({}_{3}\) transitions to one of the two metastable phases, \(\beta\)- or \(\gamma\)-Bi\({}_{2}\)O\({}_{3}\), at ~650 °C or ~640 °C, respectively.51 Upon further cooling, \(\beta\)-Bi\({}_{2}\)O\({}_{3}\) and \(\gamma\)-Bi\({}_{2}\)O\({}_{3}\), respectively, return to the \(\alpha\)-phase at ~300 °C and at a temperature dependent on the cooling rate.51 Importantly, unlike ZrO\({}_{2}\), 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 ZrO\({}_{2}\), but with the ‘intermediate’ settings for the basis set. The ‘intermediate’ settings and basis sets in FHI-aims increase the accuracy of the default numerical settings, but more importantly add an \(f\)-orbital to the ‘tier 1’ basis set for oxygen, which is necessary for describing monoclinic structures containing oxygen.

Figure 2b shows the relaxation of the \(\alpha\)-, \(\beta\)-, and \(\gamma\)-bismuth oxide phases. Multiple experimental crystal structure refinements exist51 for both the \(\beta\)- and \(\gamma\)-phases, so we take the relaxed structures from the Materials Project,54 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\)-Bi\({}_{2}\)O\({}_{3}\) 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 geometry 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 to unknown and potentially physically unrealizable ones. While the relaxation of ZrO\({}_{2}\) does eventually reach one of its known phases, bismuth oxide remains in an incorrect structure. Any further calculations on those structures would be erroneous and could lead to both false positives and skewed property descriptors during both high-throughput and machine learning studies. 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.

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 explicitly listed in the Supplementary Information (Supplementary Tables 1, 2). 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 AFLOW55 or the Materials Project54 database and converted into the right format using the Atomic Simulation Environment (ASE).56 All relaxations are done in FHI-aims using the same settings as the zirconia calculations, using both the PBE and PBEsol functionals.

Table 3 Summary of materials test set.
Full size table

Table 4 and Fig. 3 illustrate the largest benefit of using the constraints: the preservation of the initial crystallographic prototypes. Approximately, 6% of all materials tested relax to a different structure according to the AFLOW-XTAL-MATCH tool.57 This tool measures the similarity between two materials using techniques similar to those of Burzlaff, et al.58 and produces a misfit value, \(m\),

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

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:

$$\begin{array}{ll}m\,\le\, 0.1:\ \ &{\rm{structures}}\ {\rm{are}}\ {\rm{similar}}\\ 0.1\,<\,m\le\, 0.2:\ \ &{\rm{structures}}\ {\rm{are}}\ {\rm{within}}\ {\rm{same}}\ {\rm{family}}\\ m\,>\,0.2:\ \ &{\rm{structures}}\ {\rm{are}}\ {\rm{not}}\ {\rm{compatible}}.\end{array}$$

All of the 26 materials with divergent relaxations in the data set follow the same pattern as cubic ZrO\({}_{2}\) and \(\gamma\)-Bi\({}_{2}\)O\({}_{3}\). Therefore, constraining these material’s relaxation is vital because without them all further calculations on these materials would no longer be physically relevant.

Table 4 Summary of constraint performance.
Full size table
Fig. 3: Constrained and unconstrained relaxations go to the same structure.
figure 3

The misfit between the fully relaxed and constrained structures using AFLOW-XTAL-MATCH. All structures below the horizontal line are considered matching. The structures are ordered by space group and then by \(m\).

Full size image

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 toward 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 and free relaxations are acting on qualitatively different PESs. Here we define the savings, \(S\), as

$$S=\frac{{N}_{\text{free}}-{N}_{\text{constrained}}}{{N}_{\text{constrained}}}\times 100 \%$$
(4)

where \({N}_{\text{free}}\) is the number of steps needed to converge the free relaxation and \({N}_{\text{constrained}}\) is the number of steps needed to converge the constrained relaxation. The seemingly better performance of the constraints when using 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\,\ne\, 1\)) and 0.1322 (0.02038 when \(m\,\ne\, 1\)) for the PBE and PBEsol calculations, respectively. 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 are 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 4 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 increases the number of steps needed for the relaxation to converge, but in all of these cases, for PBEsol calculations, the relaxation trajectories take unproductive steps shown in the inset of Fig. 4 for PtS\({}_{2}\). 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 the convergence criteria. These results suggest that for lower symmetry structures, the potential savings from constraining the relaxation can be considerably smaller.

Fig. 4: Constraints accelerate relaxation trajectories.
figure 4

The total number of steps needed for the constrained (green) and free (purple) relaxations. Negative step numbers represent the cases where the constrained relaxations take longer than the free relaxation. The bar for orthorhombic N2O is not shown in the figure, since the free relaxation took 310 steps, and the constrained one required 296 steps to converge. The inset shows the constrained (green circles) and free relaxation (purple squares) trajectory of platinum (IV) sulfide, which is one of the materials where the constrained relaxation takes more steps than the free relaxation. The zero of the energy scale is set to 1 meV below the energy of the relaxed structure (horizontal dashed line).

Full size image

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.59 The package 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 is preserved under 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 4), consistent with findings reported by Hicks, et al.52 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, for phonon calculations using the finite difference approach, 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. The package phonopy, a very popular Python package for calculating phonon spectra with finite differences,60 uses the default spglib settings to calculate the symmetry of a material, which would misrepresent all of the tested materials obtained via a free relaxation. Commensurate phonon calculations can be achieved by performing symmetry constrained relaxations and using more robust symmetry tools, such as AFLOW-SYM.52

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 a structure’s symmetry. However, many of these defects exhibit some type of short range order, such as Jahn-Teller-type lattice distortions,61,62,63 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 calculations 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.64 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.64,65 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 catalysis66,67,68 to thermoelectricity.65 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 map all other interactions onto the main Coulombic distortion, we expect those interactions to have a minor effect on the final structure. Moreover, the restrictiveness of the constraints can be decreased further on subsequent calculations until the desired level of accuracy is achieved.

Figure 5 illustrates the uncorrected polaron energy for the free and constrained relaxation trajectories for these supercells.64 Because charge localization is necessary when studying polarons, the HSE 06 functional is used with a screening parameter of 0.11 Bohr−1 with a fraction of exact exchange, α, of one. 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. Five empty states per atom are used, 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 Coulombic interaction, reducing the amount of fine-tuning it needs to do to balance out the weaker perturbation in the other interionic interactions. Using the constraints does cause the relaxation to converge to a structure with slightly higher energy, but this is because the chosen constraints approximate 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 2 × 2 × 2 supercell and 69.1 meV higher in the 3 × 3 × 3 supercell.

Fig. 5: Constraints quickly find distorted geometries.
figure 5

The convergence of the free (squares) and constrained (circles) relaxation of a polaronic distortion in a 2 × 2 × 2 (purple) and 3 × 3 × 3 (green) supercell of rock-salt MgO. The zero of the energy scale is set to be 1 meV below the energy of the minimal energy structure. The horizontal dashed lines correspond to the final energy of constrained relaxation.

Full size image

Discussion

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 searches 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 to study only 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 forms of coordinates and straightforwardly implementable in any electronic-structure theory code, and has already been implemented as constraints with ASE by the authors.

Methods

Implementation details

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

$$\frac{{\rm{dE}}}{{\rm{{{d}}}}{{\boldsymbol{\mathcal{L}}}}}={{{\boldsymbol{\mathcal{L}}}}}^{{T}^{-1}}\ V\cdot \sigma,$$
(5)

where \(V\) is the unit-cell volume and obtain the generalized forces on the lattice, \({{{\boldsymbol{\mathcal{F}}}}}_{{{\boldsymbol{\mathcal{L}}}}}\), after cleaning from the atomic contributions

$${{{\boldsymbol{\mathcal{F}}}}}_{{{\boldsymbol{\mathcal{L}}}}}=-\frac{{\rm{dE}}}{{\rm{d}}{{\boldsymbol{\mathcal{L}}}}}-{{{\boldsymbol{\mathcal{R}}}}}_{{{\boldsymbol{\mathcal{F}}}}}^{T}\ {{{\boldsymbol{\mathcal{F}}}}}_{{{\boldsymbol{\mathcal{R}}}}}.$$
(6)

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

$${\boldsymbol{r}}={{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{Rf}}}^{-1}\ \left({{{{\boldsymbol{R}}}}}_{{{{\boldsymbol{F}}}}}-{{\boldsymbol{t}}}_{{\boldsymbol{Rf}}}\right),$$
(7a)
$${\boldsymbol{l}}={{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{L}}}^{-1}\ \left({\boldsymbol{L}}-{{\boldsymbol{t}}}_{{\boldsymbol{L}}}\right),$$
(7b)
$${{\boldsymbol{{F}}}}_{{\boldsymbol{{r}}}}={{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{R}}}^{T}\ {{\boldsymbol{{F}}}}_{{\boldsymbol{{R}}}},$$
(7c)
$${{\boldsymbol{F}}}_{{\boldsymbol{l}}}={{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{L}}}^{T}\ {{\boldsymbol{F}}}_{{\boldsymbol{L}}},$$
(7d)

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

$${{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{R}}}=\left(\begin{array}{llll}{{{\boldsymbol{\mathcal{L}}}}}^{T}&0&\ldots &0\\ 0&{{{\boldsymbol{\mathcal{L}}}}}^{T}&\ldots &0\\ \vdots &\vdots &\ddots &\vdots \\ 0&0&\ldots &{{{\boldsymbol{\mathcal{L}}}}}^{T}\end{array}\right){{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{Rf}}}.$$
(8)

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

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

provided \({{\boldsymbol{\mathcal{A}}}}\) has full column rank. The transformation back to real space can then be performed by inverting Eqs (7a)–(7d). 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, \({{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{Rf}}}\), and \({{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{L}}}\) can be constructed at each step by using analytical expressions to describe each real space degree of freedom as a function of the reduced 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 nonlinear containing for example expressions like \(c\cdot \cos (\beta )\). In these cases, the easiest solution is to substitute each nonzero lattice vector component with an independent parameter.

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

$${{{\boldsymbol{\mathcal{H}}}}}_{{\boldsymbol{r}}}={{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{R}}}^{T}{{{\boldsymbol{\mathcal{H}}}}}_{{\boldsymbol{R}}}{{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{R}}},$$
(10)
$${{{\boldsymbol{\mathcal{H}}}}}_{{\boldsymbol{l}}}={{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{L}}}^{T}{{{\boldsymbol{\mathcal{H}}}}}_{{\boldsymbol{L}}}{{{\boldsymbol{\mathcal{J}}}}}_{{\boldsymbol{L}}}.$$
(11)

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

$${{\boldsymbol{\mathcal{H}}}}=\left(\begin{array}{ll}{{{\boldsymbol{\mathcal{H}}}}}_{{\boldsymbol{r}}}&0\\ 0&{{{\boldsymbol{\mathcal{H}}}}}_{{\boldsymbol{l}}}\end{array}\right).$$
(12)

Figure 6 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 Eqs. (7a–d). 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/modified BFGS optimizer. Once the optimized parameters are obtained, the full geometry is reconstructed from the parameters and a new relaxation step can begin.

Fig. 6: The Parameteric constraints workflow.
figure 6

Workflow of the relaxation constrained to the “parameter reduced space”.

Full size image