分子动力学模拟中壁面边界的模拟方法

摘要

本发明公开提出一种分子动力学模拟中壁面边界的模拟方法，该方法综合了全反射法和随机反射法优点，能够真实反应出粒子碰壁后的运动状态，也称为半反射法，其实现步骤如下：1)定义分子动力学模拟中粒子的数据结构；2)确定粒子在通道中的当前位置A(x0，y0，z0，u0，v0，w0)；3)根据粒子当前位置及粒子受到其它粒子的作用力信息，确定粒子的运动方向，并根据粒子当前位置及受力信息判断粒子下一步是否会运动到壁面边界以外位置；4)若粒子运动到壁面边界以外位置，经过壁面反射到达下一步反射位置B(x，y，z，u，v，w)，当x方向存在壁面，则其中y＝y’，z＝z’，x＝2*壁面x方向的坐标-x’，，，粒子的速度方向随机。本发明在模拟中可以体现出壁面恒温和壁面不光滑的物理条件；可以更加真实地反应出粒子反射后的位置；也可应用在不规则的纳米通道中，简化速度矢量的计算方法。

说明书

技术领域

本发明涉及一种分子动力学模拟中的方法，尤其是模拟粒子在壁面边界条件下的运动轨迹的方法。

背景技术

现在越来越多的科学家用分子动力学方法模拟平行平板中流体和纳米通道内流体流动的动力学特性。随着研究的深入，科学家们开始对非平直纳米通道进行研究和采用分子动力学方法对此通道进行模拟，如：Xi-Jun Fan等人采用非平衡分子动力学方法模拟了简单流体在周期性喷管中的流动特性；为研究粗糙度对微通道内体流动及其边界滑移性质的影响，曹柄阳等人研究了氩气在锯齿状铂通道内的流动；J.Castillo-Tejas等人分别对线性聚合链和牛顿流体在4∶1∶4收缩-扩展通道中的流动过程进行了模拟。在使用分子动力学方法模拟两平板间和纳米通道内流体的流动特性时，以前的文献中对虚拟热壁面(Virtualthermal wall)情况下壁面边界条件的模拟有两种方法，一种是Allen，M.P.and D.J.Tildesley，在文献Computer Simulation of Liquids.中提出的全反射法，全反射法将壁面看成光滑镜面，粒子在壁面上发生镜面反射，采用全反射法模拟粒子在壁面边界条件下的运动轨迹参见图1(a)，当粒子从A点运动到B’时，由于虚拟壁面的存在，粒子不会在B’的位置出现，粒子会在壁面发生全反射后运动到B点，B与B’点关于壁面对称，与壁面平行坐标方向的速度矢量保持不变，与壁面垂直坐标方向的速度矢量大小保持不变，方向反向；虽然全反射法可以真实地反应出粒子在发生发射后的位置，但是全反射法不能模拟出粒子在定壁温和壁面不光滑的条件下的运动轨迹，并且如果在非平直纳米通道中，若采用全反射法来模拟粒子的壁面边界条件下的轨迹，则在计算速度矢量方面也存在一定的难度；另一种是Raraport，D.C.在文献The art of molecular dynamics simulation.中提出了随机反射法，随机反射法与物理上的漫发射对应，采用随机反射法模拟粒子在壁面边界条件下的运动轨迹参见图1(b)，粒子反射到B’向壁面做垂线的垂点B，速度大小由壁面温度确定，方向随机。虽然随机反射法可以模拟出壁面不光滑以及定壁温的物理条件，但是随机反射法由于粒子反射后的位置在壁面垂直方向的单一性，因此在模拟粒子运动轨迹的过程中存在一定的不真实性，不能真实反映出粒子反射后的位置。

因此，如何真实模拟粒子在壁面边界条件下的运动状态，成为了本领域内的研究人员急需解决的技术问题。

发明内容

针对现有技术存在的上述不足，本发明提供一种在分子动力学模拟中，当粒子处于壁面边界条件下时，可综合全反射法和随机反射法优点，真实体现粒子在壁面边界条件下的运动状态的壁面边界模拟方法。

为了实现上述发明目的，本发明采取以下技术方案：分子动力学模拟中壁面边界的模拟方法，其特征在于，所述模拟方法包含以下步骤：

A)定义分子动力学模拟中粒子的数据结构，粒子在三维模型中由六个信息来表示粒子在分子动力学模拟通道中的运动状态，所述六个信息分别是x方向坐标、y方向坐标、z方向坐标、x方向速度、y方向速度和z方向速度；其中x方向坐标为壁面横坐标，y方向坐标为壁面纵坐标，z方向为与XY平面垂直的方向；

B)确定粒子在通道中的当前位置A(x₀，y₀，z₀，u₀，v₀，w₀)，其中x₀表示粒子在x方向上的当前位置坐标，y₀表示粒子在y方向上的当前位置坐标，z₀表示粒子在z方向上的当前位置坐标，u₀表示粒子在x方向上的当前位置速度，v₀表示粒子在y方向上的当前位置速度，w₀表示粒子在z方向上的当前位置速度；

C)根据粒子当前位置及粒子受到其它粒子的作用力信息，确定粒子的运动方向，并根据粒子当前位置及受力信息判断粒子下一步是否会运动到壁面边界以外位置；如果粒子没有运动到壁面边界以外位置，则粒子会沿着确定的运动方向运动，到达下一步指定位置B’(x’，y’，z’，u’，v’，w’)，其中x’为粒子在x方向上的指定位置坐标，y’表示粒子在y方向上的指定位置坐标，z’表示粒子在z方向上的指定位置坐标，u’表示粒子在x方向上的指定位置速度，v’表示粒子在y方向上的指定位置速度，w’表示粒子在z方向上的指定位置速度；

D)若粒子运动到壁面边界以外位置，粒子将不能到达下一步指定的位置B’，粒子根据壁面的温度，经过壁面反射到达下一步反射位置B(x，y，z，u，v，w)，当x方向存在壁面，则其中y＝y’，z＝z’，x＝2*壁面x方向的坐标-x’，粒子的速度方向随机；当y方向存在壁面，则其中x＝x’，z＝z’，y＝2*壁面y方向的坐标-y’，粒子的速度方向随机；当z方向存在壁面，则其中x＝x’，y＝y’，z＝2*壁面z方向的坐标-z’，粒子的速度方向随机。

上述公式中，a为-1～1之间的随机数(该随机数的取得采用文献The art of moleculardynamics simulation(D.C.Rapaport，Second Edition，Cambridge University Press，2004)中的程序进行计算，在此不作详细描述)，k_B为波尔兹曼常数，T为当前温度，m为质量，x表示粒子在x方向上的反射位置坐标，y表示粒子在y方向上的反射位置坐标，z表示粒子在z方向上的反射位置坐标，u为粒子在x方向上的反射位置速度，v为粒子y方向上的反射位置速度，w为粒子在z方向上的反射位置速度。

本发明分子动力学模拟中壁面边界的模拟方法，其主要原理是在镜面反射基础上考虑了壁面粗糙度对粒子速度的影响，因此，本发明方法又可以称为半反射法。

通常情况下，统计系统的宏观量，可以通过粒子的运动状态，即在本发明方法通过粒子分别在x，y，z方向上的坐标和粒子分别在x，y，z方向的速度来反映粒子的运动状态，当粒子运动状态发生变化时，x，y，z方向上的坐标和粒子分别在x，y，z方向的速度都会发生变化。当采用半反射方法对粒子处于壁面边界条件下的运动轨迹进行模拟时，其既考虑粒子的随机性，即粒子在遇到壁面进行反射后，其反射位置速度(u，v，w)会根据壁面的温度有关，而其反射位置x，y，z方向坐标(x，y，z)则与粒子未碰到壁面的预定的下一步指定位置的x，y，z方向坐标(x’，y’，z’)有关。因此本方法在模拟粒子下一步位置时，不仅考虑到壁面的温度，也考虑到了粒子经过壁面发生反射后，其位置必定与粒子的未经过反射的位置有关系。

本发明产生的有益效果：

(1)本发明方法不仅可以在分子动力学模拟中体现恒壁温条件和壁面不光滑的物理条件，而且也可以在分子动力学模拟中真实反映粒子反射后的位置。其解决了随机反射法中不能真实反应粒子反射后位置的问题。

(2)本发明方法还可以在不规则纳米通道中，简化速度矢量计算方法，从而便于计算。解决了全反射边界法不能模拟恒壁温条件以及壁面不光滑条件下粒子的真实运动状态的问题和在不规则纳米通道中，全反射法不能简化速度矢量计算方法的问题。

(3)本发明方法可以用于模拟粒子在各种不同通道中运动到虚拟壁面边界外之后的运动情况，从而为研究粒子在不同通道中运动状态，提供了技术支撑。

附图说明

图1是分子动力学模拟中的粒子运动轨迹示意图；其中(a)为采用全反射法模拟粒子在壁面边界条件下的运动轨迹示意图；其中(b)为采用随机反射法模拟粒子在壁面边界条件下的运动轨迹示意图；

图2是分子动力学模拟中的粒子在两种边界条件下反射示意图：其中(a)半反射粒子轨迹运动示意图；(b)随机反射粒子轨迹运动示意图；

图3是分子动力学模拟中的变截面通道示意图；其中(a)为三维变截面通道示意图；(b)为三维变截面通道正视示意图；

图4是分子动力学模拟中粒子数密度分布随到壁面距离的比较图；

图5是本发明方法模拟变截面纳米通道中势能随时间变化曲线图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细地说明。

本文中提出一种综合了全反射法和随机反射法优点的新的壁面边界条件处理方法——半反射法。当粒子A点运动到B’时，如图1(a)，由于壁面的存在，粒子经过壁面反射后到达B点，坐标位置与全反射法粒子反射后相同，但是速度大小由壁面温度确定，方向随机，即速度矢量的确定方法与随机反射法相同。在模拟非平直微纳米通道时，比如非45°倾斜通道，如果壁面采用全反射边界条件处理法，粒子反射后，速度矢量的确定有一定的难度。随机反射法相对于全反射法隐含了两个条件：(1)壁面恒温；(2)壁面不光滑，与真实情况更加接近。但是随机反射法在模拟纳米通道时有一定的不真实性。如图2所示，x方向为流动方向，当粒子运动到区域A时，如果采用随机反射法，见图2(b)，无论粒子B处在A中的任意位置经过坐标调节后都出现在B’，不能真实反应出粒子运动的坐标位置。但如果采用半反射边界条件模拟方法，如图2(a)，粒子B出现在区域A中的任意位置，采用半反射法经过坐标调节后可以出现在B对应任意B’点，从而与粒子的真实运动位置更为接近。

实施例1

本实施例通过在分子动力学模拟过程中，分别采用随机反射法、全反射法和半反射三种壁面边界条件模拟方法来模拟平行平板中的氩流体流动过程，从而进一步说明本发明方法与现有的随机反射法和全反射法的不同之处。具体实施如下：

采用氩原子模拟NVT系中，用速度修正法控制系统温度，原子间相互作用采用12-6LJ势能模型，运动方程采用Verlet数值积分法，时间步长为1fs，模拟时间为100000fs，模拟尺寸为6nm×6nm×1.5nm，粒子数为648，密度为1g/cm³.在本实施例中，将坐标原点定义为0点，x和y方向采用周期性边界条件，当z方向(Z方向是指三维通道中与xy平面垂直的方向，粒子只可能超出z方向的壁面边界，因为x和y方向采用周期性边界条件，所以在这个实施例中都是无限延伸的，没有壁面)为0时，xy平面为下壁面，当z为1.5nm时，xy平面为上壁面，下面分别采用三种不同的壁面边界条件模拟方法来获取粒子的运动状态：

随机反射法：x＝x’；y＝y’；z＝0；

全反射法：x＝x’；y＝y’；z＝2*壁面z方向的坐标-z’；u＝u’；v＝v’；w＝-w’

半反射法：x＝x’；y＝y’；z＝2*壁面z方向的坐标-z’；

上述公式中，a为-1～1之间的随机数，k_B为波尔兹曼常数，T为当前温度，m为质量，x表示粒子在x方向坐标上的反射位置坐标，y表示粒子在y方向坐标上的反射位置坐标，z表示粒子在z方向坐标上的反射位置坐标，u为粒子在x方向坐标上的反射位置速度，v为粒子在y方向坐标上的反射位置速度，w为粒子在z方向坐标上的反射位置速度。本文中的“*”为乘号，“-”为减号。

从采用随机反射法、全反射法和半反射法的模拟在特定条件下平行平板中的氩流体流动中碰壁的情况来看，当采用随机反射法模拟时，参见图1(b)粒子A(x₀，y₀，z₀，u₀，v₀，w₀)在运动过程中超出模拟区域边界时，粒子不会沿着既定线路到达下一步指定位置B’(x’，y’，z’，u’，v’，w’)，而会沿着粒子反射到B’向壁面做垂线的垂点B，x＝x’；y＝y’；z＝壁面z方向的坐标；速度大小由壁面温度确定，方向随机，；当采全反射法模拟式，参见图1(a)，粒子A(x₀，y₀，z₀，u₀，v₀，w₀)在运动过程中超出模拟区域边界时，粒子不会沿着既定线路到达下一步指定位置B’(x’，y’，z’，u’，v’，w’)，粒子会在壁面发生全反射后运动到B(x，y，z，u，v，w)点，B(x，y，z，u，v，w)与B’(x’，y’，z’，u’，v’，w’)点关于壁面对称，与壁面平行坐标方向的速度矢量保持不变，与壁面垂直坐标方向的速度大小保持不变，速度方向反向；当采用半反射法模拟时，粒子A(x₀，y₀，z₀，u₀，v₀，w₀)在运动过程中超出模拟区域边界时，粒子不会沿着既定线路到达下一步指定位置B’(x’，y’，z’，u’，v’，w’)，而会在虚拟的壁面边界上发生反射，到达反射位置B(x，y，z，u，v，w)，其反射后位置B既与粒子未发生反射时的指定位置B’(x’，y’，z’，u’，v’，w’)有关，也与粒子发生反射时的壁面的温度有关，即x＝x’；y＝y’；z＝2*壁面z方向的坐标-z’；粒子反射后的速度方向随机，从而得到粒子反射后的位置。由于粒子反射后的速度与三个方向的速度分量u，v和w有关，因此，粒子反射后的速度大小为而速度的方向为三个速度分量的速度方向的合成，而速度分量是有正有负的，所以合成的速度方向也是各个角度都有的，从而粒子反射后的速度方向是随机的。本说明书中，壁面x方向的坐标、壁面y方向的坐标和壁面z方向的坐标是指壁面在三维笛卡尔坐标系的坐标，由于在实际应用过程中三维坐标系的原点可以为0点，也可以不为0点，因此，壁面本身在三维坐标系中的坐标与三维坐标系原点的设置有关。

2)模拟结束后，统计了三种壁面边界条件模拟方法情况下粒子数密度分布的情况，如附图4所示。

比较结果说明：采用半反射边界条件模拟法和采用全反射边界条件模拟法模拟得到的结果100％吻合，具有很强的实用性。从而表明本发明方法在模拟粒子在壁面边界条件下的运动状态时，不仅可以体现恒壁温条件和壁面不光滑的物理条件，而且也可以在分子动力学模拟中真实反映粒子反射后的位置。

半反射边界法与随机反射边界法和全反射边界法异同见表1和表2：

表1半反射边界法与随机反射边界法和全反射边界法特性比较表

表2半反射边界法与随机反射边界法和全反射边界法相同点

  随机反射边界条件  全反射边界条件  反射后粒子位置  不同  相同  反射后粒子速度的大小  相同  不同  反射后粒子速度的方向  相同  不同

从表1和表2中可以看出，半反射边界条件模拟方法既结合了随机反射边界条件模拟中考虑壁面温度和光滑程度的优点，也结合全反射边界条件模拟中粒子实际位置的变化，因此，其与随机反射边界条件模拟和全反射边界条件模拟相比，更能反映粒子的在碰壁后的运动状态的变化。

实施例2：采用半反射边界条件模拟方法，采用分子动力学模拟变截面通道中聚乙烯流体的流动过程。具体步骤如下：

1)采用本发明模拟在变截面通道结构中进行，聚乙烯分子的流动过程参见图3。用速度修正法控制系统温度，运动方程采用Verlet数值积分法，时间步长为1fs，模拟时间为320000fs，密度为0.2588g/cm³。

2)模拟结束后，统计了不同时间点系统的势能，如附图5所示。势能即位能，是系统宏观状态的表征物理量之一。当采用常用的随机方法时，发现当粒子在下一步运动到同时超过xy和xz平面的通道外面任意位置时，粒子发生发射后都只能在固定的16个点出现，这与实际情况不符；而当采用常用的全反射方法时，当粒子下一步运动到超过通道斜面边界时，很难确定粒子反射后的速度矢量，此时，采用本发明方法(半反射方法)是一种真实有效的方式，因为采用本发明法，粒子发生发射后，可以真实的出现在通道中的任意位置，并且易于计算出速度矢量。

模拟结果说明：本发明可以应用到分子动力学模拟实例的过程中，可以真实的模拟粒子在碰到虚拟壁面后发生反射的运动状态，是一种真实有效的数值模拟计算方法。

本发明方法可以用于模拟粒子多种通道中遇到壁面后的运动状态，在模拟粒子在通道中的真实运动状态时，既需要考虑粒子的正常运动状态，即粒子未碰到虚拟壁面的运动状态，(该正常状态的模拟方法已为成熟的现有技术，在本专利也不做详细阐述)，又需要考虑粒子非正常运动状态，即粒子在碰到虚拟壁面后运动状态，在对该运动状态的模拟中，现有技术提出了全反射法和随机反射法，但是由于该两种方法自身固有的缺陷，因而不能真实的有效的模拟出粒子碰壁后的有效状态，从而本申请人提出“半反射法”，该方法克服了全反射法和随机反射法的缺陷，能够真实有效的模拟出粒子碰壁后的运动状态。该方法的提出为有效统计粒子在每一刻的位置，速度大小，并进一步统计出系统的温度，动能、势能等宏观物理参数，从而为研究不同粒子在不同通道中的运动状态技术支持。

本发明方法可以模拟不同粒子在不同通道中碰壁后的运动状态，比如：氩流体流动过程，水分子流动过程以及各种通道中聚乙烯流体流动过程等其它情况，特别需要说明的是如果模拟通道为圆形通道时，坐标转化成柱坐标或者极坐标后，本发明方法仍然适用。