分子间势能Lennard
分子动力学模拟方法的基本原理与应用
分子动力学模拟方法的基本原理与应用摘要: 介绍了分子动力学模拟的基本原理及常用的原子间相互作用势, 如Lennard-Jones势; 论述了几种常用的有限差分算法, 如Verlet算法; 说明了分子动力学模拟的几种系综及感兴趣的宏观统计量的提取。
关键词: 分子动力学模拟; 原子间相互作用势; 有限差分算法;分子学是一门结合物理，和化学的综合技术。
分子学是一套方法，该方法主要是依靠牛顿力学来模拟分子体系的运动，以在由分子体系的不同状态构成的系统中抽取样本，从而计算体系的构型积分，并以构型积分的结果为基础进一步计算体系的量和其他宏观性质。
从统计物理学中衍生出来的分子动力学模拟方法(Molecular Dynamics Simulation, MDS) , 实践证明是一种描述纳米科技研究对象的有效方法, 得到越来越广泛的重视。
所谓分子动力学模拟, 是指对于原子核和电子所构成的多体系统, 用计算机模拟原子核的运动过程, 从而计算系统的结构和性质, 其中每一个原子核被视为在全部其他原子核和电子所提供的经验势场作用下按牛顿定律运动。
它被认为是本世纪以来除理论分析和实验观察之外的第三种科学研究手段, 称之为“计算机实验”手段, 在物理学、化学、生物学和材料科学等许多领域中得到广泛地应用。
科学工作者在长期的科学研究实践中发现，当实验研究方法不能满足研究工作的需求时，用计算机模拟却可以提供实验上尚无法获得或很难获得的重要信息；尽管计算机模拟不能完全取代实验，但可以用来指导实验，并验证某些理论假设，从而促进理论和实验的发展。
特别是在材料形成过程中许多与原子有关的微观细节，在实验中基本上是无法获得的，而在计算机模拟中即可以方便地得到。
这种优点使分子动力学模拟在材料研究中显得非常有吸引力。
分子动力学模拟就是用计算机方法来表示统计力学，作为实验的一个辅助手段。
分子模拟就是对于原子核和电子所构成的多体系统，求解运动方程(如牛顿方程、哈密顿方程或拉格朗日方程)，其中每一个原子核被视为在全部其它原子核和电子作用下运动，通过分析系统中各粒子的受力情况，用经典或量子的方法求解系统中各粒子在某时刻的位置和速度，以确定粒子的运动状态，进而计算系统的结构和性质。
兰纳琼斯势基础知识
记 因此有：
另外一种表达式：
h：两表面间的空气间隙； Wad：粘性力做工。 一般，为了简便，采用下式估算来代替兰纳琼斯势。
C U= r 6 r r0 r r0
这种简化对很多重要参数， 如平衡距离和结合能等影响不大， 却有利于计算。
h
单个原子和由相同原子构成的体之间的相互作用能为：
U a s =-
h
Cn Cn dV 2 r6 6h
可以计算两个固体的平行面的想呼呼作用能，假设每个平面都由相同的原子 构成，这种相互作用能是第一个物体沿 z 轴的积分乘以物体表面积 A 和原子密 度 n 得到的。
h
单位面积的相互作用能表示为：
2 Cn US S Cn2 =- dz 3 h A 6z 12h2
根据范德华力，若中性原子之间的距离大于或者等于其本身原子直径，就会 相互吸引。 两个距离为 r 的中性原子之间的相互作用势 Lennard-Jones 势来描述：
：势能阱深度。 ：是互相作用的势能正好为零时的两体距离。
在实际应用中，ε、σ 参数往往通过拟合已知实验数据或精确量子计算结果而 确定。两个原子之间的相互作用力为：
当 r r0 1010 m 时，分子间排斥力和吸引力平衡。分子间相互作用力为零，此 时成为平衡位置。
经过推导，两个原子之间的相互作用势能可以表示为两个表面之间的相互牵 引力。由兰纳琼斯势引起的表面作用力表示为
其中： z：两个接触平面面之间的距离。 当两表面之间的作用力为 0 时，其平衡距离 h0 为：
如果将两个物体从很远的地方移动到一起“直接接触”（接近 r0=10-10m）,则单 位面积相互作用力做的功为:
2 Cn W Cn 2 =- dz h A 6z3 12r02
第五章_固体表面吸附(固-气界面吸附)
③化学吸附具有选择性
如：CO在金属表面上的吸附
O C MM
OO CC MM
NO + 1/2O2 CH2=CH2 + 1/2O2
NO2 O
NO只在Pt上吸附 只用Ag作催化剂。
无选择性；吸附热与气体凝聚热相近；吸附速度快；多层吸附。
（2）相互作用势能
设 f 代表作用力，r 为粒子间距离，U(r)为粒子间相互作用势能，有：
f U (r) r
或
r
U(r)
f dr
永久偶极矩相互作用势能为：
U (r)
2 3
12 22
3k Tr 6
诱导偶极矩相互作用势能为：
Ui(r)
( i 2 12
使体相中某些组分在表面区产生富集的现象。 其特点为组成随表面吸附质不同而变化。
如：Ag-Pd合金，吸附CO时，体相中的Pd可通过扩散到达 表面与CO形成羰基键，从而使表面富Pd。除去CO后，表面 组成由回到原来的状态。
第二节 物理吸附和化学吸附
1、物理吸附与Lennard-Jone势能曲线
（1）物理吸附特点
S
RT
则：
p
或： bp
b(1 )
1 bp
若有两种气体存于表面而发生竞争吸附时：
则：
A
பைடு நூலகம்
bA pA 1 bA pA bB
pB
B
bB pB 1 bA pA bB
pB
如果吸附解离成两个碎片，且各占据一个吸附位置，则：
f ( ) (1 )2 f '( ) 2
势能函数 截止半径 rtap rcut
势能函数截止半径rtap rcut全文共四篇示例，供读者参考第一篇示例：势能函数是描述物体系统在不同位置之间相互作用的一种数学模型，它可以帮助科学家预测系统的行为和性质。
在许多物理学和化学领域中，势能函数都起着至关重要的作用，可以帮助研究人员深入了解分子、原子和其他微观粒子之间的相互作用。
截止半径以及截止半径内以及截至半径外的作用力的计算是非常重要的参数。
在分子模拟和计算化学中，截止半径是指在势能函数中相互作用的距离上限。
当两个粒子之间的距离超过截止半径时，它们之间的相互作用会被认为是零。
这个参数在计算中起着关键作用，可以减少计算的复杂性，同时也能够更好地模拟实际系统的行为。
另一个重要的参数是rcut，它表示在模拟中考虑的最大距离。
当两个粒子之间的距离超过rcut 时，它们之间的相互作用也会被截断。
这个参数的选择取决于所研究的系统和模拟的精度要求。
通常情况下，rcut 应该被选取得足够大，以保证系统中的所有相互作用都被考虑在内，并且不会对计算结果产生显著的误差。
势能函数中的截止半径和rcut 参数通常根据所研究的系统的特性和模拟的需要来选择。
在分子动力学模拟中，这些参数的选择对系统的稳定性和计算结果的准确性有着重要的影响。
一般来说，这些参数的选择应该是一个平衡，既要考虑计算的效率，也要考虑结果的准确性以及模拟的真实性。
在实际的研究中，科学家们会根据所研究系统的特性和模拟的需求来选择截止半径和rcut 参数。
通过不断调整这些参数，他们可以优化模拟的效率和准确性，从而更好地理解和预测系统的行为。
这些参数的选择也需要结合实验数据和理论模型，以确保模拟结果与实际情况相符。
第二篇示例：势能函数在物理学中是一个非常重要的概念，它描述了粒子或系统在不同位置的势能大小。
在化学和材料科学中，势能函数常常被用来描述原子之间的相互作用。
在许多模拟和计算中，势能函数的截止半径(rt)和截止距离(rcut)是两个关键参数，它们决定了系统中粒子之间的相互作用范围。
5_分子间相互作用势
分子间相互作用势
3.3 硬球＋吸引模型
要模拟实际流体的性质，必须包含吸引作用，一个简单的 方法是在硬球模型的基础上再加上吸引项。最简单也是实际中 应用非常多的是方阱势能模型：
u(r)
0
r r r
是吸引作用强弱的度量，称为阱深， 称为阱宽，其
值一般取1.5。
分子间相互作用势
3.3 硬球＋吸引模型
r r
L-J势能函数是计算机模拟中使用最广泛的分子间相互作
用模型。有时截断的L-J势能也被用于分子间的排斥作用：
分子间相互作用势
常用的L-J势能函数：
E(r)
4
12
6
r r
分子间相互作用势
截断的L-J势能模型：
4
12
6
u(r) r r
0
r 21/6 r 21/6
rm m / n 1/(nm)
分子间相互作用势
3.4 非极性分子的排斥＋色散模型
m和n通常取为6和12，则L-J势能函数可表示为：
rm m / n 1/( nm ) rm 6 / 12 1/(126 ) rm 1/ 2 1/ 6
rm 21/ 6
u(r)
4
12
6
u(r) Aexp(Br) C6 / r6 u(r) bexp(r / a) C6 / r6 C8 / r8
u(r)
b
ex
p
r
rm
C6 r6
C8 r8
exp
4
rm r
13
b
ex
p
r
rm
C6 r6
C8 r8
r rm r rm
分子间相互作用势
3.5 静电作用和偶极作用
范德华力的公式
范德华力的公式范德华力是分子间存在的一种较弱的相互作用，虽然它不像化学键那样强而有力，但在很多物质的性质和现象中却起着不容忽视的作用。
要了解范德华力，咱们就得先聊聊它的公式。
范德华力的公式主要有两个，一个是描述吸引作用的伦敦色散力公式，另一个是描述总体范德华力的 Lennard-Jones 势能公式。
先来说说伦敦色散力公式。
它的表达式是：$F = - \frac{3}{2} \alpha I_1 I_2 / r^6$ 。
这里的 $\alpha$ 是分子的极化率，$I_1$ 和 $I_2$ 分别是两个相互作用分子的电离能，$r$ 是分子间的距离。
Lennard-Jones 势能公式则看起来稍微复杂点：$U(r) = 4\epsilon\left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$ 。
这里的 $\epsilon$ 表示势能阱的深度，$\sigma$ 是分子间距离为平衡距离时的值。
说起来，我曾经在课堂上讲这个公式的时候，发生过一件特别有趣的事儿。
当时我在黑板上写下了这些公式，下面的学生们一个个都瞪大了眼睛，满脸的困惑。
我就问：“同学们，你们觉得这公式像不像一个神秘的密码？” 有个调皮的学生大声说：“老师，这密码太难解啦！”全班哄堂大笑。
其实啊，这些公式看起来复杂，但只要咱们一点点拆解，理解其中每个符号的含义，就会发现它们也没那么可怕。
就拿极化率 $\alpha$ 来说，它反映了分子在外界电场作用下变形的难易程度。
分子越大、越容易变形，极化率就越大，产生的伦敦色散力也就越强。
好比一个大大的气球，轻轻一戳就变形了，而一个小小的乒乓球，可就没那么容易被改变形状啦。
再看电离能 $I$ ，它反映了原子或分子失去电子的难易程度。
电离能越大，说明原子或分子越不愿意失去电子，在范德华力的作用中就相对不那么活跃。
lammps模拟时金属原子的l-j势参数提取方法
lammps模拟时金属原子的l-j势参数提取方法文章标题：如何在LAMMPS模拟中提取金属原子的L-J势参数在分子动力学模拟中，LAMMPS（Large-scale Atomic/Molecular Massively Parallel Simulator）是一个常用的开源软件工具，用于模拟原子体系的运动和相互作用。
在对金属原子进行模拟时，提取准确的Lennard-Jones（L-J）势参数对于准确模拟金属原子的相互作用至关重要。
本文将探讨在LAMMPS模拟中提取金属原子的L-J势参数的方法和步骤。
1. 理解Lennard-Jones势参数的概念Lennard-Jones势能函数是描述原子或分子间相互作用的经典模型之一，通常用于描述非键相互作用。
该势能函数包含两个部分，一个描述范德华力吸引力，另一个描述范德华力斥力。
在LAMMPS模拟中，通过调整L-J势参数的数值，可以准确描述金属原子之间的相互作用，并得到比较真实的模拟结果。
2. 原子结构模型的建立在提取金属原子的L-J势参数之前，首先需要建立金属原子的原子结构模型。
根据金属类型的不同，可以采用不同的结构模型，如面心立方（FCC）、体心立方（BCC）等。
根据原子的晶格常数、结构类型等参数，可以在LAMMPS中建立金属原子的模拟体系。
3. 模拟体系的能量最小化在模拟体系建立完成后，需要进行能量最小化以达到体系的稳定状态。
通过LAMMPS中提供的能量最小化算法，可以使金属原子的位置在势能面上找到最稳定的状态，从而为后续的L-J势参数提取做好准备。
4. 提取L-J势参数在模拟体系达到稳定状态后，可以通过LAMMPS提供的计算工具，如pair_style命令和fitting工具，来提取金属原子的L-J势参数。
通过对模拟体系的势能曲线进行拟合，可以得到最优的L-J势参数数值。
需要注意的是，不同金属原子之间的相互作用可能需要不同的L-J势参数，因此针对不同金属类型需要分别进行提取。
分子动力学模拟常用基本概念(相关碳纳米管的概念为主)
分子动力学模拟常用基本概念1、势函数： （1）Tersoff 势：Tersoff 势起源于对C 原子的处理方法，是一种共价键类型的原子间作用势，它不仅可以计算相应晶格常数、键能、键角、弹性模量和空位形成能，和其它力场模型相比，可以描述系统中化学键的形成和断裂以及原子之间化合键变化的动态过程。
Tersoff 势可以很好表述碳氢分子、石墨、金刚石间相互作用能、键能，可以表示化学键的断裂和形成，比如计算金刚石C 11、C 12、C 44的弹性常数和实验结果比较接近。
通过它可对系统进行分子动力学模拟，可以计算系统中的化学键键长、键能、键角、弹性模量和空位形成能。
Tersoff 势函数被广泛用于讨论碳纳米管的稳定结构、形成机理、力学性能以及碳纳米管中碳原子的一些动态过程。
Tersoff 势成功地被用来描述石墨、金刚石的碳键相互作用。
碳纳米管中碳原子间共价键的相互作用较广泛地采用Tersoff 势来描述并取得非常大的成功。
Tersoff 势被认为是键合强度依赖于周围原子配置的势函数，可以很好的描述表面重构能，能比较好地描述碳纳米管性质而被广泛应用。
Tersoff 势总能量函数形式为：[()()]c ij r ij ij a ij ii jf a E r b E r <Φ=-∑∑其中：排斥势：()exp()r ij ij ij ij E r A r λ=-； 吸引势：()exp()a ij ij ij ij E r B r μ=-12(1)i i innn ij ij i i a εβτ-=+；2(1)i iiim n nn ij ij i i b χβξ-=+,()()ij c ik ik ijk k i jf rg τδθ≠=∑；,()()exp[()]ij c ik ik ijk ik ij ik k i jf rg r r ξϖθσ≠=-∑角函数：22222()1(cos )iiijk ii i ijk c c g d d h θθ=+-+-截断函数：11()[1cos()]20ij ij c ik ij ijr R f r S R π⎧⎪-⎪=+⎨-⎪⎪⎩式中，αij 是截断距离，一般情况下，必须将αij 式中的β的值取得充分小，使得αij ≈1，因为在第一临近之外的范围内，τij 会指数式地变大。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
容和弹性系数的表达式。波尔兹曼的工作源于他对于简单的麦克斯韦相互作用定律的怀疑，
他比较喜欢引力，所以波尔兹曼用了不同的吸引模型重复了全部的计算，也可以得到与麦克
斯韦相近的结果。
中性原子、分子间的长程引力，后来被称为范德华（van der Waals)力。这个主要是来源
于他的著名状态方程，这个方程说明了实际气体偏离理想气体的行为：
式。用实验数据来拟合势参数的半经验方法被证明是非常有用的，并在 20 世纪得到了继续
的发展。其中开创性的工作来自于 Lennard-Jones 的论文中。这就是众所周知的 Lennard-Jones
势，已在气体和凝聚态物质的研究中得到了广泛的应用。
二、 Lennard-Jones 势
1. 势能的基本描述
3. Kihara 势 在许多情况下，分子的尺寸是不可忽视的。Kihara 在考虑了分子尺寸的基础上，对
Lennard-Jones 势做了一定的更改。根据他的设想，每一个分子都可以看做是具有凸面的刚 性转子，而在势函数表达式中的分子间距用分子与分子面之间的距离来代替，得到如下势能
函数：
V(R) = ϵ [(������������0)12 − 2 (������������0)6]
������������������������ ������������
−
������������������������我们可以
������7
看到他选择了不同的 n 的取值，用数学上的回归方法进行拟合，分别得到了以下的数据
从表格中可以看出，当 n=13 时，用公式所 得到的最近的原子距离和最低的分子势能与实验
分子间势能 Lennard-Jones 势函数研究进展
一、 分子间相互作用研究背景
粒子间相互作用定律最早是由克罗地亚物理学家 Boscovich(1711~1787)提出的，他认为 所有的物体都由点粒子组成。任何两点间的振动力随着相互接近而无限上升，并且在长距离
时符合牛顿万有引力定律的渐进行为（
）（如右图 1 所示）。
−
������������������������ 6������6
拟合之后的结果如左图：
虽然只是数学上拟合的结果，但是由于这个
势能函数的简单，而且比较实用，而得到了人们
的广泛认可，在物理化学领域得到了广泛的应用。
而且上式中的排斥项的指数 n 则可根据实际问题 的需要进行变动，通常范围为 n=12~25.
������ 12 ������ 6 V(R) = 4ϵ [(������) − (������) ] 式中，ϵ为当势函数在极小值点������������ = 21/6������时的数值，而������为势能零点的R值（如图 2 所 示）。势函数中的吸引项对应于偶极-偶极的相互作用。而排斥的一项也用幂函数进行近似。
−
������
(������������������
8
)]
Maitland 和 Smith 对上式进行了修改。他们所取的势函数与上式类似，但排斥项的指数
与 R 有如下关系：
������ = ������ + ������ (������������������ − 1) 实践也证明（m-6-8）势具有良好的可调节性，可用一套参数来描述多种单原子稀有气 体的平衡和输运性质。
式中有 4 个参数：m ,A ,B 和 C。Klein 和 Hanley 也提出另外一种方法来表示上式
V(R)
=
ϵ
6 + 2������ [( ������ − 6 )
(������������������)������
−
������
− ������(������ − ������ − 6
8)
(������������������)6
而������ = 0时，变为（12-4）势。
2. （m-6-8）势 Klein 和 Hanley 提出将偶极-四极相互作用项加入 Lennard-Jones 势，并将排斥势的幂次
设为可变的参量，这就得到了 Klein-Hanley 势：
������ ������ ������ ������(������) = ������������ − ������6 − ������8
Lennard-Jones 势的一般表示形式可表示为：
V(R)
=
������������ ������������
−
������������ ������������
他最早是 Lennard-Jones 在研究稀有气体性质，尤其是计算其维里系数的时候提出的。
后来，该模型势在各种体系的计算中得到了广泛的应用。人们常说的（12-6）势可写为：
积分是有限的，当 R→∞时，这个作用势具有 1/R^n 的形式，其中 n 不小于 6.
在 19 世纪下半叶，克劳修斯、麦克斯韦和波尔兹曼发展了气体分子的运动论。克劳修
斯接受了分子近距离相互排斥、远距离相互吸引的观点。而麦克斯韦认为分子间作用力在任
何情况下都完全是斥力。他用相互作用势 V(R)=A/R^n 来表示斥力，并得到了扩散系数、热
1
1
{������������ = √2 (������1 + ������2) , ������������ = √2 (������1 − ������2)
使哈密顿量对角化：
H
=
1 [2������
������
(������
−
2������2 ������3 )
���������2��� ]
������ 12
������ 6
������ 4
V(R) = 2ϵ [(1 + ������) (������) − 2������ (������) − 3(1 − ������) (������) ]
式中，������用来描述������−6项对体系贡献的大小。如果������ = 1，（12-6-4）势变为（12-6）势，
用势的参数可以通过下公式从同分子相互作用势得出
������0������������
=
1 2
(������0������������
+
������0������������)
������ ������������ = ( ������ ������������ ������������������)1/2
2. 参数 m=6 的确定 Lennard-Jones 势的后面一项中代表的是原子间的吸引势能，现在我们只是考虑两个原
子之间的相互作用力，不考虑两个原子之间的相互关联，考虑一个比较简单的模型，如下图
3 所示。考虑两个相距为 r 的全同线性谐振子，每个振子具有电荷±e，相距为x1和x2。当两 个振子没有相互作用时，由此可得系统的哈密顿量为
按照 Boscovich 的论文中，为了解决弹性物质的形变以及气体的 物理性质，他提出了，在分子相互作用的体系中必然同时存在 引力和斥力，而且二者可以交替发挥作用。
几乎与 Boscovich 同时，法国物理学家 Clairault 于 1743 年
提出分子间相互作用力的概念，由于解释液体在毛细管中上升。
因为|������1|，|������2| ≪ r，将其展开至二次项有
H1
≈
−
2������2������1������2 ������3
系统总的哈密顿量为
引入简正坐标
H = H0 + H1
1
1
������������ = √2 (������1 + ������2) , ������������ = √2 (������1 − ������2)
������ (������ + ������2) (������ − ������) = ������������
在随后的研究中，人们用了各种不同的经验势来解释气体和液体的的物理性质。
sutherland 在一系列论文中检验了气体分子间吸引力的某些解析形式，特别是 V(R)-A/R^3，
其中的参数通过实验数据进行拟合。这种唯象方法在该领域的后续研究中成为一种典型的方
=
1 2
ℏ(ω������
+
ω������)。没有考虑相互作用前E0
=
1 2
ℏ������0
其中������0 = (������/������)1/2
∆������
=
������
−
������0
≈
−ℏ������0
∙
1 8
2������2 2 (������������3)
=
−
������������ ������6
H0
=
1 2������
������12
+
1 2
������������12
+
1 2������
������22
+
1 2
������������22
当两个振子相互作用的时候，他们的库仑相互作用势为
������2
������2
������2
������2
H1 = ������ + ������ − ������1 + ������2 − ������ + ������2 − ������ − ������1
同分子二聚体的 Kihara 势函数的参数值
分子
������0/������0
H2
2.62
N2
3.38
O2
3.12
F2
2.85
CO2
3.20
CH4
3.10
CF4
3.10
(ϵ/k)/������ 48.5 128 160 158 337 224 298