当DP遇上量子统计和量子动力学

DrugOne

发布于 2021-08-25 16:59:29

1.2K0

文章被收录于专栏：DrugOneDrugOne

对于我们研究的大部分物理化学问题来说，由于原子核的质量比电子大很多，忽略它们之间的耦合（波恩-奥本海默近似），用量子化学的方法处理分子体系中的电子结构部分，用经典力学处理其中的原子核部分已经能得到较好的结果。然而在一些特殊情况下原子核的量子效应的贡献不能忽略，例如含有轻原子（如H）的体系、温度较低的情形或者涉及到多电子态的非绝热过程，此时需要使用量子力学的手段研究分子体系的统计或动力学性质。

统计热力学

当我们研究一个（粒子数）、（体积）、（温度）恒定的体系的统计力学性质时，经典统计力学中玻尔兹曼分布占据了重要的地位

\rho(x,p) \propto e^{- \beta U(x) - \beta p^T M^{-1} p / 2}

其中，是该体系的势能面，来自第一性原理计算或DeePMD-kit拟合得到的机器学习势能面，而动量满足高斯分布。只要确定下来体系的势能面的形式，就可以借助蒙特卡洛（MC）或分子动力学（MD）的手段对玻尔兹曼分布进行采样，而热力学量可以表达成某些物理量在玻尔兹曼分布下的平均：

\langle B \rangle = \int \rho (x,p) B(x,p) \text{d}x \text{d}p

而在量子统计力学中，由于海森堡不确定性关系，没有唯一确定的构型和动量的联合分布，类似的对应形式是玻尔兹曼算符，而热力学量可以表达成物理量与玻尔兹曼算符在求迹意义下的平均：

\langle B \rangle = \text{Tr} [\hat{\rho} \hat{B}] / Z

其中是配分函数。上面的求迹可以表达在构型表象下，例如当是构型的本征算符时：

\langle B \rangle = \int \rho(x,x) B(x) \text{d} x / Z

而对于一般情形的算符：

\langle B \rangle = \int \rho(x,y) B(x,y) \text{d} x \text{d} y / Z

得益于费曼提出的路径积分，玻尔兹曼算符在构型表象下的矩阵元可以通过对端点为的各种路径进行积分得到。离散情况下相当于在间插入若干体系的拷贝，除了每个体系受到势能面的作用力外，相邻两个体系间还受到谐振子弹簧力的约束，把所有辅助变量积分掉即可得到的边缘分布，随着插入辅助变量的个数趋于无穷多，路径积分的结果逐渐收敛于量子力学精确结果：

\rho(x,y) \propto \lim_{P \to \infty} \int e^{ - \beta U_{PI} (x_1 \equiv x, x_2, \cdots, x_P, x_{P+1} \equiv y)} \text{d} x_2 \cdots \text{d} x_P

通常，人们会对带辅助变量的扩展体系进行MC或MD采样，称为路径积分蒙特卡洛（PIMC）或路径积分分子动力学（PIMD）。也可以定义关于的有效势能或关于的有效势能，它们可以看成平均力的势能（potential of mean force）或粗粒化势能。可以通过DeePCG的手段对它进行机器学习建模，再将PIMC/PIMD中采样得到的瞬时力作为数据对模型进行训练。有了优化收敛的模型，仅需对一个变量或两个变量在有效势能下进行MC/MD采样就能得到体系的量子结果。下图验证了丙二醛分子在有效势能下采样可以获得与原始路径积分相同的分布。

动力学

许多动力学性质可以表达成时间关联函数的形式，例如热导率与能流关联函数相关，反应速率与流关联函数相关，红外/拉曼光谱与偶极/极化率关联函数相关等。经典力学中关联函数具有如下形式：

C_{AB} (t) = \int \rho (x_0,p_0) A (x_0,p_0) B(x_t,p_t) \text{d} x_0 \text{d} p_0

其中是玻尔兹曼分布，是以为起点按照经典牛顿力学（实时间）演化至时刻的结果，演化过程中的（粒子数）、（体积）、（能量）保持不变。

时间关联函数在量子力学中的对应形式是：

C_{AB} (t) = \text{Tr} [\hat{\rho} \hat{A} \hat{B} (t)] / Z

与热力学性质不同，研究动⼒学性质更为复杂，实时间的路径积分由于相位震荡引起的“符号”问题难以应⽤于⼤体系，能得到精确结果的体系不超过⼏个原⼦。要研究⼤体系的动⼒学性质就不得不转⽽发展计算代价可以接受的近似⽅法。一类近似方法是从量子力学的相空间表象出发，仿照经典关联函数的形式，用轨迹的系综平均代替算符求迹：

C_{AB} (t) \approx \int \rho_W (x_0,p_0) f_A^W (x_0,p_0) B_W(x_t,p_t) \text{d} x_0 \text{d} p_0

其中与经典的玻尔兹曼分布有所不同，是量子相空间分布函数。是以为起点在某种动力学的支配下演化至时刻的结果。然而上式始终只能是近似，无论初始分布如何选择，动力学如何选取，都无法对于一般情况得到精确的量子力学结果。

既然无法得到完全精确的结果，我们转而考虑在一些极限下能够逼近它。两个重要的极限包括，平衡态应该保证任何热力学量统计平均的动力学演化保持守恒，而对谐振子体系（即体系势能面是构型的二次函数）应该获得精确的关联函数。一类满足这两种要求的动力学方法是刘剑和W. H. Miller等人提出的平衡连续性动力学（equilibrium continuity dynamics, ECD）1，用轨迹动力学所满足的连续性方程作为限制，来保证密度分布随动力学演化不变。在ECD的运动方程中，除了涉及上面提到的有效势能外，还有一个重要元素是热质量矩阵，它是一个的矩阵，作为构型的函数，它满足平移不变性，旋转协变性和交换协变性。如果在动力学演化过程中每个时间步长都实时地（on the fly）用路径积分计算和，代价十分高昂，引入机器学习对它们进行高效建模可以把计算量降低6个数量级以上2。

为了构造热质量矩阵的模型，需要在Deep Potential基础上稍加推广，核心是从每个原子的环境中提取出的张量特征，该特征满足对构型的平移不变性和旋转协变性，再将第个和第个原子的张量特征经过线性变换获得热质量矩阵的第行第列的块。相比Deep Potential中的标量特征描述符，张量特征描述符仅需将向量内积换成向量外积：

\overleftrightarrow{D}_i^{lm} = \sum_{j,k} G_{\alpha_i \alpha_j}^l(\left|r_{ij}\right|) G_{\alpha_i \alpha_k}^m(\left|r_{ik}\right|) (r_{ij} r_{ik}^T + r_{ik} r_{ij}^T) / 2

然而在应用中我们发现，标量特征描述符中局域嵌入网络只保留径向信息作为输入而忽略角度信息的做法对张量特征描述符可能造成较大影响，例如对3原子分子水分子来说，当两个O-H键长相等即时，用上式表达的张量特征为：

\overleftrightarrow{D}_O^{lm} = G_{HH}^l(r) G_{HH}^m(r) (r_1 r_1^T + r_1 r_2^T + r_2 r_1^T + r_2 r_2^T) \\ = G_{HH}^l(r) G_{HH}^m(r) (r_1 + r_2) (r_1+ r_2)^T

可见其秩退化为了1，引入了多余的对称性，如果用这样的描述符来构造张量性质，拟合结果会随着模型参数的增加变得越来越不光滑。此时运动方程积分时的离散误差会快速累积导致动力学轨迹无法继续演化。在局域嵌入网络中加入角度信息能缓解这个问题

\overleftrightarrow{D}_i^{l} = \sum_{j,k} G_{\alpha_i \alpha_j \alpha_k}^l(\left|r_{ij}\right|, \left|r_{ik}\right|, \cos \theta_{ijk}) (r_{ij} r_{ik}^T + r_{ik} r_{ij}^T) / 2