To express the elastic properties of a mineral, the bulk modulus or compressibility is the mostly used parameter, however, for rare minerals it is not well known. During the search for another parameter, the question raised, which structural units determine the stiffness of a mineral. For minerals containing an oxygen sublattice, its packing fraction is likely to be a good proxy for the elastic properties since it dominates the occupancy of the unit cell volume. The oxygen packing fraction (OPF) is a parameter that has the great advantage that it is easily available for poor investigated minerals because OPF can be defined as

OPF = nO *V*O/(*V*UC), 1)

where nO is the number of oxygens in the unit cell, VO is the volume of the oxygen anion calculated from the data of Shannon (1976), and VUC is unit cell volume. The correlation between OPF and bulk modulus for different minerals is shown in Fig. 1 demonstrating that OPF increases from 51 to 71 % during a bulk-modulus increase from 300 to 2500 kbar. The data are from the compilation of Holland and Powell (2011) and range from tridymite with low bulk modulus over all mineral groups that contain oxygen to stishovite with a high bulk modulus (including 95 minerals).

Δ*H*mix is defined as the difference between the enthalpy of a solid solution (*H*AB) at a given composition and the enthalpy of its mechanical mixture (*X*A *H*A + *X*B *H*B), i.e.,

Δ*H*mix = *H*AB – (*X*A *H*A + *X*B *H*B), 2)

where *X*A, *X*B and *H*A, *H*B represent the mole fractions and the enthalpy of the A and B component, respectively. Δ*H*mix of different substitutions in different minerals was investigated by DFT using the single defect method (Sluiter and Kawazoe 2002). The data were then described by a symmetric Margules model,

Δ*H*mix = *X*A *X*B *w*H, 3)

where *w*H is the interaction energy between elements on a single crystallographic site (micro *w*H). This parameter is listed in Tab. 1 (given in the supplementary material) and plotted against OPF of the Mg-end members in Fig. 2.

The Mg-Al substitution was investigated in 16 minerals (ordered by increasing *w*H): melilite, chlorite, biotite, cordierite, amphibole, talc, pseudobrookite, pyroxene, sulphate (4 x hydrated), olivine, wadsleyite, sulphate (1 x hydrated), ilmenite, Mg-wolframite, anhydrous sulphate, and ringwoodite that has the spinel structure. As can be seen from Fig. 2, a linear positive correlation of *w*H versus OPF was found (fit parameters are given in the supplementary material). The Mg-Al substitution in biotite, as an example, has an energetic effect that is less than 1/2 of that in pyroxene, and less than 1/6 of that in spinel. This behaviour is expected to be a consequence of the different oxygen packing fractions, which are 53.5 % in biotite, 59.1 % in pyroxene, and 65 % in spinel. In addition to the Mg-Al substitution, the Si-Al substitution has been investigated in 8 minerals. Here, two trends were observed. One, where the tetrahedra are connected with each other (chain and sheet silicates) investigated in 4 minerals (ordered by increasing *w*H): biotite, amphibole, pyroxene, perovskite. The second trend is from island and double island silicates defined by 4 minerals, i.e., melilite, olivine, wadsleyite, ringwoodite. This group is characterised by a much flatter positive correlation. MgSiO3 perovskite, with *w*HSiAl/OPF = 577.4/77.1, formed an outlier and was not included in the fit. The interaction energies of the Mg-Ca substitution were investigated in 7 minerals (ordered by increasing OPF): MgO-CaO solid solution, amphibole, olivine, pyroxene, carbonate, garnet, and perovskite). This substitution is characterised by a flat *w*H/OPF correlation, whereas the Mg-Ti4+ substitution shows a steep slope that is defined by 4 minerals (ordered by increasing *w*H): biotite, pyroxene, olivine, ringwoodite. Finally, 7 minerals were studied for the Mg-Fe substitution finding ideal to slightly negative D*H*mix behaviour (ordered by increasing OPF): biotite, brucite, pyroxene, olivine, spinel, garnet, and perovskite.

Using the interaction energies of the separate Mg-Al and Si-Al substitutions (micro *w*´s, i.e., *w*HMgAl and *w*HSiAl) in combination with the underlying cross-site interaction (*w*Hcross), the macroscopic interaction energy for the Tschermak´s substitution (*W*HTS) can be calculated (e.g., Powell et al., 2014):

*W*HTS = ¼ (4 *w*HMgAl + *w*HSiAl – 2 *w*Hcross), 4)

where the cross-site interaction is given by the energies of the reciprocal reaction:

*w*Hcross = *H*(MgMSiT) + *H*(AlMAlT) – *H*(MgMAlT) – *H*(AlMSiT). 5)

Such calculations were performed for the Mg-Al biotite in order to verify the DFT results of charged unit cells. Using eq (4), *w*HMgAl = 82.5 kJ/mol, *w*HSiAl = 95.6 kJ/mol, and *w*Hcross of 175.1 kJ/mol, the macroscopic interaction energy of *W*HTS = 18.8 kJ/mol is obtained. It represents the interaction energy with ordering of Al on M1 and T1 sites, but without strict local-charge balance (without short-range ordering). To compare it with DFT values obtained without using charged unit cells, a Mg-biotite (KMg3[(OH)2AlSi3O10], phlogopite), a Tschermak substituted biotite (KMg2Al[(OH)2Al2Si2O10], eastonite) and intermediate solid-solution biotites were investigated. Using 10 different super cells for the intermediate biotites (50:50 composition with randomly distributed Al on M1 and T1 with sizes up to 16 formula units), a *W*HTS of 19.0 +/- 3.0 kJ/mol was obtained, which is in good agreement with the *W*HTS calculated via micro *w*´s (18.8 kJ/mol). This value is slightly lower than an experimentally derived value of 25.4 kJ/mol (Dachs and Benisek 2019) because of different ordering schemes. The same calculations were performed for the diopside – CaTs solid solution, (*w*HMgAl = 213.6 kJ/mol, *w*HSiAl = 285.6 kJ/mol, and a *w*Hcross = 422.0 kJ/mol) yielding a *W*HTS of 74.0 kJ/mol. The uncharged cells delivered 74.8 +/- 3.7 kJ/mol, also in good agreement with *W*HTS obtained via micro *w*´s. This value is larger than the experimentally derived value of ~ 24 kJ/mol (Benisek et al. 2007), if a symmetric fit is applied to their data. The difference comes mainly from the fact that the experimental data were derived using a partly disordered CaTs end-member (Benisek and Dachs 2020).