24小时热门版块排行榜

返回列表

【奖励】本帖被评价13次，作者yahoohoo增加金币 11.5 个

本帖产生 1 个模拟EPI ，点击这里进行查看

yahoohoo

铁杆木虫 (著名写手)

[资源] 【原创】分子模拟实例与挑战专题之二加速能量与力的计算

分子模拟实例与挑战专题第二期
加速能量与力的计算

本帖内容，非注明参考文献出处的，版权归本人所有，请勿以任何形式转载或传播。对于以直接或间接方式使用文中涉及代码而引起的损失，本人概不负责。

在分子模拟中，能量（Monte Carlo）与力（Molecular Dynamics）的计算最为耗时，而非键合（non-bonded）能量与力的计算占据其中最为主要的部分。我们通常采用截断势能来表示这种非键合的相互作用。也就是说，对于两体（two-body）相互作用，我们只需要考虑当两两距离小于截断半径$r_c$的所有粒子对即可。最为直观的做法就是遍历所有的粒子对：假设总的粒子数目为$N$，那么我们总共需要考虑$N(N-1)/2$对（即便这其中很多对的距离大于$r_c$）。伪代码如下：

  for (int k = 0; k < N - 1; ++k)
for (int n = k + 1; n < N; ++n)
   calcPairForce(k, n);

需要注意的是，虽然这种所谓的”Brute force”方法在实际中几乎不被使用，但它仍然是我们用来判断其他技巧（例如Cell linked list或Verlet list）正确与否的有力手段。

通过上述的讨论，我们知道”Brute force”方法计算能量或力的成本为$\mathcal{O}(N^2)$。这表明该方法用于大量数目粒子的模拟是不可行的。我们这里介绍两种经典的、计算成本为$\mathcal{O}(N)$的方法：（1）Cell linked list；（2）Verlet list。

(1) Cell linked list

Cell linked list的基本思路就是将模拟盒子划分成若干个小的元胞（cell），而这些元胞的边长不小于截断半径$r_c$。根据粒子的位置，我们可以将所有的粒子分配到对应的元胞中。这样做的好处是对于任一粒子$i$（处于元胞$\mathcal{C}_\alpha$中），我们现在只需要考虑元胞$\mathcal{C}_\alpha$及其近邻的26个元胞$\mathcal{C}_\beta$中所有粒子与$i$的相互作用即可。计算成本包括（i）分配粒子到对应的元胞中（$\mathcal{O}(N)$）和（ii）搜索近邻元胞（$\mathcal{O}(N)$）。那么总的计算成本正比于粒子数目$N$。我们以代码的形式介绍下建立列表的方法：

  for (int n = 0; n < NCELL; ++n)
Head[n] = -1; // Initialize the head list
  for (int k = 0; k < N; ++k) {
int celli = CellNo(k); // #cell of particle k
LList[k] = Head[celli];
Head[celli] = k;
  }

这里我们只需要两个一维数组Head[]和LList[]。其中Head[]对应所有元胞的序号最大的粒子，比如Head[1] = 9表示第2个元胞中序号最大的粒子为9，或者说第10个粒子（注：c++中以0为offset，不同与fortran中以1为offset），Head[2]=-1则表示第3个元胞中没有粒子。而LList[]则是粒子与粒子连接的列表，比如LList[9] = 5表示第10个粒子与第6个粒子同处于一个元胞中，LList[5] = -1表示第6个粒子为其所在元胞的最后一个粒子。下面我们给出利用这Head[]和LList[]这两个列表计算力的代码：

  for (int n = 0; n < NCELL; ++n) {
int o = Head[n];
while (i > -1) {
   // Loop over the particles in the cell where particle i is
   int j = LList[o];
   while (j > -1) {
      calcPairForce(o, j);
      j = LList[j];
   }
   // Loop over particles in the neighboring cells
   for (int k = 0; k < NEIGHCELL; ++k) {
      int cell = CLIST[n][k];
      j = Head[j];
      while (j > -1) {
      calcPairForce(o, j);
      j = LList[j];
      }
   }
}
i = LList[o];
}

上述代码中，CLIST[n][k]对应第n+1个元胞的第k+1个近邻元胞的序号。我们采用int类型（4个字节）的数组，那么Cell linked list方法所需要的内存为 NCELL * 4 Byte + N * 4 Byte + NCELL * 26 * 4 Byte。对于一百万个粒子的体系，如果平均每个元胞中的粒子数为27，那么需要的内存为8 MiB。

在分子动力学模拟中，由于粒子遵守牛顿第三定律，因此对于每一个元胞，我们只需要考虑其近邻的13个元胞（2D是4个）。我们在这里略去如何建立元胞近邻列表CLIST[][]，请读者自行思考解之。

(2) Verlet list

与上述的Cell linked list不同，Verlet list的基本思路是：我们只需要考虑以每个粒子$i$为球心、半径为$r_c$的球内包括的所有粒子$j$与$i$的相互作用。因此，使用Verlet list，我们一共需要考虑$N \times 4/3 \times \pi r_c^3 \rho$个粒子对，而使用Cell linked list则需要考虑$N \times 27 \times r_c^3 \rho$。在计算力这一点上，Verlet list的效率大约是cell linked list的7倍。不过，为了建立Verlet list列表，我们通常需要扫描所有的粒子对$ij$以确定粒子$j$是否是$i$的近邻，这就意味者建立Verlet list的成本为$\mathcal{O}(N^2)$，相比Cell linked list 的$\mathcal{O}(N)$而言，这是非常痛苦的。另外，很显然，根据$r_c$为半径所建立的列表是没有任何实用价值的。因为如果这么做，那每一步我们都需要重新更新列表。实际上我们以$r_c+r_s$为半径建立Verlet list，好处是不需要每步更新，缺点是在列表中引入了一些根本不需要考虑的作用对。我们可以通过选择最佳的$r_s$使计算成本最低。我们以具体的代码来演示如何建立并使用Verlet list：

  // Generate the list
  for (int k = 0; k < N – 1; ++k)
for (int n = k + 1; n < N; ++n)
   if (DistSQ(k, n) <= Rv2) VList[k][NList[k]++] = n;

这里，$Rv2=(r_c+r_s)^2$。NList[k]对应第k+1个粒子的近邻数目，而VList[k][j]对应第k+1个粒子的第j+1个近邻粒子的序号。不难发现，对于VList[N][MAXNEIGH]二维数组，我们需要使用一个相对较大的值MAXNEIGH 表示其第二维，一般可以选择$MAXNEIGH = 2 \times 4/3 \pi (r_c+r_s)^3 \rho$。这就说明在内存使用上Verlet list有一定的浪费。有了列表以后，我们便可以计算力或能量：

  for (int k = 0; k < N; ++k) {
int i = NList[k];
for (int n = 0; n < i; ++n)
   calcPairForce(k, VList[k][n]);
  }

目前的一个趋势是将上述的两种方法联合使用，其思路是结合Cell linked list的快速建立列表与Verlet list中大幅减少两两作用对数目的优点，根据元胞列表搜索近邻的作用对，将搜索的作用对数目从$N(N-1)/2$减少到$N\times 27 r_c^3 \rho$或$N\times 14 r_c^3 \rho$。

做分子模拟时，大家关心的一个问题是何时选用哪种方法？这里，我们不妨先考虑下每个粒子的平均作用对数目 $\mathcal{N} \equiv 4/3 \pi r_c^3 \rho$。如果$\mathcal{N}$较小时，可以考虑Cell linked list；反之，则可以考虑Verlet list。影响选择的另外一个重要因素是我们所研究体系中粒子本身的扩散系数。如果粒子扩散很快，那么我们就需要频繁的更新近邻列表。这种情况下，Verlet list很可能就不如Cell linked list有优势；即便是联合使用两种列表方法，也并不保证比单独使用Cell linked list更快。

Figure 1. Efficiency of Verlet list combined with Cell linked list in NVE MD simulation for Lennard-Jones particles at temperature $k_BT = 1.5$ and density $\rho=0.7\sigma^{-3}$ with total simulation time $t=20000 dt$. The time step for MD is $dt = 0.004662\sqrt{m\sigma^2/\epsilon}$. The simulation box is $16\sigma \times 16\sigma \times 16\sigma$. The computations were performed on a Intel Core 2 E8300 (3.16GHz) CPU.

Figure 2. Comparison of Cell linked list and the combined algorithm in NVE MD.

[总结]

Cell linked list 和Verlet list是分子模拟中常用的两种加速能量与力的计算的方法。其根本思想是在限定的局部空间而非整个模拟盒子范围内搜索可能产生相互作用的粒子对。对于简单LJ流体的MD模拟而言，联用的方法最为有效，其计算成本与粒子数成正比。

[ Last edited by yahoohoo on 2010-9-15 at 19:15 ]

24小时热门版块排行榜

[资源] 【原创】分子模拟实例与挑战专题之二加速能量与力的计算

» 收录本帖的淘帖专辑推荐

» 猜你喜欢

» 本主题相关商家推荐: (我也要在这里推广)

» 本主题相关价值贴推荐，对您同样有帮助:

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★ 三星级,支持鼓励

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

把自己的编程经历

★★★ 三星级,支持鼓励

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐

★★★★★ 五星级,优秀推荐