第13章分子动力学模拟的CUDA程序的开发

最新推荐文章于 2024-04-29 19:05:34 发布

布林组-？

最新推荐文章于 2024-04-29 19:05:34 发布

阅读量384

点赞数

分类专栏： CUDA从入门到实践文章标签：人工智能算法数据结构 c++

本文链接：https://blog.csdn.net/qq_52299902/article/details/132347752

版权

CUDA从入门到实践专栏收录该内容

13 篇文章 4 订阅

订阅专栏

在前面的章节中，已经用一些简短的程序较为系统的介绍了CUDA编程的基础知识。为了能改有更加深入地讨论CUDA程序的优化策略，将讨论一个出具规模的CUDA程序的开发与优化。该程序可实现一个简单地分子动力学模拟。分子动力学模拟的算法涉及本科水平的物理学(经典力学和经典统计力学)。

13.1 分子动力学模拟的基本算法和C++实现

本章的程序有多个版本，具有一个C++版本盒及格CUDA版本，每个版本都有几百行的代码。

13.1.1 程序的整体结构

在C++程序中没有使用类，整个程序的风格时C语言风格的，该程序由若干头文件(扩展名为.cuh)和源文件(扩展名为.cu)构成，大多是两两配对的。头文件中放置结构体定义、常量定义、宏定义及函数声明，源文件中放置有函数定义。每个源文件是一个编译单元，在连接之前会被编译成为一个目标文件。整个程序的计算流程主要体现在main.cu的main函数中。程序文件的依赖关系如下：
主函数(main.cu)
--内存的分配和释放(memory.cuh和memory.cu)
--模拟系统的初始化(initialize.cuh和initialize.cu)
--构件邻居列表(neighbor.cuh和neighbor.cu)
*实施最小镜像约定(mic.cuh)
--运动方程的积分(integrate.cuh和intergate.cu)
*求力和势能(force.cuh和force.cu)
$实施最小镜像约定

13.1.2 分子动力学模拟的基本流程

分子动力学模拟是一种数值计算方法。在这种方法中，对一个具有一定初始条件和边界条件且具有相互作用的多粒子系统的运动方程进行数值积分，得到系统在相空间中的一条离散的轨迹，并利用统计学的方法从这条相轨迹中提取出有用的物理结果。
相空间是经典力学和统计力学中常用的概念，他是坐标空间的推广。例如：对于一个三维空间中的例子(质点)，我们可以用一个具有三个分量x,y,z的矢量表示其在坐标空间的一个位置(一个点)。而该例子的相空间是一个六维的空间，其中的一个点(相空间点)可以用3个坐标分量(x,y,z)和三个动量分量(px,py,ppz)合起来描述。如果系统中有N个粒子，那么就要用3N个坐标和3N个动量来描述全部粒子的微观状态，这6N个量可以看做一个6N维空间的分量。这样构成的6N维空间就是这N个粒子的相空间，一般来说每个粒子的坐标和动量是随时间变化的，那么N个粒子的6N维相空间中的点也是随时间变化的。随时间变化的相空间点就构成了相轨迹。在分子动力学模拟中，我们只能得到离散的相轨迹。因为我们对粒子的运动方程进行数值积分时用了有限大小的时间间隔，而不是无限小的时间间隔。一条很长(很多时刻的)相轨迹包含系统的很多微观状态，而通过统计力学的方法就可以根据这些微观状态计算出有用的宏观物理量，如系统的温度和压强。
根据上述定义，可以设想一个典型的、简单地分子动力学的模拟有如下的计算流程：
(1)设置系统的初始条件，具体包括。
1)初始化各个粒子的位置矢量。
2)初始化各个粒子的速度矢量
(2)根据系统中的粒子所满足的相互作用规律，由牛顿定律确定所有粒子的运动方程(二阶常微分方程组)，并对运动方程进行数值积分，即不断的更新每个粒子的坐标和速度。最终得到一系列离散时刻系统在相空间的位置，即一条离散的相轨迹。
(3)用统计力学的方法分析相轨迹所蕴含的物理规律。

13.1.3 初始条件

初始化是指给定一个初始的相空间点，包含各个粒子初始的坐标和速度。在分子动力学模拟中，需要对3N个二阶常微分方程进行数值积分，因为每一个二阶常微分方程的求解都需要两个初始条件，所以我们需要确定6个初始条件：3N个初始坐标分量和3N个初始速度分量。
坐标的初始化指的是为系统中每个粒子选定一个初始的位置坐标。分子动力学模拟中如何初始化位置主要取决于所要模拟的体系。例如，如果要模拟固态氮的，就需要让各个氮原子的位置按照面心立方结构排列。如果要模拟液态或者气态，那么初始坐标的选取就比较随意了。重要的是，在构造的初始结构中，任何两个粒子之间的距离都不能太小，因为距离太小会使得这些粒子受到非常大的力，以致于让后面的数值积分变得非常不稳定。坐标的初始化也常被称为建模，往往需要用到一些专业知识，如固体物理学中的例子。在本例的模拟过程中，只需要考虑固态氮的模拟，所以不需要太多的专业知识，前面提到的面心立方体结构是一种晶体结构，读者若不熟悉也不用担心。
我们知道，在任何经典热力学系统中，平衡时各个粒子的速度满足麦克斯韦分布。然而，作为初始条件，我们并不一定要求粒子的速度满足麦克斯韦分布。最贱的速度初始化方法是产生3N个在某个区间均匀分布的随机速度分量，再通过如下两个基本条件对速度分量进行修正。
第一个条件是让系统的总动量为零。也就是说，不希望系统的质心在模拟过程中跑动。分子间作用力是所谓的内力，不会改变整体动量，即系统的整体动量是守恒的。只要初始的整体动量为零，在分子动力学模拟的时间演化过程中整体动量将保持为零，如果整体动量明显偏离零，则说明模拟出现了问题。这正是程序判断是否有误的标准之一。
第二个条件是系统的总动能应该与所选定的初始温度对应。我们知道，再经典统计力学中，能量均分定理成立，即粒子的哈密顿量中每一个具有平方形式的能联想的统计平均值都等于KBT/2.其中KB是波尔兹曼常数，T是系统的绝对温度。所以将在质心的动量取为零之后就可以对每个粒子的速度进行一个标度变换，使得系统的初始温度与所设定的温度一致。我们在大学物理或普通物理中学过，对于拥有N个粒子的系统，其温度T可以用每个粒子的质量mi和速度vi表达为(本章用黑体表示三维空间中的矢量)：
$\frac{3N}{2}k_{B}T=\sum_{i=1}^{N}\frac{1}{2}m_{i}V_{i}^{2}$

其中矢量Vi的平方Vi2代表它的各个分量(vix,viy,viz)的平方和：
$V_{i}^{2}=V_{ix}^{2}+V_{iy}^{2}+V_{iz}^{2}$
假设我们设置的目标温度为T0，那么需要对各个粒子的速度做怎样的操作才会让系统的温度从T变为T0呢，很简单只需要做以下变换即可：
$V_{i}\rightarrow V_{i}^{'}=V_{i}\sqrt{\frac{T0}{T}}$
文件initialize.cu中的initialize_position()函数和initialize_velocity()函数分别负责坐标和速度的初始化。其中initialize_velocity()函数调用了同文件的scale_velocity()函数，用以得到一个确定的初始温度。首先第一个需要包含的文件时时common.cuh包含了原子结构体，其中存储了所有原子的位置和动量信息：

#pragma once

#ifdef USE_DP
    typedef double real;
#else
    typedef float real;
#endif

const real K_B = 8.617343e-5;
const real TIME_UNIT_CONVERSION = 1.018051e+1;

struct Atom
{
    int *NN;
    int *NL;
    real *m;
    real *x;
    real *y;
    real *z;
    real *vx;
    real *vy;
    real *vz;
    real *fx;
    real *fy;
    real *fz;
    real *pe;
    real *ke;
    real *box;
};

接下来是initialize.cuh，作为初始化操作的头文件：

#pragma once
#include "common.cuh"

void initialize_position(int nx, real ax, Atom *atom);
void initialize_velocity(int N, real T_0, Atom *atom);

#include "initialize.cuh"
#include <stdlib.h>
#include <math.h>

static void scale_velocity(int N, real T_0, Atom *atom)
{
    real *m = atom->m;
    real *vx = atom->vx;
    real *vy = atom->vy;
    real *vz = atom->vz;
    real temperature = 0.0;
    for (int n = 0; n < N; ++n) 
    {
        real v2 = vx[n]*vx[n] + vy[n]*vy[n] + vz[n]*vz[n];     
        temperature += m[n] * v2; 
    }
    temperature /= 3.0 * K_B * N;
    real scale_factor = sqrt(T_0 / temperature);
    for (int n = 0; n < N; ++n)
    { 
        vx[n] *= scale_factor;
        vy[n] *= scale_factor;
        vz[n] *= scale_factor;
    }
}

void initialize_position(int nx, real ax, Atom *atom)
{
    atom->box[0] = ax * nx;
    atom->box[1] = ax * nx;
    atom->box[2] = ax * nx;
    atom->box[3] = atom->box[0] * 0.5;
    atom->box[4] = atom->box[1] * 0.5;
    atom->box[5] = atom->box[2] * 0.5;
    real *x = atom->x;
    real *y = atom->y;
    real *z = atom->z;
    real x0[4] = {0.0, 0.0, 0.5, 0.5};
    real y0[4] = {0.0, 0.5, 0.0, 0.5}; 
    real z0[4] = {0.0, 0.5, 0.5, 0.0};
    int n = 0;
    for (int ix = 0; ix < nx; ++ix)
    {
        for (int iy = 0; iy < nx; ++iy)
        {
            for (int iz = 0; iz < nx; ++iz)
            {
                for (int i = 0; i < 4; ++i)
                {
                    x[n] = (ix + x0[i]) * ax;
                    y[n] = (iy + y0[i]) * ax;
                    z[n] = (iz + z0[i]) * ax;
                    n++;
                }
            }
        }
    }
}
  
void initialize_velocity(int N, real T_0, Atom *atom)
{
    real *m = atom->m;
    real *vx = atom->vx;
    real *vy = atom->vy;
    real *vz = atom->vz;
    real momentum_average[3] = {0.0, 0.0, 0.0};
    for (int n = 0; n < N; ++n)
    { 
        vx[n] = -1.0 + (rand() * 2.0) / RAND_MAX; 
        vy[n] = -1.0 + (rand() * 2.0) / RAND_MAX; 
        vz[n] = -1.0 + (rand() * 2.0) / RAND_MAX;    
        
        momentum_average[0] += m[n] * vx[n] / N;
        momentum_average[1] += m[n] * vy[n] / N;
        momentum_average[2] += m[n] * vz[n] / N;
    } 
    for (int n = 0; n < N; ++n) 
    { 
        vx[n] -= momentum_average[0] / m[n];
        vy[n] -= momentum_average[1] / m[n];
        vz[n] -= momentum_average[2] / m[n]; 
    }
    scale_velocity(N, T_0, atom);
}

13.1.4 边界条件

在我们对分子动力学模拟的定义中，除了初始条件外，还提到了边界条件。边界条件对常微分方程的求解不是必要的，但在分子动力学模拟中通常会根据所模拟的物理体系来选取合适的边界条件，以期达到更合理的结果。边界条件的选取对粒子见的作用力的计算也是有影响的。常见的边界条件有好几种，但我们这里只讨论其中的一种：周期边界条件。在计算机模拟中，模拟的系统尺寸一定是有限的，通常比实验中对应的体积的尺寸小很多。选取周期边界条件通常可以让模拟的体系更加接近于实际的情形，因为原本有边界的系统在应用了周期边界条件之后，似乎没有辩解了。当然，并不能说应用了周期边界条件的系统就等价于无限大的系统，只能说周期边界条件的的应用可以部分地消除边界效应，让所模拟系统的性质更加接近于无限大系统的性质。通常，在这种条件下，我们要模拟几个不同大小的系统，分析所得结果对模拟尺寸的依赖关系。
再计算两个粒子，如粒子i和粒子j的距离时，就要考虑周期边界条件带来的影响。举个一维的例子，假设模拟在一个长度为Lx的模拟盒子中进行，采用周期性边界条件时，必须将一维的盒子想象成一个圆圈。假设Lx=10(任意距离)，第i个粒子的坐标xi=1,第j个粒子的坐标xj=8,则这两个粒子之间的距离是多少呢？如果忽略周期边界条件，那么答案是|xi-xj|=7,而j粒子在i离子的左边，且坐标值可以平移至8-10=-2.这样，j与i的距离时|xj-xi|=3,比平移粒子之前两个粒子之间的距离要小。在我们的模拟过程中，总是采取最小镜像约定：在计算两个粒子的距离时，总是取最小的可能值。定义
xj-xi=xij
则这个约定等价于如下规则：
(1)如果xij<-Lx/2,则将xij换为xij+Lx
(2)如果xij>Lx/2,则将xij换为xij-Lx
最终结果就是让变换后的xij的绝对值小于Lx/2.
很容易将上述讨论推广到二维和三维的情形。例如，在二维的情形中，将要将一个周期的模拟盒子想象成一个环面(torus)，就像一个甜甜圈或一个充了气的轮胎。在三维的情形下，就要将一个盒子想象为一个三维环面，而最小镜像约定可以做如下表述：
(1)如果xij<-Lx/2,则将xij换为xij+Lx
(2)如果xij>Lx/2,则将xij换为xij-Lx
(3)如果yij<-Ly/2,则将yij换为yij+Ly
(4)如果xij>Ly/2,则将yij换为yij-Ly
(5)如果zij<-Lz/2,则将zij换为zij+Lz
(6)如果zij>Lz/2,则将zij换为zij-Lz
这里，我们假设了三维模拟盒子中三个共点的边的长度分别为Lx,Ly,Lz,且两两互相垂直(所谓的正交盒子)。如果有任意两个共点的边不是相互垂直的，情形就要复杂一些。在通用的分子动力学模拟程序中，必须考虑这种复杂的情形，但本书只讨论正交模拟盒子的情形。
文件mic.cuh中的apply_,ic()函数负责实现最小镜像约定。该文件被定义为一个头文件，但其中定义了一个静态函数apply_mic()。该文件将仅被源文件neighbor.cu和force.cu包含。这就相当于将静态函数apply_mic()写在源文件neighbor.cu和force.cu的开头部分。这样，在编译两个源文件时，静态函数apply_mic()将被反复调用，使用内联函数更高效。若将appluy_mic()放在一个单独的编译单元，则不可能将其编译成内联函数。以下是mic.cu的代码

#pragma once

static void apply_mic(real *box, real *x12, real *y12, real *z12)
{
    if      (*x12 < - box[3]) { *x12 += box[0]; } 
    else if (*x12 > + box[3]) { *x12 -= box[0]; }
    if      (*y12 < - box[4]) { *y12 += box[1]; } 
    else if (*y12 > + box[4]) { *y12 -= box[1]; }
    if      (*z12 < - box[5]) { *z12 += box[2]; } 
    else if (*z12 > + box[5]) { *z12 -= box[2]; }
}

13.1.5 相互作用

宏观物质的性质在很大程度上是由微观粒子之间的相互作用力决定的。所以，对粒子间相互作用力的计算在分子动力学模拟中是至关重要的。粒子间有何种相互作用不是分子动力学模拟本身所能回答的；它本质上是一个1量子力学的问题。在经典分子动力学模拟中，粒子间的相互作用力常常由一个或多个经验势函数描述。经典势函数能后在某种程度上反映出某些物质的某些性质。现在已发展出很多这样的势函数。这里，我们只介绍本书将要用到的一种势函数-Lennard-Jones势函数。考察系统中的任意粒子对i和j，他们之间的相互势能可以写成：
         $U_{ij}(r_{ij})=4\epsilon (\frac{\sigma ^{12}}{r_{ij}^{12}})-\frac{\sigma ^{6}}{r_{ij}^{6}})$
其中 $\epsilon$ 和 $\sigma$ 是势函数中的参量，分别具有能量和长度的量纲；rij=|rj-ri|是两个粒子间的距离。
Lennard-Jones势比较适合描述惰性元素组成的物质。它是最早提出的两体势函数之一。两体势是指两个粒子i和j之间的相互作用势仅依赖于他们之间的距离rij，不依赖于系统中其他粒子的存在与否及具体位置。对于这样的势函数，我们将整个系统的势能U写为
         $U=\sum_{i=1}^{N}U_{i}$
   $U_{i}=\frac{1}{2}\sum_{j\neq i}^{}U_{ij}(r_{ij}))$
将以上两式合起来，可以写成：
$U=\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}^{}U_{ij}(r_{ij}))$
上面的Ui可以称为粒子i的势能。上式中的因子1/2的作用是防止将一对相互作用的势能重复计入。例如假如系统中只有三个粒子，那么总势能就是
$U=\frac{1}{2}[U_{12)}(r_{12})+U_{13)}(r_{13})+U_{21}(r_{21})+U_{23}(r_{23})+U_{31)}(r_{31})+U_{32)}(r_{32})]$
对于两体势，我们有：
$U=U_{12}(r_{12})+U_{13}(r_{13})+U_{23}(r_{23})$
如果没有那个因子1/2,那么得到的势能可以写为如下形式：
$U=\sum_{i=1}^{N}\sum_{j> i}^{}U_{ij}(r_{ij}))$ 这种写法省略的一般的计算，正是C++版本的程序中使用的公式。
能够由相互作用势能描述的粒子系统，其粒子间的相互作用力是所谓的保守力，能够表达为势能对坐标的负梯度。具体来说，粒子i所受的合外力可以表达为:
$F_{i}=-\bigtriangledown _{i}U$
其中， $\bigtriangledown _{i}$ 表示针对粒子i的梯度符号，可以写为
$\bigtriangledown _{i}=\frac{\partial }{\partial x_{i}}e_{x}+\frac{\partial }{\partial y_{i}}e_{y}+\frac{\partial }{\partial z_{i}}e_{z}$
这里的ex、ey、ez是是三个方向矢量。通过推导，我们可以得到:
$F_{i}=\sum_{j\neq i}^{}\frac{\partial U_{ij}(r_{ij}))}{\partial r_{ij})}$
我们定义一个表示粒子间相对位置的符号，我们也可以将Fi写为如下形式
$F_{i}=\sum_{j\neq i}^{}F_{ij}$
也就是说，粒子i受到总的力，等于其他粒子对他的作用力Fij的矢量和。这里Fij要理解为j粒子对i粒子的作用力，或者说i粒子受到来自于j粒子的作用力且有Fij=-Fji,这正是牛顿第三定律。在C++版本中，将利用牛顿第三定律节省一半的计算。
对于Lennard-Jones势，我们可以推导出如下的表达式：
   $F_{ij}=(\frac{24\epsilon \sigma^{6 }}{r_{ij}^{8}}-\frac{48\epsilon \sigma^{12 }}{r_{ij}^{14}})$
计算原子的受力及势能的函数find_force()定义在文件force.cu中，该函数在结构上类似为第9章讨论过的求邻居列表的函数，具体算法如下:
(1)第16-26行计算一些常量，包括阶段距离的平方 $r_{c}^{2}$ (cutoff_squere)、 $24\epsilon \sigma ^{6}(e24s6)$ 等。在循环之前尽可能多地计算常量可以省去很多不必要的计算。
(2)第27-30行将各个粒子的力和势能初始化为零。
(3)第31行对各个粒子进行循环。
(4)第33行对粒子i的所有邻居进行循环
(5)第35行找到粒子i的一个邻居j.这里所用的邻居列表由neighbor.cu中的find_neighbor()函数构建，具体算法在第9章构建过。
(6)第36行排除j小于i的情况。
(7)第37-39行计算从i到j的位置矢量rij
(8)第40行对位置矢量rij实施最小镜像约定。
(9)第41行计算i和j两个粒子之间的距离平方 $r_{c}^{2}$
(10)第42行排除 $r_{ij}^{2}>r_{c}^{2}$ .
(11)第43-48行用经济的方式计算 $r_{ij}^{-6}$ 等，尽量有惩罚代替除法。
(12)第49-53行对每个粒子的力和势能进行累加

以下是neighbor.cu

#pragma once
#include "common.cuh"

void find_neighbor(int N, int MN, Atom *atom);

#include "neighbor.cuh"
#include "mic.cuh"
#include <stdio.h>
#include <stdlib.h>

void find_neighbor(int N, int MN, Atom *atom)
{
    int *NN = atom->NN;
    int *NL = atom->NL;
    real *x = atom->x;
    real *y = atom->y;
    real *z = atom->z;
    real *box = atom->box; 
    real cutoff = 11.0;
    real cutoff_square = cutoff * cutoff;

    for (int n = 0; n < N; n++)
    {
        NN[n] = 0;
    }

    for (int n1 = 0; n1 < N - 1; n1++)
    {  
        for (int n2 = n1 + 1; n2 < N; n2++)
        {   
            real x12 = x[n2] - x[n1];
            real y12 = y[n2] - y[n1];
            real z12 = z[n2] - z[n1];
            apply_mic(box, &x12, &y12, &z12);
            real d_square = x12*x12 + y12*y12 + z12*z12;

            if (d_square < cutoff_square)
            {        
                NL[n1 * MN + NN[n1]++] = n2;
                NL[n2 * MN + NN[n2]++] = n1;
            }
        }
    }

    for (int n1 = 0; n1 < N - 1; n1++)
    {
        if (NN[n1] > MN)
        {
            printf("Error: MN is too small.\n");
            exit(1);
        }
    } 
}

以下是force.cu

#include "force.cuh"
#include "mic.cuh"

void find_force(int N, int MN, Atom *atom)
{
    int *NN = atom->NN;
    int *NL = atom->NL;
    real *x = atom->x;
    real *y = atom->y;
    real *z = atom->z;
    real *fx = atom->fx;
    real *fy = atom->fy;
    real *fz = atom->fz;
    real *pe = atom->pe;
    real *box = atom->box;
    const real epsilon = 1.032e-2;
    const real sigma = 3.405;
    const real cutoff = 10.0;
    const real cutoff_square = cutoff * cutoff;
    const real sigma_3 = sigma * sigma * sigma;
    const real sigma_6 = sigma_3 * sigma_3;
    const real sigma_12 = sigma_6 * sigma_6;
    const real e24s6 = 24.0 * epsilon * sigma_6; 
    const real e48s12 = 48.0 * epsilon * sigma_12;
    const real e4s6 = 4.0 * epsilon * sigma_6; 
    const real e4s12 = 4.0 * epsilon * sigma_12;
    for (int n = 0; n < N; ++n) 
    { 
        fx[n] = fy[n] = fz[n] = pe[n] = 0.0; 
    }
    for (int i = 0; i < N; ++i)
    {
        for (int k = 0; k < NN[i]; k++)
        {
            int j = NL[i * MN + k];
            if (j < i) { continue; }
            real x_ij = x[j] - x[i];
            real y_ij = y[j] - y[i];
            real z_ij = z[j] - z[i];
            apply_mic(box, &x_ij, &y_ij, &z_ij);
            real r2 = x_ij*x_ij + y_ij*y_ij + z_ij*z_ij;
            if (r2 > cutoff_square) { continue; }
            real r2inv = 1.0 / r2;
            real r4inv = r2inv * r2inv;
            real r6inv = r2inv * r4inv;
            real r8inv = r4inv * r4inv;
            real r12inv = r4inv * r8inv;
            real r14inv = r6inv * r8inv;
            real f_ij = e24s6 * r8inv - e48s12 * r14inv;
            pe[i] += e4s12 * r12inv - e4s6 * r6inv;
            fx[i] += f_ij * x_ij; fx[j] -= f_ij * x_ij;
            fy[i] += f_ij * y_ij; fy[j] -= f_ij * y_ij;
            fz[i] += f_ij * z_ij; fz[j] -= f_ij * z_ij;
        }
    }
}

13.1.6 运动方程的数值积分

在经典力学中，粒子的运动方程可以用牛顿第二定律进行表达。例如，对于第i个粒子，其运动方程为：
         $m_{i}\frac{d^{2}r_{i}}{dt^{2}}=F_{i}$
其中,Fi是该粒子受到的总的力；d2ri/dt2是位置矢量对时间的二阶导数，即加速度。这是一个二阶常微分方程，我们可以把它改写成两个一阶常微分方程：
         $\frac{dr_{i}}{dt}=v_{i}$
         $\frac{dv_{i}}{dt}=\frac{F_{i}}{m_{i}}$
其中上面第一个式子就是速度的定义。
对于弄懂方程进行积分的目的就是在给定的初始条件下找到各粒子在一系列离散的时间点的坐标和速度值。我们假设每两个离散的时间点之间的间隔是固定的，记为 $\Delta t$ ，称为时间步长.在分子动力学模拟中使用的数值积分方法有很多种，我们这里只介绍“速度-Verlot”积分方法(verlet是发展分子动力学模拟方法的先驱之一)，而且不去追究其推导过程。在该方法中，第i个粒子在时刻t+ $\Delta t$ 的速度vi(t+ $\Delta t$ )和位置ri(t+ $\Delta t$ )分别由以下两个式子进行给出：
         $v_{i}(t+\Delta t)=v_{i}(t)+\frac{1}{2}\frac{F_{i}(t)+F_{i}(t+\Delta t)}{m_{i}}\Delta t$
         $r_{i}(t+\Delta t)=r_{i}(t)+v_{i}(t)\Delta t+\frac{1}{2}\frac{F_{i}(t))}{m_{i}}\Delta t^{^{2}}$
由上式可以看出，t+ $\Delta t$ 时刻的坐标仅依赖于t时刻的坐标、速度和力，但t+ $\Delta t$ 时刻的速度依赖于t时刻的速度、力及t+ $\Delta t$ 时刻的力，以上两式应该对应如下的计算流程：
   $v_{i}(t)\rightarrow v_{i}(t+\Delta t/2)=v_{i}(t)+\frac{1}{2}\frac{F_{i}(t)}{m_{i}}\Delta t$
         $r_{i}(t)\rightarrow r_{i}(t+\Delta t)=r_{i}(t)+ v_{i}(t+\Delta t/2)\Delta t$
         $F_{i}(t)\rightarrow F_{i}(t+\Delta t)$
         $v_{i}(t+\Delta t/2)\rightarrow v_{i}(t+\Delta t)= r_{i}(t+\Delta t/2)+\frac{1}{2}\frac{F_{i}(t+\Delta t)}{m_{i}}\Delta t$
注意：我们引入了一个中间的速度变量vi(t+ $\Delta t/2$ ),上面第三个式子的意思是用t+ $\Delta t$ 时刻的坐标计算新的力Fi(t+ $\Delta t$ ),替换掉老的力Fi(t)。完成上述计算之后，离子的坐标、速度、和力都从t时刻更新为t+ $\Delta t$ 时刻的。这就是一个时间步的计算，反复执行这样的计算过程，系统的微观状态就会随着时间而变化，从而得到一条相空间的轨迹。系统的所有宏观性质都包含在相轨迹中。
第一个式子、第二个式子、第四个式子由文件integrate.cu中的integrate()函数实现。第3个式子对应文件force.cu中的find_force()函数。我们分子动力学模拟中有两个演化过程，一个是平衡阶段，另一个是产出阶段，他们分别对应文件integrate.cu中的equilibration()函数和production()函数，其中第一个函数会调用同文件的scale_velocity()函数控制体系的温度，而后者会调用同文件的sum()函数计算体系的总动能和总势能。其中integrate.cu代码如下：

#pragma once
#include "common.cuh"

void equilibration
(
    int Ne, int N, int MN, real T_0, 
    real time_step, Atom *atom
);

void production
(
    int Np, int Ns, int N, int MN, real T_0, 
    real time_step, Atom *atom
);

#include "integrate.cuh"
#include "force.cuh"
#include "error.cuh"
#include <stdio.h>
#include <math.h>
#include <time.h>

static real sum(int N, real *x)
{
    real s = 0.0;
    for (int n = 0; n < N; ++n) 
    {
        s += x[n];
    }
    return s;
}

static void scale_velocity(int N, real T_0, Atom *atom)
{
    real temperature = sum(N, atom->ke) / (1.5 * K_B * N);
    real scale_factor = sqrt(T_0 / temperature);
    for (int n = 0; n < N; ++n)
    { 
        atom->vx[n] *= scale_factor;
        atom->vy[n] *= scale_factor;
        atom->vz[n] *= scale_factor;
    }
}

static void integrate
(int N, real time_step, Atom *atom, int flag)
{
    real *m = atom->m;
    real *x = atom->x;
    real *y = atom->y;
    real *z = atom->z;
    real *vx = atom->vx;
    real *vy = atom->vy;
    real *vz = atom->vz;
    real *fx = atom->fx;
    real *fy = atom->fy;
    real *fz = atom->fz;
    real *ke = atom->ke;
    real time_step_half = time_step * 0.5;
    for (int n = 0; n < N; ++n)
    {
        real mass_inv = 1.0 / m[n];
        real ax = fx[n] * mass_inv;
        real ay = fy[n] * mass_inv;
        real az = fz[n] * mass_inv;
        vx[n] += ax * time_step_half;
        vy[n] += ay * time_step_half;
        vz[n] += az * time_step_half;
        if (flag == 1) 
        { 
            x[n] += vx[n] * time_step; 
            y[n] += vy[n] * time_step; 
            z[n] += vz[n] * time_step; 
        }
        else
        {
            real v2 = vx[n]*vx[n] + vy[n]*vy[n] + vz[n]*vz[n];
            ke[n] = m[n] * v2 * 0.5;
        }
    }
}

void equilibration
(
    int Ne, int N, int MN, real T_0, 
    real time_step, Atom *atom
)
{
    find_force(N, MN, atom);
    for (int step = 0; step < Ne; ++step)
    { 
        integrate(N, time_step, atom, 1);
        find_force(N, MN, atom);
        integrate(N, time_step, atom, 2);
        scale_velocity(N, T_0, atom);
    }
}

void production
(
    int Np, int Ns, int N, int MN, real T_0, 
    real time_step, Atom *atom
)
{
    float t_force = 0.0f;

    clock_t t_total_start = clock();

    FILE *fid = fopen("energy.txt", "w");
    for (int step = 0; step < Np; ++step)
    {
        integrate(N, time_step, atom, 1);

        clock_t t_force_start = clock();

        find_force(N, MN, atom);

        clock_t t_force_stop = clock();

        t_force += float(t_force_stop - t_force_start) / CLOCKS_PER_SEC;

        integrate(N, time_step, atom, 2);

        if (0 == step % Ns)
        {
            fprintf(fid, "%g %g\n", sum(N, atom->ke), sum(N, atom->pe));
        }
    }
    fclose(fid);

    clock_t t_total_stop = clock();

    float t_total = float(t_total_stop - t_total_start) / CLOCKS_PER_SEC;
    printf("Time used for production = %g s\n", t_total);
    printf("Time used for force part = %g s\n", t_force);
}

13.1.7 程序中使用的单位制

我们的分子动力学模拟程序中只涉及经典力学和热力学，所以只需用到四个基本物理量的单位。我们选择如下四个基本单位来确定各个物理量的数值：
(1)能量:电子伏特(记号为eV),约为1.6×10-19J、
(2)长度：埃(angstrom,记号为 $\dot{A}$ )即10-10m。
(3)质量：原子质量单位(atomic mass unit,记号为amu)，约为1.66×10-27kg。
(4)温度：开尔文(记号为K).
用这样的基本单位，可使程序中大部分物理量的数值接近1.我们称这样的单位为程序的“自然单位”。
从以上基本单位可以推导出程序中其他相关物理量的单位：
(1)力。因为力乘以距离等于功(能量)，所以力的单位是能量单位除以长度单位即eV $\cdot \dot{A}^{-1}$ .
(2)速度。因为动能正比于质量乘以速度的平方，所以速度的单位是能量单位除以质量单位再开根号，即eV $^{1/2}amu^{-1/2}$ .
(3)时间。因为动能正比于质量乘以速度的平方，所以速度的单位是能量单位除以速度单位，即 $\dot{A}amu^{1/2}eV^{-1/2}$ 约为1.018051×10fs(飞秒，即10-15s)。在程序的初始化阶段，我们需要设置一个积分的步长。一般来说，我们习惯用飞秒转换为该程序的自然单位，这就需要除以以下常数：
const real TIME_UNIT_CONVERSION=1.018051e+1;
该常数在头文件common.cuh中定义。
(4)玻尔兹曼常数kB。这是一个很重要的常数，他在国际单位制中约为1.38×10-23J.K-1，对应于程序自然单位制8.617343×10-5eV.K-1.所以，在开头的common.cuh还定义了一个常数：
const real K_B=8.617343e-5;

13.1.8 程序的编译和运行

我们在windows使用从网站GNU make for Windows
下载得到的make使用以下命令进行编译：make -f makefile.windows编译出ljmd.exe,后使用
ljmd nx Ne进行运行，这里，nx代表我们所模拟的固态氮模型在每个方向的晶胞个数，Ne代表平衡过程中演化的步数。

13.2 CUDA版本的分子动力学模拟程序开发

在13.1节中，我们介绍了一个完整的C++版本的分子动力学模拟程序，本节将逐步移植到CUDA程序并优化程序的性能。

13.2.1 仅加速求力和能量的部分

从13.1节的结果可知，在阐述阶段，求力的部分是最耗时的，所以很自然的想到用一个线程处理一个粒子的受力，可以将原来的求力函数分解成一个核函数和一个包装函数。核函数专注于力的计算，包装函数负责准备核函数需要的参数。因为我们这里仅仅加速求力的部分，所以在包装函数中需要进行数据传输：在调用求力的核函数之前将坐标从主机传输到设备；在调用求力的核函数之后将力与势能的数据从设备传递到主机。包装函数中还要准备核函数的执行配置。例如体系中由256个粒子，而核函数调用时指定的线程块大小为128，则每个线程块处理128个粒子，一共需要两个线程块。可见，在该并行方案中，需要较多的粒子才能达到较高的并行度。根据上述方案，写出以下求力和势能的核函数和包装函数。

#include "error.cuh"
#include "force.h"
#include "mic.h"

struct LJ
{
    real cutoff2;
    real e24s6;
    real e48s12;
    real e4s6;
    real e4s12;
};

static void __global__ gpu_find_force
(
    LJ lj, int N, int *g_NN, int *g_NL, Box box,
    real *g_x, real *g_y, real *g_z,
    real *g_fx, real *g_fy, real *g_fz, real *g_pe
)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N)
    {
        real fx = 0.0;
        real fy = 0.0;
        real fz = 0.0;
        real potential = 0.0;
        int NN = g_NN[i];
        real x_i = g_x[i];
        real y_i = g_y[i];
        real z_i = g_z[i];
        for (int k = 0; k < NN; ++k)
        {
            int j = g_NL[i + N * k];
            real x_ij  = g_x[j] - x_i;
            real y_ij  = g_y[j] - y_i;
            real z_ij  = g_z[j] - z_i;
            apply_mic(box, &x_ij, &y_ij, &z_ij);
            real r2 = x_ij*x_ij + y_ij*y_ij + z_ij*z_ij;
            if (r2 > lj.cutoff2) { continue; }

            real r2inv = 1.0 / r2;
            real r4inv = r2inv * r2inv;
            real r6inv = r2inv * r4inv;
            real r8inv = r4inv * r4inv;
            real r12inv = r4inv * r8inv;
            real r14inv = r6inv * r8inv;
            real f_ij = lj.e24s6 * r8inv - lj.e48s12 * r14inv;
            potential += lj.e4s12 * r12inv - lj.e4s6 * r6inv;
            fx += f_ij * x_ij;
            fy += f_ij * y_ij;
            fz += f_ij * z_ij;
        }
        g_fx[i] = fx;
        g_fy[i] = fy;
        g_fz[i] = fz;
        g_pe[i] = potential * 0.5;
    }
}

void find_force(int N, int MN, Atom *atom)
{
    const real epsilon = 1.032e-2;
    const real sigma = 3.405;
    const real cutoff = 10.0;
    const real cutoff2 = cutoff * cutoff;
    const real sigma_3 = sigma * sigma * sigma;
    const real sigma_6 = sigma_3 * sigma_3;
    const real sigma_12 = sigma_6 * sigma_6;
    const real e24s6 = 24.0 * epsilon * sigma_6;
    const real e48s12 = 48.0 * epsilon * sigma_12;
    const real e4s6 = 4.0 * epsilon * sigma_6;
    const real e4s12 = 4.0 * epsilon * sigma_12;
    LJ lj;
    lj.cutoff2 = cutoff2;
    lj.e24s6 = e24s6;
    lj.e48s12 = e48s12;
    lj.e4s6 = e4s6;
    lj.e4s12 = e4s12;

    Box box;
    box.lx = atom->box[0];
    box.ly = atom->box[1];
    box.lz = atom->box[2];
    box.lx2 = atom->box[3];
    box.ly2 = atom->box[4];
    box.lz2 = atom->box[5];

    int m = sizeof(real) * N;
    CHECK(cudaMemcpy(atom->g_x, atom->x, m, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(atom->g_y, atom->y, m, cudaMemcpyHostToDevice));
    CHECK(cudaMemcpy(atom->g_z, atom->z, m, cudaMemcpyHostToDevice));

    int block_size = 128;
    int grid_size = (N - 1) / block_size + 1;
    gpu_find_force<<<grid_size, block_size>>>
    (
        lj, N,  atom->g_NN, atom->g_NL, box,
        atom->g_x, atom->g_y, atom->g_z,
        atom->g_fx, atom->g_fy, atom->g_fz, atom->g_pe
    );

    CHECK(cudaMemcpy(atom->fx, atom->g_fx, m, cudaMemcpyDeviceToHost));
    CHECK(cudaMemcpy(atom->fy, atom->g_fy, m, cudaMemcpyDeviceToHost));
    CHECK(cudaMemcpy(atom->fz, atom->g_fz, m, cudaMemcpyDeviceToHost));
    CHECK(cudaMemcpy(atom->pe, atom->g_pe, m, cudaMemcpyDeviceToHost));
}

我们先来看包装函数：
(1)包装函数的开头63-73行计算Lennard-Jones势中的一些常数。我们可以将这些常数直接传递给核函数，但为了简洁期间，我们在第74-79行把这些常数封装成一个结构体变量再传递给核函数。该结构体LJ在本文件开头的5-12行定义。类似地，第81-87行将CPU中的数组box[6]包装成一个结构体变量，用于传值给核函数。将结构体传值给核函数，对结构体中的数据的读取将通过常量内存缓存。也可以将这些常数放在一个全局内存数组，但没有使用常量内存高效。读者可以记住一个结论：只要数据量在编译期间就确定且不大(明显少于4KB),在核函数中仅仅被读取，而且一个线程中的所有线程在某个时刻访问同一个地址，就适合用传参的方式使用常量内存。
(2)第89-92行将核函数需要的坐标数据从主机传输到设备。
(3)第94-101行将核函数需要的坐标数据从主机传输到设备。
(4)第103-106行将核函数求出的力和势能数据从设备传输到主机。
再看核函数的实现：
(1)第21行定义了从线程指标到粒子指标i的映射关系，即一个线程对应一个粒子。
(2)第22行将原来的C++函数中对粒子i的遍历for(int i=0;i<N;++i)改成了判断if(i<N).
(3)第22行后的代码与原来C++的版本的代码几乎一模一样。对于每一个粒子i,我们继续对他所有邻居j进行遍历，计算i与j之间的力和势能，并累加起来，将最终结果保存到全局内存。
因为求力的函数中需要用到邻居列表，我们也在GPU中建立了邻居列表。我们在第9章中已经讨论过构件邻居列表的CUDA实现。本章对应程序中的CUDA实现与第9章给出的一模一样，区别仅在于第9章没考虑周期边界条件，而分子动力学模拟中，施加了最小镜像约定，从而考虑了周期边界条件。

13.2.2 全部加速计算

#include "integrate.h"
#include "error.cuh"
#include "force.h"
#include "reduce.h"
#include <stdio.h>
#include <math.h>
#include <time.h>

static void __global__ gpu_scale_velocity
(
    int N, real scale_factor, 
    real *g_vx, real *g_vy, real *g_vz
)
{
    int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    { 
        g_vx[n] *= scale_factor;
        g_vy[n] *= scale_factor;
        g_vz[n] *= scale_factor;
    }
}

static void scale_velocity(int N, real T_0, Atom *atom)
{
    real temperature = sum(N, atom->g_ke) / (1.5 * K_B * N);
    real scale_factor = sqrt(T_0 / temperature);

    int block_size = 128;
    int grid_size = (N - 1) / block_size + 1;
    gpu_scale_velocity<<<grid_size, block_size>>>
    (N, scale_factor, atom->g_vx, atom->g_vy, atom->g_vz);
}

static void __global__ gpu_integrate
(
    int N, real time_step, real time_step_half,
    real *g_m, real *g_x, real *g_y, real *g_z,
    real *g_vx, real *g_vy, real *g_vz,
    real *g_fx, real *g_fy, real *g_fz, 
    real *g_ke, int flag
)
{
    int n = blockDim.x * blockIdx.x + threadIdx.x;
    if (n < N)
    {
        real mass = g_m[n];
        real mass_inv = 1.0 / mass;
        real ax = g_fx[n] * mass_inv;
        real ay = g_fy[n] * mass_inv;
        real az = g_fz[n] * mass_inv;
        real vx = g_vx[n];
        real vy = g_vy[n];
        real vz = g_vz[n];

        vx += ax * time_step_half;
        vy += ay * time_step_half;
        vz += az * time_step_half;
        g_vx[n] = vx;
        g_vy[n] = vy;
        g_vz[n] = vz;

        if (flag == 1) 
        { 
            g_x[n] += vx * time_step; 
            g_y[n] += vy * time_step; 
            g_z[n] += vz * time_step; 
        }
        else
        {
            g_ke[n] = (vx*vx + vy*vy + vz*vz) * mass * 0.5;
        }
    }
}

static void integrate(int N, real time_step, Atom *atom, int flag)
{
    real time_step_half = time_step * 0.5;

    int block_size = 128;
    int grid_size = (N - 1) / block_size + 1;
    gpu_integrate<<<grid_size, block_size>>>
    (
        N, time_step, time_step_half,
        atom->g_m, atom->g_x, atom->g_y, atom->g_z,
        atom->g_vx, atom->g_vy, atom->g_vz,
        atom->g_fx, atom->g_fy, atom->g_fz, 
        atom->g_ke, flag
    );
}

void equilibration
(
    int Ne, int N, int MN, real T_0, 
    real time_step, Atom *atom
)
{
    find_force(N, MN, atom);
    for (int step = 0; step < Ne; ++step)
    { 
        integrate(N, time_step, atom, 1);
        find_force(N, MN, atom);
        integrate(N, time_step, atom, 2);
        scale_velocity(N, T_0, atom);
    } 
}

void production
(
    int Np, int Ns, int N, int MN, real T_0, 
    real time_step, Atom *atom
)
{
    float t_force = 0.0f;
    CHECK(cudaDeviceSynchronize());
    clock_t t_total_start = clock();

    FILE *fid_e = fopen("energy.txt", "w");
    for (int step = 0; step < Np; ++step)
    {  
        integrate(N, time_step, atom, 1);

        CHECK(cudaDeviceSynchronize());
        clock_t t_force_start = clock();

        find_force(N, MN, atom);

        CHECK(cudaDeviceSynchronize());
        clock_t t_force_stop = clock();

        t_force += float(t_force_stop - t_force_start) / CLOCKS_PER_SEC;

        integrate(N, time_step, atom, 2);

        if (0 == step % Ns)
        {
            fprintf(fid_e, "%g %g\n", sum(N, atom->g_ke), sum(N, atom->g_pe));
        }
    }
    fclose(fid_e);

    clock_t t_total_stop = clock();

    float t_total = float(t_total_stop - t_total_start) / CLOCKS_PER_SEC;
    printf("Time used for production = %g s\n", t_total);
    printf("Time used for force part = %g s\n", t_force);
}

布林组-？

关注

0
点赞
踩
4

收藏

觉得还不错? 一键收藏
1
评论
第13章分子动力学模拟的CUDA程序的开发

在前面的章节中，已经用一些简短的程序较为系统的介绍了CUDA编程的基础知识。为了能改有更加深入地讨论CUDA程序的优化策略，将讨论一个出具规模的CUDA程序的开发与优化。该程序可实现一个简单地分子动力学模拟。分子动力学模拟的算法涉及本科水平的物理学(经典力学和经典统计力学)。
复制链接

扫一扫