全文HTML
-
从微观尺度的原子分子到宏观尺度的宇宙天体,从无机矿物到有机生命体,物质世界可谓缤纷多彩,变化万千。然而大部分材料学、化学和凝聚态物理学覆盖的物质科学领域中,我们所关心的研究体系在原子分子层次的运动都遵循一个形式简单却又难于求解的基本量子力学方程:在非相对论情形下,这个方程是相互等价的薛定谔波动力学方程[1]或海森堡矩阵力学方程[2];而在相对论情形下,这个方程是狄拉克方程[3]。其中薛定谔方程的意义和重要性已广为人知,但只有结合全同费米子(如电子)满足泡利不相容原理,量子力学方程才能解释原子壳层的电子结构及化学键合、电子跃迁、光子辐射/吸收等现象[4]。
全同费米子间的量子关联使得求解电子的多体问题变得十分困难。电子与电子之间的库仑相互作用和全同粒子关联使得即使仅包含两个电子的简单体系(如氦原子或氢分子),严格的解析求解也不可能实现。量子力学变分原理为处理这一问题提供了另一种思路,即不再寻求问题的严格解,转而仅寻求系统能量最低时的数值/近似解。这从理想主义向实用主义妥协的办法所换来的回报是巨大的,我们获得了在单粒子薛定谔方程的框架下近似处理量子多体问题的能力。能量变分原理自然推导出一个等效的薛定谔方程,其中电子与电子间的直接库仑相互作用消失了,取而代之的是一个等效的电子所感受到的平均外场,即自洽场。自洽场是一种平均效应,能够近似描述多电子体系中电子所处势场的变化,却不能保证基于其获得的准粒子能级的准确性,但可以严格保证基态能量的正确性。早期变分原理的应用产生了Hartree方法[5]和包含全同粒子交换的Hartree-Fock(HF)方法[6],这为求解原子和分子体系奠定了基础,并最终发展出目前广为流行的密度泛函理论[7–8]。这些方法均属于平均场方法,由于方程中只出现单粒子占据的电子轨道,亦称单电子近似。本文将简要介绍基于变分原理的量子多体方法(特别是密度泛函方法),并对基于这些方法的第一性原理计算与模拟及其在复杂材料高压性质研究领域中的部分进展进行简要评述,最后展望这一充满活力领域的未来发展趋势。
-
虽然平均场理论的发展极大地简化了多电子体系的理论处理,但在应用到大体系中仍十分困难。主要障碍在于:包含N个电子的波函数是一个3N维空间中的复变量,对其进行直接数学处理的计算量十分巨大[9]。为了解决上述问题,Kohn和Hohenberg[7]建议对电荷密度而非波函数进行变分,其优势是可将3N维的物理问题简化为三维问题。为此,证明了基态能量是电荷密度的泛函,且外势场唯一地确定了基态电荷密度分布,反之亦然,至此严格地建立了电子多体问题的密度泛函理论(DFT)[7]。结合平均场思想,并引入在Hartree-Fock方法中被忽略掉的电子关联贡献,自然地获得了等效的单体薛定谔方程—Kohn-Sham方程[8]。Kohn-Sham方程的求解比Hartree-Fock方程的计算量小很多,精度更高。虽然DFT在原理上是严格的,但由于包含平均场近似下未知的有效势,在具体操作层面上,DFT不可避免地引入相互作用势近似[8–12]。在Kohn-Sham框架下(等价于将一个多体相互作用系统映射到在有效外势场中运动的理想费米气体系统),该有效势被分解为电子间静电相互作用的Hartree项、电子关联和交换项。需要指出,其中的交换-关联项还包含因Kohn-Sham单电子轨道映射而被忽略的动能贡献,而非纯粹的“交换-关联效应”。DFT的各种近似处理事实上是对交换-关联项的不同处理,包括局域的LDA近似[8]、半局域的包含一阶电荷密度梯度的GGA近似[10]和二阶电荷密度Laplacian的meta-GGA[11–12]、非局域的针对满壳层体系和软物质的范德瓦尔斯vdW泛函[13]、为修正自相互作用和带隙而引入部分轨道HF贡献的EXX和杂化泛函等[14–17]。此外,结合特定物理模型,还发展出针对电子强关联体系的DFT+U方法[18–19],以及其他的高阶方法,如GW方法[20–22]、ACFDT方法[23–24]等。
上述方法的发展逐步提高了DFT的精度和可解决问题的范围,使DFT成长为目前最流行的第一性原理计算手段。可利用DFT解决的物理问题主要有:总能、电子能带结构和电荷/自旋密度分布等,以及与之相关联的粒子和能量输运问题。其中最为成功的当属总能计算以及由能量计算衍生的结构优化、力学性质、原子运动和晶格动力学计算等。高压下材料的大部分响应和性质均可由总能计算确定,其中最为典型的即为状态方程。作为材料热力学性质中的重要内容之一,DFT方法最早应用于计算材料的状态方程。对于单质材料来说,这一应用已经十分成熟:零温压缩曲线是总能计算最简单直接的应用,被广泛用于研究轻元素材料(氢、氦、锂等)、简单金属与过渡金属、双原子分子气体和惰性气体等材料的压缩行为[25–29]。由总能计算还能得到每个原子所受的力(通过Hellmann-Feynman定理),从而计算出原子的运动和声子谱,结合晶格能量和声子自由能,利用DFT方法很容易得到有限温度下单质的压缩曲线、相图和热力学性质等,并在轻元素和典型金属等材料中得到重要应用[30–34]。目前,除了少数的几个电子强关联体系(这类特殊体系的处理需要用到动力学平均场理论)[35–36],利用DFT方法计算单质和化合物的热力学性质和状态方程不再存在原则性困难。
一个往往被忽略但很基本的问题是:利用DFT及其他第一性原理方法计算材料的热力学性质和状态方程存在重要假定,即材料的成分、原子位置和分布是确定的。对于单质材料和化合物,这个假定是成立的;但对于合金以及一些非标准化学比化合物来说,这种假定往往不能满足[37]。在这类复杂材料中,材料的化学成分、原子位置和分布是不确定的,而DFT只能对其中的单个构型(Configuration)进行计算,不能获得体系的总体物理性质。因此DFT的直接应用失效了,这类复杂材料状态方程的第一性原理计算成为一个很大的挑战。
-
解决材料化学成分、原子位置和分布的不确定性困难,关键在于统计分析化学成分的可能变化和原子位置的可能分布,因此需要将DFT计算与统计力学相结合。最简单的统计模型是无相互作用的理想气体模型,相应最简单的合金状态方程即为混合物状态方程[38]。顾名思义,混合物模型简单地将不同元素的单质状态方程按化学比例均匀混合(线性叠加),忽略了不同种类原子间的相互作用和混合熵的贡献。从实验的角度,另一个简单的处理方法是将合金当作一个等效的单质或化合物,并忽略原子的不同占位和组态熵的贡献,这种方法在DFT计算中存在原则性困难,因为很难寻找到一个组态,使其他的状态方程与合金的状态方程在物理上“等效”。图1(a)给出了采用第一性原理计算Ni-Al合金在不同组分和不同原子占位情形下的压缩曲线,其中明显的结构相关性源自不同原子间的相互作用,揭示了混合物模型和等效单质/化合物模型的局限性。很明显,即使在相同组分下不同组态(即不同原子占位情况)的可压缩性差别也很大。图1(b)给出了实验结果与混合物模型的比较[39],混合物模型完全偏离了实验的状态方程数据,二者间不可忽略的偏离说明混合物模型或等效单质模型的适用范围十分有限。更重要的是,在有限温度下组态熵的贡献变得越来越重要,但在混合物模型或等效单质模型中却被完全忽略了。事实上,大多数合金的高温高压行为完全不同于混合物,这进一步凸显了建立完整的合金状态方程理论的必要性和重要性[40]。
合金体系材料物性第一性原理计算的所有困难均汇聚于如何构建一个好的统计物理模型,从而使其中的构型分布能够准确地描述真实体系中原子的占位分布情况。对于替位合金体系,基于格子气模型(Ising模型),Kikuchi[42]提出了最早的合金热力学性质的计算框架和组态熵的解析计算方法—集团变分法,并推动了合金相图研究进入了定量描述的崭新阶段。Kikuchi的最初方法建立在格点几何结构的物理分析之上,组态熵的导出步骤较为复杂,很难推广到任意格子体系。集团变分方法的数学基础—包括其中作为展开基矢完备性和组态熵表达式核心的Mobius变换—直到20世纪70年代才完全建立起来[43–47]。集团变分法的建立使得计算合金组态熵不再是一个困难,但第一性原理方法在这个问题上的应用还是障碍重重—人们不知道如何去计算合金的内能贡献,因此当时流行的合金相图与物性的理论研究是结合拟合实验或半经验的内能项和集团变分的组态熵项而获得的,虽然得到的理论相图与实验结果符合得相当好,但作为一个理论框架其预测能力和外推能力显然是残缺的[48–49]。20世纪80到90年代的一个进步是研究人员逐步意识到对组态熵的集团展开同样也适用于能量[50–52]。由于在有限格子上展开基矢的完备性,事实上集团展开也适用于定义在格子上的任意函数(例如Garbulsky等[53–55]提出单独对振动自由能进行集团展开,Geng等[56]提出对热电子的贡献部分进行单独展开等),而且由于基于格点占位的合金Ising模型的自身特征,这种展开是有限且快速收敛的[47, 57]。例如FCC晶格中的四面体展开或四面体-八面体展开就能给出很好的结果,如图2所示。自此,能量的集团展开方法获得快速发展,扫平了第一性原理方法在合金体系中应用的最后一道障碍。采用第一性原理方法对几个特定组态的总能进行计算,再利用集团展开方法,我们就能得到相同格子上任意组态的内能,结合集团变分法,自然地得到热力学平衡条件下的自由能[58],进而给出第一性原理的合金相图(参见IMR-CVM程序)。整个计算过程不依赖任何的经验参数或实验值,仅有的假定是格子气模型,其精确度取决于集团展开的收敛程度和第一性原理总能计算的精度,是一个典型的从头算方法。除了结合集团变分法计算组态熵,利用Monte Carlo模拟亦可确定原子占据组态的分布情况,再结合热力学积分即可得到组态熵,沿此路线,Walle等[55]开发了能量集团展开和原子占据组态分布的Monte Carlo模拟程序包ATAT,推进了合金性质和相图的第一性原理计算模拟。
结合集团展开法和集团变分法(或利用Monte Carlo模拟计算组态熵),基于第一性原理方法计算零压下合金及金属间化合物的相图不再是一个难题:自由能在化学组分空间的极小值给出了体系的能量凸图(Convex Hull),进而确定了可能出现的合金相,而其拐点则给出了旋节线分解点和失稳区域[59]。在有限压力下,该方法需要进行适当拓展才能应用。组态熵部分在形式上与压力无关,但能量和振动熵则随压力而变化,因此需要引入一个新的维度。Sluiter等[60]对此进行了直接而自然的拓展,在假定格子类型不发生改变的前提下,将集团展开的系数(即集团相互作用强度)从原来的常数拓展为温度和体积的函数,以描述有限压力下单个给定构型的自由能对温度和体积的依赖关系。其问题的核心在于,由于最终获得的模型需要应用到有限的压力范围,如何保证在整个应用的压力区间集团展开的收敛性均满足要求。为此,Geng等[41]阐明了只要在最小体积(即最高密度)下集团展开是收敛的,则通常在更大的体积(即更低密度)下同一个集团展开必然也是收敛的。该法则基于相互作用随距离增大而减弱这一物理事实,简化了对不同体积(密度)下进行集团展开收敛性的检验,促进了这一方法的发展,使得计算有限温度下的合金相图有了可靠保证[61]。
由于获得了完整的密度(体积)-温度空间的自由能,该技术进步还使得计算不同化学组分下的状态方程和组态熵的贡献成为可能:在给定化学组分下,自由能随密度变化的梯度给出了对应的压力[41]。在此之前,合金状态方程的计算主要基于理想混合物模型[62],既忽略了不同成分间的相互作用,还完全不考虑组态熵的影响。图1为Geng等[41]以Ni-Al合金为例,展示了不同成分间的相互作用对合金冷压曲线的影响可达10%以上,不能忽略。此外,还详细讨论了组态熵和振动熵的相对大小,并发现高温下组态熵驱动的有序-无序相变导致在冲击Hugoniot曲线上出现明显拐折[63],如图3所示,其中在低温下有序态和亚稳的无序态分别有不同的Hugoniot曲线,相变点处的体积跃变是一级相变的典型特征。这一现象在以前的理论和实验研究中均未被讨论过,是合金体系在动态压缩下冲击行为的一个新的理论预测。需要指出,合金有序-无序相变是一种典型的扩散型相变,与位移型相变相比,其时间尺度慢了几个量级,大约在微秒至亚毫秒时间尺度,甚至更慢。一般平板撞击实验的冲击波平台时间尺度为亚微秒,激光冲击实验的持续时间更短,因此很难观测到这种较慢的动力学转变。然而炸药爆轰和斜坡加载的时间尺度大约在微秒至几十微秒量级,可以与这种合金相变动力学过程耦合,从而导致复杂的物理和力学行为。这个主题与合金材料的冲击动力学和损伤破坏行为密切相关,在工程应用中有重要价值,但相关研究目前还处于起始阶段,还有大量的实验和理论工作需要开展。例如,合金在有限温度-压力下有序-无序相变中成核生长的基本物理参数仍未知,不仅不知道相变的势垒高度(因而无法估计出相变的孕育时间),甚至连原子的扩散速率和路径、机制均未知,极大地制约了相关动态实验的发展。
合金化对体系热力学性质的影响不仅体现在晶格有序占位相变和有序-无序相变上,在同一个有序相内,也会出现极为反常的热力学性质变化,而这些行为几乎无一例外地与组态熵相关。在FCC晶格的Ni-Al体系中,Geng等[61]的理论计算显示大部分的热力学性质(比热、热胀系数、Grüneisen系数等)均会在有序相的标准化学比处展现一个尖峰,并随化学成分变化展现“W”形(或其他奇特形状)的特征,如图4(a)所示。这种独特的合金行为是由组态熵驱动的在标准化学比附近的错位占据(Anti-Sites)直接导致的,是合金体系的指标式特征之一,在混合物模型中无法描述这种现象[62]。它的一个物理推论是,如果合金发生了偏析,则不仅会出现化学成分空间上的不均匀性,由于“W”形特征的存在,热力学性质的不均匀性有可能更显著。在零压下比热的这种反常变化已经得到了实验证实[62],但高压下的“W”形特征以及与状态方程密切相关的Grüneisen系数的反常行为还有待实验的进一步证实。理论分析还指出,这种反常行为不仅出现在晶格热振动的Grüneisen系数上,电子Grüneisen系数也会出现反常的变化[56]。
-
上述关于合金热力学及状态方程的理论框架构建在给定的格点上(格子气模型),最大缺点是无法直接计算不同类型格子上有序相之间的相界。例如Ni3Al和NiAl3是稳定在FCC格子上的有序相,而NiAl-B2结构则稳定在BCC格子上,无论集团展开、集团变分或是Monte Carlo模拟均无法直接处理这种跨格点类型的相图计算。虽然可以计算出在各自格子上自由能随化学组分的变化,并确定出哪一个相是热力学稳定的,但各相间过渡衔接区域的处理却一直是一个难题,强行光滑连接方法引入了一定的人为性。此外,利用各自格子上的自由能拐点给出的旋节线分解存在高估各有序相稳定区域的可能。如何解决跨格子的集团展开和组态熵计算问题可以说是合金模型目前面临的最大挑战。对此,Geng等[65]于2006年提出利用长程的对相互作用势描述晶格畸变能的扩展集团展开方法,并将局部畸变的格子映射到理想格子上,以保持多体的集团展开和组态熵形式不变,从而放宽了格子必须是理想晶格的严格要求。更进一步,他们还提出了扩展格子的方法,即将给定的格子进行扩展,以包含其他的格子类型[66]。例如:假定初始格子是BCC,在其6个面的面心引入新的格点,将获得同时包含BCC、FCC和简单立方的拓展格子。在这一框架下,必须引入白原子(即空位)以实现表示不同的格子:若体心格点被白原子占据将给出FCC格子,若面心格点被白原子占据将得到BCC格子,若体心和面心格点同时被白原子占据,则导致简单立方格子,如图5所示。在扩展格子上,由于不同类型的格子共存,必然导致大的晶格畸变,因而必须利用上述的扩展集团展开法来描述晶格畸变[64–65, 67]。这种方法还处在发展的初期阶段,尚未有实际材料的应用报道[66]。其面临的主要困难:(1)N元合金中,由于白原子的引入,原本的N元Ising模型将变为N+1元Ising模型,使得计算更为复杂;(2)由于白原子对能量没有贡献,大量白原子的引入将降低集团展开的收敛性,但这尚缺少系统的分析研究。基于扩展格子的扩展集团展开方法虽然建立在严格的物理原理之上,但用它来解决实际的物理问题还有很长的路,还有大量细节问题需要解决。首先需要系统评估白原子对集团展开收敛性的影响,给出最佳的集团基矢;其次需要分析研究扩展格子上的晶格畸变情况,找到合适的长程相互作用,以精确描述晶格畸变能;再次,多元体系的集团展开非常繁琐复杂,需要开发专门的计算程序,而扩展集团展开方法计算程序的欠缺,是导致这一方向发展迟缓的一个主要因素。
利用现有的扩展集团展开计算方法,结合格子气模型概念的变体,也推动了这一领域的发展。白原子和扩展格子的引入使得处理间隙原子和空位变得十分自然,从而将原本只适用于替位合金的集团展开理论自然地推广到间隙合金和含缺陷的更一般的复杂系统。Sluiter等[68]采用了扩展集团展开的一个简单近似,只引入了原子扩散路径上的一个代表点(一般为晶格的间隙位)以及与之紧邻的一个空位(即白原子),即ABC展开,这一处理简化了扩展集团展开的复杂性,但同时保留了其中物理的基本特质—原子扩散过程的过渡态也被近似地纳入到集团展开中,从而允许利用Kinetic Monte Carlo来描述合金中原子的扩散迁移行为,这在原格子气模型中是完全做不到的。更进一步,Geng等[67]利用扩展格子气模型,将点缺陷(空位和间隙原子)纳入到基于扩展集团展开的合金理论中统一描述。为了避开复杂的集团展开计算,他们采用了点近似,即假定间隙原子或白原子的占据是均匀无序的,而缺陷间的相互作用通过平均场的形式近似引入,从而忽略了这种作用的距离依赖性。虽然这一近似很简单,甚至有点粗糙,但仍反映了含缺陷体系与理想化合物明显不同的物理行为。以非标准化学比化合物为例,他们的计算显示缺陷可导致标准化学比附近反常的热力学性质变化,十分类似于合金中的“W”形反常,如图4(b)所示。基于这一近似模型,他们还获得了含缺陷体系的自由能,导出了相应的状态方程,并发现非标准化学比化合物允许存在一个十分独特的物理现象,即在有限温度压力下其中占主导的缺陷类型可能发生变化,从而导致自由能在两种缺陷间过渡,并在热力学性质上呈现反常的大幅度改变,甚至在体积(密度)上亦可表现出表观的“坍缩”行为(实际仍是连续的,只是在过渡区的变化较大,从而看起来类似跃变)[69–72]。这种“赝相变”现象是这类体系不同于理想化合物的一个重要特征,其幅值与化学比的偏离程度紧密相关。与合金的有序-无序相变一样,非标准化学比化合物的这一赝相变转变受限于原子扩散速率,是一个相对慢的过程,因而在较快的动态冲击实验中难以观测,但可以与相对慢的动态过程耦合,如很长的斜波加载或冲击波在材料中的多次反射,以及相应的加-卸载过程等,从而导致复杂的材料动态响应行为,这一领域目前尚属未开垦的处女地。
-
几乎在所有的量子力学计算中均存在这样的隐含假定,认为只要能将Hamiltonian对角化,则所获得的解便是物理的,其中能量最低的解即对应于基态。在平均场理论和DFT理论中,Hamiltonian是通过能量极小化的方法获得,因此该隐含假定可以解释为:在能量泛函的确定过程中,优化算法给出的残差最小的解即为基态解。诚然,若Hamiltonian能被对角化,相应的解必然为该Hamiltonian的本征态,其完备正交性保证其中能量最低的解必然为该Hamiltonian的基态解。但在平均场理论中,为了处理电子间的多体关联,通过能量极小化得到的Hamiltonian并不是完全确定的,而是依赖于电子的密度分布(即Hamiltonian的解),因此必须引入一个自洽的循环迭代过程,使得进入Hamiltonian的电子密度分布与该Hamiltonian的解完全一致[73]。整个求解过程由Hamiltonian对角化和对Hamiltonian进行修正(密度混合)两部分组成。这一自洽循环的过程在数学上等价于全局优化问题,即寻找能使能量泛函全局最小的电子密度分布[8]。
原理上,这种方法并没有大的问题,但在技术层面上同时继承了优化问题的困难和弊端:目前还没有一个行之有效的完善的全局优化算法,通用的方法基本上都是局部优化算法。事实上也确实如此,目前第一性原理计算程序采用的各种方法中,包括最陡下降法、共轭梯度法、残差法(RMM-DIIS)等均为局部优化方法。据我们所知,还没有一个第一性原理计算程序尝试过具有全局优化特性的算法,因而其获得的结果有落入局部极小的可能。这些局部极小解满足Hamiltonian对角化和电荷密度自洽这两个先决条件,但对应的能量只是局部极小而非全局最小(可视为Density-Driven Error的一种)。在通常情况下,特别是主族元素和金属中,价电子充分展开,所对应的能量曲面较为光滑平坦,很少遇到电子陷入亚稳态的情况。但在窄带体系和电子强关联系统中,能量曲面变得十分复杂,电子陷入亚稳态的几率大大增加。一个典型的例子是由5f电子主导的锕系化合物,最近的研究发现:利用DFT+U方法提高f电子的局域化程度后,电子亚稳态问题成为困扰相关材料第一性原理计算的重要障碍[74],这个问题在其他窄带体系中也可能存在。
对电子强关联材料的微观描述需要引入半经验的Hubbard在位库仑排斥修正,以此为基础的DFT+U形式的密度泛函理论在描述此类强关联电子系统上取得了很大的进步[18–19]。伴随而来的代价是同时引入了能量曲面上的大量极小值,这干扰了能量最小化算法的运行,使得电子结构优化过程常常终止在非物理的亚稳态。对局域的f轨道的占据状态进行干预(OMC方法)可以部分解决这一问题[75],即适当地改变占据矩阵的初始值来搜寻能量最低的电子态,但是这在大多数的密度泛函计算软件包中很难进行相应的操作,而且该方法是一种典型的试错法,并不能保证所获得的解就是全局能量最低的态,除非对所有可能的占据状态进行遍历。OMC方法的另一个困难在于所需要的计算量。如果不考虑自旋自由度,将n(n≤7)个电子填入7个f轨道上有
$ C_7^n$ 种可能的方式(只考虑对角情形)。非对角的占据数远大于此,但我们可以粗略地假定其只是对角情形下的数倍,即认为其他的那些占据状态并不重要。即便这样,每个原子将有$ mC_7^n$ 种不同的方式占据f轨道(假定m<10)。这只是在所模拟的原胞中只有一个不等价的重金属原子(最简单的)情形下所需要的重复计算次数。不幸的是缺陷结构具有很低的对称性,往往包含多个(例如k个)不等价的含部分局域的f轨道的重金属原子,在这种情形下$ (mC_7^n)^k$ 个独立的电子结构计算将是必须的。以二氧化铀为例,n等于2,m取值3,因此对每个不等价的铀原子需要搜寻大约60个不同的占据状态。如果考虑的是点缺陷,那么至少有两个不等价的铀原子,则总的计算次数将增加到3600。对于缺陷团簇来说,k一般大于3,因而需要进行超过百万次的电子结构计算才能获得最终的结果,即便是通过目前性能最好的超级计算机实现计算,仍然是个沉重的负担。为了解决这一困难,我们注意到对于一个经典的系统,类似的亚稳态问题可以通过模拟退火方法获得满意的解决,即缓慢地将粒子动能(对应于热运动)从系统中移走,从而使粒子逐渐收敛到基态结构上。我们发现类似的概念也可以应用到电子系统,基本思路是通过利用虚拟的能量振荡摇晃或者加热电子系统,使得有更高的几率跃过亚稳态间的势垒,从而达到获得真正基态解的目的,我们把这一方法称之为“准模拟退火”(Quasi-Annealing,QA)技术。这等价于在局域优化算法的基础上增加对全局最小的搜寻,使电子能够跳出局部极小陷阱。
该方法的理论基础在于电子系统的总能是电荷密度
$ n(r)$ 的泛函,其反过来又是由外部势场$ v(r)$ 完全确定的泛函[7]。因此可以利用粒子系统将虚拟的扰动经$ v(r)$ 传递给电子系统,总扰动为$ \int n(r) \Delta v(r) {\rm d} r$ ,其中$ \Delta v(r)$ 为外部势场的振荡。通过缓慢地移除这一虚拟的扰动,我们可以像经典的模拟退火方法一样将电子结构逐步收敛到基态。另外,可以通过跟踪电子优化过程来理解准模拟退火方法中移除亚稳态的机制:外部势场的振荡反复不停地改变极小化算法的路径和所有势垒的形状和高度,这样任何一个势垒都有可能被绕过。在实际操作中,可以利用能量优化过程中的“多余能量”模拟虚拟扰动:不收敛的自洽场导致的能量不确定范围δE会引起粒子所受力的准随机振荡,而这种力的不确定性将使得粒子在其相应的物理坐标附近形成类似于高斯分布的位置分布;反过来粒子运动的这种非物理漂移将改变电子系统所受的外场并导致δv的势场振荡,这最终提升了电子系统的总能量,即被“加热”了。通常可以用标准的结构优化算法来产生δv,我们可以在电子自洽场的冗余度为δE的条件下反复地弛豫原子结构。在这一实现中,唯一的可调参数(电子自洽场的剩余能量)控制了电子系统中的虚拟噪声。他的值在刚开始时应该足够大,这样电子系统才能在相空间中自由移动,但又要足够小以保证算法不发散。通过缓慢地降低δE的数值,可以把电子系统优化到基态。
准模拟退火的优点不仅在于它可以解决亚稳态问题,通过与粒子结构优化计算相结合,还可以同时优化电子结构和原子几何构型,因此大幅降低了计算所需的资源。这一特性在自洽场收敛特别慢的场合(例如DFT+U)中显得十分的重要。为了避免粒子偏离目标构型太远,需要在几个或者几十个粒子步长后将结构恢复到目标状态(或初始状态)。通常如果允许元胞体积和形状改变将提高准模拟退火的效率,因为这不仅扩充了搜寻的相空间区域,同时还对系统有整体的摇晃作用,并移除了Bravais格子的对称性限制,而后者往往是强关联系统优化过程中电子陷入高能量的亚稳态的主要因素之一。
准模拟退火的流程可以总结如下,其中每一步的计算需要从上一步产生的波函数开始:(1)移除所有对称性,设置适当的自洽场收敛冗余度δE和粒子移动步长δr;(2)利用标准的结构优化算法来演化结构;(3)适当减少δE和δr,重置结构,返回第2步,重复整个过程直到δE小于预设的收敛精度;(4)执行一个标准的电子结构优化过程,获得完全收敛的波函数,至此一个完整的准模拟退火已经完成;(5)对初始结构进行随机扰动,返回第2步,重复整个过程直到没有新的更低能量的态出现。这里执行第5步的目的是为了对前面的结果进行检验/验证,以免一些参数的不当设置降低了算法的效用而陷入亚稳态。
为了验证QA方法的有效性,我们首先讨论其在完整的二氧化钚晶体中的表现。其原子结构为包含12个原子的氟化钙立方胞,反铁磁相结构为1 k。交换关联能量泛函采用GGA+U及PBE近似,赝势采用PAW赝势,Hubbard模型参数U取4.0 eV,而J为0.7 eV。平面波基矢的截断动能设为400 eV,同时用63个不可约的k点对布里渊区进行采样积分。晶格常数固定在5.45 Å。选择二氧化钚是因为该系统具有十分稳定的亚稳态(相较于二氧化铀而言),例如如果对该结构施加相应的立方对称性(Oh),那么无论初始状态如何,总会陷入到总能为−125.282 eV的高位亚稳态之中。这完全不同于二氧化铀,那里对称性导致的亚稳态很容易被移除[67]。在二氧化钚中,只有把所有对称性都去掉才有可能克服该亚稳态。如图6所示,去掉对称性使能量降低了1.69 eV,但由此得到的态并非基态。绝热缓慢地引入Hubbard在位相互作用,即十分缓慢地从零增加U和J参数的数值(即后来美国西北大学倡导的U-Ramping方法[76])进一步降低了总能大约0.16 eV。而准模拟退火方法则给出了更低的能量。我们共进行了6轮独立的准模拟退火优化,所得的结果十分接近,标准偏差
$ \sigma $ (数据的分散性)为0.006 eV。另一方面,直接的电子结构优化给出的标准偏差比准模拟退火的大2~3个数量级,说明利用准模拟退火方法可以有效地剔除电子系统的高能亚稳态,并将之收敛到基态,所获结果的分散性大幅降低,提高了计算结果的可靠性,为开展精确而可靠的重金属元素及其化合物的性质计算提供了坚实的技术支撑。目前消除电子亚稳态的方法除了QA方法,还有法国原子能委员会(CEA)提出的OMC方法[75],以及美国西北大学改进的U-Ramping方法[76]。OMC方法对电子状态进行严格调控,理论上可以获得真正的电子基态,但在大体系、低对称结构中可能存在效率下降问题;Geng等[74]提出的QA方法是一种普适方法,对结构无特殊要求,在大体系、低对称系统中有独特优势,但单次模拟不一定能获得真正基态,需要多次模拟进行交叉验证。U-Ramping方法实质上是对DFT+U方法中的U值进行缓慢调增的Adiabatic Switch方法(图6),易于操作、使用方便,但对部分体系不一定能获得真正的电子基态。
-
我们对基于第一性原理的量子力学计算与模拟在复杂材料体系高压性质研究中的应用进行了简要回顾,其中特别强调了在合金、含缺陷材料以及强关联材料等复杂体系中,结合集团展开、格子气模型、准模拟退火等物理模型发展的第一性原理研究方法,而未对直接的第一性原理计算的应用(如总能计算、压缩曲线和声子谱的计算等)予以特别强调。这并不是说这些应用不再重要,相反,在地球、行星物理以及材料冲击动态响应的工程实践中有着巨大的应用价值,只是这些方面的应用和计算方法已经广为熟知,在此并未给予详细的阐述。
基于密度泛函的方法已经取得了巨大成功,推动并加速了凝聚态物质领域的发展,特别是材料设计方向的快速发展。DFT的缺点显而易见:一方面,其精度不足以描述精细的物理化学过程;另一方面,它所允许的原子分子体系尺度限制了其在大尺度复杂体系中的应用。针对精度不足的问题,目前已经开发了更精确的计算方法,并逐渐成熟,其中除了化学领域广泛使用的高阶量子化学计算方法外,针对扩展体系的基于DFT的微扰方法、准粒子方法以及量子Monte Carlo方法等也在迅速发展;针对尺度限制问题,为了提高计算效率,半经验的紧束缚方法以及等效Hamiltonian方法正在逐步复兴。DFT成功的关键因素在于,相较其他方法,该方法在牺牲了部分理论精确性的前提下大幅提升了计算效率。这种模式同样适用于紧束缚方法、等效Hamiltonian方法,甚至原子间相互作用势方法:在适当牺牲理论的严格性后,这些古老的方法有望在复杂体系的介观尺度(甚至宏观尺度)的模拟中得到复兴。特别是将这些方法与第一性原理计算相结合发展起来的多尺度计算方法,其通过在涉及原子分子尺度行为的区域采用计算精度较高但相对昂贵的DFT计算方法,而在仅关心宏观性质的区域采用连续介质模型,可以用接近DFT方法的精度模拟宏观大体系的各类性质,这一新兴方向在位错成核、裂纹传播[77]、化学催化机理[78]等多尺度现象的研究中已经展现了强大的发展潜力,将是未来第一性原理方法发展的主流方向。
纵观自然科学的发展史,我们已经发展到了这样一个路口:对控制单个粒子或少量粒子的自然规律已经非常清楚,但对大量粒子综合体所形成的集体行为的理解却很不充分。基础和前沿物理的进一步发展仍将会改写我们的世界观和宇宙观,但短期内从电子、夸克等基本粒子到分子、凝聚介质,直至星球、星系等这些尺度范围内的基本物理规律将很难再有所突破。可以这么说,目前除了继续探索超越标准模型的高能现象和宇宙规律外,自然科学发展的重心将逐渐转到利用已知的基本物理规律,认识和重构我们所处的复杂现实世界中来。基于第一性原理的量子力学计算和模拟,以及基于此衍生的一系列方法,使我们向这个宏伟目标迈出了很大的一步。随着超级计算机硬件和计算算法的快速发展,这种在虚拟世界中计算、模拟、设计并重构现实世界的潮流将不断从凝聚态物理、材料科学拓展到化学、医学、生命科学及其他领域。本文简述的部分内容,仅是从简单体系向复杂体系迈出的一小步,亦或是这一快速发展的洪流中一个微小涟漪,但都具有一定的代表性,希望该总结和回顾能对发展更先进、高效的具有预测能力的多尺度方法提供有益的参考。