| 查看: 5351 | 回复: 28 | ||||
| 【奖励】 本帖被评价13次，作者yahoohoo增加金币 11.5 个 | ||||
| 本帖产生 1 个 模拟EPI ，点击这里进行查看 | ||||
[资源]
【原创】分子模拟实例与挑战专题之二加速能量与力的计算
|
||||
|
分子模拟实例与挑战专题第二期 本帖内容，非注明参考文献出处的，版权归本人所有，请勿以任何形式转载或传播。对于以直接或间接方式使用文中涉及代码而引起的损失，本人概不负责。 在分子模拟中，能量（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人回复
-
一志愿中国石油大学（华东） 本科齐鲁工业大学
已经有3人回复
-
332求调剂
已经有4人回复
-
265求调剂
已经有9人回复
-
328求调剂，英语六级551，有科研经历
已经有8人回复
-
求调剂
已经有3人回复
-
278求调剂
已经有6人回复
-
307求调剂
已经有10人回复
-
279分求调剂 一志愿211
已经有12人回复
-
297求调剂
已经有9人回复
» 本主题相关商家推荐: (我也要在这里推广)
» 本主题相关价值贴推荐，对您同样有帮助:
- 分子模拟计算机配置选择
已经有42人回复
- 关于计算平均能量
已经有3人回复
- 如何计算金属离子的能量
已经有3人回复
- 怎么计算一个分子的能量
已经有4人回复
- 高斯计算能量问题
已经有10人回复
- Latex如何使用
已经有8人回复
- 关于波尔兹曼常数的理解
已经有10人回复
- 有关castep能量计算，请高手赐教！
已经有12人回复
- 如何能计算的更快？
已经有5人回复
- 氧分子、Ti、Bi能量计算
已经有7人回复
- 如何通过能量的Boltzman统计平均计算出分子的性质？
已经有11人回复
- 请教分子力学与分子动力学的区别
已经有14人回复
- SIGMA的单位是什么
已经有7人回复
- 分子能量最小化的软件有哪些？
已经有6人回复
- 吸附时，单独计算小分子能量时的设置
已经有6人回复
- 【求助】分子动力学模拟源代码
已经有11人回复
- 【求助】单个氧原子的能量E是多少？
已经有4人回复
- 【求助】请问分子能量是什么呀？用什么软件计算？
已经有4人回复
- 【求助】分子模拟总能量是否一定小于0？
已经有13人回复
- 【讨论】MC模拟团簇过程中能量升高了
已经有16人回复
★ ★ ★ ★ ★ ★ ★ ★ ★ ★
ghcacj(金币+10):谢谢分享 2010-09-16 08:52:05
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 ] |
3楼2010-09-16 02:37:36
★★★★★ 五星级,优秀推荐
|
做分子模拟时，大家关心的一个问题是何时选用哪种方法 本文来自: 小木虫论坛 http://muchong.com/bbs/viewthread.php?tid=2402779&fpage=1 －－－－－－－－－－－－－－－ 楼主讲的很好，很有条理，学习中！ |
4楼2010-09-16 08:37:44
5楼2010-09-16 14:44:04
6楼2010-09-16 14:51:22
7楼2010-09-16 14:57:33
8楼2010-09-16 15:24:40
9楼2010-09-16 15:26:54
10楼2010-09-16 15:27:46
11楼2010-09-16 15:45:33
12楼2010-09-16 15:53:52
13楼2010-09-16 16:10:43
14楼2010-09-26 21:04:51
15楼2010-09-26 21:09:33
16楼2010-09-26 21:11:56
17楼2010-09-28 15:53:07
18楼2010-11-07 15:24:42
19楼2010-11-08 10:16:26
20楼2010-11-15 18:40:03
把自己的编程经历
|
你可以写下来，给我看看，看过哪些书，哪些开源软件的源代码。。。。 本文来自: 小木虫论坛 http://muchong.com/bbs/viewthread.php?tid=2402779 |
21楼2010-12-10 17:14:14
22楼2011-04-01 20:07:16
23楼2011-04-01 20:46:36
24楼2011-04-01 21:27:47
25楼2011-04-08 14:20:20
26楼2011-04-08 17:13:10
27楼2014-02-14 15:55:43
28楼2014-02-16 00:34:16
29楼2014-06-11 21:53:52
简单回复
zyj81192楼
2010-09-15 16:59
回复
顶一下！
