（机械设计及理论专业论文）并行分子动力学优化与应用及固体热整流现象探讨.pdf

摘要 分子动力学并行算法的优化与应用以及固体热整流 现象探讨 摘要 姓名刘正华指导老师陈云飞教授 东南大学机械工程学院 为提高分子动力学模拟的计算效率，采用并行计算代替单机计算。作为并行算法之一的原子分 解法，由丁- 其简单性，在小型机群中较高的并行效率，和一些具体问题的适 j 性如k = 程力的计算。 因此比较常用。本文研究了原于分解方法，优化了备进程粒于的受力计算，尤其对严重影响并行效 率的全局通信进行了优化，使用最佳实现全局通信的点对点非阻塞通信替代了直接调用全收集语句 m p ia l l g a t h e r 。在进行非阻塞通信同时计算各进程粒子与本地进程粒子间及上次通信获得的其他 进程粒子问的相互作用，最终使得并行效率的提高在3 0 以上。并且使用该优化算法对同态氟的 导热系数进行了计算，模拟结果和实验值比较吻合，而且时间大大缩短。证明了该优化并行算法正 确性和有效性。此外对于大规模的粒子系统原于分解已不能满足要求，这种情况宜采用区域分解 法替代之。本文以碳纳米管的计算机模拟计算为例，对其采用区域分解算法实现并行化，结果证明 了随着粒子数量的增加并行效率是在上升的。且证实了该算法在碳纳米管理论研究中的重要地位。 本文中，使用分子动力学模拟方法对双层纳米薄膜界面处的热整流进行了计算。该双层纳米薄 膜由两种不同材料组成。模拟结果表明，界面热阻与流过界面的热流方向有关。界面处热阻的热整 流随着温度的上升而减小，但是整流效果不能仅仅由界面处的温度差解释，界面处声子的非弹性散 射可能对这种非对称传播有一定贡献。最后对于界面处的相互扩散对热整流的影响也作了一定的讨 论。通过将界面附近两薄膜层的任意原子相互交换实现界面的掺杂，分子动力学模拟结果表明了界 面处相互扩散将会降低整流效果。 关键词：分子动力学，并行算法，原子分解法，m p i 域分解法，热整流 东南大学硕士论文 t h eo p t i m i z a t i o na n da p p l i c a t i o no fp a r a l l e lm o l e c u l a rd y n a m i c s s i m u l a t i o na n dt h es t u d yo ft h e r m a lr e c t i f i c a t i o ni ns o l i d s a b s t r a e t z h e n g h u al i u ，j u e k u a ny a n g , a n dy u n f e ic h e n ( d e p a r t m e n to f m e c h a n i c a le n g i n e e r i n g , s o u t h e a s tu n i v e r s i t y , n a n j i n g2 1 0 0 9 6 ) a b s t r a e t ：t oe n h a n c et h ec a t e u l a t i o ne 币c i e n c yo fm o l e c u l a rd y n a m i c ss i m u l a t i o n w eu s ep a r a l l e l a l g o r i t h mt or e p l a c ec a l c u l a t i n go ns i n g l ep c a so n eo ft h ea l g o r i t h m so fp a r a l l e lm o l e c u l a rd 3 ，n a m i c s s i m u l a t i o n a t o m d e c o m p o s i t i o na l g o r i t h mi so f t e na d o p t e db e c a u s eo fi t ss i m p l i c i t ya n dh i g hp a r a l l e l e 币e n c yj ns m a l l s c a l ec l u s t e ra n da p p l i c a b i l i t yt os o m ep r o b l e m ss u c ha st h ec a l c u l a t i o no fl o n g - r a n g f o r c e s w eh a v es t u d y e da t o m d e c o m p o s i t i o na l g o r i t h mi nt h i sp a p e o p t i m i z e dt h ef o r o ec a l c u l a t i o no f 恍r yp r o e e s s e ra n dg l o b a lc o m m u n i c a t i o n w h i c ha f f e c t st h e p a r a l l e le f f i c i e n c ys e r i o u s l y u s e p o i n t t o p o i n tn o n b l o c k i n gc o m m u n i c a t i o nt h a tr e a l i z e st h eg l o b a lc o m m u n i c a t i o nb e s tt or e p l a c e e m p l o y i n g “m p la n g a t h e r o fg a t h e r - t o a l if u n o t i o nd i r e c t l y w h i l eb e i n gc o m m u n i c a t i n g t h ep a r t i c l e s o fe a c hp r o c e s sw i l lb ec a t c u l a t e dt h e i ri n t e r a c t i o n sw i t ht h e m s e l v e sa n dt h ep a r t i c l e st h a th a db e e n r e c e i v e da tl a s ts t e p w h i c hw i l le n h a n c et h ep a r a l l e le f f i c i e n c yt om o r ei b a n3 0 w ec a l c u l a t et h e t h e r m a lc o n d a c t i v i t yo f s o l i da r g o nb yt h eu s eo f t h eo p t i m i z e da l g o r i t h m t h es i m u l a t i o nr e s u l t sa r ev e r y c o n s i s t e n tw i t ht h ee x p e r i m e n t a lo n e s a n dt h es i m u l a t i o nt i m eh a sr e d u c e dm u c h i ti sp r o v e dt h a tt h e o p t i m i z e dp a r a l l e la l g o r i t h mi sc o r r e c ta n de f f i c i e n t i na d d i t i o n a t o m - d e c o m p o s i t i o na l g o r i t h mw i l lb e n o te f l i c i e n te n o u g hi n l a r g e - s c a l ea t o ms y s t e m s w h i c hs h o u l db er e p l a c e db vd o m a i n - d e c o m p o s i t i o n a l g o r i t h m 1 nt h i sp a p e lw et a k e t h es i m u l a t i o no fc a r b o nn a n o t u b e sf o re x a m p l ea n de m p l o y d o m a i n d e c o m p o s i t i o na l g o r i t h mt or e a l i z ei t sp a r a l l i z a t i o n t h er e s u l t ss h o wt h a tt h ep a r a l l e le f f i c i e n e y i so i lt h er i s ea st h en u m b e ro fp a r t i c l e si n c r e a s e s a n da l s oc o n f i r m e dt h a tt h ea l g o r i t h mi nt h et h e o r e t i c a l s t u d yo f c a r b o nn a n o t u b e si si na ni m p o r t a n tp o s i t i o n i nt h i sp a p e r , m ds i m u l a t i o ni su s e dt oe v a l u a t et h e r r n a lr e c t i f i c a t i o na tt h ei n t e f f a c eo f d o u b l e - l a y e r e dn a n o f i l m t h ed o u b l e 1 a y e r e dn a n o f i l mi sc o m p o s e do ft w od i s s i m i l a rm a t e r i a l s 0 u r s i m u l a t i o nr e s u l t si n d i c a t et h a tt h ei n t e r f a c i a lt h e r m a lr e s i s t a n c ed e p e n d so nt h ed i r e c t i o no fh e a tf l o w a c r o s st h ei n t e r f a c e t h er e c t i f i c a t i o no fi n t e r f a c i a lt h e r m a lr e s i s t a n c ed e c r e a s e sw i t hi n c r e a s i n g t e m p e r a t u r e b u tt h er e c t i f y i n ge f f e c t sc a nn o tb eo n l yi n t e r p r e t e db yt e m p e r a t u r ed i f i e r e n c ea tt h e i n t e f f a c e t h ej n e l a s t i cs c a t t e r i n go fp h o n o n sa ti n t e r f a c em a y b ec o n t r i b u t e st ot h i sa s y m m e t r i ct r a n s p o r t t h ee f f e c to fi r i t e r d i f i t i s i o na ti n t e r f a c eo nt h et h e r m a lr e c t i f i c a t i o na l s oh a sb e e ns t u d i e d t h es t a g g e r e d i n t e f f a e ei sr e a l i z e db yr a n d o m l yi n t e r c h a n g i n ga t o m so ft h et w ol a y e r si ni n t e r f a c i a lr e g i o n n l em d s i r e u l a t i o nr e s u l t sd e m o n s t r a t et h a ti n t e r d i f r u s i o na ti n t e r f a c ew i l ld e c r e a s et h er e c t i f y i n ge r i e c t s k e y w o r d s ：m o l e c u l a rd y n a m i c s ，p a r a l l e la l g o r i t h m ，a t o m - d e c o m p o s i t i o n ，m p l l l g a t h e r , n e i g h b o r e x c h a n g ea l g o r i t h m ，n o n b l o c k i n gc o m m u n i c a t i o n ，d o m a i n - d e c o m p o s i t i o n ，t h e r m a lr e c t i f i c a t i o n 东南大学学位论文独创性声明 本人声明所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的研究 成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包古其他人已经 发表或撰写过的研究成果，也不包含为获得东南大学或其它教育机构的学位或证书而 使用过的材料。与我一同工作的同志对本研究所做的任何贡献均己在论文中作了明确 的说明并表示了谢意。 研究生签名：盟日期：野! ! i 东南大学学位论文使用授权声明 东南大学、中国科学技术信息研究所、国家图书馆有权保留本人所送交学位论文 的复印件和电子文档，可以采用影印、缩印或其他复制手段保存论文。本人电子文档 的内容和纸质论文的内容相一致。除在保密期内的保密论文外，允许论文被查阅和借 阅，可以公布( 包括刊登) 论文的全部或部分内容。论文的公布( 包括刊登) 授权东 南大学研究生院办理。 研究生签名：兰l 兰芏导师签名： 牲日期：迸旦 第一审绪论 第一部分分子动力学并行算法 第一章绪论 1 1 纳米科技现状 纳米科学技术m a n o s c i e n c ea n dt e c h n o l o g y , n a n o - s t ) 是指在纳米尺度( 1n m 1 0 0 n m ) 上研究物 质( 包括原子、分子的操纵) 的特性和相互作用，以及利用这些特性的多学科交义的科学和技术， 它是当前国际上最活跃的研究领域之一。 纳米技术的基本概念是1 9 5 9 年由美国著名的物理学家、诺贝尔奖获得者r i c h a r df e y m a n 1 首先 提出的，他在美国物理协会年度会议上预言：如果我们能按照白己的意志去安排一个个原子，并把 这种排列控制在很小的范用内，我们将会获得具有奇特性质的物质。由此便揭开了纳米科技研究的 序幕。 1 9 8 1 年美国麻省理t 学院d r e x l e r i ”发表了第一篇关于分子纳米技术的论文，提出人造分子机器 的构想。1 9 8 2 年b j n n i n z 等首先成功研制出扫描隧道显微镜( s c a n n i n g t u n n e l i n gm i c r o s c o p y , s t m ) ， 为科学家们在原子尺度上研究物质结构开辟了新途径。1 9 8 4 年德国萨尔兰人学( u n i v e r s i t yo f s a a r l a n d e s ) g l e i t e r 等1 5 1 首次成功制备y p d 、c u 和f e 金属纳米粉末，并通过原付成形烧结得到金属 纳米多晶块体材料，开创了人类纳米材料制各和研究的历史。1 9 9 0 年i b m 公司的e i g l e r 等”1 利用扫 描隧道显微镜直接操纵原子，在金属镍板上按自己的意志安排原子组合成“i b m ”字样，如图1 1 所 示，这是人类历史上第一次直接操纵原子。1 9 9 0 年7 月第一届国际纳米科学技术会议在美国的巴 尔的摩与第五届国际扫描显微镜会议同时举办 n a n o t e c h n o l o g y ) 和 n a n o b i o l o g y 两种国际专 业期刊也在同年相继问世。自此，- - f l 崭新的科学技术纳米科技，经过3 0 年的酝酿，终于诞生 了。 图1 1 原子字样“i b m ” 图1 2 用微犁生物马达组成的纳米器 近年来自然科学和工程技术发展的一个重要趋势就是朝微型化迈进，其中微电子机械系统 f m i c r o - e l e c t r o m e c h a n i c a ls y s t e m s - - m e m s ) 尤其取得了巨大成功，并正被拓展应用于各种工业过程。 这类系统指的是那些特征尺寸在l m m 以下但x 人于1 岬的器件，它集电子及机械组件于一身，并 且要通过集成电路制造中所采用的批簧加工方法制成。纳米器什则进一步推进了微电子机械系统的 小型化。由丁现代制造与应川技术的持续进展，“微机械”或“纳机槭”的尺寸正以超乎寻常的速度降 低，而同时其性能却得到了保持甚至更好取得了诸多令人惊讶的成就。随着微机械加工技术的快 速发展，一些真正分子水平上的机器，如转子、齿轮、开关、闸门、转栅、马达等的制作甚至已经 成为可能【7 】。图i - 2 便是用微型生物马达组成的纳米器件。 纳米技术主要是研究纳米尺度( 从0 1 1 0 0 n m ) 物质的特性和相互作用，以及利用这些特性 设计和制造出具有独特功能、尺寸微小的产品。冈此，纳米科学技术使人类认识和改造客观物质世 界的能力和手段延伸剑原子或分子水平。这不但改变了人们的思维模式，而且将对科学技术与国民 经济的发展产生深远的影响。 1 2 分子动力学仿真技术现状 由于纳米科技的发展，众多机械电子等产品的微观化愈演愈烈，一些在理论分析和实验观察 上难以了解的现象就会早现山来，这就需要进行物质宏观特性和微观结构的综合考虑。分子动力学 仿真技术i 目”i 就是在揭示各种现象本质和探索新规律的微观研究中的主要科学手段，在物理、化 l 东南夫学硕j ：论文 学、生物学、材料学等诸多领域都得i j 1 泛应用，被称为“计算机实验手段”，是一种新型的计算机 模拟技术。 但是，分子动力学方法属于典型的计算瓶颈问题：( 1 ) 空间尺度小、粒子数晕大：( 2 ) 时间 步疑小，典型时间步长在秒级。任何物理现象的宏观特性都是由组成该系统的人彗原子、分子或 者电子等的运动状态所决定的。真实的物理系统所包含的微观粒子数域往往高达1 0 ”数量级，这 比_ u j 现有的电子计算机能够进行模拟计算的粒子数量人得多。模拟系统粒子数量上的巨大著别，将 导致计算结果的偏著。另外，由于分子动力学模拟的是分子级别的动力学情况，分子运动过程的相 互影响是极端复杂的，并且分子本身的状态参数也是变化不定的，所以需要采样时间间隔很小，典 型为e 秒级。 采州分子动力学进行计算机仿真时，模拟系统的粒子数越多，结果将越接近真实值。这同时 就对计算机的性能提出了更高的要求。 1 3 并行计算技术， 随着信息时代的到来，需要处理的信息量越来越庞大，需要解决的问题越来越复杂，这些都 使的计算颦剧增。通过提高单个处理器的计算速度和采用传统的串行计算技术已难以胜任。因此， 需要有功能更加强大的计算机系统和计算机技术来支撑。并行计算机及并行计算技术应运而生了。 随着并行计算机及并行计算技术的出现，大人提高了计算的速度，突破了串行计算机的物理 极限的限制，更为复杂的问题得到了解决。并且随着网络高速化及低成本集群系统的出现，使得一 般实验室拥有属于自己计算的并行计算机成为了可能。使得并行计算技术得到了普及。 目前，全球研制山的并行计算机亦即超级计算机，已经拥有相当惊人的计算速度和数据处理 能力。今年最新出炉的全球超级计算机5 0 0 强名单中，i b m 公司和美国能源部( d e p a r t m e n to f e n e r g y ) 共同开发的i b m 蓝色基因，l 如| 兰| 1 3 所示，以峰值为每秒4 7 8 万亿次浮点运算的运算速 度成为世界上运算速度最快的超级计算机系统。这个速度比普通的台式电脑速度快数万倍! 图1 - 3i b m 蓝色基因l 并行计算已经给许多研究领域带来了新的活力，且在各方面如材料、化工、物理、环境和能 源科学等领域得到了广泛的应用1 1 ”】推动了各学科的发展。此外，并行计算能够帮助人们求解 挑战性的问题。例如天气预报问题，这是在过去很好能够较为准确同答的科学计算问题。过去一段 时间取得的一项重要成果就是由n s f 的暴风雨分析与预报科学技术中心开发的地区性高级预报系 统( a r p s ：a d v a n c e dr e g i o n a lp r e d i c t i o ns y s t e m ) 1 1 3 1 , 已经能够对极端恶劣的天气进行准确预报。 在分子动力学并行算法研究方面，随着并行计算机的发展，加上结合实际模型的高敛的并行算 法，在大规模的粒子系统中并行效率达到了非常理想的效果1 1 4 1 1 15 。这其中。由曹小林，莫则尧等 开发的一种基于“块单元”数据结构的可扩展并行算法，以实现大规模、非均匀分于动力学模拟， 2 第一章绪论 采用空间填充曲线将三维区域问题转换为一维负载平衡问题，然后用基丁| 实测的多层均权法求解， 以保持处理机间负载均衡。结果在一个m p p 并行机的5 0 0 个c p u 上，模拟包含个粒子的三维金属 微喷射模型，该算法获得了4 2 0 倍的加速比。 1 4 本文研究的内容 本文主要针对当前采用分子动力学方法作为手段的计算机仿真来研究物体宏观性质时遇到的 问题，采用较为流行的分子动力学并行算法，对分子动力学并行算法中比较常_ i j 的原子分解法和区 域分解法进行了讨论，并结合实际的理论研究模艰以及并行环境，对原子分解法进行了一定的优化 并应用在了崮态氩导热系数计算中，对碳纳米管的计算应用了区域分解法。结果都得到了非常理想 的效果。 首先第二章介绍了计算机模拟方法之一的分子动力学方法的实现过程和关键技术：在第三章 介绍了并行计算机硬件以及软件的发展和现状；第四章结合两者讨论了分子动力学并行算法的主要 特点和具体实现流程；第五章针对小型的模拟对象( 粒子总数在一万以内) ，采用晟易实现的原子 分解法，并进行一定的优化措施，最终获得了较理想的并行效率；第六章针对碳纳米管的计算机模 拟，采刚区域分解法测试结果证实了该算法的高效性。 由于采_ j 分子动力学方法进行计算机模拟时，结果的真实性与模拟的粒子数是成正比的，但 同时计算时间也在相麻增加，规模人者可以周甚至月来计算。此时，应_ l j 合适并行算法，可以大大 减少计算时间，有助于开展更大规模的计算机模拟以研究诸多新物质的宏观性质。 东南大学顾i ：论文 2 1 引言 第二章分子动力学模拟方法 从物质的内在微观结构及其原子尺度成分出发米预测材料的宏观性能，并按预期目标设计材料 是材料科学的终板目标。计算材料学( c o m p u t a t i o n a m a t e r i a l s s c i e n c e ) 是实现这一目标的有效工具。 随着计算机的软硬件和计算技术的k 足发展，利用计算机对物质的微观结构和运动进行数值模 拟计算得到了迅遮发展，并在此基础上发展了川数值运算统计求和的方法：分子动力学( m d ) 方法p 1 和蒙特卡罗( m c ) il 6 1 方法。这两种方法可以按照原子( 或分子) 的排列和运动的模拟以实现宏观物理量 的计算，不仅可咀直接模拟许多物质的宏观凝聚特性，得出与实验结果相符合或可比拟的计算结果， 而且可以提供微观结构、粒子运动以及它们和物质宏观性质关系的明确图像，有利于从中提取新的 概念和理论，这两种方法又可视为“计算机实验”，利用粒子之间的真实相互作片i 势也可用半经验 势计算凝聚态物质的结构和热力学性质，从而很方便的对相互作用与宏观性质之间的关系进行考察， 据此提出有关物理现象的理论。从而，这种“计算机实验”建立了理论与实验的桥梁对新的理论和 新的实验的建立将起愈来愈人的作用。本章将详细介绍分子动力学( m o l e c u l a r d y n a m i c s ，以下记为 m d ) 基本方法。 2 2 分子动力学方法简介 1 9 5 7 年，a l d e r 和w a i n w r i g h t t ”l 首先在硬球模型下，采用经典m d 方法研究气体和液体的状态方 程，取得了尝试性的成功，其开g l t 用分子动力学模拟方法研究物质宏观性质的先例。1 9 7 2 年。l c c s 和e d w a r d s i ”等人发展了该方法并扩展到了存在速度梯度的非平衡态系统，进一步扩展了m d 方法的 应用范围。1 9 8 3 年，g i l l a n i ”i 等人将该方法推广到具有温度梯度的分平衡态系统，从而形成了非平衡 系统分子动力学体系。一直到今天，m d 方法一直在发展着，并且越来越成熟，也越来越能真实地反 映一些常见材料的宏观性能。 m d 方法的基本思想是通过原子间相互作用势，求出每一个原子所受到的力，在选定的时间步长， 边界条件、初始位置和初始速度下，对有限数目的分子( 原子) 建立其牛顿动力学方程组，用数值 方法求解，得到这些原子经典运动轨迹和速度，然后对足够长时间的结果求统计平均，从而得到多 需要的宏观物理量和力学量。 m d 方法的山发点是物理系统的确定的微观描述，也就是说它是确定性方法。而蒙特卡罗方法是 依赣于物理的概率模拟方法，它的实现依赖于均匀随机数的产生，是随机方法。 分子动力学方法的两个基本假设： n所有粒子的运动都遵循经典牛顿运动定律 2 1粒子间的相互作用满足叠加原理 这意味着m d 虽然是在原子层次研究问题，但它忽略了量子效应，仍然是一种近似计算模型。 2 3m d 模拟的初始条件 模拟时首先要初始化系统的位型分布，即首先要给定微元中分子的初始位置和初始速度。合理 的初始位置和初始速度将使得计算系统迅速松弛到所要求的平衡状态，减少计算次数分子初始位 置晟好与实际情况相类似，常用的分布有简立方晶格分布、面心立方品格分布、体心立方晶格分布 和金刚石品格分布等i ”i 对于各个原子的速度分布可按照m “w e l l b o l t z m a n n 分布来随机选取。 m a x w e l l 分布形式如下i “i ： 4 第一二章分子动力学模拟方法 川= 面m i 一1 1 2p 【三磊】 ( 2 1 ) 这个公式提供了一个质量为m ，的原子f ，温度t 时，在x 方向的速度的波尔兹曼分布是一种 高斯分布，它的实现可以利用计算机上的随机数产生器产生，具体关系见图2 1 。 p ( v 。 2 4 原子间势函数 图2 1 速度分布图 势能模型的建立是进行m d 模拟最为重要的一个环节。势能模型是对体系分子之间相互作用势 的反映，模拟能否成功取决于能否准确地选择势能模型。由于研究范围的广泛性和研究对象的复杂 性，分子间相互作片j 模型也不完全相同。 常用势能模型有对势：l e o n a r d - j o n e s ( l j ) 势1 8 1 、m o r s e 势口“、j o h n s o n 势”“等模型：多体势： 嵌入原子法( e a m ) 1 2 4 1 、s t i l l i n g e r _ w 曲e “s w ) 势”“( s w 作用势为两体一三体作用势一般用于描述 s i 、g e 等原子间的作用力) 、t e r s o f r 势i “1 ( 一般用来描述碳管原于间的作用力) 等模型。 描述含有n 个原子作爿j 的普遍的作用势包括单体、两体、三体和多体，其基本形式如下 2 7 1 ： 九( 1 2 ，) = v ，( ) + v ：( ，o ) + v s ( r , ，r j ，r d + ，j f，j t ，) ( 2 2 ) + v ( ，0 ，) 此处，f 代表第i 个原子的位置，”代表n 体作用势。在实际中，当n 增大时，”随着n 的增 大而迅速趋向于零。一阶势能项一般用来描述固壁情况或者外加作用力。在许多情况下，这一项是 不存在的。第二项代表两个粒子问的相互作用。它是最简单的可能模型。最出名的两体作用势是u 作用势 咖( ) ：拈【( 三) ”一( 旦) 6 】 ( 2 3 ) r hr h 其中为粒子间的间距，6 为l j 势阱常数，口为l j 平衡常数。l j 作用势是短距离作用关系， 它经常用来描述液体或气体间的相互作用。 2 5 分子运动方程的建立 东南人学硕i ：论文 m d 方法的出发点是对物理系统确定的微观描述。这种描述可以是哈密顿描述或拉格c ! j j 日描述， 也可以直接用牛顿运动方程来描述。每种描述都将给出组运动方程，运动方程的具体形式由分子 间的相互作片j 势，即势能模型所确定。m d 方法的具体做法就是在计算机上求分子运动方程的数值 解。采_ j 适当的格式对方程进行近似，即以离散替代连续、以差分替代微分，建立一个有限著分格 式对该差分格式方程纽进行求解。可以在相空间中生成一条路径，沿这条路径可计算出所期望的 体系的各种性质。 运动方稃的迭代方法及离散方法很多，其中以v e r l e t 算法、g e a rp r e d i c t o r - c o t r e c t o r 算法”最 为常见。算法的选择极为重要，一个好的m d 程序需要一种好的对运动方程积分的算法。一般来说， 好的积分算法应当满足以下几个方面： 1 )为了使得模拟的时间尽量的少，算法应当能够快速，同时需要的内存量要比较少： 2 ) 准确度。对较大的时间步长来说更重要，冈为所能使用的时间步长越陈，单位模拟时间内 力的计算量越少。因此，这意味着采用允许长时间步长的成熟算法是有利的： 3 )能量守恒是一个重要的判据，但时间上需要区分两种能量守恒，即短时间和长时间的。辅 助的高阶算法通常在短时间内有很好的能量守恒性，然而该方法通常会有所不期望的特征 即时间较长总能量漂移； 4 1 牛顿运动方程是时间可逆的，算法也应当如此。然而由于种种原因( 如由于运算精度的限 制) ，我们会发现在数值实现时几乎不可能是真正时间可逆的。 5 1算法应当简单而且便于编程。 本文中均采用v e r l e t 算法的改进形式速度v e r i e r 算法( v e l o c i t y - v e r i e r a l g o r i t h m ) “”对运动 方程进行积分，计算过程中原于的受力和运动关系如下： f ( ，) = 坍口 ( f + f ) = r a t ) + v 。( r ) f + ( 1 2 ) a ，( f ) 址2 v ( ，+ a t l 2 ) = v ( f ) + ( i 2 ) a ，( f ) f ( 2 4 ) ( 2 5 ) ( 2 6 ) h o + ，) = v 。o + a t 2 ) + ( 1 2 ) a ，( r + a t ) a t ( 2 7 ) 事实上，v e n e t 算法早在1 9 6 7 年即被提出，此后出现了v e r l e t 算法的一些变化形式。利用v e r i e t 算法，可以在笛卡儿坐标系中计算第i 个原子在t 时刻的位置r i 和速度t 。由泰勒( t a y l o r ) 展开式得到t 时刻附近的位移公式为： ( | + 出) = i ( f ) + t ( ，) 出+ 三( f ) f 2 + 刍巧( r ) f 3 + ( 2 8 ) ( f 一f ) = ( ，) 一t ( f ) 缸+ 吾( ，) ，2 一去耳( ，) f 3 + ( 2 9 ) 将式( 2 8 ) 和式( 2 9 ) 相加，可得到第i 个原子的位置作为其加速度函数的表达式： ( ，+ r ) = 2 ( f ) 一i o 一出) + t ( f ) 出2 + 去( f ) f 4 + ( 2 1 0 ) * 2 ( f ) 一o 一血) + e ( ，) f 2 根据保守力f i 、原子质量m 。，可计算出所要求的第i 个原子的加速度。作用力可由作用势的导 数求出。由式( 2 8 ) 减去式( 2 9 ) 可得原子速度的表达式为： 6 墨三童坌王垡当兰堡丝塑婆 v a t ) z r , ( t + a 1 t ) - i r , ( t - t ) ( 2 ，1 1 ) 从上面的计算式可看出，速度计算的截断误筹是0 ( a t2 ) 而位置计算的截断误差是o ( a t4 ) 。为 克服这一困难，发展了一些v er l e t 算法的变种，其中之一就是著名g l e a p f r o g ) 算法，其公式如 下： t ：( t + a t 2 ) = i ( ，- a t 2 ) + f ( t ) a t ，( f + f ) = ，( f ) + i ( ，+ ，2 ) 7 ( 2 1 2 ) e ( t ) = 去p ( ，一，2 ) + i o + f 2 ) 】 但蛙跳算法也存在一个缺点，就是计算过程中坐标、速度的计算不同步。给程序计算带来不便。 而v e r i e r 算法的速度形式v e l o c i t y v e r l e t 算法解决了这一问题，计算过程中坐标、速度、力( 加速度) 可 同步获得。另外，相比其它积分方法速度v e r l e t 算法使编程实现简单直接，计算速度快，而且所需 的存储空间小，因此在目前的模拟程序中，该方法被广泛应用。本文也是采用速度v e r l e t 算法对运动 方程进行积分的。 2 6 周期性边界条件 正确处理边界和边界效应对模拟方法是至关重要的，因为它是从模拟相对较少的粒子来计算物 质的宏观性质的。为了减小有限尺寸的影响，在分子动力学模拟中采用了各种有效的边界条件。通 常用于分子动力学模拟中的边界条件有三种：周期性边界条件、对称边界条件和同擘边界条件”。 物体的宏观性质是大量粒子的统计行为，模拟系统的粒子数必须非常大才能准确地再现系统的 行为。例如，边长为l u 的立方体的水就含约2 1 0 ”个水分子。另一方面，现有的计算机能力只能 处理的粒子数有限。这样少的粒子组成的集团和真实体系相比，其表面原子与体内原子数之比显然 过大，这必然造成表面效戍。为了解决上述矛盾，采用周期性边界条件的假定使中心元胞在三维 空间上重复排列，于是系统粒子i 的象粒子将在三维空间周州性出现。对于二维的例子，每个包有8 个近邻，而对于三维有2 6 个最近邻包。象粒子的坐标可以通过加上原胞边长的正整数倍得到。如果 一个粒子在模拟中离开这个包，就等于它的象从反方向进入这个包。应用周期性边界条件后，整个 体系变得粒子数无穷大了。 粒子象粒子坐标为 图2 - 2 周期性边界条件示意图 7 东南人学颂i ：论文 n = h l = ( ) 【i4 - m l ，y l + n l ，z i + k l ) ( 2 1 3 ) 其中，h = ( m ，1 1 k ) n l ，n ，k ；0 ，士1 ，土2 、3 i 一原胞的边长。 周期性边界条件已经被广泛用丁计算机模拟中，但它也有一些缺点。周期性原胞明显的限制是 它不可能得到波长大于原胞长度的波动。这在一定情况下引起些问题。图2 2 为周期性边界条件 的示意凹。 2 7 分子动力学模拟的系综 为了准确反映实际的物理过程，m d 模拟总是在一定的系综下进行。 2 7 1 系综的分类 微正则系综( n v e ) ，是孤立的、保守的、系统的统计系综。在该系综中，系统地原子数 r 、 体积v 和能量e 都保持不变。 w _ * t l 系综n v t ，在该系综中，系统的原子数、体积v 和温度r 都保持不变。且总动量为零。 等温等压系综n p t 在该系综中，系统的原子数、压力p 和温度r 都保持不变。 等压等焓系综n p h ，在该系综中，系统的原子数、压力p 和焓值h ( h = e + p v ) 都保持 不变 巨i e , l l 系综u v t ，该系综用于处理粒子数变化的系统。在该系综中，系统的化学势“、体积矿 和温度丁保持不变，但系统的粒子数n 与外界进行交换而不断变化。 2 7 2 系综的控制技术 系综的控制技术常见有调温技术和调压技术，在此介绍几种具体方法。 l 、调温技术 系统的温度t 与动能k 有关，其关系为 置：兰车：_ ( 3 n - n c ) k 日r ( 2 “) 一11 一。 式中，n 和 c 分别为系统的原子数和约束数；屯是b o l t z m 柚n 常数；e 是原子f 的速度。 ( i ) 速度标定法p i 】 速度标定法是保持系统温度恒定虽简单的方法，其具体做法是每隔一定的模拟步数，将原子运 动的速度乘以修正系数c ，c 表达式为 o 匡 t ( t ) ( 2 1 5 ) 式中t o 和t ( t ) 分别为系统的参考温度和在t 时刻的瞬时温度。 ( 2 ) b e r e n d s e n 热浴法【3 2 i b e r e n d s e n 热浴法假想系统与个温度为期望值得虚拟热浴相接触。该方法在每一步都对速度进 行标度，标度因子为 c =履t t 蕊l t ( t ) 8 ( 2 1 6 ) 第一章分了动力学模拟方法 式中瓦m 和r ( ，) 分别为系统的期望温度和在，时刻的瞬时温度；a t 为计算步长；t r 为耦合参数 通常取为0 4 p s 。 2 、调压技术 ( 1 ) b e r e n d s e n 压浴法1 3 2 】 b e r e n d s e n 压浴法假想系统与一压浴相接触。其做法为对模拟系统的体积乘以压力标度冈子o ， 对原于质心坐标乘以c 7 3 。 q ；1 一竽【b 。一p ( f ) 】。 l 一 ( 2 1 7 ) 式中，k _ 乖l l t e 为耦合参数；和p ( t ) 分别为系统的期望压力和在f 时刻的瞬时压力 ( 2 ) p a 币n e l l o - r a h m a n ( p r ) 方法 p - r 方法用于处理品格的形状和体积都发生变化的情况，可实现对原胞施加拉伸剪切以及混合加 载情况的模拟在材料科学中得到了广泛应用。具体实现见文献，此处省略。 2 8 宏观统计量的特性统计 在求出模拟系统的n 个粒子在相空间的轨迹之后，就可以应用统计物理原理计算该系统的宏观物 理特性和机构特性，从而完成分子动力学模拟的全过程。 由统计物理学可知，定义在相空间各点上的物理量的统计平均值都是系综平均( e n s e m b l ea v e r a g e ) ， 在通常条件下，它也是时间平均，即 4 = ( ) 。= ( 一) ，= ! l 叫i m l ： 一 ，协p 4 ( r ) p ( 2 1 8 ) 式中，( 一) 。为系综平均值；印) ，为时间平均值；r 为位置：p 为动量；t 为时间；f 为时间参变量。 由于不同的研究目的所需要计算的统计参数不同，这里参考a l l e n 的著作8 】列出几个与本文研究 有关的物理产量的统计形式： ( 1 ) 系统内能e e = ( ) = ( k ，+ ，+ “f 【o o ) ) ， ( 2 1 9 ) 式中，右端第一项表示动能，第二项为势能。 ( 2 ) 温度t r = 瓣= 击t ；w 仁2 。， 式中，k 8 为b o l t z i n a n n 常数。 2 9 计算中的一些技巧 2 9 1 无量纲化过程 9 东南人学硕十论文 在m d 模拟计算中，由于涉及到大晕的浮点运算和指数运算，茯i 此一般情况下，需要对计算公 式进行无量纲化处理，以便提高计算机运算效率。无龟纲化以原子质量川、势函数中的势阱常数8 和 平衡常数仃为基准进行无量化。 首先获得原子位移和能量的的无量纲值： 原子位移：r = r d 能量： e + = e 然后推导出其它各量的无量纲值： 帆r = 古后 温度：t = r ， 速率： 厩 ”i 然后推导出其它各量的无量纲值： 时间： 。t 压 h i 、i 温度： ，= 屯t i z 速率 儿v 店 作用力：厂= f a z 2 9 2 位能截断 ( 2 2 1 ) ( 2 2 2 ) ( 2 2 3 ) ( 2 2 4 ) ( 2 2 5 ) ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) ( 2 2 9 ) 对于分子数为n 的模拟体系，原则上任何两个分子间都存在相互作用，那么在计算体系位能时 须进行n ( n 1 ) 2 次运算，一般情况下要占总模拟时间的8 0 左右。非常消耗机时。为提高计算效率， 在实际模拟过程应进行势能截断，晟为常用的方法是球形截断法，截断半径一般取2 ，6 d ( d 为分子 的直径) ，这样对截断距离之外分子问的相互作用就可以忽略，模拟过程中减少了计算量9 ”。 2 9 3 邻居列表的处理 对于大型粒子系统或者复杂结构，m d 模拟是非常耗机时的，有时候出一个数据需要几天的时间。 所以当系统中粒子数量很大时，可以利用邻域列表法来判断粒子的分布情况，进一步节省计算机时。 v e r l h 算法的邻域列表法如图2 3 所示。对于系统中的每个粒子都可以建立如图的邻域表，图中- 为 截断半径，为邻域半径。如图示，计算l 粒子受力时，只计算裁断半径以内的2 、3