分子动力学计算常用于精细分析材料的性质,往往采用周期边界条件。对于多尺度计算和其他一些力学研究,周期边界条件不适用,需要设计人工边界条件^[1]。本文将针对一维线性原子链,运用机器学习设计人工边界条件。该方法需要的先验知识很少,且可以快速产生一系列边界条件。

传统的获得人工边界条件的方法有时间历史积分^[2]与在给定波数附近进行泰勒展开的方法^[3]等。本文提出的方法主要是受后者的启发。

考虑一维无穷长原子链,每个原子与其邻近原子线性相互作用。归一化控制方程为

其中$u_{j}$表示第$j$个原子的位移。

方程(1)的解可表示为模态

的线性叠加,其中

为色散关系。$\xi \in [-\pi ,0]$为右行波,$\xi \in [0,\pi ]$为左行波。

如果数值截断模拟有限长度内的原子链,需要在数值边界上设计人工边界条件。以左边界$u_{0}$为例,我们设计以下形式的线性边界条件

$N$表示边界条件用到的原子数,我们取$N=6$, $c_{0}=1$。其中$c_{j}$及$b_{j}$ ($j=0,1,\cdots,6$) 即为待定的边界条件系数。

定义残差函数

可以发现如果波数为$\xi $的左行波能够完全满足边界条件,则$\varDelta \left( \xi \right)=0$。而如果不能完全满足边界条件,$\varDelta \left( \xi \right) \ne 0$,就会产生反射。可以知道反射系数$R(\xi )=- {\varDelta (\xi )}/{\varDelta (-\xi )}$。

前馈神经网络是神经网络的一种,由一个输入层、若干个隐藏层和一个输出层组成。每一层由一组权重和一个激活函数组成,上一层的输出经过权重的线性组合和激活函数的映射之后输出给下一层隐藏层,直到最终输出。

最简单的前馈神经网络结构如图1所示(这是没有隐藏层的情况)：

神经网络学习从输入

到输出的映射,即学习函数$h_{W,b}({x})$满足

其中$\theta $为参数,代表${W},b$。

为对应每个$x_{i}$的权重,$b$为常数项。$f$称为激活函数,在得到输入加权的结果后,用激活函数将该值映射到最终的输出。常见的激活函数有：线性函数(相当于不使用激活函数),ReLU,sigmoid(用于二分类问题),softmax(用于多分类问题),Leaky ReLU,tanh等。

两层的前馈神经网络如图2所示。

更一般地,输出也可以是一个向量。

训练神经网络一般使用反向传播算法,本质是一种梯度下降算法,即通过求损失函数$J$关于每个参数的偏导数,然后用梯度下降更新参数。损失函数$J$可以定义成均方误差

梯度下降算法即

$\alpha $是学习率,一般设置为一个较小的正数。其中$h(X;\theta)$表示当前参数下神经网络的输出,$\hat{h}(X;\theta)$表示训练的目标值。

对多层的前馈神经网络,使用链式法则求损失函数$J$关于各个参数的梯度。记$W_{ij}^{k}$是第$k-1$层第$j$个神经元连接到第$k$层的第$i$个神经元的权值,$b_{j}^{k}$表示第$k$层第$j$个神经元的偏置,$z_{j}^{k}$表示第$k$层第$j$个神经元的输入,$a_{j}^{k}$表示第$k$层第$j$个神经元的输出,即$a_{j}^{k}= f(z_{j}^{k})$,$f$为激活函数。

第$k$层第$j$个神经元的误差(实际值与预测值的误差)定义为

应用链式法则,有

其中

对上式两边微分

代入$\delta_{j}^{k}$的表达式,即得

再对每个参数求梯度即可得到需要更新的梯度大小。

整体梯度下降以及随机梯度下降的收敛速度通常较慢,可使用Adam算法^[4]进行优化。

对于给定初始条件$u=f_{i}(x)$,设$f_{i}(x)$支集为闭区间$E=[-{L}/{2}, {L}/{2}]$。在区间$E_{0}=[-{L_{\rm c}}/{2}, {L_{\rm c}}/{2}]$ 上用二阶Verlet算法对时间离散进行计算,将总长度取得充分长($L_{\rm c}>5L$)并计算足够步数(使波前不超出边界),选取计算区间中一个点作为左边界点,记录对应的$u_{0},\cdots,u_{N-1}\dot{u}_{0},\cdots,\dot{u}_{N-1}$,即得参考数值解,用于训练一个单层的前馈神经网络。该神经网络以$\dot{u}_{0}$为输出,以$u_{0},\cdots,u_{N-1}\dot{u}_{1},\cdots,\dot{u}_{N-1}$为输入。由线性约束假设(4),函数应该是线性的,因此选取激活函数为线性激活函数,即$f(x)=x$。因为区间上原子位移全部为零对应于原子链静止的情形,边界原子的速度也一定是零,因此零输入一定对应零输出,所以偏置$b$取为零。于是训练得到的神经元的权值就是边界条件(4)中的系数。

以$N=6$为例,取$L=120$, $L_{\rm c}=1200$,令$j_{lb}=-90$处作为左边界点。取初值为

$A$为任意正常数,$M$为正整数,这里令$A=50$,$M=17$。离散时间步长$\Delta t=0.005$,则单位时间有200步,计算400个单位时间。将$u_{j}^{n}$, $j=j_{lb}+1,j_{lb}+2,\cdots ,j_{lb}+6$及$\dot{u}_{j}^{n}$, $j=j_{lb}+0,j_{lb}+1,\cdots,j_{lb}+6$作为训练集的输入,$u_{j_{lb}}^{n}$作为训练集的输出。训练200轮之后得到系数$c_{j}b_{j}$, $j=0,1,2,\cdots,6$(如表1所示)。 与文献^[3]提出的边界条件(表2所示)对比,两者差别很大。为了验证这组系数的正确性,在区间[-50,50]上取初值

其中$h_{0}=\dfrac{50}{1-{e}^{-6.25}}$,离散时间步长$\Delta t=0.005$,计算符合人工边界条件的数值解并在[-500,500]上计算参考数值解,结果如图3所示。

然后分析机器学习方法得到的边界条件的残差。计算$\varDelta $关于$\xi$的关系,结果如图4所示。

从图4中可见,文献^[3]提出的人工边界条件的残差在$\xi\to \pi $,即短波处发散到无穷。而机器学习得到的人工边界条件虽然在$\xi \to 0$,即长波处残差较大,但是短波处残差是有界的。

为进一步考察机器学习得到的人工边界条件在能量意义下的精确度,计算边界条件对应的反射系数$R$关于$\xi$的关系,结果如图5所示。这里有反射系数大于1的波段,但前述算例仍是稳定的;此外,还可以通过在训练数据中增加更多这个波段的模态来改善。

对于线性原子链的人工边界条件这一问题,本文在基于泰勒展开的匹配边界条件方法基础上,提出了基于前馈神经网络的机器学习算法。该方法不需要任何对方程(1)的细致分析就能得到一组较好的边界条件。虽然该方法在一维线性原子链模型上的应用仍然存在残差和反射系数较大的问题,但是应用该方法需要的先验知识很少,编程方便且容易实现,因此具有较高的实用价值,可以进一步推广到其他人工边界条件的设计。

(责任编辑: 胡漫)