Compressional behavior of the aragonite-structure carbonates to 6 GPa

The behaviors of aragonite (CaCO3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}), strontianite (SrCO3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}), cerussite (PbCO3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}), and witherite (BaCO3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}) at increasing pressure have been studied up to 6 GPa using density functional theory with plane waves. A parallelism of the orthorhombic carbonates with the closed-packed AsNi structure is considered in our analysis, being the CO32-\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3^{2-}$$\end{document} groups not centered in the interstice of the octahedron. The decomposition of the unit-cell volume into atomic contributions using the Quantum Theory of Atoms in Molecules has allowed the analysis of the bulk modulus in atomic contributions. The bulk, axes, interatomic distances, and atomic compressibilities are calculated. The largest compression is on the c crystallographic axis, and the c linear modulus has a linear function with the mineral bulk modulus (K0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_0$$\end{document}). Many of the interatomic distances moduli of the alkaline earth (AE) carbonates show linear functions with the bulk modulus; however, the whole series (including cerussite) only gives linear functions when K0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_0$$\end{document} is related either with the CC distances modulus or the modulus of the distances of the C to the faces of the octahedron perpendicular to c. These last distances are the projections of the Metal–Oxygen (MO) distances to the center of the octahedron. K0AE\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$K_{0AE}$$\end{document} carbonates also show linear functions with the atomic moduli of their cations. However, the whole series show a linear relation with the atomic modulus of C atoms. Therefore, the whole series highlight the importance of the C atoms and their interactions in the mechanism of compression of the orthorhombic carbonate series.


Introduction
Carbonates are an important group of minerals that are present in the Earth's crust forming mainly part of sedimentary and metamorphic rocks. They are important in atmospheric and earth surface processes such as the global carbon cycle. The saturation of the oceans in carbonate minerals as calcite, aragonite, and dolomite has significant implications in the climate change. Carbonate rocks can also be involved in the subduction zones and may persist in the slab at depths greater than those where dehydration reaction occurs. Moreover, carbonate minerals are even found in peridotite, kimberlite, and xenoliths (Berg 1986;van Achterbergh et al. 2002;Zanetti et al. 1999;Ducea et al. 2005). Magnesite and some stable phases of CaCO 3 are present in the Earth's lower mantle down to at least 1500 km depth (Katsura and Ito 1990;Santillán and Williams 2004). The results of experiments conducted in a laser-heated diamond anvil cell indicate the stability of carbonates to 50 GPa (Biellmann et al. 1993) and especial phases are even stable to more than 200 GPa (Pickard and Needs 2015). The presence of carbon in different chemical phases in the Earth's lower mantle indicates that one of the important fate of CO 2 could be the Earth's mantle and to be recycled to surface throughout the volcanic processes. For these reasons, the effect of pressure on the carbonate structures has been widely investigated [ (Martens R 1982;Catti et al. 1993;Martinez et al. 1996;Chung-Cheng and Lin-Gun 1997b, a;Holl et al. 2000;Ono et al. 2005;Wang et al. 2015;Merlini et al. 2016;Palaich et al. 2016) and references therein].
Aragonite (CaCO 3 ), strontianite (SrCO 3 ), cerussite (PbCO 3 ), and witherite (BaCO 3 ) are naturally isostructural orthorhombic group that belong to the space group Pmcn. This mineral group is also known as the aragonite group. Among these minerals, aragonite is the most frequent in nature, thermodynamically stable at high pressures (Wang et al. 2015;Carlson 1983), forming the first high-pressure polymorph of CaCO 3 . Aragonite is more stable than calcite at approximately 0.3 GPa (Johannes and Puhan 1971). Calcium carbonate precipitates as metastable aragonite polymorph in marine environments, rather than the thermodynamically stable phase calcite. Aragonite is also found in the Earth's mantle and in the subducted slabs (Wang et al. 2015;Ono and Mibe 2013). At high pressure (40 GPa), a phase transition is found going to post-aragonite, a trigonal structure (Santillán and Williams 2004). By ab initio calculations, different phases at higher pressures were found (Pickard and Needs 2015;Oganov et al. 2006Oganov et al. , 2008Arapan et al. 2007;Arapan and Ahuja 2010). Regarding the effects of pressure, the isotropic compression of a system always produces an increase of the density if the system's mass remains constant. Bond length shortening is expected when pressure is applied to all type of crystals, although its magnitude varies with the different polyhedrons formed by cations and their anions. Some authors have remarked that the compressibility of a polyhedron is directly proportional to the cube of the distance between ions and therefore approximately to the volume of the polyhedron (Hazen and Finger 1985). The marked shortening of the longest ionic Ba-O bonds in the witherite structure with pressure (3.2% between room pressure and 7.05 GPa) compared with the other orthorhombic carbonates has been considered as the cause of its higher compressibility along the c crystallographic axis (5.14%, from 0.2 to 5.1 GPa) (Holl et al. 2000;Wang et al. 2015). On the other hand, the a axis is only reduced by 1.47% and the b axis by 1.68% in this range of pressures (Wang et al. 2015). These different linear incompressibility moduli indicate this series is highly anisotropic in front of the pressure, especially that concerning to the c axis. Lin and Liu (Chung-Cheng and Lin-Gun 1997b, a) found phase transitions for the strontianite, cerussite, and witherite to a high-pressure phase at 35 GPa, 17 GPa, and 8 GPa, respectively. Holl et al. (Holl et al. 2000) found a high-pressure phase for witherite, being the phase transition pressure of the first order at 7.2 GPa, and the new phase was trigonal ( P − 31c).
Previous works performing ab initio calculations on aragonite have been done. Oganov et al.(Oganov et al. 2006) described theoretically a phase transition from aragonite (Pmcn) to a postaragonite form (pyroxene-type structure C222 1 ), where the carbonate group transforms from a planar triangular structure (CO 2− 3 ) to a tetrahedral arrangement (CO 4− 4 ) above 137 GPa. Other calculations made by Arapan et al. (Arapan et al. 2007) by means of ab initio calculations found sp 3 hybridized bonds and a especial stability of CaCO 3 at very high pressure.
Taking into account the previous references, we have revisited the aragonite-type carbonates, also including cerussite by performing DFT calculations. The effects of pressure to 6 GPa on structure, atomic groups and atomic parameters are studied. Besides, in order to give a deeper insight into the compressibility mechanism, to understand and systematize the aragonite group minerals, the Quantum Theory of Atoms in Molecules [QTAIM, (Bader 1990)] is applied by performing a partition of the unit-cell volume in terms of atomic contributions, and from these atomic volumes the atomic bulk moduli are calculated. This type of partition has been successfully used by some authors in the analyses of spinels (Pendás et al. 2000) or polymorphs of silica (SiO 2 ) (Morales-García et al. 2014).

Crystal model
The structure of these orthorhombic carbonates consists of different layers of cations that alternate with corrugated layers of CO 2− 3 along the c axis ( Fig. 1). This structure was early related to the closed-packed structures of the NiAs (Ewald and Hermann 1931;O'Keeffe and Hyde 1985) (Fig. 2a) by assuming a replacement of As by metals (M) and Ni by C (the oxygen atoms bonded to C are considered as a stuff (O'Keeffe and Hyde 1985). Therefore, it can be envisaged as closed-packed sheets of M and CO 2− 3 . Comparing the structure of aragonite-type minerals with the NiAs structure, it is remarkable that whereas the Ni atoms are centered in the octahedral interstices between the As atoms in the hexagonal closest-packed arrangement the CO 2− 3 anions are not (Jarosch and Heger 1986;Bevan et al. 2002). They are shifted from the center of the octahedron formed by the metal atoms ( Fig. 2b) toward one of the two octahedral faces along the c axis. Thus, the CO 2− 3 layers appear corrugated. The displacement along the c axis of the C atoms from the center of the octahedral interstices allows to define two distances from the carbon atom to the two faces of octahedron (see Fig. 2b): (i) one long distance from the carbon atom to the upper face defined by the three cations directed toward the three oxygens bonded to C ( d a ); and (ii) one short distance from the carbon atom to the bottom face defined by the three cations directed between the three oxygens atoms bonded to C ( d b ). The distances from the O's to the M's in the octahedral vertexes on the side of d a are notated as d a1 and d a2 and a similar notation with d b (Fig. 2b). There are also two different C-O distances in the Fig. 1 Aragonite 2x3x3 unit cell. Oxygen atoms in red, carbon atoms in gray and calcium atoms in green. CO groups are bonded by sticks Fig. 2 a NiAs structure (Nickel atom in green, Arsenide atoms in purple). The Ni atom is located at the center of the trigonal prism. b NiAs-type structure of the orthorhombic carbonates. The anion CO 2− 3 is not located at the center of the distorted trigonal prism. The internal distances are defined inside the octahedron. c Witherite 0 GPa internal distances. d Witherite 6 GPa internal distances slightly non-planar CO 2− 3 group of the aragonite-type structures, which becomes more symmetrical and planar with the increasing of the cation radii of the isostructural orthorhombic carbonates (Antao and Hassan 2009). In opposite to the NiAs structure, the behavior of the structure of aragonitetype minerals with pressure can not only be studied by analyzing the compression of the distorted trigonal prism, but by analyzing the changes of bond distances and atomic volumes. The crystal structure of these compounds can be studied of different structural points of view, such as cations coordinate to nine oxygens coming from the carbonates, as a layered structure by stacking of cations and anions, such as Fig. 1 suggest, however, the NiAs structure includes the least extension with the maximum information, especially from the symmetry/asymmetry point of view, which is very suitable from the computational point of view.
The initial crystal structure of the aragonite, cerussite, strontianite, and witherite minerals were obtained from the American Mineralogist Crystal Structure Database (Downs and Hall-Wallace 2003)

Atomic compressibilities
In order to look for a partition of the bulk modulus, it is necessary to take into account that the crystal volume is the sum of the atomic volumes of the atoms in the unit cell: The compressibility, , of the crystal at constant temperature is defined by: being K the bulk modulus. Introducing Eq. 1 into Eq. 2: From the above, the atomic compressibility is defined as: So, the compressibility of the crystal is decomposed into atomic contributions as it is determined by the expression: Where f i is the fractional volume occupancy of an atom in the unit cell. Moreover, taking into account the inverse relationship between compressibility and incompressibility/bulk modulus, the previous equation can also be decomposed in atomic incompressibility moduli contributions: The above equations together with the QTAIM theory show that the compressibility and the bulk modulus of a crystal can be decomposed in a volume-weighted sum of atomic contributions. Therefore, there are two factors that determines the compressibility, the relative volume that an atom occupies and its local compressibility.

Introduction to QTAIM
The topological analysis of the electronic density (r) according to the QTAIM (Bader 1990) provides an accurate mapping of the chemical concepts of atoms, chemical bonding, and structure. It provides a unique partitioning of the physical space of a molecule or a solid in subsystems where any quantum chemical observable is univocally well defined. Every subsystem (or basin) is bounded by the presence of zero-flux surfaces of the gradient of the electron density that envelop one local maxima in (r) . The definition of an atom inside a molecule emerges when it is considered that the nuclear positions behave topologically as local maxima in (r) . As a result, either a molecule or a solid property associated to a quantum chemical observable can be partitioned into atomic basins contributions using the QTAIM. Therefore, the unit-cell volume can be divided into atomic volumes, allowing such partition for the analysis of the bulk properties of a solid in terms of atomic contributions.

Computational details
Density functional theory calculations were performed with the Quantum-Espresso package, (Giannozzi et al. 2017(Giannozzi et al. , 2009) using the PAW method (Blöchl 1994) in a plane wave (PW) basis set. The presence of an attractive and weak carbon-carbon interactions in this orthorhombic carbonates has been described previously by some authors (Nelyubina and Lyssenko 2012;Vidal and Navas 2014), so the selection of a density functional capable of take into account weak interactions is important. To this end, a combination of the B86b (Becke 1986) exchange functional with PBE correlation (Perdew et al. 1996) was used together with the exchange-hole dipole moment model (Becke and Johnson 2007; to incorporate dispersion effects. Such combination has been proven as a good choice to describe weak interactions . Some calculations were previously done to obtain an optimal PW kinetic energy cutoff, density cutoff and k-point grid, resulting in 90 Ry, Page 5 of 15 13 800 Ry, and 5x3x5, respectively. The electronic densities obtained from the calculations were analyzed by means of the QTAIM theory (Bader 1990), using the CRITIC2 program (de-la Roza et al. 2009Roza et al. , 2014. This program is able to calculate the atomic charge and the atomic volume by integration over the basin of every atom in the QTAIM framework using several algorithms. In this work the integration were performed using the algorithm of Yu and Trinkle (Yu and Trinkle 2011).

Birch-Murnaghan equation of state fitting
The geometries and the unit-cell volumes of the minerals were optimized at several target pressures, which go from 0 GPa to 6 GPa, with pressure increments of 1 GPa. This compression limit looks like suitable to get results in agreement with pressure in the Earth's upper mantle, without to undergo any phase transition. By parametrically changing the target pressure a set of (V,P) values can be obtained; fixing the volume at a place of the pressure a set of (E,V) can be also obtained (Blanco et al. 2004). Any of the two sets can produce equivalent equation of states (Blanco et al. 2004). The equation of state (EoS) was determined by fitting the pressure-volume (P,V) values to a third-order Birch-Murnagan EoS (Murnaghan 1944;Birch 1947 where f v is the Eulerian finite strain, K 0 is the bulk modulus at a 0 GPa pressure, and K ′ is the first derivative of K 0 with respect of pressure. The EosFitCalc (Angel et al. 2021) and EosFit7.6c (Angel et al. 2014) were used for the least square fitting. The cell parameters and distances moduli were fitted with the same program by building a cubic polyhedron, from which the incompressibility linear moduli were calculated.

Crystal structure
Cell parameters of this series at a pressure of 0 GPa are shown in Table 1. Theoretical results are in agreement with the experimental results. If we take into account the average of the deviations of the theoretical values with respect to the experimental values a maximum value of 2.2% is obtained. The b axis of cerussite and the c axis of witherite show the largest absolute average deviations in percentage, 1.5% and 2.2%, respectively. It should be noted that other ab initio calculations at PBE level for witherite are also in agreement with our values (Table 1).
Cell parameters are plotted as a function of the Shannon 9-fold coordination radius of the cation (SIR), (Shannon and Prewit 1969;Shannon 1976) and linear equations are fitted by the least square method (Fig. 3) In Fig. 3a, the experimental values of volume are also included. Most of them are close to the calculated values with the exception of Ono  for the witherite. The volume shows a similar behavior to the cell axes with respect to the SIR (Fig. 3a). However, Li and Liu (Chung-Cheng and Lin-Gun 1997b) showed a curved trend between the molar volume and the SIR. For those authors, the aragonite is the most outsider of the linear trend, doing a curve with respect to the other three carbonates which followed a linear trend.
For the cell axes (Fig. 3b), the linear fittings to the ionic radii yield square correlation coefficients. The a axis looks like to be the less affected by the SIR of the cation, and the b axis shows the largest slope of the three parameters, being approximately three times the slope of the a axis, and c axis present near twice the slope of a axis.
In Table 2, several geometric parameters are compared to the experimental ones coming from different references. Most of our calculated distances are very close to the experimental distances, with a 2% as the largest deviation from the average of experimental values. The C-O distances are also divided into two sets: (i) two O's with a symmetry plane between them, being symmetrical in the CO 3 group; and (ii) one single O being asymmetric in the CO 3 group. They also have small deviations from the known experimental values, being our calculated values slightly larger than the experimental values. In Table 2, we have also added the CC distances between the two closest stacked CO 3 groups, whose stacking is quasi parallel to the c axis, which are the largest of all chosen distances in Table 2; they increase as a function of the cations' SIR: The cerussite decreases the linear trend, and meanwhile, the alkaline earth cations follow a best linear trend.
Vidal and Sánchez-Navas (Vidal and Navas 2014) showed the presence of an attractive interaction between C atoms in cerussite. This CC long-range interaction must be taken into account as well for the understanding of the behavior under pressure of this group of minerals. Table 2 also reports the distances from the carbon atom to the closest and farthest bases of the octahedral polygons, indicated as d a and d b (depicted in Fig. 2b). Both distances and the difference between them, d a − d b , increase with the cation size; the latter is a measure of the asymmetry of the location of the CO 2− 3 groups within the octahedrons and it also reveals the deviation compared to the NiAs structure model (Table 2). Therefore, our calculations rightly describe the crystal structure and the atomic group interactions in concord with the experimental data at room pressure.

Atomic charges and atomic volumes
In    (Shannon 1976;Shannon and Prewit 1969), and our calculated values come from non-spherical volumes but from the atomic basins. These atomic basins are quite spread, and they are limited by the boundaries (QTAIM zero-flux surfaces) of the other surrounding atomic basins. The volumes of the C atoms increase with the volume of the metals, and a linear relationship is also found: although the slope is very small. This fact could indicate the importance of the carbon atom in the system (vide supra). The difference of the distances of the CO 2− 3 to the bases of the octahedron (Table 2) is a quadratic function of the atomic volumes of C: d a − d b (measures the octahedron asymmetry) as a function of the volumes of the CO 2− 3 groups of the alkaline earth carbonates also show a linear trend, but cerussite goes out of the linear trend. Nonetheless, the atomic volumes of carbon atoms yield a good fitting with the cations volume (Eq. 10), meanwhile the oxygens are considered like stuff (O'Keeffe and Hyde 1985) in the crystal structure of carbonate series. From these results, it is also remarkable the CO 2− 3 is not a rigid group as it has generally been considered.

Crystal structure
In general, the increasing pressure yields compression to atomic distances, axes and volumes. The DFT values of the volumes as a function of P are depicted in Fig. 4a along with BM3 EoS and the experimental values from different sources. The volumes and their slopes are quite close to experimental values (Fig. 4a), nonetheless, our values are slightly lower than the experimental ones. From BM3 least square fitting, the bulk moduli of the compounds are obtained, which are given in Table 4 in comparison with experimental values. Table 3 Integrated atomic charges Q ( e − ), atomic volumes V (Å 3 ) and CO 3 volumes of the aragonite-structure group at zero pressure. O 1 is the asymmetric O and O 2 is one of the two symmetric Oxygens of CO 3 group. The cell volumes are calculated with Z=4  Bulk moduli as a function of the volume of the compounds can be studied following to Anderson and Nafe(Anderson and Nafe 1965). They found a dependence of the bulk modulus upon volume of different compounds especially in oxides. The values of the bulk moduli are taken from the elastic constants and they are converted to isotropic values from the averaged method of Reuss-Voigt-Hill. In order to study the variation of the bulk modulus at room pressure as a function of the volume of the compounds by number of atoms in the system, they deduce the following equation: where C is the integration constant and K � = dK dP . This equation is exact considering K � = Constant . The two first values of K' of our series are 4.6, which fulfills 12; however, from strontianite going on, the series K' values decrease. Most of the experimental K ′ values are around 4, which come from the BM2 fitting. Theoretical and experimental values low down as a function of the V 0 (Table 4). Our DFT values go from 4.6 (aragonite) to 1.6 (witherite) (Fig. 4b and Table 4 ), which follow increasing volume. A decreasing linear function is found from the strontianite to witherite, indicating that in this series, when the pressure increases, the sequence of K 0 will go increasing more slowly with the increasing cell volume. The lack of constancy of the BM3 fitting K ′ DFT values can be expanded in Taylor's (see supplementary material) series around the maximal value and taking only the first derivative of the series, and after integration this (12) ln(K) = −K � ln(V) + C expansion yields a modified Anderson and Nafe equation (if it is compared to Eq. 12): The first term of the second member is that coming from the Anderson and Nafe's equation, which is constant, the second one is a correction to the first one and C is a constant adding the integration constants. This can be seen in (Fig. 4b).
DFT results give a linear trend for Eq. 13, and most of the experimental results cluster along this linear trend, although they are calculated using K � 0 = K � DFTSr (strontianite) and V 0 = V DFTSr in Eq. 13. This equation is linearly fitted by the least square method. The integration constant goes from the intercept of the fitting.
When the pressure increases, a axes decrease producing, at 6 GPa, a maximum deformation of 1.55% in strontianite and a minimum deformation of 0.82 in witherite (Fig. 5a). The majority of deformations do not clearly show any linear trend, although aragonite looks like to be the most linear and the witherite is the most bend. For the nonlinear trends, we can at least make two regions out: i) a linear function from 0 GPa to a relatively low pressure value; and ii) the previous trend bends at 2 GPa, with the exception of cerussite that bends at 1 GPa. These two regions can be also seen in the experimental results depicted in Fig. 5a. Results of Palaich et al. (Palaich et al. 2016) and Wang et al. (Wang et al. 2015) for the aragonite and strontianite, respectively, are quite close to ours, especially in the first region; results from Holl Table 4 Calculated and experimental volumes V 0 (Å 3 ), bulk moduli K 0 (GPa) and its first derivative with respect to pressure K ′ 0 of the Birch-Murnagham equation fitting of aragonite-structure carbonates.
SSD is the sum of the square differences between the objective P and the calculated P from the fitting  (Antao and Hassan 2009) et al. (Holl et al. 2000) for witherite are also quite close to our results; however, those of Wang et al. (Wang et al. 2015) for the same mineral are far away from our results having a larger compression than our results. It could be possible that the silicon oil compression media used in their experiment could affect to their results. The second region, in general, shows lower slope than the first region, which looks like to indicate a certain trend to the softening of the axis. Besides, reductions of the compression's slope of the a axis for the strontianite, cerussite and witherite are also observed in our calculations (Fig. 5a). The values of a 0 , K 0a , K ′ a and K ′′ a for aragonite (Table S1) are quite close to the experimental values of Palaich et al. (Palaich et al. 2016). The b axis has a similar behavior than a axis (Fig. 5b). In witherite, however, around 3 GPa, a softening is observed, which goes on to 6 GPa. This softening is also found in the experimental results of Holl et al. (Holl et al. 2000). However, the results of Wang et al. (Wang et al. 2015) do not reproduce the softening behavior and they are far away from Holl et al.'s (Holl et al. 2000) and ours results. The Holl et al.'s (Holl et al. 2000) deformation is approximately parallel to our results. This softening in witherite looks could be the previous physical transformations in the crystal structure to a phase transition from the orthorhombic phase to the trigonal/post-witherite phase (Holl et al. 2000). The K 0b values for aragonite and strontianite are larger than the known experimental values (Table S1). For the witherite no value has been calculated because of the softening from 3 GPa.  (Wang et al. 2015) using silicon oil compression media; Wang15b for strontianite means Wang et al. (2015) (Wang et al. 2015) using methanol and ethanol compression media; Holl00 is ref. (Holl et al. 2000) Finally, the most important compression is for the c axis, which shows a maximal deformation of 9% for witherite, and a minimal for aragonite, 3% (Fig. 5c). Which seems to indicate the compression of the c axis is related with the size of the cations' SIR, which follows a linear behavior. The linear equation for the alkaline earth carbonates is: The cerussite is located slightly out of this linear trend. Compression of the aragonite in reference (Palaich et al. 2016) goes to 35 GPa and the c/c0 values follow a curve trend, our linear functions is the first part of these authors' trend. In general, our theoretical values follow quite close to experimental values with the exception of those of Wang et al. (Wang et al. 2015) for witherite. However, Holl et al's (Holl et al. 2000) results follow quite close our results for witherite. Our K 0c values which go from 154 GPa for the aragonite to 57.9 GPa for witherite (Table S1), giving a linear trend with K 0 : All these compressional results for the c axis are the most meaningful, indicating the structure will compress mainly along the c axis. Thus, the crystallographic direction [001] must be considered the softest direction along the crystal structure, and this is quite close in direction with the CC long-range attractive interaction. This fact was previously described by some other authors for orthorhombic carbonates (Nelyubina and Lyssenko 2012;Vidal and Navas 2014). The fitting of the cell volumes in J∕(mol * bar) units as a function of cell compressibility ( Mbar −1 ) yields a linear function: so, the slope turns out to be in energy units (MJoule/mol). This energy per mol could be considered as the necessary one of the compound to change the cell volume to reach an increase of the compressibility of one Mbar −1 , which could be thought about as a property of the particular series.

Interatomic distances
The CC distances are quasi parallel to the c axis and are clearly affected by compression. The linear moduli of these interatomic/intergroup distances, being their values related to K 0c (Table S1): This shows K 0c depends mainly on the compression of CC distances, considering Eq. 15, where K 0 is related with K 0c , so K 0 is also related with K CC . The variation of the longest distance of the CO 2− 3 group (Fig. 2c, d) to the base of the octahedron ( d a ) as a function of pressure for the orthorhombic carbonates clearly decreases up to 6 GPa. The largest linear modulus corresponds to the aragonite while the smallest modulus to the witherite (Table 5).
K 0da is much higher than the bulk moduli, yielding a linear trend with low correlation. d b 's are shown as a function of the pressure in Fig. 6a. Each mineral shows different slope and a crossing point is found around 4 GPa, indicating that d b at this pressure is independent of the system. d b incompressibility moduli are given in Table 5. These moduli are of the same order than the bulk moduli, which shown a linear trend with the bulk moduli with low correlation.
The lowest linear incompressibility modulus is for witherite, which is near 1.6 times smaller than its bulk modulus (Table 5); however, K db and K ′ db of the aragonite have similar values than the bulk values. Due to the low value of d b incompressibility modulus of witherite, and the largest ionic radius of the Ba 2+ , at a relatively small pressure change, the CO 2− 3 group can be too close to the metal octahedron base, and a structural instability could be arisen, and consequently, a change of phase could be provoked. It is well known, witherite shows an early pressure in the series mineral to change from the orthorhombic phase to the trigonal phases [7 GPa (Holl et al. 2000)]. The difference between d a and d b distances indicates the asymmetry of the CO 2− 3 position in the interstices of the octahedral polygons, in other words, the displacements of the CO 2− 3 groups from the center of the octahedrons. These asymmetries linearly increase with the (17) K 0c = 32 + 0.78K 0CC ;R 2 = 0.998 Table 5 Atomic distances incompressibility moduli (GPa) and their first derivatives (K') with respect to P of the largest C ⋯ Oct distance ( d a ) and the shortest C-Oct distance ( d b ) (See Fig. 2 (Fig. 6b), and, so, a large value can induce a change of phase, and the witherite having the largest slope could be the first compound of the series to change the phase. The pressure of phase transition (PPT) of the alkaline earth carbonates turns out to be a linear function of the slopes of d a − d b (SDAB) in Fig. 6b: using as phase transition pressure 45 GPa for the aragonite (an average between 40 GPa (Palaich et al. 2016) and 50 GPa (Santillán and Williams 2004), 35 Gpa for strontianite (Chung-Cheng and Lin-Gun 1997b) and 7 Gpa for witherite (Holl et al. 2000). The phase transition will be reached when a certain value of the asymmetry parameter or a certain value of the relation (d a − d b )/c is reached. In this meaning, the slope SDAB could be considered as an indicative parameter. The longest distance, d a , has a larger modulus ( K da ) than d b ( K db ) ( Table 5). In principle, such behavior could be thought contradictory, large distances should give small linear moduli, but it should be considered that the distances d ai are the shortest and accordingly they give the strongest OM interactions (see Fig. 2c and d, internal distances for the witherite at 0 and 6 GPa, respectively) and therefore K da is the largest. So, d a and d b are the projections of d ai and d bi , where the MO interactions occur.

Atomic compressibilities
In Table 6, atomic moduli, compressibilities and atomic volume fractions are shown. Atomic moduli of cations decrease along the series but although follows the sequence of the cell bulk moduli a very low correlation is found, being K Pb quite close to K Ba . However, K 0 of the alkaline earth cations The slope indicates the alkaline earth cell bulk moduli varies close to the atomic moduli variation. This out-of-linear-trend of Pb 2+ is brought about in most of our relationships, indicating once again the importance of identical atomic electronic structure in many of the cell properties explanations from the atomic properties. Carbons moduli follow linear relationships with K 0 including the C atom of cerussite in the least square fitting: The linear trend including all minerals of the series indicates that this relationship determines the rationale of the compressibility of the whole series.
The largest atomic oxygen moduli correspond to the asymmetric oxygen, being the symmetric oxygens slightly smaller than the asymmetric ones. The bulk moduli increase as K Oasym increase; however, no right linear relationship is found if the cerussite oxygen is included in the function, but the alkaline earth oxygens give a linear relationship. The K CO 3 follows a similar trend to the previous O.
The gist of this analysis is the weak bond interactions, volumes, distances and atomic moduli of the carbon atoms along the c axis. If cerussite is removed from the series, and we keep only the alkaline earth carbonates, many other magnitudes will look like to determine compressibility, previously pointed out in literature, especially Metal ⋯ Carbonate. Our least square fitting should be considered interpretative because of shortness of series of the orthorhombic carbonates studied, and more systems should be consequently (20) K 0cell = 21.9 + 0.50K C ;R 2 = 0.999 included in the series. Therefore, it would be very interesting to study a broader series including many other orthorhombic carbonates pertaining to different columns of the periodic table. However, the other orthorhombic carbonates are the rutherfordine and widenmannite (Dana 1997) with UO 2+ 2 and Pb 2+ cations, respectively, but their spatial groups are not the same than the carbonates group studied in this work. From these results, it can be highlighted the role of the carbon atom in the CO 2− 3 group, not being a passive role at all as it could be thought at first sight, even more its compressibility increases along the series. The anisotropic character of the aragonite-type structure (and of the NiAs type structure) leads to compression perpendicular to the layers. Diverse layer distances may all well correlate to lattice behavior mathematically. However, the change of the distance between the O 3 C ⋯ CO 3 layer stacking distances with pressure would explain perfectly the compression on the c crystallographic axis because they are quasi-parallel.
If we take into account the alkaline earth atomic volumes in J/mol·bar units at 0 GPa and their compressibilities in Mbar −1 units, a linear correlation can be found: This slope turns out to be in MJ/mol and it is of the same nature of Eq. 16, that is, it is the necessary energy to change the atomic volume of the cation to increase of the compressibility in one Mbar −1 . The atomic volume of C as a function of its atomic compressibility (in the previous units) yields: In this event, the carbon atom of cerussite is included in the fitting, so, the square correlation coefficient lows. The slope for C's is the lowest found in this analysis; the carbon slope turns out to be 57 times lower than the slope of alkaline earth cations.

Conclusions
Theoretical crystal structure of orthorhombic carbonates has been determined with DFT at increasing pressure. The calculated cell parameters and geometries of atomic groups agree with the experimental results. Even the lengthening of the b axis of witherite measured during the experimental compressions from 3.5 GPa is reproduced in our theoretical work, a fact that is consistent with an early change of phase from orthorhombic to a trigonal phase. The bulk moduli also agree with the average of the experimental values. With respect to the axial compression of the orthorhombic carbonates' unit cells, the c crystallographic parameter is by far the most compressible. The C ⋯ C moduli are related to K c , which also is related to K 0 . These relationships indicated the important role of the c axis and the compressibility of atoms and distances along this direction in the compressional behavior of orthorhombic carbonates.
Moreover, the atomic charges are approximately constant along the series and the atomic volume of cations are approximately twice the volume coming from the SIR, indicating the different atomic basins coming from the QTAIM analysis with respect to the SIR.
On the other hand, it has been found that the bulk modulus K 0 is a linear function of the metallic atomic bulk modulus K M , but only for the alkaline earth metals, although the Pb 2+ can also get into the function, but decreasing the correlation coefficient, indicating the importance of the different electronic structure of the cation.
The distances from carbon atom to the faces of octahedron measured quasi parallel to the c axis show different compressibilities and an isocompression point. The difference of distances to the faces of the octahedron gives the asymmetry/corrugation of the CO 3 groups in the octahedral structure, which show different behavior as a function of the pressure. The slopes of the asymmetric distances ( d a − d b ) of the alkaline earth cations could be related with the pressure of phase transition from the orthorhombic phase to trigonal phase.