数值优化和几何优化（密度泛函理论）

最新推荐文章于 2024-04-03 22:14:37 发布

hhjc

最新推荐文章于 2024-04-03 22:14:37 发布

阅读量892

点赞数 1

分类专栏： MS CASTEP的理论基础-密度泛函理论文章标签：数值优化几何优化

本文链接：https://blog.csdn.net/weixin_46465199/article/details/129928073

版权

MS CASTEP的理论基础-密度泛函理论专栏收录该内容

7 篇文章 8 订阅

订阅专栏

1数值优化

为了能够在实际操作中实现数值收敛的DFT计算，我们需要借助于一些数值优化方法，来有效处理具有多个自由度的问题。如同理解倒易空间一样，充分理解优化的核心概念对有效进行DFT计算十分重要。下面将主要了解相关内容。

1.1一维空间中的优化

先看一个运行DFT计算比较关键的问题：从数学问题入手：找到 $f\left ( x \right )$ 的一个极小值，其中 $f\left ( x \right )$ 是由单值确定的某个函数。假设 $f\left ( x \right )$ 是一个“光滑函数”，也就是说该函数的导数是存在并且连续的。可以认为 $f\left ( x \right )$ 所表征的是一个非常复杂的函数，无法用代数问题对其进行求解。换言之，我们必须使用计算机来求得该问题的数值解。我们所要求的极值是局域的，并不是整个区间的极值。我们主要是求解一个局域化问题，并非一个全局优化问题。

(极值处理问题在DFT计算中十分重要。对于确定FCC金属的晶格常数这一问题，就可以将其作为一个最小化问题。如果我们并不知道超晶胞中N个原子的精确排布构型，那么这个超晶胞的总能可以写成 $E\left ( x \right )$ ，其中是定义原子位置的维矢量。找到这些原子在物理上最为有限的排布构型就相当于将这个能量函数最小化。)

继续上面的数学问题，假如能够找到满足 $f{}'\left ( x \right )=df/dx=0$ 这一条件的值，那么这个点就确定了该函数的一个极大值或者极小值点，所以能够把上述问题重新定义为：找到一个满足 $f{}'\left ( x \right )=0$

的值。下面采用两种方法来求解这一问题。

1二分法

找函数零点的话首先要确保存在 $f{}'\left ( x_{1} \right )f{}'\left ( x _{2}\right )< 0$ 。只要 $f{}'\left ( x \right )$ 是光滑的，那么就一定能在 $x_{1},x_{2}$ 之间找到一点使 $f{}'\left ( x \right )$ 等于0。但是由于不知道这一点在哪，需要对两者的中点考察一下 $f{}'\left ( x^{\ast } \right )$ 。此时，可能会出现三种情况，要么 $f{}'\left ( x^{\ast } \right )=0$ ,要么 $f{}'\left ( x_{1} \right )f{}'\left ( x^{\ast } \right )< 0$ ,要么 $f{}'\left ( x^{\ast } \right )f{}'\left ( x_{2} \right )< 0$ 。除了第一种情况，剩余两种情况都将继续寻找零点下去，并且只要有足够长的时间重复计算，这个 $x^{\ast }$ 就一定会越来越逼近所需要寻找的零点。

2牛顿法

该方法思路如上图所示。定义 $g\left ( x \right )=f{}'\left ( x \right )$ ,那么根据泰勒展开式，则有 $g\left ( x+h \right )\approx g\left ( x \right )+hg{}'\left ( x \right )$ 。这是忽略高阶项的近似，也很好理解，其实这个近似在h无限小的情况下可以看成梯形图解法求积分的近似，加上h(宽）乘以（ $g{}'\left ( x \right )$ ）(高)。当在 $g\left ( x \right )\neq 0$ 的位置考察 $g\left ( x \right )$ 时，就可以根据上述近似方程，得到该函数为0时的一个良好估测值： $x^{\ast }=x+h=x-g\left ( x \right )/g{}'\left ( x \right )$ 。这个式子如果按上面近似替换的话并不能连等，可以按这种思路来思考：牛顿法求零点、极值点(66条消息) 牛顿法求零点、极值点_sunlanchang的博客-CSDN博客

因为泰勒展开式舍去了高阶项，所以并不能精确确定零点位置。就如同二分法那样，也需要进行大量重复的计算来逼近这个值。

举个栗子：

使用牛顿法从x=1.8开始寻找 $f\left ( x \right )=e^{-x}cosx$ 的极小值，可以生成以下一系列近似解： $x^{\ast }=1.8,2.183348,2.331986,2.355627,2.356194$ （保留六位小数）。

对于二分法和牛顿法，二者显著区别在于收敛速度。下图进行比较：

现在值了解了在足够多次数计算后越来越逼近真实解，那么如何判断何时停止计算？在实际应用中，我们并不确定真实解，因此无法用 $\varepsilon _{i}=\left | x_{i} -x\right |$ 来表述当前估值多么精确。为了能有一个合理的判断条件，需要定义一个计算终止标准。常用方法就是：在连续两次迭代之间，如果两者差值小于某个容限值时 $\left | x_{i+1}-x_{i} \right |<\delta$ ，就终止迭代计算。

上面栗子体现了数值优化的几个普遍性特点：

$\left ( 1 \right )$ 这些算法都是迭代的，因此他们无法给出精确解，但是能提供一些列精确解近似值

$\left ( 2 \right )$ 使用这些算法时，都必须提供一个真实解的初始估值。

$\left ( 3 \right )$ 所运算的迭代次数由容限值参数控制，当然也取决于初始估值，这个容限值估测了当前值与真实解的接近程度。

$\left ( 4 \right )$ 使用不同的容限值或者不同的初始估值，采用任何一种算法进行多次迭代计算都最终逼近与真实解的近似，虽然可能形成的一系列解并不相同。

$\left ( 5 \right )$ 使用不同的算法，收敛于真实解的速度也不相同。因此选择合适的算法十分重要。

但是：

$\left ( 1 \right )$ 如果只进行一次计算，这能找到对于这次计算的区间的极值，无法判断函数是否存在多个极小值。即使采用多次设置初始值进行迭代，也不能找到全局的多个极小值。

$\left ( 2 \right )$ 对于大多数方法，无法保证该方法能从任意一个初始估值收敛于真实解。

1.2大于一维的优化

回到前面的问题：通过不断改变超晶胞的原子位置，对超晶胞中一组原子的总能 $E\left ( x \right )$ 进行最小化。这个例子将会凸显出在试图求解多变量优化时候一些复杂难题，这些难题在一维优化并没有体现。如果定义 $g_{i}\left ( x \right )=\partial E\left ( x \right )/\partial x_{i}$ ,那么最小化 $E\left ( x \right )$ 就相当于找到一组原子位置，使 $g_{i}\left ( x \right )=0$ 对于 $i=1,\cdots ,3N$ 都同时成立。

把一维的对分法推广用于求解这个多为问题很难。但是牛顿法推广可以应用于这个情形。一维牛顿法是根据泰勒展开式推导得到，对于多维问题也可以用相同思路处理。如果结果中含有 $3N\times 3N$ 的导数矩阵,其元素为 $J_{ij}=\partial g_{i}/\partial x_{i}$ 。请注意，该矩阵的元素是我们真正所需要在意函数 $E\left ( x \right )$ 的二阶偏导数。牛顿法通过下式定义了一系列迭代，即：

$x_{i+1}=x_{i}-J^{-1}\left ( x_{i} \right )g\left ( x_{i} \right )$ （1）

这个应用范围更广。带入一个情形：假如有十个原子，在立体三维中就有30个变量。每次迭代中，都必须计算矩阵的55个不同元素；其他45个矩阵元素随之而定。然后就是求解 $30\times 30$ 矩阵的线性代数问题。计算量相比于一维骤增。可以看出牛顿法的计算量增大并不适合进行多维问题的求解。

一维牛顿法和多维牛顿法：(69条消息) 一维牛顿法与多维牛顿法_多维度牛拉法_一条兔子的博客-CSDN博客

下面两种方法是常用的多维问题数值方法：拟牛顿法 Quasi-Newton Method 和共轭梯度法（ Conjugate-Gradient Method ）。

拟牛顿法

拟牛顿法和只要求每一步迭代时知道目标函数的梯度。通过测量梯度的变化，构造一个目标函数的模型使之足以产生超线性收敛性。这类方法大大优于最速下降法，尤其对于困难的问题。另外，因为拟牛顿法不需要二阶导数的信息，所以有时比牛顿法(Newton's Method)更为有效。如今，优化软件中包含了大量的拟牛顿算法用来解决无约束，约束，和大规模的优化问题。

拟牛顿法是解非线性方程组及最优化计算中最有效的方法之一.它是一类使每步迭代计算量少而又保持超线性收敛的牛顿型迭代法。

拟牛顿法：

(69条消息) 拟牛顿法_一条兔子的博客-CSDN博客

共轭梯度法：

(69条消息) 共轭梯度法(Conjugate gradient)详解_bitcarmanlee的博客-CSDN博客

关于数值优化最主要还是先理解上面总结7点，然后具体问题具体分析。

2几何优化

2.1内部自由度

试想一下，假如我们所感兴趣的是与氮有关的一组反应，那么必定有一件事是我们想知道的，就是气相 $N_{2}$ 分子的能量，为了计算这个能量，我们需要找到其分子能量最小化的 $N_{2}$ 构型。

模拟一个气相分子，一个简单方法就是建立一个边长为的立方超晶胞。然后把两个原子放在这个超晶胞的分数坐标 $\left ( 0,0,0 \right )$ 和 $\left ( +d/L,0,0 \right )$ 位置上。只要显著大于的长度，这个超晶胞就表示了一个键长为的一个隔离态 $N_{2}$ 分子。采用拟牛顿法或者共轭梯度法，通过最小化上述超晶胞的总能（在最小化过程中，固定超晶胞的尺寸和形状，但是允许两个原子的分数坐标自由变动）就可以得到 $N_{2}$ 分子 DFT 优化后的键长。此外，我们必须重新定义一个终止标准，用以决定是否这些迭代已经收敛到最小值。由于在计算的 $N_{2}$ 原子构型中，两个原子受力平衡，所以设定两个原子上的力数量级小于 0.001eV/A ,就终止迭代计算。在这个例子中，由于对称性，在两个原子上的力数量级始终是完全相同的，当然，两个原子上面的力分别指向两个相反的方向。这个标准设置是否合理？如果我们把一个原子的位置稍稍改变一个很小了量（ $\Delta r$ ）,由于这一变化所导致的总能量变化就可以估算为 $\left | \Delta E \right |\approx \left | F \right |\left | r \right |$ ,其中 $\left | F \right |$ 是作用在原子上面的力。如果两个原子上面的力都小于 0.01eV/A ,那么，当把任何一个原子移动 0.1A 时（相对于化学键长，这是一个较大的距离），则总能变化都小于 0.001ev ，这是一个很小的能量。可以看出，这个迭代终止标准设置是合理的。

所需要考虑的第二个因素就是：在初始放置两个原子时，两者应该距离多远。这一项似乎并不是很重要，因为在进行能量最小化，只有一个键长对应于体系最小能量。但是，对于键长初始估值的选择，也会在很大程度上影响计算是否能顺利完成。就比如，我们使用共轭梯度法进行计算，使用两个不同初始距离带入。当把原子间距初始值设置成时，计算进行11次迭代之后，力已经满足了上述收敛标准，并且最终所得键长为 1.12A 。实验上观测到的键长为 1.10A ，可以看出DFT求解虽然不是精确解，但是和实际差别也很小。在第二次计算时，初始原子间距设置为 0.7A 。此时，经历共次迭代之后，两个原子之间的间距为 2.12A ，这个迭代情况很不理想，并且能量还没有进行收敛。为什么会这样？这是由于这个初始状态，当两个间距设置成 0.7A 时，就相当于物理上对其键长进行了极大压缩（相比于 $N_{2}$ 分子键长的正常值）。这个值意味着存在会存在一个极大的排斥力，把两个原子互相推开。这会产生很大影响，因为在这个优化方法中，要根据某个位置上的能量导数估测能量变化有多快。所编写的优化程序在进行判断时首先把它判定为需要将原子移开，这回导致程序初始步长会较大，可能直接超过了最小能量体系的原子间距，出现震荡现象，导致长时间不收敛。这说明此次计算没有成功。

关于怎样从一些细节判断是否已经出现了这类失误，与所使用的优化算法相关；但是，对于任意一种几何优化算法，如果对原子构型使用了很小的初始值，就会出现上述这种情况。

给出一个经验：花费较大的计算工作量，对原子坐标给出一个良好的初始估值，可以极大的提高DFT计算速度，这个值对最终的收敛也会有影响。

第二个例子:

优化一个 $CO_{2}$ 分子的几何构型。再次使用一个边长为L的立方超晶胞，我们可以通过下列方法创建一个 $CO_{2}$ 分子：把原子放置在分数坐标 $\left ( 0,0,0 \right )$ ，两个氧原子放置在 $\left ( +d/L,0,0 \right )$ 和 $\left ( -d/2,0,0 \right )$ 。从 d=1.3A 的初始状态开始优化该分子的能量，并采用与 $N_{2}$ 分子相同的迭代终止标准，最终得到优化后的结果是 C-O 键长是 1.17A ,且 O-C-O 键角是 $180^{\circ}$ 。

这个结果看似很合理，但是考察一下这个原子构型中作用在原子上面的力，如果我们把这个力写成 $f=\left ( f_{x} ,f_{y},f_{z}\right )$ ，那么根据对称性，无论我们怎么对取值，都有 $f_{y}=f_{z}=0$ 。这说明：在能量最小化过程中，分子的几何构型虽然被迭代，但是 O-C-O 键角始终保持在 $180^{\circ}$ ，也就是我们在初始构型中所给定的角度值。这存在一个问题就是，这个键角很有可能并不是真实体系能量最低对应的键角，并非一个真实收敛结果。

优化这一个分子更加可靠的方法就是设置初始键角的时候，不要使两边的分量由于对称性而导致分量为0。

2.2具有约束原子的几何优化

优化一部分原子位置，其他原子原值保持不变，在上述前提下进行超晶胞能量最小化，这种计算也很有现实意义。相比于无约束优化，牵扯到约束的数值优化问题在计算上更具有挑战性。基于力的几何优化，可以很容易的扩展到一个或者多个原子都被约束的情形。对于由约束的体系，在优化计算的每次迭代过程中，只有那些未约束的原子的位置会不断更新。