二维自旋转变模型中自旋转变温度的分子动力学模拟

Molecular Dynamics Simulation of the Spin Transition Temperature in a Two Dimensional Spin Transition Model

摘要

在自旋转变材料中，分子大小随自旋态的改变而改变，从而形成晶格畸变。受晶格畸变影响的分子间相互作用会导致自旋态的协同行为，进而改变温度可以发生高低自旋态之间的突变。我们考虑了一个包括分子间相互作用和分子内相互作用的二维自旋转变模型，基于Nosé-Hoover的分子动力学方法研究了分子内势能参数对于高低自旋态转变温度及高低自旋态平均占据数的影响。计算结果表明，在温滞范围内自旋转变温度会随着分子内高低自旋态的耦合强度和分子内高低自旋态的能量差的增大而增大，而随着自旋态简并度比值的增大而减小；温滞宽度会随着高低自旋态的耦合强度和分子内高低自旋态的能量差的增大而减小，而随着自旋态简并度比值的增大而增大。如果要进一步扩大自旋转变的温滞范围，分子内高低自旋态的耦合强度和分子内高低自旋态的能量差有一个下限值，而高低自旋态简并度比值有一个上限值。

Abstract

In the spin-transition materials, the size of each molecule changes with the spin state, which cause lattice distortion. The intermolecular interaction affected by the lattice distortion will lead to the synergistic behavior of spin states, in turn, cause the abrupt transition between the low spin state and the high spin state. Here we consider a two-dimensional spin-transition model incorporating the intermolecular and the intramolecular interactions and investigate the effect of the intramolecular interaction parameters on the spin-transition temperature and the average occupation number based on the Nosé-Hoover molecular dynamics method. Numerical results show that during the hysteresis the spin-transition temperature increases with the increase of the intramolecular coupling parameter and the energy detuning between the high spin state and the low spin state, yet decreases with the increase of the ratio of the degeneracies between the high spin state and the low spin state. The width of the hysteresis decreases with the increase of the intramolecular coupling parameter and the energy detuning, yet increase with the increase of the ratio of the degeneracies between the high spin state and the low spin state. There is a lower limit for the intramolecular coupling parameter and the energy detuning and an upper limit for the ratio of the degeneracies between the high spin state and the low spin state. when the hysteresis can exist.

Graphical abstract

关键词

自旋转变 / 分子内相互作用 / 温滞效应 / Nosé-Hoover方法

Key words

spin transition / intramolecular interaction / temperature hysteresis effect / Nosé-Hoover method

0 引言

在过渡金属化合物形成的自旋转变材料中，当核外电子配对能大于晶体场劈裂能，过渡金属离子会呈现高自旋（High Spin， HS）态；反之，当电子配对能小于晶体场劈裂能时，过渡金属离子会呈现低自旋（Low Spin， LS）态。自旋转变材料受到光照^［1-3］、温度^［4-10］、压力^［11-14］、配体^［15-18］等外界因素调制时，会展现出HS态和LS态的相互转化。伴随着高低自旋态的转变，该材料还会发生如导电性、磁性、发光性质、介电常数及机械性能等物理化学性质的转变。自旋态转变及伴随性质的转变作为热力学相变的一类受到理论和实验的广泛关注。

早期，Cambi等在研究过渡金属离子时实验观察到离子的磁矩会随着温度的变化而产生变化^［19］。Spiering等提出了包括非简谐效应的自旋转变分子系统的自由能，并解释了自旋转变温度对压强的依赖作用^［20］。Nishino等基于分子动力学的理论方法研究了由分子间相互作用所引起的自旋态间协同作用现象，通过温度调控，发现热滞回线的宽度会随分子间相互作用强度的增加而增加^［21］。

近来，Lorenc等实验上采用超短光脉冲诱导Fe（Ⅲ）化合物发生自旋转变并进行光学时域分辨研究，给出了各种不同时间尺度的转换动力学^［22］。Paez-Espejo等同时考虑自旋转变时体积变化和弹性阻挫的自洽弹性相互作用，通过调制弹性阻挫强度，自旋转变材料中可以出现单步或多步转变现象^［23］。Cruddas等在笼目结构中采用弹性模型预测了自旋转变材料中存在新相^［24］。Ndiaye等从弹性阻挫的角度分析了自旋转变材料中加热和冷却时产生的不对称热滞回线^［25］。Cruddas等在自旋转变材料中证明了分子间的相互作用可以产生三种不同的自旋态相，通过温度调节熵与焓之间贡献实现不同相之间的转变^［26］。Affes等采用蒙特卡罗的方法研究了单晶自旋转变材料的热和结构性质，特别考虑了掺杂浓度效应，结果表明随掺杂浓度变化，自旋转变会从一阶突变转化为二阶连续转变^［27］。Chen等实验上通过改变\mathrm{F}{\mathrm{e}}_{\mathrm{4}}\mathrm{C}{\mathrm{o}}_{\mathrm{4}}配体的方法实现了自旋态转变的调控^［28］。

之前的理论研究主要集中在调制分子间相互作用参数（分子间参数）^［21-26］，采用方法包括Nishino等运用的分子动力学方法^［21］和蒙特卡罗模拟方法^［27］。近年来人们理论上通过改变掺杂浓度的方法^［27］或者实验上通过改变配体的方法^［28］来研究自旋态转变的调控。而改变掺杂浓度或者改变配体都会显著影响分子内相互作用，进而影响自旋态转变。因此，研究分子内相互作用参数（分子内参数）的调控对于自旋态转变的研究有重要的科学意义。在本文中，我们主要基于分子动力学方法研究分子内势能参数的调控对自旋转变温度的影响，特别是研究分子内势能参数如何影响温滞行为。

1 理论模型与方法

在过渡金属化合物形成的自旋转变材料中，由于过渡金属的电子配对能与晶体场相互竞争，过渡金属会呈现高低自旋态。对应不同的自旋状态，由于电子振动耦合，即电子畸变使得过渡金属与周围的配体之间的距离会不同。把过渡金属与配体整体简化为一个球形分子，这样分子大小会随着自旋态的不同而变化。同时，这种畸变会引起邻近分子间的自旋相互作用，进而这种畸变的相互作用会导致自旋态间的长程有效相互作用。这里主要考虑由于自旋态不同导致的分子大小不同所产生的晶格畸变。本文考虑了一个16×16的简单正方形晶格模型^［21］。在过渡金属化合物中，如典型的过渡金属离子\mathrm{C}{\mathrm{o}}^{\mathrm{3}+}\left({\mathrm{d}}^{\mathrm{6}}\right)和\mathrm{F}{\mathrm{e}}^{\mathrm{2}+}\left({\mathrm{d}}^{\mathrm{6}}\right)，其HS态对应电子构型为{\mathrm{t}}_{\mathrm{2}\mathrm{g}}^{\mathrm{4}}{\mathrm{e}}_{\mathrm{g}}^{\mathrm{2}}\left(S=\mathrm{2}\right)，而LS态对应电子构型为{\mathrm{t}}_{\mathrm{2}\mathrm{g}}^{\mathrm{6}}\left(S=\mathrm{0}\right)（如图1（a）和图1（b））。由于HS态分子半径大于LS态分子半径，图1（c）与1（d）分别给出了分子均处于LS态和HS态对应二维晶格的示意图。

我们基于Nishino等的工作考虑了分子内相互作用和分子间相互作用 ^［21］。假设分子仅具有一个对称振动模式的主活跃模式（且对于该振动模式分子有各向同性的体膨胀系数），关于依赖于分子大小的分子内势能，采用双阱谐振子势V_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}\left(x\right)，

V_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}\left(x\right)=\frac{A}{\mathrm{2}}\{d+b{(c-x)}^{\mathrm{2}}+a{x}^{\mathrm{2}}-\sqrt[]{\mathrm{4}{J}^{\mathrm{2}}+[d+b{(c-x)}^{\mathrm{2}}-a{x}^{\mathrm{2}}{]}^{\mathrm{2}}}\}\mathrm{。}（1）

这里x=r_i-r_{\mathrm{L}\mathrm{S}}是分子半径对理想的LS态半径的偏离，r_{\mathrm{L}\mathrm{S}}为LS态分子半径，r_i为分子半径大小。该双阱谐振子势V_i^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}\left(x\right)可看作是由LS态谐振子势V(x)=a{x}^{\mathrm{2}}与HS态谐振子势V(x)=d+b{(c-x)}^{\mathrm{2}}经由非对角项J耦合形成的最低势能。a，b分别代表分子处于LS态和HS态附近时的耦合强度，d可以调节高低自旋态势能极小值之差，A可以整体调节分子内势能强度。为了简化，选取x=\mathrm{0}时作为分子的LS态，取x=\mathrm{1}\left(c=\mathrm{1}\right)时作为分子的HS态（如图1（a）和1（b））。图2为单个分子内势能随x的变化曲线，其中参数分别为A=\mathrm{10}，J=\mathrm{0.05}，a=\mathrm{10}，b=\mathrm{0.1}，c=\mathrm{1.0}，d=\mathrm{0.1}，且采用无量纲化处理。之所以采用上述取值一方面是因为可以保证势能函数存在两个极小值，分别对应于高低自旋态，另一方面可以保证高低自旋态有足够高的简并度比值。对于第i个分子，与其半径r_i相对应的共轭动量记为p_i，其约化质量为m。因此N个分子的总分子内动能和分子内势能构成的哈密顿量为

H_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}={\displaystyle \sum _{i=\mathrm{1}}^N}\frac{{p_i}^{\mathrm{2}}}{\mathrm{2}m}+{\displaystyle \sum _{i=\mathrm{1}}^N}{V_i}^{intra}(x)。（2）

除了分子内相互作用，还考虑了分子间的相互作用，包括分子间最近邻和次近邻相互作用，对于最近邻的两个分子i和j，其相互作用记为V_{ij}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}({\mathit{R}}_i,{\mathit{R}}_j,r_i,r_j)=f\left(d_{ij}\right)，d_{ij}=|{\mathit{R}}_i-{\mathit{R}}_j|-(r_i+r_j)，{\mathit{R}}_i和{\mathit{R}}_j分别是分子i和j的二维坐标，r_i和r_j分别是分子i和j的半径。选取函数f为

f(d_{ij})=D({\mathrm{e}}^{a^{\text{'}}(d_{ij}-u_{\mathrm{0}})}+{\mathrm{e}}^{-b^{\text{'}}(d_{ij}-u_{\mathrm{0}})})，（3）

D为分子间相互作用强度常数。u_{\mathrm{0}}是一个常数，选取u_{\mathrm{0}}为对应d_{ij}=\mathrm{0}时使f\left(d_{ij}\right)取极小值时对应的值，a^{\mathit{\text{'}}}，b^{\mathit{\text{'}}}均为常数，选取a^{\mathit{\text{'}}}=\mathrm{0.5}和b^{\mathit{\text{'}}}=\mathrm{1.0}。

为了保持二维结构稳定，还需考虑分子之间的次近邻相互作用，与最近邻相互作用类似，对于次近邻的两个分子i和k，其相互作用记为V_{ik}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}({\mathit{R}}_i,{\mathit{R}}_k,r_i,r_k)=f\left(d_{ik}-\mathrm{\Delta }r\right)，d_{ik}=|{\mathit{R}}_i-{\mathit{R}}_k|-(r_i+r_k)，相互作用函数f与最近邻相互作用有相同的表达式，但不同于最近邻相互作用，在次近邻作用中，选取a^{\mathit{\text{'}}}=\mathrm{0.1}和b^{\mathit{\text{'}}}=\mathrm{0.2}。这样可以保证次近邻相互作用远小于最近邻相互作用。次近邻相互作用的分子由于不接触，因而它们之间有一定间隙，假定间隙\mathrm{\Delta }r=(\sqrt[]{\mathrm{2}}-\mathrm{1})(r_{\mathrm{L}\mathrm{S}}+r_{\mathrm{H}\mathrm{S}})时，次近邻相互作用是最稳定的。实际上，这个间隙\mathrm{\Delta }r会随着温度而改变，但后面数值模拟忽略了这个效应。对于每个分子，其质心坐标为{\mathit{R}}_i=\left(X_i,Y_i\right)，与其对应的共轭动量记为{\mathit{P}}_i=\left(P_{X_i},P_{Y_i}\right)，分子质量为M，则N个分子总质心动能和分子间势能构成的哈密顿量为

\begin{array}{c}{H}_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}={\displaystyle \sum _{i=\mathrm{1}}^N}\frac{{{\mathit{P}}_i}^{\mathrm{2}}}{\mathrm{2}M}+\\ \frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^N}\left(\sum _j{V}_{ij}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}(d_{ij})+\sum _k{V}_{ik}^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}(d_{ik})\right),\end{array}（4）

因此，N个分子构成的总系统的哈密顿量为H_{\mathrm{s}\mathrm{y}\mathrm{s}}=H_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}+H_{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}。

为了研究N个分子的系统其自旋态对温度的依赖性，对于给定温度T采用Nosé-Hoover方法^［29-30］来产生正则系综。为此引入热库的哈密顿量为

H_{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}}=\frac{{P_s}^{\mathrm{2}}}{\mathrm{2}Q}+\mathrm{3}N{k}_{\mathrm{B}}T\mathrm{l}\mathrm{n}s，（5）

其中s称为标度因子，{\mathit{P}}_s为s共轭动量，Q为与s有关的有效质量，N是总粒子数。热库的势函数选择为\mathrm{3}N{k}_{\mathrm{B}}T\mathrm{l}\mathrm{n}s是为了保证该系统能够产生正则系综。系数选择\mathrm{3}N是整个系统包括描述质心的两个自由度和分子半径的一个自由度，此外它还与动力学演化方程是实时还是虚时的具体表达方式有关（实时选\mathrm{3}N，虚时选\mathrm{3}N+\mathrm{1}）。为了后面方程的简化，引入\mathrm{\zeta }=\frac{{\mathit{P}}_s}{Q}。

于是，包括热库的系统总哈密顿量为H_{\mathrm{t}\mathrm{o}\mathrm{t}}=H_{\mathrm{s}\mathrm{y}\mathrm{s}}+H_{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}}。基于Nosé-Hoover方法可以给出系统的演化动力学方程组为

\frac{\mathrm{d}{r}_i}{\mathrm{d}t}=\frac{p_i}{m}，（6）

\frac{\mathrm{d}{p}_i}{\mathrm{d}t}=-\frac{\partial V^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}}}{\partial r_i}-\frac{\partial V^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}}{\partial r_i}-\mathrm{\zeta }{p}_i，（7）

\frac{\mathrm{d}{\mathit{R}}_i}{\mathrm{d}t}=\frac{{\mathit{P}}_i}{M}，（8）

\frac{\mathrm{d}{\mathit{P}}_i}{\mathrm{d}t}=-\frac{\partial V^{\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}}}{\partial {\mathit{R}}_i}-\mathrm{\zeta }{\mathit{P}}_i，（9）

\frac{\mathrm{d}s}{\mathrm{d}t}=s\mathrm{\zeta }，（10）

\frac{\mathrm{d}\mathrm{\zeta }}{\mathrm{d}t}=\frac{\mathrm{1}}{Q}\left({\displaystyle \sum _{i=\mathrm{1}}^N}\frac{{p_i}^{\mathrm{2}}}{m}+{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{{{\mathit{P}}_i}^{\mathrm{2}}}{M}-\mathrm{3}N{k}_{\mathrm{B}}T\right)。（11）

上述方程组中，式（6）和式（7）描述了分子半径和与其共轭动量随时间的演化，分子半径只与其对应的共轭动量有关，而其共轭动量不仅依赖于分子内势能，还依赖于分子间相互作用。类似的，式（8）和式（9）描述分子质心坐标和与其共轭动量随时间的演化，分子质心坐标只与其对应的共轭动量有关，而其共轭动量依赖于分子间相互作用。式（10）和式（11）描述标度因子和与其共轭动量随时间的演化，标度因子与自身和相对应的共轭动量有关，而其共轭动量依赖于与分子半径相对应的共轭动量和质心坐标所对应的共轭动量，且与系统设定温度及系统总自由度有关。

2 数值结果与分析

在数值计算中，我们采用分子半径x_i=r_i-r_{\mathrm{L}\mathrm{S}}来表征分子的自旋态，考虑了一个16×16个分子组成的二维系统。选取开放边界条件，首先对此二维系统加热，使系统温度从0.1上升到2.0，然后再从2.0降温到0.1，每次温度变化步长为0.1。在每一个温度做了60 000步分子动力学时间演化，舍弃前面40 000步暂态演化，之后的20 000步用来计算x_i的平均值，然后再计算N=\mathrm{256}个分子的平均值来表征系统的自旋态，此处用〈x〉代表分子半径偏离理想LS态的系综平均值，由于各态遍历假说，实际的平均包括对每个分子的时间平均和对所有分子的平均。分子动力学模拟采用的是自编的Fortran程序，时间演化的步长选取\mathrm{\Delta }t=\mathrm{0.01}，采用四阶龙格-库塔法来实现系统动力学演化，该方法可以精确到{\left(\mathrm{\Delta }t\right)}^{\mathrm{4}}。为了计算方便，选取m=\mathrm{1.0}，M=\mathrm{1.0}，Q=\mathrm{1.0}，k_B=\mathrm{1.0}。其他分子内势能参数选择为A=\mathrm{10}，J=\mathrm{0.05}，a=\mathrm{10}，b=\mathrm{0.1}，d=\mathrm{0.1}。涉及能量的参数，选取\mathrm{25}\mathrm{m}\mathrm{e}\mathrm{V}作为基准对其做无量纲化处理。

首先检验了不同的分子间相互作用强度下，表征分子自旋态〈x〉随温度的变化情况（如图3）。从图中可以看出，在较小的分子间相互作用下（如图3（a）），系统自旋态随着温度的升高从LS态逐渐转化为HS态，而随着温度的下降系统又从HS态逐渐回到LS态，转变过程是连续的。随着分子间相互作用的增强（如图3（b）），系统自旋态随着温度的升高从LS态突变为HS态，而随着温度的下降系统又从HS态突变回LS态，但是突变温度不同于升温情况。把从LS突变为HS的温度记为T_{\uparrow }，把从HS态突变为LS态的温度记为T_{\downarrow }，温度滞回宽度记为T_{\uparrow }-T_{\downarrow }。进一步，可以发现随着分子间相互作用的继续增大（如图3（c）），温度滞回宽度会进一步增大。这些结论与已有结论是一致的^［21］。关于分子间相互作用，真实实验上可以采用核磁共振、X射线等方法来测量分子间相互作用强度，为了调节其大小，一方面可以更换配体，调节晶体场大小；另一方面也可以通过外加压强调节。

其次，我们分析分子半径、质心动量在不同温度下的分布情况，并计算转变温度的尺寸效应。由图3可以看到在高温时，〈x〉约为0.7。为此计算了在D=\mathrm{35}，温度分别为T=\mathrm{0.5}和T=\mathrm{1.9}时，系统在稳定演化时，截取了100步时间步长，统计了256个粒子，对25 600演化半径做了统计分析（如图4（a）和4（b））。由图4（b）可以看出高温时系统绝大部分电子自旋都呈现HS态，而还有少部分自旋处于LS态，即粒子自旋不会全部发生转化。进一步截取了T=\mathrm{1.9}和时间演化到59 900步时粒子的位置构型图4（c），可以看出，高温时仍有部分分子处于LS态。此时也截取了100步时间步长，对25 600个质心演化的二维动量进行了统计（如图4（d）和4（e）），其分布近似为高斯分布；关于二维动量的空间分布，其每个位置的分布应都与这个分布相似。为了对比，增加了低温下T=\mathrm{0.5}时截取的25 600个质心演化的二维动量分布（如图4（d）和4（e）），由图可以看到，温度越高，分布越宽。进一步考虑了系统的尺寸效应，构建了\mathrm{8}\times \mathrm{8}、\mathrm{12}\times \mathrm{12}、\mathrm{16}\times \mathrm{16}、\mathrm{20}\times \mathrm{20}、\mathrm{24}\times \mathrm{24}个分子组成的二维系统，计算了转变温度随粒子数的变化情况（如图4（f）），由图可知，当粒子数较少时，系统存在突变行为，但没有回线；而当粒子数增加时，系统的高低自旋态转变会出现温滞回线；且粒子数增加时T_{\uparrow }将会增大，而T_{\downarrow }将会下降。由此可见，尺寸效应对于温滞现象有一定影响。

然后，为了进一步考虑自旋态随温度变化情况，特别是温度升降时高低自旋态突变温度及滞回宽度对分子内相互作用参数的依赖情况。我们考虑了分子自旋态随分子内高低自旋态的耦合强度J的变化情况（如图5）。由图5（a）知，随着分子内高低自旋态的耦合强度J的增大，从LS突变为HS的温度T_{\uparrow }从有限值缓慢增加，而从HS态突变为LS态的温度T_{\downarrow }只有当耦合强度J达到一定值时才会存在，且该温度随J的增大而增大。这样分子自旋态随分子内高低自旋态的耦合强度J变化图中可以分成三个区域（如图5（b）—图5（d））。当高低自旋态的耦合强度较小时，系统随着温度的上升可以从LS态转变为HS态，而随着温度的下降，系统在研究的温度范围内没有观测到自旋态的转变，正如图5（b）中所示，从LS突变为HS的温度T_{\uparrow }存在，而从HS态突变为LS态的温度T_{\downarrow }不存在。当高低自旋态的耦合强度J增大时，从HS态突变为LS态的温度T_{\downarrow }进入本文的计算范围，这时整个系统的自旋态展现出温滞行为（如图5（c））。进一步可以看到随着高低自旋态的耦合强度J增大，温滞回线的宽度\left(T_{\uparrow }-T_{\downarrow }\right)会减小。随着分子内高低自旋态的耦合强度J的进一步增大，T_{\uparrow }与T_{\downarrow }靠近变成了在同一点发生转变，即高低自旋态转变的温滞行为消失，但此时转变仍是突变的（如图5（d））。随着分子内高低自旋态的耦合强度J的继续增加，当其超过某个临界值时，自旋态转变从突变行为变为渐变行为。

之所以存在自旋态的温滞行为，物理上可以简单理解为，在上述的分子内高低自旋态的耦合强度J较小时，系统的自由能在不同温度下可以出现两个极小值点，即稳定态和亚稳态；随着温度升高，HS亚稳态会转变为稳定态，由于分子间的协同相互作用，系统中大部分保留在初始的LS态；当温度继续升高，HS稳定态变为唯一极值点，即唯一稳定态，分子间的协同相互作用会使整个系统中的大部分分子才会转化为HS态。反之，当温度下降时，系统中大部分保留在初始HS态，由于分子间的协同相互作用，即使HS态变为亚稳态，高自旋态不会立即转变为LS态，只有当LS态变为唯一稳定态，此时由于分子间协同作用，系统中大部分分子才会变为LS态。当分子内高低自旋态的耦合强度J较小时，随着温度降低，在该文的温度范围内系统的自由能始终都有两个极值点，所以降温时由于分子间协同相互作用，系统一直保持在HS态。分子内高低自旋态的耦合强度J增大到一定程度时，随着温度的上升系统自由能从单个极值点（低温时）转化为两个极值点（适当温度时），然后继续升高温度，又从两极值点转化为一个极值点（更大温度时）。

接着考虑了分子自旋态随分子内势能高低自旋态两个极小值点的能量差的变化情况，这个能量差可由参数d来调节（如图6）。

由图6（a）知，随着分子内势能曲线中高低自旋态的能量差\mathrm{\Delta }E的增大，T_{\uparrow }从有限值缓慢增加，只有当能量差\mathrm{\Delta }E达到一定值时才会存在T_{\downarrow }，且该温度随能量差\mathrm{\Delta }E的增大而增大。分子自旋态随高低自旋态的能量差\mathrm{\Delta }E的变化图中可以分成三个区域（如图6（b）—图6（d））。当分子内势能高低自旋态的能量差\mathrm{\Delta }E较小时，系统随着温度的上升可以从初始的LS态转变为HS态，而随着温度的下降，系统在研究的温度范围内没有观测到自旋态的转变，正如图6（b）中所示，这时转变温度T_{\uparrow }存在，而从HS态突变为LS态的转变温度T_{\downarrow }不存在。当分子内势能曲线中高低自旋态的能量差\mathrm{\Delta }E增大时，在研究的温度范围内T_{\downarrow }出现，这时整个系统的自旋态随能量差\mathrm{\Delta }E的变化展现出温滞行为（如图6（c））。进一步可以看到，随着分子内势能的能量差\mathrm{\Delta }E的增大，温滞回线的宽度逐渐减小。随着分子内势能的能量差\mathrm{\Delta }E的进一步增大，T_{\uparrow }与T_{\downarrow }在计算精度范围内变成了在同一点发生转变，此时自旋态转变温度随能量差\mathrm{\Delta }E的温滞行为消失，但此时的转变仍保持是突变的（如图6（d））。随着分子内势能能量差\mathrm{\Delta }E的继续增加，当其超过某个临界值时，自旋态转变从突变行为变为渐变行为。

关于自旋态随能量差\mathrm{\Delta }E出现的温滞行为，与前面类似，当分子内势能的能量差\mathrm{\Delta }E较小时，自由能在不同温度下仍然可以出现两个极小值点，即稳定态和亚稳态；随着温度增加，起始稳定的LS态变为亚稳态，而HS亚稳态会转变为稳定态，但分子间的协同相互作用，使系统中大部分分子保留在初始的LS态。继续增大温度，HS稳定态变为唯一极值点，即唯一的稳定态，分子间的协同相互作用会使整个系统中的分子大部分转化为HS态。反之，当温度下降时，系统中大部分分子保留在初始HS态，由于分子间的协同相互作用，即使分子HS转变为亚稳态，HS不会立即转变成LS态，只有当LS态变为唯一稳定态时，由于分子间协同作用，系统中大部分分子才会变为LS态。当分子内势能的能量差\mathrm{\Delta }E较小时，随着温度降低，在模拟的温度范围内系统的自由能始终都有两个极值点，所以降温时由于分子间协同相互作用，系统一直保持在HS态。分子内势能的能量差\mathrm{\Delta }E增大到一定程度时，温度的上升使系统自由能从单极值点（低温时）转化为双极值点（在一定温度范围内），进一步升高温度，自旋系统又从双极值点转化为单极值点（更高温度时）。

最后，考虑了分子自旋态随高低自旋态的简并度比值的变化情况，可以通过改变参数a(b)，从而改变高低自旋态的简并度比值g，这里简并度比值定义为g=\sqrt[]{a/b}，改变简并度比值g一方面可以改变高低自旋态的熵差，另一方面会影响分子内势能曲线的弯曲程度，从而影响不稳定点（极大值点）所处的位置以及能量大小。因此，需要考虑分子自旋态随高低自旋态的简并度比值g的变化情况（如图7）。

由图7（a）知，随着分子内势能曲线中高低自旋态的简并度比值g的增大，T_{\uparrow }从有限值缓慢减小，T_{\downarrow }只有当简并度比值g达到一定值时才会出现，且该温度随自旋态的简并度比值g的增大而减小。选取不同的简并度比值g，做出系统自旋态随温度的变化图。当分子内高低自旋态的简并度比值g较小时，系统随着温度的上升表现为从初始的LS态逐渐转变为HS态，但不存在突变，此时不存在T_{\uparrow }和T_{\downarrow }（如图7（b））。当分子内高低自旋态的简并度比值增大超过某一临界值时，在这里模拟的温度范围存在T_{\uparrow }和T_{\downarrow }及突变温度。这时整个系统的自旋态随自旋态简并度比值g的变化有温滞行为（如图7（c））。进一步可以看到随着自旋态的简并度比值g的增大，温滞回线的宽度逐渐增大。随着分子内自旋态的简并度比值g的继续增加，当其超过某个临界值时，系统随着温度的上升可以从初始的LS态转变为HS态，而随着温度的下降，系统在模拟的温度范围内没有观测到自旋态的转变，正如图7（d）所示，这时转变温度T_{\uparrow }存在，而从HS突变为LS的转变温度T_{\downarrow }不存在。

关于自旋态随自旋态简并度比值g出现的温滞行为，当分子内自旋态简并度比值g较小时，自由能在不同温度下只有单极小值点；随着温度增加，起始只有单稳定态的LS态，而逐渐转变为HS的单稳态，由于分子间的协同相互作用，系统从初始LS态逐渐变为HS态。继续增大温度，HS稳定态变为唯一极值点，即唯一的稳定态，分子间的协同相互作用会使整个系统中的分子大部分转化为HS态。当分子内自旋态简并度比值g增大超过某临界值时，自由能在不同温度下仍然可以出现两个极小值点，即稳定态和亚稳态；与前面自由能双极值点的机理一样，同样可以解释此时的温滞行为。当自旋态简并度比值g增大到一定程度时，温度的上升使系统自由能从单极值点（低温时）转化为双极值点（在一定温度范围内），进一步升高温度，自旋系统又从双极值点转化为单极值点（更高温度时）。反之，当温度下降时，系统自由能从单极值点（高温时）转化为双极值点（在模拟的温度范围内），这时转变温度T_{\uparrow }存在，而模拟的温度范围内转变温度T_{\downarrow }不存在。

3 结论

本文考虑了一个简单的二维自旋转变模型，该模型包括了分子间相互作用和分子内相互作用。基于Nosé-Hoover的分子动力学方法研究了高低自旋态转变温度及高低自旋态平均占据数。首先检验了分子间相互作用参数对自旋转变温度的影响，得到了与已有结果一致的结论。接着重点探讨了分子内相互作用参数对自旋转变温度的影响。这里主要研究了分子内高低自旋态的耦合强度、分子内高低自旋态两个极小值点的能量差、自旋态的简并度比值对于自旋转变温度的影响。在温滞范围内自旋转变温度会随着分子内高低自旋态的耦合强度和分子内高低自旋态的能量差的增大而增大，而随着自旋态简并度比值的增大而减小；温滞宽度会随着高低自旋态的耦合强度和分子内高低自旋态的能量差的增大而减小，而随着自旋态简并度比值的增大而增大。通过减小分子内高低自旋态的耦合强度和分子内高低自旋态两个极小值点的能量差，或者增加自旋态的简并度比值可以更有利于观测自旋转变行为。如果要进一步观测自旋转变的温滞回线且有尽可能大的温滞范围，分子内高低自旋态的耦合强度和分子内高低自旋态两个极小值点的能量差存在一个下限值，而自旋态的简并度比值存在一个上限值。本文没有讨论发生自旋态突变时精确确定突变温度和相关参数的方法，这个问题可以采用更精确的基于平行回火算法的分子动力学模拟来解决。

参考文献

[1] OHKOSHI S I, IMOTO K, TSUNOBUCHI Y, et al. Light-induced Spin-crossover Magnet[J]. Nat Chem, 2011, 3(7): 564-569. DOI: 10.1038/nchem.1067 .

[2] NAKAYA M, OHTANI R, LINDOY L F, et al. Light-induced Excited Spin State Trapping in Iron(Iii) Complexes[J]. Inorg Chem Front, 2021, 8(2): 484-498. DOI: 10.1039/D0QI01188F .

[3] EKANAYAKA T K, MAITY K P, DOUDIN B, et al. Dynamics of Spin Crossover Molecular Complexes[J]. Nanomaterials, 2022, 12(10): 1742. DOI: 10.3390/nano12101742 .

[4] BAKER W A, BOBONICH H M. Magnetic Properties of some High-spin Complexes of Iron(II)[J]. Inorg Chem, 1964, 3(8): 1184-1188. DOI: 10.1021/ic50018a027 .

[5] FÉLIX G, NICOLAZZI W, SALMON L, et al. Enhanced Cooperative Interactions at the Nanoscale in Spin-crossover Materials with a First-order Phase Transition[J]. Phys Rev Lett, 2013, 110(23): 235701. DOI: 10.1103/PhysRevLett.110.235701 .

[6] BROOKER S. Spin Crossover with Thermal Hysteresis: Practicalities and Lessons Learnt[J]. Chem Soc Rev, 2015, 44(10): 2880-2892. DOI: 10.1039/c4cs00376d .

[7] SATO O. Dynamic Molecular Crystals with Switchable Physical Properties[J]. Nat Chem, 2016, 8(7): 644-656. DOI: 10.1038/nchem.2547 .

[8] ZHANG X, COSTA P S, HOOPER J, et al. Locking and Unlocking the Molecular Spin Crossover Transition[J]. Adv Mater, 2017, 29(39): 1702257. DOI: 10.1002/adma.201702257 .

[9] VALVERDE-MUÑOZ F J, SEREDYUK M, MUÑOZ M C, et al. Thermochromic Meltable Materials with Reverse Spin Transition Controlled by Chemical Design[J]. Angew Chem Int Ed Engl, 2020, 59(42): 18632-18638. DOI: 10.1002/anie.202006453 .

[10] 康婷, 焦虎军. 自旋转变材料中温度调制的自旋泵浦[J]. 山西大学学报(自然科学版), 2021, 44(1): 73-79. DOI: 10.13451/j.sxu.ns.2020044 .

[11] KANG T, JIAO H J. Temperature Modulated Spin Pumping in the Spin-transition Materials[J]. J Shanxi Univ Nat Sci Ed, 2021, 44(1): 73-79. DOI: 10.13451/j.sxu.ns.2020044 .

[12] GUTLICH P, KSENOFONTOV V, GASPAR A. Pressure Effect Studies on Spin Crossover Systems[J]. Coord Chem Rev, 2005, 249(17/18): 1811-1829. DOI: 10.1016/j.ccr.2005.01.022 .

[13] GASPAR A B, MOLNÁR G, ROTARU A, et al. Pressure Effect Investigations on Spin-crossover Coordination Compounds[J]. Comptes Rendus Chim, 2018, 21(12): 1095-1120. DOI: 10.1016/j.crci.2018.07.010 .

[14] LIU Z H, SAKAI Y, YANG J Y, et al. Sequential Spin State Transition and Intermetallic Charge Transfer in PbCoO₃ [J]. J Am Chem Soc, 2020, 142(12): 5731-5741. DOI: 10.1021/jacs.9b13508 .

[15] 宿凤琴, 焦虎军. 自旋转变材料中压强调制的自旋泵浦[J]. 山西大学学报(自然科学版), 2021, 44(6): 1123-1130. DOI: 10.13451/j.sxu.ns.2021024 .

[16] SU F Q, JIAO H J. Pressure Modulated Spin Pumping in the Spin-transition Materials[J]. J Shanxi Univ Nat Sci Ed, 2021, 44(6): 1123-1130. DOI: 10.13451/j.sxu.ns.2021024 .

[17] HALDER G J, KEPERT C J, MOUBARAKI B, et al. Guest-dependent Spin Crossover in a Nanoporous Molecular Framework Material[J]. Science, 2002, 298(5599): 1762-1765. DOI: 10.1126/science.1075948 .

[18] BAO X, SHEPHERD H J, SALMON L, et al. The Effect of an Active Guest on the Spin Crossover Phenomenon[J]. Angew Chem Int Ed, 2013, 52(4): 1198-1202. DOI: 10.1002/anie.201207218 .

[19] LI J Y, HE C T, CHEN Y C, et al. Tunable Cooperativity in a Spin-crossover Hoffman-like Metal–Organic Framework Material by Aromatic Guests[J]. J Mater Chem C, 2015, 3(30): 7830-7835. DOI: 10.1039/C5TC00432B .

[20] ZHANG C J, LIAN K T, HUANG G Z, et al. Hysteretic Four-step Spin-crossover in a 3D Hofmann-type Metal-organic Framework with Aromatic Guest[J]. Chem Commun, 2019, 55(74): 11033-11036. DOI: 10.1039/c9cc06017k .

[21] CAMBI L, SZEGO L. Uber Die Magnetische Susceptibilitat Der KomplexenVerbindungen[J]. Ber Dtsch Chem Ges,1931,64(10): 2591-2598.DOI:10.1002/cber.19310641002 .

[22] SPIERING H, BOUKHEDDADEN K, LINARES J, et al. Total Free Energy of a Spin-crossover Molecular System[J]. Phys Rev B, 2004, 70(18): 184106. DOI: 10.1103/physrevb.70.184106 .

[23] NISHINO M, BOUKHEDDADEN K, KONISHI Y, et al. Simple Two-dimensional Model for the Elastic Origin of Cooperativity Among Spin States of Spin-crossover Complexes[J]. Phys Rev Lett, 2007, 98(24): 247203. DOI: 10.1103/PhysRevLett.98.247203 .

[24] LORENC M, BALDE C, KASZUB W, et al. Cascading Photoinduced, Elastic, and Thermal Switching of Spin States Triggered by a Femtosecond Laser Pulse in an Fe(III) Molecular Crystal[J]. Phys Rev B, 2012, 85(5): 054302. DOI: 10.1103/physrevb.85.054302 .

[25] PAEZ-ESPEJO M, SY M, BOUKHEDDADEN K. Elastic Frustration Causing Two-step and Multistep Transitions in Spin-crossover Solids: Emergence of Complex Antiferroelastic Structures[J]. J Am Chem Soc, 2016, 138(9): 3202-3210. DOI: 10.1021/jacs.6b00049 .

[26] CRUDDAS J, POWELL B J. Spin-state Ice in Elastically Frustrated Spin-crossover Materials[J]. J Am Chem Soc, 2019, 141(50): 19790-19799. DOI: 10.1021/jacs.9b09191 .

[27] NDIAYE M, ISLAM BELMOURI N EL, LINARES J, et al. Elastic Origin of the Unsymmetrical Thermal Hysteresis in Spin Crossover Materials: Evidence of Symmetry Breaking[J]. Symmetry, 2021, 13(5): 828. DOI: 10.3390/sym13050828 .

[28] CRUDDAS J, POWELL B J. Multiple Coulomb Phases with Temperature-tunable Ice Rules in Pyrochlore Spin-crossover Materials[J]. Phys Rev B, 2021, 104(2): 024433. DOI: 10.1103/physrevb.104.024433 .

[29] AFFES K, SINGH Y, BOUKHEDDADEN K. Electro-elastic Modeling of Thermal Spin Transition in Diluted Spin-crossover Single Crystals[J]. Int J Mol Sci, 2022, 23(22): 13854. DOI: 10.3390/ijms232213854 .

[30] CHEN Z Y, LIU Q, CHENG Y, et al. Manipulating Electron-transfer Events in [Fe₄Co₄] Cubes via a Mixed-ligand Approach: The Impact of Elastic Frustration[J]. Angew Chem Int Ed Engl, 2023, 62(29): e202301124. DOI: 10.1002/anie.202301124 .

[31] NOSÉ S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods[J]. J Chem Phys, 1984, 81(1): 511-519. DOI: 10.1063/1.447334 .

[32] HOOVER W G. Canonical Dynamics: Equilibrium Phase-space Distributions[J]. Phys Rev A, 1985, 31(3): 1695-1697. DOI: 10.1103/physreva.31.1695 .

基金资助

国家自然科学基金(11574186)