计算物理复习[ep.4|分子动力学、第一性原理与密度泛函理论、机器学习基础]

本文含有大量数学公式，请使用电脑或平板电脑浏览以获得最佳阅读体验

第八章 分子动力学

短程相互作用力的计算方法

二体势与三体势

二体势、三体势（多体势）$U_{2}(\mathbf{r}_{i},\mathbf{r}_{j})\quad U_{3}(\mathbf{r}_{i},\mathbf{r}_{j},\mathbf{r}_{k})$ 来自粒子之间的相互作用，决定了整个系统的性质，是分子动力学准确性的先决条件。
实践中通常把这些作用势表达为一个解析函数，函数形式受物理原理启发 (比如中远程吸引、近程排斥)，其参数则通过拟合实验结果和第一性原理计算结果来得到，因此称作经验作用势 (empirical potential) 或者直接简称为势 (potential)。

常见二体势

• 兰纳-琼斯势 Lennard-Jones Potential

  用来模拟两个电中性的分子或原子间相互作用势能的一个比较简单的数学模型，定义两粒子间的距离为 r：

  $$\boldsymbol{r}=\|\mathbf{r}_{i}-\mathbf{r}_{j}\|_{2} \tag{1}$$

  则LJ势能表达式为：

  $$U_{\text{LJ}}(\mathbf{r}_i,\mathbf{r}_j)=4\epsilon\bigg[\bigg(\frac{\sigma}{r}\bigg)^{12}-\bigg(\frac{\sigma}{r}\bigg)^6\bigg] \tag{2}$$

  作用力表达式为：

  $$f_{ij}=-\frac{\mathrm{d}U_{\mathrm{LJ}}}{\mathrm{d}r}=4\epsilon\bigg(12\frac{\sigma^{12}}{r^{13}}-6\frac{\sigma^6}{r^7}\bigg) \tag{3}$$

  

• 纯排斥势 Weeks-Chandler-Andersen Potential

  $$V(r)=\begin{cases}{4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right]}&{r\leq2^{1/6}\sigma}\\ {0}&{r>2^{1/6}\sigma}\end{cases} \tag{4}$$

• Morse 势

  Morse 势是一种对于双原子分子间势能的简易解析模型。 一方面，对 Morse 势求解薛定谔方程具有解析解，方便分析问题；另一方面，由于它隐含地包括了键断裂这种现象，对于分子振动的微细结构的具有良好的近似。

  $$U_{\text{Morse}}(r)=-D_e+D_e\Big(1-e^{-a(r-r_e)}\Big)^2 \tag{5}$$

  

• 模拟共价键连接：弹簧势能 Harmonic Spring Potential

  $$E_{\mathrm{bond}}=K_{\mathrm{bond}}\bigl(r-r_0\bigr)^2 \tag{6}$$

热浴

两粒子足够接近时，会在短时内受到巨大的力，使速度巨量增加，温度无法短时弛豫，造成系统崩溃。

• 简单标度：

  $$\begin{aligned} & \frac{1}{2}\langle m v_{\mathrm{cur}}^{2}\rangle=\frac{3}{2}k_{B}T_{\mathrm{cur}} \\ & \frac{1}{2}\langle m v_{\mathrm{tar}}^{2}\rangle=\frac{3}{2}k_{B}T_{\mathrm{tar}} \end{aligned} \quad \Rightarrow \quad v_\mathrm{tar}=v_\mathrm{cur}\sqrt{\frac{T_\mathrm{tar}}{T_\mathrm{cur}}} \tag{7}$$

• Berendsen 标度：

  $$v_{\mathrm{tar}}=v_{\mathrm{cur}}\sqrt{1+\frac{\Delta t}{\tau}\left(\frac{T_{\mathrm{tar}}}{T_{\mathrm{cur}}}-1\right)} \tag{8}$$

• Andersen 热浴：

  1. 从一组位置和动量开始计算
  2. 整合运动方程
  3. 选择粒子进行与热浴的碰撞；选择的概率为 $\nu\Delta t$
  4. 选定用于碰撞的粒子的速度取自所需温度下的玻尔兹曼分布

• Nosé-Hoover 热浴

  • 修改运动方程，使分子动力学 Molecular Dynamics(MD) 系统能够与恒温的热储器相互作用

    $$\left\{\begin{array}{l} \frac{d\mathbf{p}_i}{dt}=\mathbf{F}_{i}\\ \frac{d\mathbf{r}_{i}}{d t}=\frac{\mathbf{p}_{i}}{m_{i}} \end{array}\right . \quad \Rightarrow \quad \begin{cases} \frac{d\mathbf{p}_i}{dt}=\mathbf{F}_{i}-\eta\mathbf{p}_{i}\\ \frac{d\mathbf{r}_{i}}{d t}=\frac{\mathbf{p}_{i}}{m_{i}} \end{cases} \tag{9}$$

    其中 $\eta$ 是阻尼系数

  • $\eta$ 的变化根据是

    $$\frac{d\eta}{dt}=\frac{1}{\tau_T^{2}}\bigg(\frac{T}{T_{set}}-1\bigg) \tag{10}$$

    • $T_{\textit{set}}$ 预期温度

    • $\tau_T$: 时间常数，描述了系统和热储层的耦合强度。时间常数越长意味着 MD 系统对储层的反应越慢。

周期性边界条件

• 计算位置：

  if periodic_x
      if x < -x_size*0.5
          x = x + x_size;
      end
      if x >= x_size*0.5
          x = x - x_size;
      end
  end

• 计算距离：

  if periodic_x
  	dx = x(j) - x(i);
      if dx > x_size*0.5
          dx = dx - x_size;
      end
      if dx <= -x_size*0.5
          dx = dx + x_size;
      end
  end

Q：请问图中红圈标出的两个粒子的距离是否是红线的长度？

A：不是，考虑周期性边界条件，真实距离如绿线所示。

分子动力学模拟过程

关于分子动力学的模拟练习，可以参考这篇文章，该文章用使用 GROMACS 软件进行 Lennard-Jones 势能下的液氩模拟。

第九章 第一性原理与密度泛函理论

第一性原理（ab-initio）

在物理学中，第一性原理，或称从头算，指从基本的物理学定律出发，不外加假设与经验拟合的推导与计算。例如利用薛定谔方程在一些近似方法下解电子结构，但不从实验数据得到拟合参数的从头计算法。

广义从头计算 – 密度泛函理论

用电子密度取代波函数做为研究的基本量

Hohenberg – Kohn 定理

• 定理一：体系的基态能量仅仅是电子密度的泛函。

  $$E_0=E_0[\{\rho_0\}] \tag{11}$$

• 定理二：以基态密度为变量，将体系能量最小化后就得到了基态能量。

根据这两个基本定理，系统的基态能量可以分为多电子系统的相互作用能 (U) 、多电子在外势场中的能量 (V ) 和动能 (T )，即：

$$E=T+V+U \tag{12}$$

该理论对任意相互作用的多粒子系统有效。

Kohn – Sham方程

该方程在密度泛函理论里面指的是与真实体系相关的虚拟体系所满足的薛定谔方程。该虚拟体系中的粒子（通常是电子）在无相互作用的有效势场中运动。虚拟系统中的粒子是彼此无相互作用的费米子，因此 Kohn – Sham 方程的精确解为单个斯莱特行列式，行列式中的轨道则称为 Kohn – Sham 轨道，每一个 Kohn – Sham 轨道都可以表示为原子轨道的线性组合，也可以按照基函数展开。方程形式为：

$$\bigg(-\frac{\hbar^2}{2m}\nabla^2+v_{\textrm{eff}}\left(\mathbf{r}\right)\bigg)\chi_i(\mathbf{r})=\varepsilon_i\chi_i(\mathbf{r}) \tag{13}$$

式中 $\varepsilon_i$ 为 $K-S$ 轨道 $\chi_i$ 的轨道能。含有 $N$ 个粒子的 $K-S$ 系统的电子密度则由下式给出：

$$\rho(\mathbf{r})=\sum_i^N|\chi_i(\mathbf{r})|^2 \tag{14}$$

用无相互作用的粒子模型代替有相互作用粒子的能量：

$$E[\rho]=T_0[\rho]+E_{\mathrm{ext}}[\rho]+E_{\mathrm{H}}[\rho]+E_{\mathrm{XC}}[\rho] \tag{15}$$

其中：

• $E_{\mathrm{ext}}[\rho]$：电子在外势场中的能量
• $E_{\mathrm{H}}[\rho]$：Hartree 项，描述电子之间由于库仑力产生的势能
• $E_{\mathrm{XC}}[\rho]$：电子之间的交换关联能，用来代替电子之间的相互作用

交换关联能近似

• 局域密度近似（Local Density Approximation, LDA）

  近似认为 只与 相关但与电子密度梯度无关。

  $$E_{\mathrm{XC}}^{\mathrm{LDA}}[\rho]=\int\rho(r)\varepsilon_{\mathrm{XC}}(\rho)\mathrm{d}r \tag{16}$$

  $\rho$ 为电子密度， $\varepsilon_{\mathrm{XC}}$ 为交换相关能量密度，它仅仅是电子密度的函数。

• 广义梯度近似（Generalized Gradient Approximation, GGA）

  交换关联能密度不仅与该体积微元内局域电荷密度有关,还与邻近体积微元内电荷密度有关,这种情况下就要考虑整个空间电荷密度的变化,也就是电荷密度梯度。

  $$E_{\mathrm{XC}}^{\mathrm{GGA}}[\rho]=\int\varepsilon_{\mathrm{XC}}[\rho_{\alpha}(r),\rho_{\beta}(r),\nabla\rho_{\alpha}(r),\nabla\rho_{\beta}(r)]\mathrm{d}r \tag{17}$$

第十章 机器学习基础

监督学习与非监督学习

• 监督学习 (Supervised Learning)

  给定数据 $(X_1,Y_1), (X_2,Y_2), … ,(X_n,Y_n)$

  对新的 $X_i$，预测其 $Y_i$

  分类，回归

• 非监督学习 (Unsupervised Learning)

  给定数据 $X_1, X_2, … ,X_n$

  $f(X_i), f(X_i, X_j)$

  降维，聚类

监督学习和非监督学习的核心差异在于数据是否存在标签

监督学习经典模型与算法

1. 线性回归
2. 决策树
3. 逻辑回归
4. 朴素贝叶斯
5. 支持向量机
6. 神经网络

分类问题评价指标（二分类）

• 平均绝对误差 MAE(Mean Absolute Error)

  $$MAE=\frac{1}{n}\sum_{i=1}^n|f_i-y_i| \tag{18}$$

• 均方误差 MSE(Mean Square Error)

  $$MSE=\frac{1}{n}\sum_{i=1}^n(f_i-y_i)^2 \tag{19}$$

• 方根误差 RMSE(Root Mean Square Error)

  $$RMSE=\sqrt{MSE} \tag{20}$$

• R 平方

  $$r^2=1-\frac{SS_{res}}{SS_{tot}}=1-\frac{\sum(y_i-f_i)^2}{\sum(y_i-\bar{y})^2} \tag{21}$$

• 分类问题-模型评估参数：混淆矩阵

  

  混淆矩阵各个部分含义

  True Positive（TP）：真正类。样本的真实类别是正类，并且模型识别的结果也是正类。

  False Negative（FN）：假负类。样本的真实类别是正类，但是模型将其识别为负类。

  False Positive（FP）：假正类。样本的真实类别是负类，但是模型将其识别为正类。

  True Negative（TN）：真负类。样本的真实类别是负类，并且模型将其识别为负类。

  • 精确率（Accuracy）：精确率是最常用的分类性能指标。可以用来表示模型的精度，即模型识别正确的个数/样本的总个数。

    $$Accuracy = \frac{TP+TN}{TP+FN+FP+TN} \tag{22}$$

  • 精度（Precision）：表示在模型识别为正类的样本中，真正为正类的样本所占的比例。

    $$Precision = \frac{TP}{TP+FP} \tag{22}$$

  • 召回率（Recall）：召回率表现出在实际正样本中，分类器能预测出多少。

    $$Recall = \frac{TP}{TP+FN} \tag{23}$$

  • F-β Score：使模型更倾向于 Precision 或 Recall。F-β 的物理意义就是将精度和召回率的一种加权平均，在合并的过程中，召回率的权重是精度的 β 倍。

    $$F_{\beta} =(1+\beta^2)\cdot\frac{\text{Precision}\cdot\text{Recall}}{(\beta^2\cdot\text{Precision})+\text{Recall}} \tag{24}$$

  • F1 Score：它被定义为精度和召回率的调和平均数。在 β=1 的情况，F1-Score的值是从0到1的，1是最好，0是最差。

    $$\operatorname{F1}\operatorname{Score}=2\cdot\frac{\operatorname{Precision\cdot Recall}}{\operatorname{Precision+Recall}} \tag{25}$$

逻辑回归原理、最大似然估计与交叉熵

逻辑回归（线性边界函数）

逻辑回归是一种以决策边界函数为依据进行分类概率计算的算法，以二分类为例，设有数据集 $\{(x_1^1,x_2^1),(x_1^2,x_2^2),\ldots,(x_1^n,x_2^n)\}$ 并存在一条直线 $h(\mathbf{x})=\mathbf{w}^T\mathbf{x}+b=w_1x_1+w_2x_2+b$ 可以将数据分为正类和负类。

决策边界可以表示为 $\mathbf{w}^T\mathbf{x}+b=w_1x_1+w_2x_2+b=0$，如果有某个点 $\mathbf{x^*}(x^*_1,x^*_2)$ 代入边界函数：$w_1x^*_1+w_2x^*_2+b>0$ 那就可以把它判断为正类。

找到判断依据并没有结束，完整的逻辑回归还需要找出判断为正类的概率和输入向量 $\mathbf{x^*}$ 之间的关系。为做判断，引入 Logistic 分布函数：

$$f(x)=\frac{ke^{-k(x-x_0)}}{(1+e^{-k(x-x_0)})^2} \tag{26}$$

它的累积函数（分布函数）为：

$$F(x)=\frac{H}{1+e^{-k(x-x_0)}} \tag{27}$$

概率密度函数和分布函数长这个样子👇：

如果令 $H=1,k=1,x_0=0$ 就可以得到 sigmoid 函数：

$$\sigma(x)=\frac{1}{1+e^{-x}} \tag{28}$$

显然，sigmoid 函数值域为 $[0,1]$ ，$x$ 越大越接近1。把 $h(\mathbf{x})$ 代入得到：

$$y(\mathbf{x}) = \sigma(h(\mathbf{x}))=\frac{1}{1+e^{-(\mathbf{w}^T\mathbf{x}+b)}} \tag{29}$$

求反函数：

$$ln\frac{y}{1-y}=\mathbf{w}^T\mathbf{x}+b \tag{30}$$

这里将 $y$ 视为判断 $\mathbf{x}$ 为正例的概率，则 $1-y$ 为判断 $\mathbf{x}$ 为反例的概率，两者比值称为**几率（odds）**指该事件发生与不发生的概率比值。以 $y=1$ 为判断标准有：

$$\begin{gathered} \mathbf{w}^T\mathbf{x}+b=l n\frac{P(y=1|\mathbf{x})}{1-P(y=1|\mathbf{x})}\\ P(y=1|\mathbf{x})={\frac{1}{1+e^{-(\mathbf{w}^T\mathbf{x}+b)}}} \end{gathered} \tag{31}$$

由（31）式可以看出，通过 sigmoid 函数在输出 $y=1$ 和决策边界函数 $h(\mathbf{x})$ 之间建立起了函数关系，因此根据一个特定的决策边界函数在输入一个向量 $\mathbf{x}$ 后就可以得到一个判断 $\mathbf{x}$ 为正类的概率，然后我们就能根据这个概率来进行判断。

最大似然估计

上一小节告诉我们，逻辑回归是一种以决策边界函数为依据进行分类的算法，它只负责分类。那么有人问了，计算分类概率的方法我会了，接下来如何得到分类的依据即决策边界函数呢？接下来要介绍的最大似然估计就可以解决这个问题。

最大似然估计一言以蔽之：通过观测值求已知模型中的未知参数。对于逻辑回归模型来说，未知的参数就藏在决策边界函数 $h(\mathbf{x})$ 里面，即 $\vec{\theta}=\begin{pmatrix}\mathbf{w}\\ b\end{pmatrix}$。

现有一组观测值 $\{\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n\}$ （可以通过做实验之类的得到），设：

$$\begin{array}{l} P(y=1|\mathbf{x},\vec{\theta})={\frac{1}{1+e^{-\vec{\theta}\mathbf{x}}}}\\ P(y=0|\mathbf{x},\vec{\theta}) = 1-P(y=1|\mathbf{x},\vec{\theta}) \end{array} \tag{32}$$

有：

$$P(y|\mathbf{x},\vec{\theta})=P(y=1|\mathbf{x},\vec{\theta})^yP(y=0|\mathbf{x},\vec{\theta})^{1-y} \tag{33}$$

则似然函数为：

$$L(\vec{\theta})=\prod_{i=1}^NP(y=1|\mathbf{x}_i,\vec{\theta})^{y_i}P(y=0|\mathbf{x}_i,\vec{\theta})^{1-y_i} \tag{34}$$

从（34）可以看出，似然函数把所有给定观测值 $\{\mathbf{x}_1,\mathbf{x}_2,\ldots,\mathbf{x}_n\}$ 对决策边界函数的影响都考虑在内了，它表示在当前参数 $\vec{\theta}$ 下，模型和观测值的匹配程度。为了得到最大的似然函数并防止浮点数溢出，我们要将似然函数（34）取对数变成对数似然函数之后令其梯度为零（如果参数是一个标量就令导数为零）。

取对数：

$$\begin{aligned} l(\vec{\theta})&= \mathrm{ln}L(\vec{\theta}) \\ &=\sum_{i=1}^N\big[y_i\mathrm{ln}P(y=1|\mathbf{x}_i,\vec{\theta})+(1-y_i)\mathrm{ln}P(y=0|\mathbf{x}_i,\vec{\theta})\big] \end{aligned} \tag{35}$$

令梯度为零（这一步可以用之前提到的梯度下降法来计算）：

$$\nabla l(\vec{\theta})=0 \tag{35}$$

满足（35）的参数 $\vec{\theta}$ 即为所求。

逻辑回归的损失函数（交叉熵）

交叉熵损失函数（Cross Entropy Loss Function）

$$J(\vec{\theta})=-\frac{1}{N}\sum_{i=1}^N\big[y_i\mathrm{ln}P(y=1|\mathbf{x}_i,\vec{\theta})+(1-y_i)\mathrm{ln}P(y=0|\mathbf{x}_i,\vec{\theta})\big] \tag{36}$$

（36）是一个典型的二分类交叉熵损失函数

多层感知器网络（神经网络）

这部分内容参考了这位大佬的文章，感谢大佬提供的优秀思路 Orz…

神经网络主要用于数据的分类，以下面这个神经网络为例：

该神经网络主要分为：

• 输入层 X $(1\times2)$
• W1 $(2\times50)$
• 隐藏层 H
• W2 $(50\times4)$
• 输出层 Y $(1\times4)$

下面将逐步对神经网络的正向传播进行讲解：

输入层

输入层每一次会输入一个向量 $\mathbf{x}$ 经过矩阵 W1 的运算进入隐藏层：

$$\mathbf{H=X\ W1+b1} \tag{37}$$

隐藏层（激活层）

如果隐藏层直接连接 W2 然后就输出的话，我们就没必要费老大劲搭建神经网络了，直接拿一个线性式 $\mathbf{Y=X(W1+b1)(W2+b2)=X\ M}$ 表达就完了。真正让神经网络成为神经网络的是隐藏层中的激活层，就是这层激活层中的非线性函数让输入和输出之间的关系变为非线性的

假如经过公式 $\mathbf{H=X\ W1+b1}$ 计算得到的 H 值为：$(1,-2,3,-4,7...)$，那么经过阶跃函数激活层后就会变为 $(1,0,1,0,1...)$，经过 ReLU 激活层之后会变为 $(1,0,3,0,7...)$。

需要注意的是，每个隐藏层计算（矩阵线性运算）之后，都需要加一层激活层，要不然该层线性计算是没有意义的。

正规化处理

经过矩阵 W2 的运算进入输出层 Y 的值可能会是 $(3,1,0.1,0.5)$ 这样的矩阵，虽然可以通过比较这四个数的大小来找到 $\mathbf{x}$ 所属的分类，但是还不够直观，Y 的结果需要经过正规化处理即通过“Softmax”层。

简单来说分三步进行：（1）以e为底对所有元素求指数幂；（2）将所有指数幂求和；（3）分别将这些指数幂与该和做商。这样求出的结果中，所有元素的和一定为1，而每个元素可以代表概率值。

衡量输出的好坏

通过 Softmax 层之后，我们得到了I，II，III和IV这四个类别分别对应的概率，但是要注意，这是神经网络计算得到的概率值结果，而非真实的情况。

比如，Softmax输出的结果是 $(90\%,5\%,3\%,2\%)$，真实的结果是 $(100\%,0,0,0)$。虽然输出的结果可以正确分类，但是与真实结果之间是有差距的，一个优秀的网络对结果的预测要无限接近于100%，为此，我们需要将 Softmax 输出结果的好坏程度做一个“量化”。

一种直观的解决方法，是用1减去 Softmax 输出的概率，比如 $1-90\%=0.1$。不过更为常用且巧妙的方法是，求对数的负数。

还是用 90% 举例，对数的负数就是：$-\mathrm{ln}0.9=0.046$

可以想见，概率越接近100%，该计算结果值越接近于0，说明结果越准确，该输出叫做“交叉熵损失（Cross Entropy Error）”。（是的，就是3.3.3小节的那个）

我们训练神经网络的目的，就是尽可能地减少这个“交叉熵损失”。

至此，神经网络的正向传播过程介绍完毕，不知道俺说清楚没有，反正之前一听到神经网络就不明觉厉，但是在看完这位大佬的文章之后就两个字“通透”。。。

后日谈

啊~计算物理这个系列终于写完了呢，坚持下来真是太好了呢，如果能帮到大家就更好了呢~(*°ω°)ノ♡

其实计算物理是一个很有意思的学科（能敲代码的都有意思。。。吧？），而这门课程有趣不仅仅是因为敲代码，而是由于计算物理及其广泛（还会越来越广）的应用面，使得许多学过的知识都能和它串联起来，比如大一大二的微分方程、概率论的内容（老师讲的很好，是我太笨没听懂，责任全在我），能亲手把它们从课本上的公式变成一行行可编译执行的机器语言，我收获的不仅仅是运行程序的成就感，还有对各种方法计算原理的更深入的了解。

由于时间不够&我太菜了，有很多知识点我其实是没有写出来的 反正也不考~(*°▽°)，比如插值和拟合章节的最小二乘法拟合，还有蒙特卡洛方法里的马尔科夫链等，都是些很有意思的算法，有关文章有可能我以后会补上（坑挖好了嗷）。

总之，写完这些文章真是不容易呢，感谢俺自己（午饭还没吃呢）；感谢在我写这些文章的过程中提供帮助的老师和网页；感谢俺的播放器和耳机，感谢你们为我播放优美的音乐；最后，谢谢你能看到这里，感谢陪伴，完结撒花！ヾ(≧ v ≦*)ノ〃

有误之处欢迎反馈！