（机械设计及理论专业论文）纳米单晶铜杆的分子动力学模拟与性能研究.pdf

西北工业人学硕l 论史摘要 摘要 本文主要应用分子动力学方法对纳米单晶铜杆的表面效应及拉伸、弯曲性 能做了分析和研究。 对分子动力学方法的基本原理进行了总结，研究了分子动力学模拟方法及 程序编制，并以l j 势为例给出了分子动力学方法的具体实现。 运用分子动力学方法对纳米单晶铜在接近绝对零度下的表面效应进行了分 析，给出了表面原子位置变化、能量大小、应力的大小。 研究了纳米单晶铜杆在接近绝对零度下的拉伸下的应力应变曲线，并计算 了杨氏模量。拉伸作用下纳米杆的截面形状变化明显。对其拉伸变形的内部机 理进行了说明。研究了铜纳米杆在高低两种应变率的变拉伸以及纳米单晶铜杆 的弯曲，观察到了位错的产生。 研究表明，纳米材料表面效应不可忽略，未受外力作用下，材料表面内部 有应力存在。拉伸作用下的应力应变曲线变化趋势并不与宏观介质连续力学下 完全相同，不同的应变率会对材料变形造成不同的影响。弯曲下的纳米铜杆具 有与宏观尺度下不同的力学特性。 关键词：分子动力学嵌入原子法弛豫 纳米 两北t 业大学砸上论文a b s t r a c t a b s t r a e t s u r f a c ee f f e c t ，t e n s i o na n db e n d i n go fn a n os i n g l e - c r y s t a lc o p p e ra r er e s e a r c h e d o nt h eb a s i so fm o l e c u l a rd y n a m i c su s i n gt h ee a m ( e m b e d d e da t o mm e t h o d ) p o t e n t i a l t h eb a s i cm e t h o d so fm o l e c u l a rd y n a m i c sa r es u m m a r i z e d i m p l e m e n tw a y so f m o l e c u l a rd y n a m i c sa r es t u d i e d as i m p l em o l e c u l a rd y n a m i c sp r o g r a mi sg i v e n s o m ee x a m p l e sa r es h o w no nt h eb a s i s o fm o l e c u l a rd y n a m i c s u s i n g t h el - j p o t e n t i a l ar e l a x a t i o ns u r f a c ee f f e c ti ss i m u l a t e do nt h eb a s i so fam o l e c u l a rd y n a m i c s u s i n gt h ee a m ( e m b e d d e d a t o mm e t h o d ) p o t e n t i a lb ya p p l y i n gap e r i o d i cb o u n d a r y c o n d i t i o nf o rc o p p e r s u r f a c ea t o m sp o t e n t i a l ，s t r e s s a n dp o s i t i o n c h a n g e s a r e p r e s e n t e d t e n s i o ns t r e s sa n d s t r a i nr e l a t i o na l o n gt h ezd i r e c t i o ni sr e s e a r c h e da b o u t s i n g l ec o p p e rp o l en e a rt h ez e r ot e m p e r a t u r e t h ee l a s t i cm o d u l u s i s c o m p u t e d t h e s i z eo fc r o s sc h a n g ei sr e m a r k a b l eu n d e rt e n s i o n s o m ee x p l a i n sa r eg i v e na b o u t i n t e r n a ld e f o r m a t i o n su n d e rt e n s i o n ，t h et w os t r a i nr a t e sa r ea p p l i e d t h eb e n d i n g o f n a n o s i n g l ec o p p e rp o l ei s s t u d i e da l s oi t p r o v e st h a tl l a n os i n g l e c r y s t a lo fc o p p e r p o l eh a st h ed i f f e r e n c ed e f o r m a t i o nc o m p a r e d w i t hm a c r o s c a l e t h ed i s l o c a t i o n sa r e s e e nf r o mt h ea t o ms n a p s h o t i ti sf o u n dt h a ts u r f a c ee f f e c ti so b v i o u si nt h en a n o s c a l e n a n o m e t e rm a t e r i a l s t r e s si s n o tz e r ow i t h o u tf o r c e t h e s t r e s s s t r a i nc u r v e sa r ed i f f e r e n tf r o m m a c r o s c a l es i t u a t i o n v a r i e ds t r a i nr a t e sh a v ed i f f e r e n te f f e c t st on a n o m e t e rm a t e r i a l d e f o r m a t i o n c o m p a r e d w i t h m a c r o s c a l e ，n a n o s c a l e m e c h a n i c a lb e h a v i o u ro f b e n d i n g i sd i s s i m i l a r i t y k e y w o r d s ：m o l e c u l a rd y n a m i c se m b e d d e d 。a t o mm e t h o d r e l a x a t i o nn a n o 西北一 业大学顽j 一论文第一帝绪诠 第一章绪论 1 1 纳米科技与分子动力学 1 9 9 0 年被认为是纳米科技正式诞生的年代。纳米科技是现代科学( 量子力 学、介观物理、混沌物理等) 以及先进技术( 包括微电子、技术、电子计算机、 扫描隧道显微技术等) 相结合的产物。它是在纳米尺度( 0 1 n m m 0 0 n m ) 上研 究自然界现象中原子、分子行为和相互作用规律，旨在深化对客观世界的认识， 并在此基础上创造出性能独特的产品。纳米使得人类在认识和改造自然方面进 入到一个新的层次，能够进一步开发出物质的潜在能力，因此，它的发展无疑 地将深刻影响国民经济和现代科学技术的未来。纳米尺度的新现象和新工艺的 发现，新实验工具和模拟工具的兴起，为纳米技术的进一步发展提供了新的契 机。纳米技术向人类揭示出了一个可见的原子、分子世界，人们可按自己的意 愿操纵单个原子和分子，实现对微观世界的有效控制，生产出具有特定功能的 产品，纳米技术向我们展示了一幅幅诱人的蓝图。 纳米技术不是小尺寸技术的延伸，不能简单地将宏观连续介质力学套用到纳 米材料上，也不能简单地认为是尺度的减小。在纳米尺度下，材料表现了与宏 观材料截然不同的性质，表现出了许多优异的特性。纳米材料微粒尺寸相当或 小于光波波睦，周期性的边界条件被破坏，光、电、热力学性能呈现明显的小 尺j r 效应，如纳米微粒的熔点远低于块体，纳米单晶材料的力学性能显著提高。 表面与界面效应是纳米材料的另外一个显著特点，纳米材料的表面积相当大， 当粒径降至1 1 m a 时，表面原子数比例已达到9 9 ，增加了表面能，从而使纳米材 料具有很高的化学活性。纳米粒子还具有贯穿势垒的能力，人们称之为隧道效 应，扫描隧道显微镜s t m 正是基于量子的隧道效应，而于1 9 8 2 年由b i n n i g 和 r o h r e r 1 宣布研制成功。 对于纳米材料的研究形成了两种研究手段，一种是采用实验手段，如扫描探 针显微镜、原子力显微镜( s t m ) 的出现提供了在原子尺度上观察表面及其变化的 有效手段。另一种是采用数值模拟方法，它主要以计算机作为模拟手段，又可 分为分子动力学方法( m o l e c u l a rd y n a m i c sm e t h o d ) 和蒙特卡罗方法( m o n t ec a r l o m e t h o d l 。这两种方法都是基于原子尺度上的模拟方法，蒙特卡罗方法作为一种 两北丁业人学颂 论文 笫一章绪论 统计学的概率性方法，可以深入到相空间研究m a r k o v 链的随机行为，而分子动 力学作为一种确定性方法，它可以跟踪每个粒子的运动，给出材料的动力学量。 本论文中只考虑分子动力学方法，故仅对分子动力学方法加以介绍。 分子动力学计算机模拟是研究复杂的凝聚态系统的有力工具，这一技术既能 得到原子的运动轨迹，还能像做实验一样做各种观察。分子动力学模拟具有沟 通宏观特性与微观的作用，对于许多在理论分析和实验观察上都难以观察的现 象可以做出一定的解释，因此它在物理学、化学、生物学和材料学等许多领域 中得到广泛引用。在揭示各种现象的本质和探索新规律的微观研究中取得很大 的成功，被认为是2 0 世纪以来除理论分析和实验手段之外的第三种科学研究手 段，人们称之为“计算机实验”手段【2 j 。 基于粒子之间基本作用的分子动力学方法是近年来模拟纳米尺度力学行为 的重要数值手段，分子动力学计算方法的关键步骤是确定原子间的作用势，然 后通过求解经典牛顿力学体系的正则哈密顿方程组口】求解问题。采用不同的原 子间相互作用势，可衍生出不同的分子动力学分支。 经典分子动力学通过实验结果或经验模型确定原子间作用势，其优点是计 算量较小，适用于解决较大规模的问题，缺点是可移植性差，针对不同的问题， 需要确定不同的经验参数。 最早采用的原子间作用势是适用于原子间共价键作用的对势模型( p a i r p o t e n t i a l ，比较经典的对势模型有b u c k i n g h a m 势【4 1 、m o r s e 势口1 以及附加势】 等。对势模型仅考虑原子对之间的作用力，无法考虑化学键的方位交化。鉴于 此，人们开始考虑粒子间的多体作用，构造出多体势结构。多体势模型不仅通 过对势项考虑r 原子排斥作用、静电作用、极化作用等，还通过三势项考虑了 化学键阃的弯曲以及多体耦合等作用。该模型被成功地应用于碳纳米管，纳米 电子器件等的模拟。考虑到原子所处位置对电子密度的影响，人们先后提出了 确定原子间作用势的原子嵌入法1 7 j ( e m b e d d e da t o mm e t h o d ，e a m ) 和修正内嵌 原子方法【8 i ( m o d i f i e de m b e d d e da t o mm e t h o d ，m e a m ) 。e a m 和m e a m 方法特 别适于解决界面和表面问题。在很多情况下，原子局部环境( 比如微缺陷、表 面、界面等) 往往不同于整体环境，为了考虑这种差异的影响，人们提出了可 变电荷模型 9 l ，这种模型非常适用于解决边界问题。 为了克服传统分子动力学可移植性差这一缺陷，人们考虑直接从量子力学 轨道理论出发获取原子间作用势，基于量子力学的分子动力学称之为现代分子 动力学。现代分子动力学包括紧束分子动力学( t b m d ) 、密度泛函分子动力学 两北 _ 业人学硕f ：论文 第一章绪论 ( d f m d ) 和第一原理分子动力学( f p m d ) 。紧束模型通过考虑电子的波动性 而提高了其精度和普适性。t b m d 主要运用于硅、无定性碳的结构和动力学特 性等的研究。近几年t b m d 得到了广泛的发展和应用，是目前应用最为广泛的 分子动力学。d f m d 基于量子力学密度泛函理论，直接从基本原理考虑电子云 结构，模拟更为精确，可移植性更好，同时计算量更大，对计算机要求更高。 d f m d 仅能应用于简单的材料系统，而且计算规模仅有几十到数百个原子【 o , i i i ， 但模拟结果表明d f m d 的计算非常精确。计算精度和计算规模永远是一对矛盾， 为了协调这个矛盾，人们开始考虑综合t b m d 和d f m d 各自的优势，形成所 谓的f p m d 。f p m d 已经被广泛应用于固体微结构的演化及纳米力学行为分析 等方面【j “，报道的结果充分浼明了其模拟的精确性。 1 2 分子动力学模拟对机械设计理论研究的必要性和推动作用 随着纳米科技的到来，许多新的学科产生了，例如纳米电子学、纳米生物 学、纳米材料学、纳米机械学等。人们的注意力逐渐从宏观物体转向小尺度及 相应的器件，其中微机械系统( m i c r om a c h i n e ) 称微型机电系( 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 ，m e m s ) ：j f ：其取得了成功，并正被拓展应用于各种工业 过程。由图1 1 可知，分子动力学正是处于m 尺度下的研究方法。 长度( m ) 图1 - 1不同模拟方法所对应的空间和时间尺度 西北工业夫学顾1 论义第一章绪沧 由于微型机械在尺度、材料、构造、制造方法和工作原理上与传统机械截 然不同，因而在研究方法和研究内容上也与传统机械不同。研究在纳米尺度下 的运动变换、材料的力学性能以及在原子、分子尺度上接触的相互作用、摩擦 润滑等尤为重要。如在微机械中的润滑膜处于纳米范围，此时宏观摩擦学已不 适用，必须从原子、分子的相互作用来研究；还有微观接触上必须考虑到表面 与界面的相互作用，正如清华大学的余寿文提出：对于纳米尺度下研究。必须 重新来审视不同作用力在纳米下起的作用，一些以前忽略的力，如表面作用， 必须重新考虑。 利用分子动力学，人们己在纳米摩擦学中取得了可喜的成果。b u l l o u g h 【i 驯 以及c o t t e r i l l 和d o y a m a h 】分别采用b o r n m a y e r 和m o r s e 势对位错核心下的晶 格构形及位错的扩展进行了分析和模拟，d e c e l i s 【l5 j 等人采用j o h n s o n 和m o r s e 势分别对d f e 和铜中裂纹尖端构形在加载过程中的演化进行了模拟，z h o u 等 人采用多体势对f c c 中晶体中位错的运动及加氢后对位错运动的影响进行 了模拟，c a g i n 等利用分子动力学对纳米齿轮和转子进行了研究。 对于表面和界面问题。b a s k e s 工作皇且【”i l 剐运用分子动力学分别研究了n i ，s i ， a l s i 等界面的断裂、剥离等纳观力学行为。l i 等1 1 9 j 应用e a m 分子动力学模拟 了压应力或热载作用下铜晶体中的微裂纹愈合，清华大学的杨卫1 2 0 j 等人提出了 三重嵌套模型。关于固体表面粘着微观接触问题，l a n d m a n 和l u e d t k e 等人提 出2 1 i 用分子动力学模拟了硬的镍探针与软的金基片之间趋近与分离过程，随后 他们对于离子型晶体c a f 2 2 _ i n t 22 1 、半导体硅材料之1 b l t 2 3 1 之间的粘着及接触分离 的过程也得到同样的结果，清华大学的冷永胜1 2 4 】运用分子动力学做了二维网4 性 平冲头与c u 的接触，给出了接触表面在零作用下的压力分布。 实践证明，分子动力学模拟的特点是能够有效地进行微观的动态过程分析， 计算出其空间尺度可咀达到纳米量级，时间尺度甚至可达1 0 “5 s ，因此特别适合 于纳米摩擦学研究【2 l 。迄今为止，纳米摩擦学涉及的各个研究领域几乎都有分 子动力学模拟计算的研究报道。在固体表面接触与分离、表面牯着行为与材料 转移、微观摩擦与粘滑现象、微磨损与微切削加工，以及固体与液体界面作用、 纳米薄膜流变特性、有序分子膜的构性关系等方面，分子动力学模拟研究都取 得了重大进展，其中不少研究结果已为实验所验证。发展至今，分子动力模拟 技术不仅成为纳米摩擦学理论分析的基本手段，而且有望用于进行纳米结构的 设计和性能评定。 西北d 世大学顺 ：论文 第一章绪论 1 3 分子动力学方法的发展 分子动力学作为一种计算机模拟手段，己在材料设计、纳米摩擦学、固体 力学、高分子等领域得到应用。分子动力学模拟已从传统的经典分子动力学发 展为多体分子动力学( m a n y b o d y m o l e c u l a rd y n a m i c s ，m b m d ) ，可变电荷分子动力 学( v a r i a b l e - c h a r g em o l e c u l a rd y n a m i c s ，v c m d ) ，紧束分子动力学( ( t i g h tm o l e c u l a r d y n a m i c s ，t b m d ) ，密度泛函分子动力学( d e n s i t y - f u n c t i o n a ld y n a m i c s ，d f m d ) 等。 为了克服传统分子动力学可移植性差这一缺陷，人们直接从量子力学出发来获 取原子间的相互作用，如密度泛函分子动力学( d f m d l ，第一性原理分子动力学 等。 由于分子动力学方法本身的限制，再加上计算机硬件的发展赶不上人们对 运算速度的要求，要求解复杂的问题会遇到很多困难，要想像有限元数值模拟 方法一样得到广泛应用，成为一种得到工程实际应用的并为普通人所能掌握的 一种数值运算手段还有一段距离，还要很多工作要做。 寻找扩大计算规模的运行算法是以后重要的研究方向，人们已在降低运算 规模的线性规划算法( 1 i n e a rs c a l i n g a l g o r i t h m s ，o ( i v ) ) 的研究方面取得了重要进展 1 2 6 - 2 9 。还有采用并行计算来提高求解问题的规模，即通过并行处理器解决计算 效率的问题。随着计算机技术的发展，廉价的高性能微处理器和高速网络技术 以及开放的系统并行软件的出现，采用微机集群进行计算，可显著提高运算速 度，速度最高的微机集群已达每秒l l 万亿次浮点运算口。对算法效率进行优化 的只有一些探索性工作，如冷冻原子法，但该方法目前用于材料的变形分析， 还不能用于热力学分析。 由于分子动力学本身为所模拟的时间尺度、空间尺度所限，为了扩大求解 问题的规模，人们提出了多尺度的概念，提出构建基于各种尺度耦合的方法， 即将适用于各种尺度的方法组合在一起，这是个急需解决的课题，如基于有限 元与分子动力学相结合的方法。杨卫建立了宏观细观微观三重嵌套模型用于描 述裂尖断裂。a b r a n a m j 结合从头算法( a bi n i t i oc a l c u l a t i o n ) ，提出宏观原子从 头分子动力学方法( m a c r o a t o m i s t i c a b i n i t i od y n a m i c s ，m a a d ) 。西北工业大学的 崔俊芝提出用窗函数的概念来构建多尺度方法。 有效地对分子动力学的数据进行后处理也是一个重要的研究方向。对数掘 进行后处理，获取有用的信息也是一个很烦琐的工作。现有的一些分子动力学 西北：c 业大学硕l 论文 第一章绪论 软件功能不够强大，只是偏重于某一行业， 软件，可维护性不强。还没有出现集建模、 软件。 通用性不高，还有一些软件为自山 求解、后处理于一一体的分子动力学 不能实现所见即所得。将虚拟现实技术( v i n u a lr e a l i t y ，简写为v r ) 引入分子动 力学当中，可以对模拟结果实现进行动态显示及结构预测。应用v r 技术，将有 助于深入理解固体中界面的结合强度运动、流体中原子的运动。 l 。4 本文的研究内容 研究分子动力学方法，首先应弄清它的基本原理、实现过程、应用范围等。 目前，分子动力学方法在许多领域都有应用。对于在纳米尺度下材料的力学性 能有必要从理论上给出一定的解释。本文主要研究在纳米尺度下单晶铜在接近 绝对零度下的表面弛豫、表面原子重新排布、周期性边界条件的影响，分析了 在受力下的应力变化，模拟了受拉与弯曲情况下的变形情况。 本论文主要将分子动力学方法应用到微机械中材料的力学性能分析上，并 给予一定的合理解释，从而加深对纳米材料力学性能的了解。 第二章分了动力学模拟技术 2 1 基本原理 第二章分子动力学模拟技术 随着纳米科学技术的发展，应用连续介质力学来处理处于原予、分子状态 的固体或液体己不太合适。在经过长时间的探索之后，人们发现从统计物理学 中衍生出来的分子动力学模拟，实践证明是一种描述微观现象的有效方法f 3 3 1 。 这种方法已得到人们越来越广泛的应用，并应用于各种领域。分子动力学根据 粒子间相互作用势，计算多体系统的结构和动力学过程，并可以计算物质的结 构和性质。分子动力学的基本原理是建立一个粒子系统来模拟所研究的系统， 系统中各粒子之间的相互作用根据量子力学来确定，人们通过求解粒子的动力 学方程组来得到各粒子在相空间的运动规律和轨迹，然后按照统计物理原理得 出该系统相应的宏观物理特性。 由于分子动力学模拟具有沟通宏观特性与微观结构的作用，特别对于许多 在理论分析和实验观察上都难以了解的现象可以作出一定的微观解释，并可以 模拟一些极端条件下的微观现象，并得到广泛应用。被公认为是理论和实验观 察相联系的第三种科学手段。 分子动力学计算中基于以下假设口“。 ( 1 ) 所有粒子之间的的运动都遵守经典牛顿力学定律； ( 2 ) 粒子之间的相互作用满足叠加原理。 对于粒子数为1 1 的物理系统，其力学描述的哈密顿形可写为p 5 j 口= 票 f ：1 “2 一，” ( 2 1 ) 9 2 面 净1 “一” ( 2 1 ) 哈密顿函数为 0 h p j 。一- = 一 o q i h ：y 盟+ u 争2 ( 2 2 ) ( 2 3 ) 两北t 业大学硕卜论文 第二章分了动力学模拟技术 其中，g 。，p 。分别为第i 个粒子的广义坐标和动量；m ，为该粒子的质量；u 为系绕 的总势能。 该系统力学的牛顿方程研式为 ，= f ( ) m 。 i = l ，2 ，n ( 24 ) f = 一v 。u ( r 。，t ，)( 25 ) 式中，f 是第i 个粒子所受到的合力，是该原子的位置坐标。 粒子之间的相互作用是通过势函数来计算的，而势函数的选择非常重要， 它决定了计算量和所计算模型与真实系统的近似程度。人们因此提出了若干势 函数形式，并且势函数已从最初的两体势函数发展到多体势函数，人们已经提 出了很多势函数形式，并在实际中使用。碱金属盐、碱土金属盐、硅氧化合物、 蛋白质、半导体、过度金属等多种物质的相互作用势已建立了，并取得了应用， 随着研究工作的不断深入，尤其是精度高的相互作用势函数的确立，可以极大 地提高运算精度。可以根据所研究的不同对象，选用不同的势函数。 由所选定势函数，即可求出系统的总能量，即 u = u ，( ，) ( 2 6 ) 这样，粒子之间的相互作用力可由下式得出 f , = - f j ，= 掣 眨， r q 2 r i r = 巴 ( 2 8 ) ( 2 9 ) 根据粒子之间的相互作用，对运动方程进行积分，可以得到粒子的运动细 节，并由统计物理学原理得出系统的宏观量，并可对此进行分析，就好比在做 实验一样，可碍到系统的物理特性。 由统计物理学可知，定义在相空间各点上物理量的统计平均值都是系统 平均( e n s e m b l ea v e r a g e ) ，在通常条件下，它也是时间平均，即 1o 矿= ( 妒) 。= ) 。= ；i m 【 r ”( f ) ，p ”( f ) d f ( 2 10 ) f “ ( 2 1 0 ) 式中，( 妒 一。为系统平均值， ，+ ( “。啄1 ) 】) ， ( 2 - 1 1 ) 式中第一项表示动能，第二项表示势能。 ( 2 1 温度t 卜壶( 莩哪 但1 2 。为b o l t z m a t m 常数。对于含有n 。个内部约束的n 个粒子，则系统 温度的统计形式为 r 。面孤1 ( 莩) ( 2 j 3 ) ( 3 ) 压力p p v = 础月7 + ( f ) ， ( 2 - 1 4 ) 其中，f 是作用在i 粒子上的台力，v 为体积 ( 4 ) 径向分布函数g ( r ) 径向分布函数表示任一给定原子在距离为r 处的平均数密度p ( ，) 与在同 样的总密度情况下的理想气体中某原子在距离为r 处的密度之比。通过测量 它，我们可以知道结构的无序化程度，并可以与中子衍射作对比。在理想 气体中，g ( r ) = 1 ，对理想晶体的原子密度为p 。，则径向分布函数为 g ( ，) = 盟p , v ；孝等 ( 2 1 5 ) ( 5 ) 静态结构因子【3 7 】 静态结构因子是表征结构无序化程度的物理量，其表达式为 9 两北工业大学硕i 论文第二章分子动力学模拟挫术 蹦h 帅隆p o x b 1 6 ) s ( k ) = ( 1 a ，) e xo ) j ( 2 1j = jl 是总原子数，k 为倒格矢，f 原子的位置矢量，对理想晶体其结构因 子为l ，对理想流体，则为0 2 2 分子动力学模拟的实现及程序编制 2 2 1 初始化条件的设定 要对一个系统进行模拟，应对体系中所有粒子赋予初始位置和速度，所选粒 子位置应与所模拟的对象相适应。一般可以采用随机分布的方式，但很可能会 造成粒子之间的重叠。通常，将粒子初始位置放在一个品格上，常见晶格有简 单立方、体心立方、面心立方等，分别如图2 1 、2 - 2 、2 - 3 所示。在对于大 量流体分子系统的模拟过程中，品格结构会消失，成为液体结构。 图2 - 1 简单立方图2 - 2 体心立方 图2 - 3 面心立方 初始速度的选择通常为一随机分布，但也需要满足一定的温度条件或者指 定的动量为零的要求，即 p = ，q = 0 ( 21 7 ) 对于平衡分子动力学系统，动量为零可以作为系统模拟是否正确的一个判 断条件，下面以液态氩为例加以说明，选用l - j 势【3 4 。 叫巾。州2 删 仁旧 西北工业夫学顺卜论义 第二章分子动= ：j j 学模拟技术 上式中，和盯分别是能量参数和长度参数，对于液态氩来说，常取心； 1 1 98 k ，口= o 3 4 n m 。设对比温度t = 1 ，对比密度p 。= o 8 。这里的单位均采用 l j 势的对比单位，取长度单位盯、质量单位肌、能量单位为基本单位，由这 些基本荦位，可得到其它单位。例如势能的对比单位可取为中：，= 中- ，弦。初始 构形为2 0 * 2 0 的平面二维结构，初始构形图如图2 - 4 所示。时间步长为0 0 0 5 ， 每1 0 0 步作一次平均，总步长数为1 0 0 0 0 步，下面为模拟过程的系统动量随时 问的变化图。图2 5 的横坐标为系统模拟的步长数，纵坐标为系统的动量( 采用 对比单位) 。 1 00 0 j y | 口 i 。 oo o ! l - 1 00 00 2 00 0 _ 一 - - ， 一一一 2 00 0 1 00 0 00 01 00 02 0 0 0 x | o 图2 4 初始构形图 10 0 05 0 一 p1 0 ，0 0 一 - 05 【) l 一1 0 0 一一 一一 rr i 1 02 0 0 0 4 0 0 0 8 0 0 08 0 0 010 0 0 0 n u m b e r 图2 - 5 动量随时间变化图 西= | h - l l 业大学硕二l 论史 第二市分了动力学模拟技术 在热平衡状态下，有下式成立，系统中的粒子速度分布满足 m a x w e l l b o l t z m a n n 分布，对于原子体系 p ( v 。) = ( 川，2 础。t ) 1 2e x p ( 一三m v ：。r )( 2 1 9 ) 这里，p ( v 。) 为速度分量的概率密度函数，同样，对y 、z 方向也一样。 图2 - 6 速度概率分布 图2 6 为系统在对比湿度丁+ = l ，根据( 2 。1 9 ) 式所产生的在x 方向的平衡 速度概率分布，x 轴代表速度，y 轴代表速度的概率分布。 因此，从系统的速度分布可以检验系统是否处于平衡态，这是判断系统是 否处于平衡状态的重要判别方法。对系统中粒子初始速度全赋零也可以，我们 只是给予初始条件，在系统的模拟中，还要对温度进行控制。 2 2 2 力的计算 根据给定的初始条件和边界条件，就可以求解运动方程组，来求得粒子在 相空间的运动轨迹，进而由统计物理学得出该系统的宏观特性。 而在分子动力学模拟中对力的计算最费时的部份是计算每一粒子上的力， 并考虑所有邻近原子作用于i 粒子上的力。如果仅考虑计算与另一粒子的最近 邻的作用，则我们必须计算n f n 一1 ) 2 个成对作用【3 引。如果不使用任何技巧的话， 力的计算时间与n 2 成正比。有几种有效的方法可以用来提高计算力的效率。 最常用的是截断半径法，粒子删的作用力随着它们之间的距离增大而减小， 两北工业大学倾【。论义 第二草分了动力学模拟技术 在所选的截断半径r c 以外，粒子列作n j i j , j 零，所以只计算r c 以内粒子f 自j 的相 互作用力。当粒子数目很大时，可采用v e r l e t 列表法或近邻列表法，在该方 法中，第二截断半径r l r 。被引入，且在我们计算相互作用之前，制作以粒子l 为中一1 5 , 的v e r l e t 链表，如图2 7 所示。 7 图2 - 7v e r l e t 表 计算粒子l 的受力时，只需计算截断半径以内的2 、3 、4 三个粒子对l 粒 子的作用力之和，其它粒子对l 粒子的作用力忽略不计。计算过程中每隔一定 步数要更换此表。一般更新邻域表的步数为1 0 步至2 0 步。在此期间，较小的 r 1 r 。能够保证在模拟过程中位于邻域表以外的原子进入邻域表内。当5 、6 粒 子与1 粒子4 距离小于r 。或者7 粒子与1 粒子距离小于r i 时，则需要更新邻域 表。即在每涉中控制速度最大的粒子满足下述条件时f 3 9 】： ( m 脚v ，i ) 杀 ( 2 2 0 ) s t e p s 则需要更新邻域表。这个公式相当于检查原子的位移。 运行时间与n 正比的算法是元胞列表法，元胞列表法的思想如图2 - 8 所示。 模拟盒子被分成若干大小相等或大于r 。的元胞。给定元胞中某一个粒子只与相 同元胞或相邻的原胞相互作用。元胞列表法所用计算时间与n 成正比。对于二 维而言，粒子i 与9 个元胞有相互作用，对于三维而言，则与2 7 个元胞有相互 作用( 这里也包括i 粒子自身所在的元胞) 。模拟元胞被分成大小为r c r 。的元胞。 两北_ _ 【二业大学硕i j 论文 第一二章分子动力学模拟技术 ：i 图2 8 元胞列表法 对这些算法的效率进行比较是有意义的，有助于我们在模拟中找到合适的 算法，并提高计算效率。在v e r l e t 列表中，在三维情形下需计算距离的粒子的 数目为 4 。= - - = ：w r ； ( 2 2 1 ) 对于元胞列表法相应的数目为 = 2 7 p r ( 2 2 2 ) 将我们通常所用的参数值l j 势代入，则v e r l e t 计算的对距离数比元胞列表 法少，v e r l e t 方法在计算相互作用中更有效的事实促使a u e r b a c h 等口8 】使用两个 列表的联合方法：用元胞列表法来制作v e r i e r 列表。关于力的不同求解方法计 算的比较见a l l e n r 的著作。西安交大的吴江涛【4 0 j 则对短程力的上述各种方法作 了比较，并给出了计算效率与所给粒子数、密度之间的关系。 2 2 3 积分求解方法 在分子动力学中常用积分算法有两类：g e a r 预测校正算法和v e r l e t 算法。 其基本思想是给定时刻t 的位置、速度，通过求解运动微分方程组，给出下一 时刻f + 的位置和速度。 已知t 时刻的位置，( f ) 、加速度口( ，) 和前一时刻的位置，p 一，) ，由泰勒展 开式可得下一步的位最。 ，f + a t ) = r ( t ) + a t v ( t ) + 0 5 a t2 口o )( 2 2 3 ) r r f a t ) = r ( t ) 一z x t v ( t ) + o 5 a t 2 口( f )( 2 2 4 ) 对两式相加，则得出下一步的位置如下： 两北t 业大学硕| ：论义第二章分予动力学模拟技术 r ( t + a t ) = 2 r ( t ) 一r ( t a t ) + a t ：口o )( 2 2 5 ) 从此式就可以得到粒子在相空间的轨道，并臣该式的截断误差是o ( a t 4 ) 。 速度并不能直接得出，但速度可以计算系统的动能，速度由下式得出： v ( f ) = 坐学 ( 2 2 6 ) 在这罩，下一步的位置没有得到之前，难以得到速度项。另外，它不是一 个自启动算法。在t = 0 时，只有一组位置，所以必须通过其它方法得到卜出的 位爨。获得f 一出的位黄的方法这一是采用近似式 4 1 1 r ( - a t ) = r ( o ) 一a t v ( o )( 2 2 7 ) 速度v e r l e t 算法可以同时给出位置、速度与加速度，它的公式如下： ，o + a r ) ：，( ，) + p ( ，) 十二；盟，2 ( 2 2 8 ) z m v ( ，础) ：v ( 卅丛兰生盟，( 2 2 9 ) 这种算法要求的空间比较大，需要总共9 n 个变量存储空间 h o c k n e y 【4 2 1 提出的l e a p f r o g 算法是y e r l e t 算法的另一种形式。与y e r l e t 算 相比它有两个优点( 1 ) 包括显速度项，( 2 ) 计算量较小，但它的位置与速度不是同 步的，不可能在位置一一定时，同时计算动能对总能量的贡献。它的计算公式如 下： r ( f + a t ) = r ( t ) + a t v ( t + 亡f ) ( 2 3 0 ) v ( h i la t ) = v ( 卜i 1a t ) + z x t a ( ，) ( 2 3 1 ) 所需存储的变量有当前的位置r ( t ) 、加速度“t ) 、半步的速度v ( t o 5 a t ) 先求出半步的速度v o + 0 , 5 a t ) ，则当前的速度可n t 式求出 v ( t ) = o 5 ( v ( t + 0 5 a t ) + v o o 5 a t ) ) ( 2 3 2 ) g e a r 预测一校正算法4 3 l 首先根据泰勒展歼式，预测新的位置、速度与加速 度，然后根据新的位置计算得到的力计算加速度口( f + f ) ，再将此加速度与 t a y l o r 展开的所预测的加速度进行比较，两者之差用来在校正步里校正位置与 两北1 二业人学坝l 论文 第一章分子动力学模拟控术 速度。 对于四阶预测校正算法，我们定义位黄对时间的系列导数如下，0 、 = f ( a , o a t ) 、= o 5 ，2 ( d 2 础2 ) 、_ = 1 6 ，3 ( d 3 , o a , 3 ) ，则相应的泰勒展 开式预测项为： 瞄0 + 出) c ( t + a t ) 哆( ，+ a t ) ( f + a t ) 1l 01 o0 o0 鹾 这个矩阵为帕斯卡上三角阵，相应的校正步矩阵为 ( f + r ) ( f + a t ) ( r + f ) 巧( f + a t ) 喈o + f ) r ( f + a t ) 哆o 十血) 瞄o + ，) + c 0 c 1 c 2 q f 2 3 3 ) ( 2 3 4 ) 校正系数的值依赖于所要求解的差分方程，对于一阶运动方程 ，= 厂一)( 2 3 5 ) ，= 叶一r ，r 为预测的值，寸为校正的值。表2 1 为一阶运动方程的校 正系数。 表2 - 1 一阶运动方程的预测系数 v a l u e c 0 c lc 2 c 3c 4c 5 35 1 2 11 2 43 813 41 6 52 5 l 7 2 0l1 1 1 21 31 2 4 69 5 2 8 8l2 5 2 43 5 7 25 4 81 1 2 0 对于二阶运动方程 ，= 厂( ，)( 2 3 6 ) 西北工业夫学碗l 论文第二章分子动力学模拟挫求 a r = 嘭一哆，哆为预测值，是将瞄代入运动方程中得到的。二阶运 动方程的系数由表2 2 给出。 表2 - 2 二阶运动方程的校正系数 v a l u e c o c 1 c 26c 5 3ol 1 41 65 611 3 51 9 1 2 03 4l1 ，21 1 2 63 2 02 5 1 3 6 011 1 1 8l 6 0 对于如下的二阶运动方程 i = f ( r , k 、 ( 2 3 7 ) 对于c 。的值要作一调整，对于值5 要调整为1 9 9 0 ，对于值6 要调整为3 1 6 。对 于四阶预测校正法而言，需要1 5 n 的存储空间，对于5 阶预测校正算法则要 1 8 n 的存储空间。一般说来，展开的阶数越高，精度越好。但阶数取得过高并 不一定能提高运算精度1 3 “。 当然一个好的积分方法是我们所想要的。那么什么是一个好的算法，主要 从以下几方面去考虑，速度、较大的时间步长、能量是否守恒。最重要的是一 个算法能同时在短时间和长时间内预测所有粒子轨迹【3 引，事实上这样的算法并 不存在。牛顿运动方程是时间可逆的，我们的算法也该如此。可事实上许多算 法并不是时间可逆的。这就是说在某一瞬间将所有粒子的动量反向，体系并不 在相空间上反演。只有在有限短的时间步长情况下，这样的算法才是可逆的。 即使我们能找到积分可逆的方法，在数值上也不是真正可逆的。这是因为我们 所用的计算机也存在截断误差。 2 2 4 系统控制方法 1 调温方法 在n v t 系统中或n p t 系统中，常需要将系统温度调到期望值。常见的有 速度标定法【4 4 】、a n d e r s o n t 4 引、g a u u s s i a n 热浴1 4 6 1 、n o s e h o o v e r m 4 8 】热浴。 速度标定直接对速度进行标度。如果t 时刻的温度为t ( t ) ，速度乘以因子五 酣北工业大学顿七论文 笫一章分予动力学模拟技术 则温度的变化为 a t = ( 五2 一1 ) t ( t ) = ( 2 3 8 ) ( 2 3 9 ) 这里，一。是期望的参考温度。 a n d e r s o n 热浴中，体系与一热浴相耦合。与热浴的耦合由偶尔作用于随机选 择粒子上的随机力表示。与热浴耦合的强度由随机碰撞的频率来决定。 g a u u s s i a n 热浴法是在运动方程中加上一项“摩擦力”，其运动方程由下式给 出。 m 出, d ：2 r , = ，f 一笋v 。 ( 2 4 0 ) 笋= ( v ，一) ( 仇咋- ) 。1 ( 2 4 1 ) n o s e h o o v e r g 与g a u u s s i a n 方法采用相同的运动方程，但摩擦因子不同 警一c p 似一芳，q b a ：， 2 调压方法 在等压模拟中，可以通过改变模拟原胞的三个方向的尺度来实现体积的变 化，计有约束法、扩展系统法、c a r - p a r r i n e l l o 方法。 b e r e n d s e n 方法 4 9 是假想系统与一压浴相耦合，模拟原胞体积的标度因子为 兄，而原子坐标的标度因子为a n l 3 。则标度因子为： 五：i 一声兰( 尸一p 6 。m )( 2 4 3 ) p 这里，r ，为耦合参数，p 。为热浴的压力，p 为真实的压力。 扩展系统法中，个额外的自由度加入，相当于一个活塞作用在一个真实的系 统上，这个自由度的动能为去q 矿2 ，q 为活塞的质量，v 为系统的体积，通过 活塞的体积来调节内外压力平衡来保持系统平衡。系统的运动方程为 吝= f j ( m v 们) 一( 2 3 ) s 1 2 v( 2 4 4 ) 互 两北丁业人学颂 论文 第二章分了动力学模拟挫术 v = ( p p ) q( 2 4 5 ) p a r r i n e l l o r a h m a n 方法考虑晶体单胞的形状及体积大小都发生改变的情况 f ，在外加均匀静压，时，对于允许单胞形状及大小发生改变和起伏的等压方 法。考虑如图2 - 9 所示的单胞，其形状为三个基矢a ，b ，c 所围成的平行六面 图2 - 9 单胞 体，粒子的坐标可以表示为 = h t a + h 2 b + h