《npj Computational Materials》：Enhancing the efficiency of time-dependent density functional theory calculations of dynamic response properties
编辑推荐：
研究人员针对极端条件下X射线汤姆逊散射（XRTS）谱模拟时，含时密度泛函理论（TDDFT）计算量大、效率低的瓶颈，结合从动态结构因子（DSF）到虚时密度-密度关联函数（ITCF）的拉普拉斯变换映射，在虚时域进行展宽参数η的收敛性检验，并开发了一种基于约束的窄带波动衰减滤波方法。该方法在保证结果无明显偏差的前提下，将TDDFT计算速度提升了一个数量级，显著降低了极端条件物质诊断与性质模拟的计算成本，为高效、精准的从头算材料模拟提供了新工具。
随着高亮度、高分辨率、高重复率的X射线自由电子激光器（XFELs）如欧洲XFEL、LCLS和SACLA的出现，人类探测极端条件（如高压、强激光加热）下物质性质的能力得到了前所未有的提升。X射线汤姆逊散射（XRTS）已成为诊断此类极端条件物质状态的核心技术。为了从实验测量中提取准确的物理信息，如温度、密度、电子结构等，必须对实验测得的散射谱进行高精度、从头算的理论建模。其中，含时密度泛函理论（TDDFT）是计算电子动态结构因子（DSF）S(q, ω)等动态响应性质最精确的从头算方法之一。然而，与常温常压条件相比，极端条件下的强热激发、以及实际分析中需要考虑温度密度变化、探测器有限尺寸等因素，都使得TDDFT模拟的计算成本急剧增加，成为阻碍其广泛应用的主要瓶颈。计算中通常引入一个有限的展宽参数η来抑制数值波动，理论上η应趋于零以获得精确响应，但实践中η必须设为有限值以获得平滑的DSF曲线。减小η会放大由有限k点网格采样引起的伪波动（表现为曲线“振铃”），而为了在更小的η下获得平滑曲线，则需要更密的k点网格，这将导致计算成本以三次方甚至更快的速度增长。对于具有强热激发或无序原子结构的体系，此问题尤为严重。因此，发展能够在不显著增加计算量的前提下，高效、可靠地获得收敛、平滑的TDDFT动态响应结果的方法，对于支持前沿XRTS实验、深入理解极端条件物质性质、以及为辐射流体力学等应用领域提供可靠输入模型，都具有迫切的需求和重要的意义。本研究发表于《npj Computational Materials》期刊。
本研究主要采用了两个关键技术方法：1. 虚时域η收敛性检验：利用动态结构因子S(q, ω)与虚时密度-密度关联函数F(q, τ)之间通过双边拉普拉斯变换存在的一一映射关系，在虚时域对TDDFT计算中使用的展宽参数η进行系统性收敛测试。通过分析ITCF的最小值、面积等特征量对η的依赖关系，确定在给定计算参数下不引入明显偏差的最大η值，此η对应的DSF结果被定义为“参考DSF”。2. 约束性窄带波动衰减滤波：针对参考DSF中残存的、由有限k点网格引起的准周期性窄带数值波动，本研究开发了一种基于约束的Savitzky-Golay滤波方法。该方法结合修正的Akaike信息准则进行滤波器参数优化，并在频域施加惩罚以消除目标频率的波动贡献，同时确保滤波后的DSF所对应的ITCF与参考ITCF高度一致，从而最大程度地保留物理信息并抑制伪波动。
研究结果
方法一：虚时域η收敛性检验
研究人员首先阐述了所提出方法的理论基础，即S(q, ω)与F(q, τ)通过方程(1)关联。他们以固态密度氢（T=4.8 eV, ρ=0.08 g/cc, q=0.946 ?-1）为例，展示了不同η值（0.05 eV, 0.1 eV, 0.5 eV）下的DSF和ITCF（图1）。在频率域，直接判断DSF的η收敛性较为困难且模糊；而在虚时域，计算得到的ITCF能清晰地区分收敛与未收敛的结果。计算表明，当η ≤ 0.1 eV时，ITCF的特征量（如平移后的ITCF最小值、面积）基本收敛，而η=0.5 eV的结果则明显偏离。因此，η=0.1 eV被确定为该体系在当前计算参数下的参考值。进一步系统性的η扫描（图2）定量证实了ITCF在η ≤ 0.1 eV时达到收敛，从而为可靠地选择参考DSF提供了客观判据。
方法二：约束性窄带波动衰减
研究人员将参考DSF分解为物理部分Sp(q, ω)和非物理的窄带波动部分Sn(q, ω)。通过分析参考DSF导数的功率谱（图3a, b），识别出波动的主要频率成分。随后，他们采用结合了频域惩罚和二阶导数惩罚的修正Akaike信息准则，对Savitzky-Golay滤波器的多项式阶数和窗口长度进行网格搜索和自动优化（图3c），以获得能有效滤除窄带波动同时最大限度保持ITCF不变的滤波器参数。应用优化后的滤波器得到的不同参数组合的Sp(q, ω)高度一致（图3d）。滤波后残差的ITCF Fn(q, τ)的均值和均方根偏差相对于参考ITCF的最小值可以忽略不计（图4），满足方程(5)和(7)提出的约束条件，表明滤波过程引入的偏差极小。整个工作流程总结于图5。
应用实例一：固态密度氢
将上述方法应用于固态密度氢体系。结果显示，基于约束滤波得到的Sp(q, ω)与参考DSF（η=0.1 eV）的ITCF高度吻合，且与精确的路径积分蒙特卡洛（PIMC）基准结果符合得很好（ΔF?(q, τ) < 0.75%）（图6a-c）。相比之下，简单地增大η（0.3 eV, 0.5 eV）进行平滑虽然能得到光滑的DSF曲线，但其对应的ITCF明显偏离参考值及PIMC结果，表明引入了显著的物理偏差（图6a, b, e, f）。滤波产生的波动残差Sn(q, ω)围绕零值对称振荡，而增大η产生的“残差”则呈现出系统的、非对称的偏离，证明了滤波方法在去除伪信号的同时更好地保留了物理信息（图6d-f）。研究还通过对比不同k点网格（10×10×10 vs. 20×20×20）的计算结果，证明准周期性波动直接源于有限的k点网格。在η=0.2 eV时，对10×10×10网格的结果进行约束滤波后，与20×20×20网格的直接计算结果高度一致（图7a）。在η=0.1 eV时，即使使用20×20×20网格，结果仍存在小幅度波动，而约束滤波能有效消除这些波动，并使两种网格的结果趋于一致，表明滤波后的结果在η=0.1 eV时对k点密度基本收敛（图7b）。利用涨落-耗散定理和Kramers-Kronig关系，可以从Sp(q, ω)进一步推导出其他动态介电性质，如动态介电函数的虚部、动态电导率等。图8显示，基于Sp(q, ω)的结果与基于参考DSF（η=0.1 eV）但未滤波的波动结果，以及基于增大η（0.3 eV, 0.5 eV）的有偏平滑结果，在峰值等关键特征上存在明显差异。
应用实例二：等容加热铝
第二个应用实例是面心立方结构的铝。研究比较了常温（T=0.025 eV）和电子等容加热（T=25 eV）两种情况。虚时域η收敛性测试表明，在T=25 eV时，η ≤ 0.5 eV达到收敛（图9b, c）。在T=25 eV时，由于大量部分占据的能带和需要高截断能来描述L壳层（2s, 2p）热激发，计算面临巨大挑战。使用10×10×10 k点网格和η=0.5 eV得到的参考DSF在等离子体激元和L壳层特征区存在显著波动。将η增大至1.5 eV可平滑波动，但严重压制了关键的2s-2p束缚-束缚跃迁特征。而约束滤波在有效衰减窄带波动的同时，完好地保留了这一特征（图9a）。对于大动量转移（q=7.75 ?-1，对应于背散射几何）的情况，研究展示了更大的计算加速效益。对T=25 eV的铝，使用10×10×10 k点网格并应用约束滤波得到的结果，与使用32×32×32 k点网格（未滤波）的结果质量相当，均能清晰显示尖锐的2s-2p特征以及与L壳层特征重叠的宽康普顿信号，而计算时间降低了超过50倍（从约8.8×105CPU小时降至约1.6×104CPU小时）（图10）。
研究结论与讨论
本研究成功开发了一种普适性强的高效TDDFT动态响应性质计算新策略。该策略的核心在于利用动态结构因子与虚时密度-密度关联函数之间严格的数学映射关系，将评估展宽参数η收敛性的难题从嘈杂的频率域转移至平滑的虚时域，从而能够客观、可靠地确定“参考DSF”。在此基础上，针对有限k点网格引入的、具有特定频率的窄带数值波动，发展了一种基于物理约束（即要求滤波对ITCF的影响极小）的滤波衰减方法。两者结合，使得研究人员能够在无需持续加密k点网格（即不显著增加计算成本）的情况下，获得平滑、可靠且无显著偏差的动态响应性质计算结果。
该方法的价值和意义体现在多个层面。首先，它直接解决了极端条件下XRTS光谱模拟计算量巨大的瓶颈问题，通过消除对过度密集k点网格的依赖，将计算效率提升了一个数量级，使得对复杂体系（如强热激发、无序材料）进行系统性的TDDFT研究变得更为可行。其次，该方法为“真正”的从头算模拟提供了更严谨的框架。它用客观的、基于ITCF收敛的判据，取代了传统上依靠视觉主观判断选择平滑参数的做法，提升了计算结果的可靠性和可重复性。最后，本工作所倡导的“虚时域思考”范式，不仅适用于DSF和XRTS建模，还可以便捷地扩展到其他动态介电性质（如密度响应函数、介电函数、电导率、吸收系数等）的计算中。因此，该工作有望成为在广阔条件和研究领域内，估算材料性质、为实验测量提供理论支持的新标准工具，并对惯性约束聚变能源、高能量密度物理、材料科学等前沿领域产生积极的推动作用。
