分子动力学模拟(MD)

分子动力学模拟是一种通过模拟体系原子和分子运动轨迹研究体系微观和宏观性质的计算机模拟方法，是研究分子体系的重要手段。根据经典力学，只要给定分子体系中原子和分子之间的相互作用，就可以利用牛顿运动方程可以求解分子或原子的运动轨迹，然后使用一定的统计方法计算出系统的力学、热力学、动力学性质。在分子动力学中首先将由N个粒子构成的系统抽象成N个相互作用的质点，每个质点具有坐 标(通常在笛卡尔坐标系中)、质量、电荷及成键方式，按目标温度根据Boltzmalm分布随机指定各质点的初始速度，然后根据所选用的力场中的相应的成键和非键能量表达形式对质点间的相互作用能以及每个质点所受的力进行计算。接着依据牛顿力学计算出各质点的加速度及速度，从而得到经一指定积分步长(通常1fs)后各质点新的坐标和速度，这样质点就移动了。经一定的积分步数后，质点就有了运动轨迹。设定时间间隔对轨迹进行保存。最后可以对轨迹进行各种结构、能量、热力学、动力学、力学等的分析，从而得到感兴趣的计算结果。分子动力学模拟方法已被广泛应用于化学、物理、材料科学、生物学等学科领域。

MD模拟的统计力学基础

统计力学利用统计学方法和微观粒子的运动规律，研究多粒子体系的宏观性质，是现代物理学的支柱之一。利用统计力学解决热力学问题时，统计力学被称为统计热力学。宏观体系处于宏观静止状态，但其组成原子、分子处在不断的运动状态。在热力学，常用一组宏观可观测量描述宏观体系的状态。Gibbs提出，可以利用宏观可观察量作约束条件，定义具有不同统计特性的热力学体系，即热力学系综。重要的热力学系综包括：微正则系综（microcanonical ensemble, NVE），正则系综（canonical ensemble, NVT），巨正则系综（grand canonical ensemble, µVT），等温等压系综（isothermal–isobaric ensemble, NPT），等焓等压系综（isoenthalpic-isobaric ensemble, NPH）。

在NVE系综中，体系的粒子数（N）、体积（V）、总能量（E）均保持守恒，对应与环境没有能量交换、粒子交换的孤立体系。在NVT系综中，体系的粒子数（N）、体积（V）、温度（T）保持守恒，在MD模拟中通过体系与恒温热浴的能量交换实现温度守恒。常用控温算法有Nosé–Hoover控温、Berendsen控温、Andersen控温、Langevin控温等。在NPT系综中，体系的粒子数（N）、压强（p）、温度（T）均保持守恒，但体系体积可发生变化。实验室化学反应常在恒压下进行，NPT系综常被用于与化学反应相关的研究。µVT系综对应开放体系，不但与热浴存在能量交换，还与环境进行物质交换，体系的能量和粒子数不守恒，但系综与环境保持统计平衡。在NPH系综中，保持系统的粒子数N、压力P和焓值H都不变。然而，模拟时要保持压力与焓值为固定值,有一定难度。事实上,这种系综在实际的分子动力学模拟中很少见。

根据统计热力学，可以利用经典力学或量子力学计算系综平均，得到体系的热力学性质等。

分子内和分子间相互作用

分子间力是分子间相互作用产生的力，包括分子与近邻粒子(原子或离子)之间相互作用产生的引力或斥力，影响体系的物理化学性质。分子内力是同一分子中的原子通过共用电子对产生相互作用力，影响体系的化学性质。分子间力通常小于分子内力，共用电子引起的共价键远强于分子间相互作用力。分子体系相互作用由分子间力和分子内力组成。

在MD模拟中，常用分子间相互作用势U(r)来表示分子间相互作用，r为分子间距离。因此，分子间相互作用势决定了物质的结构性质和物理化学性质。Van der Waals在研究实际气体和液体性质时，考虑了分子间相互作用力，提出在分子间距离r较小时分子间存在斥力，随着距离r的增大分子间的斥力逐渐被引力所取代，推导得到实际气体van der Waals状态方程。

分子间相互作用力分为分子间排斥力和吸引力两个部分。分子间斥力的本质是遵循Pauli不相容原理的波函数发生重叠时的排斥力。分子间引力由分子电子分布中心与原子核正电中心分离、色散等效应引起。在Lennard-Jones势函数中，近似于一对非键合原子或分子之间的非静电相互作用的势能，分子间引力随 $r^{-6}$ 变化，分子间斥力随$r^{-12}$ 变化，数学函数可表示为

$U(r)=4{\phi }_0[(\frac{\sigma }{r}{)}^{12}-(\frac{\sigma }{r}{)}^6]$

其中，r=σ对应体系势能为0的位置，${\phi }_0$是势阱的深度。对于复杂体系，难于利用量子力学理论精确计算U(r)，也难于用实验测量分子间相互作用力。

MD模拟的积分算法

实际体系中，描述粒子间相互作用的势能函数复杂，难以求解牛顿运动方程的解析解。因此，在MD模拟过程中采用数值积分的方法求解体系的运动方程，有关常用数值积分算法如下，

1. Verlet算法：利用Taylor展开粒子的位置坐标， $$\begin{array}{rl}& r(t+\delta t)=r(t)+v(t)\delta t+\frac{1}{2}a(t)\delta t^2\\ & r(t-\delta t)=r(t)-v(t)\delta t+\frac{1}{2}a(t)\delta t^2\end{array}$$ 上述两式相加得到Verlet算法的基本公式，

$$\begin{array}{rl}& r(t+\delta t)=2r(t)-r(t-\delta t)+a(t)\delta t^2\end{array}$$

利用Verlet算法计算$t+\delta t$时刻的粒子位置，需要t时刻的粒子位置和加速度、以及 $t-\delta t$时刻的粒子位置。该算法计算简单，不直接计算粒子的速度，但算法精度不高。

2.Leap-frog算法：首先计算$t+\frac{1}{2}\delta t$时刻的粒子速度，

$$\begin{array}{r}v(t+\frac{1}{2}\delta t)=v(t-\frac{1}{2}\delta t)+a(t)\delta t\end{array}$$

然后，计算 $t+\delta t$时刻的粒子位置，

$r(t+\delta t)=r(t)+v(t+\frac{1}{2}\delta t)\delta t$

Leap-frog算法方法虽然直接计算体系的粒子速度，但体系的粒子速度和位置不同步。 t时刻的粒子速度近似为，

$v(t)=\frac{1}{2}[v(t-\frac{1}{2}\delta t)+v(t+\frac{1}{2}\delta t)]$

3.Velocity Verlet算法：具有更高的计算精度，体系的粒子位置、速度分别表示为，

$$\begin{array}{rl}& r(t+\delta t)=r(t)+v(t)\delta t+\frac{1}{2}a(t)\delta t^2\\ & v(t+\delta t)=v(t)+\frac{1}{2}[a(t)+a(t+\delta t)]\delta t\end{array}$$

4.Beeman’s算法：基于Verlet算法改进体系的粒子位置和速度分别为， $$\begin{array}{rl}& r(t+\delta t)=r(t)+v(t)\delta t+\frac{2}{3}a(t)\delta t^2-\frac{1}{6}a(t-\delta t)\delta t^2\\ & v(t+\delta t)=v(t)+v(t)\delta t+\frac{1}{3}a(t)\delta t+\frac{5}{6}a(t)\delta t-\frac{1}{6}a(t-\delta t)\delta t\end{array}$$ Beeman’s算法的计算精确度得到了极大的提高，但计算成本也相应提高

分子力场类型

分子力场根据量子力学的波恩-奥本海默近似，一个分子的能量可以近似看作构成分子的各个原子的空间坐标的函数，简单地讲就是分子的能量随分子构型的变化而变化，而描述这种分子能量和分子结构之间关系的就是分子力场函数。一般而言，分子力场函数由以下几个部分构成：

键伸缩能：构成分子的各个化学键在键轴方向上的伸缩运动所引起的能量变化

键角弯曲能：键角变化引起的分子能量变化

二面角扭曲能：单键旋转引起分子骨架扭曲所产生的能量变化

非键相互作用：包括范德华力、静电相互作用等与能量有关的非键相互作用

交叉能量项：上述作用之间耦合引起的能量变化

不同的分子力场会选取不同的函数形式来描述上述能量与体系构型之间的关系。

实例：

GAFF

GAFF(Generation Amber Force Field) 力场是普适型有机小分子力场，函数形式和AMBER力场相同，与AMBER力场完全兼容。 $$\begin{array}{rl}{E}_{pair}=& \sum _{bonds}{K}_r(r-r_{eq}{)}^2+\sum _{angles}{K}_{\theta }(\theta -{\theta }_{eq}{)}^2+\sum _{dihedrals}\frac{V_n}{2}[1+\cos (n\varphi -\gamma )]+\sum _{i<j}[\frac{A_{ij}}{R_{ij}^{12}}\\ & -\frac{B_{ij}}{R_{ij}^6}+\frac{q_i{q}_j}{\epsilon R_{ij}}]\end{array}$$

AMBER

AMBER 力场是传统力场之一，主要适用于蛋白质和核酸体系、多糖。 $$\begin{array}{rl}{E}_{pot}=& \sum _b{K}_2(b-b_0{)}^2+\sum _{\theta }{H}_{\theta }(\theta -{\theta }_0{)}^2+\sum _{\varphi }\frac{V_n}{2}[1+\cos (n\varphi -{\varphi }_0)]\\ & +\sum \epsilon [(\frac{r^{\ast }}{r}{)}^{12}-2(\frac{r^{\ast }}{r}{)}^6]+\sum \frac{q_i{q}_j}{{\epsilon }_{ij}{r}_{ij}}+\sum [\frac{C_{ij}}{r_{ij}^{12}}]+\sum [\frac{C_{ij}}{r_{ij}^{12}}-\frac{D_{ij}}{r_{ij}^{10}}]\end{array}$$

MMFF

MMFF(Merck Molecular Force Field)是小分子力场。 $$\begin{array}{rl}{V}_{total}=& \sum _{bonds}{K}_{bond}(r-r_{eq}{)}^2(1+cs(r-req)+\frac{2}{7}(c{s}^2(r-r_{eq}{)}^2))\\ & +\sum _{angle}{K}_{\theta }(\theta -{\theta }_{eq}{)}^2(1+cb(\theta -{\theta }_{eq}))+\sum _{angle,linear}{K}_{al}(1+\cos (\theta ))\\ & +\sum _{stretch,bend}(K_{ijk}(r_{ij}-r_{eq})+K_{kji}(r_{kj}-r_{eq}))(\theta -{\theta }_{eq})+\sum _{outofplane}{K}_{OOP}(\chi {)}^2\\ & +\sum _{dihedrals}\frac{V_1}{2}[1+\cos (\varphi )]+\frac{V_2}{2}[1+\cos (2\varphi )]+\frac{V_3}{2}[1+\cos (3\varphi )]\\ & +\sum _{i<j}[ϵ(\frac{1.07\sigma }{r_{ij}+0.07\sigma }{)}^7(\frac{1.12{\sigma }^7}{r_{ij}^7+0.07{\sigma }^7}-2)-\frac{q_i{q}_j}{D(r_{ij}+\delta )}]\end{array}$$

UFF

UFF(Universall Force Field) 力场覆盖了周期表中所有元素，应用最为广泛。

Bond:

1.Harmonic

$E_R=1/2K_{ij}(r-r_{ij}{)}^2$

2.Morse

$E_R=D_{ij}[e^{-\alpha (r-r_{ij})}-1{]}^2\alpha =[\frac{k_{ij}}{2D_{ij}}{]}^{1/2}$

Angle:

$E_{\theta }=\frac{K_{ijk}}{n^2}[1-\cos (n\theta )]$

Torsion:

$E_{\varphi }=1/2V_{\varphi }[1-\cos n{\varphi }_0\cos n\varphi ]$

LJ:

$E_{vdw}=D_{ij}\{-2[\frac{{\chi }_{ij}}{\chi }{]}^6+[\frac{{\chi }_{ij}}{\chi }{]}^{12}\}$

DPMD

深度势分子动力学（Deep Potential Molecular Dynamics，DPMD）是一种基于深度学习技术的分子动力学模拟方法，它可以高效地模拟分子间的复杂相互作用。相较传统的分子动力学模拟方法，DPMD不需要手动构建势函数，而是使用深度神经网络来学习分子间的相互作用势能函数。这样可以大大提高模拟的准确性和效率。DPMD的核心思想是将分子势能函数表示为一个多层的神经网络模型，输入是分子的结构信息，输出是势能函数。通过在大量的分子结构数据上进行训练，神经网络可以学习到分子间的相互作用势能函数。在模拟过程中，DPMD使用这个学习得到的势能函数来计算分子间的相互作用，从而推动分子动力学模拟的进行。DPMD的优点包括：

可以高效地模拟分子间的复杂相互作用。

不需要手动构建势函数，节省了大量的人力和时间成本。

模拟结果更加准确，可以更好地反映实际分子系统的行为。

可以应用于复杂的分子体系，如蛋白质、DNA等。

通过深度学习技术的不断发展，进一步提高模拟的准确性和效率。

运行实例（基于LAMMPS）

LAMMPS

LAMMPS程序（Large-sale Atomic/Molecular Massively Parallel Simulator）是由美国Sandia国家实验室开发的开源经典力学MD模拟程序，侧重于材料领域的模拟研究。LAMMPS程序在模拟固态材料（金属、半导体）、柔性物质（生物分子、聚合物）、 粗粒度介观体系等方向具有广泛的应用。持续内置多种原子间势（力场模型），可以实现原子、聚合物、生物分子、固态材料（金属、陶瓷、氧化物）、粗粒度体系的建模和模拟。该程序既可以模拟二维体系，也可以模拟三维体系，可以模拟多达数百万甚至数 十亿粒子的分子体系，具有模拟效率高、计算时间短等优点。LAMMPS程序具有良好的用户界面，用户可以自由修改或扩展新的力场模型、原子类型、边界条件等以满足课题研究需求。LAMMPS可以在单个处理器的台式机和笔记本本上运行且有较高的计算效率， 但是它是专门为并行计算机设计的。他可以在任何一个安装了C++编译器和MPI的平台上运算，这其中当然包括分布式和共享式并行机和Beowulf型的集群机。LAMMPS的部分功能还支持OpenMP多线程、矢量化和GPU加速。通常意义上来讲，LAMMPS是根据不同的边界条件和初始条件对通过短程和长程力相互作用的分子，原子和宏观粒子集合对它们的牛顿运动方程进行积分。

下面以计算Cu熔点为例进行演示

LAMMPS输入文件

Cu.in

对其中的部分参数进行解释：

variable t equal 1300：是对体系温度的设置，由于Cu的熔点在$1357\text{K}$，故选取的t取值为$1200$~$1400\text{K}$，间隔$20\text{K}$作为粗略观测。

velocity all create $t 23456789使原子运动速度随机分布。

timestep 0.001：设置时间步长为$1,\text{fs}$

run 5000：运行总步数，为保证$5 10,\text{ps}$的模拟时长，则总步数至少为$5000$

Cu.data

盒子长度选为两倍晶格长度，原子数量取为32，在保证计算精度的同时使得$Density=\frac{Volume}{Totalmass}$接近铜的实际密度8.9$g/c{m}^3$

输入运行脚本：

绘制密度关于温度变化的图像如下：

密度在1300K左右发生突变，故熔点在1300K附近,与Cu实际熔点接近