运用机器学习方法设计原子链人工边界条件1)

张 慊 乔 丹 唐少强,2)

*(北京大学工学院力学与工程科学系，北京100871)

†(北京大学数学科学学院概率统计系，北京100871)

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

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

1 原子链的人工边界条件

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

其中uj表示第j个原子的位移。

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

的线性叠加，其中

为色散关系。ξ∈[-π,0]为右行波，ξ∈[0,π]为左行波。

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

N表示边界条件用到的原子数，我们取N= 6,c0= 1。其中cj及bj(j= 0,1,···,6) 即为待定的边界条件系数。

定义残差函数

可以发现如果波数为ξ的左行波能够完全满足边界条件，则Δ(ξ) = 0。而如果不能完全满足边界条件，Δ(ξ) /= 0，就会产生反射。可以知道反射系数

2 前馈神经网络

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

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

图1 单层神经网络

神经网络学习从输入

到输出的映射，即学习函数hW,b(x)满足

其中θ为参数，代表W,b。

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

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

图2 多层神经网络

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

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

梯度下降算法即α是学习率，一般设置为一个较小的正数。其中h(X;θ)表示当前参数下神经网络的输出，ˆh(X;θ)表示训练的目标值。

对多层的前馈神经网络，使用链式法则求损失函数J关于各个参数的梯度。记是第k-1 层第j个神经元连接到第k层的第i个神经元的权值表示第k层第j个神经元的偏置表示第k层第j个神经元的输入表示第k层第j个神经元的输出，即为激活函数。

第k层第j个神经元的误差(实际值与预测值的误差)定义为应用链式法则，有

其中

对上式两边微分

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

3 数值算例

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

以N= 6 为例，取L= 120,Lc= 1200，令jlb=-90处作为左边界点。取初值为

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

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

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

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

表1 N =6 时机器学习得到的人工边界条件的系数

表2 N =6 时文献[3]提出的人工边界条件的系数

图3 N =6 时机器学习方法与参考解的数值结果比较，实线代表参考解，圆圈代表采用人工边界条件的数值解。由于解的对称性，只画出了通过右侧边界的波形

图4 N =6 时Δ 关于ξ 的关系，其中点线与ξ 轴重合，表示精确的边界条件，点划线代表文献[3]提出的人工边界条件的残差，实线代表机器学习得到的人工边界条件的残差

图5 N =6 时R 关于ξ 的关系，其中点线与ξ 轴重合，表示精确的边界条件，点划线代表文献[3]提出的人工边界条件的反射系数，实线代表机器学习得到的人工边界条件的反射系数

4 结论

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