Tight-binding approximation enhanced global optimization

Subscriber access provided by - Access paid by the | UCSB Libraries

Structure Prediction

Tight-binding approximation enhanced global optimization Maxime Van den Bossche, Henrik Grönbeck, and Bjørk Hammer J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00039 • Publication Date (Web): 28 Mar 2018 Downloaded from http://pubs.acs.org on March 28, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Tight-binding approximation enhanced global optimization Maxime Van den Bossche,∗,†,‡ Henrik Grönbeck,¶ and Bjørk Hammer§ †Department of Chemistry, Brown University, Providence, RI, United States ‡Science Institute and Faculty of Physical Sciences, University of Iceland, 107 Reykjav´ık, Iceland ¶Department of Physics and Competence Centre for Catalysis, Chalmers University of Technology, 412 58 Göteborg, Sweden §Department of Physics and Astronomy, Aarhus University, DK-8000 Aarhus, Denmark E-mail: maxime van den [email protected] Abstract Solving and predicting atomic structures from first-principles is limited by the computational cost of exploring the search space, even when relatively inexpensive density functionals are used. Here, we present an efficient approach where the search is performed using density functional tight-binding, with an automatic adaptive parametrization scheme for the repulsive pair potentials. We successfully apply the method to the genetic algorithm optimization of bulk carbon, titanium dioxide, palladium oxide and calcium hydroxide, and assess the stability of the unknown crystal structure of palladium hydroxide.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

Finding global energy minima is important in computational chemistry and material science, as thermodynamically stable structures are frequently the most abundant and therefore determine the chemical and physical properties of the material. Following the increase in computational power over the last decades, such global optimization (GO) investigations are increasingly being performed using first-principles electronic structure methods such as Kohn-Sham density functional theory (DFT). 1,2 Recent examples include metal clusters on oxide surfaces, 3 metal oxide clusters, 4 tin dioxide surface terminations, 5 superhard transition metal borides, 6 and stepped titania surfaces. 7 However, such studies have yet not become routine due to the large computational cost of performing on the order of hundred to several thousand geometry optimizations, with each typically requiring some 50 to 500 force calls. Even with comparatively inexpensive density functionals, the total cost grows rapidly to excessive levels with increasing system size. This problem is even more pronounced in the case of crystal structure prediction, where also the cell vectors should be optimized. Three main approaches have emerged to increase the overall computational efficiency. The first is the development of more efficient local optimization algorithms, among which the universal preconditioned optimizers are a prominent recent example. 8 The second approach is the formulation of more efficient GO algorithms, among which genetic algorithms (GAs) have arguably seen most widespread use and with major development done by e.g. the Oganov, 9–11 Zurek, 12,13 Johnston, 14,15 and Hammer 16–18 groups. Adaptive approaches are a third, less used strategy where a cheaper (semi-)empirical method is trained using first-principles data and refined during the GO. 19 Such approaches have also been considered for uses other than structure prediction, such as examining the robustness of interatomic potentials. 20 To our knowledge, its first example was the GAGA approach by Hartke and coworkers, which was applied to small silicon clusters. 21,22 The empirical method in Ref. 21 utilized modified Stillinger-Weber analytical potentials, which have often been used in modeling of group IV semiconductors. More recent examples are work by the Ho group, 23–25 where the embedded 2


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


atom method (EAM) is chosen as the empirical potential. Even though such potentials can be used for certain metal alloys, (per)oxides and silicates, the shortcomings are apparent. Firstly, the initially obtained EAM parametrizations are of low quality, so that many DFTbased reparametrizations are required before the GM structure is eventually found. Secondly, the approach is expected to have a limited applicability, considering the inappropriateness of EAM potentials for describing e.g. directional bonds and charge redistribution. It should be noted that also machine learning using neural networks has recently been proposed as a suitable empirical method. 26 Although neural networks are known to be adaptive, this approach is currently not able to provide a large advantage over regular DFT optimizations, due to the large computational cost of both training and running the neural network. In this work, we show that using density functional tight-binding (DFTB) 27 as the semiempirical method holds several key advantages compared to previously reported alternatives. Thanks to the ‘built-in’ chemical knowledge present in the electronic parameters (orbital eigenvalues, Hubbard values, and Slater-Koster integrals), only a limited amount of DFT input is needed to obtain a sufficiently accurate DFTB model for exploring the search space. Additionally, the parametrization procedure is straightforward in comparison to the fitting of interatomic potentials. Even though DFTB calculations are significantly slower than calculations using interatomic potentials, a speedup of circa two orders of magnitude is retained compared to typical DFT calculations. A DFTB-based strategy is, furthermore, expected to possess greater reliability compared to interatomic potentials. We have adopted the TANGO acronym (Tight-binding Approximation eNhanced Global Optimization) for this approach and its Python implementation. The remainder of the paper is structured as follows. First, we briefly describe the relevant aspects of DFTB (Section 2). The TANGO method is then described in detail in Section 3. Next, the general features of the GA searches are outlined (see Section 4). The associated DFT setup can be found in Section 5. The performance of the method is illustrated for a series of materials in Section 6. The studied examples with known crystal structures include

3



Page 4 of 36

bulk carbon, titanium dioxide, palladium oxide, and calcium hydroxide. We then turn to the unknown crystal structure of palladium hydroxide and discuss its stability with respect to hydrated palladium oxide. Benefits and potential improvements of the method are discussed in Section 7 before summarizing.

2

Density functional tight-binding

DFTB exists in different versions, and we have chosen to use self-consistent charge (SCC) DFTB, 28 also referred to as DFTB2. 29 By approximating the DFT energy as a second-order perturbation with respect to the reference charge density distribution, the SCC-DFTB total energy is obtained as a sum of the band structure energy, a charge fluctuation term, and repulsive interactions: Etot = Eband + Efluc + Erep .

(1)

Eband is represented as the sum of the eigenvalues of the lowest occupied orbitals,

Eband =

X

σ σ wk fn,k n,k ,

(2)

n,k,σ

with wk being the k-point weight, f the orbital occupation and n, k and σ the band, k-point and spin indices, respectively. The eigenvalues are obtained by solving the secular equation for each spin component: [H(k) − k S(k)]C(k) = 0

(3)

with H the Hamiltonian, S the overlap matrix and C the expansion coefficients. The (molecular) orbitals ψ are expanded in a minimal basis set consisting of the (confined) valence orbitals φ of the isolated atoms: ψn =

X

Cn,m φm

(4)

m

In a one-center approximation, the diagonal elements of the H are taken equal to the (Kohn4


Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Sham) DFT eigenvalues of the corresponding orbitals φ of the isolated atoms. The offdiagonal Hamiltonian elements are evaluated using a two-center approximation, and are as a result related to the corresponding orbitals and Kohn-Sham potentials of the isolated atoms and the interatomic distance. The overlap integrals and Hamiltonian matrix elements can, therefore, be computed in advance. Using the Slater-Koster transformation, 30 the number of orbital orientations which need to be evaluated is drastically reduced, leading to the so-called Slater-Koster tables. The density fluctuation term Efluc accounts for the effect of charge transfer as a result of differences in electronegativities. This includes both the individual charge transfer energies as well as the resulting Coulombic interactions. By representing the density change using atom-centered Gaussian distributions, this energy can be written as

Efluc =

1X γij (rij )∆qi ∆qj 2 i,j

(5)

where the sums run over the atomic indices. The γij function depends on the interatomic distance rij and the Hubbard parameters U of the elements involved. The U -values are in general taken equal to those of the isolated atoms, i.e. as the difference between the ionization potential and electron affinity. It is furthermore possible to use different U parameters for the different orbital angular momenta `. The charge fluctuations ∆q are typically calculated using a Mulliken population analysis. 31 This entails that Equation 3 needs to be solved iteratively since the density fluctuation term in H should be consistent with the resulting occupied electronic states. The repulsive term Erep includes core-core repulsion and exchange-correlation effects, and is generally treated empirically to improve the agreement with the parent density functional result. A simple, yet effective approach is to express Erep as a sum of short-ranged pairwise repulsive interactions: Erep =

1X Vrep (rij ). 2 i6=j 5


(6)


Initial DFTB parametrization Atomic parameters εℓ, Uℓ

SlaterKoster tables

Dimer curves E(rij)

Initial Vrep(rij) fitting

Circa 20 short DFT relaxations

Relax random structures Single-point DFT calculations on selected structures

Global optimization (GO) phase

Perform DFTB-GO searches

Final refinement

Select best unique structures

Refine Vrep(rij)

DFT relaxation

Figure 1: Overview of the work flow in the TANGO approach. The repulsive potentials Vrep are set up for each element pair (A,B) and are required to AB smoothly decay to zero at the pair-specific cutoff radii rcut . These cutoff radii are regularly

chosen in between the typical nearest-neighbor and next-nearest-neighbor distances. What are now the variables that can be adjusted to make the DFTB parametrization reproduce the DFT results? By design, the main degrees of freedom are the form and range of the repulsive potentials. Additionally, it is possible to tune electronic aspects of the Hamiltonian, namely the orbital eigenvalues , the Hubbard values U and the orbital confinement procedure 32–36 Here we have chosen to follow standard procedures for the electronic parameters and to only optimize the repulsive potentials. This is both out of convenience and out of lack of necessity, as previous work has shown that the repulsive potentials are most influential in describing geometries and energy differences. 35,37 The `-dependent eigenvalues and Hub6


Page 6 of 36

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


bard parameters are calculated here using the atomic DFT code available through the CP2K suite 38–40 with relativistic effects included through the Douglas-Kroll-Hess method. 41–43 The Slater-Koster tables are generated using the Hotbit code 44 with a parabolic confinement potential for the atomic orbitals. Standard confinement radii are used which amount to 1.85 times the covalent radius of the element in question. 27,44 All DFTB calculations are carried out with DFTB+. 45

3

The TANGO approach

We start with a couple of general considerations. The overall objective of the method is to accelerate the search towards low-lying minima of the DFT potential energy surface (PES), most notably the global minimum (GM). To this end, the method aims to generate a suitable DFTB parametrization with which the PES can be searched, resulting in valuable suggestions for the low-lying minima of the DFT PES. The generated DFTB parametrization should therefore mimick the chosen density functional in terms of geometries and total energy differences with sufficient accuracy. The method should use as few DFT force calls as possible and be largely automated. These considerations have led us to the scheme outlined in Figure 1, which broadly consists of two main stages followed by a final refinement step. In the first stage, the electronic parts of the DFTB Hamiltonian are constructed from the properties of the isolated atoms, as described in Section 2. This includes the the SlaterKoster tables, the (`-dependent) eigenvalues and Hubbard values for chosen electronic configurations of the different elements. An initial set of repulsive potentials is then constructed on the basis of random structures generated using the same random structure generator as employed by the GO algorithm (see Section 4). We have found it sufficient to generate only about 20 of such random structures, which are subjected to very short DFT variable-cell relaxations (circa 25 force calls each). Only the energies and forces of the initial and partially relaxed structures are included in the training data, as these samples are least correlated.

7



Page 8 of 36

The fits of Vrep are thereafter obtained following the procedure described below (paragraph 3.1). The obtained potentials are quickly tested by relaxing a number (e.g. 20) of new random structures with the resulting DFTB parametrization and performing single-point DFT calculations on the relaxed structures, followed by refitting the repulsive potentials. For best results, these random tests should be run more than once (e.g. five times). The main benefit from performing these tests is to ensure that the repulsive potentials are sufficiently accurate at short bond distances, which may be underrepresented in the initial training set. The accuracy with which the DFTB model describes the low-lying DFT minima will, up to a certain limit, improve upon reparametrization. However, this improvement represent a further refinement rather a than qualitative change. In the next stage, several GO runs are performed with the resulting DFTB parametrization. From the generated list of structures, the most stable, unique structures are selected. The search can be aborted at this step, after which the selected structures can, for example, be reoptimized at the DFT level and subjected to further analysis. Alternatively, only single-point DFT calculations are conducted and the results are added to the training set, followed by refitting the repulsive potentials and resuming or restarting the GO searches with the updated DFTB parametrization. The best unique DFTB structures are collected and subjected to further analysis, generally including a full DFT reoptimization.

3.1

Repulsive potential fitting

We have chosen to represent each repulsive potential Vrep (r) by the following piece-wise function:

    e−a1 r + a2 + a3     Pp i Vrep (r) = i=2 ci (rcut − r)       0

8

0 ≤ r < rexp rexp ≤ r < rcut rcut ≤ r


(7)

Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


At the shortest interatomic distances (r < rexp ), an exponentially decaying function is used, containing three free parameters (a1 , a2 , a3 ). The separation point rexp is determined as the gas phase dimer bond length with a potential energy 5 eV above the energy of the isolated, non spin-polarized atoms. The ai parameters are fitted to reproduce the DFT cohesive energies of the gas phase dimers (at separations below rexp ) using the Levenberg-Marquardt algorithm with Boltzmann-distributed weights (kB T = 50 eV). This choice performs satisfactorily as interatomic distances below rexp are not present in relaxed structures. It furthermore avoids the need to sample such high-energy configurations of the studied material with DFT. At the more relevant interatomic distances (rexp ≤ r < rcut ), we employ a polynomial function, which is a common choice in DFTB. The reason is that the optimization of the pair potentials with this choice becomes an ordinary linear regression problem in the polynomial coefficients ci . The cutoff radii rcut are simply taken equal to 1.5 times the sum of covalent radii of the elements being considered. The degree p of the polynomial is typically chosen between 3 and 9. The polynomial coefficients are optimized through a least-squares fit following a similar approach as previous studies on automatic DFTB parametrization strategies. 46,47 For a training set consisting of M structures with N atoms each, the total residual to be minimized is given by

RSSQ =

M X

ωi

i=1

+ωf

DFTB, no rep. DFT Ecoh, −C i − Ecoh, i

N X j=1

2 #

DFTB, no rep. 2 DFT ||F¯i,j − F¯i,j ||

(8)

i2 X h pair pair pair pair +ωd Vrep, poly (rexp ) − Vexp (rexp ) . pairs

Ecoh stands for the cohesive energy1 and F are the forces; we have not found it necessary to also include the stress tensor. The energy and force discrepancies of each structure with 1

For simplicity Ecoh is calculated with respect to the non spin-polarized atoms.

9



Figure 2: Flowchart of the GA strategy employed here as the GO algorithm. respect to DFT are weighted by a factor ωi which we take to be Boltzmann distributed, i.e. DFT ωi = exp −(EiDFT − Emin )/kB T , with kB T = 10 eV. The fitting parameter C furthermore shifts all DFTB cohesive energies by a constant. This is commonly employed to correct systematic over- or underbinding tendencies. 36,46,48,49 An additional weighting factor ωf is ˚2 in order to adhere most weight to required for the forces, which is taken equal to 10−3 A the description of the energetics. The last term in Eq. 8 penalizes discontinuities between the polynomial and exponential functions at the separation point rexp . A weight ωd of 106 ensures the match is sufficiently close in all cases.

4

Genetic algorithm

Our GO method employs GA searches using the framework 51 in the Atomic Simulation Environment (ASE). 52,53 In connection to this work, we have implemented in the ASE-

10


Page 10 of 36

Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: The different genetic operators used in the GA. Operator

Description

Cut-and-splice pairing

Combines the (scaled) atomic positions from two parent structures lying on opposite sides of a cutting plane with normal along a randomly chosen unit cell direction. 9 The new unit cell is created as a linear combination of the parent cells with random weighting. Applies a (symmetric) strain matrix to the atomic positions and the unit cell of the parent structure. 9 The randomly generated strains are normally distributed with standard deviation of 0.7.

Strain mutation Soft mutation

Displaces the atoms along the vibrational mode with lowest nonzero frequency. 50 The normal mode analysis is carried out based on a pairwise harmonic model. The force constants k are calculated as a function of the distance r as k = k0 (r0 /r)2 , with k0 and r0 the force constant and bond distance of the gas-phase dimer (calculated with PBE).

Rattle mutation

Introduces random displacements of the circa 80% of the atoms with amplitudes up to 2.5 ˚ A.

Rattlerotation mutation

In addition to rattle mutation, 30% of the hydroxide moieties (if present) are rotated over angles between π/2 and 3π/2 along randomly chosen axes.

GA code several genetic operators for bulk crystal optimization which have been previously developed by the Oganov group (see description below). For the studied metal hydroxides, a constrained GA approach is used 54 where the hydroxide bond distances are preserved throughout the entire GA search (i.e. including local optimizations and genetic operations). The characteristics of the applied GA strategy are shown in Figure 2. First, an initial set of structures is generated at random within appropriate limits for the cell shape and interatomic distances. Cell vector lengths and cell angles are allowed to vary between 2 and 8 ˚ A and π/5 and 4π/5 radians, respectively. We furthermore employ a standard cell splitting strategy 50 with splitting factors of 2 (50 %), 2-2 (25%) and 1 (no splitting, 25%). In addition, the cells are scaled to yield densities of circa 0.1 atoms/˚ A3 . All interatomic distances are restricted to be larger than 0.6 times the sum of covalent radii. The generated structures are in the next step subjected to variable-cell local minimization using the limited-memory BFGS 55 and FIRE 56 algorithms with Armijo linesearches 57 and a 11



universal, exponential-type preconditioner. 8 In our experience, the preconditioning reduces the required number of force calls by a factor of 2 to 5, and we therefore strongly recommend its use. The actual GA loop starts by ensuring that the population contains the Npop most stable, unique structures. We have set Npop to 20 and evaluate uniqueness using fingerprint functions (see Ref. 11). Unless the maximum number of GA iterations Imax has been reached, a new candidate is created by one of the genetic operators which are described in Table 1. For carbon, titanium dioxide and palladium oxide, the operators are chosen with probabilities 25/30/30/15 (pairing, soft mutation, strain mutation, rattle mutation). For calcium and palladium hydroxide, the probabilities are 25/25/25/25 (pairing, soft mutation, strain mutation, rattle-rotation mutation). The required parents (1 or 2 depending on whether a mutation or pairing operator is selected) are randomly chosen from the population with a bias towards the more stable and less frequently selected candidates, as described in Ref. 16. As the genetic operations regularly produce configurations which contain short interatomic distances, the atoms in the generated structure are expanded so that all distances exceed a specified minimum bond length (e.g. half the sum of the covalent radii). After this step, the local minimization is carried out until the largest forces fall below 0.1 eV/˚ A and all stress components are less than 5 · 10−3 eV/˚ A3 . The resulting structure is thereafter added to the GA database.

5

DFT setup

The DFT calculations for the TANGO searches are carried out using the Quickstep module in CP2K 58 (version 5.0, git revision 4f7b2ef). We employ GTH pseudopotentials 59–61 with standard valence states for C, H and O and additional s and p semicore states for Ca, Ti and Pd. The Kohn-Sham wave functions are expanded in a shorter-ranged double-ζ basis set

12


Page 12 of 36

Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with polarization functions from the molecularly optimized (MOLOPT) basis set series. 62 The electron density is represented using plane waves with kinetic energies up to 600 Ry. A relative energy cutoff of 60 Ry is used to control the accuracy with which the Gaussian functions are mapped onto the real-space multi-grids. Such relatively high cutoff energies are required for adequate convergence of total energy differences, as the plane wave basis set changes with the size and shape of the unit cells. Monkhorst-Pack grids 63,64 are used for sampling the first Brillouin zone with k-point ˚−1 ) of 2.5 points carbon, 1.5 for titanium dioxide, 3.0 densities dk (in units of points per A for palladium oxide, 1.5 for calcium hydroxide and 2.5 for palladium hydroxide. The number of k-points along a given reciprocal cell vector (associated with a unit cell vector of length L) is then given by 2π dk /L. The exchange-correlation energy is approximated using the Perdew-Burke-Ernzerhof (PBE) functional. 65 It should be noted that generalized-gradient functionals such as PBE may yield an incorrect energetical ordering of the polymorphs, as will be apparent in the case of titanium dioxide. For the purpose of testing the present approach, however, the choice of the functional is not crucial. In the case of palladium hydroxide, we have also considered the GPAW method 66,67 where the plane-wave basis set and PAW potentials in principle offer increased technical accuracy. An 800 eV plane-wave energy cutoff is used in the GPAW calculations and both the PBE and BEEF-vdW functionals are evaluated.

6

Results

In the following paragraphs, the TANGO approach is applied to a series of crystal structures, namely carbon, titanium dioxide, palladium oxide, calcium hydroxide and palladium hydroxide. The considered cells include between 4 and 12 formula units, yielding a total number of atoms between 12 and 20. We have used sixth-order polynomials for all element pairs,

13



Page 14 of 36

Table 2: Electronic configurations for the different elements as well as their `-dependent eigenvalues ` and Hubbard parameters U` (in eV) calculated with the PBE functional. The rightmost column contains the covalent radii rcov in ˚ A (from Ref. 68) which are used to calculate e.g. the parabolic confinement radii (as 1.85 × rcov ) and the cutoff radii for the A B repulsive interactions (as 1.5 × (rcov + rcov )). Element Core Valence H 1s1 C [He] 2s2 2p2 O [He] 2s2 2p4 Ca [Ar] 4s2 2 Ti [Ar] 4s 3d2 4p0 Pd [Kr] 4d9 5s1 5p0

`=0 -6.493 -13.748 -23.956 -3.766 -4.466 -4.259

` `=1 `=2 -5.284 -9.028 -1.393 -4.264 -0.407 -6.605

`=0 14.339 17.428 28.310 5.548 6.927 8.794

U` `=1 `=2 11.872 19.461 44.903 10.175 59.027 14.069

rcov 0.31 0.76 0.66 1.76 1.60 1.39

except for the metal hydroxides where cubic polynomials are employed for H-H, O-O and H-O. For these pairs, higher-order polynomial fits would suffer from overfitting as a result of the relatively scarcity and high energy of configurations with short bond distances between these elements. Additionally, configurations with O-O and H-H bonds less than 1 ˚ A are discarded during the associated GA searches. For each material, 20 random structures (partially relaxed with DFT) and 5 batches of 20 random structures (relaxed with DFTB) are used to obtain the initial repulsive potentials (see Section 3 and Figure 1). Twenty independent GO runs are then performed using the GA described in Section 4 in 2 steps with each 50 GA iterations. In between these steps, the repulsive potentials are further refined using the 100 best unique DFTB structures found thus far. The final set of coefficients for the repulsive interactions are shown in Table S1 and Figure S2 of the Supporting Information. The employed electronic DFTB parameters for the various elements are listed in Table 2.

6.1

Carbon

Carbon compounds have frequently been analyzed using tight-binding models, and many bulk polymorphs are known. Panel (a) in Figure 3 shows a parity plot of the final 100 best unique structures gathered from the TANGO calculation. The relative ordering according to 14


Page 15 of 36

(a) Carbon (C12)

(b) Titanium dioxide (Ti4O8)

Anatase (A)

Graphite (G) GD

A R Rutile (R)

Diamond (D) (d) Calcium hydroxide (Ca4(OH)8)

(c) Palladium oxide (Pd6O6)

P

}

P }

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Palladinite (P)

Portlandite (P)

Figure 3: Parity plots and relevant crystal structures for the four different test cases. The plots show the relative DFT and DFTB energies per formula unit for the 100 best, unique structures from the DFTB GA runs. The reference energies correspond to the structure with the lowest DFT energy. The models have been rendered with VESTA. 69 The H, C and O atoms are colored white, black and red, respectively. The different metal atoms (Ca, Ti, Pd) are colored blue-green. the DFT energies is well reproduced by the final DFTB parametrization. Graphite (with various stackings of the graphene sheets) and diamond are the two types of structures with the lowest DFT (single-point) energies. These indeed correspond to the PBE GM and the second-best structures. In addition, several other relevant carbon polymorphs are found among the lowest-energy structures, such as lonsdaleite, graphite with 7- and 5-membered rings, and a chiral pentagon-only diamond structure. 70 Figure 4 displays the repulsive potential between the carbon atoms at two different stages of the TANGO run. The first parametrization in this case already produces a qualitatively representative repulsive potential, which is then further refined during the subsequent reparametrizations. Note that the repulsive potential is allowed to adopt moderately 15



First fit Last fit

40 Vrep, C−C [eV]

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

30 20 10

rcut

rexp

0 1

1.5 ˚ rC−C [A]

2

Figure 4: The C-C repulsive potential at two different stages of the C12 TANGO run. The training set for the first fit (dashed red line) only includes the 20 partially relaxed DFT structures. In the last fit (filled blue-green line) also 5 batches of 20 random structures (relaxed with DFTB) are included in the training set, as well as the best 100 unique structures obtained from 20 GA runs with 50 iterations each. negative values, as it is treated in an empirical manner.

6.2

Titanium dioxide

Titanium dioxide is an often-used test case for the development of GO methods for reason of its rich polymorphism. The corresponding parity plot is shown in panel (b) of Figure 3, with indication of the anatase and rutile phases. Even though both structures are recovered among the configurations with lowest DFT energies, in this instance the TiO2 (B) polymorph is the best structure found at this point. This indicates that DFT reoptimization is necessary in this case to obtain the true PBE ordering, where anatase and rutile are the preferred phases2 . Other TiO2 polymorphs are also encountered, such as brookite, columbite and lepidocrocite-like nanosheets. 2

GGA as well as hybrid functionals wrongly predict anatase to be more stable than rutile. 71

16


Page 16 of 36

Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6.3

Palladium oxide

Palladium oxide represents an interesting case, as so far no DFTB parametrizations have been reported for Pd-O compounds. It should be noted that the PBE functional predicts PdO to be metallic 72 whereas experimentally PdO is a semiconductor with a band gap of around 1 eV. 73 The tetragonal palladinite structure is recovered as the best structure from the TANGO run, according to both DFTB and the single-point DFT energies (see panel (c) in Figure 3). Rocksalt-like polymorphs have been proposed to form under high pressure 74 or as supported ultrathin layers. 75 Such structures are also retrieved here, however not among the set of 100 best unique structures. This appears to be consistent with the high relative energy with respect to the tetragonal structure (circa 0.85 eV per formula unit at the GGA-DFT level 76 ).

6.4

Calcium hydroxide

The results for calcium hydroxide are shown in panel (d) of Figure 3. Both according to the DFTB energies and the DFT single-point energies, the known portlandite mineral is obtained as the preferred candidate. Each of the five best structures indicated in the parity plot correspond to this brucite-type structure, with either perfectly aligned or slightly perturbed orientations of the hydroxide groups. Other polymorphs, such as the two highpressure phases (see Ref. 77) are absent in this set. As with rocksalt PdO, finding these structures may require the addition of a pV term to the potential energy.

6.5

Palladium hydroxide

The final application focuses on palladium hydroxide, for which no crystal structure has yet been proposed. Its structure is relevant to palladium catalysis, in particular at high H2 O chemical potentials. In several studies, the formation of unreactive Pd(OH)2 has been suggested to inhibit methane oxidation over PdO catalysts. 79–84 Other work has found the

17



Page 18 of 36

HP3

HP5

CP2K PBE GPAW PBE GPAW BEEF-vdW

Figure 5: Panel (a): parity plot comparing the relative DFT and DFTB energies per formula unit of the best 100 unique Pd4 (OH)8 structures gathered from the GA runs. The reference energies correspond to the structure with the lowest DFT energy. Panel (b): DFT energy differences of the 25 best structures from single-point calculations at the DFTB geometries (black line) and after DFT reoptimization (blue-green line). The energy of the first relaxed structure is used as the reference. Panel (c): layered representations of hydrated PdO (based on Ref. 78). Panel (d): energetic comparison of the HP3 and the 25 best structures with respect to HP5 using different computational setups. activity of PdO for CO oxidation to increase upon hydroxide formation. 85 In organic chemistry, Pearlman’s catalyst is frequently described as Pd(OH)2 on carbon. 86 However, there exists substantial evidence that Pd(OH)2 consists of hydrated palladium oxide rather than a bulk hydroxide. 78,87 One useful piece of information in this context would be to know the relative stability of Pd(OH)2 and hydrated PdO. To this end, we compare a set of Pd(OH)2 structures obtained from a TANGO search with previously proposed structures for hydrated PdO. 78 The results using a Pd4 (OH)8 stoichiometry are shown in Figure 5. The parity plot

18


Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: Structural models of the most stable variations of three selected structural classes (I, II, II) of palladium hydroxide. The atomic color coding is as in Figure 3. Hydrogen bonds are indicated as dashed gray lines. in panel (a) displays significantly more scatter than the previous test cases, meaning that obtaining an accurate parametrization is more challenging for this compound. Nevertheless plausible configurations are recovered among the final 100 best unique structures. The 25 best candidates are reoptimized at the DFT level, and the corresponding energy changes are shown in panel (b). Understandably, this reoptimization leads to a certain reordering of the relative stabilities. Panel (c) displays the hydrated PdO models (based on Ref. 78) with which the hydroxide structures will be compared. These consist of hydroxylated PdO(001) slabs with additional H2 O layers. The stoichiometries correspond to Pd3 (OH)6 (‘HP3’) and Pd5 (OH)10 (‘HP5’). It should be reminded that the present TANGO calculations are performed with a Pd4 (OH)8 stoichiometry and that the search space is explicitly confined to hydroxide compounds, thereby excluding the disproportionation of OH groups and the formation of hydrated PdO phases. In panel (d), the stabilities of the 25 Pd(OH)2 structures with respect to HP3 and HP5 are evaluated after additional reoptimizations using different theoretical descriptions3 . Using the PBE functional with the employed CP2K setup, the HP5 structure is found to represent the energetically most stable phase, and the most stable hydroxide structures are competitive with the HP3 model. We have then reoptimized each structure using the GPAW code, with the same choice of functional (PBE). The ordering of the hydroxide structures is not We have increased the k-point density to 4 points per ˚ A−1 and the (relative) energy cutoff to 1000 (100) Ry to ensure convergence of the total energy. 3

19



Page 20 of 36

significantly perturbed, though the HP phase is destabilized and is no longer the preferred structure. The presence of separate water layers, however, suggests that dispersion interactions can be relevant to the relative stability of hydrated PdO. Using the BEEF-vdW functional 88 in GPAW, the ordering is indeed shifted to the extent that the HP3 and HP5 configurations become most energetically favored. From the calculated single-point energies Esp and relaxation energies ∆Erelax , it is furthermore possible to estimate the probability P that other, higher-energy structures would surpass the currently obtained GM structure after DFT relaxation. Based on this information, one may judge whether the initial batch of relaxations (here consisting of 25 runs) is sufficient or additional DFT relaxations are needed to reach a certain level of certainty. Under the assumption that the relaxation energies are normally distributed, this probability P is given by: P =1−

unrelaxed Y i

Φ

Esp,i − EGM − µ ˜ σ ˜

,

(9)

with Φ the cumulative distribution function of the standard normal distribution, and µ ˜ and σ ˜ the estimated mean and standard deviation of the available relaxation energies. In this example, we obtain P = 0.06, meaning that there is a very low (6 %) probability to improve on the current GM structure by relaxing the 75 other best unique structures. The most competitive hydroxide phases can be grouped in three classes (I, II, III), of which the most stable representatives are shown in Figure 6. The coordinates and structural models of each structure are included in the Supporting Information. Structures within each class display relatively minor changes, such as different orientations of the hydroxide groups. Interestingly, the most stable motif (I) consists of linear Pd(OH)2 chains with hydrogen bonding between adjacent chains. This arrangement is in fact very similar to that of αPdCl2 , which is consistent with the notion of OH as a pseudohalogen. The second motif (II) is of a more conventional type, reminiscent of the brucite structure common among metal hydroxides. Motif III is also composed of two-dimensional layers, however with a square arrangement of the metal ions (as opposed to the hexagonal ordering in brucite). In all 20


Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


cases, the palladium atoms are coordinated in a square-planar fashion to the neighboring oxygen atoms, which is indeed a common occurrence in palladium compounds. Within the present approximations, the assignment of Pd(OH)2 to nanocrystalline, hydrated PdO is hence preferred over a bulk hydroxide compound. Presumably, a 7-layer slab (with Pd7 (OH)14 stoichiometry) would be an even more stable representation of hydrated PdO. It remains to be seen, however, to what extent e.g. vibrational effects and higher-level functionals affect the free energy difference between these two forms.

7

Discussion

In this section we provide additional discussions of the merits and possibilities of the TANGO method, as well as ideas for further improvements.

7.1

Benefits of the approach

The main merit of the approach is that it allows to locate the low-lying local minima (including the GM) of the first-principles energy landscape at much reduced computational cost. While the overall speedup is difficult to evaluate explicitly, we estimate the total CPU time to be reduced by a factor of at least fifty to hundred when using semi-local functionals, as compared to GO at only the DFT level. It should be noted that GO problems cannot be efficiently tackled in a massively parallel fashion, meaning that the time-to-solution cannot be strongly reduced by simply using a large amount of computational resources. Parallelism can for example be introduced by evaluating more than one candidate simultaneously. Efficient speedup, however, only allows limited parallelism because typically an increasingly fit population is needed to generate increasingly better candidates. Efficient parallelization of the electronic structure calculations themselves is also limited, for different reasons. Enhancement by low-cost semi-empirical methods such as DFTB is therefore valuable even when vast computational resources are

21



available. While the structures obtained from TANGO calculations should, in general, be reoptimized at the first-principles level, these relaxations are entirely independent and can therefore be performed simultaneously. Naturally, the search would become even more cost effective if an appropriate DFTB parametrization is already available from the start, which, however, is generally not the case. Transition metal compounds, in particular, have only in a few cases been investigated with DFTB. Additionally, certain parametrizations may have been generated on a substantially different class of compounds and, as a result, the transferability to the problem at hand may be limited. Training sets focusing on (transition) metal complexes, 48,89 for example, may yield parameters unfit for more bulk-like compounds. Similarly, the recently developed GFN-xTB approach 90 is intended for molecular chemistry and does not (as of yet) handle periodic boundary conditions or metal compounds other than complexes with few metal atoms. Likewise, parameters in the QUASINANO set 32,33 are only available up to Ca and are unlikely to offer the accuracy needed for reliable GO searches (though it can possibly be used to generate the initial training set). Even though we have chosen bulk crystals as the focus of our tests, this method is equally applicable to other types of materials. Relevant applications include, for example, the structure of surface terminations and other interfaces, (supported) nanoparticles, or substances inside mesoporous materials such as zeolites and metal-organic frameworks. Even though further testing will be required to ascertain the suitability of the method for such applications, we expect adequate performance as we consider the present test cases to be fairly challenging due to a large variety of bonding configurations (coordinations, distances and bond angles). Lastly, it should be noted that the present approach is not limited to GA as the search algorithm. Other GO strategies such as particle-swarm optimization, minima hopping, or simulated annealing could, in principle, be used.

22


Page 22 of 36

Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


7.2

Potential improvements

Although the method is found to perform adequately for the studied materials, the parity plots in Figures 3 and 5 suggest that it may be beneficial to reach more accurate DFTB parametrizations. As only the repulsive potentials are fitted in the present implementation, a higher accuracy may be attained by also tuning the atomic eigenvalues, Hubbard parameters, confinement radii, and/or cutoff radii. It would also be a possibility to employ other functional forms of the repulsive pair potentials, such as B-splines. Furthermore, the repulsion may be described more accurately through the use of environment-dependent screening functions, as has been done for molybdenum 91 and high-pressure carbon. 92 Other options include adding third-order corrections to DFTB 93 or adding polarization functions. In case the employed first-principles method is capable of describing van der Waals interactions, naturally also dispersion corrections should be added to the DFTB parametrization. Another path to further improving the efficiency of the method would be to perform most of the initial parametrization using smaller, yet still representative systems. In the studied examples, one could initially include less formula units (e.g. Pd(OH)2 or Pd2 (OH)4 instead of Pd4 (OH)8 ). Lastly, as for other adaptive approaches, there may be a risk of ‘agnosis’. In principle, structures at and/or near the GM may display a very different but favorable bonding pattern which is unstable in the initial DFTB parametrization. Such patterns are then unlikely to appear in stable DFTB structures, and would then never be sampled with DFT. We expect this to be a minor risk, in particular if the initial set of random structures (partially relaxed with DFT) is relevant to the problem at hand. This risk can, in principle, be lowered by including additional random or even custom-made structures to the initial training phase.

23



8

Conclusions

We have developed and tested an adaptive global optimization strategy involving iterative refinement of an SCC-DFTB description of the system under study. The approach strongly reduces the computational cost of first-principles structure prediction while providing adequate reliability, at least for the benchmarks investigated thus far (bulk carbon, titanium dioxide, palladium oxide and calcium hydroxide). This has been possible to achieve despite the fact that only the shape of the repulsive potentials are empirically adjusted. Applying this method to the crystal structure of Pd(OH)2 suggests an α-PdCl2 -type structure as the ground state. The formation of this structure is, within the applied approximations, slightly endothermic with respect to that of nanocrystalline hydrated PdO. We furthermore identify possible avenues for further methodological improvement. We hope this approach will prove to be useful for various problems in material science and development and will encourage the use of global optimization techniques. The TANGO code will be made freely available at http://gitlab.com/mvdb/tango and the implemented genetic operators will be incorporated into the ASE trunk.

Acknowledgement The authors thank the Swedish Research Council (2016-05234), Chalmers Areas of Advance Nano and Transport and COST action CM1104 for financial support. The calculations were performed at C3SE (Göteborg) via a SNIC grant. The Competence Centre for Catalysis (KCK) is hosted by Chalmers University of Technology and is financially supported by the Swedish Energy Agency and the member companies: AB Volvo, ECAPS AB, Haldor Topsøe A/S, Scania CV AB, Volvo Car Corporation AB, and Wärtsilä Finland Oy. The authors furthermore acknowledge support from the Danish Council for Independent Research Natural Sciences (grant no. 0602-02566B) and from VILLUM FONDEN (Investigator grant, project no. 16562). 24


Page 24 of 36

Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Supporting Information Available Structural models of the Pd(OH)2 structures as well as plots and parameters of the final repulsive potentials obtained for the studied materials. Coordinate files for the different variations of the class I, II and III Pd(OH)2 structures and the hydrated PdO models.

References (1) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864– B871. (2) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133–A1138. (3) Heard, C. J.; Heiles, S.; Vajda, S.; Johnston, R. L. Pdn Ag(4-n) and Pdn Pt(4-n) clusters on MgO (100): a density functional surface genetic algorithm investigation. Nanoscale 2014, 6, 11777–11788. (4) Trinchero, A.; Klacar, S.; Paz-Borbón, L. O.; Hellman, A.; Grönbeck, H. Oxidation at the Subnanometer Scale. J. Phys. Chem. C 2015, 119, 10797–10803. (5) Merte, L. R.; Jørgensen, M. S.; Pussi, K.; Gustafson, J.; Shipilin, M.; Schaefer, A.; Zhang, C.; Rawle, J.; Nicklin, C.; Thornton, G.; Lindsay, R.; Hammer, B.; Lundgren, E. Structure of the SnO2 (110) – (4x1) Surface. Phys. Rev. Lett. 2017, 119, 096102. (6) Zhang, R. F.; Wen, X. D.; Legut, D.; Fu, Z. H.; Veprek, S.; Zurek, E.; Mao, H. K. Crystal Field Splitting is Limiting the Stability and Strength of Ultra-incompressible Orthorhombic Transition Metal Tetraborides. Sci. Rep. 2016, 6, 23088. (7) Martinez, U.; Vilhelmsen, L. B.; Kristoffersen, H. H.; Stausholm-Møller, J.; Hammer, B. Steps on rutile TiO2 (110): Active sites for water and methanol dissociation. Phys. Rev. B 2011, 84, 205434. 25



(8) Packwood, D.; Kermode, J.; Mones, L.; Bernstein, N.; Woolley, J.; Gould, N.; Ortner, C.; Csányi, G. A universal preconditioner for simulating condensed phase materials. J. Chem. Phys. 2016, 144, 164109. (9) Glass, C. W.; Oganov, A. R.; Hansen, N. USPEX – Evolutionary crystal structure prediction. Comput. Phys. Commun. 2006, 175, 713–720. (10) Oganov, A. R.; Lyakhov, A. O.; Valle, M. How Evolutionary Crystal Structure Prediction Works – and Why. Acc. Chem. Res. 2011, 44, 227–237. (11) Lyakhov, A. O.; Oganov, A. R.; Stokes, H. T.; Zhu, Q. New developments in evolutionary structure prediction algorithm USPEX. Comput. Phys. Commun. 2013, 184, 1172–1182. (12) Lonie, D. C.; Zurek, E. XtalOpt: An open-source evolutionary algorithm for crystal structure prediction. Comput. Phys. Commun. 2011, 182, 372–387. (13) Lonie, D. C.; Zurek, E. Identifying duplicate crystal structures: XtalComp, an opensource solution. Comput. Phys. Commun. 2012, 183, 690–697. (14) Lloyd, L. D.; Johnston, R. L.; Salhi, S. Strategies for increasing the efficiency of a genetic algorithm for the structural optimization of nanoalloy clusters. J. Comput. Chem. 2005, 26, 1069–1078. (15) Heiles, S.; Johnston, R. L. Global optimization of clusters using electronic structure methods. Int. J. Quantum Chem. 2013, 113, 2091–2109. (16) Vilhelmsen, L. B.; Hammer, B. A genetic algorithm for first principles global structure optimization of supported nano structures. J. Chem. Phys. 2014, 141, 044711. (17) Jørgensen, M. S.; Groves, M. N.; Hammer, B. Combining Evolutionary Algorithms with Clustering toward Rational Global Structure Optimization at the Atomic Scale. J. Chem. Theory Comput. 2017, 13, 1486–1493. 26


Page 26 of 36

Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Jacobsen, T.; Jørgensen, M.; Hammer, B. On-the-Fly Machine Learning of Atomic Potential in Density Functional Theory Structure Optimization. Phys. Rev. Lett. 2018, 120, 026102. (19) Hartke, B. Global optimization. WIREs Comput Mol Sci 2011, 1, 879–887. (20) Tipton, W. W.; Hennig, R. G. A grand canonical genetic algorithm for the prediction of multi-component phase diagrams and testing of empirical potentials. J. Phys.: Condens. Matter 2013, 25, 495401. (21) Hartke, B. Global geometry optimization of clusters guided by N-dependent model potentials. J. Chem. Phys. Lett. 1996, 258, 144–148. (22) Tekin, A.; Hartke, B. Global geometry optimization of small silicon clusters with empirical potentials and at the DFT level. Phys. Chem. Chem. Phys. 2004, 6, 503–509. (23) Wu, S. Q.; Ji, M.; Wang, C. Z.; Nguyen, M. C.; Zhao, X.; Umemoto, K.; Wentzcovitch, R. M.; Ho, K. M. An adaptive genetic algorithm for crystal structure prediction. J. Phys. Condens. Matter 2014, 26, 035402. (24) Zhao, X.; Cuong Nguyen, M.; Wang, C.-Z.; Ho, K.-M. Structures and stabilities of alkaline earth metal peroxides XO2 (X = Ca, Be, Mg) studied by a genetic algorithm. RSC Adv. 2013, 3, 22135–22139. (25) Zhao, X.; Shu, Q.; Nguyen, M. C.; Wang, Y.; Ji, M.; Xiang, H.; Ho, K.-M.; Gong, X.; Wang, C.-Z. Interface Structure Prediction from First-Principles. J. Phys. Chem. C 2014, 118, 9524–9530. (26) E. L. Kolsbjerg, A. A. Peterson, B. Hammer. Neural Network Enhanced Evolutionary Algorithm Applied to Supported Metal Clusters. To be published . (27) Porezag, D.; Frauenheim, T.; Köhler, T.; Seifert, G.; Kaschner, R. Construction of

27



tight-binding-like potentials on the basis of density-functional theory: Application to carbon. Phys. Rev. B 1995, 51, 12947–12957. (28) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties. Phys. Rev. B 1998, 58, 7260–7268. (29) Gaus, M.; Cui, Q.; Elstner, M. Density functional tight binding: application to organic and biological molecules. WIREs Comput Mol Sci 2014, 4, 49–61. (30) Slater, J. C.; Koster, G. F. Simplified LCAO Method for the Periodic Potential Problem. Phys. Rev. 1954, 94, 1498–1524. (31) Mulliken, R. S. Electronic Population Analysis on LCAO-MO Molecular Wave Functions. I. J. Chem. Phys. 1955, 23, 1833–1840. (32) Wahiduzzaman, M.; Oliveira, A. F.; Philipsen, P.; Zhechkov, L.; van Lenthe, E.; Witek, H. A.; Heine, T. DFTB Parameters for the Periodic Table: Part 1, Electronic Structure. J. Chem. Theory Comput. 2013, 9, 4006–4017. (33) Oliveira, A. F.; Philipsen, P.; Heine, T. DFTB Parameters for the Periodic Table, Part 2: Energies and Energy Gradients from Hydrogen to Calcium. J. Chem. Theory Comput. 2015, 11, 5209–5218. (34) Chou, C.-P.; Nishimura, Y.; Fan, C.-C.; Mazur, G.; Irle, S.; Witek, H. A. Automatized Parameterization of DFTB Using Particle Swarm Optimization. J. Chem. Theory Comput. 2016, 12, 53–64. (35) Gaus, M.; Chou, C.-P.; Witek, H.; Elstner, M. Automatized Parametrization of SCCDFTB Repulsive Potentials: Application to Hydrocarbons. J. Phys. Chem. A 2009, 113, 11866–11881.

28


Page 28 of 36

Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(36) Loureno, M. P.; Silva, M. C. d.; Oliveira, A. F.; Quinto, M. C.; Duarte, H. A. FASP: a framework for automation of Slater-Koster file parameterization. Theor. Chem. Acc. 2016, 135, 250. (37) Elstner, M.; Seifert, G. Density functional tight binding. Phil. Trans. R. Soc. A 2014, 372, 20120483. (38) Sch¨ utt, O.; Messmer, P.; Hutter, J.; VandeVondele, J. In Electronic Structure Calculations on Graphics Processing Units; Walker, R. C., Götz, A. W., Eds.; John Wiley & Sons, Ltd, 2016; pp 173–190. (39) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. CP2K: atomistic simulations of condensed matter systems. WIREs Comput Mol Sci 2014, 4, 15–25. (40) Borˇstnik, U.; VandeVondele, J.; Weber, V.; Hutter, J. Sparse matrix multiplication: The distributed block-compressed sparse row library. Parallel Comput. 2014, 40, 47– 58. (41) Hess, B. A. Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators. Phys. Rev. A 1986, 33, 3742–3748. (42) Jansen, G.; Hess, B. A. Revision of the Douglas-Kroll transformation. Phys. Rev. A 1989, 39, 6016–6017. (43) Douglas, M.; Kroll, N. M. Quantum electrodynamical corrections to the fine structure of helium. Ann. Phys. 1974, 82, 89 – 155. (44) Koskinen, P.; Mäkinen, V. Density-functional tight-binding for beginners. Comp. Mater. Sci. 2009, 47, 237–253. (45) Aradi, B.; Hourahine, B.; Frauenheim, T. DFTB+, a Sparse Matrix-Based Implementation of the DFTB Method. J. Phys. Chem. A 2007, 111, 5678–5684. 29



(46) Bodrog, Z.; Aradi, B.; Frauenheim, T. Automated Repulsive Parametrization for the DFTB Method. J. Chem. Theory Comput. 2011, 7, 2654–2664. (47) Goldman, N.; Fried, L. E.; Koziol, L. Using Force-Matched Potentials To Improve the Accuracy of Density Functional Tight Binding for Reactive Conditions. J. Chem. Theory Comput. 2015, 11, 4530–4535. (48) Zheng, G.; Witek, H. A.; Bobadova-Parvanova, P.; Irle, S.; Musaev, D. G.; Prabhakar, R.; Morokuma, K.; Lundberg, M.; Elstner, M.; Köhler, C.; Frauenheim, T. Parameter Calibration of Transition-Metal Elements for the Spin-Polarized SelfConsistent-Charge Density-Functional Tight-Binding (DFTB) Method: Sc, Ti, Fe, Co, and Ni. J. Chem. Theory Comput. 2007, 3, 1349–1367. (49) Dolgonos, G.; Aradi, B.; Moreira, N. H.; Frauenheim, T. An Improved Self-ConsistentCharge Density-Functional Tight-Binding (SCC-DFTB) Set of Parameters for Simulation of Bulk and Molecular Systems Involving Titanium. J. Chem. Theory Comput. 2010, 6, 266–278. (50) Lyakhov, A. O.; Oganov, A. R.; Valle, M. How to predict very large and complex crystal structures. Comput. Phys. Commun. 2010, 181, 1623–1632. (51) Vilhelmsen, L. B.; Hammer, B. Systematic Study of Au6 to Au12 Gold Clusters on MgO(100) F Centers Using Density-Functional Theory. Phys. Rev. Lett. 2012, 108, 126101. (52) Bahn, S. R.; Jacobsen, K. W. An object-oriented scripting interface to a legacy electronic structure code. Comput. Sci. Eng. 2002, 4, 56–66. (53) Larsen, A. H.; Mortensen, J. J.; Blomqvist, J.; Castelli, I. E.; Christensen, R.; Duak, M.; Friis, J.; Groves, M. N.; Hammer, B.; Hargus, C.; Hermes, E. D.; Jennings, P. C.; Jensen, P. B.; Kermode, J.; Kitchin, J. R.; Kolsbjerg, E. L.; Kubal, J.; Kristen Kaasbjerg,; Lysgaard, S.; Maronsson, J. B.; Maxson, T.; Olsen, T.; Pastewka, L.; Andrew 30


Page 30 of 36

Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Peterson,; Rostgaard, C.; Schiøtz, J.; Sch¨ utt, O.; Strange, M.; Thygesen, K. S.; Tejs Vegge,; Vilhelmsen, L.; Walter, M.; Zeng, Z.; Jacobsen, K. W. The atomic simulation environmenta Python library for working with atoms. J. Phys. Condens. Matter 2017, 29, 273002. (54) Zhu, Q.; Oganov, A. R.; Glass, C. W.; Stokes, H. T. Constrained evolutionary algorithm for structure prediction of molecular crystals: methodology and applications. Acta Cryst. B 2012, 68, 215–226, arXiv: 1204.4756. (55) Liu, D. C.; Nocedal, J. On the limited memory BFGS method for large scale optimization. Math. Program. 1989, 45, 503–528. (56) Bitzek, E.; Koskinen, P.; Gähler, F.; Moseler, M.; Gumbsch, P. Structural Relaxation Made Simple. Phys. Rev. Lett. 2006, 97, 170201. (57) Armijo, L. Minimization of functions having Lipschitz continuous first partial derivatives. Pacific J. Math. 1966, 16, 1–3. (58) VandeVondele, J.; Krack, M.; Mohamed, F.; Parrinello, M.; Chassaing, T.; Hutter, J. Quickstep: Fast and accurate density functional calculations using a mixed Gaussian and plane waves approach. Comput. Phys. Commun. 2005, 167, 103–128. (59) Goedecker, S.; Teter, M.; Hutter, J. Separable dual-space Gaussian pseudopotentials. Phys. Rev. B 1996, 54, 1703–1710. (60) Hartwigsen, C.; Goedecker, S.; Hutter, J. Relativistic separable dual-space Gaussian pseudopotentials from H to Rn. Phys. Rev. B 1998, 58, 3641–3662. (61) Krack, M. Pseudopotentials for H to Kr optimized for gradient-corrected exchangecorrelation functionals. Theor. Chem. Acc. 2005, 114, 145–152. (62) VandeVondele, J.; Hutter, J. Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases. J. Chem. Phys. 2007, 127, 114105. 31



(63) Monkhorst, H. J.; Pack, J. D. Special points for Brillouin-zone integrations. Phys. Rev. B 1976, 13, 5188–5192. (64) Pack, J. D.; Monkhorst, H. J. ‘Special points for Brillouin-zone integrations’ – a reply. Phys. Rev. B 1977, 16, 1748–1749. (65) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (66) Mortensen, J. J.; Hansen, L. B.; Jacobsen, K. W. Real-space grid implementation of the projector augmented wave method. Phys. Rev. B 2005, 71, 035109. (67) Enkovaara, J.; Rostgaard, C.; Mortensen, J. J.; Chen, J.; Duak, M.; Ferrighi, L.; Gavnholt, J.; Glinsvad, C.; Haikola, V.; Hansen, H. A.; Kristoffersen, H. H.; Kuisma, M.; Larsen, A. H.; Lehtovaara, L.; Ljungberg, M.; Lopez-Acevedo, O.; Moses, P. G.; Ojanen, J.; Olsen, T.; Petzold, V.; Romero, N. A.; Stausholm-Møller, J.; Strange, M.; Tritsaris, G. A.; Vanin, M.; Walter, M.; Hammer, B.; Häkkinen, H.; Madsen, G. K. H.; Nieminen, R. M.; Nørskov, J. K.; Puska, M.; Rantala, T. T.; Schiøtz, J.; Thygesen, K. S.; Jacobsen, K. W. Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method. J. Phys. Condens. Matter 2010, 22, 253202. (68) Cordero, B.; Gómez, V.; Platero-Prats, A. E.; Revés, M.; Echeverr´ıa, J.; Cremades, E.; Barragán, F.; Alvarez, S. Covalent radii revisited. Dalton Trans. 2008, 0, 2832–2838. (69) Momma, K.; Izumi, F. VESTA 3 for Three-Dimensional Visualization of Crystal, Volumetric and Morphology Data. J. Appl. Crystallogr. 2011, 44, 1272–1276. (70) Zhu, X.; Su, H. Chiral Pentagon Only Diamond-like Structures. J. Phys. Chem. C 2017, 121, 13810–13815.

32


Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(71) Labat, F.; Baranek, P.; Domain, C.; Minot, C.; Adamo, C. Density functional theory analysis of the structural and electronic properties of TiO2 rutile and anatase polytypes: Performances of different exchange-correlation functionals. J. Chem. Phys. 2007, 126, 154703. (72) Seriani, N.; Harl, J.; Mittendorfer, F.; Kresse, G. A first-principles study of bulk oxide formation on Pd(100). J. Chem. Phys. 2009, 131, 054701. (73) Nilsson, P. O. Optical Properties of PdO in the Range of 0.5-5.4 eV. J. Phys. C: Solid State Phys. 1979, 12, 1423. (74) Christy, A. G.; Clark, S. M. Structural behavior of palladium (II) oxide and a palladium suboxide at high pressure: An energy-dispersive x-ray-diffraction study. Phys. Rev. B 1995, 52, 9259–9265. (75) Kumar, J.; Saxena, R. Formation of NaCl- and Cu2 O-type oxides of platinum and palladium on carbon and alumina support films. J. Less Common Met. 1989, 147, 59–71. (76) Derzsi, M.; Piekarz, P.; Grochala, W. Structures of Late Transition Metal Monoxides from Jahn-Teller Instabilities in the Rock Salt Lattice. Phys. Rev. Lett. 2014, 113, 025505. (77) Iizuka, R.; Komatsu, K.; Kagi, H.; Nagai, T.; Sano-Furukawa, A.; Hattori, T.; Gotou, H.; Yagi, T. Phase transitions and hydrogen bonding in deuterated calcium hydroxide: High-pressure and high-temperature neutron diffraction measurements. J. Solid State Chem. 2014, 218, 95–102. (78) Parker, S. F.; Refson, K.; Hannon, A. C.; Barney, E. R.; Robertson, S. J.; Albers, P. Characterization of Hydrous Palladium Oxide: Implications for Low-Temperature Carbon Monoxide Oxidation. J. Phys. Chem. C 2010, 114, 14164–14172.

33



(79) Cullis, C. F.; Nevell, T. G.; Trimm, D. L. Role of the catalyst support in the oxidation of methane over palladium. J. Chem. Soc. Faraday Trans. 1 1972, 68, 1406–1412. (80) Cullis, C. F.; Willatt, B. M. The inhibition of hydrocarbon oxidation over supported precious metal catalysts. J. Catal. 1984, 86, 187–200. (81) Ribeiro, F. H.; Chow, M.; Dallabetta, R. A. Kinetics of the Complete Oxidation of Methane over Supported Palladium Catalysts. J. Catal. 1994, 146, 537–544. (82) Burch, R.; Urbano, F.; Loader, P. Methane combustion over palladium catalysts: The effect of carbon dioxide and water on activity. Applied Catalysis A: General 1995, 123, 173–184. (83) Hayes, R. E.; Kolaczkowski, S. T.; Li, P. K.; Awdry, S. The palladium catalysed oxidation of methane: reaction kinetics and the effect of diffusion barriers. Chem. Eng. Sci. 2001, 56, 4815–4835. (84) Ciuparu, D.; Katsikis, N.; Pfefferle, L. Temperature and time dependence of the water inhibition effect on supported palladium catalyst for methane combustion. Appl. Catal. A 2001, 216, 209–215. (85) Oh, S.-H.; Hoflund, G. B. Low-temperature catalytic carbon monoxide oxidation over hydrous and anhydrous palladium oxide powders. J. Catal. 2007, 245, 35–44. (86) Wang, Z. Comprehensive Organic Name Reactions and Reagents; John Wiley & Sons, Inc., 2010. (87) Albers, P. W.; Möbus, K.; Wieland, S. D.; Parker, S. F. The fine structure of Pearlman’s catalyst. Phys. Chem. Chem. Phys. 2015, 17, 5274–5278. (88) Wellendorff, J.; Lundgaard, K. T.; Møgelhøj, A.; Petzold, V.; Landis, D. D.; Nørskov, J. K.; Bligaard, T.; Jacobsen, K. W. Density functionals for surface science:

34


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Exchange-correlation model development with Bayesian error estimation. Phys. Rev. B 2012, 85, 235149. ˇ aˇc, J.; Elstner, M. Parameterization of the (89) Kubillus, M.; Kubaˇr, T.; Gaus, M.; Rez´ DFTB3 Method for Br, Ca, Cl, F, I, K, and Na in Organic and Biological Systems. J. Chem. Theory Comput. 2015, 11, 332–342. (90) Grimme, S.; Bannwarth, C.; Shushkov, P. A Robust and Accurate Tight-Binding Quantum Chemical Method for Structures, Vibrational Frequencies, and Noncovalent Interactions of Large Molecular Systems Parametrized for All spd-Block Elements (Z = 1-86). J. Chem. Theory Comput. 2017, 13, 1989–2009. (91) Haas, H.; Wang, C. Z.; Fähnle, M.; Elsässer, C.; Ho, K. M. Environment-dependent tight-binding model for molybdenum. Phys. Rev. B 1998, 57, 1461–1470. (92) Goldman, N.; Goverapet Srinivasan, S.; Hamel, S.; Fried, L. E.; Gaus, M.; Elstner, M. Determination of a Density Functional Tight Binding Model with an Extended Basis Set and Three-Body Repulsion for Carbon Under Extreme Pressures and Temperatures. J. Phys. Chem. C 2013, 117, 7885–7894. (93) Gaus, M.; Cui, Q.; Elstner, M. DFTB3: Extension of the Self-Consistent-Charge Density-Functional Tight-Binding Method (SCC-DFTB). J. Chem. Theory Comput. 2011, 7, 931–948.

35



Graphical TOC Entry

36


Page 36 of 36

Tight-binding approximation enhanced global optimization

Recommend Documents