菱镁矿MgCO₃属于三角晶系，其空间群为$R\bar 3c$，原胞中包含10个原子（2个Mg原子、2个C原子和6个O原子）。为了讨论铁及其在高温高压下的自旋转变对菱镁矿物理性质的影响，本研究选择用1个Fe原子替换1个Mg原子，得到Fe和Mg的摩尔比为1∶1，分子式为（Mg_0.5Fe_0.5）CO₃。

几何结构优化和相关能量的计算^[34]采用基于密度泛函理论（Density of functional theory, DFT）的第一性原理分子动力学计算软件VASP完成。由于结构中有铁存在，需要考虑强关联作用，因此计算中考虑了Hubbard参数U的影响，即LDA+U。已有的研究表明，U值会随自旋转变而改变，本研究选取Tsuchiya等^[18]、Krukau等^[35]通过线性响应理论计算得到的U值，分别为U_LS = 5.3 eV，U_HS = 4.0 eV。K点采用Monkhorst-Pack方法生成以Γ点为中心的15 × 15 × 15网格，截断能设置为1 000 eV^[36]。总能量的收敛阈值设置为1 × 10⁻⁶ eV/cell，原子力的收敛阈值为1 × 10⁻³ eV/cell。采用基于准谐近似的PHONOPY软件计算得到二阶原子间力常数（IFCs）、声子谱和其他热力学参数^[37]。构造求解高对称q点声子频率本征值的动力学矩阵，然后利用倒空间中动力学矩阵的傅里叶变换，计算其他一般q点的声子频率。所有计算中，以Γ为中心的q点网格选取为30 × 30 × 30。计算热力学参数之前，用VASP软件完成不同位移或不同体积下结构的能量计算，计算过程中的参数设置与上述几何优化参数设置相同。

根据已有的研究结果，MS态时的吉布斯自由能可表示为

式中：n为LS态在MS态中所占的百分比，p为压力，T为温度，G_HS、G_LS分别为HS态和LS态的吉布斯自由能，G_mix为HS-LS混合态的吉布斯自由能。G_mix可以表示为

式中：X_Fe为铁在含铁菱镁矿中的摩尔分数，k_B为玻尔兹曼常数。在给定的压力和温度条件下，通过求最小化MS态吉布斯自由能得到LS态的占比n

式中：ΔG_LS−HS为(Mg_0.5Fe_0.5)CO₃在LS态和HS态的吉布斯自由能之差。对于Fe²⁺来说，自旋和轨道简并量子数分别为s = 2，m = 3。得到n后可以推导出MS态的热膨胀系数α(n)和等温体积模量K_T(n)分别为

式中：${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$和${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$为LS态占比的导数。根据式（4）和式（5）导出其他参数的表达式

式中：K_S为绝热体积模量，${C_{V,\rm{m}}}$为定容比热容，ρ为密度。根据上述公式可以计算出MS态下的体积V、速度v和格临爱森常数γ。

含铁菱镁矿(Mg_0.5Fe_0.5)CO₃的LS态与HS态的焓差ΔH（H_HS–H_LS）随压力的变化趋势如图1所示。由图1可知：在低压时，ΔH < 0，HS态的焓值小，因此HS态更稳定；在高压时，ΔH > 0，LS态的焓值变小，LS态变得更稳定。本研究计算得到的自旋转变压力为44.5 GPa，与实验观察得到的结果（40～52 GPa）很好地符合^{[5, 10-11, 29-30]}，与Hsu等^[32]的计算结果（48 GPa）也符合较好，较小的差异源于铁在含铁菱镁矿中的摩尔分数不同。

图2(a)给出了菱镁矿含铁前后的体积对比。由图2(a)可以看出，菱镁矿含铁后HS态对应的体积比LS态对应的体积大，在30 GPa、300 K的温压条件下，(Mg_0.5Fe_0.5)CO₃在HS和LS态下的体积分别为77.138 10 ų和73.366 07 ų，自旋转变后体积减小5.0%左右。对比含铁菱镁矿与不含铁菱镁矿的体积可知，含铁菱镁矿LS态的体积减小，HS态的体积则在低温端增大，在高温端减小，说明含铁菱镁矿的HS态对应的热膨胀系数小于不含铁菱镁矿。本研究中自旋转变导致的体积变化幅度比Liu等^[11]实验测量得到的变化幅度稍小，这是由于实验所用样品为(Mg_0.35Fe_0.65)CO₃，体积变化幅度的差异源于测量样品中铁的摩尔分数不同。

热膨胀系数是研究矿物热力学性质的重要参数，本研究计算了含铁对菱镁矿热膨胀系数的影响。图2(b)分别给出了含铁菱镁矿的HS、LS态在30 GPa和60 GPa时的热膨胀系数，为了便于比较，图2(b)中给出了对应压力下MgCO₃的热膨胀系数以及部分前人的结果^{[11, 32]}。由图2(b)可知，在相同的温压条件下，(Mg_0.5Fe_0.5)CO₃的HS与LS两种自旋态的热膨胀系数均小于MgCO₃，即含铁会导致菱镁矿的热膨胀系数急剧减小。在30 GPa、300 K的温压条件下，计算得到MgCO₃与(Mg_0.5Fe_0.5)CO₃的热膨胀系数分别为8.05 × 10⁻⁶ K⁻¹和3.13 × 10⁻⁶ K⁻¹，减小幅度为61.0%；(Mg_0.5Fe_0.5)CO₃的HS和LS态的热膨胀系数分别为3.13 × 10⁻⁶ K⁻¹和4.22 × 10⁻⁶ K⁻¹，铁的自旋转变导致的热膨胀系数增加幅度为34.8%。此外，热膨胀系数会随着压力的增加而减小，并且对温度的依赖性减弱。这与含铁方镁石在低压下HS态的热膨胀系数比LS态大，而在高压下HS态的热膨胀系数反而比LS态小的变化趋势有所差别^[28]。根据格临爱森定律

由于影响热膨胀系数的主要物理量为格临爱森常数γ、定容比热容C_V,m、体积模量Κ_T及体积V，因此可以从上述热力学参数的变化来解释含铁菱镁矿热膨胀系数的变化机制。

图3给出了(Mg_0.5Fe_0.5)CO₃的热力学参数随温度和压力的变化关系。由图3可知，定容比热容在低温段先随温度升高而急剧增大，随后趋于饱和，MgCO₃的定容比热容为231.75 J/(mol·K)，(Mg_0.5Fe_0.5)CO₃在HS、LS态时的定容比热容分别为232.99 J/(mol·K)和231.58 J/(mol·K)，彼此之间的差异很小，因此定容比热容对热膨胀系数的影响可以忽略不计。由1.2节的研究结果可知，菱镁矿含铁及铁的自旋转变引起的体积变化在5.0%左右，无法成为热膨胀系数剧烈变化的决定性因素。由图3(b)和图3(c)可知，菱镁矿含铁后格临爱森常数明显减小，由HS态转变为LS态后明显增大，与热膨胀系数的变化趋势很好地符合。需要注意的是，含铁菱镁矿LS态的格临爱森常数γ比不含铁时大，而热膨胀系数却比不含铁时小，由此可以推断，体变模量K_T对热膨胀系数的变化也起到至关重要的作用。因此，菱镁矿含铁及铁的自旋转变导致的热膨胀系数变化是由格临爱森常数γ与体变模量K_T共同决定的。

Liu等^[11]、Hsu等^[32]、Fu等^[33] 的研究表明，含铁菱镁矿的MS态在高温高压下会引起弹性及声速的异常变化。研究MS态热力学参数的前提是计算出自旋共存时LS态的占比n。根据式（1）～式（8）计算出MS态下LS的占比n，n对p和T的一阶导数${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$、${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$，以及α、V、v等热力学参数。

图4给出了不同温度下n随压力的变化关系。由图4可知：在低温时，n急剧增大到1，即HS态与LS态共存的压力区间很窄；随着温度升高，HS态与LS态共存的压力区间明显增大；当温度为300 K时，HS态与LS态共存的压力区间约为5 GPa；温度为1 500 K时，两者共存的压力区间增大到约15 GPa，与Liu等^[11]、Hsu等^[32]通过实验和理论计算得到的趋势一致。同时，通过计算发现：自旋转变的压力随温度的升高而增加；温度为300 K时，自旋转变发生在约40 GPa；而温度为3 000 K时，自旋转变发生在约60 GPa。该自旋转变行为在铁方镁石的理论计算和实验研究中也有报道^{[18, 20, 22, 38-40]}。HS态与LS态共存区间和自旋转变随压力的变化关系决定了n对温度和压力的一阶导数，如图5所示。计算表明：当压力为50 GPa时，${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$异常变化的起始温度为500 K左右；70 GPa时，${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$异常变化的起始温度为1 500 K左右。同时${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$峰的宽度也由50 GPa时的1 000 K增大到70 GPa时的1 700 K，由此可知，自旋转变引起的${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$峰随着压强的增加向高温方向移动，同时峰的宽度也随之增加。图5(b)给出了${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$随压强的变化关系，其峰高和峰宽随压力的变化趋势与${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$类似。

因为${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$与${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$可以不等于零，由式（4）和式（5）可知，等温体弹模量K_T(n)和热膨胀系数α(n)受${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$和${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$的影响。本研究计算了(Mg_0.5Fe_0.5)CO₃在MS态下的体积、热膨胀系数和声速等热力学参数，并与MgCO₃、(Mg_0.5Fe_0.5)CO₃的HS态的参数及前人的研究结果进行了对比。首先是体积的变化，由于$V(n) = {\left. {\dfrac{{\partial G(n)}}{{\partial p}}} \right|_T} = {\left. {\dfrac{{\partial G(n)}}{{\partial p}}} \right|_{T,n}} + {\left. {\dfrac{{\partial G(n)}}{{\partial n}}} \right|_{T,p}}{\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$，在平衡时，${\left. {\dfrac{{\partial G(n)}}{{\partial n}}} \right|_{T,p}} = 0$，因此，(Mg_0.5Fe_0.5)CO₃在MS态的体积为V(n) = nV_LS(p,T) + (1－n)V_HS(p,T)，其中V_LS和V_HS分别为(Mg_0.5Fe_0.5)CO₃在LS态与HS态时的体积。由于在相同压力下，(Mg_0.5Fe_0.5)CO₃的HS态体积比LS态体积大，因此，由图6(a)可以看出，在压力转变区间内，MS态的体积比HS态体积小，与前人的理论和实验结果是相符的。

图6(b)和图6(c)分别给出了热膨胀系数α和速度v随压力的变化关系。由MS态时热膨胀系数的计算公式(4)可知，α主要受到LS态的占比n对温度的一阶导数${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$的影响，其中α与${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$成反比，由图5(a)可知，随着HS态向LS态发生自旋转变，在自旋共存区间内，${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$的变化趋势为先突然减小然后增大，因此热膨胀系数α的变化趋势为先突然增大然后减小，这与体积的变化趋势不同；由式（5）、式（6）、式（8）可知，MS态的速度主要由LS态的占比n对压强的一阶导数${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$决定，且v与${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$成反比，由图5(b)可知，在自旋共存区间内，${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$的变化趋势为先突然增大后减小，因此速度v在自旋共存区间的变化趋势为先突然减小然后增大。在300 K时，MS态(Mg_0.5Fe_0.5)CO₃在约40 GPa时产生异常变化，而在1 500 K时，这类异常变化出现在约50 GPa，该变化趋势与${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$和${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$的变化一致，表明热力学参数的变化受n的影响。α、V、v由于自旋转变引起的异常随温度的增大向高压方向移动，并且变得更加平滑，但在2 500 K左右时，α、V、v的异常仍然比较明显。本研究计算结果与前人的实验研究均符合较好，在相同温度下，仅峰的位置存在较小的偏差，这是由于实验中Fe²⁺的浓度较高，其自旋转变的压力也随之相应地提高。(Mg_0.5Fe_0.5)CO₃的MS态体积和速度的异常变化也会影响其在地球深部条件下的热传导特征。

采用第一性原理计算方法研究了含铁菱镁矿(Mg_0.5Fe_0.5)CO₃中铁的掺入以及自旋转变对其热力学性质的影响，得出以下主要结论。

（1）含铁后菱镁矿晶体的体积发生明显变化，HS态在低温端体积增大，而在高温端体积比不含铁时减小，LS态的体积明显减小，表明铁的自旋转变会导致含铁菱镁矿晶体体积减小。

（2）含铁后菱镁矿的热膨胀系数显著减小，但是铁的自旋转变导致热膨胀系数有所增大，但仍然比不含铁时要小；格临爱森常数与体变模量的变化是热膨胀系数变化的主要决定因素。

（3）含铁菱镁矿的HS态与LS态两种自旋共存时，${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$、${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$在自旋共存区间内不等于零，并在自旋过渡区域出现异常变化，此特征会导致含铁菱镁矿的体积在自旋共存压力区间突然减小至接近LS态时的体积；而热膨胀系数、速度也因受${\left. {\dfrac{{\partial n}}{{\partial T}}} \right|_p}$、${\left. {\dfrac{{\partial n}}{{\partial p}}} \right|_T}$的影响，分别在自旋转变过程中出现异常增大峰和异常减小峰，并且这些异常变化峰随着温度的升高向高压方向移动。

若要进一步明确含铁菱镁矿的热力学性质对地幔在碳循环中的作用，还需全面考虑地幔中其他矿物对菱镁矿的影响，如铁方镁石和钙钛矿等在地球深部同时存在时所表现出的物理化学性质。