分子动力学中优化算法和积分算法

文章目录

- 1.优化算法
- - 1.优化算法
  - 2.最速下降法
  - 3.一维搜索算法
  - 4.共轭梯度法
  - 5.多原子体系构型优化
- 2.积分算法
- - 1.积分算法
  - 2.前向Euler算法
  - 3.Verlet算法
  - 4.Velocity Verlet算法
  - 5.Leapfrog算法
  - 6.积分算法的稳定性

1.优化算法

优化分为局部最优化和全局最优化(网格搜索、模拟退火、遗传算法等)，而凸函数的局部最优解就是全局最优解
局部优化分为：
间接法，利用目标函数的一阶或二阶导数，如最速下降法、共轭梯度法、牛顿法
直接法，利用目标函数值，如坐标轮换法、鲍威尔法
间接法的收敛速率更高，直接法的可靠性较高
分子的势能面通常是一个复杂的依赖于坐标的多维函数，需要通过优化算法计算得到分子的最优构型 $E_{total}=E(x_{1},y_{1},z_{1},x_{2},y_{2},z_{2},...,x_{N},y_{N},z_{N})$
优化算法目标是找到函数的最大值/最小值，一个系统通常有很多的局部最小值，优化将前往最近的局部最小值
优化算法流程：给定初始点 $x^{0}$
1.在迭代点 $x^{k}$ 处，按一定规律确定搜索方向 $p^{k}$
2.沿搜索方向确定适当的搜索步长 $\alpha_{k}$
3.令 $x^{k+1}=x^{k}+\alpha_{k}p^{k}$ ，
a.若 $x^{k+1}$ 满足某种终止条件，则停止迭代，得到近似最优解 $x^{*}$
b.否则，重复以上步骤
初始点 $x^{0}$ 、搜索方向 $p^{k}$ 、迭代步长 $\alpha_{k}$ 称为优化算法的三要素
其中搜索方向最为重要和突出，确定搜索方向是研究优化算法的最根本任务之一
负梯度方向是函数值在该点附近的范围内下降最快的方向

2.最速下降法

二次型是 n 个变量上的二次齐次多项式 $f(x)=\frac{1}{2}x^{T}Ax-b^{T}x+c$ 梯度表达： $\nabla f(x)=Ax-b$
二次型是除线性函数外最简单的函数，具有对称矩阵表示形式，二次型的优化是研究优化问题的基础
二次型的最速下降法
正定二次型的最优步长可以使用解析法得到 $\alpha_{k}=\underset{\alpha}{argmin}\{ \varphi(\alpha)=f(x^{k}+\alpha p^{k})\}$ 根据方向导数定义： $\nabla_{p}f=\frac{\partial f}{\partial x}p_{x}+\frac{\partial f}{\partial y}p_{y}+\frac{\partial f}{\partial z}p_{z}$ $f(x)=\frac{1}{2}x^{T}Ax-b^{T}x+c$ 在 $(x^{k}+\alpha p^{k})$ 点沿着负梯度方向 $p^{k}$ 的导数为： $\varphi^{'}(\alpha)=(x^{k}+\alpha p^{k})^{T}Ap^{k}-b^{T}p^{k}$ 根据极值条件， $\alpha$ 使得该函数取极小值
令 $\varphi^{'}(\alpha)=0$ ，即： $\alpha(p^{k})^{T}Ap^{k}=-((x^{k})^{T}A-b^{T})p^{k}$ 而 $x^{k})^{T}A-b^{T}=(-p^{k})^{T}$ ，得最优步长 $\alpha$ $\alpha=\frac{(p^{k},p^{k})}{(Ap^{k},p^{k})}$

Python代码实现

def steepest_descent(x0,tol=1e-5,max_iter=1000)
    x = x0
    for i in range(max_iter):
        p = -grad_f(x)
        if np.linalg.norm(p) < tol:
             break
        alpha = p.T @ p / (p.T @ A @ p)
        x = x + alpha * p
    return x

最速下降法的缺点：随着迭代步的继续，搜索方向呈现震荡的行为 $x^{k+1}=x^{k}+\alpha p^{k}$ $\frac{d}{d\alpha}f(x^{k+1})=f^{'}(x^{k+1})^{T}\frac{d}{d\alpha}x^{k+1}=f^{'}(x^{k+1})^{T}\cdot p^{k}\equiv 0$ 在最速下降法中，相邻两个迭代点上的函数梯度相互垂直。而搜索方向是负梯度方向，相邻两个搜索方向互相垂直

3.一维搜索算法

对于一般的函数形式，往往依赖数值解法来确定最佳步长 $\alpha_{k}=\underset{\alpha}{argmin}\{ \varphi(\alpha)=f(x^{k}+\alpha p^{k})\}$ 分类：精确一维搜索(进退法、分割法、插值法)、非精确一维搜索(Backtracking-Armijo)
一维搜索的步长确定：
1.确定搜索区间，使函数 f(x) 存在一个局部最小值
a.图行方法：绘制函数图像，通过观察来估计最小值可能存在的区间
b.分析方法：使用函数的导数或其他数学工具进行分析，以确定可能存在最小值的区间
c.进退法(Bracketing Method)：从一个初始点出发，通过增加或减小 x 的值来找到一个包含最小值的区间
2.缩小搜索区间，确定最小值的精确位置
a.黄金分割搜索：使用黄金比率来选点，并根据两个点的函数值来更新区间
b.二分法：取区间的中点，然后根据中点的函数值或导数值来决定保留哪半个区间
c.插值方法：使用插值多项式来逼近 f(x)，然后选择插值多项式的最小值作为新的点
d.牛顿法：如果函数的导数信息可用，可以使用牛顿法来快速找到最小值
精确的求解每个一维子问题往往需要较大的计算量，更实际的方法是使用非精确的线搜索，以最小的代价获得使目标函数充分降低的步长
回溯线搜索法(Armijo准则)： $f(x^{k}+\alpha_{k}p_{k})-f(x^{k})\le\alpha_{k}\beta g_{k}^{T}p_{k}$

4.共轭梯度法

为避免锯齿的发生，取下一次的迭代搜索方向直接指向极小点，如果选定这样的搜索方向，对于二元二次函数只需进行两次直线搜索就可以求到极小点 $x^{1}=x^{0}+\alpha_{0}d^{0}$ $\frac{\partial f}{\partial \alpha}|_{x^{1}}=\nabla f(x^{1})^{T}d^{0}=0$ $x^{*}=x^{1}+\alpha_{1}d^{1}$
$d^{1}$ 应满足什么条件？
对于二次函数 f(x) 在 x* 处取得极小点的必要条件 $\nabla f(x^{*})=Ax^{*}-b=0$ $\nabla f(x^{*})=A(x^{1}+\alpha_{1}d^{1})-b=Ax^{1}-b+\alpha_{1}Ad^{1}$ $\nabla f(x^{*})=\nabla f(x^{1})+\alpha_{1}Ad^{1}=0$ 等式两边同乘 $d^{0})^{T}$ ，得 $d^{0})^{T}Ad^{1}=0$ $d^{0}$ 、 $d^{1}$ 是对 $A$ 的共轭方向，如果 $A = I$ ，则 $d^{0}$ 与 $d^{1}$ 正交
共轭梯度法：共轭方向由迭代点的负梯度构造出来 $f(x)=\frac{1}{2}x^{T}Ax-b^{T}x+c$ 从点 $x^{k}$ 出发，沿 $A$ 某一共轭方向 $d^{k}$ 作一维搜索，到达 $x^{k+1}$ $x^{k+1}=x^{k}+\alpha_{k}d^{k}$ 而在点 $x^{k}$ 、 $x^{k+1}$ 处的梯度分别为： $g^{k}=Ax^{k}-b\quad\quad g^{k+1}=Ax^{k+1}-b$ $d^{0}=-g^{0}=-Ax^{0}+b$ ，根据一维精确搜索， $g^{1})^{T}d^{0}=0$
令 $d^{1}=-g^{1}+\beta_{0}d^{0}$ ，则 $(d^{1})^{T}=-(g^{1})^{T}+\beta_{0}(d^{0})^{T}$ 选择 $\beta_{0}$ 使得 $d^{1})^{T}Ad^{0}=0$ ，上式两边同乘 $Ad^{0}$ ，可得： $(d^{1})^{T}Ad^{0}=-(g^{1})^{T}Ad^{0}+\beta_{0}(d^{0})^{T}Ad^{0}=0$ $\beta_{0}=\frac{(g^{1})^{T}Ad^{0}}{(d^{0})^{T}Ad^{0}}$
共轭梯度法流程
1.沿 $d^{k-1}$ 方向一维搜索 $x^{k}=x^{k-1}+\alpha_{k-1}d^{k-1}$ 2.构筑新的共轭方向 $d^{k}=-g^{k}+\beta_{k-1}d^{k-1}$ $\beta_{k-1}=\frac{(g^{k})^{T}Ad^{k-1}}{(d^{k-1})^{T}Ad^{k-1}}$ 简化，Fletcher-Reeves公式 $\beta_{k-1}=\frac{(g^{k})^{T}g^{k}}{(g^{k-1})^{T}g^{k-1}}$ 3.沿新的共轭方向 $d^{k-1}$ 一维搜索 $x^{k+1}=x^{k}+\alpha_{k}d^{k}$ Shewchuk, Jonathan Richard. “An introduction to the conjugate gradient method without the agonizing pain.” (1994): 7.

Python代码实现

def linear_conjugate_gradient(A,b,x0,tol=1e-5,max_iter)
# 初始化
    x = x0; r = b - A @ x; p = r
    rsold = r.T @ r
    for i in range(max_iter)
        Ap = A @ p
        alpha = rsold / (p.T @ Ap)
        x += alpha * p
        r -= alpha * Ap 
        rsnew = r.T @ r
        # 检查收敛性
        if np.sqrt(rsnew) <  tol:
             break  
        p = r + (rsnew / rsold) * p
        rsold = rsnew
   return x

5.多原子体系构型优化

在原子间相互作用势的作用下，通过改变粒子的几何构型，以能量最小为依据，得到体系的最优构型
以约化的LJ势为例，构建一个包含5个原子的二维初始构型，使用优化算法得到能量最低构型 $U(r_{ij})=4[(\frac{1}{r_{ij}})^{12}-(\frac{1}{r_{ij}})^{6}]$ $X=(x_{1},y_{1},x_{2},y_{2},x_{3},y_{3},x_{4},y_{4},x_{5},y_{5})$

定义Reduced LJ势函数

def LJ_reduced_pot(r):
    return 4*((1.0/r)**12-(1.0/r)**6)

定义导数

def LJ_reduced_force(r):
    return 24*(2*(1.0/r)**12-(1.0/r)**6)/r

生成随机构型

def genPts(N):
    np.random.seed(10)
    pts = np.zeros(N*2)
    for i in range(N*2):
        pts[i] = np.random.uniform(-2,2)
    return pts

定义目标函数

def E(x):
    N = int(len(x)/2)
    eng_pot = np.zeros(N)
    R = x.reshape(N,2)
    for i in range(N-1):
        for j in range(i+1,N):
            Rij = R[j]-R[i]
            r = np.linalg.norm(Rij)
            eng_pot[i] += 0.5*LJ_reduced_pot(r)
            eng_pot[j] += 0.5*LJ_reduced_pot(r)
    return np.sum(eng_pot)

2.积分算法

1.积分算法

N 个原子组成的分子体系，第 i 个原子的坐标、速度、动量及其作用力分别用 $r_{i}(t)$ ， $v_{i}(t)$ ， $p_{i}(t)$ ， $f_{i}(r,t)$ 表示，其初始值为 $r_{i}(0)$ ， $v_{i}(0)$ ， $p_{i}(0)$ ， $f_{i}(0)$ 。则第 i 个原子运动的牛顿方程为 $f_{i}=m_{i}\frac{d^{2}r_{i}}{dt^{2}}=m_{i}\ddot{r}_{i}\quad\quad (i=1,2,3,...,N)$
无约束条件下，3N 个自由度，3N 个分量微分方程
根据初始时刻系统状态求解 3N 个牛顿运动方程，通过积分算法可以获得系统所有原子的轨迹
第 i 个粒子的运动轨迹： $r_{i}(t)=\int v_{i}(t)dt$ 未知 $v_{i}(t)$ 解析形式，计算机程序算法： $r_{i}(t)=\sum_{0}^{t}v_{i}(t)\Delta t$ 使用计算机离散化求解数学问题近似解
分子动力学积分算法要求：
1.速度快，内存小
2.可以使用一个较长的时间步 $\Delta t$
3.可以接近和重复经典的路径
4.满足已知的能量和动量守恒以及时间反演性
5.形式简洁，易编程

2.前向Euler算法

一维谐振子运动方程的解析解
运动方程： $\frac{md^{2}x}{dt^{2}}=f=-kx$ 初始条件： $x(0)=0,\quad v(0)=v_{0}$ 解析解： $x(t)=\frac{v_{0}}{\omega}sin(\omega t)\quad\quad\omega=\sqrt{\frac{k}{m}}$ 体系总能： $E(t)=K(t)+V(t)=\frac{1}{2}mv_{0}^{2}=E(t=0)$
Euler算法基本思想
采用有限差分法中的一阶微分形式表示，即展开到 $\Delta t$ 的一阶泰勒展开公式： $f(t+\Delta t)=f(t)+\Delta t\frac{df}{dt}+O(\Delta t^{2})$ 将 $r (t)$ 、 $v (t)$ 泰勒展开后，得到前向Euler算法： $r(t+\Delta t)=r(t)+v(t)\Delta t+O(\Delta t^{2})$ $v(t+\Delta t)=v(t)+\frac{F(t)}{m}\Delta t+O(\Delta t^{2})$

Python代码实现

class HarmonicOscillator:
    def __init__(self,k,m):
        self.k = k
        self.m = m
     
    def acceleration(self,x):
        return -self.k * x / self.m

class ForwardEulerIntegrator:
    def step(self,system,x,v)
        # 执行一个时间步
        # x:当前位置  v:当前速度
        a = system.acceleration(x)
        x_new = x + self.dt * v
        v_new = v + self.dt * a
        return x_new,v_new

实际求解结果：
在分子动力学用向前Euler法计算体系能量时，经常会发生过热现象，即计算出的能量大于体系的真实能量，甚至可能导致体系的能量不收敛
真实位置： $r^{'}(t+\Delta t)=r(t)+\bar{v}\Delta t$ 预测位置： $r(t+\Delta t)=r(t)+v(t)\Delta t$ 预测位置大于真实位置(势能偏大)
真实速度： $v^{'}(t+\Delta t)=v(t)+\frac{\bar{F}}{m}\Delta t$ 预测速度： $v(t+\Delta t)=v(t)+\frac{F(t)}{m}\Delta t$ 预测速度大于真实速度(动能偏大)
体系过热现象的理论证明： $r(t+\Delta t)=r(t)+v(t)\Delta t+O(\Delta t^{2})$ $v(t+\Delta t)=v(t)+\frac{-dV(t)}{mdr}\Delta t+O(\Delta t^{2})$
忽略二阶小量，代入到体系总能量表达式中，计算相邻两步的能量差 $E(t+\Delta t)-E(t)=\frac{1}{2}m[v(t)+\frac{-dV(t)}{mdr}\Delta t]^{2}+[V(t)+\frac{dV(t)}{dt}\Delta t]-\frac{1}{2}mv(t)^{2}-V(t)$ $E(t+\Delta t)-E(t)=-v(t)\frac{dV(t)}{dr}\Delta t+\frac{dV(t)}{dt}\Delta t+\frac{1}{2m}[\frac{dV(t)}{dr}\Delta t]^{2}>0$ 总能增大，体系过热，前向Euler法无法保证能量守恒

3.Verlet算法

在 t 时刻求解 $t+\Delta t$ 时刻的坐标和力，分别向前和向后Taylor展开 $r(t+\Delta t)=r(t)+\dot{r}(t)\Delta t+\frac{\ddot{r}}{2}\Delta t^{2}+\frac{1}{6}\overset{...}{r}\Delta t^{3}+O(\Delta t^{4})$ $r(t-\Delta t)=r(t)-\dot{r}(t)\Delta t+\frac{\ddot{r}}{2}\Delta t^{2}-\frac{1}{6}\overset{...}{r}\Delta t^{3}+O(\Delta t^{4})$ 位置更新： $r(t+\Delta t)=2r(t)-r(t-\Delta t)+\frac{F(t)}{m}\Delta t^{2}+O(\Delta t^{4})$
在Verlet算法中，位置的更新不需要知道速度信息。然而想要得到动能信息就必须计算速度
速度更新： $v(t)=\frac{r(t+\Delta t)-r(t-\Delta t)}{2\Delta t}+O(\Delta t^{2})$
Verlet算法的启动 $r_{n+1}=2r_{n}-r_{n-1}+(\frac{F_{n}}{m})\Delta t^{2}+O(\Delta t^{4})...$
算法启动时， $n = 0$ ，已知 $r_{0}$ ，计算 $r_{1}$ 必须需要知道 $r_{-1}$ 才能开始更新。为避免计算 $r_{-1}$ ，通常通过Taylor展开来更新第一步： $r(t+\Delta t)=r(t)+v(t)\Delta t+\frac{a(t)}{2}\Delta t^{2}$ 经过这一步后，就可以使用计算得到的 $r_{0}$ 和 $r_{1}$ ，根据Verlet算法持续更新位置

Python代码实现

def step(self,system,x,v):
    if self.previous_x is None:
        self.previous_x = x-v0*self.dt+0.5*system.acceleration(x)* self.dt ** 2
    # 使用Verlet算法计算新位置
    a = system.acceleration(x)
    new_x = 2*x - self.previous_x + a * self.dt ** 2
    # 为当前位置计算速度
    if self.previous_x is not None:
         current_v = (new_x - self.previous_x) / (2*self.dt)
    else:
         current_v = 0
     # 为下一步更新previous_x
     self.previous_x = x
     return new_x,current_v

Verlet算法优点：1.位置精确度高，误差 $O(\Delta t^{4})$ 2.每次积分只计算一次力 3.时间可逆
Verlet算法缺点：1.速度有较大误差 2.轨迹与速度无关，无法与热浴耦联

4.Velocity Verlet算法

将 $r (t)$ 、 $v (t)$ 泰勒展开到二阶项： $r(t+\Delta t)=r(t)+v(t)\Delta t+\frac{1}{2}a(t)\Delta t^{2}+O(\Delta t^{3})$ $v(t+\Delta t)=v(t)+a(t)\Delta t+\frac{1}{2}\dot{a}(t)\Delta t^{2}+O(\Delta t^{3})$ 引入 $\dot{a}(t)=\frac{a(t+\Delta t)-a(t)}{\Delta t}$
位置更新： $r(t+\Delta t)=r(t)+v(t)\Delta t+\frac{1}{2}a(t)\Delta t^{2}$ 速度更新： $v(t+\Delta t)=v(t)+\frac{1}{2}[a(t)+a(t+\Delta t)]\Delta t$
引入半步速度
位置更新： $v(t+\frac{\Delta t}{2})=v(t)+\frac{\Delta t}{2}a(t)$ $r(t+\Delta t)=r(t)+v(t+\frac{\Delta t}{2})\Delta t$ 速度更新： $a(t+\Delta t)=\frac{F(r(t+\Delta t))}{m}$ $v(t+\Delta t)=v(t+\frac{\Delta t}{2})+\frac{\Delta t}{2}a(t+\Delta t)$

Python代码实现

def step(self,system,x,v):
    a = system.acceleration(x)
    x_new = x + self.dt * v + 0.5 * self.dt**2 * a
    a_new = system.acceleration(x_new)
    v_new = v + 0.5 * self.dt * (a + a_new)
    return x_new,v_new

5.Leapfrog算法

Leapfrog算法 $r(t+\Delta t)=r(t)+\Delta t\cdot v(t+\frac{\Delta t}{2})$ $v(t+\frac{\Delta t}{2})=v(t-\frac{\Delta t}{2})+\frac{\Delta t}{m}\cdot F(t)$
Leapfrog算法的自启动 $r(\Delta t)=r(0)+\Delta t\cdot v(\frac{\Delta t}{2})$ $v(\frac{\Delta t}{2})=v(-\frac{\Delta t}{2})+\frac{\Delta t}{m}\cdot F(0)$ $v(-\frac{\Delta t}{2})=v(0)-\frac{F(0)}{2m}\Delta t$
Leapfrog算法优点：节约内存，准确性和稳定性较高，速度、计算精度提高
Leapfrog算法缺点：速度和位置总是差半个时间步，不在同一时间定义，因此动能和势能也不能同时定义，无法直接计算总能

6.积分算法的稳定性

粒子演化涉及到位置对时间的积分，数值计算中必须考虑 $\Delta t$ 的选取
1. $\Delta t$ 的选取是一个折中的方案。 $\Delta t$ 太小，会导致有限长的时间内无法遍历所有的相空间； $\Delta t$ 太大，会导致积分算法的不稳定
2.分子动力学中常用的 $\Delta t$ 一般介于0.1-10 fs之间(10^-15 s)
3.对于液体的模拟，时间步长需要远小于两个液体分子的平均碰撞时间
4.对于固体或者分子振动，时间步长通常不大于最短的运动周期的1/10
牛顿运动方程的时间可逆性
$F=m\frac{d^{2}r}{dt^{2}}$ 时间反演 $t'=-t\quad\quad\frac{dr}{dt'}=-\frac{dr}{dt}\quad\quad F=m\frac{d^{2}r}{dt'^{2}}$ 牛顿运动方程形式不变
前向Euler算法的可逆性 $r(t+\Delta t)=r(t)+v(t)\Delta t$ $v(t+\Delta t)=v(t)+\frac{F(t)}{m}\Delta t$ 时间反演 $r'(t)=r(t+\Delta t)-v(t+\Delta t)\Delta t=r(t)-\frac{F(t)}{m}\Delta t^{2}$ $r'(t)\ne r(t)$ $v'(t)=v(t+\Delta t)-\frac{F(t+\Delta t)}{m}\Delta t=v(t)-\frac{F'(t)}{m}\Delta t^{2}$ $v'(t)\ne v(t)$
Verlet算法的可逆性
从 $r(t-\Delta t)$ ， $r (t)$ 推导 $r(t+\Delta t)$ $r(t+\Delta t)=2r(t)-r(t-\Delta t)+\frac{F(t)}{m}\Delta t^{2}$ 时间反演后，从 $r(t+\Delta t)$ ， $r (t)$ 推导 $r(t-\Delta t)$ $r'(t-\Delta t) = 2r(t)-r(t+\Delta t)+\frac{F(t)}{m}\Delta t^{2}=r(t-\Delta t)$
Velocity Verlet算法的可逆性
正向推导： $r(t+\Delta t)=r(t)+v(t)\Delta t+\frac{1}{2}a(t)\Delta t^{2}$ $v(t+\Delta t)=v(t)+\frac{1}{2}[a(t)+a(t+\Delta t)]\Delta t$ 逆向推导： $r'(t)=r(t+\Delta t)-v(t+\Delta t)\Delta t+\frac{1}{2}a(t+\Delta t)\Delta t^{2}=r(t)$
Leapfrog算法的可逆性
正向推导： $r(t+\Delta t)=r(t)+v(t+\frac{\Delta t}{2})\Delta t$ $v(t+\frac{\Delta t}{2})=v(t-\frac{\Delta t}{2})+\frac{F(t)}{m}\Delta t$ 逆向推导： $r'(t)=r(t+\Delta t)-v(t+\frac{\Delta t}{2})\Delta t=r(t)$ $v'(t-\frac{\Delta t}{2})=v(t+\frac{\Delta t}{2})-\frac{F(t)}{m}\Delta t=v(t-\frac{\Delta t}{2})$
前向Euler算法不满足时间可逆性，而Verlet、Velocity Verlet和Leapfrog算法满足时间可逆性
前向Euler算法不满足辛结构(Symplectic structures)，而Verlet、Velocity Verlet和Leapfrog算法满足辛结构