说明:本文华算科技主要介绍原子间势函数的能量定义、势能与受力的关系、双体和多体作用,以及势函数对材料性质的描述范围。


势函数的定义:从原子构型到势能
给定一块铜晶体中全部原子的位置,拉伸晶格、移走一个原子或改变表面排列,体系的能量都会变化。原子间势函数 U(R1, …, RN)用原子种类和坐标计算这组构型的势能,Ri 是第 i 个原子的位置。本文讨论材料模拟中的这种势函数;晶胞、体系电荷等条件在比较时也须确定。
在玻恩—奥本海默近似下,固定原子核坐标并求解所选电子态,电子能量与核间排斥能之和形成原子核运动的势能面。DFT 可以为不同构型提供能量和力,原子间势则用解析表达式、数值表格或拟合模型近似这张势能面,计算每个新构型时通常无需重新求解电子波函数。

对于孤立双原子,整体平移和转动不改变内部能量,势函数可简化为只随间距 r 变化的 U(r)。图1左侧把不同二聚体的距离和结合能归一化,极小值均位于 r/r0=1,两侧曲线仍有不同形状。相同的归一化势阱深度,并没有消除不同成键作用随距离衰减的差异。
多原子体系具有大量构型自由度;移动一个原子时,其余原子是否固定、是否允许弛豫,会得到不同的能量曲线。论文里的 U(r) 或二维等高线往往只是高维势能面的一条路径或一个截面。路径上的最低点只针对指定的约束,沿未扫描的位移方向仍可能存在降低能量的构型。
势能 U 不含原子核运动的动能。在组成固定的体系中,给所有构型加上同一个能量常数,不改变受力或能量差;比较不同组成时则要保持元素参考能一致。力场通常还包括原子类型、参数和作用项的组织方式,势函数是其中的能量模型。电子结构计算中的赝势描述电子与离子实的等效作用,所作用的对象与 U(R) 不同。


势阱的深度、斜率和曲率能控制什么?
原子间保守力由Fi=−∇RiU给出。对双原子的径向运动,Fr=−dU/dr;取增大间距为正方向,势阱左侧常有正的排斥力,右侧有负的吸引力。能量用 eV、距离用 Å 时,力的单位为 eV·Å−1。
图2中 H2 的几种电子结构处理给出不同的势阱深度和平衡键长,归一化后曲线却较接近。拟合能重现给定参考曲线,也会保留参考方法的偏差。实线的 Morse 拟合接近对应数据,三次多项式在较大拉伸下明显偏离;势阱附近拟合良好,仅约束了平衡键长附近的能量变化。

常见 Lennard-Jones 势写为 U(r)=4ε[(σ/r)12−(σ/r)6]。其中 ε 是以无限远为零点时的势阱深度,σ 是势能过零的距离;平衡距离为 21/6σ。r−6 项描述色散吸引,r−12 项近似短程排斥。它适用于氩等近似非键相互作用模型,不能用同一套参数完整描述铜的金属键或碳的共价网络。
极小值处 dU/dr=0,附近曲率 κ=d2U/dr2 决定小位移的恢复力。在双原子的谐振近似下,角频率 ω=√(κ/μ),μ 为约化质量。势阱深度控制解离能尺度,曲率控制小振动频率;图1右侧几种势在最低点具有相同曲率,远离平衡位置后的吸引和排斥仍不同。晶体声子则由所有原子位移的二阶导数组成力常数矩阵。
实际模型常在距离 rc 以外省略短程作用。直接截断会使势能跳变,减去常数只能消除能量跳变;若截断处斜率非零,力仍有突变。图3左侧比较了势能、力及曲率的平滑处理,右侧压力曲线在邻居壳层进出截断范围时出现不同响应。

压力来自能量对体积的导数,弹性常数又涉及更高阶导数。图3中仅平移势能的蓝色曲线出现压力跃变,力平滑后压力连续,压力对体积的导数仍可能不连续。平滑处理也改变了截断范围以内的函数形状,会移动势阱并影响状态方程。改变截断距离或过渡函数后,原有参数拟合的平衡体积和弹性常数可能随之偏移。


为什么材料势函数要描述邻近原子?
双体模型把总能写为 U=Σiφ(rij),每一项只依赖一对原子的距离。同种原子、相同距离在双体项中总有相同的能量贡献。若第三个原子改变了原来的成键电子分布,单个 φ(r) 缺少这种环境变量;增加可调参数也无法消除这一函数形式的限制。
铜和碳的 DFT-LDA 结果中,配位数改变时,平衡近邻距离和平均每键结合能随之变化。图4包括二聚体、链、二维排列及多种晶体,其中含有假想结构;铜与碳的每键结合能随配位增加总体减小,而每原子内聚能还计入邻居数量。每键能是按配位数归一化的量,并非直接测出的独立化学键能。
嵌入原子法 EAM 常用于金属,总能包含双体项与嵌入能 ΣiF(ρi),其中 ρi=Σj≠if(rij)。F 对邻居贡献的非线性响应产生多体作用,使同一对原子的有效作用随周围配位变化。这里的 ρi 是模型中的局域密度变量,通常由参数化函数构造,不能据其数值绘制 DFT 电荷密度图。

共价材料对键的方向有强烈响应。Stillinger-Weber 势含显式三体键角项,Tersoff 势通过邻居环境改变键级;固定拓扑的分子力场则预先指定键、键角和二面角。研究断键或成键时,能量表达式必须允许键连接改变。ReaxFF 等反应力场采用连续键级与电荷模型,其反应能垒仍须用相应化学反应的数据检验。
机器学习势把元素种类、邻居距离及几何信息映射为原子能量,再求和得到总能。“机器学习”描述拟合方式,不规定只含双体作用;模型可以包含高阶多体信息。对无外场的孤立体系,总能应在整体平移、旋转和同元素原子置换后保持不变,力则随整体旋转改变方向。

图5的局域检验把中心原子及内区全部原子的坐标固定,仅改变外区原子的位置;若参考受力发生变化,局域输入相同的模型便无法区分这些力。原子能截断半径为 Rc 时,求导后的力依赖范围可达 2Rc,因为邻近原子的能量也贡献中心受力,局域检验所固定的内区应包含这一范围。离子材料的长程静电和界面的电荷响应,常需显式静电或全局电荷处理。


怎样确认势函数能描述的材料性质?
经验势的参数可以来自实验晶格常数、弹性常数或量子化学能量;机器学习势通常拟合大量构型的能量、力,有时还包括应力。训练数据限定了模型受到约束的构型范围。只包含平衡晶体小位移的数据,对液体中的短距离碰撞、表面配位不饱和或缺陷迁移鞍点的约束往往不足。
硅的几套 ACE 势采用不同数量的原子环境进行拟合,图6顶部的能量—体积曲线较接近,声子分支和热膨胀却表现出更大的模型差异。图例中的 3%、4%、5% 指原数据库中参与拟合的原子环境比例,第二幅直方图表示其体积分布;这些百分比没有表示计算误差。

硅的热容曲线几乎重合时,Grüneisen 参数和低温热膨胀仍能分辨几套势的差异。后两者对振动频率随体积的变化敏感。能量拟合、受力拟合和性质验证约束不同的导数或构型:声子检验平衡结构附近的力常数,缺陷形成能检验低配位能量,扩散能垒检验迁移路径中鞍点附近的能量。
正则系综中,构型的概率密度与 exp[−U/(kBT)] 成正比,kB 是玻尔兹曼常数。局部出现过深的人工势阱,会指数式提高该区域的统计权重,即使多数测试构型的平均误差很小,轨迹仍可能反复进入错误结构。对这些异常构型计算参考能量,可以检验新增势阱是否来自拟合误差。

图7的局部高估使对应构型的概率降低;局部低估则产生尖锐的概率峰,并减少其他区域的占据。这种能量误差改变的是采样分布本身。用于高温熔化、表面重构或化学反应的势,应在相应构型中核验能量与受力,测试样本还应与训练轨迹分开,避免相邻帧几乎相同而掩盖外推误差。
常规基态势把电子自由度消去,通常不随每个构型输出能带、磁激发或电荷转移量。强激光激发、电子温度变化或自旋态切换可能要求额外电子态变量或多势能面描述。粗粒化势若由有限温度的构型分布反演,常近似平均力势,已包含被消去自由度的熵贡献;其温度、密度或组成改变后,原有势函数可能失效。
