分子动力学模拟的优化与并行研究

《分子动力学模拟的优化与并行研究》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

·84·计算机应用研究　2005年分子动力学模拟的优化与并行研究张勤勇1,2,蒋洪川2,刘翠华1(1.西华大学材料学院,四川成都610039;2.电子科技大学微电子与固体电子学院,四川成都610054)摘　要:分析讨论了分子动力学模拟的算法特征和计算特点,对串行程序作了优化,并使之适合于作并行化。对模拟体系使用区域分解的方法,在计算节点间保留了部分重叠区域,采用基于消息传递的MPI设计平台,在可扩展机群上实现了并行化,获得了90%以上的并行效率。关键词:分子动力学模拟;优化;并行中图法分类号:TP30116　　　文献标识码:A　　　文章编号:100123695(2005)0820084202ResearchofOptimizationandParallelizationforMolecularDynamicsSimulationZHANGQin2yong1,2,JIANGHong2chuan2,LIUCui2hua1(1.SchoolofMaterial,XihuaUniversity,ChengduSichuan610039,China;2.SchoolofMicroelectronics&SolidElectronics,UniversityofE2lectronicsScience&Technology,ChengduSichuan610054,China)Abstract:Analyzesthealgorithmiccharacterofmoleculardynamicssimulation,andoptimizestheserialprogram.Inthetask decomposeforparallelization,spatialdecompositionmethodisadopted,andpartofareaiskeepingasoverlapinneighborsub2 area.OnthebasisofthemessagepassingparallelprogrammingplatformMPI,andonthenetworkofworkstation,paralleliza2 tionisimplemented.Intheparallelenvironment,theparallelefficiencyisabove90%.Keyword:MolecularDynamicsSimulation;Optimization;Parallelization　　分子动力学(MD)是一种物质微观领域的模拟方法[1,2],解上述微分方程的方法很多,本文采用Verlet法[2],该方这种方法通过计算机模拟物质体系中微观粒子(主要是原子、法简便易懂,稳定性也比较好。①通过式(2)计算得到全部粒分子)之间的相互作用及运动来得到体系中粒子的运动轨迹,子的受力;②解式(1),得到粒子在该时刻的运动速度v和位再按照统计物理的方法计算得出物质的宏观性能等。自此方置r;③时间步长加1。重复上述过程,直到计算完预定的时间法提出以来,尤其是随着计算机技术的进步,发展十分迅速,逐步长数。通过这个过程,得到了体系内全部粒子的微观运动轨迹,然后运用统计物理的方法可以计算诸如分布函数、扩散系渐成为模拟微观系统、预测宏观性能的主要方法之一。数、各种相关函数、自由能、热传导系数等,从而达到与试验对分子动力学是典型的计算瓶颈问题。①组成物质的微观粒子的数量是十分庞大的,如固体的典型密度为每0135μm的照检验、指导进一步的试验研究、预测物质性能等目的。在立方体拥有109个原子;②为了精确跟踪粒子的运动,保证数MD模拟中,改变物质密度、组分以及温度高低、升降温速度等试验条件是十分方便的,能够完成现实中很多根本无法完成的 值模拟的稳定性,必须使用很小的时间步长。分子动力学中典试验,甚至假想试验。进行系列MD模拟,就相当于进行了系 型的时间步长为飞秒,即使是微秒级的现实模拟也需要10万列现实试验,对研究成本降低、研究周期缩短都是有极大好处 步,显然计算量是十分庞大的。所以分子动力学模拟对计算机的。相比于其他模拟方法,MD模拟得到了系统任意时刻的位 性能的需求可以说是无限的,并行计算是必由之路。形,方便研究人员进行各种研究,尤其重要的是这种位形是系1　串行算法及优化统在内禀动力学驱动下形成的,不具有任何随机性。所以,分111　基本串行算法子动力学模拟越来越得到研究人员的重视,其发展相当迅速。在上述计算过程中,计算量最大的是对式(2)的计算,即系统中每个粒子满足牛顿运动方程:力的计算,约占到全部计算时间的90%左右。为了达到模拟d2rm=F　i=1,2,⋯,N(1)dtri-rjFi=∑fij,　fij=-φ′(|ri-rj|)(2)j≠i|ri-rj|实际宏观系统的目的,一般采用周期性边界条件[1],以一个较小系统在三维方向周期性复制,形成较大的系统来代表实际宏观系统,计算时只需要计算小系统。这样做可以大大减少计算量,但小系统不能太小,否则界面效应将十分明显,从而使偏离 式中,N为要模拟的全部粒子数,Fi为每个粒子受到的系统内实际系统的程度过大。为了尽量使模拟贴近实际试验,则需要 全部相互作用力的总和,fij则为i粒子与j粒子的相互作用力,模拟尽可能多的粒子数,这显然将加大计算量,甚至有时在串 φ(x)为粒子相互作用的势函数。势函数是系统内粒子之间所行情况下根本就无法完成。 有相互作用的一种近似。 112　算法特征分析及优化收稿日期:2004209201;修返日期:2004211228在上述计算过程中,存在如下特征: 第8期张勤勇等:分子动力学模拟的优化与并行研究·　85　·　(1)一般来说,粒子的相互作用中,以短距离作用为主,长要建立一个存储各元胞链头的链表。若采用近邻表法,则要对距离作用可以忽略。在式(2)的计算中,可以认为i粒子只与每个粒子建立链表,用以存储其近邻表。其相距Rc距离内的粒子发生相互作用,从而只计算在Rc距离经过优化后,程序用多种体系进行试验,运行效率提高约内的粒子间的相互作用,称Rc为截断距离。这样一来,计算量在40%～45%,而且模拟系统规模越大,优化效果越明显。可大幅度减少。当然,Rc的选取要靠试验来确定,因为Rc对系统的演化是有影响的,一般以保证能量守恒为条件来确定。2　并行化(2)可以对粒子进行分区,形成小的计算元胞。计算式从设备情况来看,大型的并行机是十分昂贵的,一般科研(2)时,只计算元胞内及相邻元胞间的相互作用,而不再计算单位不具备这个条件。并行可以在共享存储池[4]和分布式与之相隔较远的元胞间的相互作用。此即元胞法[1],如图1[1],如图1池[5,6]两方面应用,文献[7]对并行化作了比较全面的评述。所示。但相对而言,比较可行的方法是用普通PC机作计算工作站,假设在二维情况下,在计算A元胞中粒子的受力时,只需若干台PC机构成机群,采用基于消息传递的MPI计程序,使之具有可扩展性,是既经济又可行的办法。所以本文主要考虑在可扩展的机群上实现并行化。并行化主要是并行方案的选择,选择的依据是负载平衡、通信量大小、同步开销等。从MD算法分析可以知道,并行方案可以采用两种途径:①粒子分解,即将所有参加模拟的粒子编号,平均分配给所有计算节点,这种方法的好处是负载平衡容易保证,同步开销小。但粒子在空间是随机分布的,且在不停运动,为了使计算能够进行,每运算完一个时间步后,需要把所以粒子位置、速度信息进行交换,这样通信量会很大,而且不具有可扩展性。②区域分解,即像元胞法那样,将模拟系统平均分成若干个计算区域,每个区域对应一个计算节点,每个区域内当然还可以采用元胞法、近邻表法等。由于粒子运动毕竟是比较慢的,在一个时间步内,运动出本计算区域的粒子数是比较少的,所以通信量比较小;可以根据模拟系统的大小决定需要计算的节点数,具有良好的可扩展性;但由于每个区域的粒子数不同,负载不容易平衡,同步开销比较大。本文基于MPI消息传递并行程序设计平台,采用区域分解法,用C++语言编程。设计程序时,区域划分、负载平衡、通信量控制等一并考虑,其中主要考虑通信量控制。从计算过程看,通信的主要构成有:计算力时为了完成边界粒子的受力计算,需要其他计算区域内部分粒子的位置信息;经过一个时间步后,本计算区域内的粒子产生了运动,可能有的粒子运动出了本计算区域,需要把这部分粒子的位置、速度信息发给相应的计算节点,同时接收其他计算区域中运动到本区域的粒子的信息。在计算分区边界粒子时,需要邻近分区的粒子的位置信息,一般采用把所有邻近分区粒子信息存储在本节点,这样做既消耗内存,同时还使通信量大大增加。本文在区域分解时按平均分解,但相邻区域间存在一个2×Rc厚度的信息重叠层(),仅包含位置信息通信量减少、内存消耗减小、效率大大提高,如图4所示()仅画出两个分区之间的重叠情况。由图4可以看出,计算节点p包含了计算区p+1中邻近计算区p的厚度为Rc的粒子的信息,而计算节点p+1包含了计算区p中邻近计算区p+1的厚度为Rc的粒子的信息。如此,则计算节点p中包含了要完成区域p所需要的全部粒子信息。整个过程为:当所有节点完成一个时间步后,节点p向其余节点发送运动出本区域的粒子的信息及本区域边界厚度为Rc的所有粒子的信息,然后接收邻近区域发送来的包含运动到本区域的粒子信息和邻近区域厚度为Rc的重叠层的粒子信 设计平台设要计算A元胞中的粒子及与A元胞近邻的B～I元胞中的粒子的相互作用,共需要计算九个元胞。进一步考虑到牛顿第三定律,实际计算时只需要计算五个元胞,如A～E,A元胞内的粒子间受力是必须计算的,而把与B计算的结果加到B中粒子受力结果中,与C计算的结果加到C中的粒子受力结果中,依此类推。当然,在计算F元胞中粒子受力时,同样要把与A计算的结果加到A中粒子受力结果中;计算G～I元胞中粒子受力时,同理都需相应加到A中。这样,既保证了A中粒子受力计算的完整性,又减少了一半的计算量。在三维情况下,可以同样处理,对27个元胞,只需要计算14个元胞。(3)可以考虑将粒子i近邻一定半径R范围内的粒子构i[1],在计算i粒子受力时,只需要查对应造成一个表,即近邻表的近邻表,计算近邻表中的粒子,可以进一步减少计算量。此即近邻表法,如图2所示。当然,由于粒子是在运动的,所以粒子的近邻表需要随时更新。为了减少更新的次数,可以将近邻半径设得大一点,则特定步长H内粒子不会运动出近邻区,从而只需在特定步长H后再更新近邻表,提高效率。受力计算在分子动力学模拟中消耗的时间最多,所以一般只对受力计算作优化。由上述分析还可以看到,分子动力学模拟尤其是受力计算完全可以进行并行化,让多个处理器同时运算,然后汇总,这对提高计算速度提供了良好的先决条件。113　串行优化算法的实现　　从上述分析可以确定串行程序的流程如图3所示。其中,系统驰豫是指系统从不平衡到平衡的过程,通常就是一定步数的MD模拟。由于系统初始化几乎可以任意,此时必然是不平衡的,而且不一定要在模拟的相空间内,经过一定时间后,系统会自动恢复到平衡状态。只有在平衡态下开始投产模拟,才能得到有用的结果。本文采用元胞法。由于元胞中的粒子数不相等,所以用链表结构存储粒子信息是合适的。对每个元胞建立链表,同时还息,更新重叠层信息,更新本区域的邻近表链,(91)下转第页 第8期张　晓等:复制子系统性能评估模型的研究与实现·　91　·　　　试验结果和根据I/O路径长度推测的结论一致。另外从试验结果可以看出,三种备份方式速度分别相差一个数量级,这说明基于SAN的备份方式由于直接在磁盘间进行数据的复制,不需要占用CPU和网络的带宽,缩短了I/O路径,可以获得很好的性能。6　结论5　实验2:不同实现方式的复制性能的比较本文提出了一种备份系统性能分析的模型,利用该模型可对利用不同实现方式实现的复制子系统进行复制性能的比较,备份从I/O路径来说有三种方式,分别是通过网络备份、也可以利用该系统进行不同复制模式下复制性能和应用程序通过主机备份以及使用SAN的备份方式。三种实现方式中I/性能的试验。该系统是进一步进行性能研究的基础,利用该系O路径的示意图如图5～图7所示。从备份路径可以看出三统可对影响复制性能的各个因素进行调整,研究不同因素对于种备份方式的性能应该是:SAN的备份方式>通过主机的备复制性能的影响,包括网络速度(10Mbps/100Mbps/1GbpsIP份方式>通过LAN的备份方式。下面使用我们提出的模型进网络或FC网络),不同文件系统及文件系统分块大小等。行验证。参考文献:[1]RichardBarker,PaulMassiglia.StorageAreaNetworkEssentials[M].2001.2712293.[2]EMCMirrorViewversion1.xAdministrator’sGuide[R].EMC试验环境:操作系统Solaris5.9,复制300MB的数据。通过LAN方式的备份使用FTP来代替。由于FTP本身提供传输速率的数据,网络备份的性能数据采用FTP提供的数据。试验结果如表2所示。Corp.,2002.[3][美]MarcFarley.SAN存储区域网络[M].孙功星,等.北京:机械工业出版社,2001.1662204.[4]张晓,陈建全,郑炜.基于SAN的高可用性存储管理复制子系统设计与实现[J].计算机工程,2004,30(6):70272.[5]BethHamel.IntroductiontoReplicationPerformanceTuning[R].IBM,1997.[6]http://www.sql2server2performance.com[EB/OL].作者简介:表2　试验结果(单位:KB/s)12345平均张晓(19782),男,博士研究生,主要研究领域为存储区域网相关软件LAN2Backup941.841065.921114.351118.541113.331070.796Direct2Copy7475.217419.977507.387671.027428.477500.41SAN2Backup53901.5154058.5354084.2253526.8353849.9853884.21的设计及优化、存储系统相关算法的设计与优化、参与异构系统间共享文件系统、存储管理系统、备份系统的设计与实现;李战怀,教授,博士生导师;贾瑞勇,博士研究生,研究方向为基于网络存储的文件系统;陈建全,博士研究生,研究方向为存储集群。(上接第85页)开始下一个时间步的计算。计算节点程序流程图如图5所示。4　结论本文讨论了分子动力学算法优化与并行,并采用基于MPI程序设计平台实现了分子动力学程序并行化,取得了较好的性能结果。程序PMD具有良好的可扩展性,但由于粒子是运动的,估计对高温情况(粒子运动速度快)会出现负载不平衡,导致性能下降,需要进一步研究。参考文献:3　并行计算数据及分析基于MPI消息传递并行程序设计平台和以上并行算法及其优化技术,我们设计了并行程序,程序名叫做PMD,在并行机P2IICluster上试验Cu薄膜晶化问题,粒子数量为6万个。试验了不同数量的计算节点情况下每时间步消耗的时间,结果如表1所示。表1　程序PMD的并行计算性能(Cu薄膜晶化问题)计算节点数12468并行时间(s/步)14.877.924.082.582.06加速比1.001.883.65.767.22并行效率1.000.930.910.960.90　　从表1可以看出,并行效果是明显的,并行效率均保持在90%以上。预计随着粒子数量增加,效果会更明显。[1]AllenMP.ComputerSimulationofLiquids[M].Oxford:ClarendonPress,1987.17235.[2]VerletI.Computer"Experiments"onClassicalFluids,theThermody2namicalPropertiesofLennard2JonesMolecules[J].PhysicalReview,1967,159(1):982103.[3]PhlimptonS.FastPrarllelAlgorithmsforShort2RangeMolecularDy2namics[J].JournalofComputationalPhysics,1995,117(5):125.[4]董科军,刘让苏,郑采星,等.液态金属原子大系统凝固过程模拟的并行算法研究[J].自然科学进展,2001,11(2):1872193.[5]莫则尧,张景琳,等.二维分子动力学(MDP)的并行与优化[J].计算物理,2000,17(122):1932195.[6]朱宇翔,张景琳.在可扩展机群系统上二维分子动力学问题的并行计算[J].计算物理,1999,16(4):4222426.[7]PriyaV,AiichiroN,etal.MolecularDynamicsSimulationofCovalentAmorphousInsulatorsonParallelComputers[J].JournalofNon2crys2tallineSolids,1995,182(3):59262.作者简介:张勤勇(19722),男,四川仁寿人,讲师,硕士,研究方向为纳米材料的第一原理与分子动力学模拟;蒋洪川(19712),男,四川成都人,讲师,硕士,研究方向为磁性材料及器件、MEMS;刘翠华(19692),男,四川资中人,讲师,学士,研究方向为材料的分子动力学模拟。