基于几何边界运算的分子动力学边界条件快速施加方法
阅读说明：本技术 基于几何边界运算的分子动力学边界条件快速施加方法 (Method for quickly applying molecular dynamics boundary conditions based on geometric boundary operation ) 是由 黄志刚 葛露明 黄亚军 何斌 覃圣超 史纪童 于 2021-01-29 设计创作，主要内容包括：本发明公开了一种基于几何边界运算的分子动力学边界条件快速施加方法,包括：根据分子动力学仿真模型,提取其几何边界,通过几何边界运算找出边界粒子,获取边界域内每个边界粒子的位置信息；为每个边界粒子赋予合理的速度和作用力,使得边界粒子的宏观物理量与微观状态相一致,从而得到完整的分子动力学边界条件。本发明给出了在复杂分子动力学系统上快速施加边界条件的一般方案；施加了平滑过渡的速度边界条件和应力边界条件,边界粒子的速度分布和势能分布,与系统局部的微观动力学状态一致,将边界附近的状态波动限制在较低的水平,解决施加边界条件任务繁重的问题,同时极大减少边界振荡。(The invention discloses a method for quickly applying a molecular dynamics boundary condition based on geometric boundary operation, which comprises the following steps: extracting a geometric boundary of the molecular dynamics simulation model according to the molecular dynamics simulation model, finding out boundary particles through geometric boundary operation, and acquiring position information of each boundary particle in a boundary domain; and (3) endowing each boundary particle with reasonable speed and force, so that the macroscopic physical quantity of the boundary particle is consistent with the microscopic state, and the complete molecular dynamics boundary condition is obtained. The invention provides a general scheme for rapidly applying boundary conditions on a complex molecular dynamics system; the velocity boundary condition and the stress boundary condition of smooth transition are applied, the velocity distribution and the potential energy distribution of boundary particles are consistent with the local micro dynamic state of the system, the state fluctuation near the boundary is limited to a lower level, the problem of heavy task of applying the boundary condition is solved, and the boundary oscillation is greatly reduced.)
技术领域
本发明涉及分子动力学计算仿真技术领域，具体涉及一种基于几何边界运算的分子动力学边界条件快速施加方法。
背景技术
随着计算机技术和数值仿真理论的发展，分子动力学(Molecular Dynamics,MD)方法越来越为人们所关注。分子动力学是通过数值求解分子体系力学运动方程得到体系的相轨迹，并统计体系结构特征与性质的一种数值仿真方法。广泛应用于生物计算、微流体和微纳工程问题的模拟。
分子动力学是一种基于粒子的计算方法，边界条件的施加是一项非常困难的任务。传统基于网格的计算方法，比如有限元法等，边界条件施加于几何边界上，可方便地通过图形交互等方式进行操作。然而分子动力学作为离散的粒子系统，必须要将边界作用施加于所有的边界粒子上。对于数目巨大的边界粒子的选择和操作，将是非常繁重的任务。特别是对于复杂形状的计算区域，边界粒子的快速操作，目前仍缺少有效的通用处理方法。
此外，分子动力学系统中，粒子的运动状态，包括位置、速度和势能等，会形成特定的分布。边界粒子的运动状态必须与局部的分布函数相适应，而不能简单地将边界作用均分给边界粒子，否则将引起系统状态的振荡，加大计算误差，甚至造成计算不收敛。
目前分子动力学计算规模越来越大，计算区域形状越来越复杂，对非稳态、非均匀系统的计算需求越来越多。如何用可视化方法快速建立计算模型，并施加相应的分子动力学边界条件，是现代分子动学仿真面临的重要挑战。主要需解决如下困难：
(1)边界粒子的快速选取和操作，对于大规模分子动力学计算而言，边界粒子数目非常巨大，可以多到百万级别，这将是一项非常繁重的任务。
(2)赋予每个边界粒子合适的速度和势能，分子动力学粒子之间的作用非常强，如果粒子运动状态与分布函数不相应，有可能发生剧烈的振荡，局部条件和目标状态之间会出现瞬时偏差，密度涨落等边界效应容易发生。
发明内容
本发明的目的是提供一种基于几何边界运算的分子动力学边界条件快速施加方法，用以解决现有技术施加边界条件任务繁重的问题。
为了实现上述任务，本发明采用以下技术方案：
一种基于几何边界运算的分子动力学边界条件快速施加方法，包括：
根据分子动力学仿真模型，提取其几何边界，通过几何边界运算找出边界粒子，获取边界域内每个边界粒子的位置信息；
为每个边界粒子赋予合理的速度和作用力，使得边界粒子的宏观物理量与微观状态相一致，从而得到完整的分子动力学边界条件。
进一步地，所述通过几何边界运算找出边界粒子，获取边界域内每个边界粒子的位置信息，包括：
依据边界的几何特征构建几何边界的实体模型；
根据所述的实体模型，获取边界抽象的几何信息，根据边界模型的几何信息对实体模型进行离散操作，得到离散粒子；根据边界的几何特征采用相交法判别离散粒子与边界的位置关系，从而确定边界粒子和模拟域粒子，最终得到完整的边界粒子的位置信息。
进一步地，所述依据边界的几何特征构建几何边界的实体模型的过程中，对于较为规则的几何边界，采用偏置法得到，对于较为复杂的边界，使用布尔运算获得。
进一步地，所述为每个边界粒子赋予合理的速度，包括：
对于所有边界粒子，基于速度修正的方式，采用松弛因子对边界粒子的速度进行修正，使得边界粒子的速度向着模拟域的方向缓慢变化，速度修正式为：
Vi'＝Vi+h(xi)(u-Vi) 式1
式1中，i表示第i个边界粒子；Vi'为边界粒子修正后的速度；Vi为边界粒子的初始速度；h(xi)为松弛因子，其值在0到h0区间内变化；xi为粒子在沿着模拟域法线方向从0到A变化的某一坐标值；u为边界粒子在分子动力学仿真系统目标状态时的期望速度，A为边界域沿法向量n方向的宽度。
进一步地，基于边界的几何特征，确定边界域朝着模拟域方向的法向量n，依据边界粒子的位置信息，沿着法向量方向，通过所述速度修正式的数值计算方式，依次为所有边界粒子赋予速度。
进一步地，所述为每个边界粒子赋予合理的作用力，包括：
实际模拟中，要在边界域内施加边界力，假定在边界域外存在一些外部粒子，他们对边界粒子产生作用势，外部粒子在截断半径RC范围内对边界粒子存在有效作用势，这些作用势的总和构成了有效边界作用力；假定任意边界粒子在边界域外周围包裹有一层薄壳，首先确定单层薄壳对边界域内任意边界粒子的作用力，继而确定外部粒子施加在边界域内任一边界粒子的总作用力。
进一步地，所述确定单层薄壳对边界域内任意边界粒子的作用力，包括：
壳内的粒子数可以表示为：ρg(r)S(r)dr，壳内的粒子对边界粒子存在作用力，沿着模拟域法线方向，任意一层薄壳内粒子对边界域内任一粒子的作用力为：其中，ρ为分子动力学仿真系统中的粒子密度；
g(r)为粒子径向分布函数；S(r)为壳的面积；T(r)为薄壳面积在边界域外层上的投影；V(r)为单个壳内粒子对边界域内任一粒子的作用势；dr为壳的厚度。
进一步地，所述外部粒子施加在边界域内任一边界粒子的总作用力表示为：
上式中，Φ为应力张量；P为维里压强；T(r,x)为边界粒子在坐标为x时，薄壳面积在边界域外层上的投影。
进一步地，对于边界粒子i，在其截断半径Rc范围内，边界域外任一外部粒子j对其作用力为Fij，通过假定任意边界粒子在边界域外周围包裹有一层薄壳，利用式2计算得到边界域外截断半径范围内所有粒子施加在边界粒子i上的有效作用力Fm，最终实现为所有边界粒子施加合适的边界作用力。
与现有技术相比，本发明具有以下技术特点：
1.本发明方法可以迅速找出边界粒子并施加合理边界条件，解决施加边界条件任务繁重的问题，同时极大减少边界振荡。
2.本发明给出了在复杂分子动力学系统上快速施加边界条件的一般方案；施加了平滑过渡的速度边界条件和应力边界条件，边界粒子的速度分布和势能分布，与系统局部的微观动力学状态一致，将边界附近的状态波动限制在较低的水平。
附图说明
图1为本发明方法的流程示意图；
图2本发明获取边界粒子原理示意图；
图3本发明赋予边界粒子速度原理图；
图4本发明施加边界粒子作用力原理图。
具体实施方式
参见图1，本发明提供一种基于几何边界运算的分子动力学边界条件快速施加方法，该方法对于分子动力学仿真模型，首先通过几何边界运算，快速找出边界粒子，获取边界域内每个边界粒子的位置信息；其次为每个边界粒子赋予合理的速度和作用力，使得边界粒子的宏观物理量与微观状态相一致，从而得到完整的分子动力学边界条件。
具体地，本发明的分子动力学边界条件快速施加方法包括三个过程，具体说明如下：
1.根据给定的几何边界获取边界粒子
该步骤是施加边界条件时最基本的步骤，主要是依据分子动力学仿真模型中的几何边界，离散得到分子动力学边界粒子信息，构建边界模型。具体步骤如下：
1.1构建几何边界的实体模型
根据分子动力学仿真模型，提取其几何边界，依据边界的几何特征构建几何边界的实体模型；在这过程中，对于较为规则的几何边界，可以采用偏置法得到，对于较为复杂的边界，则可以使用布尔运算获得；其中，对于实体模型，其二维模型中的直线或三维模型中的平面可认为是较为规则的几何边界。
如图2的(a)所示，在图2给出的示例中，对于较为规则的分子动力学模型边界的几何特征，可以直接采用偏置法，将模型沿着X方向偏置b，沿着Y方向偏置a，得到边界实体模型；如图2的(b)所示，当仿真模型边界不规则或为复杂曲面时，边界模型可以通过两个或多个模型进行布尔运算得到。
1.2建立几何边界粒子模型
完成步骤1.1后，根据几何边界的实体模型，获取边界抽象的几何信息，根据边界模型的几何信息对实体模型进行离散操作；得到离散粒子后，根据边界的几何特征采用相交法判别离散粒子与边界的位置关系，从而确定边界粒子和模拟域粒子，如图2所示，最终得到完整的边界粒子的位置信息，并构建边界粒子模型。
2.为边界粒子赋予合理的速度
施加在分子动力学仿真系统的边界条件主要分为两部分，速度和作用力。该步骤中为边界粒子赋予速度，其速度需要满足狄利克雷边界条件，同时使得分子动力学仿真系统粒子的微观状态与宏观物理量相符。为避免边界域内和模拟域内的粒子速度产生较大偏差，引发系统震荡，必须让边界域内的粒子速度缓慢变化，即速度“软过渡”；该步骤的具体实现方式如下：
对于所有边界粒子，基于速度修正的方式，采用一个“钟形函数”形式的松弛因子，对边界粒子的速度进行修正，使得边界粒子的速度向着模拟域的方向缓慢变化，速度修正式为：
Vi'＝Vi+h(xi)(u-Vi) 式1
式1中，i表示第i个边界粒子；Vi'为边界粒子修正后的速度；Vi为边界粒子的初始速度；h(xi)为松弛因子，其值在0到h0区间内变化，h0的具体值取决于具体的分子动力学仿真模型，例如一般情况下h0的取值可以为1；xi为粒子在沿着模拟域法线方向从0到A(边界域沿法向量n方向的宽度)变化的某一坐标值；u为边界粒子在分子动力学仿真系统目标状态时的期望速度。
边界粒子由边界开始，向着模拟域法向量n方向的速度值从u开始缓慢变化，引入的松弛因子函数应满足h'(0)＝0，h'(A)＝0，而且松弛因子在边界区域内紧支，从而将边界作用限制在边界区域内，松弛因子函数的引入取决于具体的分子动力学仿真模型。
如图3所示，基于边界的几何特征，确定边界域朝着模拟域方向的法向量n，依据步骤1.2得到的边界粒子的位置信息，沿着法向量方向，通过式1的数值计算方式，依次为所有边界粒子赋予速度。
引入松弛因子和采用速度修正式后，可以快速得到合理的边界粒子速度，使其符合狄利克雷边界条件，并满足仿真系统的微观状态，减小系统震荡。
3.为边界粒子施加合适的作用力
为保证边界粒子在微观状态下与仿真系统内应力及微观作用势相一致，引入边界作用力，确保系统内应力的连续性，并与局部粒子分布函数相适应。根据局部粒子分布函数，引入虚拟的“外部粒子”，重建外部边界对内部边界粒子的作用，基于维里定理推导施加于边界粒子上的额外作用力。具体实现方式如下：
实际模拟中，要在边界域内施加边界力，可以假定在边界域外存在一些外部粒子，他们对边界粒子产生作用势，外部粒子在截断半径RC范围内对边界粒子存在有效作用势，这些作用势的总和构成了有效边界作用力。由于分子动力学仿真模型中所有粒子间的作用势始终限制在截断半径RC范围内，因此，可以假定任意边界粒子在边界域外周围包裹有一层薄壳，薄壳的半径r在x(边界域内粒子沿法线方向到边界的距离)到RC范围内变化。
壳内的粒子数可以表示为：ρg(r)S(r)dr，壳内的粒子对边界粒子存在作用力，沿着模拟域法线方向，任意一层薄壳内粒子对边界域内任一粒子的作用力为：其中，ρ为分子动力学仿真系统中的粒子密度；
g(r)为粒子径向分布函数；S(r)为壳的面积；T(r)为薄壳面积在边界域外层上的投影；V(r)为单个壳内粒子对边界域内任一粒子的作用势；dr为壳的厚度。
单层薄壳对边界域内任一边界粒子的作用力确定后，可以得到外部粒子施加在边界域内任一边界粒子的总作用力为：
上式中，Φ为应力张量；P为维里压强；T(r,x)为边界粒子在坐标为x时，薄壳面积在边界域外层上的投影。
如图4所示，对于边界粒子i，在其截断半径Rc范围内，边界域外任一外部粒子j对其作用力为Fij，采用前面所述的薄壳法，通过式2积分计算得到边界域外截断半径范围内所有粒子施加在边界粒子i上的有效作用力Fm，最终实现为所有边界粒子施加合适的边界作用力，使得边界粒子的宏观物理量与微观状态相一致，从而得到完整的分子动力学边界条件。
以上实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的精神和范围，均应包含在本申请的保护范围之内。