新手写MD程序的一点心得体会
心得体会, 程序
本帖最后由 chenys 于 2009-6-20 08:12 编辑
前段时间开始着手自己编MD程序，就写点心得体会吧。呃，首先是我的启蒙帖，非常感谢cuiloveyang版主！
1。单位制
宏观的SI单位制不适用于原子尺度的MD，所以选择一套单位很重要，因为我对Gromacs比较熟悉所以就直接采用了它的单位体系，即以nm，ps，kJ/mol，bar等为基础的一套单位。这套单位的好处在于计算加速度时不用乘以额外的常数，力的单位kJ mol^-1 nm^-1，质量单位是原子单位制（即碳12原子的1/12质量），加速度单位nm ps^-2，量纲刚好满足a=f/m。唯一需要注意的就是压力换算的常数，1 bar = 1/16.6054 kJ mol^-1 nm^-3，以及计算库仑力的常数f = 138.9354859 kJ mol^-1 nm e^-2。
2。程序语言
最开始我用fortran，因为用得比较多，后来发现用c语言算起来更快一些就把程序改成了c语言，刚开始的时候犯了很傻的错误就是数组越界，orz。如果不局限于科学计算的话，我觉得c/c++的应用范围更广，所以也希望借着写MD程序来学一点c的知识和技巧。
3。体系自由度
对于三维的MD体系，自由度N_free = 3 * N - N_c - N_COM，其中N是原子个数，N_c是constraints的数目，N_COM是质心自由度，周期性体系N_COM=3，真空中（非周期性）N_COM=6。N_free出现在计算温度的公式里，当然它会影响thermostat，同时也会影响与thermostat相耦合的barostat。
同时在MD的每一步要消除质心的平动，如果是真空体系还要消除质心转动。这部分代码我参考了Tinker源码，消除平动很容易了，而消除转动的代码要复杂一些，我基本照抄过来，未求甚解，惭愧。
4。积分算法和constraint算法
我采用了velocity Verlet积分算法，暂时没有尝试其它算法。与velocity Verlet积分算法相对应的就是RATTLE constraint算法了。RATTLE分为两步，第一步是在t+0.5dt的时候（dt为步长）实施SHAKE算法调整原子坐标和速度，第二步是在t+dt的时候用RATTLE算法调整速度。源码参考Tinker。
5。力和能量的计算
短程力和能量的计算参考Computer simulations of liquids的随书程序，很简单。长程库仑力的计算我采用了Wolf算法，这个算法非常简洁而且有效，cuiloveyang版主做过详细介绍。传统的Ewald和PME算法我大概知道其原理但是数学不过关所以自己编不起来，orz。
6。热浴和控压
传统Nose-Hoover算法。这部分我主要参考了DL_POLY 3的源码，其描述为“velocity Verlet with Nose-Hoover thermostat and barostat (isotropic pressure control) (symplectic)”，参考文献Mol. Phys., 1996, Vol. 87 (5), p. 1117。这里的控压算法大概是Anderson？
参考文献里有一个Nose-Hoover chain的概念，通俗的讲大概可以理解为“热浴的热浴”，Tinker的热浴源码采用了数目为2的Nose-Hoover chain。
7。邻位算法
对于大体系，需要采用Verlet list或者linked-cell list等算法来加速计算。Verlet list算法因为每隔十步左右就需要更新列表，而更新列表的计算量为O(N*N)，感觉效率不够高，我就直接采用了linked-cell list算法。源码参考Computer simulations of liquids的随书程序。
8。并行
使用cuiloveyang版主推荐的openmp。一顿折腾之后才明白并行不是简单的cpu越多越快，而是整个程序的数据处理方式都要有很大改变，要消除数据之间的依赖性。还有个很尴尬的问题是经常碰到并行后反而比串行更慢的情况。我自己改了很多次代码之后终于实现了并行计算力和能量，但是我不知道是如何成功的，因为之前都是并行以后反而更慢，改着改着忽然变快了，我也很orz。现在的并行效率还比较可怜，还需要继续摸索，不过openmp真的很方便，对于shared memory计算机（也就是个人机或小型工作站？）来说是不二选择。
总结
自己写MD程序的一点心得体会，现在能够做rigid SPC/E水的NPT，尚
本文来源于分子模拟论坛 ，原文链接：
阅读(2518) | 评论(0) | 转发(2) |
