原文链接:
看完《分子模拟从算法到应用》那本书的第四章,不用全看完,而且起码要对分子动力学模拟过程有一个了解。试着根据书的过程做个Ar的NVE,虽然Ar和离子晶体以及其它的任何材料的差异仅仅是势函数的问题,即使由势函数带来了一些问题,并且这种都不是本质问题。从初始化的原子数,原子位置,初始速率,时间步长,初始体温等等这种初始化结束了之后,选择一个简单的积分算法,如6阶的Gear预测校准,不要控控温压,就是一个简单的NVE,不要考虑任何的提升效率的邻位算法,由于这个时侯我们可以选择5×5×5的超原胞,总共的原子数也就500个,不须要考虑邻位算法。开始循环估算:预测----估算原子的力和能量—校正—输出能量。
这样最简单的NVE就编成了,总共也就1000多行,是个很小的程序。自己先试着体会一下。虽然当这样的小程序完成之后,你会感觉分子动力学编程也很简单,这么接出来的复杂的分子动力学也不会是哪些问题。
做完了这一些,你须要晓得的是这些是和材料无关的东西,这么就尽量的分离,开始使用一个个的函数。诸如,原子的位置是和具体的材料相关,并且初始速率却和材料一点关系都没有,同样的数值积分中的预测和校准也是和材料无关的,之后的温控和控压算法也是和材料无关的。当规模大了之后,邻位算法也是和材料无关的,像这种和材料无关的部份最好自己弄成小函数,选择调用。之后换材料的时侯程序也不会有太大的改动。
编程小方法
1.选择用intel编译器,个人喜欢用10.1或则9.1的版本,打开优化选择,类似的/QaxS/QxS/Qipo/Qprec-div-等等,之后可以使用Openmp的并行估算(具体的可以参考intel编译器的帮助指南)
2.尽量的简化估算,比如2×a就要写成a+a,在计算机中,加减是一个数目的估算,乘是一个,除是一个。估算量逐步降低,所以选择a+a来取代2×a会减轻一些估算量
3.选择数据来取代结构体,结构体看的比较便捷,而且估算效率要低,尤其是在编译过程中的矢量化的时侯,而链表则可以挺好的矢量化,也更适宜并行。结构体的估算是先找到结构体表针,之后再找上面的参数,当并行的时侯,多线程同时找结构体表针,会很大的减少速率。
4.乘法尽量用加法来代替,有了程序之后,自己可以仔细的剖析上面的估算消耗,这样可以更好的优化估算,其实最主要的消耗是在力和能量的估算中,可以选择离散-配准的方式,促使不同的势函数有相同的估算效率,也可以选择邻位算法来分块估算。
NVT程序
当简单的NVE做完之后,可以直接做NVT,这仅仅是降低一个简单的温控,是很简单的,初期编程的时侯,建议选择Nose-Hoover温控,你会发觉NVT也是这么的简单。接出来,你可以考虑做做NPH,然而因为引入了浮力算法,原先的原子的位置和速率等等和笛卡尔系有关的一切内容都要发生改变,这将让你重新写程序,可以说是一个完全不同的程序。不过,幸而我们早已有了上面的一些经验,NPH即使复杂,而且并不是不可能完成的,注意的是和笛卡尔系相关的量要变化。初期的估算,可以选择Anderson控压算法。在完成了这种,可以说你早已把最主要的分子动力学程序都完成了。
分子动力学算法
这么,接出来就是复杂的分子动力学算法了,试着选择PR算法来替代Anderson浮力算法,之后选择Metric-tensor来取代PR算法;选择Nose-Poincare来取代Nose-Hoover算法,选择Generalized-Leap-Frog算法来代替预测校准算法,用Wolf来代替Ewalds算法,这种一步步的改进就会让你有好多新的发觉。你会感觉你的程序达到了如今大部份软件包没有的功能。随着规模的减小,你可以选择Verlet列表,原胞列表、结合法以及快速排序来实现邻位算法,假如再有了openmp并行估算的加入,你会发觉你的程序早已有了一个质的突破。
END
当那些都完成了之后,接出来就可以改变势函数来模拟不同的材料了,但是你对分子动力学的理解也会上升到一个新的高度,再者,其余相关的结果也不会是问题了,比如,弹性挠度的估算,MSD,径向分布函数等等。你可以直接的加入代码来实现你想要的过程。