Deformation twinning and the role of amino acids and magnesium in calcite hardness from molecular simulation

We employ classical molecular dynamics to calculate elastic properties and to model the nucleation and propagation of deformation twins in calcite, both as a pure crystal and with magnesium and aspartate inclusions. The twinning is induced by applying uniaxial strain to the crystal and relaxing all stress components except the uniaxial component. A detailed analysis of the atomistic processes reveal that the twinning mechanism involves small displacements of the Ca ions and cooperative rotations of the CO 3 ions. The volume of the twinned region expands under increased uniaxial strain via the propagation of steps along the twin boundaries. The energy cost of the twin boundaries is compensated by the reduced hydrostatic stress and strain energy. The presence of biogenic impurities is shown to decrease the strain required to induce twin formation in calcite and, thus, the yield stress. This increased propensity for twinning provides a possible explanation for the increased hardness and penetration resistance observed experimentally in biominerals.


Introduction
5][6][7] One of the most widely studied biominerals is calcite, the most stable polymorph of CaCO 3 .Apart from being a rock-forming mineral abundant in the Earth's crust, calcite is also biogenic.It often forms part of the skeleton or exoskeleton of invertebrates, where exceptional hardness and toughness are required for protection and support.3][14] One way that this is achieved is deformation twinning, as was recently shown using nanoindentation on a biomineralised calcite exoskeleton. 15eformation twinning is a mechanism which occurs in some crystals when they are subjected to a sufficiently large shear stress.Crystalline twins are characterised by the formation of a mirror plane which separates the original ''parent'' and ''twinned'' regions.The highly co-ordinated atomic movement that leads to the twin boundary formation has been extensively studied for metals and minerals. 16,17Twinning has a very important role in a range of physical properties of crystals, as well as growth and phase transformations.It has been shown that, at the nanoscale, twin boundaries increase hardness in the same way as grain boundaries 18,19 by inhibiting dislocation motion, while gliding of dislocations along densely distributed twin boundaries increases fracture toughness.This was recently exemplified by the creation of nanotwinned diamond with exceptional hardness and thermal stability. 202][23][24] In this study we examine twinning in calcite using molecular dynamics.We analyse the formation and development of the twin boundaries that are created by the application of a uniaxial strain and calculate the resulting stress-strain curve.We repeat the simulations in the presence of biogenic impurities known to affect mechanical properties (Mg ions and aspartate), [25][26][27] in order to investigate the superior mechanical properties of biogenic calcite.deformation or as defects during crystal growth.e and r are the most common deformation twinning planes, depicted in Fig. 1.It should also be noted that the r-plane is the usual slip plane.The e-twin is the one most commonly observed experimentally, and it can be easily introduced by indentation. 15,280][31][32] Other carbonate minerals, such as dolomite CaMg(CO 3 ) 2 , have also been analysed in terms of their twin structure and necessary shearing directions to produce twins. 33,34However, the atomistic mechanism of nucleation and propagation of twin boundaries in calcite has not been studied computationally and this is one of the aims of the current study.

Computational method
Molecular dynamics calculations were performed in the DL_POLY Classic code (v.1.9), 35 using a rigid ion version of the interatomic potential by Pavese et al. 36,37 This potential was found to reproduce the mechanical properties of calcite better than the more recent potential by Raiteri et al., 38,39 which is unsurprising since it was specifically fitted to elastic and vibrational data.
The periodic supercell used in the simulations was a 7 Â 4 Â 3 multiplication of the conventional orthorhombic cell of calcite, containing 1008 CaCO 3 units, and with the z-axis corresponding to the c-axis of the crystal.The electrostatics were handled by the Ewald summation method and the short-range potentials had a cutoff of 10.1 Å.The simulation timestep was 1 fs, using the Nose ´-Hoover thermostat at 300 K with a 100 fs relaxation time.The initial system was relaxed to equilibrium using the NsT ensemble with a 1 ps barostat, for a duration of 1 ns.
Small concentrations of impurities were added to the supercell to examine their effect on the elasticity and deformation twinning.We added 5 and 10 mol% Mg in our supercell by random substitution of Ca 2+ ions, and modelled the Mg interactions with the potential by de Leeuw for MgCO 3 . 40We repeated the simulations with 0.4 and 0.8 mol% aspartate (Asp).The zwitterionic form of aspartate (Asp 1À ) was used ( + H 3 N-CH-(CH 2 CO 2 À )-CO 2 À ), and its potential was generated with AMBER. 41The interactions between aspartate and calcite were fitted using the methods described in ref. 42.The configuration, in this case, was obtained by using metadynamics to drive the crystallisation of disordered calcium carbonate that contained the aspartate molecules.The resulting structure is shown in Fig. 2 where it can be seen that the amino group in the aspartate takes the place of a Ca ion in the calcite lattice, while the two carboxylic acid groups replace two carbonate groups, hence the aspartate fits in the calcite lattice with minimal disruption.

Elastic properties
We used the small strain method to calculate the six independent elastic constants (C 11 , C 33 , C 44 , C 12 , C 13 , C 14 ) of calcite, with and without Mg and Asp impurities.Small (up to AE1%) tensile, compressive and shear strains were applied to the equilibrium lattice vectors, and the resulting stresses calculated. 43This was repeated for 20 strain increments close to the equilibrium vectors, with each strain direction corresponding to one or more elastic constants.Specifically, given an initial orthorhombic calcite supercell with axes a, b and c, the three matrices chosen to define the axes of the deformed supercell were the following, with an increment d = 0.001 Å: a da 0 from which all the elastic constants were obtained.The elastic constants were extracted by a linear fit to the stress-strain data (s i = C ij e j ).This is possible because the strain is sufficiently small that the deformation stays in the linear regime.The stresses were calculated by sampling the canonical ensemble for 0.2 ns after each strain increment.The bulk and shear moduli were calculated from the elastic constants using the Voigt-Reuss-Hill averaging scheme. 44

Deformation twinning
A similar method was used to deform the supercell along the z-direction, i.e. the c-axis of the unstrained and relaxed unit cell.The implementation of the NsT ensemble in DL_POLY was modified in order to constrain the deformation in the z direction only.All other components of the cell vectors were relaxed to relieve all stress components except s z .Thus, in general, all components of the strain vector are finite and all components of the stress vector are zero except s z .We are effectively simulating an experimental setup in which only e z is constrained.
The supercell, including the atomic coordinates, was stretched in small incremental steps at a rate of 0.15 m s À1 , up to B30% strain.The strain increase at each step was 0.28% and a 100 ps constant strain simulation at 300 K was carried out after each step.The final configurations of these 100 ps simulations were stretched and used as the starting configuration for the  This journal is © the Owner Societies 2015 next step.The z-component of the stress -the only non-zero component -was calculated after each strain increment by averaging over the final 98 ps of the 100 ps constant strain simulations.The calculated values were used to plot the uniaxial stress-strain curve.The atomistic configuration was also captured after each strain increment and the atomic motion involved in twin formation and growth was studied in detail.

Calculation of local stress field
In addition to calculating the average uniaxial stress for the simulation cell, we calculated the local hydrostatic stress distribution throughout the cell.This gives information about the stress distribution near the twin boundaries and the stress introduced by the impurities.We also use this to compare the stress distribution in the simulation cell before and after twinning.
A method for calculating stress fields in ionic crystals from molecular dynamics was introduced by Branicio and Srolovitz 45 and applied by Darkins et al. 46 The continuous stress field is derived from the time-averaged atomic virials, hW i abi, using where all species X and atoms i are summed, and P(r) is a threedimensional normalised Gaussian of the form where R parametrises the size of the Gaussian and is equivalent to three standard deviations.In our calculations R was set to 4 Å. a and b take on values x, y, z to generate the components of the virial stress tensor.The supercell is divided into a grid, and the hydrostatic stress field Tr(P(r))/3, which is a scalar field equivalent to the negative of the pressure field, is calculated at every grid point using eqn (1).By using a dense grid of 1 Å divisions, we can examine the stress distribution throughout our system, both before and after deformation.The calculations were performed at a very low temperature (10 K) to reduce thermal noise.

Elastic properties
The calculated elastic constants of calcite, with and without Mg and Asp impurities, are presented in Table 1.The calculated values are in good agreement with both the values calculated previously with a different methodology 37 and the measured values.These values give us confidence that the potentials we use are appropriate for modelling the mechanical properties of calcite.
The addition of Mg makes calcite stiffer -the 10 mol% impurity caused a 2% increase in C 11 .This increase in stiffness is in agreement with studies on dolomite, which has considerably higher elastic moduli compared to calcite. 47Asp has the opposite but larger effect, achieving approximately the same change by including only 0.8 mol% of impurity.As the concentrations used here were small, the symmetry of the elastic matrix remained unaffected.

Deformation twinning
The calculated stress-strain curve for pure calcite up to a uniaxial strain of 25% is shown in Fig. 3.The curve is linear up to 2% strain and the gradient gives a tensile modulus of 66.4 AE 0.2 GPa, which is in excellent agreement with the value calculated from the elastic constants in section 4.1 above.The curve is close to linear up to a strain of 10.8%, at which point there is a sharp drop in stress, signalling the formation of a pair of {104} twin boundaries separated by 2-3 atomic layers.
The stress at which the stress-strain curve deviates significantly from linear behaviour, the yield stress, was found to be 7.7 GPa.We note that the calculated yield stress and strain are much higher than would be achieved in a real calcite crystal since, in our model, there are no surfaces, cracks or other flaws that would usually serve as nucleation sites for such failure.
The atomistic configuration of the simulation cell immediately after the formation of the pair of deformation twins is shown in the second snapshot of Fig. 3.The twinned volume between the two twins is two or three atomic layers thick and has a different crystallographic orientation to the original (parent) crystal.The change in the number of layers is the result of ''steps'' in each of the twin boundaries.There can be one or more steps in each twin, and their number increases with strain.We repeated this experiment using a simulation cell that was 1.5 times smaller in each direction, and we found that r-twins formed at the same stress/strain values as in the original cell; the only difference was that the twin boundaries spanned the entire cell, contrary to the ones shown in Fig. 3, which end part way through.

Twin propagation
The stages of deformation were identified by analysing the trajectories of each simulation.The motion of the twin boundary steps results in widening of the twinned region, with increased tensile strain.As the lattice is strained further, the steps move along the boundary, making the twinned region wider with respect to the parent region.When the step reaches the end of the twin boundary, a step nucleates on the opposite side of the twinned region and continues to propagate, further widening the twin region.Steps may also form spontaneously at the end of the twin boundaries.When two steps meet, the two twinned areas join, and in this instance, due to the periodic boundary conditions, the ''parent'' region is finally consumed and the entire structure adopts the crystallographic orientation of the ''twin'' region.In the pure calcite lattice this happens at B25% strain.The a angle of the supercell also gradually changes to accommodate the new periodic structure, from 901 to 1011.Snapshots of the process are illustrated in Fig. 3, together with the corresponding stress-strain relationship.Whenever the steps move a large distance or annihilate, we notice a corresponding drop in stress.Interestingly, near 25% strain, when the parent regions disappear completely, the structure is in a compressed state (negative stress).The energy cost of the twin boundaries is larger than that of a twin-free lattice with stored elastic energy.Specifically, using the initial unstrained lattice as reference, the energy per unit volume drops from 237 MJ m À3 to 64 MJ m À3 , before and after the disappearance of the twins.We repeated the simulations with low concentrations of Mg (5 mol%) and Asp (0.4 mol%) impurities.The twin boundaries that formed following the initial drop in stress are shown in Fig. 4.
In the Mg case, the twin boundaries started forming in the region of the supercell where the Mg concentration was higher than the cell average, around 9%.Similarly, in the simulation with Asp impurities, the twinned region initially formed close Fig. 3 Deformation twinning and widening of the twinned region caused by tensile strain along the z-axis of the calcite supercell.The stress-strain relationship and four configurations are shown between 11% and 25% strain.The final configuration has entirely adopted the crystallographic orientation of the twinned region.The dashed blue arrows show the direction of step propagation, which results in widening of the twinned volume.The twin boundaries are drawn as black lines.The thin dashed lines denote the end of the twinned region, as well as the steps in the boundaries.
to the amino acid molecules.We attribute this to high localised stress fields around the impurities and we discuss this further in Section 4.4.
The calculated stress-strain curves for calcite with Mg and Asp impurities are shown in Fig. 5.The curve for pure calcite is included for comparison.We note that the initial twin formation occurs at lower strain (9.7% for Mg, 9.4% for Asp) and stress (6.7 GPa for Mg, 6.1 GPa for Asp) values in the presence of impurities.This is an important result as it implies that impurities enhance boundary formation and, consequently, reduce the yield stress of calcite.

Stress field
We calculated the hydrostatic stress fields in our supercells with both Mg and Asp impurity ions.A volumetric representation of the stress field in the 5 mol% Mg supercell, just before and just after the twin boundaries form, are shown in Fig. 6.
As expected, the volume with the highest local stress before deformation is the volume with the largest impurity concentration.
The Mg 2+ ions, being smaller than the Ca 2+ ions they substituted, cause a compressive stress around them in the lattice.When there is a high localised impurity concentration, the stress in that volume is higher and the deformation twin nucleates close to that region in order to relieve the stress.The compressive (positive) stress around the Mg 2+ ions at the moment just before twinning is visible in Fig. 6.After the twin has formed, the previously-compressed region is now replaced by a pocket of negative stress.The same effect was observed when Asp was the impurity; the twins form near the impurities and act to relieve the local accumulation of stress.

Conclusions
We have used molecular dynamics to calculate the stress-strain curve for calcite and to observe the process of deformation twinning induced by uniaxial strain in the (001) direction.We observe elastic deformation up to a strain of around 10%.We note that this is much higher than would be expected in an experimental sample where grain boundaries and other defects introduce stress concentrations that initiate failure.Our simulations are instead representative of a perfect calcite crystal and provide detailed information about intrinsic deformation mechanisms.
At a strain of 10.8%, and a stress of 7.7 GPa, a pair of {104} twin boundaries form by small displacements of the Ca ions and displacement and cooperative rotation of the CO 3 ions.The twin boundaries are initially separated by two to three crystal layers with a different orientation to the original crystal, with steps in the twin boundaries separating the two-and three-layer regions.As the strain is increased further, the volume of the Fig. 5 Comparison of the stress-strain relationship between the system with 5% Mg, 0.4% Asp, and the pure calcite system.Twinning occurs at lower stress when impurities are introduced in the calcite lattice.We repeated the simulations with Mg and aspartate impurities, which are additives commonly found in biogenic calcite.The calculated stress-strain curves for the impure calcite crystals were similar to that for pure calcite.However, the initial drop in stress, corresponding to the formation of the pair of twin boundaries, occurred at significantly lower strains and stresses (9.7% and 6.7 GPa for Mg and 9.4% and 6.1 GPa for Asp).The twin boundaries formed initially close to the regions with the highest impurity concentration and analysis of the local stress distribution around the impurities before and after twin formation showed that twin formation affects the local stress associated with the defects.
The fact that impurities enhance deformation twinning may explain the enhanced hardness and penetration resistance of biogenic calcite in nanoindentation experiments. 15,48The inclusion of impurity ions, amino acids and proteins will result in a higher density of twin boundaries in response to the application of stress.The twin boundaries, in turn, will inhibit dislocation motion and nano-crack propagation, increasing both hardness and penetration resistance.Deformation twinning may therefore prove to be a key mechanism in the effort to engineer biomimetic minerals with superior mechanical properties.
In summary, we have used molecular dynamics to investigate the response of calcite to uniaxial strain, with and without added impurities.We have found that plastic deformation occurs via the formation of pairs of {104} twin boundaries and proceeds via the motion of steps along these boundaries.Added impurities reduce the stress and strain required for twin formation, which may contribute to the enhanced hardness of biogenic calcite.

Fig. 2
Fig. 2 Snapshots before (left) and after (right) the inclusion of Asp in the calcite lattice.The positions of the substituted ions are shaded.

Fig. 4
Fig. 4 (a) {104} twins in calcite with 5 mol% Mg.The shaded regions shows the area with the largest Mg 2+ concentration.The twin nucleates from that area.Mg 2+ ions are coloured yellow.(b) {104} twins with Asp 1À ; here the shaded circles mark the position of the Asp 1À molecules.

Fig. 6
Fig. 6 Volume rendering of hydrostatic stress (a) before and (b) after deformation twinning.Calcium ions are not shown for clarity; magnesium ions are shown in yellow.The deformed structure on the right corresponds to Fig. 4(a).

Table 1
47lculated elastic constants of pure calcite and with Mg and Asp impurities.The values in brackets are the experimental values for pure calcite from Chen et al.47