（凝聚态物理专业论文）固体表面动力学过程的第一原理理论研究.pdf

摘要 本论文首次向国内计算物理学界详细地介绍了目前国际上一种最先 进的计算方法，从头计算分子动力学方法( a bi n i t i o m o l e c u l a rd y n a m i c s m e t h o d ) 【赵字军等，物理学进展，第1 8 卷，第1 期，4 7 7 5 ( 1 9 9 8 ) 。 同时对硅表面、i i i - v 簇半导体表面、金属表面等几种典型体系的表面动 力学过程，用从头计算分子动力学方法和离散变分x a 【( d v - x a ) 方法进行 了第一原理的模拟计算。 从头计算分子动力学方法有机地结合了第一原理计算和分子动力学 技术，在真实的离子动力学中引入了表征在不同能级间跃迁速率的电子虚 拟动力学，使基于局域密度泛函理论的第一原理计算直接应用于统计力学 模拟成为可能，是第一原理研究动力学过程的理想工具。本文详细介绍了 这一方法的原理，范数不变( n o r m - c o n s e r v i n g ) 赝势和超软( u l t r a s o f l ) 赝势的 构造，以及在平面波基和上述两种赝势基础上的具体实现的过程和程序框 吼 所谓表面复位( s u r f a c er e s t o r a t i o n ) 是指表面重构得到消除，结构回 复到理想晶格表面的过程。我们用从头计算分子动力学通过完全的第一原 理动力学对s 和s e 在s i ( 1 0 0 ) 面的吸附以及吸附引起的表面复位作了较为 系统的研究，发现s 和s e 在s i ( 1 0 0 ) 表面的吸附结构对覆盖度很敏感，从 而导致了实验结果的复杂性，给出了以前实验结果之间不一致的原因。对 于s e 的吸附结构，在0 2 5 m l 和0 5 0 m l 的时候，为( 2 x 2 ) 的吸附结构，在 1 o m l 时确证能很好地消除重构，形成( 1 x 1 ) 结构。s 在s i ( 1 0 0 ) 面吸附时 结构就显得更为复杂，在o 2 5 m l 时，s 在桥位的吸附不稳定，在体系的 弛豫中自发地转移到穴位附近，形成2 x 2 结构。在o 5 0 m l 时，s 也自发 地从桥位迁移到弯曲d i m e r 结构的顶位附近，同样形成2 x 2 结构。这一结 果表明弯曲d i m e r 结构由于高低原子间的电荷转移，使其硅原子具有一定 浙江大学博士学位敝 的离子性，形成了0 5 0 m l 下s 在顶位附近吸附的结构，改变了认为二度 位才是最佳吸附位的观点。在覆盖度增至1 0 m l 时，验证了s 能破坏硅 d i m e r 结构，使s i ( 1 0 0 ) 2 x l 复位至( 1 x 1 ) 结构。 同时我们抓住吸附能在g a n 薄膜生长中对m g 掺杂的决定作用，用 从头计算分子动力学和离散变分方法计算了m g 和g a 在g a n 薄膜( o o o l l 面的两种不同极性下的吸附，间接地给出了m g 的掺杂与g a n 薄膜的极 性的关系。计算发现+ c 结构下g a 的吸附能比m g 要大约2 0 9 e v ，在g a n 薄膜生长条件下估算。m g 的吸附几率仅为g a 的1 0 气1 0 一。对于c 结构 而言，m g 与g a 的吸附能差别不大，仅为0 5 0 e v ，m g 的吸附几率为g a 的1 0 一1 0 。通过价电荷分布的分析知，由于c 极性中，表面的n 只有 一个悬挂键，而+ c 极性下有三个悬挂键，因此g a 的吸附能在+ c 下比较 大，比c 下大1 5 l e v 。对于m g 来说，由于其离子性很强，而且只有两 个价电子，对衬底的n 的悬挂键是一个还是三个影响不大，吸附能变化 很小，因此导致了m g 在c 结构下掺杂相对容易。 另外，我们用离散变分方法通过坐标优化的方法研究了乙烯在n i ( 1 1 1 ) 表面脱氢解离的可能的反应途径。给出了乙炔在各个吸附位的吸附结构， 并在假定乙烯在脱氢前没有迁移的基础上，对各种可能的反应途径作了模 拟研究，发现乙炔在平行桥位的吸附结构不稳定，通过1 2 4 a 的平移和3 0 度的旋转的耦合，迁移到垂直桥位，迁移过程中只有0 0 2 e v 的势垒。同 时给出了迁移途径中乙炔分子的碳碳键长，碳氢键角等结构参数在迁移过 程中的变化：。工一 a b s t r a c t a bi n i t i om o l e c u l a rd y n a m i c s ( m d ) m e t h o d ，o n eo f 也em o s ti m p o r t a n ta n d a d v a n c e ds c h e m e si nc o m p u t e rs i m u l a t i o n i sf i r s t l yi n t r o d u c e dt ot h es c i e n t i s t s i nc h i n g tu n d e rw o r k so f t h ea u t h o rf y u j u n z h a o e t c p r o g r e s s 加p h y s i c s v 0 1 1 8 ，n o 1 ，1 9 9 8 ( 4 7 7 5 ) 1 f u r t h e r m o r e ，s o m ed y n a m i cp r o c e s s e so ns u r f a c eo f t y p i c a ls y s t e m s 。s u c h 淞s i l i c o n ，i i i v n i t r i d es e m i c o n d u c t o ra n dt r a n s i t i o n m e t a l a r es t u d i e db y ：a bi n i t i om d m e t h o d c o u p l e d w i 血d v x qm e t h o d a bi n i t i om dm e t h o dc o m b i n e sm o l e c u l a rd y n a m i c sa n dd e n s i t yf u n c t i o n a l t h e o r y ，a n dp e r m i t st h ea p p l i c a t i o no fd e n s i t y - f u n c t i o n a lt h e o r y t om u c h l a r g e r s y s t e m st h a np r e v i o u sf e a s i b l e 。i ti sa ni d e a lt o o lt os i m u l a t ed y n a m i cp r o c e s s b yf i r s tp r i n c i p l e t h ef o u n d a t i o n s a n dp r a c t i c a li m p l e m e n t a t i o no f t h em e t h o d i sr e v i e w e d a sw e l la si t sd e v e l o p m e n ti nr e c e n ty e a r s t h ei d e ao fs u r f a c er e s t o r a t i o nb ya d s o r b a t e sw a sn a s c e n ti ns t u d i e so fm o r e t h a na q u a r t e ro f ac e n t u r ya g o h e r e ，t h ea d s o r p t i o na n ds u r f a c er e s t o r a t i o no f s i ( 10 0 ) ：sa n ds i ( 1 0 0 ) ：s ea r es t u d i e db y a bi n i t i om dm e t h o d i ti sf o u n dt h a t t h ea d s o r p t i o ns t r u c t u r eo ft h e s es y s t e m sa r ev e r ys e n s i t i v et ot h ec o v e r a g e ， w h i c hm a yc o n t r i b u t et ot h ec o m p l e x i t yi nt h ep r e v i o u se x p e r i m e n t a lr e s u i t s i nt h ec 黯eo f s i ( 1 0 ：s e ，i tt a m st ob ea ( 2 x 2 ) r e c o n s t r u c t e ds t r u c t u r ea tt 1 1 e c o v e r a g eo f0 2 5 m l _ o 5 0 m l w h i l ei t r e a c h e s 也e1 x lb u l ks t r u c t u r ea tt h e c o v e r a g eo f1 om l a st os i ( 1 0 0 ) ：s ，i ta p p e a r sa l i t t l em o r ec o m p l e x i t y a t 0 2 5 m lsa t o mt r a n s l a t e st ot h ec a v es i t ef r o mb r i d g es i t es p o n t a n e o u s l y ，a n d f o r m sa2 x 2s t r u c t u r e a t0 5 0 m l sa t o m sp r e f e rt oa d s o r b e dn e a rt h et o p s i t e so fb u c k l e ds i l i c o nd i m e r s w h i c hp r e s e n t sn e wa d s o r p t i o ns i t e st ot h e c o n v e n t i o n a lo p i n i o n a st h ec o v e r a g ei n c r e a s e st o1 o m l i ti sc o n f i r m e dt h a t sc a na l s or e s t o r et h es i ( 1 0 0 ) 2 x lt oi x lb u l ks u r f a c e i na d d i t i o n i ti sn o t i c e dt h a tt h ea d s o r p t i o ne n e r g i e so fm ga n dg ap l a ya n i m p o a a n tr o l ei nt h ed y n a m i cp r o c e s so fm gi n c o r p o r a t i o no n t h eg a nf i l m t h ea d s o r p t i o ne n e r g i e so fm ga n do ao nd i f f e r e n tp o l a r i t i e so fo a n ( 0 0 0 1 ) s u r f a c ea r ec a l c u l a t e db ya bi n i t i om dm e t h o da n dc h e c k e db yd v x am e t h o d ， 浙江大学博士学位论文 a n dt h ed e p e n d e n c e 。f _ i 曲o r 两妇c i i = i | ；j = 誓ot h ep o l a r i t yi sd i s c u s s e d i ti sf o u n d t h a tm gh a sl i t t l ec h a n c et ob ea d s o r b e do i l + cs t r u c t u r eb e c a u s eo ft h el a r g e d i f f e r e n c eo fa d s o r p t i o ne n e r g yb e t w e e ng aa n dm g i ti ss u g g e s t e dt h a tt h e m u l t i - d a n g l i n g - b o n do n - cs t r u c t u r er e d u c e st h ea d s o r p t i o ne n e r g yo fg a ， w h i l ea f f e c t sl i t t l eo nm gb e c a u s eo ff f si o n i z a t i o n f u r t h e r m o r e t h ea d s o r p t i o no f a c e t y l e n eo nt h en i ( 1 “1s u r f a c ei si n v e s t i g a t e d b yd v - x c tm e t h o d ，a sw e l la si t sm i g r a t i o np a t h w a yf r o mt h ea l i g n e db r i d g e s i t e ，t h ef a v o r i t ea d s o r p t i o ns i t ef o re t h y l e n e t o 也ec r o s s b r i d g es i t e a f t e r s t u d y i n gs e v e r a lr e a c t i o n p a t h w a y s o f c 2 h 2o n n i ( 1 1 1 ) ，w e f m d t h a ta c e t y l e n e ， a st h ep r o d u c to fd e h y d r o g e n a t e de t h y l e n e ，w i l l m i g r a t ef r o mt h ea l i g n e d b r i d g es i t et ot h ec r o s s - b r i d g es i t ew i t hav e r ys m a l le n e r g yb a r r i e r 0 0 2 e v t h ef a v o r i t ep a t h w a yo ft h em i g r a t i o ni sa 仃a n s l a t i o nt ot h en e a r e s tc r o s s - b r i d g es i t e w i t har o t a t i o nt h r o u g h 3 0 0s i m u l t a n e o u s l y 浙江大学博士学位论文 致谢 y 3 0 3 1 2 9 本论文是在我的导师曹培林教授的悉心指导、热情关怀和具体帮助下 完成的。导师渊博的知识、开放的思维方式、科学严谨的治学态度，给作 者留下了深刻的印象，并将对作者今后的学习和工作产生深远的影响。尤 其是导师平易近人，乐观开朗的生活态度，处处为学生着想的无私精神， 为我们创造了榆快的学习工作环境。能有这样的导师，我感到无比的荣幸。 在此，向我的导师谨致以诚挚的感谢和衷心的祝福! 感谢物理系的鲍世宁教授、唐景昌教授、张训生教授、徐伦彪教授多 年来对我的关心和帮助，是他们引导我进入表面科学领域的。特别是鲍世 宁教授，在本论文工作中给作者以具体的指导和帮助，在此致以衷心的谢 意。 感谢研究生院的周煦炎副教授，他几年来对作者的关心和帮助，使作 者能顺利完成论文工作。 感谢日本静冈大学的来关明博士在作者留目期间给予的生活和学习上 的帮助，特别是为部分工作提供了优越的计算机条件。同时也感谢m s u m i y a 博士对作者的论文工作的热忱帮助。 感谢本实验室的姜明、仇曼，几年来，与他们在生活，学习上的台作 愉快，特别在学术上的讨论，常常使作者受益非浅，衷心希望今后能有更 多的合作。 感谢我的父母、亲友在我近二十年的学习生涯中对我无私的关怀和支 持。 感谢所有给予我关心和帮助的老师、同学和朋友。 赵宇军 一九九八年 一二月 于求是园 第一章绪论 1 1 动力学过程第一原理研究的重要性和困难 表面动力学，包括表面重构，弛豫，分子吸附，表面复位( r e s t o r a t i o n ) ， 表面扩散，晶格生长等，涉及表面物理的各个方面，对基础理论研究和工 程技术的应用都有着极其重要的作用。在实验研究之外，分子动力学计算 机模拟是研究复杂的凝聚态系统的有力工具| 1 ，2 lo 这一技术既能得到原子 的运动轨迹，还能象做实验一样作各种观察。对于平衡系统，可以在一个 分子动力学观察时间( o b s e r v a t i o nt i m e ) 内作时间平均来计算一个物理 量的统计平均值，对于一个非平衡系统过程，只要发生在一个分子动力学 观察时间内( 一般为l l o p s ) 的物理现象也可以用分子动力学计算进 行直接模拟。可见数值实验是对理论和实验的有力补充，特别是许多与原 予有关的微观细节，在实际实验中无法获得，而在计算机模拟中可以方便 地得到。这种优点使分子动力学在表面动力学现象的研究中显得非常有吸 引力，因为原则上上面提到的现象都可以进行计算机模拟研究。 分子动力学假定原子的运动是由牛顿运动方程决定的，这意味着原予 的运动是与特定的轨道联系在一起的。当核运动的量子效应可以忽略，以 及绝热近似严格成立时，分子动力学的这一假定是可行的。绝热近似也就 是要求在分子动力学过程中，每一时刻电子均要处在相应原子结构的基 态。大多数情形下，这一条件都是满足的，要进行分子动力学模拟需要知 道原子间正确的相互作用势，从而必须知道相应的电子基态。电子基态的 计算是一个非常困难的量子多体问题。好在密度泛函的引入使这方面的计 算有了很大的简化，意味着我们可以把这一多体问题转化为一组自治的单 粒子轨道方程【3 ，4 】。对交换相关势采用局域密度近似，这组方程就实际可 解了【4 ，5 。这是目前凝聚态物理电子结构计算中普遍采用的方法。对手非 强相关联系统，局域密度泛函非常有效。在局域密度泛函基础上的第一原 理性计算，在研究原子、分子、和晶体的结构中取得了巨大的成功【5 】。 然而由于计算上的复杂性，基于局域密度泛函的第性原理计算长期以来 认为是不可能直接用于统计力学模拟，因为上百个原子的模拟就需要对 1 0 4 1 0 6 种原子构型作电子结构计算。因此，在分子动力学模拟中一般只 能采用经验势来代替原子间实际作用势。一般就用l c n n a r d - j o n e s 势之类 的二体作用势来替代。这种替代可能对惰性气体之类的系统进行模拟有 效，但许多实际材料科学中感兴趣的问题都是多体效应的结果，不可能光 用二体势代替能解决的。尽管人们对金属和共价系统的模型势作了许多尝 试f 6 ，7 】，并取得了一些成果，但一般一种模型势只能适用于一种原子，不 能满足实际模拟的需要。 另外经验势近似有一个本身的局限性，它丢失了局域电子结构之间存 在着的强相关作用信息，也就是说，不能得到成键性质( b o n d i n g p r o p e r t i e s ) ，以及原子动力学过程中的电子性质。尽管在某些特殊情况 下，电子性质也能通过经验势近似来得到，例如用经验势的分予动力学来 计算原子结构时，选取某几种原子构型来作电子结构计算 8 1 ，这样不仅 耗时，而且原子动力学和电子结构计算成了相互独立的过程。综上所述， 用第一原理来计算各种动力学现象是计算机模拟实验的一个必然趋势。 1 2 实现动力学过程第一原理研究的可行性 1 9 8 5 年c a r 和p a r r i n c l l o 在传统的分子动力学中引入了电子的虚拟动 力学，把电子和核的自由度作统一的考虑，首次把密度泛函理论与分子动 力学有机地结合起来，提出了从头计算分子动力学方法( 也称c p 方法) f 9 ，在不少体系中使基于局域密度泛函理论的第一原理计算在直接用于 统计力学模拟成为可能，极大地扩展了计算机模拟实验的广度和深度。 到目前为止，从头计算分子动力学方法主要是基于赝势和平面波实 现，尽管对赝势的构造和平面波的收敛性有作了改进，但对过渡金属元素 和第一行元素的计算仍然比较困难。有时对这些体系目前无法直接应用从 2 浙江大学博士学位论文 头计算分子动力学方法实现对其动力学过程的第一原理研究。面对各种复 杂的表面动力学现象，我们认为有条件对其作从头计算分子动力学模拟 时，最好用它来模拟，因为这种方法得到的信息最多。在条件不充足时， 可以通过抓动力学过程的主要矛盾，简化物理过程，甚至通过非动力学方 法坐标优化等途径来实现动力学现象模拟的第一原理研究。总之，在对不 同的体系进行第一原理计算时，需要用不同的方法以及针对性的简化处理 才能得到满意的结果。本论文的工作中分别就几种典型的系统，即硅表面 的吸附和表面复位，含第一行元素的i i i v 族半导体g a n 薄膜的m g 掺杂， 过渡金属n i ( 1 l1 ) 表面乙炔的吸附和反应途径等问题用不同的处理或计算 方法实现对其动力学过程的第一原理理论研究，解决了上述体系动力学过 程中一直没有解决的问题，同时说明随着计算机条件和计算方法的改进， 表面动力学过程的研究已经步入第一原理研究阶段。 1 3 本论文的安排 本论文系统地介绍了从头计算分子动力学方法，围绕对表面动力学过 程第一原理理论研究，分别对三种典型的体系，用不同的处理方法或计算 方法进行研究，很好地实现了这些体系动力学过程的第一原理理论研究。 第二章，我们较为详细地介绍了从头计算分子动力学方法的原理以及 基于范数不变赝势 1 0 ，l1 1 和平面波的实现。 第三章，介绍了用从头计算分子动力学方法模拟硫族元素在s i o o o ) 表面的吸附以及吸附引起的表面复位。讨论了吸附结构与覆盖率的依赖关 系，对这一体系以往实验为什么给出互不相同的结果作出了解释。 第四章，介绍了g a n 薄膜极性与m g 掺杂关系的研究，这儿我们抓 住在掺杂中起关键作用的吸附能，通过从头计算分子动力学方法对g a ， m g 在不同g a n 薄膜极性下的吸跗能的计算，并通过离散变分方法验证， 给出了m g 掺杂与g a n 薄膜的关系。 第五章，介绍了用离散变分方法，通过对乙炔所有坐标的优化研究乙 烯在n i ( 1 11 ) 表面脱氨解离的化学过程，讨论了乙，炔从乙烯的最佳吸附位 到乙烯脱氢后的最佳吸附位的化学反应途径。 最后，在附录中给出了超软( u l t r a s o 哟赝势 1 2 】的构造和从头计算 分子动力学方法基于超软赝势和平面波的实现。 浙江大学博士学位论文 参考文献 【l 】a r a h n l a n ，i nc o r r e l a t i o nf u n c t i o n sa n dq u a s i p a r t i c l ei n t e r a c t i o n s 加 c o n d e n s e d m a t t e r ，e d i t e db yj w h a n e y ( p l e n u m , n y 19 7 8 ) f 2 f o r ar e c e n to v e r v i e ws e ea l s o m o l e c u l a r - d y n a m i c s ? s i m u l a t i o no f s t a t i s t i c a l - m e c h a n i c a ls y s t e m s ，e d i t e db yg c i c c o t t ia n dw g h o o v e r ( n o r t h - h o l l a n d ，a m s t e r d a m ，1 9 8 6 ) 3 】p h o h e n b e r ga n dw k o h n ，p l a y s r e v 13 6 ，b 8 6 4 ( 1 9 6 4 ) 【4 】w k o h na n dl j s h a m ，p h y s r e v 1 4 0 ，a 11 3 3 ( 1 9 6 5 ) 【5 】o p s f i v a s t a v aa n dd w e a i r e ，a d v p h y s 3 6 ，4 6 3 ( 1 9 8 7 ) 【6 】f e r e o l e s s i ，m p a r r i n e l l oa n de t o s a t t i ，p h y s r e v l e f t 5 7 ，7 1 9 ( 1 9 8 6 ) 7 】m w f i n n i sa n dj e s i n c l a i r ，p h i l m a g a 5 0 ，4 5 ( 1 9 8 6 ) ： 【8 】j q b r o u g h t o na n dp b a l l e n ，m a t l e t t r e s ：s o c s y m p p r o c 6 3 ，1 3 ( 1 9 8 6 ) ( 9 】r c a ra n dm p a r r i n e l l o ，p h y s r e v l e r 5 5 ，2 4 7 1 ( 1 9 8 5 ) 【1 0 d r h a m m a n ，m s c h l u e t e r ，a n dc c h i a n g ，p h y s r e v l e f t 4 3 ，1 4 9 4 ( 1 9 7 9 ) 11 】g b b a c h e l e t ，d r h a m m a na n dm s c h l u e t e r ，p h y s r e v b2 6 ，4 1 9 9 ( 1 9 8 2 ) 1 2 d v a n d e r b i l t ，v h y s r e v b4 1 ，7 8 9 2 ( 1 9 9 0 ) 5 第二章从头计算分子动力学方法介绍 2 1 离子势的第一原理计算 2 1 1 条件与假设 第一原理推导离子间相互作用力主要基于两个假定。首先认为离子实 可以看作经典粒子来处理，其次假定系统满足b o r n o p p e n h e i m e r ( b o ) 绝热近似，即电子的量子行为与离子的经典行为具有相对的独立性。这就 要求电子都处在单能态，特别地，要求电子的基态是非简并的。对于半导 体与绝缘体，由于基态与第一激发态之间存在着一个较宽的能隙，这个能 隙一般大于离子运动有关的能量，所以b o 近似一般都能满足。对于金属， b o 近似的有效性就显得有点微妙了【1 。不过一般地，只要金属的等离激 元的振动能量远大于离子的热力学能量，我们也可以认为b o 近似对金属 适用【2 。 对于满足b o 近似的系统，离子问相互作用势v ( j ，) ) 可以从第一原 理导出： v ( r ，) ) = ( o 矧中o ) 。 ( 2 1 ) 其中h 是相应离子坐标下的系统哈密顿量，o 。为相应的电子基态。为了 在实际分子动力学中应用( 2 1 ) 式，由于模拟中离子坐标詹，不断变动，电子 基态中。的计算次数不会少于1 0 4 。 由上面的讨论可以知道，要使( 2 1 ) 式在分子动力学模拟中实际可行， 有两个最基本的要求。首先要有一套处理系统的多体哈密顿量h 的理论 方法，使得电子基态能量的计算在不同原子构型下有基本相同的精度，否 则由于分子动力学模拟过程中原子坐标变动很大，就会影响计算结果。另 外还需要一种当代计算机能够忍受的具体可行的计算方法，而不能在分子 动力学过程中，对每一种原子构型的电子基态都需要从头计算。 从第二条，即计算机条件来考虑，就把许多第一原理方法给排除了。 尽管它们在大量的原子结构计算中取得了成功，但由于其计算量一般随着 系统的电子数n 的增加按指数增长( o ( n “) ，m 4 ) ，即使是从头计算 h a r t r e e f o c k 方法，计算量也是以o ( n 4 ) 增长，更重要的是，这些方法对 不同原子构型的电子基态都需重新迭代计算，而在m d 模拟中，将要计算 1 0 0 左右的原子构型的电子基态，因此这些第一原理方法直接用于分子动 力学模拟是不现实的。十年前c a r 、p a r r i n e l l o 提出的从头计算分子动力 学方法【3 】，在实际的动力学物理系统中引入了电子的虚拟动力学，基本 上从方法上解决了这一问题，在介绍这一方法之前，在这一部分先扼要地 讲一下密度泛函理论和求解k s 方程的一些方法。 2 1 2 密度泛函理论和局域近似 对一个多原子系统，当应用b o r n - o p p e n h e i m e r 近似以后，问题化成 了一个多电子系统，即一个存在相互作用的非均匀电子气。1 9 6 5 年， h o h e n b e r g 和k o l l i l 等人建立的密度泛函理论【4 】不但为单电子近似建立 了坚实的理论基础，也为处理非均匀爹电子体系的基态性质提供了可行的 办法。而在此基础上k o h n 、s h a m 导出的单电子轨道理论更清楚地给出 了单电子近似的理论基础【5 。到目前，这一理论已被广泛应用，包括晶 体、表面、分子和原子集团等方面的计算。 密度泛函理论有两个基本假定：( i ) 对于一个费米子系统，任何一 个可观察置( 如h 。) 在基态的期望值都是电子基态p 。的一个泛函： 【户。1 皿j 【i d 。】 = 【p 。】 ( 2 2 ) 泛函是不依赖于外场的普适函数，也就是说，对于处在原子、分子、 以及固体等不同环境中的电子的密度泛函都是一样的。( i i ) 电子系统 的基态能量对应于能量泛函e ，( 力的极小值，即 e o = m i n e ， 纠( 2 3 ) 其中e ，【纠= 【纠+ f 布矿( f ) 烈芦) ，y ( f ) 为电子系统所受的外场( 如核对 电子的作用势) 。k o h n ，s h a m 单粒子轨道理论的中心假定是，对于任何 一个相互作用的电子系统，存在一个电子基态密度与之相同的无相互作用 粒子的等效系统( 等效系统的哈密顿量为：胃，= t + v ，总能为 e ；= 正+ i d r y 。p ，) 即p o = 风，这样p o 具有单一的表达式： p 。( ，) = ，i y ，p ) | 2 ( 2 4 ) 单粒子正交轨道妒。可由下面的s h r o d i n g e r 方程得到： 日n 妒，= ( 一 v 2 + p 0 ) y = s 。y ， ( 2 5 ) 本文中均采用了原子单位( c = 幌= a = 1 ) ，简单起见，电子占有数 在推 导过程中认为是相等的。这样外场v 作用下的相互作用系统的能量泛函 e v 可以写为： e p o 】= r 。 p o 】+ 缈( 力p o + ；f 黼眢岷【p o ( 2 6 ) 由( 2 6 ) 式可以看出，能量e 被分成四个部分，非相互作用系统的动能、 与外场的相互作用能、经典电子电子相互作用库仑能，以及包含量子效 应的交换相关能。相应的有效势f 0 为： 川卅p 高4 。丽a e o n 厂 ( 2 7 ) 方程( 2 4 ) ，( 2 5 ) ，以及有效势表达式( 2 7 ) ，就构成了k s 方法。 遗憾的是，它还是形式上的东西，因为包含电子多体效应的交换相关 部分e 。的具体表达式至今还不清楚。这样又不得不求助与某种近似，现 在广泛应用的就是局域密度近似，即把非均匀电子系统分割成一些小块， 在这些小块中，认为电子气是均匀的，这样尹处的子块中的交换相关能密 度s 。( ，) ，只取决于该点的厦芦) ，整个系统的交换相关能为： e 。= i 6 ：厦尹) p 。 厦芦) 】 ( 2 8 ) 至于交换相关势p ，的具体表达式不断有作者给出，目前用得比较广 浙江大学博士学位论文 泛的是c e p e l e y 和a l d e r 在八十年代用m o n t e c a r l o 方法导出的解析表达式 【6 】。从方法推导中的近似来看，局域密度近似只局限于那些电子密度变 化平缓的系统，但实际上，在一些不满足这个局限条件的系统中局域密度 近似也取得了较大的成功( 如具有很强方向性的半导体) 。总之，局域密 度泛函理论，自六十年代提出以来经过许多入的完善和发展，已广泛地应 用于各个领域的多原子体系电子结构计算。由于在此理论近似下单电子运 动方程中的交换相关势形式简单，使计算工作量大为减少，因此，现在大 多数能带或原子集团模型中的计算方法，都是在局域密度泛函理论近似的 基础上建立起来的。 2 1 3k o h n s h a m 方程的求解 在传统的固体、分子的局域密度泛函电子结构计算【7 中，k s 自洽方 程组一般是这样求解的：先给定一个初始的电荷风( 尸) ，譬如可以选用 各原子电荷密度的叠加，然后代入( 2 ，7 ) 式计算得到有效势，再通过对 角化来求解( 2 5 ) 式，得到矗。得本征矢，从而可以得到一个新的密度 分布函数风( ，) ，继续进行迭代，直至自洽。由于对角化的计算量以 m3 增长( m 为k s 轨道展开所需的基函数个数) ，这样对于较大的系统 工作量非常庞大。 k s 方程组也可以通过其它方法求解。为求得能量泛函对电子密度变 分的最小值，我们可以在电子自由度妒。的空间中引入一个虚拟的动力学 过程，来使泛函达到极小值。最简单的动力学考虑悬速降法【8 】( s t e e p e s t d e c e n t ) ： 脯 州 f ) _ 一葫蒜 2 9 ) 这里，点是代表对时间t 的微分，变分罟= 妒。，依赖于y ，。 o w i 求解( 2 9 ) 式时，为保证波函数的正交性，还得加上一个约束条件，虚 9 浙江大学博士学位论文 拟动力学方程就变为： 厢 一 少一( f ，) 一赢+ ；a ( 2 1 0 ) 1 一 j j 实际上，先用无约束的动力学方程( ( 2 9 ) 式) 计算妒，然后通过适当 的正交化过程对y ，进行修正，并得到系数人护当系统达到极小值时， 少。= 0 ，即z k 。= a f y ，约束矩阵a f 对角化后，就得到了k ，s ， 方程组的本征矢和本征值。经验表明。对于一定的原子坐标，e 只有一个 极小值，这样速降法也能得到最小值。值得注意的是，用这种动力学方法 自洽地得到能量最小值，与传统方法相比省去了大量的对角化工作。 速降法的有效性取决于波函数达到收敛的有效步数，i 这卞步数可能需 要很多，特别是在低对称性的情形下。人们也做了许多改善速降法的尝试， 如用更高次的积分方程，以及改进的速降法【9 】。摄近还发现采用共轭梯 度法来计算hk 能量泛函的极值还是比较有效的 1 0 1 。 前面已经提到，象速降法和共轭梯度这样的算法在哈密顿量宣。作用 到矿，的每步计算中都需要适当的正交化操作。如果妒是用平面波展开的 ( 设n 个单电子轨道均用m 个平面波展开) ，则日。y 的计算就需要n m 2 次浮点计算，加上正交化需要n 2 m 次操作，计算量很大，为减少计算量 我们可以把算符h k s 分成两部分，即动能项( 一 v 2 ) 习订势能项( 圪，) ， 动能项在倒格矢空间是对角的，作用到波函数| ；f ，。上只需o ( n m ) 次操作。 如果势能项k ，是局域的( 1 0 c a l ) ，那么作用至渡函数矿，上也可方便地 使用f f t 技术在实空间中计算，也只要o ( n m l n 次操作。因此随着m 的增大，计算量是以线性或o ( n m i nm ) 增长，这与传统的标准对角化技 术相比，工作量有显著减少。 1 0 渐江大学博士学位论文 2 1 4 局域密度近似下的离子闯相互作用势的计算 传统地，如果对给定的原子坐标 j i ，1 由从头计算方法推导离子间的 相互作用势矿( j ，) ) ，然后进行第一原理分子动力学模拟，一般有三个基 本步骤，即在每个m d 过程中：( i ) 对于给定的 j 2 ，) 求解自洽的k s 方 程，( i i ) 根据h e l l m a n f e y m a n n 定理【1 1 计算每个离子受到的力，( i i i ) 求解牛顿运动方程鸩豆，= 一v m v 。十年前，c a r 、p a r r i n e l l o 给出了一种 较先进的方法来进行从头计算分子动力学模拟，这种方法的重要特点是能 同时处理电子和离子的系统。也就是说用不着在每一步m d 中去重复自治 地求解k s 方程组，就能同时得到离子的轨迹和相应的电子基态。这个方 法以及全局优化技术将在下一章中详细讨论。在这之前，我们先来考虑怎 样同时对电子和离子的自由度进行优化，得到能量的极值。 同时对电子和离子的自由度进行优化的想法是由b e n d t 和z u n g e r 最 先提出来的1 1 2 。利用速降法，通过下面的方程组就可简单实现这想法： 死沙，= 一熹+ a 口妒j 1 ) v r ， j 1 厨，豆= 一娄 ( 2 1 2 ) 。 觎， 其中“质量”刀。和羁的引入可以调节与两套参数 妒。) 和 j 矗相关的时 间标度。方程组( 2 1 1 ) 、( 2 1 2 ) 就确定了能量面上从初始点到附近 极小值之间的一条轨迹。在 j i ， 不变时，通过变化 矿 ，就可找到能量 面e 上唯一的极小值，相应地，在b o 势能面上也有与之对应的一个极小 值：矿( 詹， ) = 耳卸研( 妒。) ， 孟，) 】，这样能量面e 上的局域极小值也就成了 t r ， 势能面v 上的局域极小值。通过同时对 妒j 和 夏； 进行变分，方程组 ( 2 1 1 ) 、( 2 1 2 ) 就可以用于复杂的原子结构计算，与传统的第一原 理计算方法直接用于分子动力学模拟相比，这一方法在计算量上就大为减 少。方程( 2 11 1 2 ) 决定的动力学虽然很简单，也成功地解决了许多复 杂原子结构的局域优化问题 1 3 - 1 5 】。 渐江大学博士学位论文 不过，方程( 2 1 1 1 2 ) 决定的动力学毕竟还有很大的局限性，它所 得到的结果只是局域优化的结果，而不是全局优化的。实际上，势能面v 一般有多个极小值。为了解决这个问题，我们将在下章介绍一个比较满 意的全局优化算法，即c a r 、p a r r i n e l l o 提出的基于统计力学的从头计算 分子动力学f 3 j 。 2 2 从头计算分子动力学( c p 方法) 的原理 8 5 年，c a r 、p a r r i n e l l o 提出的从头计算分子动力学方法 3 】，最重 要的一点是在真实的物理系统中引入一个虚拟的电子动力学系统。这样 组成的新系统的势能面e 是离子和电子自由度的一个总泛函，恰好这个 虚拟系统产生的轨迹与具有势能面v 的实际物理系统的轨迹非常一致。 实际物理系统的经典拉格朗日量为： = 鸩豆卜j i 朋 ( 2 1 5 ) ， 而虚拟系统的广义经典拉格朗日量为： = j 谚扩) 2 + m ，霹一研 晨， ， 妖 】+ 一 7 ( 2 1 6 ) 一 、7 乞a 口( j 方? ( ，) ( 尹) 一岛) f ，j 是二套自由度 y ，和 j l ，) 的泛函，它本身不显含时间r ，但；f ，和晨， 与时间有关。“为任意参数( 单位为：质量长度2 ) ，从式子中看， u ，相当于电子的“质量”，实际上它起着调节电子运动时间标度的作用。 为简单起见，我们可以令朋与电子波函数。无关，即2 。；，( 2 1 6 ) 式的拉格朗臼量中，第一项和第二项分别是电子和离子的动能疋和置， e 孟，) ， l l f ，。) 】是电子和离子耦合虚拟系统的势能，拉格朗日乘子a 。是为 了保证，的正交性而引入的，在经典力学中，它们就是一个完整约束。 1 2 浙江大学博士学位论文 从( 2 ，1 6 ) 式我们容易得到相应的e u l e r 方程为： 乃妒，= 一丽8 e + 莩a 鸬 砑，盖，；一粤 。 积， ( 2 1 7 ) ( 2 1 8 ) 我们注意到通过改变离子的速度豆，和电子的“速度”舻。，可以控制这一 虚拟系统的温度，从而可以对它们进行各种热力学处理，如退热、冷却 等。在这些热力学处理过程中，可以同时处理离子和电子的自由度，特 别地，我们可以设计这样一个热力学过程，即慢慢地冷却离子的温度， 最后使系统在t = 0 时达到其平衡状态，这一方法称为动力学退热模拟【3 】。 其实这只是k i r k p a t f i c k 等提出的退热模拟概念 1 6 】的特殊应用，所不同的 只是这儿用了m d 技术代替原先的m o n t ec a r l o 技术。 用动力学退热模拟来作全局优化计算的优点是明显的。在有限温度 下，各自由度的速度是不为零的，而且由于系统的复杂性，会产生随机 的热力学运动。由各态历经知，只要初始温度足够高，体系就能到达势 能面上的各个点。通过降温，体系在势能面上的随机运动慢慢局域于几 个极小值之间的跳跃，而且跳跃的频率逐渐减小。只要温度冷却得足够 慢，系统最终一定会平衡在势能面的最低点。值得注意的是，有些情况 下，这种全局优化