过程工程学报
The Chine Journal of Process Engineering
第21卷第3期2021年3月
Vol.21 No.3
Mar. 2021
.-------
1流动与传递匸
DOI: 10.12034/j.issn.l 009-606X.220107
¢35-
Application of cutoff distance lection in molecular dynamics
simulation of LJ argon system
Chenyang SUN 1'2, Chaofeng HOU 1*, Wei GE 1-2
1. Institute of Process Engineering, Chine Academy of Sciences, Beijing 100190, China
2. School of Chemical Engineering, University of Chine Academy of Sciences, Beijing 100049, China
Abstract: In Lennard-Jones (LJ)
potential argon system investigated by
molecular dynamics simulation, the
different cutoff distances are frequently
applied to calculate the interactive force between atoms, and some reports have
given out the influence of the cutoff
distances on the simulation systems.
More and more calculations suggest to u 4.5 cr or even larger truncation
distances (a is the diameter of argon atom) to obtain the more accurate thermodynamic properties of the
1.0 1.1
0.9Temperature
Mapping between thermodynamic states under different r c .
s
c
tu o
d o E o l
e -I o c I
-
1.5
-2.0
-2.5-3.0
-3.5
-4.0-4.5
Lennard-Jones potential
c w
.u
-g
p u n j
u o u n q
'E
-
.s
p e y
0；8 '=
"a u o 3o l n e
Comparison of system properties under different r c .
systems. In this work, a simple method was propod to solve the problem of superheating encountered in the calculation of melting point by direct heating, where an independent track and enmble at each temperature point are run. And then, the effect of different cutoff distances on the pha diagrams of the melting and boiling points of argon
system in the NPT enmble was studied. The melting point was in good agreement with th
e experimental and theoretically calculated values when the cutoff d istance of 2.5crwas ud. However, the deviation from the experimental
家常红烧肉做法melting point became more evident when the larger cutoff distances were employed. In order to find out the underlying mechanism behind the deviation, the radial distribution functions and velocity autocorrelation function at the melting
points and different thermodynamic states of the liquid argon with different cutoff distances were analyzed. It was
found that the same thermodynamic properties can be obtained at the corresponding thermodynamic state points under
different truncation distances. The mapping between the different thermodynamic state points was understandable due
to the different thermodynamic states at the same temperatures under the varied truncation distances, and was beneficial to significantly reduce the computational workload at the smaller cutoff distance. This work propod an exploratory
way for the lection of the cutoff distance in the simulation of liquid argon, where the truncation distance of 2.5 cr can
抗日英雄有哪些
meet the requirements of computational accuracy and performance in the simulations.
Key words: molecular dynamics simulation; LJ potential; cutoff distance; liquid argon
收稿：2020-03-27.修回：2020-05-10，网络发表：2020-05-25, Received: 2020-03-27, Revid: 2020-05-10, Published online: 2020-05-25
基金项目：国家自然科学基金(编号：21776280; 91834303; 91934302)；中国科学院资助项目(编号：XDC01040100; QYZDJ-SSW-JSC029)；北京市自然
科学基金委员会■■北京市教育委员会联合资助项目(编号：KZ201910017019)；多相复杂系统国家重点实验室自主研究课题(编号：MPCS-2019-A-10)
作者简介：孙晨阳(1993-)，男，山东省滨州市人，硕士研究生，化学工程专业，E-mail: ************* ：侯超峰，通讯联系人，E-mail:
************ .
引用格式：孙晨阳，侯超峰，葛蔚.LJ 势氮系统分子动力学模拟中截断半径的选择.过程工程学报,2021,213:259-264.
Sun C Y, Hou C F, Ge W. Application of c utoff distance lection in molecular dynamics simulation of LJ argon system (in Chine). Chin. J. Process Eng., 2021, 21(3): 259-264, DOI: 10.12034/j.issn. 1009-606X.220
107.
260过程工程学报第21卷LJ势氮系统分子动力学模拟中截断半径的选择
孙晨阳12，侯超峰"，葛蔚12
1.中国科学院过程工程研究所，北京100190
2.中国科学院大学化学工程学院，北京100049
摘要：近年来在分子动力学方法研究LJ(Lennard-Jones)势氨系统时，越来越多的计算建议釆用4.5cr甚至更大的截断半径为氨原子直径)。本工作研究了不同截断半径对等温等压NPT系综氮系统熔沸点相图的影响，分析了不同截断半径下在氨系统的熔点及液相区域不同热力学状态点的径向分布函数和速
度自相关函数。结果表明，以熔点为基准，在距离其相同液态温区分率的热力学温度点能获得相同的热力学性质，2.5a的截断半径在模拟的准确性和计算性能上均能满足模拟的要求，本工作为液氮模拟中截断半径的选择指明了一个新的思考方向。
关键词：分子动力学模拟；LJ势；截断半径；液氫
中图分类号：TQ460.6+4文献标识码：A文章编号：1009-606X(2021)03-0259-06
1前言
分子动力学模拟(Molecular Dynamics Simulation, MD)是一种在原子分子尺度的计算机模拟方法，主要依靠对经典牛顿运动方程的求解模拟体系内粒子的运动，得到粒子的位置信息和动量信息，再根据经典力学和统计物理理论获得体系的热力学量和其他宏观性质，是研究凝聚态平衡和非平衡热力学过程的重要理论手段山。随着计算机计算能力的迅猛发展，分子动力学模拟技术逐步成熟。目前，分子动力学模拟已广泛应用于化学、物理、材料、生物、医药等领域，在解释物质和热力学过程的微观机理方面发挥了不可替代的作用。
在氨原子体系的MD模拟中，截断半径的合理使用对模拟的准确性和积分算法的选取有重要作用R現为减少计算量，节约计算资源，通常计算时只考虑与原子相距小于截断半径r c(Cutoff Distance)的近邻
原子的贡献。但统计研究发现模拟过程中原子间力的截断取舍导致的能量损耗，可能会对系统的热力学性质有明显影响。近年来由于计算能力提升，对于采用LJ势函数的模拟系统，越来越多的计算建议采用4.5c甚至更大的截断半径。Rutkai等⑷考虑，由于完全忽略了超过r c=2.5aM远程相互作用，使用带有截断的LJ势在计算惰性元素液体密度、等压热容、声速等方面与通过状态函数理论计算的结果存在明显偏差。此外，在气液两相界面动力学性质和界面密度同、纳米流体表面力场的计算同和高频率相互作用诱导散射谱⑴等研究方面，均建议使用相对较大的截断半径来避免大的模拟偏差。然而对于径向分布函数和速度自相关函数的计算，虽然己有文献研究截断半径的影响，但并未分析影响背后的底层物理机制。
本工作通过研究不同截断半径下NPT系综氮系统的熔沸点相图及径向分布函数和速度自相关函数的变化，提出了一种氮体系分子动力学模拟中截断半径选择的新的思考方向。
2模拟方法
在分子动力学模拟中，原子间作用势函数的选取对模拟结果起决定性作用。Lennard-Jones(LJ)12-6势在氨原子体系的模拟中广泛应用，是准确描述氮等惰性气体原子之间相互作用的最重要的势函数之一。常用的Lennard-Jones12-6势函数的表达式如下：
式中，参数&与o分别为能量特征常数与长度特征常数，巾为粒子i与丿之间的距离。
对所有物理量以体系中原子质量皿、Lennard-Jones 势中的能量特征参数£、长度特征参数理行无量纲化处理，分子间距r=r/a，长度L=Lla，温度T=kT/s，能量时间t*=gm巧叭氮的基本参数歸］为：单位长度0.3405nm，单位质量6.636x10_26kg，单位温度119.8 K，单位时间2.157x10-%，单位能量1.655x10-21J。
应用开源软件LAMMPS(22Augl8)进行算例设计与计算。模拟体系包含2048个原子，采用周期性边界条件，模拟系综为NPT，使用Verlet速度时间积分和
第3期孙晨阳等:LJ势氫系统分子动力学模拟中截断半径的选择261
No-Hoover控温控压方法，统计的时间步长为0.001（约化单位制）。分别在截断半径为2.0cr,2.5o;3.50；4.5<T, 5.0cr下分析NPT系综氮体系各自熔点和液体区域2个对应温度（相对熔点增加的相对温度分率为50%，见表1）的径向分布函数和速度自相关函数。粒子的初始位置按面心立方FCC晶格构型排布，初始速度分布满足麦克斯韦-玻尔兹曼分布。
表1p=101kPa时NPT系中综氯原子系统不同截断半径下的液氟热力学状态温度范围
Table1Thermodynamic temperature range of liquid argon at different cutoff distances at0=101kPa in the NPT enmble Cutoff distance Melting point,T m&+50%厶Boiling point，兀Liquid temperature range,L
2.0cr0.630.720.810.18
2.5(y0.740.850.960.22
3.5cr0.790.935 1.080.29
4.5cr0.810.965 1.120.31
5.0cr0.810.97L130.32
Note:Experimental melting point and boiling point of argon system are82.95K and86.45K|191,the corresponding dimensionless values are0.69 and0.72,respectively(all physical quantities in the table are in reduced units).
3分析与讨论
3.1不同截断半径下氮原子体系NPT系综熔沸点变化
考察了不同截断半径对氨体系熔沸点相图的影响。通常计算物质熔点有滞后法I"。］、两相固液共存法".⑷和自由能法""I等。两相固液共存法通过降低成核自由能的势垒模拟固液共存状态来计算熔点，该方法计算较准确，但模拟温度的扰动需限制在熔点附近。自由能方法计算较复杂，计算量大，尤其液体的自由能很难准确计算。滞后法直接简单，通过模拟类似实验加热或冷却过程计算熔点，但由于连续升温或冷却不能保证体系在较短模拟时间内达到热力学平衡状态，存在过热和过冷的现象（滞后现象），观测到的体系融化和凝固的温度与理论值存在一定偏差。对于加热过程，融化滞后的程度与模拟中使用的原子间作用势函数以及体系的加热速度有关，计算的熔点误差约为0AT m~0.5T^2K加热速率越快，体系过热程度越大，误差越大。为消除加热过程伴随的误差，在每个温度点均采用单独算例进行计算。在外压为一个大气压的NPT系综下，每原子平均总能（势能+动能）随温度的变化曲线如图1所示。每个模拟运行100万步，取最后达到热力学平衡状态的30万步计算单原子总能的平均值。通过一系列不同温度下氨体系热力学平衡状态的计算，根据文献［17,18］给出的判据，获取了体系在不连续跳变点的熔点和沸点信息。可看出，随溫度的上升模拟体系由固体逐渐转变为液体，在熔点时每原子平均总能明显增大，继续升温至沸点时，原子平均总能突然跳变，体积急速膨胀达到气态，完全进入气相状态，模拟过程复现了固、液、气三相转变过程。氮原子系统在不同截断半径下熔沸点的数据对比见表1。截断半径2.5耐，计算的熔点与实验值（82.95K）问的误差约为0.069心（误差由计
算得到的无量纲熔点减去真实的无量纲熔点再除以该真实的无量纲熔点得到），远小于滞后法加热测得的熔点，与文献［20,21］模拟的结果一致，验证了模拟计算的可靠性。截断半径为2.0册熔点的计算误差约为0.09&，显著大于截断半径2.5册计算误差。对于直接加热法，从初始晶格点阵排布的固态融化获取熔点会因过热滞后现象导致计算的熔点偏高，而截断半径2.0耐计算的熔点己显著小于氮的真实熔点，实际计算的熔点应该更低，因此，不建议釆用2.0o M截断半径。
2
-2
E
o
l
e
3d
A
m
u
u
u
水仙花什么时候开
一00
1
-6
-8
0.40.60.8 1.0 1.2 1.4 1.6
Temperature
图1p=101kPa时NPT系综中氮原子系统不同截断半径下
的每原子平均总能随温度的变化
Fig」Variations of the average total energy per atom with temperature at different cutoff distances of the argon system at
/7=101kPa in the NPT enmble
从表1和图1可以看出，不同截断半径下氨系统的熔沸点存在显著变化。增大截断半径，体系达到相变点的温度越高，熔点与实验值偏离越大，且截断半径为4.5cr和5.0耐，熔沸点相变温度相差很小，表明NPT系综在截断半径为4.5c以上时截断半径的影响己经很弱。
262
过程工程学报第21卷
截断半径越大，液态区域所覆盖温度的范围越大，不同 的截断半径导致了模拟的氮体系存在不同的热力学液
态温度范围。
3.2不同截断半径下氮原子体系NPT 系综熔点及液态
热力学状态点的热力学性质
利用LAMMPS 复现了 Anikeenko 等凹在截断半径
2.56统计约化温度0.85下液态氮体系的扩散系数。根
据Lee 〔23］研究的体系规模对扩散系数的影响，用2048个 粒子在NVT 系综下得到体系扩散系数£>=2.385x10-3
nm 2/ps ，与文献［22］—致。另外，复现了 Huang 等卩勺在 高温高压条件下研究的NPT 系综下截断半径对径向分
布函数（Radial Distribution Function, RDF ）的影响的部分 工作，如图2所示。第一、二峰的出现位置与文献［24］
基本一致。但由于模拟体系原子数目和径向统计网格宽
度不同，第一峰的强度存在细微差别，该结果验证了本 模拟工作的可靠性。虽然Huang 等使用了相同温度下不 同截断半径的RDF 对比，但由于不同截断半径会造成
氮原子体系熔沸点相图的变化，此时相同的温度并不代 表二者体系处于相同的热力学状态。例如模拟温度为
0.8，在截断半径为2.5耐体系为液态，截断半径为4.5c
时体系为固体。因此，图2中相同温度下RDF 第一峰
峰高的差异可能是体系处于不同的热力学状态造成的。
3
2
cow ury
U C H 3q E S _p S P E H
assists
4 6 8 10
Radial distance/x 10 1 nm
图2 NPT 系综、7=133 K 、p=20 MPa 液氮体系不同截断半
径下径向分布函数与文献结果a 】的比较
Fig.2 Comparison of RDF between the system and the literature
results 1241 under different cutoff distances of liquid argon system
in NPT enmble of 7M33 K, p=20 MPa
为了分析不同截断半径下液体热力学性质的变化，
以不同截断半径下的熔点为基准温度点，考察熔点及离 熔点相同温度比例的液相区域两个热力学状态点的径 向分布函数和速度自相关函数（Velocity autocorrelation
function, VACF ）（详细参数见表1）。计算了体系在不同截 断半径下各自熔点处的RDF 和VACFo 体系的模拟时长
为100万步（归约时间步长为0.01）, RDF 的计算基于最
后30万步体系的平均结果。在100万步之后使用0.001 的时间步长再运行1万步统计体系的VACF，结果如图
3和4所示，不同截断半径下对应各自熔点的径向分布
函数和速度自相关函数高度一致。
c o 二。unj
uo-lnqunslp
12
69
Radial distance/xlO~ nm
图3 p=101 kPa, NPT 系综液氨体系不同截断半径下各自熔
点径向分布函数的比较
Fig.3 Comparison of the RDF of the respective melting points at different cutoff distances in the NPT enmble at p= 101 kPa
0 5
L
o.a u o _Q u a
uo 二
eRJJOOOlne
U >
0.0 0.5
1.0
Integrating time
图4 p=\0\ kPa,NPT 系综液氮体系不同截断半径下各自熔
点速度自相关函数的比较
营销电话Fig.4 Comparison of VACFs at the respective melting points
with different cutoff distances in the NPT en s emble atp=I01 kPa
香蕉鸡蛋饼虽然不同截断半径造成了体系熔点温度不同，但熔
点处表达的热力学状态相同。RDF 和VACF 的高度一致 表明不同截断半径下各自熔点处体系的平均配位原子 分布相同，液体的结构和热力学性质近似。以此为基础,
分别统计了距离熔点相同温度分率的温度点（r ra +50%Z ）
所对应的热力学状态的径向分布函数和速度自相关函 数，如图5所示。从图中可以看出，在对应热力学状态
一致的状态下，即便截断半径和对应温度不同，体系的
径向分布函数和速度自相关函数未发生根本改变，统计
曲线基本一致。图中出现细微差距的原因可能与热力学
第3期孙晨阳等:LJ 势氨系统分子动力学模拟中截断半径的选择
263
状态相同点的温度精度、统计时域和数值波动等的影响 有关。所以在模拟液体的结构和热力学性质时，不必过
分追求大的截断半径，可考虑以熔点为基准点计算其在 相同的液态温区长度分率的状态点来表达近似相同的
热力学状态(图6，其中Lx 和厶2分别为对应2.50-和4.5<t 截断半径下液态温区的范围，几和兀为每种模拟状况
下的熔点和沸点)，因为此时它们具有近似相同的体系 空间结构稳定性和时间关联性。
o
o o
.5
L 0.U O R O U n J
u o p e a x lo o o ln e
A lo O R A
3
2
o
012
6
3
o (、)"
c O I Q u r y
u o p n q -b s -p
-.s p e E
Radial distance/xlO - nm Integrating time
图5 p=101 kPa, NPT 系综液氮体系不同截断半径下距离熔点50%液态温区长度的温度点（几+50%厶）
的径向分布函数和速度自相关函数的比较
Fig.5 Comparisons of the RDF and VACF at the liquid temperature （Tm+50%厶）under different cutoff distances
in the NPT enmble at p=101 kPa
l
u s e 」u d A e>D
」u u u 一s o l
o 45 o
4 5.5.0
-1.27
a 9
o o
Temperature
狗狗聪明排行榜图6不同截断半径下液氨体系不同热力学状态
点之间的映射
Fig.6 Mapping between different thermodynamic state points
under different cutoff distances
算性能、力计算部分计算能力和邻居列表更新部分计算 性能分别为截断半径为4.5。时的4.43倍、4.71倍和4.06 倍。
S B U n l
p u s d e E
150
)()
50
2.50
3.5
Cutoff distance/cr
3.3不同截断半径下的相同规模算例的计算性能
图7为基于不同截断半径下的一个NPT 系综热力 学状态点(r m +50%£)的计算时长分析，体系规模同表1,
每2步检查一次邻居列表，运行1万步，并统计程序的 总计算时长、力计算耗时、更新邻居列表耗时等。测试
运行环境为 CentOS relea 6.5 (Final) Kernel 2.6.32-
431.el6.x86_64, CPU 为 Intel(R) Xeon(R) CPU E5-2630 v4 (2.20 GHz)，内存为 DDR3 64GB, LAMMPS 使用 icpc
编译加-03优化选项。由图7可知，截断半径越大，计
算耗时越长，计算性能显著下降。当选用截断半径2.5c 时，在不影响RDF 和VACF 计算结果的条件下，总计
图7不同截断半径下液氨体系相同热力学状态点(F m +50%Z.)
模拟计算的性能对比
Fig.7 Performance comparison of simulation calculation for the
same thennodynamic state point (7^+50%厶)of liquid argon
system under different cutoff distances
肉糜蒸蛋4结论
研究了不同截断半径下NPT 系综氨原子系统熔沸 点相图的变化，统计分析了熔点处及液相区域不同热力
学状态点的径向分布函数和速度自相关函数，得到如下 结论：
(1)以不同截断半径下的熔点状态为基准，在距离
