版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

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 ]

回复此楼

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

资源收集

分子模拟

» 猜你喜欢

读博已经有5人回复
到新单位后，换了新的研究方向，没有团队，持续积累2区以上论文，能申请到面上吗已经有13人回复
博士申请都是内定的吗？已经有6人回复
之前让一硕士生水了7个发明专利，现在这7个获批发明专利的维护费可从哪儿支出哈？已经有5人回复
博士读完未来一定会好吗已经有29人回复
投稿精细化工已经有4人回复
高职单位投计算机相关的北核或SCI四区期刊推荐，求支招！已经有4人回复
导师想让我从独立一作变成了共一第一已经有9人回复
心脉受损已经有5人回复
Springer期刊投稿求助已经有4人回复

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

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

1楼2010-09-15 16:55:31

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

yahoohoo

铁杆木虫 (著名写手)

★ ★ ★
御剑江湖(金币+3): 谢谢 2011-04-02 12:11:05

引用回帖:

Originally posted by Wuhe at 2011-04-01 20:07:16:
师兄好，我也是只做分子模拟的，看了师兄临近列表的讲解受益匪浅，自己现在也开始在做相关方面的工作。我一个问题想咨询下师兄，您文章图一（Figure 1）下模拟温度k_BT=1.5,那么温度是多少？是按k_B=1.3806505(24 ...

这里的$k_BT$是约化温度。在模拟中我们使用$\sigma$作为基本的能量单位，那么根据量纲，温度的基本单位$T_0$应该为$T_0 \equiv \sigma / k_B$。我们便用$T_0$来约化温度，讲义中$k_B T = 1.5$是不严格的，严格的描述应该是温度为$T=1.5\sigma / k_B$或$T/T_0=1.5$。

希望这能解释你的问题。

赞一下(1人)

回复此楼

24楼2011-04-01 21:27:47

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

查看全部 29 个回答

yahoohoo

铁杆木虫 (著名写手)

★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢分享 2010-09-16 08:52:05

Cell linked list：

有些文献报道，使用更小边长的元胞（例如$r_c/2$、$r_c/3$），可以减少邻近元胞搜索的粒子对数目。实际上，当元胞边长为$r_c/2$时，需要搜索125个元胞，而这其中的若干元胞很可能是空的。换言之，尽管搜索的粒子对数目减少，但我们需要对更多的近邻元胞（相对于原来的27个）进行循环，而且元胞列表CLIST[][]的大小与边长的立方成反比，因此使用更小的元胞在内存上也是不利的。综合起来，使用更小的元胞并不能提高Cell linked list的效率。

感兴趣的读者可以思考以下两个问题：
(i)将原来的元胞(边长等于$r_c$)划分成若干子元胞（subcell），对于某个元胞及其近邻元胞，我们可以事先确定这些subcell的距离，如果大于$r_c$，那么处于这些subcell中的粒子对则可以不予考虑。这样做是否可以提高Cell linked list的效率？如果可以，需要怎样的条件？

(ii)在每个元胞中按照粒子的坐标对粒子进行排序。这样做会导致建立列表的时间有所增加，那我们可以减少邻近元胞搜索粒子对的数目吗？

Verlet list:

判断更新列表的条件是什么？

[ Last edited by yahoohoo on 2010-9-16 at 02:39 ]

赞一下(1人)

回复此楼

3楼2010-09-16 02:37:36

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖