First-principles calculated decomposition pathways for LiBH4 nanoclusters

We analyze thermodynamic stability and decomposition pathways of LiBH4 nanoclusters using grand-canonical free-energy minimization based on total energies and vibrational frequencies obtained from density-functional theory (DFT) calculations. We consider (LiBH4)n nanoclusters with n = 2 to 12 as reactants, while the possible products include (Li)n, (B)n, (LiB)n, (LiH)n, and Li2BnHn; off-stoichiometric LinBnHm (m ≤ 4n) clusters were considered for n = 2, 3, and 6. Cluster ground-state configurations have been predicted using prototype electrostatic ground-state (PEGS) and genetic algorithm (GA) based structural optimizations. Free-energy calculations show hydrogen release pathways markedly differ from those in bulk LiBH4. While experiments have found that the bulk material decomposes into LiH and B, with Li2B12H12 as a kinetically inhibited intermediate phase, (LiBH4)n nanoclusters with n ≤ 12 are predicted to decompose into mixed LinBn clusters via a series of intermediate clusters of LinBnHm (m ≤ 4n). The calculated pressure-composition isotherms and temperature-pressure isobars exhibit sloping plateaus due to finite size effects on reaction thermodynamics. Generally, decomposition temperatures of free-standing clusters are found to increase with decreasing cluster size due to thermodynamic destabilization of reaction products.

To predict decomposition paths and reaction thermodynamics, we first need to predict the ground state geometries of LiBH 4 nanoclusters and candidate reaction products. We considered (LiB) n , (LiH) n , Li n , B n and Li 2 B n H n as possible reaction products for all n. The method for predicting the geometries of these nanoclusters will be described below. Since small lithium and boron clusters have been investigated theoretically by several authors, we have recalculated the boron clusters B n (n = 2 to 12) based on the the structures reported in the literature. [1][2][3] The structures of the ionically bonded (LiBH 4 ) n and (LiH) n clusters were obtained via prototype electrostatic ground state search (PEGS), 4 while the metallic clusters (Li n and (LiB) n ) were obtained via genetic algorithm (GA). Since (LiH) n clusters were found to form cutouts of the rocksalt structure, possible cluster topologies were enumerated manually and relaxed using the conjugate-gradient method and forces obtained from firstprinciples DFT calculations. We adopted the conformations of Li 2 B n H n clusters with n ≥ 6 from a previous study. 5

PEGS for (LiBH 4 ) n clusters
The reactant clusters, (LiBH 4 ) n were generated using the PEGS method described in ref. 4 , which has been used successfully to find several ground-state structures, high-temperature polymorphs and metastable structures in ionic materials. [6][7][8][9] In PEGS, one minimizes the electrostatic energy of a system of positively and negatively charged ion complexes, with no a priori symmetry constraints. Since covalent and metallic interactions are not explicitly included, PEGS per se is not expected to provide accurate energetics, and first-principles DFT calculations are always used to calculate the energies of the candidate ground state structures predicted by PEGS. The Hamiltonian is given by where Q i are the ionic charges, r ij are the distances between ions, and a soft-sphere repulsion term is nonzero when the inter-ionic distance is smaller than the sum of the ionic In order to find the structures of (LiBH 4 ) n clusters, our implementation utilizes the uniform histogram Wang-Landau (WL) algorithm 10 to efficiently search the configuration space. Due to the fact that the WL method uses an estimate of the density of states to accept or reject trial moves, it can avoid being trapped in local energy minima and hence explore all possible configurations with energies below a given upper bound. A detailed implementation of PEGS using the WL algorithm is given in ref. 11 , which describes application of the method to obtain ionic (NaH) n , (AlH 3 ) n , and (NaAlH 4 ) n nanoclusters with n ≤ 8. Each WL run produces many topologically distinct candidate structures with low electrostatic energies; these structures are relaxed using DFT and the lowest-energy structure is selected as the ground state configuration.

Genetic algorithm for metallic clusters
Since the PEGS Hamiltonian is not expected to provide a reasonable description of the energetics of metallic clusters, a different approach was adopted for Li n and Li n B n clusters, where we resorted to a genetic algorithm. Our genetic algorithm has been successfully used in global structural optimizations for C 60 , 12 medium-sized silicon clusters, 13 metal clusters, 14 and the reconstructions of crystal surfaces. [15][16][17][18] In the present study, our implementation is similar to our previous study for pure Al clusters. 14 Three types of operations including mating, mirror, and mutation, have been used to generate new structures in our GA optimization. In the mating operation, two parent structures A and B are randomly picked from the population pool and randomly rotated about their centers of mass which are set to be the origin. Then, a child structure is generated by mating the lower-half of parent A with the upper-half of parent B. If the number of atoms in the new structure does not match with that in the parent, parents A and B are shifted along the z -axis in order to achieve the correct match. In the mirror operation, a parent structure is picked from the population pool and rotated randomly about its center of mass which is set to be the origin.
This parent is shifted along the z-axis in order to balance the number of atoms in both halves. An offspring is created simply by duplicating the lower-half into the upper-half.
The new upper-half is shifted along the z-axis to reach an appropriate distance between the two halves. Finally, in the third operation, mutation, a parent structure is picked from the population pool and the outmost atoms of this cluster are randomly removed and then recapped on the surface of the cluster. The mutation operation was used in our previous study 19,20 and is useful for studying doping in the cluster.
In order to maintain diversity, it is important to apply a structural comparison algorithm when updating the population pool. Our structural comparison algorithm is based on a contact matrix 21 rather than the computationally demanding method used previously in ref. 14 . Note that similar but different isomers may have the same contact matrix which will be excluded from the generation pool. Thus, we imposed a condition that if the child structure has an energy difference greater than 30 meV/atom compared to a structure with the same contact matrix, we still insert the new child structure into the populations provided that the child structure has a sufficiently low energy to be included in the pool.
The genetic algorithm optimizations in the present study were performed as follows. In generation zero, a population pool of p = 50 structures for a given cluster size n were generated by putting n atoms randomly in a cubic simulation box of moderate side length (see details in the next section). The structures were relaxed to their local minima using a first-principles DFT calculation with the steepest descent method and a lower cut-off energy. In each subsequent generation, one new structure was generated from three types of operations (i.e., mating, mutation, and rotation as described above) with a given probability distribution: 80% of chance for mating, 10% of chance for mutation, and 10% for rotation operations. These new structures were relaxed to a local minimum configuration by a firstprinciples DFT calculation. The population pool was updated if the newly relaxed structures were not identical to any isomers inside the population pool and have lower energies. Genetic algorithm optimizations were produced for up to 1000 generations.
At the end of the global minimum structure search, as many as p = 50 structures remain as possible candidates. All candidates are further optimized using accurate first-principles DFT calculations with a high cut-off energy and a larger simulation box (see below).

First-principles DFT calculations
First-principles DFT calculations were carried out within the generalized gradient approximation (GGA) to the exchange-correlation functional 22-24 using projector-augmentedwave pseudopotentials 25 as implemented in the Vienna Ab-initio Simulation Package. 26,27 Two different settings were used for different accuracies. When the VASP is coupled with our genetic algorithm, the kinetic energy cut-off is set to be 239 eV and non-spin-polarized calculation is performed. The length of the simulation box is set to the diameter of a cluster plus 8Å with a smaller box size. The structural optimization is done with the steepest descent method and without symmetry until the forces on the atoms are less than 0.01 eV/Å.
At the end of a GA run, the kinetic energy cut-off is increased to 398.3 eV and spinpolarized calculation is performed. For clusters containing less than 8 formula units, the length of simulation cell is set to 15Å. For n > 9, the supercell length is 20Å. After the structures are determined, we use a higher kinetic energy cut-off equal to 875 eV, and the structural optimization is done with the conjugate gradient algorithm and without symmetry until the forces on the atoms are less than 0.001 eV/Å. This setting with a high accuracy was used for all the clusters in this study (shown in Fig. 1 from main text) including those obtained from PEGS and from the literature. After precise electronic structure calculations, we performed vibrational frequency calculations using the PHONOPY package. 28 We ensure the mechanical stability of all the clusters by checking the phonon frequencies for imaginary values. Temperature dependent vibrational free energies are obtained by summing contributions from all phonon modes, excluding rigid translations and rotations.
We further examined whether the candidate structures from the GA coupled with VASP with a lower accuracy setting underwent substantial relaxations upon first-principles optimizations. We found that the calculations using both low and high accurate settings retain the same energy trends and atomic structures. In order to ensure that the ground state structures for large clusters are not missed, we have chosen, for each cluster size, all 50 candidate structures in the final GA pool as the initial structures for our first-principles optimization using a higher accuracy.
ensuring that the total number of all non-hydrogen atoms is fixed. Feasible solutions to Eq. S1 form a finite set, which can be enumerated using integer programming techniques (i.e., using Wolfram Mathematica's built-in function FindInstance). The probability of finding any given configuration f is then given by where β = 1/k B T and F i (T ) is the free energy of cluster i, including vibrational contributions (pressure dependence of cluster free energies can be neglected). The number of hydrogen atoms in the cluster i is given by n H i , and the chemical potential of gaseous hydrogen as a function of pressure and temperature is given by µ H (p, T ). The partition function Z(p, T ) is calculated by summing over all configurations f that are consistent with the constraints given by Eq. S1: