Molecular simulations of crystal growth: From understanding to tailoring

ARTICLE IN PRESS

Molecular simulations of crystal growth: From understanding to tailoring Patrick Duchstein, Philipp Ectors, Dirk Zahn* Lehrstuhl f€ ur Theoretische Chemie/Computer Chemie Centrum, Friedrich-Alexander Universit€at Erlangen-N€ urnberg, Erlangen, Germany *Corresponding author: e-mail address: [email protected]

Contents 1. Introduction 2. Fundamental theory of crystal nucleation and growth 2.1 Classical nucleation theory and the role of crystal shape 2.2 Step nucleation and crystal growth 2.3 Wulff construction of crystallite shape 3. Explicit molecular simulations of crystal growth 3.1 Straight-forward molecular dynamics simulations of crystal nucleation and growth 3.2 Kinetic Monte-Carlo simulations of crystal nucleation and growth 3.3 Kawska-Zahn simulations of crystal nucleation and growth 4. Conclusions and outlook Acknowledgments References

2 3 3 5 7 13 13 16 17 21 22 22

Abstract Modern approaches to molecular simulations allow an increasingly accurate account of modeling crystal growth. Starting with easy-to-use modeling tools of limited predictive power, we review the recent progress in the field. The complexity of crystal growth calls for scale-bridging methods to provide high accuracy of describing atomic interactions where crucially needed, while effectively crossing time and length scales to ensure convergence from the viewpoint of statistical significance. By picking landmark case studies rather than claiming completeness, we illustrate the current forefront of modeling crystal growth. This includes shape prediction from connecting molecular scale simulations to μm sized crystallite models, the role of ripening reactions and of surfactant molecules. While there is still much room for improving computational efficiency and userfriendliness of the methods, the in-depth mechanistic understanding at reach makes molecular simulations an increasingly attractive tool for tailoring crystal growth.

Advances in Inorganic Chemistry ISSN 0898-8838 https://doi.org/10.1016/bs.adioch.2018.11.004

#

2019 Elsevier Inc. All rights reserved.

1

ARTICLE IN PRESS 2

Patrick Duchstein et al.

1. Introduction Crystals have always fascinated for their (apparent) purity and the puzzling manifold of shapes with well-defined geometry. Indeed, the presence of sharp faces and symmetric habit already hinted at mathematical principles of ordering, long before the concept of periodically arranged unit cells was demonstrated. With the beginning of modern crystallography—which we date back to Max von Laue’s seminal work on X-ray diffraction in 19121 and the development of Bragg’s law2—also approaches to understanding the shape of crystals improved swiftly. Crystal surfaces were then considered as cuttings of the crystal structure, in the simplest case exposing well-defined crystallographic planes as crystal faces. While we by now know that this may occur to be a drastic simplification, this picture is still used—at least as a starting point—for today’s concepts of constructing polyhedral models to describe crystallite shape. Provided that we know the relative preferences of the different crystal growth fronts, a fast and well-established route to geometric models may be derived from Wulff’s construction.3 The key question, however, remains how to assess this preference. To our knowledge, the first approach based on information of the crystal structure was taken by associating the different spacing of differently oriented lattice planes to their stability when exposed as crystal surface.4 While intuitively reflecting different trends in layerby-layer packing, molecular simulations offer the calculation of association energy per surface area, hence lifting the surface layer stability concept to a more quantitative level. In the past decade, molecular dynamics simulations moreover paved the way to assess growth mechanisms and to calculate the speed of different growth fronts of forming crystals. Wulff’s construction based on growth front dynamics, nowadays represents the most accurate, but also a computationally challenging approach to predicting crystal shape. The aim of the present review is to describe the fundamental theory of crystal growth and to summarize the molecular models in use. Depending on the simulation approach, a careful choice between quantum chemical descriptions, molecular mechanics (force-fields) and combined models have to be made. The choice is also system dependent, as the accuracy in calculating energy and forces, which is best for quantum methods, must be traded off for limited sampling of configurations. The simulation approach itself may vary from simple energy calculations, geometry optimization from energy

ARTICLE IN PRESS Molecular simulations of crystal growth

3

minimization to molecular dynamics and other kinetic approaches—which will all be discussed in this review. Finally, we illustrate the current perspectives of using molecular simulation to understand crystal growth from selected case studies.

2. Fundamental theory of crystal nucleation and growth 2.1 Classical nucleation theory and the role of crystal shape Classical nucleation theory provides a general framework to rationalizing nucleation, crystal growth/decay and crystal stability on the basis of thermodynamics. Despite a series of recent extensions and re-formulations, which shall be discussed later, the key ideas were already provided by Gibbs who contrasted the size-dependent interplay of two opposing driving forces to crystal growth or dissociation.5,6 The most obvious requirement to crystal growth is given by a gain in free energy ΔGbulk. Hence, crystal formation will only take place if ΔGbulk < 0 ΔGbulk ¼ ΔEpot  T  ΔS + p  ΔV

(1)

where the first term reflects the change in potential energy upon formation of a bulk crystal. In order to exist as a stable form, atomic/molecular packing in a crystal lattice must invoke a gain in potential energy plus volume work as compared to the solution or the melt (ΔEpot + p ΔV < 0). On the other hand, the entropy term TΔS is positive as entropy is highest for the disordered melt/solution and packing in periodic lattices implies a loss of entropy. Balancing ΔEpot + p ΔV ¼ TΔ S implies the transition point, i.e., provides a theoretical definition of the saturation concentration as a function of temperature T and pressure p—or of transition temperature and pressure in case of constant concentration. However, even when going beyond the crystallization conditions described by ΔGbulk < 0, a common observation is that crystal formation will not take place immediately, but is instead confined by the crossing of an activation barrier (or even several barriers). For the case of crystal formation from solution, the barrier may be intuitively be rationalized by the need to disconnect solute-solvent interactions. For crystallization from the melt, a somewhat less intuitive explanation is given by the loss of entropy by establishing local order in an elsewise disordered phase. The interface of the forming crystal to its parent phase (solution/melt) implies a considerable disfavoring to crystal growth. This, however, only applies to the early stage of

ARTICLE IN PRESS 4

Patrick Duchstein et al.

crystal formation, as only small crystals exhibit a significant surface to bulk ratio. In other terms, the thermodynamics of a forming crystal must be considered as size-dependent. ΔGcryst ¼ ΔGbulk + ΔGsurface

(2)

The most simple formulation of the interplay of surface/interface terms and bulk free energy gain assumes spherical crystals (with radius r) and thus denotes ΔGbulk and ΔGsurface as functions of the sphere volume and surface, i.e., ΔGbulk  r3 and ΔGsurface  r2, respectively. In what follows, we discourage this notation, and instead advocate a shape-independent formulation of classical nucleation theory based on the number of formulae units N comprised in the forming crystal. This is possible as any regular geometric shape implies that the surface area is proportional to the enclosed volume to the power of 2/3. The size-dependent formulation of Eq. (2) is then given by ΔGcryst ðN Þ ¼ Δμ  N + σ surface  N 2=3

(3)

Along those lines, the volume term in free energy is replaced by the N times the free energy change per formulae unit Δμ N. In other words, we directly denote the chemical potential difference of the compound in its crystalline form with respect to the solution/melt. Provided that ΔGbulk in Eq. (1) is taken per formulae units, it is identical to Δμ. While Δμ must be negative to allow crystal formation, the surface/interface term is positive and hence opposes crystal formation. Therein, Δσ surface is a materials constant that describes the overall crystal surface to solution/melt interactions at a given shape. In case that the shape is known, we may particularize the interface free energy as σ surface  N 2=3 ¼ σ ð100Þ  Að100Þ + σ ð110Þ  Að110Þ + ⋯

(4)

where A refers to all relevant faces observed and σ reflects the corresponding free energy density of forming such interfaces. A series of free energy profiles as function of the number of precipitated formulae units according to Eq. (3) is illustrated in Fig. 1. The different curves illustrate different chemical potential changes upon crystal formation—thus reflect different thermodynamic driving forces as resulting from different concentrations, temperature, etc. For Δμ < 0, crystal formation is exothermic, however requires the crossing of an energy barrier ΔG# ¼ G(Ncrit).

ARTICLE IN PRESS 5

Molecular simulations of crystal growth

G(N) = Δ μ⋅N+σ ⋅N2/3 Dm1 = 0

Dm2 < Dm1

ΔG

#

Ncrit

N Dm3 < Dm2

Fig. 1 Thermodynamic description of crystal nucleation and growth by classical nucleation theory. The free energy profiles are shown as functions of the number of precipitated formulae units N and illustrate the role of the chemical potential difference Δμ between the crystalline phase and the melt/dissolved state, while keeping the prefactor of the surface term σ constant. The location of the energy barriers and thus the size of critical nuclei Ncrit ¼ (⅔ σ/Δμ)3 depend on inverse of the cube of the degree of thermodynamic driving Δμ and thus approaches infinity at vanishing Δ μ (as for example observed for precipitation from solution until saturation concentrations are reached).

After exceeding critical size Ncrit, nuclei grow increasingly exothermic. For N ! ∞ the free energy per formulae unit converges to lim(N ! ∞): G(N)/N ¼ Δμ. Classical nucleation theory hence provides a thermodynamic mainframe to rationalize the existence of a nucleation barrier, but may also be extended to multi-step nucleation routes and Ostwald ripening of competing nuclei.7 The concept of relating interface and bulk terms to free energy may furthermore be transferred to describe crystal growth by step formation and propagation as discussed in the following.

2.2 Step nucleation and crystal growth While the most commonly used model of a crystal surface is that of a sharp plane, the description of crystal growth calls for more detailed models. Indeed, the most common route of crystal growth is that of step propagation. Hence, a given imperfection of the previously mentioned plane acts as an attractor for associating further material to amend the imperfection and drive the crystal surface toward the perfect plane—which would reflect the energetically favored state. The imperfections, in turn, are promoted by entropy

ARTICLE IN PRESS 6

Patrick Duchstein et al.

Fig. 2 Illustration of crystal growth as a cascade of building block association events. The most recently associated block is highlighted in blue for (i) spontaneous step formation, and (ii) propagation of the growth front. The kinetics of step nucleation and step growth are triggered by different changes in free energy Δ G(step nucleation) 6¼ ΔG (step propagation). The completion of layer deposition (iii) may hence reflect a particularly stable intermediate.

or the kinetics of precipitation process. In a simplistic picture, akin to the previously discussed rationalization of nucleation, crystal growth can be described as: (i) Spontaneous formation of an “island,” possibly a small patch of a single adlayer on an elsewise planar crystal surface (Fig. 2, left). Such asperity exposes considerable interface domains and its formation energy per formulae unit is much less as compared to the energy gain in forming the bulk crystal. Under mild crystallization conditions (moderate supersaturation, etc.) the free energy change in forming the bulk crystal Δμ(bulk) is slightly negative, whereas the free energy change for forming an adlayer island is positive, ΔG(step nucleation) > 0. Such situation may be considered as the nucleation of an adlayer in full analogy to the homogenous nucleation discussed in Section 2.1; free energy as a function of the number of precipitated formulae units would hence initially increase, as additional interface zones need to be formed. (ii) Once such an asperity, however, exceeds a critical size, the association of further solutes is preferred over asperity dissociation. Further growth is then directed by propagation of the step over the whole basal layer in order to “complete” the adlayer. The kinetics of this propagation process are typically much faster than those of step nucleation. As a rationalization of this, one may consider step propagation as a re-shaping of the crystal-solvent (or melt) interface, without actually increasing the overall size of unfavorable zones (Fig. 2, center). Indeed, fully canceling out interface contributions would imply that the free energy change of step propagation becomes equivalent to that of forming the bulk, suggesting ΔG(step propagation)  Δμ(bulk).

ARTICLE IN PRESS Molecular simulations of crystal growth

7

(iii) Unlike classical nucleation theory for homogenous nucleation, crystal growth by step propagation, however, stops when a full adlayer is deposited on a given crystal face. Once this happens, further growth requires the nucleation of another asperity as denoted in (i). This picture applies to nucleation-driven mechanisms of adlayer deposition and there are a number of cases following different routes. For example, crystal growth of CaCO3 in the calcite form is known as propagation of screw-dislocations—which makes step nucleation less important (or fully irrelevant) as the step permanently persists during crystal growth.8 Moreover, crystal growth may fully ignore the layer-by-layer concept and instead associate incoming solutes wherever they approach the precipitate. For CaCO3, the precipitation of the amorphous form represents an obvious example for this. Moreover, multi-step processes like deposition of weakly ordered adlayers, followed by structural transitions such as the amorphous ! aragonite and aragonite ! calcite transformations of CaCO3, are elusive to the simple scheme discussed above. Similar to structural transformation steps, secondary reorganization processes may also result from ripening reactions (like bicarbonate ! carbonate transformation by proton transfer reactions). Reliable modeling of crystal growth mechanisms thus requires considerable insight from the experiment or explicit molecular dynamics simulations as discussed in Section 3. Nevertheless, the step nucleation and propagation model provides a valuable concept for both general understanding and numerical modeling. From a general viewpoint, it helps to rationalize the kinetics of crystal growth qualitatively. Using mechanistic insights and formation energies/rates from detailed molecular dynamics simulations, the step nucleation and propagation model may furthermore form the basis of bridging time and length scales. Examples for quantitative modeling of crystal growth up to the mm scale by kinetic Monte-Carlo simulations will be provided in Section 3.

2.3 Wulff construction of crystallite shape The foundations of assessing the equilibrium geometry of a crystallite both in terms of shape and inner structure were provided by Gibbs concept of free energy. Assuming full relaxation, Gibbs suggested that the shape of a crystallite should be such that surface energy is minimal.5,6 This led to the development of a polyhedra-based approach to create geometric models of crystallites by Wulff.3 While the “Wulff construction” method dates back

ARTICLE IN PRESS 8

Patrick Duchstein et al.

to the beginning of the 20th century, it is still widely used, and only very recent developments in molecular dynamics and kinetic Monte-Carlo simulations paved the way to more detailed assessment of crystal growth and the prediction of shape from modeling. The latter types of simulation approaches as discussed in Section 3, however, not necessarily entirely rule out Wulff’s construction, but may also be taken as a new source of profound input to the Wulff construction method, or variations thereof. Consequently, in what follows we discuss both the original (and most common) implementation of Wulff’s construction and variations thereof. The most important model assumption of Wulff is to describe crystal surfaces as ideal planes, hence ignoring any imperfections such as asperities, steps, kinks, defects, etc. In principle, any type of surface planes may be modelled, however, it is very reasonable to limit the considerations according to suitable cuttings of the crystal lattice. There are a few rules of thumb of how such cuttings may look like. Based on empirical observations, crystal faces typically reflect planes of low Miller indices. Moreover, crystal cuttings that would lead to charged surfaces appear unfavorable from an electrostatic viewpoint. Instead of simply discarding such surfaces, we however advocate the use of zig-zag type cuttings in order to ensure charge-neutral, but slightly rough surfaces. While surely more realistic, such deviation from strictly planar cuts comes at the price of adding an extra degree of freedom to model creation. In case that suitable cuts for a given (hkl) layer cannot be provided from intuition, the modeler may scan a series of alternatives and choose the one of lowest energy (i.e., crystal cuts which invoke smallest energy costs for taking the 3D-periodic model apart along the [hkl] direction). Moreover, the relevant crystal faces and their local structure may also be derived from explicitly studying atom-by-atom growth from molecular simulation as discussed in Section 3. To predict the shape of a crystallite formed from the vapor or solution, Wulff’s construction uses the specific surface σ surf(hkl) or interface σ int(hkl) tension of the corresponding face. The latter are defined as the energy of forming the surface or interface, respectively, per surface area. Planes of high energy are typically ruled out and the construction of a polyhedron is then based on a small set of (hlk) faces of low surface/interface tension. For a given polyhedron volume, the ratio of the different surface areas is determined by Gibbs concept of lowest total surface energy. In other terms, the area A(hlk) of an (hlk) face is expected to be inversely proportional to its surface tension σ surf(hkl) (or σ int(hkl) for interfaces to a solution). To construct such geometries, Wulff defined a distance constant d to scale the overall volume of the

ARTICLE IN PRESS Molecular simulations of crystal growth

9

polyhedron, and an orientation-dependent distance constant d(hkl) ¼ d  σ surf(hkl).3 The latter is suggested as the center-of-mass to (hlk) surface distance. Hence, the extension of the polyhedron along the [hlk] direction (face-to-face distance) is given by

 dðhkl Þ + d hkl ¼ d  σ surf ðhkl Þ + d  σ surf hkl

(5)

which simplifies to 2  d  σ surf(hkl) for symmetric crystal faces. Mathematical proofs of the validity of this approach date back to von Laue,1 but extensions and generalization work are continuing up to present time.9 For a recent review of applying Wulff’s construction, we direct the reader to the more extensive article of Barmparis et al. who also discussed the modeling of nanoparticles and surfactants using extensions of Wulff’s approach.9 While there are numerous successful applications of Wulff’s construction, there are also cases of utter failure. Therefore, we feel it is crucial to discuss the limitations of this in principle quite powerful method. An obvious limitation to routine applications is provoked by the earlier mentioned possible problems in creating realistic crystal faces. To illustrate this effect, we used different cuttings of the wurtzite lattice of CdSe and different techniques of estimating the aspect ratio of nanorods. CdSe reflects a prominent case study example, as control of size and shape is of imminent importance for tuning the optical properties.10–13 Moreover, there are tailor-made molecular interaction models available that were proven to reproduce the different polymorphic forms.14 A CdSe nanorod model as employed in our calculations is illustrated in Fig. 3. According to Eqs. (3) and (4) the total energy of the nanorod comprising N formulae units is given by: E ðN Þ ¼ Δμ  N + σ f0001g  Af0001g + σ f1010g  Af1010g 2 pffiffiffi 2 3a Af0001g ¼ 2  3 Af1010g ¼ 6  a  c

(6)

where a and c are the dimensions of the hexagonal rod, while A{1010} and A{0001} denote the total surface area of all lateral and top/bottom faces, respectively. In the gas phase, Δμ ¼ 11.48 eV is taken as the lattice energy per formulae unit of the wurtzite crystal. For precipitation from solution, this energy must be reduced by the free energy of solvation per formulae unit. Using a small set of nanorods of different sizes, we can identify the

ARTICLE IN PRESS 10

Patrick Duchstein et al.

Fig. 3 Center: illustration of the 7  7  10 hexagonal nanorod model of CdSe in water. Left and right: 10  10  10 and 7  7  20 models as used for comparing total energies and identifying σ {0001} and σ {1010} from Eq. (6). Colors: Cs (yellow), Se (brown), water (green).

different surface/interface energy densities σ {10–10} and σ {0001} from calculating E(N) (resp. G(N) for crystal growth from solution). Based on 5  5  10, 7  7  5, 7  7  10, 7  7  20 and 10  10  10 hexagonal a  a  c models and the force-field of Rabani,14 we provided estimates of the nanorod shape as preferred in the gas phase. The nanorods were then immersed into solution (mimicked by 3D-periodic bulk models of ten thousands of ethanol, DMSO and water molecules—hence covering a range of differently polar solvents with εr ¼ 24.6, 46.7 and 78.4, respectively). To illustrate the role of suitable cuttings of the crystal structure, the (0001)/(0001) faces were taken as (i) planar cuts, thus leading to polar nanorod models, and (ii) zig-zag type cuts with 50% occupancy

ARTICLE IN PRESS 11

Molecular simulations of crystal growth

Table 1 Calculated surface/interface energy densities for the two types of faces of a CdSe nanorod in various solvents and shape prediction based on Wulff construction. Model σ {10210}/eV Å22 σ {0001}/eV Å22 c:a ratio

Polar cut (vacuum) Zig-zag (vacuum)

0.11 0.13

0.40 0.43

3.6:1 3.3:1

Polar cut (ethanol) Zig-zag (ethanol)

(0.19) (0.19)

(0.10) (0.10)

1.9:1 1.7:1

Polar cut (DMSO) Zig-zag (DMSO)

(0.06) (0.05)

(0.06) (0.07)

3.6:1 2.4:1

Polar cut (water) Zig-zag (water)

(0.24) (0.24)

(0.18) (0.17)

2.7:1 2.3:1

For the polar models, the average of the two different interface energy densities according to Cd/Se-terminated (0001) and (000–1) surfaces is shown. In all solvent simulations, the embedding liquid is chosen as pure ethanol/DMSO and water at 1 atm. and 300 K, respectively. Concentration-dependent chemical potential differences are ignored by only sampling averages of the solvent-surface and solvent-solvent interactions (σ values in brackets). The resulting interface terms are then used for re-scaling the nanorod shape predicted for the gas phase. On this basis, solvent effects are estimated for the approximation of infinitely diluted solutions.

of the outermost surface layer implemented by randomly removing the corresponding atoms. The resulting surface/interface energy densities obtained for vacuum, ethanol, DMSO and water, along with estimates of the c:a aspect ratio as obtained from Wulff’s construction are denoted in Table 1. It is educative to compare these findings to the attachment energy approach.15 For this, we first attempted shape prediction in the gas phase based on studying the association of a formulae unit to the different faces of the nanorod. To assess the role of different solvent environment, we then calculated the association energy for attaching a single solvent molecule to the nanorod surface. This provides hints toward shape changes as compared to the gas phase structures calculated before. We note that the attachment of single solvent molecules ignores cooperative phenomena of crystal-solvent interfaces and uses the corresponding energies for qualitative shape estimates, only (Table 2). Indeed, comparing the predicted nanorod shapes from attachment energy to the c:a ratios obtained from interface energy density, we find quite remarkable discrepancies. Based on the established experimental evidence, the shape of wurtzitetype CdSe ranges from nearly spherical particles to elongated rods.10–13 As a simple classification, Alivisatos and coworkers suggested the shape as

ARTICLE IN PRESS 12

Patrick Duchstein et al.

Table 2 Calculated attachment energies for the two types of faces of a CdSe nanorod. Model Eattach{10 210}/eV Eattach{0001}/eV c:a ratio

Polar cut—CdSe pair Zig-zag—CdSe pair

3.4  3.4

4.2/4.3  35.1

1.2:1 10:1

Polar cut—ethanol molecule Zig-zag—ethanol molecule

0.6  0.6

0.9/0.7  1.2

Elongate Elongate

Polar cut—DMSO molecule Zig-zag—DMSO molecule

0.7  0.7

1.9/1.3  1.4

Elongate Elongate

Polar cut—water molecule Zig-zag—water molecule

1.2  1.2

0.8/0.7  1.5

Flatten Elongate

The gas phase study refers to association of a CdSe ion pair and is used for quantitative estimates of the a:c aspect ratio. Comparing the total energy of the polar and unpolar species, we note that the zig-zag shaped surfaces are drastically disfavored. This leads to unrealistic estimates of the attachment energy for adding the CdSe pair. For the polar models, two attachment energy values are provided to discriminate the Cd/Se-terminated (0001) and (000 1) surfaces. Moreover, a qualitative estimate of the effect expected from solvent (or surfactant) molecules in scaling the a:c ratio is deduced from single molecule attachment energies.

hexagonal rods in each case, however, with aspect ratios (dimensions along c: a axes) varying from 1 to 10, respectively.16 The underlying syntheses routes involve precursor molecules, injection into anti-solvents, reactions at elevated temperature and surfactant molecules.16 Hence, the diversity of mechanisms is immense and clearly exceeds the simplistic models applied for the beforehand energy calculations. We argue that this is actually the mayor concern about estimating crystal shape from Wulff’s construction. In case that orientation-dependent growth rates were available, the construction of realistic polyhedra would be rather straight-forward. While we discuss at least some degree of extending molecular simulations to realistic setups of crystal growth in the next section, it is still educative to interpret the simple models from a qualitative viewpoint. Comparing the attachment-energy approach to the Wulff construction from surface/interface energy densities, we indeed investigate rather different systems as prototype models to describing crystal growth. The attachmentenergy of a formulae unit to a given surface plane is often interpreted as an easy-to-get estimate of the overall energy density of a given crystal face. While this is a rather rough estimate, we however suggest that the underlying models can also (and likely better) be interpreted as proxies to step nucleation as illustrated in Fig. 2(i). As a consequence, the attachment-energy approach might actually become useful for predicting the shape of crystals which growth mechanism is triggered by the nucleation of a deposit island, followed

ARTICLE IN PRESS Molecular simulations of crystal growth

13

by comparably fast step propagation. In turn, the Wulff construction using surface/interface energy densities appears more suited for syntheses that provide high temperature or rapid (and thus indifferent) reactions of precursors. In such cases, the above kinetic argument imposed for the formation of deposit islands on planar surfaces no longer holds and step nucleation and step propagation steps would occur at nearly equal pace. When being comparably free to move along various modes of growth, the forming crystals should evolve toward the thermodynamically favored shape as suggested from balancing the surface/interface energy densities.

3. Explicit molecular simulations of crystal growth 3.1 Straight-forward molecular dynamics simulations of crystal nucleation and growth In the past decades, the progress in computational resources paved the way to direct molecular simulation of crystal nucleation and growth from solution. Initial studies were related to strongly simplified systems mimicking solutes and solvents by spherical particles.17–22 Such generic models allow insights into the fundamental physics of solute association, solvent dissociation and solute organization. For example, Anwar and coworkers used two types of soft-sphere particles to mimic a solution in contact with a crystal surface of the solute species to explore crystal growth from direct molecular dynamics simulation.23 On this basis, they demonstrated different modes of crystal growth as functions of supersaturation, classified as low, moderate and high. Low supersaturation leads to surface growth by nucleation of (multiple) steps and propagation as a rough growth front. High supersaturation instead showed rapid nucleation in the solution leading to initially dispersed crystallites. In-between, a “mixture” of both processes was observed. Here, solute clusters are formed in the solution, while heterogeneous nucleation is identified as the surface-induced transformation of weakly ordered aggregates into crystalline nuclei.23 The key motivation of such generic models is to mimic crystal/solution systems at moderate computational costs. Simple and short-ranged interaction potentials suffice to describe soft-sphere models, and hard-spheres allow even simpler treatment. Going from generic models to realistic setups of crystal nucleation and growth evokes two extra challenges to molecular simulation. First, appropriate interaction models are required to describe solutesolute, solute-solvent and solvent-solvent interactions both accurate and computationally efficient. This is typically implemented by molecular

ARTICLE IN PRESS 14

Patrick Duchstein et al.

mechanics, rather than quantum chemistry. While force-fields are usually parameterized to quantum characterization of small subsystems, their accuracy in describing solutions and crystals still requires careful benchmarking. Second, most solute-solvent systems of interest actually refer to drastically lower solute concentrations as compared to the saturated soft-sphere solutions used for the generic models simulations. While direct molecular dynamics simulations of crystal nucleation and growth can cope with solute to solvent ratios of 1:1 to 1:100, experimental setups often imply ratios of 1:104 to 1:106. For such systems of (initially) very dispersed solutes, direct molecular simulations are impractical and we suggest special techniques to implicitly model diffusion processes as discussed in Sections 3.2 and 3.3. Given the above limitations, it is not surprising that the mayor body of direct molecular dynamics simulation studies of crystal nucleation and growth from solution so far focused on NaCl/water24–30 and urea/water.31–35 Both solute species are highly solvable in water and realistic model setups of low supersaturation require only around 10 water molecules per solute. Another feature of these two systems is the availability of robust interaction potentials which suitability have been shown by confirmation of crystal structures (including polymorphs),36,37 melting point and solubility estimates38,39 in line with the experiment. Molecular simulations have contributed significantly to our understanding of NaCl from salt water solutions. The current mechanistic model as illustrated in Fig. 4 involves several stages, and nicely demonstrates the complexity of “real” nucleation and growth as compared to the theoretical concepts discussed in Section 2. An increasingly prominent issue is that of so-called two-step nucleation.7,40 Starting from slightly supersaturated salt water solution models, Chakraborty and Patey studied the nucleation and growth mechanisms of NaCl precipitates comprising up to 1000 ions.30 In agreement with earlier work on the initial stages of ion aggregation,27 the precipitation of NaCl from salt water is expected to be initiated from the association of Na+ and Cl ions in partially ordered, not-yet stoichiometric aggregates that comprise a significant amount of water molecules. Instead of directly forming solute aggregates arranged according to the structure of the final crystal, the system hence takes a “short-cut” to precipitation in terms of forming a less/differently ordered intermediates. The observed partial de-solvation along with partial ordering of solutes before crystal formation may intuitively be rationalized from Ostwald’s concept based on structural similarity41—or from a kinetic viewpoint based on the route of lowest transformation barrier. However, there is increasing

ARTICLE IN PRESS Molecular simulations of crystal growth

15

Fig. 4 Current understanding of NaCl nucleation and growth from salt water. The initial aggregate (A) and the outer boundaries (B) of infant aggregates are only partially ordered and retain some content of solvent molecules (green). Simulations however suggest local motifs of the later NaCl crystal to act as centers of stability (yellow) which determine the fate of the aggregate. Successful precipitates (C) show increasing ordering with increasing size. The snapshots in the upper panel are adapted with permission from Zahn, D. Phys. Rev. Lett. 2004, 92, 040801/1–040801/4; Chakraborty, D.; Patey, G.N. J. Phys. Chem. Lett. 2013, 4, 573–578 and refer to 20 and 500 formulae units of the solute, respectively.

evidence that the apparent cascade in free energy of a “solution ! loose agglomeration ! desolvation and ordering to crystals” process might instead reflect the minimum energy path as a function of precipitate size.7 The stability of small aggregates is significantly influenced by the interface free energy term Δ Gsurface (Eq. 2). Favorable interface structures may hence compensate the somewhat less favorable bulk energy ΔGbulk of the intermediate forms, provided that the aggregates are sufficiently small. For larger particles, the surface/interface terms, however, become less relevant and the overall thermodynamic stability is dominated by ΔGbulk. So far, even the largest simulation study on NaCl crystallization from salt water currently available did not get beyond the precipitation of 500 formulae units of NaCl as observed within 150 ns.30 The inset in Fig. 4 illustrates a final configuration from this work. It shows a single-crystalline domain with much faceting and the favoring of {100} type surfaces may be seen only vaguely in ref. 30. To predict the habit of meso- to macroscopic

ARTICLE IN PRESS 16

Patrick Duchstein et al.

crystals of NaCl, the common approach is still to anticipate that on the millisecond to hour time scales, there will be structural reorganization and Ostwald ripening much beyond what could have been observed from direct molecular dynamics simulations so far. Using ideal cuttings of the rocksalt crystal, Smith demonstrated the different solvent-surface interplay for NaCl-water and NaCl-urea interfaces.29 On this basis, the preference of {100} and {111} faces was rationalized for water and urea solutions, respectively, hence allowing to explain the cubic versus octahedral shape formation seen in the experiments.42

3.2 Kinetic Monte-Carlo simulations of crystal nucleation and growth As a stand-alone, the beforehand discussed direct molecular dynamics simulations are not yet suited for studying the whole process of forming crystals with well-defined faces from solution. With current computational resources, even the particularly efficient NaCl-water and urea-water systems are too demanding for such brute-force endeavor. However, lessons from these studies, such as the average rate of solute uptake on different crystal faces (prepared as slab models from cutting ideal crystal lattices), may be used to boost simulations considerably. The key to this is to simplify the complexity of possible association steps to a representative set of prototype growth steps as illustrated in Fig. 2. Piana and Gale outlined kinetic Monte-Carlo simulations of urea crystal growth from solution indeed bridging the atomic scale of urea/solvent solutions to meso-scale crystals.31–33 On this basis, the experimentally observed habit of μm-sized urea crystallites grown from water and from methanol could be reproduced using simulation models of equal size. The starting point to this was molecular dynamics simulations of four different types of urea slaps in contact with urea/water and urea/methanol solutions. On this basis, the association rates to planar cuts of [100], [110], [111] and [1 1 1] surfaces were characterized from direct molecular dynamics simulations. Moreover, Piana and Gale made use of the abundant interplay of association and dissociation events observed on the 50 ps scale for the chosen system. Based on the hundreds of ns simulation time, an extended characterization of the manifold of deposition $ dissociation steps going on was possible. Classification of these transformation steps and their rates was accomplished on the basis of nearest-neighbor and second-nearest neighbor coordination. This led to considering a total of 34 different processes that include local features of kinks, steps and islands—along with statistical

ARTICLE IN PRESS Molecular simulations of crystal growth

17

sampling of the respective association rates (cf. Fig. 2 which illustrates only two types of growth steps). Starting from seed crystallites of post-critical size, kinetic Monte Carlo simulations were then performed to reach much larger time and length scales. Instead of considering the full flexibility and dynamics of all molecules (as done in molecular dynamics simulations), crystal growth was reduced to the set of the beforehand characterized types of 34 growth events.32,33 To study time-resolved crystal growth, a time step of 50 ps was associated to each Monte-Carlo step. Growth steps were then implemented randomly using controlled probabilities. For the different types of growth events the individual occurrence is simply given by the corresponding association rate multiplied by the selected time step. This manner of controlling the occurrence of Monte-Carlo steps reflects a considerable speed-up over the conventional implementation of MonteCarlo moves based on first attempting a fully random move, and then accepting/rejecting according to a probability of p ¼ exp(ΔE/kBT). In the latter type of Monte-Carlo simulations, random moves often lead to unfavorable energy changes ΔE and thus a decline of the growth attempt. The approach of Piana and Gale not only benefits from larger acceptance rates, but also avoids the calculation of energy changes for individual Monte-Carlo steps at all.33

3.3 Kawska-Zahn simulations of crystal nucleation and growth Despite the impressive boost in time and length scales demonstrated for kinetic Monte-Carlo simulations, we caution that the underlying classification of possible growth steps and assessment of the respective rates remains a mayor issue. So far, we discussed systems of particularly high solute solubility for which interface models of solutions and crystal surfaces are comparably easy to handle in molecular simulations. However, contrasting NaCl precipitation from salt water to that of other compounds, drastically lower solute concentrations must be chosen to provide realistic accounts of the experiment.43 Coping with millions of solvent molecules per solute, as for example needed in CaF2 water solutions, causes direct molecular simulations impractical in terms of both time and length scales. For such systems the required information on crystal growth steps and rates can no longer be assessed from direct sampling statistics and we miss the prerequisite for the kinetic Monte-Carlo as discussed in the previous section. There are ongoing efforts aiming at employing enhanced sampling techniques to probe single solute association steps triggered by adding extra

ARTICLE IN PRESS 18

Patrick Duchstein et al.

forces to boost association/dissociation processes. A particularly illustrative study of this kind was recently presented by Gale and coworkers on the uptake of calcium carbonate to a surface terrace of calcite.44 While demonstrating insights into aspects of the growth mechanism (propagation of a screw dislocation), a full account of all relevant types of growth events akin to the beforehand discussed urea study appears as far from feasible. Slow association kinetics may stem from the energy barrier related to de-solvation, i.e., the re-arrangement of solute-solvent coordination upon solute-solute contact. However, at low concentration also the simple need of long-ranged diffusion before solutes actually reach a growing crystallite considerably affects the growth rates. Solute diffusion in bulk solution may, however, be rather easily be approximated: solutes will approach a forming crystal from random incoming direction at a rate determined by solution concentration and the diffusion constant. The Kawska-Zahn approach to crystal nucleation and growth makes use of this feature and hence implements Monte-Carlo steps as the approaching of a solute next to a forming aggregate or growing surface.43,45 As detailed knowledge on the modes of solute association and precipitate re-arrangement upon growth is typically not known, this process is studied from explicit molecular dynamics simulations. In other terms, Monte-Carlo steps and molecular dynamics are combined, in such way that Monte-Carlo steps are used for placing a solute near the precipitate, while the incorporation or rejection process is explored from unconstrained molecular dynamics simulations. This approach performs for a wide range of systems, provided that dissociation events play only a minor role to crystal growth (dissociation is considered during the molecular dynamics steps, but proper rate balance of association/dissociation events is not provided43,45). The technique may furthermore be extended to modeling proton transfer or redox-reaction steps as part of the relaxation process, by means of quantum/classical and semi-empirical approaches, respectively.46,47 Snapshots from ZnO nucleation and growth from ethanol solution are illustrated in Fig. 5. Using the Kawska-Zahn approach implies ion-by-ion association and relaxation after every growth step. Simulation studies dedicated to aggregate formation and the evolution of nuclei are hence typically limited to the assembly of a few hundreds to thousands of ions.48 On this basis, the key mechanisms of self-organization are assessable, including secondary transformations of forming nuclei (such as the ripening of Znx(OH)y aggregates to ZnO/Zn(OH)2 core-shell nanoparticles shown in Fig. 546). However, for the characterization of surface growth mechanisms, it is much more efficient to prepare different crystal-solvent interface models and

ARTICLE IN PRESS Molecular simulations of crystal growth

19

Fig. 5 Molecular simulation of zinc oxide nucleation and growth from ethanol solution. (A) The first stage of nucleation reflects the association of Zn2+ and OH ions from solution. Upon aggregation of a few tens of ions, OH  OH hydrogen bonds are observed. In case of dissimilar local ion arrangement, OH + OH ! O2 + H2O reactions become exothermic (Note the strong coordination by Zn2+ ions of the oxygen atom donating the proton and the comparably weaker coordination by the accepting hydroxide ion). Such ripening reactions lead to (B) drastic re-arrangements and the formation of a ZnO domains of wurtzite structure. (C) Similar to homogeneous nucleation, also the initial formation of an island on a ZnO surface is observed as the association of zinc and hydroxide ions. However, upon further ion deposition, ripening reactions lead to the growth of the ZnO phase while the crystal grows. Note that the initially planar surface evolves in favor of a ridged surface. The ridges were found to expose (1  100) and (0  110) faces, hence analogs of the {10–10} surface type, and comprise of Zn/OH ions at corners and edges, while less exposed surface patches are comprised of Zn/O ions. Pictures adapted with permission from references Kawska, A.; Duchstein, P.; Hochrein, O.; Zahn, D. Nanoletters 2008, 8, 2336–2340; Anwar, J.; Zahn, D. Angew. Chem., Int. Ed. 2011, 50, 1996–2013.

explore ion deposition on pre-modelled crystal faces. In full analogy to the direct crystal growth simulations discussed in Section 3.2, the starting systems need not be realistic; upon deposition of a few ion layers the system will evolve to preferred interface structures by its own (Fig. 5C49). On this basis, unprejudiced surface growth mechanisms may be studied from 10 nm scale atomistic models. Apart from the role of the solvent on stabilizing different surfaces, also the interaction with surfactants can be analyzed at the atomic level of detail. A particularly illustrative example is given by the association of surfactants to the (1010) growth front of zinc oxide

ARTICLE IN PRESS 20

Patrick Duchstein et al.

shown in Fig. 5. This was explored for acetate and citrate by means of picking a snapshot of the stepped surface and systematically probing surfactant docking to different sites. The roughness of realistic crystal growth fronts indeed implies considerable heterogeneity to binding the surfactants. Namely, the carboxyl-groups were found to favorably coordinate a single zinc ion. Such coordination strongly depends on the local configuration of the corresponding Zn2+ at the rough ZnO surface. The more exposed the zinc ion, the stronger becomes the stabilization by acetate surfactants. As a consequence, acetate selectively binds to (in this order of preference) islands, corners and edges of the stepped (10  10) surface (Fig. 6, upper panel).

Fig. 6 Association of acetate (upper panel) and citrate (lower panel) to the stepped (10  10) surface growth front illustrated in Fig. 5. Picking a snapshot of the rough surface, a series of docking attempts were performed to sample statistics of the local binding arrangements. From this, rather drastic preferences were identified. The example of citrate demonstrates the analysis of growth-control mechanisms by molecular recognition of growth steps. Pictures adapted with permission from Anwar, J.; Zahn, D. Angew. Chem., Int. Ed. 2011, 50, 1996–2013.

ARTICLE IN PRESS Molecular simulations of crystal growth

21

On the other hand, citrate exhibits three suitable binding partners to form salt bridges with zinc ions. Because of the molecular alignment of the citrate, it can however only establish a limited number of salt-bridges when binding to a flat surface patch. Still such association was found as exothermic by 2 eV. Citrate association to a surface step however allows binding arrangements that make use of all three carboxy-groups yielding an extra energy gain of around 4 eV.50 This demonstrates the pinning of growth fronts by citrate additives (Fig. 6, lower panel).

4. Conclusions and outlook Molecular simulations are advancing to an increasingly reliable and insightful tool of investigating crystal formation. This encompasses a series of current topics in solid state research. Syntheses planning may benefit from modeling-based search for possible crystal structures. Early studies in this field were dedicated to inorganic compounds and date back to the 1990s.51,52 Using inexpensive force-field model, Monte-Carlo or even fully random walkers were implemented to explore the manifold of solute arrangements within the solid phase. Ordered motifs were classified and transferred into periodic arrangements. The latter may then be subjected to high-level force-field or quantum mechanical calculations to provide an energy ranking of the suggested crystal structures. By now, computational resources advanced to a degree that crystal structure searches may be directly based on quantum calculations.53 Moreover, crystal structure prediction from theory was extended to molecular compounds54,55 and by now became an indispensable guide for the pharmaceutical industry to identify possible polymorphs relevant to drug patenting and drug formulation.56 A related field is that of studying crystal nucleation by molecular simulation.48 Atomic-scale simulations are of particular value to provide in-depth understanding of how a crystal forms. Thus, theory is not only used to suggest possible polymorphic forms of a material, but also helps in identifying possible synthetic routes to actually produce them. Such rational design of syntheses is, however, far from routine and there is still much room for advances—both in terms of efficiency and accuracy.48 Nevertheless, we feel that ongoing progress is promising and the importance of the topic justifies the wide attention from the molecular simulation community. Moreover, molecular simulation techniques initially developed for tackling crystal nucleation often also qualify for studying crystal growth. Indeed, the case studies described in Section 3 are all based on approaches that could be envisaged for exploring crystal nucleation.

ARTICLE IN PRESS 22

Patrick Duchstein et al.

From the viewpoint of theory, the approaches to crystal growth and prediction of crystallite shape experienced a paradigm change. The pioneering work of Gibbs,5,6 Wulff 3 and Donnay4 has been extended considerably. Today, Wulff’s construction of polyhedra mimicking crystal shape can be based on attachment-energies or interface energy densities from routine calculations.9 The case study of CdSe nano-rods however shows the imminent transferability problem when using simplistic models adjusted to equilibrium geometries for rationalizing actually rather complex non-equilibrium reactions.16 With new molecular simulation techniques at hand, we are gradually approaching realistic setups that at least describe key aspects of the manifold of nanoparticle syntheses routes.32–35,49 Such methods might once fully replace Wulff’s construction. However, an appealing challenge could also lie in using insights from detailed molecular simulation of solute-by-solute crystal growth to re-direct Wulff’s construction toward more realistic interface models and statistically more relevant sampling of the corresponding energy terms.

Acknowledgments The authors gratefully acknowledge the DFG for funding within the EXC315 initiative.

References 1. Friedrich, W.; Knipping, P.; von Laue, M. Sitzungsberichte der Mathematisch-Physikalischen Klasse der K€ oniglich-Bayerischen Akademie der Wissenschaften zu M€ unchen. K€ oniglich Bayerische Akademie der Wissenschaften, 1912;303. 2. Bragg, W. H.; Bragg, W. L. Proc. R. Soc. London, Ser. A 1913, 88, 428–438. 3. Wulff, G. Z. Krist. 1901, 34, 449–530. 4. Donnay, J. D. H.; Harker, D. Am. Mineral. 1937, 22, 463. 5. Gibbs, J. W. Trans. Conn. Acad. Art Sci. 1876, 3, 108–248. 6. Gibbs, J. W. Trans. Conn. Acad. Art Sci. 1876, 16, 343–524. 7. Zahn, D. ChemPhysChem 2015, 16, 2069–2075. 8. Sekkal, W.; Zaoui, A. Sci. Rep. 2013, 3, 1587. 9. Barmparis, G. D.; Lodziana, Z.; Lopez, N.; Remediakis, I. N. Beilstein J. Nanotechnol. 2015, 6, 361–368. 10. Heath, J. M. Ed.; Acc. Chem. Res. 1999, 32 (Nanoscale materials special issue), 387–387. 11. Alivisatos, A. P. Science 1996, 271, 933–937. 12. Lieber, C. M. Solid State Commun. 1998, 107, 607–616. 13. Smalley, R. E.; Yakobson, B. I. Solid State Commun. 1998, 107, 597–606. 14. Rabani, E. J. Chem. Phys. 2002, 116, 258. 15. Hartman, P.; Bennema, P. J. Cryst. Growth 1980, 49, 145–156. 16. Peng, X.; Manna, L.; Yang, W.; Wickham, J.; Scher, E.; Kadavanich, A.; Alivisatos, A. P. Nature 2000, 404, 59–61. 17. Lomakin, A.; Asherie, N.; Benedek, G. B. Proc. Natl. Acad. Sci. U.S.A. 2003, 100, 10254–10257. 18. Auer, S.; Frenkel, D. J. Chem. Phys. 2004, 120, 3015–3029. 19. van Duijneveldt, J. S.; Frenkel, D. J. Chem. Phys. 1992, 96, 4655–4668.

ARTICLE IN PRESS Molecular simulations of crystal growth

23

20. ten Wolde, P. R.; Frenkel, D. Science 1997, 277, 1975–1978. 21. Anwar, J.; Boateng, P. K. J. Am. Chem. Soc. 1998, 120, 9600–9604. 22. Moroni, D.; Wolde, P. R.; Bolhuis, P. G. Phys. Rev. Lett. 2005, 94, 235703/1–235703/4. 23. Anwar, J.; Khan, S.; Lindfors, L. Angew. Chem. Int. Ed. 2015, 54, 14681–14684. 24. Shinto, H.; Sakakibara, T.; Higashitani, K. J. Chem. Eng. Jpn. 1998, 31, 771–779. 25. Shinto, H.; Sakakibara, T.; Higashitani, K. J. Phys. Chem. B 1998, 102, 1974–1981. 26. Oyen, E.; Hentschke, R. Langmuir 2002, 18, 547–556. 27. Zahn, D. Phys. Rev. Lett. 2004, 92, 040801/1–040801/4. 28. Bahadur, R.; Russell, L. M.; Alavi, S.; Martin, S. T.; Buseck, P. R. J. Chem. Phys. 2006, 124, 154713. 29. Smith, P. E. Fluid Phase Equilib. 2010, 290, 36–42. 30. Chakraborty, D.; Patey, G. N. J. Phys. Chem. Lett. 2013, 4, 573–578. 31. Piana, S.; Gale, J. D. J. Am. Chem. Soc. 2005, 127, 1975–1982. 32. Piana, S.; Reyhani, M.; Gale, J. D. Nature 2005, 438, 70–73. 33. Piana, S.; Gale, J. D. J. Cryst. Growth 2006, 294, 46–52. 34. Salvalaglio, M.; Vetter, T.; Giberti, F.; Mazzotti, M.; Parrinello, M. J. Am. Chem. Soc. 2012, 134, 17221–17233. 35. Salvalaglio, M.; Perego, C.; Giberti, F.; Mazzotti, M.; Parrinello, M. PNAS 2015, 112, E6–E14. 36. Zahn, D.; Leoni, S. Phys. Rev. Lett. 2004, 92, 250201–250204. 37. Ozpinar, G. A.; Beierlein, F. R.; Peukert, W.; Zahn, D.; Clark, T. J. Mol. Model. 2012, 18, 3455–3466. 38. Anwar, J.; Frenkel, D.; Noro, M. G. J. Chem. Phys. 2003, 118, 728–735. 39. Benavides, A. L.; Aragones, J. L.; Vega, C. J. Chem. Phys. 2016, 144, 124504/1–124504/14. 40. Vekilov, P. G. J. Cryst. Growth 2005, 275, 65–76. 41. Ostwald, W. Z. Phys. Chem. 1897, 22, 289–330. 42. Retzel, A. J. Z. Krist. 1912, 49, 152–192. 43. Kawska, A.; Brickmann, J.; Kniep, R.; Hochrein, O.; Zahn, D. J. Chem. Phys. 2006, 124, 24513. 44. De La Pierre, M.; Raiteri, P.; Stack, A. G.; Gale, J. D. Angew. Chem. Int. Ed. 2017, 56, 8464–8467. 45. Milek, T.; Duchstein, P.; Seifert, G.; Zahn, D. ChemPhysChem 2010, 11, 847–852. 46. Kawska, A.; Duchstein, P.; Hochrein, O.; Zahn, D. Nanoletters 2008, 8, 2336–2340. 47. Milek, T.; Zahn, D. Nanoletters 2014, 14, 4913–4917. 48. Milek, T.; Duchstein, P.; Zahn, D. CrystEngComm 2015, 17, 6890–6894. 49. Anwar, J.; Zahn, D. Angew. Chem. Int. Ed. 2011, 50, 1996–2013. 50. Milek, T.; Zahn, D. Z. Anorg. Allg. Chem. 2016, 642, 902–905. 51. Pannetier, J.; Bassas-Alsina, J.; Rodriguez-Carvajal, J.; Caignaert, V. Nature 1990, 346, 343–345. 52. Sch€ on, J. C.; Jansen, M. Angew. Chem. Int. Ed. Engl. 1996, 35, 1287–1304. 53. Woodley, S. M.; Catlow, R. Nat. Mater. 2008, 7, 937–946. 54. Gavezzotti, A. J. Am. Chem. Soc. 1991, 113, 4622–4629. 55. Price, S. L. Acc. Chem. Res. 2009, 42, 117–126. 56. Neumann, M. A.; van de Streek, J.; Fabbiani, F. P. A.; Hidber, P.; Grassmann, O. Nat. Commun. 2015, 6, 7793.