使⽤Python代码实现⼀个简单的分⼦动⼒学模拟程序
1. 前⾔
理解分⼦动⼒学模拟最好的⽅法是编写⼀个分⼦动⼒学程序，换句话说，教会计算机去理解并做执⾏，⾃⼰才算理解会了。因此本⽂将从常⽤于描述分⼦间的⾮键相互作⽤中的Lennard-Jones potential讲起，最后将其应⽤到⼀个简单的分⼦动⼒程序中，来模拟⼆维空间下流体的分⼦动⼒学⾏为，为了更加直观的理解和展⽰这⼀过程，本⽂主要是Python代码实现相关的代码和程序，并使⽤jupyter notebook来进⾏编码展⽰和描述，会尽可能的解释出每⼀⾏代码的含义，让分⼦动⼒学模拟的初学者和爱好者能从0到1编写出⼀个初级的分⼦动⼒学⼩程序，以加深对分⼦动⼒学模拟过程的理解。
2. 分⼦间相互作⽤⼒
2.1 简介
根据现代物理⼒学理论可知，⼒的本质其实就是宇宙中存在的四种相互作⽤，按照从强到弱的顺序，排在第⼀的是强相互作⽤，其可以使原⼦核⾥⾯的质⼦，中⼦稳定地存在。排在第⼆的是电磁相互作⽤，它贡献了化学键，分⼦间相互作⽤⼒，以及各种宏观世界的的⼒，⽐如表⾯张⼒，静电⼒，磁⼒等等。 排在第三的是弱相互作⽤，其也只有在原⼦核内有效，⽽且强度⾮常弱，虽然没有什么存在感，但是很重要，排在最后就是万有引⼒了，其强度极弱，对物质的性质没有可观察的影响，可以忽略不计了。
在这⾥主要介绍由电磁相互作⽤引起的分⼦间相互作⽤⼒，亦称分⼦间引⼒，是介导分⼦间相互作⽤的⼒，包括作⽤于原⼦和其他类型的相邻粒⼦（例如原⼦或离⼦）之间的吸引⼒或排斥⼒。分⼦间相互作⽤⼒(⾮键相互作⽤) 相对于分⼦内相互作⽤⼒（键相互作⽤） 来说是微弱的。例如，涉及原⼦之间共享电⼦对的共价键⽐相邻分⼦之间存在的⼒强得多。这两组⼒都是分⼦⼒学中常⽤⼒场的重要组成部分。
分⼦间相互作⽤⼒，它主要包括：
范德华⼒（Van der Waals force）：起初为了修正范德华⽅程⽽提出。普遍存在于固、液、⽓态任何微粒之间，与距离六次⽅成反⽐。次级键（secondary bond）： 键长长于共价键、离⼦键、⾦属键，但短于范德华相互作⽤的微观粒⼦相互作⽤。
1. 氢键（Hydrogen bonding）：氢与氮、氧、氟所键结产⽣的作⽤⼒。
2. ⾮⾦属原⼦间次级键：存在于碘单质晶体中。python新手代码及作用
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 相互作⽤也是⾮常弱的范德华相互作⽤，不会发⽣在含有电解质的⽔溶液中。具有固有偶极矩和的分⼦i和j，其相互作⽤势能为：
表⽰随机取向的偶极相互作⽤的平均作⽤能
:⾃由空间的介电常数
: 温度
：玻尔兹曼常数，
：分⼦间距离。
2.诱导⼒（1921 Debye force）：诱导偶极与固有偶极之间的电性引⼒，如离⼦诱导偶极⼒和偶极诱导偶极⼒。(热平均固有偶极——诱导偶极相互作⽤)
第⼆个贡献为项称为感应⼒（极化）⼒或德拜⼒，这些⼒以荷裔美国物理化学家彼得德拜的名字命名。为具有极化率的分⼦在周围分⼦固有偶极矩电场的作⽤下产⽣诱导偶极矩，诱导偶极矩与周围分⼦固有偶极矩的相互作⽤。当⼀个具有固有偶极⼦的分⼦排斥另⼀个分⼦的电⼦时，就会发⽣这些诱导偶极⼦。具有永久偶极⼦的分⼦可以在相似的相邻分⼦中诱导偶极⼦并引起相互吸引。原⼦之间不会发⽣德拜⼒。因为⼀般单原⼦并⾮分⼦，⽆极化率，⽽特殊当原⼦分⼦如惰性⽓体等，为⾮极性分⼦⽆固有偶极，它们之间没法互相诱导，故德拜⼒不会产⽣与原⼦之间。诱导偶极⼦和固有偶极⼦之间的⼒不像 Keesom 相互作⽤那样依赖于温度，因为诱导偶极⼦可以围绕极性分⼦⾃由移动和旋转。其相互作⽤的势能表⽰为：
μi μj E =e −32[(4πε)K T 02B μμi 2j 2r 61
](1)
E e ε0T K B r αE =p −4πε0αμ+αμi j 2j i 2r 61
(2)
当具有永久偶极⼦的分⼦（例如 HCN）与没有分⼦偶极⼦的分⼦碰撞时，碰撞本⾝会通过分⼦内电⼦密度的变化导致偶极⼦出现。HCN 中的氮原⼦富含电⼦，分⼦偶极⼦指向该原⼦的⽅向。碰撞时，第⼆个原⼦的电⼦云会被氮上多余的电⼦密度排斥，因此带正电的原⼦核会更
靠近 N 并与其相互作⽤。
Fig 2.诱导⼒⼒⽰意图
3.⾊散⼒（London force）：瞬时偶极之间的电性引⼒；(诱导偶极——诱导偶极相互作⽤)
第三个也是主要的贡献是⾊散或伦敦⼒（波动偶极诱导偶极），由于分⼦中电⼦和原⼦核不停地运动,⾮极性分⼦的电⼦云的分布呈现有涨有落的状态，从⽽使它与原⼦核之间出现瞬时相对位移，产⽣了
瞬时偶极，分⼦也因⽽发⽣变形。分⼦中电⼦数愈多、原⼦数愈多、原⼦半径愈⼤，分⼦愈易变形.瞬时偶极可使其相邻的另⼀⾮极性分⼦产⽣瞬时诱导偶极，且两个瞬时偶极总采取异极相邻状态，这种随时产⽣的分⼦瞬时偶极间的作⽤⼒为⾊散⼒(因其作⽤能表达式与光的⾊散公式相似⽽得名)。瞬间偶极距随时间不断变化，统计平均为0，因此实验测量得到的分⼦偶极距离任然为0，⾊散⼒作为⼀种相互吸引⼒，可以使得体系能量降低，⼀般地，极化的体积分别为和的两个分⼦或离⼦间⾊散能可以写成：
虽然瞬时偶极存在暂短，但异极相邻状态却此起彼伏，不断重复，因此分⼦间始终存在着⾊散⼒.⽆疑,⾊散⼒不仅存在于⾮极性分⼦间，也存在于极性分⼦间以及极性与⾮极性分⼦间。⾊散⼒存在于⼀切分⼦之间.⾊散⼒与分⼦的变形性有关，变形性越强越易被极化,⾊散⼒也越强。稀有⽓体分⼦间并不⽣成化学键，但当它们相互接近时，可以液化并放出能量，就是⾊散⼒存在的证明.
所有分⼦通过伦敦⾊散⼒或诱导偶极相互作⽤。在下图中，⼀个2个原⼦的分⼦与⼀个3个原⼦的分⼦碰撞。第⼀个分⼦的电⼦云排斥它撞击的分⼦的电⼦云，导致⼀些电⼦密度远离原⼦核。原⼦核然后被它⾃⼰的电⼦屏蔽得很差，并吸引了第⼀个分⼦的电⼦云。
Fig 3.⾊散⼒⽰意图
4.由泡利不相容原理产⽣的排斥成分，可防⽌分⼦坍塌。
上⾯讨论的分⼦间相互作⽤均为吸引作⽤，并且势能随分⼦间距离的减⼩⽽降低。但是，这些关系并不适合分⼦相互靠得特别近的情况。当两个分⼦相互靠近，达到或接近电⼦云相互重叠的距离时，分⼦间的排斥作⽤将显著增加，引起体系势能的迅速升⾼。此外，原⼦核间的静电排斥作⽤，也引起体系的势能升⾼。完整的分⼦间相互作⽤函数必须同时考虑排斥作⽤的贡献。理论分析表明, 排斥作⽤对体系势函数的贡献与分⼦间的距离呈指数关系：
为了计算⽅便和效率，常把排斥能表⽰为幂函数的形式：αi αj E =d −23
[I +I i j I I i j ][(4πε)02ααi j r 61](3)
u (r )=Aexp (−βr )
(4)
式中，A为参数;n则在8~16中变化。同时排斥⼒也是⼀会⾛近程相互作⽤。随分⼦间距离增⼤⽽迅速衰减。当原⼦间距离低于0.4 nm 时，
此时的对范德华⼒的贡献主要来⾃泡利不相容原理，让物质内部有所排斥，表现为排斥⼒。
Fig 4. 各分⼦范德华⼒吸引⼒的分配情况
可以说范德华⼒属于库伦相互作⽤（静电相互作⽤），但是我们经常看到分⼦间相互作⽤⼤体分为两种：静电相互作⽤和范德华⼒相互作⽤（⼴义的）。前⾯说到的范德华⼒属于静电相互作⽤，其为狭义的范德华⼒它通常指伦敦⾊散⼒（London dispersion force），它的相互作⽤能与分⼦间距（r）是的关系，⽽⼴义的范德华⼒，是因为Keesom force和Debye force在热平均后都与London force⼀样具
有的关系。如果没有热平均，偶极-偶极相互作⽤是的关系也就是说Keesom force要更加短程，⽽相对短程也是范德华⼒（短程相互作⽤）区别于静电相互作⽤（长程相互作⽤）的特征。因此可以看出⾮键相互作⽤分类依据为两者相互作⽤随距离变化的衰减特点不同。范德华⼒的主要特点是：它们⽐普通的共价键和离⼦键弱。范德华⼒是可加的，不能饱和。它们没有⽅向特性。它们都是短程⼒，因此只需要考虑最近粒⼦（⽽不是所有粒⼦）之间的相互作⽤。分⼦距离越近，范德华引⼒越⼤。除了偶极-偶极相互作⽤之外，范德华⼒与温度⽆关。
在低分⼦量醇中，其极性羟基的氢键性质⽀配了其他较弱的范德华相互作⽤。在较⾼分⼦量的醇中，⾮极性烃链的性质占主导地位并决定了它们的溶解度。2.2.1 Lennard-Jones potential 模型
Lennard-Jones potential，它以John Lennard-Jones.（约翰·伦纳德-琼斯）的名字命名被称为LJ势或者12-6势。它描述了简单原⼦和分⼦之间相互作⽤的基本特征：两个相互作⽤的粒⼦在⾮常近的距离相互排斥，在中等距离相互吸引，在⽆限距离上不相互作⽤，LJ势是⼀对势，即该势不包含三体或多体相互作⽤。LJ势是研究最⼴泛和最彻底的势，被认为是简单⽽现实的分⼦键相互作⽤的半经验模型。LJ势的常⽤表达式为：
u (r )=A /r n
E d r 61
r 61r 31E =LJ 4ε()−()[r ij σ12r ij σ6](5)
：为两体之间的势能⼤⼩。
$\varepsilon $ ：为控制相互作⽤的强度（以 eV 为单位），从本质上讲，它衡量两体相互吸引的强度，通常称为势阱深度（也叫“⾊散能”）;
：为相互作⽤的势能正好为零时的两原⼦间的距离，给出了两个⾮键合粒⼦的接近程度的度量，因此被称为范德⽡尔斯半径。它等于⾮成键粒⼦之间核间距离的⼆分之⼀。
: 为i,j两体之间的距离。
在实际应⽤中，，参数往往通过拟合已知实验数据或精确量⼦计算结果⽽确定。 另⼀种写法是：
：为势阱底部时，两体间的距离。
物理背景有很多⽅程能描述范德华相互作⽤，其中最常见的为LJ势函数。Lennard-Jones 势模拟了两个最重要和最基本的分⼦相互作⽤：其中排斥
项()描述了在相互作⽤粒⼦的短距离上由于重叠电⼦轨道⽽产⽣的泡利斥⼒，吸引项()描述了长距离相
互作⽤（主要为⾊散⼒）的吸引⼒，在两个粒⼦之间的⽆限距离处消失。短距离的陡峭排斥相互作⽤导致固相和液相的低压缩性；有吸引⼒的分散相互作⽤对凝聚相起到稳定作⽤，尤其是汽液平衡。
吸引项的函数形式，即指数“6”，具有物理上的合理性，对于具有指数“12”的排斥项并不严格。简单原⼦和分⼦之间有吸引⼒的⾊散相
互作⽤是部分电荷波动的结果。量⼦化学计算表明，这种⾊散贡献必须以()的速度衰减。
() 项被使⽤，主要是因为它在计算上可以⾮常有效的实现为()的平⽅，这对于’12’以外的值不成⽴。此外，() 合理地接近了泡利斥⼒，可以使⽤任意指数⽽不是12和6来推⼴Lennard-Jones势能，在这⾥我们仅讨论LJ 6-12.
下⾯以单原⼦分⼦氩为例，应⽤LJ势能函数描述氩-氩间的相互作⽤，使⽤以下Python代码定义了LJ势的每个组成部分和相互作⽤的总能量。然后将氩-氩间LJ势能的变化绘制在⼀张图上。和的值是与氩-氩相互作⽤有关的值。
#内嵌绘图
%matplotlib inline
import  numpy as  np
import  matplotlib .pyplot as  plt
#1.LJ 势能函数主要代码：...............................................................
#定义吸引⼒项
def  attractiveEnergy  (r ,epsilon ,sigma ):
return  -4 * epsilon * np .power (sigma /r ,6)
#定义排斥项
def  repulsiveEenergy (r ,epsilon ,sigma ):
return  4 * epsilon * np .power (sigma /r ,12)
#定义LJ 6-12 势能函数
def  ljEenergy (r ,epsilon ,sigma ):
return  attractiveEnergy (r ,epsilon ,sigma ) + repulsiveEenergy (r ,epsilon ,sigma )
#.....................................................................................
#2.图像绘制....................................................
#画氩-氩原⼦间的LJ 相互作⽤图像（参数：a=epsilon; b=sigma ）
a =0.0103
b =3.4
#预设氩-氩原⼦为3到8埃，产⽣以下100个系列数据E LJ σσr ij εσE (r )=LJ ij ϵ()−2()（amber ⽤的这个）
[r j i r min
12r ij r min 6](6)
r =min 2σ=61  1.122σr ij 121r ij 61r ij 61r ij 121r ij 61r ij 121ϵσ

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。