2.10 几何优化和从头算分子动力学

在玻恩——奥本海默近似下，原子核由经典力学描述，电子根据与每个原子构型相关的基态波函数从属于原子核。总能量E（{RI}）提供了原子间势能的第一原理描述。
几何优化（或结构优化）问题旨在寻找势能面E（{RI}）上的（全局或局部）极小值，相应的极小值{R∫I}是（全局或局部）稳定构型。大多数自然存在的分子和固体都处于或接近这样的极小值。
定义第i个原子核的原子力为
F I ( { R I } ) = − ∂ E ( { R I } ) ∂ R I ; F_I(\{R_I\})=-\frac{\partial E(\{R_I\})}{\partial R_I}; FI​({RI​})=−∂RI​∂E({RI​})​;

那么，确定极小值的一个必要条件是，对于所有原子，原子力$F_I$都应该为零。
E t o t = ∑ I = 1 M M I R ˙ I 2 2 + E ( { R I } I = 1 M ) , E_{\mathrm{tot}}=\sum_{I=1}^M\frac{M_I\dot{\boldsymbol{R}}_I^2}2+E(\{\boldsymbol{R}_I\}_{I=1}^M), Etot​=I=1∑M​2MI​R˙I2​​+E({RI​}I=1M​),

其中MI是第i个原子核的质量。
原子核的动力学性质可以用牛顿方程来研究
M I R ¨ I ( t ) = − ∂ E ∂ R I ( { R I ( t ) } ) = F I ( { R I ( t ) } ) . M_I\ddot{\boldsymbol{R}}_I(t)=-\frac{\partial E}{\partial\boldsymbol{R}_I}(\{\boldsymbol{R}_I(t)\})=\boldsymbol{F}_I(\{\boldsymbol{R}_I(t)\}). MI​R¨I​(t)=−∂RI​∂E​({RI​(t)})=FI​({RI​(t)}).

这被称为从头算分子动力学（AIMD）。几何优化和分子动力学消除了构建经验原子间势的需要，并使化学和材料科学中静态和动态性质研究的广泛应用成为可能。

在几何优化和AIMD中，关键因素是评估原子力。在三维空间中，这种导数信息正式需要3M独立电子结构计算，这将是极其昂贵的。幸运的是，赫尔曼——费曼定理允许我们评估原子力，与单个电子结构计算相比，几乎没有增加成本。

原子力的评估

为了引入赫尔曼——费曼定理 ，让我们首先考虑一个参数相关的哈密顿量H（λ）。对于每个给定的λ，我们假设基态是非简并的。它的特征值ε（λ）和本征函数$\psi (r;\lambda )$可以计算为
$$H(\lambda )\psi (\lambda )=\epsilon (\lambda )\psi (\lambda )\text{\hspace{1em}(2.10.4)}$$
特征值对λ的导数是
$$\begin{array}{rl}\frac{\mathrm{d}\epsilon (\lambda )}{\mathrm{d}\lambda }& =\frac{\mathrm{d}}{\mathrm{d}\lambda }⟨\psi (\lambda )|H(\lambda )|\psi (\lambda )⟩\\ & =⟨\psi (\lambda )|\frac{\mathrm{d}H(\lambda )}{\mathrm{d}\lambda }|\psi (\lambda )⟩+\epsilon (\lambda )\frac{\mathrm{d}}{\mathrm{d}\lambda }⟨\psi (\lambda )|\psi (\lambda )⟩\\ & =⟨\psi (\lambda )|\frac{\mathrm{d}H(\lambda )}{\mathrm{d}\lambda }|\psi (\lambda )⟩.\end{array}$$
这里我们对所有λ使用了归一化条件<φ（λ）|φ（λ）>=1。因此，赫尔曼——费曼定理表明，特征值相对于外部参数扰动的一阶导数简单地由哈密顿量相对于λ的变化给出。
赫尔曼——费曼定理的另一个观点可以从变分原理中理解
$$\epsilon (\lambda )=\underset{|\psi )\ne |0⟩}{\inf }E[\psi ,\lambda ]:=\underset{|\psi ⟩\ne |0⟩}{\inf }\frac{⟨\psi |H(\lambda )|\psi ⟩}{⟨\psi |\psi ⟩}.$$
因此，
$$\frac{\mathrm{d}\epsilon (\lambda )}{\mathrm{d}\lambda }={\frac{\partial E[\psi ,\lambda ]}{\partial \lambda }|}_{\psi =\psi (\lambda )}+\int {\frac{\delta E[\psi ,\lambda ]}{\delta \psi (r)}|}_{\psi =\psi (\lambda )}\frac{\partial \psi (r;\lambda )}{\partial \lambda }\mathrm{d}r=⟨\psi (\lambda )|\frac{\mathrm{d}H(\lambda )}{\mathrm{d}\lambda }|\psi (\lambda )⟩.$$
请注意，虽然导数$\frac{\partial \psi (r;\lambda )}{\partial \lambda }$一般不会消失，但由于欧拉——拉格朗日方程，它对$\frac{\mathrm{d}\epsilon (\lambda )}{\mathrm{d}\lambda }$没有贡献

现在，对于Kohn-Sham DFT中的总能量，只有外部势Vext和离子——离子相互作用EII明确依赖于原子构型{RI}。此外，外部势一般可以分解为每个原子的贡献，如
V e x t ( r ; { R I } ) = ∑ I V e x t I ( r − R I ) . V_\mathrm{ext}(r;\{R_I\})=\sum_IV_\mathrm{ext}^I(\boldsymbol{r}-\boldsymbol{R}_I). Vext​(r;{RI​})=I∑​VextI​(r−RI​).
然后，应用赫尔曼——费曼定理，我们发现
F I = − ∂ E ( { R I } ) ∂ R I = − ∫ ρ ( r ) ∂ V e x t I ( r − R I ) ∂ R I   d r − ∂ E I I ( { R I } ) ∂ R I . \boldsymbol{F}_I=-\frac{\partial E(\{\boldsymbol{R}_I\})}{\partial\boldsymbol{R}_I}=-\int\rho(\boldsymbol{r})\frac{\partial V_\mathrm{ext}^I(\boldsymbol{r}-\boldsymbol{R}_I)}{\partial\boldsymbol{R}_I}\mathrm{~d}\boldsymbol{r}-\frac{\partial E_\mathrm{II}(\{\boldsymbol{R}_I\})}{\partial\boldsymbol{R}_I}. FI​=−∂RI​∂E({RI​})​=−∫ρ(r)∂RI​∂VextI​(r−RI​)​ dr−∂RI​∂EII​({RI​})​.

Car-Parrinello分子动力学(CPMD)

在AIMD中，电子结构问题必须完全自洽地解决，以便达到每个原子构型的基态。因此，这种类型的AIMD也被称为玻恩——奥本海默分子动力学（BOMD）。很长一段时间，尽管BOMD有明显的潜力，但它被认为过于昂贵。AIMD是由开创性的Car-Parrinello分子动力学（CPMD）[14]实现的，它引入了一个扩展的拉格朗日函数，包括原子核和电子的自由度：
L c p ( { R l } , { R ˙ l } , { ψ i } , { ψ ˙ i } ) = ∑ I M I R ˙ I 2 + ∑ i μ 2 ∫ ψ ˙ i ∗ ( r ) ψ ˙ i ( r ) d ⁡ r − E ( { ψ i } ; { R I } ) . (2.10.13) \begin{aligned}\mathcal{L}_{\mathrm{cp}}(\{R_l\},\{\dot{R}_l\},\{\psi_i\},\{\dot{\psi}_i\})\\&=\sum_I\frac{M_I\dot{R}_I}2+\sum_i\frac\mu2\int\dot{\psi}_i^*(\boldsymbol{r})\dot{\psi}_i(\boldsymbol{r})\operatorname{d}\boldsymbol{r}-E(\{\psi_i\};\{\boldsymbol{R}_I\}).\end{aligned}\tag{2.10.13} Lcp​({Rl​},{R˙l​},{ψi​},{ψ˙​i​})​=I∑​2MI​R˙I​​+i∑​2μ​∫ψ˙​i∗​(r)ψ˙​i​(r)dr−E({ψi​};{RI​}).​(2.10.13)

这里μ被称为虚拟电子质量。E（{φ i}；{RI}）是Kohn-Sham DFT中的能量泛函，而{φ i}可能是也可能不是原子构型{RI}的最小化Kohn-Sham轨道。

由拉格朗日方程（2.10.13）引起的运动方程给出了CPMD
M I R ˙ I ( t ) = − ∫ ρ ( r , t ) ∂ V e x t I ( r − R I ( t ) ) ∂ R I dr ⁡ − ∂ E I I ∂ R I ( { R I ( t ) } ) , I = 1 , … , M , μ ψ i ¨ ( t ) = − H [ ρ ( t ) ; { R I ( t ) } ] ψ i ( t ) + ∑ j ψ j ( t ) Λ j i ( t ) , i = 1 , … , N , ρ ( r , t ) = ∑ i = 1 N ∣ ψ i ( r , t ) ∣ 2 . \begin{aligned} M_{I}\dot{\boldsymbol{R}}_{I}(t)& =-\int\rho(\boldsymbol{r},t)\frac{\partial V_\mathrm{ext}^I(\boldsymbol{r}-\boldsymbol{R}_I(t))}{\partial\boldsymbol{R}_I}\operatorname{d\boldsymbol{r}}-\frac{\partial E_\mathrm{II}}{\partial\boldsymbol{R}_I}(\{\boldsymbol{R}_I(t)\}),\quad I=1,\ldots,M, \\ \mu\ddot{\psi_{i}}(t)& =-H[\rho(t);\{\boldsymbol{R}_I(t)\}]\psi_i(t)+\sum_j\psi_j(t)\Lambda_{ji}(t),\quad i=1,\ldots,N, \\ \rho(\boldsymbol{r},t)& =\sum_{i=1}^N\lvert\psi_i(r,t)\rvert^2. \end{aligned} MI​R˙I​(t)μψi​¨​(t)ρ(r,t)​=−∫ρ(r,t)∂RI​∂VextI​(r−RI​(t))​dr−∂RI​∂EII​​({RI​(t)}),I=1,…,M,=−H[ρ(t);{RI​(t)}]ψi​(t)+j∑​ψj​(t)Λji​(t),i=1,…,N,=i=1∑N​∣ψi​(r,t)∣2.​

这里，Λ是拉格朗日乘子，使得{φ i（t）}在任何时间都是一个正交函数集。与BOMD相比，CPMD使用虚拟动力学来引导电子的运动，而不需要收敛的SCF迭代。电子轨道的动力学可以粗略地看作是在每个分子动力学步骤执行SCF迭代的一种特殊方式。由于哈密顿结构，CPMD的数值模拟是稳定的，并且与没有严格SCF收敛准则的BOMD相比，能量在更长的时间内是保守的。当系统有光谱间隙时，CPMD的精度由一个参数控制，即虚拟电子质量μ。当μ趋于零时，CPMD的结果接近BOMD的结果[69，9]。然而，也有研究表明，CPMD不适用于间隙消失的系统，例如金属系统[69]。