CN103279645A - 基于gpu并行运算的碳纳米管分子动力学仿真方法 - Google Patents

Abstract

本发明公开了一种基于GPU并行运算的碳纳米管分子动力学仿真方法，包括以下步骤：1)由数据库文件构造由碳纳米管和C60分子组成的碳纳米管系统模型；2)C60分子受力，在平衡位置上来回振动，当C60分子处于平衡状态时，使碳纳米管管壁上的碳粒子和C60分子的位置和速度发生变化；3)在CUDA平台上将碳纳米管分割成多层计算单元，采用CPU遍历计算单元进行轮转计算得到可并行运算队列，调度GPU的流处理单元进行并行运算和处理；4)重复执行步骤3)，最终输出C60分子随着模拟时间的运行轨迹，并描绘碳纳米管内随着能量传递的温度变化曲线图，完成仿真过程。本发明极大地提高分子动力学仿真的运算效率。

Description

基于GPU并行运算的碳纳米管分子动力学仿真方法

技术领域

本发明涉及一种碳纳米管分析动力学仿真方法，尤其是一种基于GPU并行运算的碳纳米管分子动力学仿真方法，属于仿真技术领域。

背景技术

在碳纳米管分子动力学仿真过程中，由于其运算中涉及超大规模的分子量及大量的复杂的浮点运算，使得碳纳米管分子动力学模拟需要极强的计算能力。因为碳纳米管由大量分子构成，其模拟仿真中涉及大量的浮点运算，随着分子动力学仿真规模的扩大，仿真结果的可视化也非常重要。对于这些碳纳米管分子动力学模拟系统来说，提高仿真算法的浮点运算能力和并行计算能力具有重大作用。

在新兴的纳米工程领域，建立在连续介质基础上的宏观机理很难解释纳米工程中出现的一些特殊现象，因而分子动力学方法成为重要的研究手段之一。碳纳米管分子动力学仿真，主要通过牛顿力学模拟碳纳米管体系的运动，在其不同状态时的数据集合中抽取数据样本，从而计算碳纳米管体系的构型积分，并以构型积分的结果为基础进一步计算碳纳米管体系的热力学量和其他宏观性质，从而高效仿真碳纳米管系统。但碳纳米管分子动力学仿真算法受计算能力两大局限：其一是极小空间内包括极大量粒子，计算量极大；其二是为保证数值模拟稳定性，典型的分子动力学时间步长为飞秒(fs)级，对计算机性能要求高。

受限于计算机性能的发展，为扩大仿真算法规模，国内外许多专家学者对分子动力学仿真算法做了大量研究。由最初扩大单机规模的串行算法，发展到现在通过将计算任务分配给多个CPU而扩大模拟规模的并行算法，可以模拟的规模由最初的几千个原子提高到上百万、上千万甚至到上亿个原子的规模。虽然CPU的计算能力已经非常强大，但对于具有实际意义的超大规模的模拟系统来说，其运行速度还远不能满足实际需求，而且对于其批量处理碳，单CPU的串行实现效率相对来说也比较低。因此，大规模的碳纳米管分子动力学仿真算法需要进一步突破。

发明内容

本发明的目的是为了解决上述现有技术的缺陷，提供一种基于GPU并行运算的碳纳米管分子动力学仿真方法，该方法能大幅度提高碳纳米管的模拟规模，极大地提高分子动力学仿真的运算效率。

本发明的目的可以通过采取如下技术方案达到：

基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于包括以下步骤：

1)由数据库文件构造碳纳米管系统模型，该模型由碳纳米管和设置在碳纳米管内部的C60分子组成，对模型进行初始化，读取模型中所有粒子的速度和位置坐标；

2)C60分子受到碳纳米管管壁上碳粒子产生的力，在平衡位置上来回振动，当C60分子处于平衡状态时，在碳纳米管的前、后两端同时分别设置不同温度的条件进行模拟，并设定模拟时间，使碳纳米管管壁上碳粒子和C60分子的速度和位置不断发生变化，进而使C60分子的受力、运动状态和能量不断发生变化；

3)在CUDA平台上对碳纳米管进行分割，分割成多层相互独立的计算单元，采用CPU遍历计算单元进行轮转计算得到可并行运算队列，再调度GPU的流处理单元，采用Verlet算法进行并行运算和处理，重复执行该步骤直至模拟时间结束；

4)采用数据库文件记录相关数据，输出C60分子随着模拟时间的运行轨迹，并描绘碳纳米管内随着能量传递的温度变化曲线图，完成仿真过程。

作为一种优选方案，所述碳纳米管为层状中空结构，由多个六边形碳环结构组成；所述C60分子有一个或多个。

作为一种优选方案，所述步骤2)中，采用Gaussian热浴模拟法在碳纳米管的前、后两端分别设置不同温度的热浴条件，并设定热浴时间作为模拟时间。

作为一种优选方案，所述步骤3)中，采用CPU遍历计算单元进行轮转计算具体如下：

a)若所有粒子均已作为中心计算粒子，跳到步骤h)；

b)按次序找到第一个非竞争粒子并没有按中心计算的粒子，加入可并行运算队列；

c)标记该粒子的所有近邻为本并行队列一度不可并行粒子，若近邻已经是二度不可并行粒子，则提升为一度不可计算；

d)标记该粒子所有次邻为本并行队列二度不可并行粒子，若次邻为一度不可并行粒子，则不修改其度数；

e)标记本粒子为已按中心计算粒子；

f)是否已遍历到粒子队尾？若是继续执行，若否跳到步骤a)；

g)下一个可并行队列开始，返回步骤b)；

h)结束。

作为一种优选方案，所述步骤3)中，采用Verlet算法进行并行运算和处理具体如下：

a)按碳纳米管系统模型所有粒子的位置，计算近邻粒子以及次邻粒子间的键关系和角度关系；

b)调度GPU的流处理单元，并行计算被分割在不同计算单元里的粒子，并积累计算碳纳米管管壁上每个粒子与其近邻粒子的相互作用力；

c)根据C60分子所在区域，计算C60分子内每个粒子与其所在碳纳米管区域中的每个粒子的相互作用力；

d)根据粒子所受到的力及其速度，更新粒子位置，再执行步骤b)和c)；

e)根据粒子所受到的力，计算粒子本次速度及碳纳米管前、后两端热浴的热流值；

f)若达到循环频数，则计算结束；否则，在CPU端间隔保存粒子的数据，返回步骤d)。

作为一种优选方案，所述步骤3)中，Verlet算法具体设计如下：

a)将x(t+Δt)和x(t-Δt)进行泰勒展开如下式所示：

$\stackrel{\→}{x}(t+\mathrm{\Δt})=\stackrel{\→}{x}\left(t\right)+\stackrel{\→}{v}\left(t\right)\mathrm{\Δt}+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}+\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left(\mathrm{\Δ}{t}^4\right)---\left(1\right)$

$\stackrel{\→}{x}(t-\mathrm{\Δt})=\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{v}\left(t\right)\mathrm{\Δt}+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}-\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left(\mathrm{\Δ}{t}^4\right)---\left(2\right)$

其中，x(t+Δt)表示为前一时刻的位置，x(t-Δt)表示为后一时刻的位置；

b)将式(1)和式(2)相加，得到位置表达式如下：

$\stackrel{\→}{x}(t+\mathrm{\Δt})=2\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{x}(t-\mathrm{\Δt})+\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2+O\left(\mathrm{\Δ}{t}^4\right)---\left(3\right)$

在已知碳纳米管系统模型粒子当前时刻t的位置和加速度，以及前一时刻t-Δ的位置的情况下，推算出下一时刻t+Δt的位置；

c)将式(1)和式(2)相减，再两边同时除以2Δt，获得速度的表达式如下：

$\stackrel{\→}{v}\left(t\right)=\frac{\stackrel{\→}{x}(t+\mathrm{\Δt})-\stackrel{\→}{x}(t-\mathrm{\Δt})}{2\mathrm{\Δt}}+O\left({\mathrm{\Δt}}^2\right)---\left(4\right)$

在已知碳纳米管系统模型粒子前一时刻t-Δt的位置和后一时刻t+Δt的位置的情况下，推算出当前时刻t的速度v(t)；

d)由式(1)～(4)，在已知粒子t-2Δt时刻的位置、t-Δt时刻的位置和t-Δt时刻的加速度的情况下，启动Verlet算法进行积分：根据t-2Δt时刻的位置、t-Δt时刻的位置和t-Δt时刻的加速度，将t＝t-Δt代入式(3)，获得当前时刻t的位置；根据当前时刻t的位置，基于一定的势函数更新当前时刻t的加速度；同时，根据当前时刻t的位置和t-2Δt时刻的位置，将t＝t-Δt代入式(4)，更新t-Δt时刻的速度；即得到粒子当前时刻t的位置、t-Δt时刻的速度和t时刻的加速度，重复执行该步骤。

本发明相对于现有技术具有如下的有益效果：

1、本发明采用数据库文件构造碳纳米球C60模型和碳纳米管结构模型，从而对碳纳米管分子动力学进行仿真，不仅方便计算，也便于构造多种不同模型进行比较，可移植性较高，并利用算法可以模拟足球状结构的C60分子在碳纳米管里的运行轨迹。

2、本发明采用由NVIDIA公司发布的基于GPU的并行运算平台CUDA进行了图像处理，CUDA平台是一个基于C语言的基础架构，具备有大量高性能计算指令和良好的编程接口，可以极大地提高碳纳米管分子动力学仿真方法的效率。

3、本发明通过CUDA平台可以对大型的碳纳米管系统模型进行分割，分割成多层大小适宜、相互独立且GPU能承受的计算单元，使得各算法在CUDA中并行运行，极大提高了碳纳米管分子动力学仿真方法的速度。

附图说明

图1为本发明的C60分子示意图。

图2为本发明的碳纳米管示意图。

图3为本发明构造的碳纳米管系统模型的结构示意图。

图4为本发明在CPU与GPU中碳纳米管分子动力学仿真的数据流框图。

图5为本发明在CUDA环境下进行碳纳米管分子动力学仿真的流程示意图。

图6为本发明基于GPU并行运算与传统的串行运算的运行时间对比图。

具体实施方式

实施例1：

如图1～图5所示，本实施例基于GPU并行运算的碳纳米管分子动力学仿真方法如下：

1)在数据库文件中读取碳纳米管系统的参数条件，如初始温度、粒子数、密度时间等，构造碳纳米管系统模型MOD1及参数模型MOD2，模型MOD2可以设置各种参数在模型MOD1上使用，模型MOD1如图1～图3所示，由碳纳米管和设置在碳纳米管内部的一个足球状C60分子(由60个碳原子构成的分子，形似足球，又名足球烯)组成，所述碳纳米管为层状中空结构，管身是准圆管结构，由多个六边形碳环结构单元(由碳原子构成)组成，其直径一般在一到几十个纳米之间，长度则远远大于其直径，对该模型进行初始化，读取模型中所有粒子的速度和位置坐标；

2)C60分子受到碳纳米管管壁上碳粒子产生的力，在平衡位置上来回振动，当C60分子处于平衡状态时，采用Gaussian热浴法(约束温度调节方法)，在碳纳米管两端同时分别设置不同温度的热浴条件，并设定热浴时间，使碳纳米管管壁上碳粒子和C60分子因平衡状态被打破而发生位置和速度的变化，而这些变化是相互影响的，具有传递性的，伴随着这些粒子位置和速度的变化，能量在长管状的碳纳米管系统内发生传递，由于C60分子位于碳纳米管内，其受力会随着碳纳米管组成粒子的位置改变而改变，因此伴随这能量传递的过程，C60分子的运动状态和能量也会发生改变；

采用Gaussian热浴法的基本原理是在运动方程中加入摩擦力f_i，并将其与粒子速度v_i联系起来，其受力公式为：

f_i＝ma_i+ξmv_i

当平衡态时，系统温度不变，因此有dEk/dt＝0

即有： $\underset{i}{\mathrm{\Σ}}{v}_i{a}_i=0$

由此可得： $\mathrm{\ξ}=\frac{\underset{i}{\mathrm{\Σ}}{f}_i{v}_i}{m\underset{i}{\mathrm{\Σ}}{v}_i^2}$

3)在CUDA平台上对碳纳米管进行分割，分割成多层大小适宜、相互独立的计算单元，分割原则是在避免过多重复计算的前提上提高并行度，即分割的计算单元必须合理粗细化，因为计算单元过大，则并行不明显，计算单元过小，会导致太多不必要的重复的计算；

如图4所示，采用CPU遍历计算单元进行轮转计算得到n个可并行运算队列，如下：

a)若所有粒子均已作为中心计算粒子，跳到步骤h)；

b)按次序找到第一个非竞争粒子并没有按中心计算的粒子，加入可并行运算队列；

c)标记该粒子的所有近邻为本并行队列一度不可并行粒子，若近邻已经是二度不可并行粒子，则提升为一度不可计算；

d)标记该粒子所有次邻为本并行队列二度不可并行粒子，若次邻为一度不可并行粒子，则不修改其度数；

e)标记本粒子为已按中心计算粒子；

f)是否已遍历到粒子队尾？若是继续执行，若否跳到步骤a)；

g)下一个可并行队列开始，返回步骤b)；

h)结束；

如图4和图5所示，调度GPU的流处理单元，采用Verlet算法进行并行运算和处理，如下：

a)按碳纳米管系统模型所有粒子的位置，计算近邻粒子以及次邻粒子间的键关系和角度关系；

b)调度GPU的流处理单元，并行计算被分割在不同计算单元里的粒子，并积累计算碳纳米管管壁上每个粒子与其近邻粒子的相互作用力；

c)根据C60分子所在区域，计算C60分子内每个粒子与其所在碳纳米管区域中的每个粒子的相互作用力(范德华力)；

d)根据粒子所受到的力及其速度，更新粒子位置，再执行步骤b)和c)；

e)根据粒子所受到的力，计算粒子本次速度及碳纳米管前、后两端热浴的热流值；

f)若达到循环频数，则计算结束；否则，在CPU端间隔保存粒子的数据，返回步骤d)。

4)重复执行步骤3)，直至热浴时间结束，采用数据库文件记录相关数据，输出C60分子随着模拟时间的运行轨迹，并描绘碳纳米管内随着能量传递的温度变化曲线图，完成仿真过程。

所述步骤3)中，由于粒子间的作用力间距影响较大，必须设定一个距离临界值，以判断其位置关系，当粒子间的间距小于距离临界值，为近邻粒子，当粒子间的间距大于距离临界值，为次邻粒子，则认为其相互间的分子作用力可以忽略不计。

所述步骤3)中，Verlet算法具体设计如下：

a)将x(t+Δt)和x(t-Δt)进行泰勒展开如下式所示：

$\stackrel{\→}{x}(t+\mathrm{\Δt})=\stackrel{\→}{x}\left(t\right)+\stackrel{\→}{v}\left(t\right)\mathrm{\Δt}+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}+\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left(\mathrm{\Δ}{t}^4\right)---\left(1\right)$

$\stackrel{\→}{x}(t-\mathrm{\Δt})=\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{v}\left(t\right)\mathrm{\Δt}+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}-\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left(\mathrm{\Δ}{t}^4\right)---\left(2\right)$

其中，x(t+Δt)表示为前一时刻的位置，x(t-Δt)表示为后一时刻的位置；

b)将式(1)和式(2)相加，得到位置表达式如下：

$\stackrel{\→}{x}(t+\mathrm{\Δt})=2\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{x}(t-\mathrm{\Δt})+\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2+O\left(\mathrm{\Δ}{t}^4\right)---\left(3\right)$

在已知碳纳米管系统模型粒子当前时刻t的位置和加速度，以及前一时刻t-Δ的位置的情况下，推算出下一时刻t+Δt的位置；

c)将式(1)和式(2)相减，再两边同时除以2Δt，获得速度的表达式如下：

$\stackrel{\→}{v}\left(t\right)=\frac{\stackrel{\→}{x}(t+\mathrm{\Δt})-\stackrel{\→}{x}(t-\mathrm{\Δt})}{2\mathrm{\Δt}}+O\left({\mathrm{\Δt}}^2\right)---\left(4\right)$

在已知碳纳米管系统模型粒子前一时刻t-Δt的位置和后一时刻t+Δt的位置的情况下，推算出当前时刻t的速度v(t)；

d)由式(1)～(4)，在已知粒子t-2Δt时刻的位置、t-Δt时刻的位置和t-Δt时刻的加速度的情况下，启动Verlet算法进行积分：根据t-2Δt时刻的位置、t-Δt时刻的位置和t-Δt时刻的加速度，将t＝t-Δt代入式(3)，获得当前时刻t的位置；根据当前时刻t的位置，基于一定的势函数更新当前时刻t的加速度，从而得到其受力；同时，根据当前时刻t的位置和t-2Δt时刻的位置，将t＝t-Δt代入式(4)，更新t-Δt时刻的速度；即得到粒子当前时刻t的位置、t-Δt时刻的速度和t时刻的加速度(受力)，重复执行该步骤。

如图6所示，可以看到本发明采用的GPU并行运算，其执行时间要远远小于传统的串行运算，因此，其仿真效率较之传统的串行运算有很大的提高。

实施例2：

本实施例的主要特点是：所述步骤1)中，设置在碳纳米管内部的C60分子可以为多个，其余同实施例1。

以上所述，仅为本发明优选的实施例，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明所公开的范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都属于本发明的保护范围。

Claims (6)

1.基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于包括以下步骤：

1)由数据库文件构造碳纳米管系统模型，该模型由碳纳米管和设置在碳纳米管内部的C60分子组成，对模型进行初始化，读取模型中所有粒子的速度和位置坐标；

2)C60分子受到碳纳米管管壁上碳粒子产生的力，在平衡位置上来回振动，当C60分子处于平衡状态时，在碳纳米管的前、后两端同时分别设置不同温度的条件进行模拟，并设定模拟时间，使碳纳米管管壁上碳粒子和C60分子的速度和位置不断发生变化，进而使C60分子的受力、运动状态和能量不断发生变化；

3)在CUDA平台上对碳纳米管进行分割，分割成多层相互独立的计算单元，采用CPU遍历计算单元进行轮转计算得到可并行运算队列，再调度GPU的流处理单元，采用Verlet算法进行并行运算和处理；

4)重复执行步骤3)，直至模拟时间结束，采用数据库文件记录相关数据，输出C60分子随着模拟时间的运行轨迹，并描绘碳纳米管内随着能量传递的温度变化曲线图，完成仿真过程。

2.根据权利要求1所述的基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于：所述碳纳米管为层状中空结构，由多个六边形碳环结构组成；所述C60分子有一个或多个。

3.根据权利要求1或2所述的基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于：所述步骤2)中，采用Gaussian热浴模拟法在碳纳米管的前、后两端分别设置不同温度的热浴条件，并设定热浴时间作为模拟时间。

4.根据权利要求3所述的基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于：所述步骤3)中，采用CPU遍历计算单元进行轮转计算具体如下：

a)若所有粒子均已作为中心计算粒子，跳到步骤h)；

b)按次序找到第一个非竞争粒子并没有按中心计算的粒子，加入可并行运算队列；

c)标记该粒子的所有近邻为本并行队列一度不可并行粒子，若近邻已经是二度不可并行粒子，则提升为一度不可计算；

d)标记该粒子所有次邻为本并行队列二度不可并行粒子，若次邻为一度不可并行粒子，则不修改其度数；

e)标记本粒子为已按中心计算粒子；

f)是否已遍历到粒子队尾？若是继续执行，若否跳到步骤a)；

g)下一个可并行队列开始，返回步骤b)；

h)结束。

5.根据权利要求4所述的基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于：所述步骤3)中，采用Verlet算法进行并行运算和处理具体如下：

a)按碳纳米管系统模型所有粒子的位置，计算近邻粒子以及次邻粒子间的键关系和角度关系；

b)调度GPU的流处理单元，并行计算被分割在不同计算单元里的粒子，并积累计算碳纳米管管壁上内每个粒子与其近邻粒子的相互作用力；

c)根据C60分子所在区域，计算C60分子内每个粒子与其所在碳纳米管区域中的每个粒子的相互作用力；

d)根据粒子所受到的力及其速度，更新粒子位置，再执行步骤b)和c)；

e)根据粒子所受到的力，计算粒子本次速度及碳纳米管前、后两端热浴的热流值；

f)若达到循环频数，则计算结束；否则，在CPU端间隔保存粒子的数据，返回步骤d)。

6.根据权利要求5所述的基于GPU并行运算的碳纳米管分子动力学仿真方法，其特征在于：所述步骤3)中，Verlet算法具体设计如下：

a)将x(t+Δt)和x(t-Δt)进行泰勒展开如下式所示：

$\stackrel{\→}{x}(t+\mathrm{\Δt})=\stackrel{\→}{x}\left(t\right)+\stackrel{\→}{v}\left(t\right)\mathrm{\Δt}+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}+\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left(\mathrm{\Δ}{t}^4\right)---\left(1\right)$

$\stackrel{\→}{x}(t-\mathrm{\Δt})=\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{v}\left(t\right)\mathrm{\Δt}+\frac{\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2}{2}-\frac{\stackrel{\→}{b}\left(t\right){\mathrm{\Δt}}^3}{2}+O\left(\mathrm{\Δ}{t}^4\right)---\left(2\right)$

其中，x(t+Δt)表示为前一时刻的位置，x(t-Δt)表示为后一时刻的位置；

b)将式(1)和式(2)相加，得到位置表达式如下：

$\stackrel{\→}{x}(t+\mathrm{\Δt})=2\stackrel{\→}{x}\left(t\right)-\stackrel{\→}{x}(t-\mathrm{\Δt})+\stackrel{\→}{a}\left(t\right){\mathrm{\Δt}}^2+O\left(\mathrm{\Δ}{t}^4\right)---\left(3\right)$

在已知碳纳米管系统模型粒子当前时刻t的位置和加速度，以及前一时刻t-Δ的位置的情况下，推算出下一时刻t+Δt的位置；

c)将式(1)和式(2)相减，再两边同时除以2Δt，获得速度的表达式如下：

$\stackrel{\→}{v}\left(t\right)=\frac{\stackrel{\→}{x}(t+\mathrm{\Δt})-\stackrel{\→}{x}(t-\mathrm{\Δt})}{2\mathrm{\Δt}}+O\left({\mathrm{\Δt}}^2\right)---\left(4\right)$

在已知碳纳米管系统模型粒子前一时刻t-Δt的位置和后一时刻t+Δt的位置的情况下，推算出当前时刻t的速度v(t)；

d)由式(1)～(4)，在已知粒子t-2Δt时刻的位置、t-Δt时刻的位置和t-Δt时刻的加速度的情况下，启动Verlet算法进行积分：根据t-2Δt时刻的位置、t-Δt时刻的位置和t-Δt时刻的加速度，将t＝t-Δt代入式(3)，获得当前时刻t的位置；根据当前时刻t的位置，基于一定的势函数更新当前时刻t的加速度；同时，根据当前时刻t的位置和t-2Δt时刻的位置，将t＝t-Δt代入式(4)，更新t-Δt时刻的速度；即得到粒子当前时刻t的位置、t-Δt时刻的速度和t时刻的加速度，重复执行该步骤。

Images