Magnesium-rich nanoprecipitates in calcite: atomistic mechanisms responsible for toughening in Ophiocoma wendtii

Abstract

The brittlestar Ophiocoma wendtii is theorised to employ a technique already used in metallurgy in order to optimise the mechanical properties of calcitic microlenses within their skeletons. These microlenses contain arrays of Mg-rich nanoprecipitates, which are proposed to inhibit crack propagation through the compression of the local host lattice. Here, we employ classical molecular dynamics in order to study the effects of Mg-rich nanoprecipitates on lattice strain, stress distributions and crack propagation in calcite. Our quantitative results on lattice strain and stress induced on the host matrix are compatible with empirical estimates. Simulations of crack propagation demonstrate that the inclusion of a Mg-rich region results in an increase in stress required to fracture the crystal, as well as higher residual stress in the fractured crystal. This is the result of an inhomogeneous stress distribution causing a more disordered fracture, as well as deflections of the crack away from the lowest energy (10.4) surface. The results agree with the proposal that the compression of the host lattice inhibits propagation, and offer insight into other mechanisms through which the nanoprecipitates affect crack propagation.

---

1 Introduction

The ability to produce minerals with bespoke properties is crucial for the survival of living organisms. Materials formed through biomineralisation provide functions such as protection,¹ skeletal support² and mastication.³ Biological systems are limited to the chemical composition and ambient conditions of their environment when forming biominerals; nevertheless, evolution has provided living organisms with the means to produce a range of highly sophisticated materials with properties tailored to their respective purposes. An example of a mechanism used in both shells^4,5 and bones⁶ is the formation of complex hierarchical structures of hard and soft materials.⁷ The resulting structures have significantly improved toughness, and often it is mechanisms such as those found in nature that materials scientists aim to exploit to design materials with improved properties.

A more unusual example of a toughening mechanism has been identified in the brittlestar Ophiocoma wendtii. The arm plates of these brittlestars are covered by roughly 10-micrometre sized lenses composed of calcite.⁸ The function of these lenses is to focus light onto photoreceptor nerve bundles positioned beneath the lenses. Each lens is aligned along its optical axis parallel to the c axis of calcite in order to minimise the effects of birefringence. Calcite, being transparent and highly abundant in nature, provides an ideal material for these lenses. However, calcite is famously brittle, with a low fracture toughness. Increasing fracture toughness using the hierarchical structure mechanism described above would be inappropriate, as it would be detrimental to the optical properties of the material. However, in 2017, Polishchuk et al.⁹ identified a new toughening mechanism in the brittlestar when examining the nanostructure of the lenses. Each calcitic lens was found to contain a dense array of magnesium calcite nanoprecipitates, about 4 nm in diameter. These magnesium-rich nanoprecipitates were found to be coherent with the host matrix. This coherency allows the preservation of the optical properties of calcite. The nanoprecipitates are believed to form during crystallisation from amorphous calcium carbonate (ACC) into calcite.¹⁰ The solubility of Mg in ACC is much higher than in calcite, therefore during crystallisation Mg rich regions would be expected to form. Coherent nanoparticles, often referred to as Guinier–Preston (GP) zones,^11,12 are already well known in metallurgy, and have been the subject of previous molecular dynamics studies.^13,14 Prior to the findings of Polishchuk et al., no examples of such a phenomenon had been observed in any living organism. Whereas the formation of GP zones in metallurgy involves extensive heating,¹⁵ the brittlestar is able to achieve this process in calcite under ambient conditions. In metals, GP zones increase the tensile strength of the host material,¹⁶ generally at the expense of increased brittleness. Polishchuk et al. propose that the nanoprecipitates have a different effect on the mechanical properties of calcite: rather than increasing the tensile strength by inhibiting dislocation motion, as in metal alloys, the nanoprecipitates in calcite increase the fracture toughness, by inducing a compressive stress in the host matrix. Such a prestressing mechanism is employed in other brittle materials, such as tempered glass and prestressed concrete. In calcite, the compressive stress in the host matrix is induced by the coherent Mg-rich nanoprecipitates. Due to the small size of the Mg ion compared to the Ca ion, there is a local tensile stress within the nanoprecipitate, causing a compression of the lattice spacing. The continuity of the lattice planes ensure that, where the nanoprecipitates are under a tensile stress, the surrounding matrix is under a compressive stress. It is proposed that this compressive stress inhibits crack propagation and, therefore, increases toughness.

The mechanism proposed by Polishchuk et al. has not, as yet, been confirmed unambiguously as the stress induced on the host matrix, and the resulting effect on crack propagation is challenging to quantify experimentally.⁹ However, molecular simulation allows a direct calculation of the stress field induced by the nanoprecipitate as well as the host crystal's resilience to crack propagation. In this study, we use molecular dynamics to investigate the effects of magnesium incorporation in calcite. We examine the effects of different concentrations of magnesium on calcite lattice parameters. We investigate the hydrostatic stress field and the magnitude of the compensating stress field in the host matrix. Finally, we use crack propagation simulations to examine the effect of magnesium nanoprecipitate incorporation on calcite fracture toughness.

2 Computational method

Molecular dynamics simulations were performed using LAMMPS.¹⁷ The calcite interactions were modelled with the rigid ion version of the Pavese et al.^18,19 force field, since it was fitted explicitly to the elasticity tensor of calcite. The Mg-ion interactions were obtained from Raiteri et al.^20,21 The system timestep was set to 1 fs, and the system was kept at temperature 300 K using a Nosé–Hoover thermostat with a relaxation time of 100 fs. For constant pressure (NPT) simulations, the pressure was set to atmospheric pressure using a Nosé–Hoover barostat with a relaxation time of 1 ps.

Periodic simulation cells with edge lengths roughly equal to 8 nm in all dimensions were populated with calcite. Magnesium ions were introduced into the calcite lattice by randomly substituting for calcium ions. A probability of assignment for Ca and Mg was used to set the percentage of magnesium ions. When modelling spherical nanoprecipitates, the Mg-substitutions were confined to a spherical volume with a diameter of 4 nm, roughly the size of the observed nanoprecipitate. In all simulations, the nanoprecipitates were populated with 40 mol% Mg in order to replicate the empirically found concentration observed by Polishchuk et al.⁹

2.1 Calculating the local stress field

For simplicity, the simulation box was set up such that the x-direction corresponded to the a-axis and the z-direction corresponded to the c-axis of the crystal. The box dimensions were chosen such that the lengths corresponded to 18, 18 and 6 repetitions of the primitive hexagonal cell along the a, b and c axes respectively. All cell dimensions were relaxed under NPT at atmospheric pressure to allow for stress relaxation. The crystal was equilibrated for 100 ps, and a 500 ps simulation was carried out for the calculation of the local stress field.

The definition of stress breaks down in the atomic limit. However, since the stress tensor of a supercell is just the volume-average of the atomic virials, it is common to average the atomic virials over an appropriate local volume to produce a local stress field. Branicio and Srolovitz²² presented a general method for this, later applied to titania nanoparticles²³ and calcite defects,^24,25 that involves calculating the time averaged atomic virial tensors 〈W_αβ〉_i, and multiplying their value by a normalised smearing function P(x⃑ − x⃑_i) where x_i is the position of the atom. A continuous stress field is obtained by summing the product of the virial and smearing function over all N atoms.

Here, P(r − 〈r_i〉) is chosen to be a radial Gaussian, namely:where R is equal to three standard deviations and parametrises the distance above which any contributions to the stress field can be considered negligible (and therefore the resolution of the stress field). Throughout this paper, R was set to 9 Å. α and β label the cartesian components of the virial stress tensor. The hydrostatic stress field, which is equivalent to the negative of the pressure field, is therefore calculated by taking the average of the diagonal elements, Tr(Π_αβ(x⃑))/3. Using the calculated atomic virials from simulations, a grid of points with spacing 1 Å was constructed, and eqn (1) was evaluated at each point.

While the integral of eqn (1) over all space x⃑ recovers the correct stress tensor of the system, interpreting it at a local level is a little delicate. If the local stress tensor is uniform across a volume that exceeds the smearing volume then it is a physically meaningful measure of stress. Otherwise, its value is sensitive to the means of smearing and so it is, at best, an order-of-magnitude estimate. Crucially though, the sign is a faithful indicator of whether the local stress is compressive (negative) or tensile (positive), and it is this feature of the stress field that we are primarily interested in.

2.2 Crack propagation

In order to model a calcite crystal with a pre-existing crack, a thin slab of calcite was simulated. A monoclinic, periodic simulation cell was used to allow periodicity along two of the rhombohedral crystal's three axes. The crystal was oriented such that the free (10.4) surface of the slab was normal to the z-direction of the simulation cell. The length of the cell in the x and y-directions was 75.22 Å. Along the z-direction, the crystal accounted for 81.46 Å. The simulation cell was also periodic in the z-direction, but with a 40 Å vacuum region to emulate free surfaces. At the upper surface of the crystal, an initial crack was introduced to the slab by defining two adjacent regions of atoms. Each region was effectively infinite along the y-direction, had a thickness of 15.7 Å along the x-direction, and penetrated 15.1 Å into the crystal in the z-direction. All interactions between these two regions were disabled so as to create an effective interface that would nucleate a crack upon being strained. During simulations, the Ca ions at the lower face of the slab were constrained to their initial z-coordinates using a harmonic potential with a spring constant of 50 kJ (mol)⁻¹ Å⁻² in order to simulate the effect of bulk crystal at the base of the slab.

Simulations consisted of an equilibration process of 0.2 nanoseconds, where all cell vectors except the z-length were relaxed under NPT at atmospheric pressure. The barostat was then changed such that only the cell length in the x-direction was able to fluctuate at atmospheric pressure. The cell was deformed along the y-direction using a constant engineering strain rate of 0.01 ps⁻¹ over a period of 10 ps. During this 10 ps period, and following a 2 ps equilibration period, the yy-component of the stress tensor was evaluated every 0.1 ps and averaged over the remaining 8 ps. Simulations were repeated 10 times using different velocity seeds and nanoprecipitate configurations (i.e. different seeds for determining the Mg distribution within the nanoprecipitates).

3 Results and discussion

3.1 Force field validation: strain dependence on Mg content

To test the suitability of the force fields, we began with a uniform random distribution of Mg ions in a bulk calcite lattice, and relaxed the cell vectors to eliminate the stress. The resulting cell strain as a function of Mg concentration is shown in Fig. 1. The results are compared with the empirically derived relations for the lattice parameters in the a and c-directions as a function of Mg concentration as found by Polishchuck et al. It is clear from the results that the force fields predict a reasonably accurate relationship between strain and concentration, especially in the more elastic c-direction.

Fig. 1 Comparison of lattice strain in the a direction and c direction between molecular dynamics simulations and the empirically derived relation found by Polishchuk et al.⁹

3.2 Stress distribution for calcite with 40 mol% Mg nanoprecipitates

The local hydrostatic stress field was computed for a cross-section of the simulation cell that bisected the centre of the nanoprecipitate, shown in Fig. 2. The radial and angular components of the stress tensor field are shown in Fig. 3. Fig. 2 shows a high tensile stress within the nanoprecipitate, caused by the smaller size of the Mg ions, as can be seen in the red region encompassing the nanoprecipitate in Fig. 2. The surrounding host matrix experiences a compressive stress in order to compensate for the tensile stress, as can be seen in the blue regions outside the nanoprecipitate. The hydrostatic pressure appears roughly homogeneous throughout the bulk crystal, although Fig. 3 sheds more light on the complexity of the stress distribution: the radial stress distribution (left) shows a host matrix almost entirely under tensile stress. Furthermore, the radial stress is larger in magnitude nearer the nanoprecipitate. The tangential component, however, shows an entirely compressive stress. It is apparent, therefore, that the stress encountered by a propagating crack will be dependent on where it is relative to the nanoprecipitate. These results provide insight into how crack propagation may be inhibited in some regions, and deflected in others, as elaborated in the following section.

Fig. 2 Hydrostatic stress tensor for a periodic (bulk) calcite crystal with an incorporated 40 mol% Mg nanoprecipitate of diameter 4 nm. Atoms from a cross-section of the cell are also shown. Ca is shown in green, C in grey, O in red and Mg in yellow. Here, red represents regions under tensile stress and blue represents regions under compressive stress. The cell dimensions are 86 Å in the x-direction and 105 Å in the z-direction. The dashed circle indicates the Mg-rich domain. The hydrostatic stress field is roughly homogeneous throughout the host crystal.

Fig. 3 Radial (left) and tangential (right) components of the stress tensor for the same system as Fig. 2 in cylindrical coordinates. The smaller lattice spacing within the nanoprecipitate causes a tensile radial stress (red). This causes an angular compression (blue).

The calculation of the atomic virials also allows the average stress induced in the host matrix to be calculated. By excluding all atoms within the defined spherical region, and accounting for the resulting volume change, the stress tensor components for regions outside the nanoprecipitate (σ_αβ) can be calculated by averaging the remaining atomic virials over the remaining volume. The results of this process for the xx and zz components, along with the hydrostatic stress, are given in Table 1, where they are compared with the empirical estimates of Polishchuck et al.⁹ derived from a continuum elasticity model. While the values are consistently lower, the agreement is generally quite impressive considering the various approximations in both data sets.

Table 1 Comparison of stress tensor components for the host matrix between simulations and estimates of Polishchuck et al.⁹

	Simulations	Polishchuck et al.
σ xx (MPa)	−145	−180
σ zz (MPa)	−127	−140
σ H (MPa)	−139	−170

---

3.3 Crack propagation

Fig. 4 shows visualisations of crack propagation in calcite under several different conditions. In the control situation, where the calcite contains no Mg (Fig. 4(a)), the crack propagates cleanly down the low-energy (10.4) plane and, with the exception of the constrained region at the lower surface, the two surfaces are completely separated. The inclusion of the nanoprecipitate, however, brings about a visibly significant change in the nature of the crack propagation. When the crack bisects the nanoprecipitate (Fig. 4(b)), propagation still occurs along the (10.4) plane, but the fracture is far more disordered, with ions bridging the two resulting surfaces. The Mg-rich domain itself is also observed to fail before the crack has reached the precipitate, presumably because of the tensile stress in the nanoprecipitate. When the crack reaches the nanoprecipitate slightly off-centre, as shown in Fig. 4(c), the crack changes direction when propagating through the Mg-rich domain, producing a rough, high-energy interface. This again may be attributed to the high tensile stress within the nanoprecipitate. Where, as shown in Fig. 4(d), the crack propagates around the nanoprecipitate, the crack is deflected from the low-energy (10.4) plane towards the nanoprecipitate. In light of the stress distribution in Fig. 3, it is unsurprising that the crack would be deflected towards the nanoprecipitate, since the radial stresses are highest near the Mg-rich domain, and a crack propagating through this region is able to relieve said radial stress. This reveals an interesting consequence of the presence of a Mg-rich domain: the complex nature of the stress distribution may result in a deflection of crack propagation away from the lowest energy (10.4) surface, which will contribute to the increased toughness.

Fig. 4 Snapshots of a crack propagating under different conditions including: crack propagating through pure calcite (a); crack bisecting the nanoprecipitate (b); crack propagating slightly off-centre of the nanoprecipitate (c); crack propagating around the nanoprecipitate. (d) Here, the Mg ions have been magnified for easier visualisation.

The stress strain curves calculated during crack propagation are shown in Fig. 5. Note that the stress values are averaged over 10 simulations for each point. The points of maximum stress on these curves correspond to the initiation of the crack propagation. The curves corresponding to Fig. 4(a) and (d) are indistinguishable before the point of fracture, suggesting that the precipitates have little influence on cracks that do not approach closely. The curves corresponding to Fig. 4(b) and (c) display two important features. Firstly a higher tensile stress (an additional 223 and 180 MPa respectively) is reached before crack propagation, and secondly some stress is retained at the end of the simulation. The increase in maximum tensile stress is the result of the compressive stress in the matrix, as suggested by Polishchuk et al. The residual stress is a result of the incomplete fracture, as the fracture surfaces are bridged by a disordered region of ions at the end of the calculation. Both of these effects would increase the fracture toughness. The crack that propagates off centre (Fig. 4(d)) does not display a significant increase in the tensile stress at the initiation of crack propagation but it does show significant residual stress, partially due to bridging of the surfaces.

Fig. 5 Stress–strain plots averaged over 10 simulations. An insignificant difference with respect to pure calcite is seen when the crack is far from the nanoprecipitate (d), but the difference becomes more apparent when the crack directly bisects the nanoprecipitate (b) and (c). In each case, the residual stress is higher when a nanoprecipitate is present.

As seen in Fig. 4: rather than the crack propagating cleanly down a (10.4) surface, the breaking of the crystal is more complex; this is clearly an effect of the crack being deflected towards, and even through, the Mg-rich domain. This result demonstrates a possible inhibitory mechanism for crack propagation not previously reported to our knowledge, and may suggest a novel mechanism for improving fracture toughness. Polishchuk et al. proposed a toughening by which cracks are deflected on a macroscale due to variations in the material density.⁹ Our observations suggest that a nanoscale deflection may also be an important part of the story.

It is worth emphasising that the stress and strain at the point of crack propagation is only increased when the crack directly bisects the nanoprecipitate, as can be observed in Fig. 5. This can be explained by the results in Fig. 3, which demonstrate that different regions of the crystal will be under different stresses. As a crack propagates, the stress opposing propagation will be the stress perpendicular to the direction of propagation. If the crack propagates towards the precipitate, then it will encounter the compressive tangential stress observed in Fig. 3, whereas if the crack propagates around the precipitate, the relevant component of the stress tensor field becomes the radial component. This is why the stress at the point of propagation of the crystal is only increased when the crack bisects the nanoprecipitate. The difference of 223 MPa observed in Fig. 5 is of the same order of magnitude as the tangential stress observed in the host matrix in Fig. 3, indicating that the prestressing of the host matrix does contribute to the increased fracture toughness. This relies on the assumption that the crack approaches the nanoprecipitate, although it can be assumed that, given the high precipitate density, a propagating crack through a real crystal would, at some point, approach a nanoprecipitate.

4 Conclusions

We have used molecular dynamics to investigate the effects of magnesium-rich nanoprecipitates on the stress distribution and crack propagation in calcite domains. Such nanoprecipitates have been identified in calcitic microlenses within the skeleton of the brittlestar Ophiocoma wendtii, and are observed to increase the fracture toughness in comparison to non-biogenic calcite. Our simulations provide unique insight into the atomistic mechanisms responsible for the increased toughness. We have identified three possible sources. Firstly, the small size of the magnesium ions compared to calcium, together with the coherent nature of the nanoprecipitates, induces a compressive stress in the surrounding matrix, which counteracts the applied tensile stress during crack propagation, similar to the toughening of tempered glass and prestressed concrete. This mechanism was suggested by Polishchuk et al. who also estimated the magnitude of this effect using continuum mechanics. Here, we were able to calculate the components of the stress tensor at atomic resolution and show that, although the hydrostatic stress is relatively homogeneous throughout the matrix, the radial and tangential components show separated regions of compressive and tensile stress. As well as providing an explanation for why crack propagation is inhibited only when the crack approaches the nanoprecipitate, the inhomogeneous stress distribution is significant because it induces deflections of the cracks away from the lowest energy surfaces. This second source contributes to increased energy adsorption during fracture and, consequently, increased fracture toughness.

The third contribution to increased fracture toughness identified by the simulations is the mechanism by which the cracks propagate through the nanoprecipitates. The highly inhomogeneous stress distribution within the precipitates result in strongly disordered fracture surfaces which form bridges between the crack surfaces that support residual stress at the termination of the simulation. Such crack bridges would adsorb energy during fracture and may be a significant contribution of the nanoprecipitates to increased toughness.

In summary, our atomistic simulations of stress and fracture in calcite with embedded Mg-rich nanoprecipitates have identified novel mechanisms by which such nanoprecipitates may increase the toughness of calcite.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by an Engineering and Physical Sciences Research Council Programme grant (EP/R018820/1). AB acknowledges support of an EPSRC Doctoral Training Partnership studentship (EP/R513143/1).

References

1. J. C. Weaver, G. W. Milliron, A. Miserez, K. Evans-Lutterodt, S. Herrera, I. Gallana, W. J. Mershon, B. Swanson, P. Zavattieri, E. DiMasi and D. Kisailus, Science, 2012, 336, 1275–1280.
2. A. H. Heuer, D. J. Fink, V. J. Laraia, J. L. Arias, P. D. Calvert, K. Kendall, G. L. Messing, J. Blackwell, P. C. Rieke, D. H. Thompson, A. P. Wheeler, A. Veis and A. I. Caplan, Science, 1992, 255, 1098–1105.
3. L. M. Gordon, M. J. Cohen, K. W. MacRenaris, J. D. Pasteris, T. Seda and D. Joester, Science, 2015, 347, 746–750.
4. J. D. Currey, Proc. R. Soc. B, 1977, 196, 443–463.
5. F. Barthelat, H. Tang, P. D. Zavattieri, C. M. Li and H. D. Espinosa, J. Mech. Phys., 2007, 55, 306–337.
6. S. Weiner and H. D. Wagner, Annu. Rev. Mater. Sci., 1998, 28, 271–297.
7. P. Fratzl and R. Weinkamer, Prog. Mater. Sci., 2007, 52, 1263–1334.
8. J. Aizenberg, A. Tkachenko, S. Weiner, L. Addadi and G. Hendler, Nature, 2001, 412, 819–822.
9. I. Polishchuk, A. A. Bracha, L. Bloch, D. Levy, S. Kozachkevich, Y. Etinger-Geller, Y. Kauffmann, M. Burghammer, C. Giacobbe and J. Villanova, et al. , Science, 2017, 358, 1294–1298.
10. E. Seknazi, S. Kozachkevich, I. Polishchuk, N. B. Stein, J. Villanova, J.-P. Suuronen, C. Dejoie, P. Zaslansky, A. Katsman and B. Pokroy, Nat. Commun., 2019, 10, 1–9.
11. A. Guinier, Nature, 1938, 142, 569–570.
12. G. D. Preston, Proc. R. Soc. A, 1938, 167, 526.
13. W. Verestek, A.-P. Prskalo, M. Hummel, P. Binkele and S. Schmauder, Phys. Mesomech., 2017, 20, 291–304.
14. C. Singh and D. Warner, Acta Mater., 2010, 58, 5797–5805.
15. S. P. Ringer and K. Hono, Mater. Charact., 2000, 44, 101–131.
16. D. M. Duffy, Science, 2017, 358, 1254–1255.
17. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
18. A. Pavese, M. Catti, G. D. Price and R. Jackson, Phys. Chem. Miner., 1992, 19, 80–87.
19. A. Pavese, M. Catti, S. C. Parker and A. Wall, Phys. Chem. Miner., 1996, 23, 89–93.
20. P. Raiteri, J. D. Gale, D. Quigley and P. M. Rodger, J. Phys. Chem. C, 2010, 114, 5997–6010.
21. P. Raiteri, R. Demichelis and J. D. Gale, J. Phys. Chem. C, 2015, 119, 24447–24458.
22. P. S. Branicio and D. J. Srolovitz, J. Comput. Phys, 2009, 228, 8467–8479.
23. R. Darkins, M. L. Sushko, J. Liu and D. M. Duffy, Phys. Chem. Chem. Phys., 2014, 16, 9441–9447.
24. A. S. Cote, R. Darkins and D. M. Duffy, Phys. Chem. Chem. Phys., 2015, 17, 20178–20184.
25. Y.-Y. Kim, R. Darkins, A. Broad, A. N. Kulak, M. A. Holden, O. Nahi, S. P. Armes, C. C. Tang, R. F. Thompson and F. Marin, et al. , Nat. Commun., 2019, 10, 1–15.

---