1.分子动力学方法

1.1 简单介绍

分子动力学模拟是一种在计算机上进行的虚拟实验:在分子动力学模拟中,原子被视为质点,质点的位置和速度是分子动力学处理的变量.在给定原子初始坐标和初始速度的前提下,通过求解运动方程得到原子随时间的运动轨迹,在不同的宏观热力学控制条件下,分子动力学求解不同的运动方程获得对相应的系综的采样.例如微正则系综需要求解 Hamilton 系统,正则系综需要求解 Nose-Hoover 动力学或 Langevin 动力学,假设我们可以求解无穷长的原子运动轨道,且在遍历性假设成立的条件下,这些轨道经历的体系状态对应相应系综分布的一个采样.体系的大量可观测性质可以写为相应系综下的期望或时间关联函数,因此可以通过分子动力学模拟进行计算.

1.2 verlet算法

首先列出分子满足的哈密顿正则方程：

$$\frac{d{q}_i}{dt}=\frac{\partial H}{\partial p_i},\frac{d{p}_i}{dt}=-\frac{\partial H}{\partial q_i}$$

在非相对论情况下做泰勒展开，我们得到：

$$r(t+\Delta t)=2r(t)-r(t-\Delta t)+\frac{f(t)}{m}\Delta t^2+O(\Delta t^4)$$

这里的$r$可以精确到$\Delta t^4$量级。然而，随着模拟的推进，$r$的算法由于引进了之前累积下来的误差，其误差会越来越大。可以证明，如果模拟的时长$T=n\Delta t$，那么$r(t+T)$的累积误差可以达到$\Delta t^2$量级。

同样对于速度的递推关系式有：

$$v(t)=\frac{r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}+O(\Delta t^2)$$

由于速度误差不累计，$v$可以精确到$\Delta t^2$量级。

以上两式合称为Verlet算法。

通过上式我们可以得到系统的速度等微观物理量，再经由一定的热力学关系得到宏观物理量。如温度可以通过$E_k=\frac{3}{2}(N-1)k_{B}T$来计算。

除此之外，进行MD模拟时也要考虑系综的选择，其主要作用是在模拟过程中进行系综调节。系综调节主要是指在进行分子动力学计算过程中，对温度和压力参数的调节，包括调温技术和调压技术。

1.3 势能面与深度学习理论

在上述理论中，还缺失了一个关键的参数$f(t)$，这个参数体现了不同分子间的本质区别，即分子间的相互作用。这可以通过势能面来描述.

势能面对应一个系统在给定参数下的势能，这个系统通常是原子的集合，而这些参数通常是原子的坐标.

描述势能面的方法有多种。一是给定一个经验的解析的表达式,计算效率高，但结果不一定和实验吻合.二是用第一性原理,即密度泛函理论来计算,这需要求解薛定谔方程,增加了运算量的同时提高了准确度。

除此之外较为先进的方法是机器学习方法。先用第一性原理精确建模,再用深度学习方法对高维的势函数进行拟合，以此实现计算速度和精度的优化平衡。

关于深度学习构造势函数的原理，参考以下文献说明王涵.机器学习原子间相互作用建模[J].计算数学,2021,43(3):18.DOI:10.12286/jssx.j2021-0833.。

假设体系的原子坐标为$\mathcal{R}=\{r_1,r_2,\dots ,r_N\}$，通过第一性原理精确计算得到体系能量$E=E(\mathcal{R})$和第i个原子的受力$F_i$，以及体系的virial张量${\Xi }_{\alpha \beta }=-\frac{\partial E}{\partial h_{\gamma \alpha }}{h}_{\gamma \beta }$。利用这三个量得到损失函数 $$\begin{array}{r}\mathcal{L}=\frac{1}{|\mathcal{B}|N^2}\left\{\sum _{k\in \mathcal{B}}\right(p_e\Vert E({\mathcal{R}}^{(k)})-E^w({\mathcal{R}}^{(k)}){\Vert }^2+p_f\Vert \mathbf{F}({\mathcal{R}}^{(k)})-{\mathbf{F}}^w({\mathcal{R}}^{(k)}){\Vert }^2+p_v\Vert \mathbf{\Xi }({\mathcal{R}}^{(k)})-{\mathbf{\Xi }}^w({\mathcal{R}}^{(k)}){\Vert }^2\left)\right\}\end{array}$$ 其中$\mathcal{B}$为批次大小，三项前各自有个系数$p$调整对损失函数的贡献。我们只要最小化损失函数就完成了训练。

在物理上考虑，能量应该是广延的，即可以写成对每个原子的分解式。 $$E^w=\sum _{i=1}^N{E}_i^w=\sum _{i=1}^N{E}_i^w({\mathcal{R}}_i)$$ 从而我们只需要对每个$E_i$分别建模。其次，体系的能量有对称性要求，包括平移，旋转，交换对称性，这体现为，能量有约束 $$E(\mathcal{R})=E(U\mathcal{R})$$ U代表某种变换。但一般的神经网络是不会满足这样的对称性的，所以我们可以引入描述子： $$\mathcal{D}({\mathcal{R}}_i)=\mathcal{D}(U{\mathcal{R}}_i),\forall U\in \mathcal{U}$$ 其为对原子坐标的映射，满足所有的对称性要求，再将建模的能量函数写成关于描述子的函数，这样能量就自然的满足了对称性。

2.DPMD方法计算举例

以下将举例利用DPMD+LAMMPS方法，采取NPT系综，模拟计算Sc元素的熔点。（所需文件见数据集）

判断熔化相变发生的主要依据为金属密度的变化，查阅资料得到Sc密度在固态常温下约为2.985g/cm^3。

首先准备Sc-molten-model.pb文件作为训练出的势能模型，然后构建INPUT文件Sc.in以及结构文件Sc.data如数据集所示。

由于Sc的熔点在1800K左右，首先选取t的范围1500-2000K，步长100K进行初步扫描。设定时间步长为1fs，总运行5000步。

根据图像确定，相变大概发生在1700K左右。进一步细化扫描得到：

关于步长的设定，需要考虑时间成本和准确性得到一个折中的大小，查阅资料可得，取分子振动周期的1/10比较合适。