使用Python代码实现一个简单的分子动力学模拟程序

1. 前言

理解分子动力学模拟最好的方法是编写一个分子动力学程序，换句话说，教会计算机去理解并做执行，自己才算理解会了。因此本文将从常用于描述分子间的非键相互作用中的Lennard-Jones potential讲起，最后将其应用到一个简单的分子动力程序中，来模拟二维空间下流体的分子动力学行为，为了更加直观的理解和展示这一过程，本文主要是Python代码实现相关的代码和程序，并使用jupyter notebook来进行编码展示和描述，会尽可能的解释出每一行代码的含义，让分子动力学模拟的初学者和爱好者能从0到1编写出一个初级的分子动力学小程序，以加深对分子动力学模拟过程的理解。

2. 分子间相互作用力

2.1 简介

根据现代物理力学理论可知，力的本质其实就是宇宙中存在的四种相互作用，按照从强到弱的顺序，排在第一的是强相互作用，其可以使原子核里面的质子，中子稳定地存在。排在第二的是电磁相互作用，它贡献了化学键，分子间相互作用力，以及各种宏观世界的的力，比如表面张力，静电力，磁力等等。 排在第三的是弱相互作用，其也只有在原子核内有效，而且强度非常弱，虽然没有什么存在感，但是很重要，排在最后就是万有引力了，其强度极弱，对物质的性质没有可观察的影响，可以忽略不计了。

在这里主要介绍由电磁相互作用引起的分子间相互作用力，亦称分子间引力，是介导分子间相互作用的力，包括作用于原子和其他类型的相邻粒子（例如原子或离子）之间的吸引力或排斥力。分子间相互作用力(非键相互作用) 相对于分子内相互作用力（键相互作用） 来说是微弱的。例如，涉及原子之间共享电子对的共价键比相邻分子之间存在的力强得多。这两组力都是分子力学中常用力场的重要组成部分。

分子间相互作用力，它主要包括：

范德华力（Van der Waals force）：起初为了修正范德华方程而提出。普遍存在于固、液、气态任何微粒之间，与距离六次方成反比。
次级键（secondary bond）： 键长长于共价键、离子键、金属键，但短于范德华相互作用的微观粒子相互作用。

1. 氢键（Hydrogen bonding）：氢与氮、氧、氟所键结产生的作用力。
2. 非金属原子间次级键：存在于碘单质晶体中。
3. 金属原子与非金属原子间次级键：存在于金属配合物中。
4. 亲金作用。
5. 亲银作用。
   此外科学家也不断研究新的分子间作用力，包括双氢键和金键等。

2.2 范德华力

范德华力（英语：van der Waals’ Forces；又名:范德瓦耳斯力）在化学中指分子或非极性原子（稀有气体）之间非定向的、无饱和性的、较弱的相互作用力，包括排斥力和吸引力。根据荷兰物理学家van der Waals（约翰内斯·范德瓦耳斯）命名，其于1910获得Noble物理奖。范德华力是分别来自两个分子中的一种中性原子之间通过瞬间静电相互作用产生的一种弱的分子之间的力，这是由于附近粒子波动极化的相关性引起的（量子动力学的结构），具体而言，电子密度可能会暂时更大地转移到原子核的一侧。这会产生一个瞬态电荷，附近的原子可以被吸引或排斥。但它比化学键或共价键弱得多，通常其能量小于5kJ/mol（0.4~4kJ/mol）。当两个原子之间的距离为它们的范德华半径之和时，范德华引力最强。强的范德华排斥作用可以防止原子相互靠近。范德华力的大小和分子的大小成正比。

在介绍范德华力的来源前，先简单了解几个概念：

固有偶极：极性分子中由于组成元素不同，其吸引电子的能力各有差异（元素周期律），这就使得分子中有电子偏移的现象，这样就产生了极性，并且偶极持续存在，称为固有偶极。

诱导偶极：是指非极性分子在电场中或者有其他极性分子在较近距离的情况下，由于电子带负电，核带正电，它们会发生偏移，这种现象称为诱导偶极。

瞬时偶极：一切分子中，不管是极性分子还是非极性分子，原子核时刻在震动，电子时刻在运动、跃迁，在它们运动的时候，偶尔离开平衡位置而产生极性，只不过这个过程持续时间很短，故称瞬时偶极。

根据范德华力的来源可分为以下几个部分：

1.取向力（1912-Keesom force）：固有偶极之间的电性引力；(热平均固有偶极——固有偶极相互作用)

范德华力中吸引力的第一个贡献项是由于固有偶极子、四极子（所有对称性低于立方的分子）和多极子之间的静电相互作用。它被称为Keesom相互作用，以Willem Hendrik Keesom 的名字命名。这些力源于固有偶极子（偶极分子）之间的吸引力，并且与温度有关。
Fig 1.取向力示意图

它们由集合的偶极子之间的吸引力相互作用组成对偶极子的不同旋转方向进行平均。假设分子不断旋转并且永远不会被锁定到位。这是一个很好的假设，但在某些时候分子确实会被锁定到位。Keesom 相互作用的能量取决于距离的倒数六次方，而两个空间固定偶极子的相互作用能量取决于距离的倒数三次方。Keesom 相互作用只能发生在具有固有偶极矩的分子之间，即两个极性分子之间。Keesom 相互作用也是非常弱的范德华相互作用，不会发生在含有电解质的水溶液中。具有固有偶极矩${\mu }_i$和${\mu }_j$的分子i和j，其相互作用势能为：

$$E_e=-\frac{2}{3}\left[\frac{{\mu }_i^2{\mu }_j^2}{(4\pi {\epsilon }_0{)}^2{K}_{B}T}\frac{1}{r^6}\right]\text{\hspace{1em}(1)}$$
$E_e$表示随机取向的偶极相互作用的平均作用能

${\epsilon }_0$ :自由空间的介电常数

$T$: 温度

$K_B$:玻尔兹曼常数，

$r$:分子间距离。

2.诱导力（1921 Debye force）：诱导偶极与固有偶极之间的电性引力，如离子诱导偶极力和偶极诱导偶极力。(热平均固有偶极——诱导偶极相互作用)

第二个贡献为项称为感应力（极化）力或德拜力，这些力以荷裔美国物理化学家彼得德拜的名字命名。为具有极化率$\alpha$​​​的分子在周围分子固有偶极矩电场的作用下产生诱导偶极矩，诱导偶极矩与周围分子固有偶极矩的相互作用。当一个具有固有偶极子的分子排斥另一个分子的电子时，就会发生这些诱导偶极子。具有永久偶极子的分子可以在相似的相邻分子中诱导偶极子并引起相互吸引。原子之间不会发生德拜力。因为一般单原子并非分子，无极化率，而特殊当原子分子如惰性气体等，为非极性分子无固有偶极，它们之间没法互相诱导，故德拜力不会产生与原子之间。诱导偶极子和固有偶极子之间的力不像 Keesom 相互作用那样依赖于温度，因为诱导偶极子可以围绕极性分子自由移动和旋转。其相互作用的势能表示为：

$$E_p=-\frac{{\alpha }_i{\mu }_j^2+{\alpha }_j{\mu }_i^2}{4\pi {\epsilon }_0}\frac{1}{r^6}\text{\hspace{1em}(2)}$$

当具有永久偶极子的分子（例如 HCN）与没有分子偶极子的分子碰撞时，碰撞本身会通过分子内电子密度的变化导致偶极子出现。HCN 中的氮原子富含电子，分子偶极子指向该原子的方向。碰撞时，第二个原子的电子云会被氮上多余的电子密度排斥，因此带正电的原子核会更靠近 N 并与其相互作用。

Fig 2.诱导力力示意图

3.色散力（London force）：瞬时偶极之间的电性引力；(诱导偶极——诱导偶极相互作用)

第三个也是主要的贡献是色散或伦敦力（波动偶极诱导偶极），由于分子中电子和原子核不停地运动,非极性分子的电子云的分布呈现有涨有落的状态,从而使它与原子核之间出现瞬时相对位移,产生了瞬时偶极,分子也因而发生变形。分子中电子数愈多、原子数愈多、原子半径愈大,分子愈易变形.瞬时偶极可使其相邻的另一非极性分子产生瞬时诱导偶极,且两个瞬时偶极总采取异极相邻状态,这种随时产生的分子瞬时偶极间的作用力为色散力(因其作用能表达式与光的色散公式相似而得名)。瞬间偶极距随时间不断变化，统计平均为0，因此实验测量得到的分子偶极距离任然为0，色散力作为一种相互吸引力，可以使得体系能量降低，一般地，极化的体积分别为${\alpha }_i$和${\alpha }_j$的两个分子或离子间色散能可以写成：

$$E_d=-\frac{3}{2}\left[\frac{I_i{I}_j}{I_i+I_j}\right]\left[\frac{{\alpha }_i{\alpha }_j}{(4\pi {\epsilon }_0{)}^2}\frac{1}{r^6}\right]\text{\hspace{1em}(3)}$$

虽然瞬时偶极存在暂短,但异极相邻状态却此起彼伏,不断重复,因此分子间始终存在着色散力.无疑,色散力不仅存在于非极性分子间,也存在于极性分子间以及极性与非极性分子间。色散力存在于一切分子之间.色散力与分子的变形性有关,变形性越强越易被极化,色散力也越强。稀有气体分子间并不生成化学键,但当它们相互接近时,可以液化并放出能量,就是色散力存在的证明.

所有分子通过伦敦色散力或诱导偶极相互作用。在下图中，一个2个原子的分子与一个3个原子的分子碰撞。第一个分子的电子云排斥它撞击的分子的电子云，导致一些电子密度远离原子核。原子核然后被它自己的电子屏蔽得很差，并吸引了第一个分子的电子云。

Fig 3.色散力示意图

4.由泡利不相容原理产生的排斥成分，可防止分子坍塌。
上面讨论的分子间相互作用均为吸引作用,并且势能随分子间距离的减小而降低。但是,这些关系并不适合分子相互靠得特别近的情况。当两个分子相互靠近,达到或接近电子云相互重叠的距离时,分子间的排斥作用将显著增加,引起体系势能的迅速升高。此外,原子核间的静电排斥作用,也引起体系的势能升高。完整的分子间相互作用函数必须同时考虑排斥作用的贡献。理论分析表明, 排斥作用对体系势函数的贡献与分子间的距离呈指数关系：
$$u(r)=Aexp(-\beta r)\text{\hspace{1em}(4)}$$

为了计算方便和效率，常把排斥能表示为幂函数的形式：

$$u(r)=A/r^n$$

式中，A为参数;n则在8~16中变化。同时排斥力也是一会走近程相互作用。随分子间距离增大而迅速衰减。当原子间距离低于0.4 nm 时，此时的对范德华力的贡献主要来自泡利不相容原理，让物质内部有所排斥，表现为排斥力。

Fig 4. 各分子范德华力吸引力的分配情况

可以说范德华力属于库伦相互作用（静电相互作用），但是我们经常看到分子间相互作用大体分为两种：静电相互作用和范德华力相互作用（广义的）。前面说到的范德华力属于静电相互作用，其为狭义的范德华力它通常指伦敦色散力（London dispersion force），它的相互作用能$E_d$与分子间距（r）是$\frac{1}{r^6}$的关系，而广义的范德华力，是因为Keesom force和Debye force在热平均后都与London force一样具有$\frac{1}{r^6}$的关系。如果没有热平均，偶极-偶极相互作用是$\frac{1}{r^3}$的关系也就是说Keesom force要更加短程，而相对短程也是范德华力（短程相互作用）区别于静电相互作用（长程相互作用）的特征。因此可以看出非键相互作用分类依据为两者相互作用随距离变化的衰减特点不同。

范德华力的主要特点是：

• 它们比普通的共价键和离子键弱。

• 范德华力是可加的，不能饱和。

• 它们没有方向特性。

• 它们都是短程力，因此只需要考虑最近粒子（而不是所有粒子）之间的相互作用。分子距离越近，范德华引力越大。

• 除了偶极-偶极相互作用之外，范德华力与温度无关。

• 在低分子量醇中，其极性羟基的氢键性质支配了其他较弱的范德华相互作用。在较高分子量的醇中，非极性烃链的性质占主导地位并决定了它们的溶解度。

2.2.1 Lennard-Jones potential 模型

Lennard-Jones potential,它以John Lennard-Jones.（约翰·伦纳德-琼斯）的名字命名被称为LJ势或者12-6势。它描述了简单原子和分子之间相互作用的基本特征：两个相互作用的粒子在非常近的距离相互排斥，在中等距离相互吸引，在无限距离上不相互作用，LJ势是一对势，即该势不包含三体或多体相互作用。LJ势是研究最广泛和最彻底的势，被认为是简单而现实的分子键相互作用的半经验模型。LJ势的常用表达式为：

$$E_{LJ}=4\epsilon \left[(\frac{\sigma }{r_{ij}}{)}^{12}-(\frac{\sigma }{r_{ij}}{)}^6\right]\text{\hspace{1em}(5)}$$

$E_{LJ}$:为两体之间的势能大小。

$\varepsilon $ :为控制相互作用的强度（以 eV 为单位）,从本质上讲，它衡量两体相互吸引的强度,通常称为势阱深度（也叫“色散能”）;

$\sigma$:为相互作用的势能正好为零时的两原子间的距离，$\sigma$给出了两个非键合粒子的接近程度的度量，因此被称为范德瓦尔斯半径。它等于非成键粒子之间核间距离的二分之一。

$r_{ij}$: 为i,j两体之间的距离。

在实际应用中，$\epsilon$,$\sigma$参数往往通过拟合已知实验数据或精确量子计算结果而确定。 另一种写法是：

$$E_{LJ}(r_{ij})=ϵ[(\frac{r_{min}}{r_{i}j}{)}^{12}-2(\frac{r_{}}{r_{ij}}$$