金属锡4种经典原子间势函数的热力学性质比较

邹继德^1,2 ✉ 王丽芳³ ✉ 王越超³刘瑜³陈欣³赵继泽^1,2宋海峰³

1. 兰州大学物理科学与技术学院，甘肃兰州 730000

2. 兰州理论物理中心和甘肃省理论物理重点实验室，甘肃兰州 730000

3. 北京应用物理与计算数学研究所计算物理实验室，北京 100094

摘要： 本文选取了4种典型锡（Sn）原子间势函数：Pot_Ravelo1997，Pot_Vella2017，Pot_Eteami2018和Pot_Ko2018，对其在计算α相和β相Sn的晶格常数、平衡体积和能量、弹性常数和模量、等温压缩、等压热膨胀等热力学物性方面的适用性进行了系统比较评估，结合已有的实验和第一性原理数据比较结果，分别给出这4种势函数的准确性和适用范围：Pot_Ravelo1997能较好地再现α-，β-Sn的晶格常数、α-Sn的弹性常数和模量，并能够定性描述较宽温压范围内的Sn的等温压缩和等压热膨胀过程；Pot_Vella2017能够较好地给出α-Sn的弹性常数和模量，计算的β-Sn平衡体积与实验和密度泛函理论（DFT）计算结果具有一致性，但其计算晶格常数、β-Sn弹性常数和模量与实验值相比具有一定偏差；Pot_Eteami2018能较好再现α-，β-Sn的晶格常数、定性描述较宽温压范围的α-，β-Sn等温压缩和等压热膨胀过程；Pot_Ko2018能较好再现α-，β-Sn的晶格常数、β-Sn的弹性常数和模量，定量描述β-Sn等温压缩和等压热膨胀过程。本工作为金属锡分子动力学模拟选择合适的势函数提供了参考和依据。

关键词： 金属锡；分子动力学；势函数；热力学物性

基金： 国家自然科学基金项目(12004046,U2230401,U1930401,12247101)

收稿日期： 2023-03-28

金属锡（Sn）由于具有丰富的固-固和固-液相变现象，且相变条件对温度和压力都较为敏感^［1–5］，因此其结构、电子性质、热力学性质、相图、动态力学响应行为等一直是金属材料领域的研究热点^［6–9］。理论和实验的研究结果表明^［1］，Sn在0~194 GPa的压强下至少存在5种固相：金刚石结构的α相（ $F d \bar{3} m$ ）、体心四方的β相（I4₁/amd）、体心四方的γ相（I4/mmm）、体心立方的δ相（ $I m \bar{3} m$ ）和密排六方的ε相（P6₃/mmc）。近年来以分子动力学为主的材料计算模拟方法的发展为研究Sn的多相热力学性质、相图和动力学响应行为等提供了更好条件^{［10–12］}。1997年，Ravelo和Baskes^［13］开发了Sn的修正嵌入原子势（modified embedded atom method，MEAM），并研究了其α相和β相的结构性质、弹性性质和热力学性质，同实验结果的一致性较好。Bernard和Maillet^［14］基于该势函数模型发展了Sn的新MEAM，并采用分子动力学方法模拟研究了Sn在压力下熔化行为和冲击雨贡纽。同时，为研究Sn的动力学响应行为，Lane和Marder^［15］模拟研究了Sn中冲击波的传播演化过程。2014年，为了更好地再现Sn的冲击雨贡纽，Sapozhnikov等^［16］开发了Sn的嵌入原子势（embedded atom method，EAM），此势函数忽略了γ相和δ相的区别，适用于研究δ相段的熔化曲线问题，模拟计算得到的Sn的冲击温度与实验符合得较好。近几年，随着高性能计算模拟能力提升，千万原子量级的Sn冲击动力学响应行为研究开始涌现。Tang等^［17］采用分子动力学模拟了Sn在激光和平板冲击加载下的表面微喷现象，Wang等^［18］以Sn为例研究了单晶金属在平面波应力加载下的层裂损伤机制。此外，另有研究学者则针对β-Sn展开了系列研究，包括其表面、晶界等界面性质，Sn的原子表面扩散和自扩散行为、β-Sn的黏性、以及β-Sn的应力应变曲线等^{［19–22］}。

尽管分子动力学方法在研究Sn的物性方面具有不可替代的作用，但其计算模拟的准确性依赖于原子间势函数。除了上述提到的1997年开发的Sn的MEAM和2014年开发的Sn的EAM，2017年，Vella等^［23］为了更好地描述液相Sn的性质，通过拟合固相的内聚能、体积模量、平衡体积和液相的密度、自扩散系数等数据对MEAM势函数参数进行了优化。2018年，Ko等^［24］针对β-Sn的热力学性质进一步优化了MEAM势函数。同时，Etesami等^［25］将该势函数模型应用于Pb-Sn合金体系中，并对相关参数做了修正和优化。目前针对Sn开发的原子间势函数较为丰富，但尚缺乏对这些势函数准确性、适用性的系统评估，这影响了这些势函数在Sn物性模拟研究中的进一步应用。因此，本文选取了4种MEAM势函数模型：Pot_Ravelo1997^［13］，Pot_Vella2017^［23］，Pot_Eteami2018^［25］和Pot_Ko2018^［24］，对其在计算模拟α相和β相Sn的晶格常数、平衡体积和能量、弹性常数和模量、等温压缩、等压热膨胀等热力学物性方面的适用性进行了系统计算评估，结合已有的实验和第一性原理数据比较结果，分别给出这4种势函数的准确性和适用范围，为Sn分子动力学模拟选择合适的势函数提供参考和依据。

1 方法

1.1 原子间势函数

评估的4种Sn势函数都采用了Baskes等^{［26–27］}提出的MEAM模型，在该势函数模型中，体系原子总能量表达为： $E_{tot} = \frac{1}{2} \sum_{i, j (j \neq i)} Φ_{s_{j} s_{i}} (r_{i j}) + \sum_{i} F_{s_{i}} (\bar{ρ_{i}})$ （1）式中，等号右边第一项为体系对势，第二项为体系原子嵌入能，代表体系的多体效应，是原子坐标和电子密度的函数。特别地，为更好的描述过渡金属体系电子的d轨道杂化，使原子键呈现一定程度上的共价性，电子密度函数被视为张量，并将 $\bar{ρ_{i}}$ 表达为张量不变量的函数，即： ${(\bar{ρ_{i}})}^{2} = {({\bar{ρ_{i}}}^{(0)})}^{2} + {({\bar{ρ_{i}}}^{(1)})}^{2} + {({\bar{ρ_{i}}}^{(2)})}^{2} + {({\bar{ρ_{i}}}^{(3)})}^{2}$ （2） ${({\bar{ρ_{i}}}^{(0)})}^{2} = {[\sum_{j \neq i} ρ_{s_{j}}^{(0)} (r_{i j})]}^{2}$ （3） ${({\bar{ρ_{i}}}^{(1)})}^{2} = \sum_{α} {[\sum_{j \neq i} ρ_{s_{j}}^{(1)} (r_{i j}) \frac{r_{i j}^{α}}{r_{i j}}]}^{2}$ （4） $\begin{array}{l} {({\bar{ρ_{i}}}^{(2)})}^{2} = \sum_{α, β} {[\sum_{j \neq i} ρ_{s_{j}}^{(2)} (r_{i j}) \frac{r_{i j}^{α} r_{i j}^{β}}{r_{i j}^{2}}]}^{2} - \\ \frac{1}{3} {[\sum_{j \neq i} ρ_{s_{j}}^{(2)} (r_{i j})]}^{2} \end{array}$ （5） ${({\bar{ρ_{i}}}^{(3)})}^{2} = \sum_{α, β, γ} {[\sum_{j \neq i} ρ_{s_{j}}^{(3)} (r_{i j}) \frac{r_{i j}^{α} r_{i j}^{β} r_{i j}^{γ}}{r_{i j}^{3}}]}^{2}$ （6）在这样的修正下，势函数的精度提高了，可以很好的模拟过渡金属体系，甚至一些共价固体体系^［28］。

1.2 分子动力学模拟

分子动力学模拟均通过LAMMPS软件实现^［29］。模拟体系均采用9×9×9的超胞和周期性边界条件。对于晶格常数、平衡体积以及弹性常数和模量计算，分子静力学能量和力收敛条件设置分别为1×10⁻¹² eV和1×10⁻¹² eV·Angstram⁻¹，经测试满足收敛要求。对等温压缩和等压热膨胀等热力学过程模拟，采用等压等温（NPT）系综进行温压控制，模拟总时长600 ps，并在200~600 ps时间内每隔1 ps进行采样和统计平均，以获得该温压条件下体系温度、压力、体积和内能等热力学量。动力学模拟时间步长为0.001 ps。

2 结果与讨论

2.1 晶格常数、平衡体积和能量

首先比较4种势函数在计算Sn两种晶格结构：α相和β相方面的表现和性质。Sn的α相和β相晶格结构如图1所示，α-Sn为金刚石结构，β-Sn为体心四方结构。4种势函数计算获得的α-Sn和β-Sn晶格常数、平衡体积以及内聚能如表1所示。可以看到，Pot_Ravelo1997势函数给出与实验一致的α-Sn晶格常数、平衡体积和内聚能，β-Sn晶格常数和内聚能比实验值略高；Pot_Vella2017势函数计算的α-Sn和β-Sn晶格常数都明显低于实验值和第一性原理计算值；Pot_Etesami2018势函数计算α-Sn晶格常数与基于Generalized Gradient Approximation （GGA）的密度泛函理论（DFT）计算结果较接近，β-Sn晶格常数与基于Local-Density Approximation （LDA）的DFT计算结果和实验值较一致，α-Sn和β-Sn的内聚能均接近实验值；Pot_Ko2018势函数给出均接近实验值的α-Sn和β-Sn晶格常数和内聚能。总体来看，4种势函数中，除Pot_Vella2017明显低估α-Sn的晶格常数和平衡体积、β-Sn的晶格常数外，其他4种势函数均能较好给出α-Sn和β-Sn的晶格常数、平衡体积和能量性质。

图1 Sn的α相和β相晶格结构

Fig.1 Lattice structure of tin with α phase （a）and β phase （b）

表1 α-Sn和β-Sn的晶格常数、平衡体积和内聚能

Table 1 Lattice constants (LC)，equilibrium volumes and cohesion energies of α-Sn and β-Sn

Structures	Properties		Interatomic potentials				DFT		Exp.
Structures	Properties		Pot_Ravelo1997	Pot_Vella2017	Pot_Etesami2018	Pot_Ko2018	LDA^［30］	GGA^［24］	Exp.
α	LC	a/nm	0.6483	0.6304	0.6581	0.6430	0.6462	0.6658	0.6483^［31］
	Eq.vol/（nm³·atom⁻¹）		0.034050	0.031308	0.035634	0.033224	0.03372		0.03405^［31］
	E_c/（eV·atom⁻¹）		–3.140	–3.220	–3.135	–3.209	–3.723	–3.154	–3.14^［32］
β	LC	a/nm	0.5920	0.5682	0.5859	0.5914	0.5800	0.5938	0.5831^［31］
		c/nm	0.3235	0.3334	0.3206	0.3234	0.3141	0.3224	0.3184^［31］
		c/a	0.546	0.587	0.547	0.547	0.541	0.543	0.546^［31］
	Eq.vol/（nm³·atom⁻¹）		0.028345	0.026911	0.027512	0.028301	0.02642		0.02707^［31］
	E_c/（eV·atom⁻¹）		–3.085	–3.115	–3.102	–3.091	–3.688	–3.113	–3.10^［33］
	ΔE/（eV·atom⁻¹）		0.055	0.105	0.033	0.118	0.035	0.041	0.03

2.2 弹性常数和模量

计算了4种势函数在α-Sn和β-Sn弹性常数和模量方面的性质和表现，如表2所示，与已有文献计算结果一致^［24］。可以看到，对于Pot_Ravelo1997势函数，其计算得到的α-Sn弹性常数和模量均与实验符合较好，但β-Sn的两个弹性常数C₁₁和C₃₃则明显高于实验值；Pot_Vella2017势函数函数计算得到的α-Sn的弹性常数和弹性模量与实验值和DFT数据较为接近，但β-Sn的弹性常数与实验值和第一性原理计算值偏差较大；对于Pot_Etesami2018势，除α-Sn的弹性模量与DFT结果较为接近外，其余的弹性常数与β-Sn的弹性模量均与实验结果存在较大偏差；对于Pot_Ko2018势函数，其α-Sn弹性常数与实验值和DFT计算结果偏离较大，弹性模量与实验值接近，而其计算得到的β-Sn弹性常数和模量均与实验值较为一致。由此，可以认为，在弹性常数和模量计算方面，对于α-Sn，Pot_Ravelo1997和Pot_Vella2017具有较准确的计算结果和性质表现，而对于β-Sn，Pot_Ko2018具有较准确的计算结果和性质表现；Pot_Etesami2018计算性质不够理想。由于4种势函数均不能同时描述和计算α-Sn和β-Sn，因此涉及α，β两相共存或α-β相变的金属Sn弹性常数和模量计算的体系，目前尚没有合适的势函数进行计算模拟，有待新的Sn势函数开发和参数优化。

表2 α-Sn和β-Sn的弹性常数和模量

Table 2 Elastic constants and modulus of α-Sn and β-Sn

Structures	Properties	Interatomic potentials				DFT^［30］	Exp.^［34］
Structures	Properties	Pot_Ravelo1997	Pot_Vella2017	Pot_Ko2018	Pot_Etesami2018	LDA	Exp.^［34］
α	C₁₁/GPa	70.358	64.884	50.413	81.949	67.6	69.1
	C₁₂/GPa	28.121	33.873	35.660	24.429	31.0	21.3
	C₄₄/GPa	36.757	42.578	10.532	94.856	44.0	42.6
	Bulk modulus/GPa	42.200	44.210	40.578	43.602	43.2	42.6
β	C₁₁/GPa	109.255	123.312	89.745	132.879	73.7	73.4
	C₁₂/GPa	62.522	34.258	46.675	46.905	61.9	59.9
	C₁₃/GPa	24.426	54.027	36.915	16.558	32.8	39.1
	C₃₃/GPa	139.634	59.179	93.729	157.116	95.3	90.7
	C₄₄/GPa	0.707	4.163	7.933	1.158	23.2	22.0
	C₆₆/GPa	22.495	21.500	10.626	28.097	24.7	23.9
	Bulk modulus/GPa	64.544	65.603	57.136	64.768	55.3	57.0

2.3 等温压缩性质

采用4种势函数分别模拟计算了α-Sn和β-Sn的等温压缩线，结果分别如图2和图3所示。其中，作为比较参考，在图3的β-Sn的等温压缩线图中给出Salamat等^［35］的实验结果和本课题组研究的零温DFT结果。对于α-Sn，4种势函数给出定性一致的等温压缩线变化趋势：不同温度条件下，α-Sn的等温压缩线呈现斜率随压强增加而逐渐减小的凹函数特征，意味着α-Sn的等温压缩系数随压力的增加而减小。另一方面，由图2可以看到，不同势函数给出不同的α-Sn构型稳定温压区间和平衡体积：对于Pot_Ravelo1997势函数，α-Sn在200~1400 K，–2~10 GPa的温压区间呈现分子动力学模拟尺度上的稳定性；对于Pot_Vella2017，α-Sn在200~1200 K，–2~10 GPa的温压区间能够稳定保持晶格结构；对于Pot_Etesami2018，α-Sn在200~2200 K，–2~10 GPa的温压区间能够保持晶格结构稳定；而对于Pot_Ko2018，α-Sn仅在200~800 K，–2~9 GPa的极窄温压范围保持晶格结构。这意味着Pot_Ko2018势函数在较宽温压区间内对α-Sn的适用性有待改进。

图2 不同势函数模拟计算的α-Sn的等温压缩线

Fig.2 Isothermal compression line for α-Sn with different interatomic potentials

（a）Pot_Ravelo1997；（b）Pot_Vella2017；（c）Pot_Eteami2018；（d）Pot_Ko2018

图3 不同势函数模拟计算的β-Sn的等温压缩线（300 K-EXP.数据来源于参考文献［35］）

Fig.3 Isothermal compression lines for β-Snwith different interatomic potentials （300 K-EXP. data are derived from Ref.［35］）

（a）Pot_Ravelo1997；（b）Pot_Vella2017；（c）Pot_Eteami2018；（d）Pot_Ko2018

由图3可以看出，对于β-Sn，4种势函数给出与实验和DFT计算结果定性一致的等温压缩线变化趋势：不同温度条件下，β-Sn的等温压缩线也呈现斜率逐渐减小的凹函数特征，意味着β-Sn的等温压缩系数也随压力的增加而减小。同时，由图3可以看到，不同势函数计算模拟下，β-Sn晶格结构的稳定温度上限，即熔化温度随压力变化：压力增加，熔化温度升高，与实验认识一致^［36］。定量上，以实验和DFT计算结果为参考，Pot_Vella2017势和Pot_Ko2018势计算获得的β-Sn平衡体积与实验和DFT结果较为接近，但Pot_Ravelo1997势和Pot_Eteami2018势都高估了β-Sn的平衡体积。另一方面，4种势函数给出定量不同的β-Sn晶格结构保持稳定的稳压区间：对于Pot_Ravelo1997势函数，β-Sn在0~1600 K，0~20 GPa的温压区间呈现分子动力学模拟尺度上的稳定性；对于Pot_Vella2017，β-Sn在0~1400 K，0~20 GPa的温压区间能够稳定保持晶格结构；对于Pot_Etesami2018，β-Sn在0~2200 K，0~20 GPa的温压区间能够保持晶格结构稳定；而对于Pot_Ko2018，β-Sn在0~1400 K，0~20 GPa的温压范围保持晶格结构。可见，4种势函数对β-Sn的等温压缩性质描述较好，都适用于较宽温压范围。综上，针对β-Sn的等温压缩性质，4种势函数都能获得与实验及DFT定性一致的计算结果，且均具有面向以β-Sn为初始状态的动力学响应行为模拟研究的应用适用性；但定量上Pot_Vella2017和Pot_Ko2018更接近实验和DFT值。

2.4 等压热膨胀性质

等压条件下体系的热膨胀过程同样可以反映势函数在热力学物性计算模拟方面的适用性。采用4种势函数分别模拟计算了α-Sn和β-Sn的等压热膨胀线，结果如图4和图5所示。对于α-Sn，4种势函数给出定性一致的等压热膨胀规律：较高压力条件下，α-Sn体积随温度的变化呈现较好的线性变化趋势；较低压力至负压条件下，α-Sn体积随温度升高而增加呈现非线性趋势，热膨胀效应随温度升高加剧。在定量规律上，由图4可以明显看到，不同势函数的适用温压范围明显不同：对于Pot_Ravelo1997势函数，α-Sn在0~1400 K，–2~10 GPa的温压区间呈现分子动力学模拟尺度上的稳定性；对于Pot_Vella2017，α-Sn在0~1200 K，–2~10 GPa的温压区间能够稳定保持晶格结构；对于Pot_Etesami2018，α-Sn在0~2200 K，–2~10 GPa的温压区间能够保持晶格结构稳定；而对于Pot_Ko2018，α-Sn仅在0~800 K，–2~9 GPa的极窄温压范围保持晶格结构。这与由等温压缩线计算获得的势函数在分子动力学尺度上的构型稳定性一致。

图4 不同势函数模拟计算的α-Sn的等压热膨胀线

Fig.4 Isobaric thermal expansion line of α-Sn with different interatomic potentials

（a）Pot_Ravelo1997；（b）Pot_Vella2017；（c）Pot_Eteami2018；（d）Pot_Ko2018

图5 不同势函数模拟计算的β-Sn的等压热膨胀线

Fig.5 Isobaric thermal expansion line of β-Sn with different interatomic potentials

（a）Pot_Ravelo1997；（b）Pot_Vella2017；（c）Pot_Eteami2018；（d）Pot_Ko2018

对于β-Sn，由图5可以看到，4种势函数同样给出定性一致的等压热膨胀曲线变化趋势：在0~20 GPa的稳定压力区间范围内，β-Sn的平衡体积随温度都呈现逐渐增加的非线性趋势，意味着β-Sn的热膨胀系数具有随温度增加的变化规律。而在定量描述方面，4种势函数同样给出定量不同的β-Sn晶格结构保持稳定的稳压区间：对于Pot_Ravelo1997势函数，β-Sn在0~1600 K，0~20 GPa的温压区间呈现分子动力学模拟尺度上的稳定性；对于Pot_Vella2017，β-Sn在0~1400 K，0~20 GPa的温压区间能够稳定保持晶格结构；对于Pot_Etesami2018，β-Sn在0~2200 K，0~20 GPa的温压区间能够保持晶格结构稳定；而对于Pot_Ko2018，β-Sn仅在0~1400 K，0~20 GPa的温压范围保持晶格结构。这也与由等温压缩线计算获得的势函数在分子动力学尺度上的构型稳定性一致。但值得注意的是，Pot_Ravelo1997和Pot_Eteami2018都给出了整体偏大的β-Sn平衡体积，意味着在采用这两种势函数进行以α-或β-Sn为初始状态的动力学响应行为模拟研究时，存在低估初始体系密度的潜在问题，建议使用相对密度进行标定。

3 结论

针对金属Sn，系统评估和比较了4种典型原子间势函数：Pot_Ravelo1997，Pot_Vella2017，Pot_Eteami2018和Pot_Ko2018，分别在计算模拟α相和β相Sn的晶格常数、平衡体积和能量、弹性常数和模量、等温压缩、等压热膨胀等热力学物性方面的适用性。结果发现，Pot_Ravelo1997能较好再现α-，β-Sn的晶格常数、α-Sn的弹性常数和模量，并能够定性描述较宽温压范围内的Sn的等温压缩和等压热膨胀过程；Pot_Vella2017能够较好地重复α-Sn的弹性常数和模量，计算的β-Sn平衡体积与实验和DFT结果一致性良好，但其处理晶格常数、β-Sn弹性常数和模量与实验值存在一定偏差；Pot_Eteami2018能较好再现α-，β-Sn的晶格常数、定性描述较宽温压范围的α-，β-Sn等温压缩和等压热膨胀过程；Pot_Ko2018能较好再现α-，β-Sn的晶格常数、β-Sn的弹性常数和模量，定量描述β-Sn等温压缩和等压热膨胀过程。综上，本文认为Pot_Ravelo1997和Pot_Eteami2018势在对α-，β-Sn的定性模拟方面具有更好的适用范围和稳定性，Pot_Ko2018势针对β-Sn具有更准确的定量计算结果。

参考文献：

［1］

Rehn D A，Greeff C W，Burakovsky L，Sheppard D G，Crockett S D. Multiphase tin equation of state using density functional theory［J］.Physical Review B，2021，103（18）：184102.

［2］

Weir S，Lipp M，Falabella S，Samudrala G，Vohra Y. High pressure melting curve of tin measured using an internal resistive heating technique to 45 GPa［J］.Journal of Applied Physics，2012，111（12）：123529.

［3］

Hafner J. Ab initio calculation of the pressure-induced A4→ A5→ A2 （distorted）→ A3 phase transitions in tin［J］.Physical Review B，1974，10（10）：4151.

［4］

Katzke H，Bismayer U，Tolédano P. Theory of the high-pressure structural phase transitions in Si，Ge，Sn，and Pb［J］.Physical Review B，2006，73（13）：134105.

［5］

Christensen N E，Methfessel M. Density-functional calculations of the structural properties of tin under pressure［J］.Physical Review B，1993，48（9）：5797.

［6］

Prasad S C，Wooster W A. The study of the elastic constants of white tin by diffuse X-ray reflexion［J］.Acta Crystallographica，1955，8（11）：682.

［7］

Groves S，Paul W. Band structure of gray tin［J］.Physical Review Letters，1963，11（5）：194.

［8］

Thewlis J，Davey A R. Thermal expansion of grey tin［J］.Nature，1954，174（4439）：1011.

［9］

Mabire C，Hereil P L. Shock Induced Polymorphic Transition and Melting of Tin［A］.AIP Conference Proceedings ［C］.American Institute of Physics，2000. 1.

［10］

郑佳琳，姜森宝，汤泽军，叶福田，李细锋.锆扩散连接行为的分子动力学研究［J］.塑性工程学报，2023，30（11）：116.

［11］

杜春志，庞帅，吴文平，刘兵飞.单晶Ni₃Al裂纹扩展行为的分子动力学模拟［J］.锻压技术，2023，48（7）：255.

［12］

肖寒，陈昊，崔鋆昕，周瑀杭，陈磊，张雄超.退火工艺对铸态ZCuSn10P1锡青铜组织性能的影响［J］.稀有金属，2023，47（6）：815.

［13］

Ravelo R，Baskes M. Equilibrium and thermodynamic properties of grey，white，and liquid tin［J］.Physical Review Letters，1997，79（13）：2482.

［14］

Bernard S，Maillet J B. First-principles calculation of the melting curve and Hugoniot of tin［J］.Physical Review B，2002，66（1）：012103.

［15］

Lane J M D，Marder M P. Molecular dynamics of shock fronts and their transitions in Lennard-Jonesium and tin［J］.arXiv Preprint Cond-Mat，2006. https：//doi.org/10.48550/arXiv.cond-mat/0607335.

［16］

Sapozhnikov F，Ionov G，Dremov V，Soulard L，Durand O. The Embedded atom Model and Large-scale MD Simulation of Tin under Shock Loading［A］.Journal of Physics：Conference Series ［C］.IOP Publishing，2014. 1.

［17］

Tang G Q，Li B，Deng X L，Wang L，Zhu W J，Hu W Y，Gao F. Comparative investigation of microjetting from tin surface subjected to laser and plane impact loadings via molecular dynamics simulations［J］.Mechanics of Materials，2020，148：103479.

［18］

Wang X X，He A M，Zhou T T，Wang P. Spall damage in single crystal tin under shock wave loading：a molecular dynamics simulation［J］.Mechanics of Materials，2021，160：103991.

［19］

Sellers M S，Schultz A J，Basaran C，Kofke D A. β-Sn grain-boundary structure and self-diffusivity via molecular dynamics simulations［J］.Physical Review B，2010，81（13）：134111.

［20］

Voter A F，Montalenti F，Germann T C. Extending the time scale in atomistic simulation of materials［J］.Annual Review of Materials Research，2002，32（1）：321.

［21］

Lee Y，Basaran C. Atomic-level shear stress-strain behavior of β-Sn［J］.Journal of Nanomechanics and Micromechanics，2013，3（4）：04013003.

［22］

Lee Y，Basaran C. Molecular dynamics of viscoplasticity in β-tin lattice and grain boundary［J］.Computational Materials Science，2013，68：290.

［23］

Vella J R，Chen M，Stillinger F H，Carter E A，Debenedetti P G，Panagiotopoulos A Z. Structural and dynamic properties of liquid tin from a new modified embedded-atom method force field［J］.Physical Review B，2017，95（6）：064202.

［24］

Ko W S，Kim D H，Kwon Y J，Lee M H. Atomistic simulations of pure tin based on a new modified embedded-atom method interatomic potential［J］.Metals，2018，8（11）：900.

［25］

Etesami S A，Baskes M I，Laradji M，Asadi E. Thermodynamics of solid Sn and PbSn liquid mixtures using molecular dynamics simulations［J］.Acta Materialia，2018，161：320.

［26］

Baskes M I. Modified embedded-atom potentials for cubic materials and impurities［J］.Physical Review B，1992，46（5）：2727.

［27］

Baskes M I，Angelo J E，Bisson C L. Atomistic calculations of composite interfaces［J］.Modelling and Simulation in Materials Science and Engineering，1994，2（3）：505.

［28］

文玉华，朱如曾，周富信，王崇愚.分子动力学模拟的主要技术［J］.力学进展，2003，1：65.

［29］

Thompson A P，Aktulga H M，Berger R，Bolintineanu D S，Brown W M，Crozier P S，Veld P J，Kohlmeyer A，Moore S G，Nguyen T D，Shan R，Stevens M J，Tranchida J，Trott C，Plimpton S J. LAMMPS-a flexible simulation tool for particle-based materials modeling at the atomic，meso，and continuum scales［J］.Computer Physics Communications，2022，271：108171.

［30］

Aguado A. First-principles study of elastic properties and pressure-induced phase transitions of Sn：LDA versus GGA results［J］.Physical Review B，2003，67（21）：212104.

［31］

Barrett C S，Massalski T B. Structure of Metals［M］.New York：McGrawHill，1966. 1.

［32］

Kittel C. Introduction to Solid State Physics［M］.New York：John Wiley，1986.185.

［33］

Ihm J，Cohen M L. Equilibrium properties and the phase transition of grey and white tin［J］.Physical Review B，1981，23（4）：1576.

［34］

Brandes E A，Brook G B. Smithells Metals Reference Book［M］.Butterworth-Heinemann：Oxford，1992. 1044.

［35］

Salamat A，Briggs R，Bouvier P，Petitgirard S，Dewaele A，Cutler M E，Corà F，Daisenberger D，Garbarino G，McMillan P F. High-pressure structural transformations of Sn up to 138 GPa：angle-dispersive synchro-tron X-ray diffraction study［J］.Physical Review B，2013，88（10）：104104.

［36］

Dudley J D，Hall H T. Experimental fusion curves of indium and tin to 105 000 atmospheres［J］.Physical Review，1960，118（5）：1211.