简介这是一份面向量子化学、凝聚态物理及材料模拟研究者的 Octopus 4.1.0 资源包聚焦基于时间依赖密度泛函理论TDDFT的电子结构计算。软件支持从基态 DFT 预计算、TDDFT 参数设置、激发态响应计算到吸收、荧光、拉曼光谱分析与可视化呈现的完整流程尤其适合研究太阳能电池材料、光催化剂及复杂分子激发态性质的科研场景。压缩包以 gz 格式打包约 4.58MB便于在 Linux 类环境中部署包内涉及的源代码、编译脚本、文档和示例输入文件可支撑使用者完成从环境配置、输入文件编写到结果解读的闭环实践省去大量摸索时间。目前已有 409 人学习下载对刚开始接触 Octopus 或希望系统理解 TDDFT 模拟流程的初学者和研究人员而言是一份不错的入门参考亦可作为搭建专项计算流程的起点。1. 用 Octopus 跑 TDDFT到底能算出什么值不值得入坑红外光谱、紫外可见吸收、电子激发态寿命、非线性光学响应、表面等离激元共振——这些听起来分属不同领域的光学性质在 Octopus 里其实共用同一条 TDDFT 技术路线先让体系在基态自洽收敛再用含时外场扰动它最后从电子密度的演化中提取响应信息。这个流程的好处是你不需要为每一种光谱单独搭一套方法学一套实空间网格上的含时演化代码就能覆盖从气相小分子到纳米团簇、二维材料、周期性表面的光物理问题。反直觉的地方在于TDDFT 最耗时的部分往往不是传播过程本身而是你把模拟盒子建多大、网格切多细、基态收敛到什么程度——这三个前置决定占了踩坑总量的八成。这篇文章面向两种读者一是已经用 VASP、Quantum ESPRESSO 这类平面波代码做过基态 DFT刚接触含时动力学的新手二是被商业化软件的限制、力场方法精度不够逼过来的从业者想在光响应模拟上找一个开源、透明、能自由改参数的方案。我会按实际跑计算的顺序展开先讲 Octopus 的实空间离散和含时演化原理再用一个可以完整照抄的最小输入算吸收谱接着分析数据怎么拿振子强度最后单列一章讲我这些年翻过车的坑。全程用原子单位制不绕弯。2. 实空间网格与含时演化为什么 Octopus 用网格而不是平面波2.1 实空间格点离散Octopus 和平面波代码的本质区别Octopus 不用平面波基组而是把模拟盒子划分成均匀的实空间格点直接在格点上表示 Kohn-Sham 波函数和电子密度。这个技术路线的直接后果是没有 G 向量截断能这个收敛参数但多了一个网格间距参数 h 和盒子形状参数。平面波代码处理周期性体系天然顺手因为 Bloch 波和周期性边界是配套的而实空间网格对非周期体系孤立分子、纳米团簇、有界表面更直观——你直接在有限区域里切一个球、圆柱或盒子把波函数掉到边界附近之前截断掉就行。网格间距 h 决定了动能算符的离散精度。Octopus 里默认用高阶级数有限差分表示拉普拉斯算子h 越小精度越高但格点数量按三次方增长。我通常的收敛测试流程是固定盒子不动把 h 从 0.4 扫到 0.3 再到 0.25 原子单位对比基态总能的变化——总能差小于 1 meV/原子 就够用没必要无脑加细。这里有个和平面波习惯不同的地方不要只盯着总能收敛要看你关心的光学响应量收敛没有。偶极矩、跃迁能对格点误差更敏感用吸收谱的峰位漂移作为收敛判据更靠谱。2.2 含时 Kohn-Sham 方程与传播子TDDFT 的核心循环加外场后电子要满足含时 Kohn-Sham 方程i ∂/∂t ψ_i(r,t) [-1/2 ∇² v_ext(r,t) v_H(r,t) v_xc(r,t)] ψ_i(r,t)v_ext 是外加激光脉冲或 delta 脉冲的势v_H 是 Hartree 势电子间库仑排斥的平均场v_xc 是交换关联势。注意 v_H 和 v_xc 都依赖当前时刻的电子密度所以每个时间步都得先由波函数造出密度、再算这两个势、再推进波函数。这就是自洽传播。Octopus 里传播子的实现不是简单的一阶 Euler 推进默认是 Taylor 展开或近似指数传播子envelope 方法级数阶数由 PropagationType 下的 Taylor 阶数控制。时间步 dt 的选取有一个硬约束它必须小于体系最高占据态与最低空态间驱动振荡的倒数否则就会出现密度振荡发散、总能量漂移。我一般推荐 dt 取 1.02.0 原子时间单位1 a.u. ≈ 2.419e-17 秒和 25 eV 量级的电子激发能量比其振荡周期大概是 0.5 a.u.所以 1.5 a.u. 的步长在工程上是一个安全又能接受的折中。2.3 交换关联泛函与赝势的选择影响跃迁能的三个变量赝势和泛函的选择会直接决定算出来的激发能准不准。LDA 和 GGA 泛函的系统性缺陷是低估 HOMO-LUMO 间隙导致吸收谱整体红移杂交泛函有 HF 交换成分间隙明显改善但计算量暴涨——实空间网格对非局域交换算符的处理不算特别高效所以大体系一般先算 LDA 看趋势再挑小模型用杂交泛函验证极值点。Octopus 里常用的赝势是 SG15 和 Hamann 格式OGM 是从 PseudoDojo 转换来的。跑 TDDFT 之前先跑一遍基态检查 HOMO-LUMO 间隙和实验值/文献值的偏差心里有数再往下走。Octopus 的赝势文件存放是用户自己控制的在 inp 文件里通过 PseudopotentialSet 指定格式、Species 条目写元素名和对应的文件名。这里有个入门常见的暗坑赝势文件路径写错时Octopus 报错不一定很直白有时只提示找不到文件有时会出现读取了错误格式文件的乱码报错。建议把赝势集中放在一个目录在 PseudoPotentialDir 里用绝对路径或相对路径统一指定避免每次换计算目录都手改。3. 最小可复现的 TDDFT 吸收谱计算从盒子和基态开始跑3.1 盒子类型与收敛测试球、最小盒子与外部区域的取舍Octopus 支持的盒子类型有 sphere、parallelepiped、cylinder、box 等。对孤立分子我一般优先选 sphere 配合 SphericalBox 的半径参数——它给各个方向的真空一致、不会引入形状伪影。radius 怎么定先粗略按原子的共价半径之和乘以 1.2 再加大约 4 Å但更稳的做法是用 Octopus 自带的基态收敛信息看如果波函数在边界处没有掉到接近零说明盒子小了边界切断会造成电子密度畸变。具体的收敛逻辑是固定 h 只加盒子半径观察吸收谱第一个峰位的移动。比如半径从 8 加到 10 再加到 12 a.u.如果峰位移动小于 0.05 eV就可以认为收敛了。盒子太大也不是好事格点数三次方膨胀内存占用涨得很快。另一个绕开大盒子的手段是 Octopus 的盒子外部区域BoxShapeOutside它只在外围用吸收边界条件absorbing boundaries抑制波函数反射中心区域保持正常精度。这个功能在需要模拟电子解离、强场电离时是主场域普通吸收谱不太用。3.2 基态计算输入文件inp 文件怎么写才不是玄学Octopus 的输入文件是 Fortran namelist 风格和 VASP 的 KPOINTS/POTCAR 体系完全不同。最省心的办法是用官网的初始化工具生成模板再手工改关键参数。下面这份是跑乙烯模拟的基态输入可以做骨架抄写CalculationMode gs ExperimentalFeatures yes TheoryLevel dft BoxShape sphere Radius 12.0*angstrom Spacing 0.25*angstrom %Species C | species_pseudo | pseudo | 6 | hgh_lda H | species_pseudo | pseudo | 1 | hgh_lda % %Coordinates C | 0.0000 | 0.0000 | 0.6710 | units | angstrom C | 0.0000 | 0.0000 | -0.6710 | units | angstrom H | 0.0000 | 0.9230 | 1.2430 | units | angstrom H | 0.0000 | -0.9230 | 1.2430 | units | angstrom H | 0.0000 | 0.9230 | -1.2430 | units | angstrom H | 0.0000 | -0.9230 | -1.2430 | units | angstrom % KohnShamConverter yes Smearing 0.005 ExtraStates 4 ConvAbsEv 1.0e-8 Output potential density OutputHow xcrysden这段配置里几个关键点需要解释一下。Spacing 的值给了两个乘积项的组合写法写成 0.25*angstrom 是让 Octopus 自动换算成原子单位比手算 0.25/0.529177 少一步。Radius 同理Smearing 用 0.005原子单位约 0.14 eV是给金属或多参考体系提供占据数展宽的平滑对闭壳层分子可以设 0 或去掉。ExtraStates 给多少个空态直接影响 TDDFT 能描述多少激发态太少会切掉高能跃迁太多则白吃网格和传播时间。通常先设 68 个再根据激发能范围增删。ConvAbsEv 是总能收敛判据1e-8 原子单位约 2.7e-7 eV对多数光学响应计算足够紧了再小纯属浪费算力。跑完基态后要检查 oct-static-info 文件里的 Fermi 能级、总能和占据数异常占据分数偏离 0 或 1 超过 0.05通常是盒子或赝势出问题的信号。3.3 线性响应吸收谱delta-kick 与传播参数吸收谱的计算采用线性响应 TDDFT 的时域方案就是给体系一个极短的动量扰动再让密度自由演化。输入文件里的变化集中在三个位置CalculationMode 改成 td加一个 delta-kick 外场定义再把传播参数补上CalculationMode td FromScratch yes %TDOutput multipoles % %MomentumKick 0.001 | 0.0 | 0.0 % %TDKickFunction delta % TdPropagator etrs Dt 1.0 TdMaxSteps 1500MomentumKick 是沿 x 方向的动量转移数值一般取 0.001 原子单位。这个值不能太大否则进入非线性区太小则会因为数值噪声淹没有效信号。常见取值 1e-3 到 1e-2换算成电场强度大约对应 10⁷ V/m 量级放心落在线性区。传播到这里 TDDFT 步数乘 Dt1500 × 1 a.u. ≈ 3.6e-14 秒对应到光谱上的能量分辨率 ΔE 2π/总传播时间≈ 658/1500 a.u.≈ 0.44 eV 的粗略估计实际分辨率还要看短时傅里叶的手法。想分辨两个距离 0.2 eV 的峰就把 TdMaxSteps 提到 30005000代价是传播时间翻倍、内存和墙钟时间同步上涨。两条经验一是基态必须用和后续 TDDFT 完全一致的网格和赝势换任何一个参数都会引入基态-传播态之间的系统性错位吸收谱会出现基础峰位偏移二是如果基态用了 SmearingTDDFT 里也要维持同样的 Smearing不然占据数分布突变会反映在偶极矩上。4. 从偶极矩到吸收光谱数据处理与振子强度提取4.1 时域偶极矩信号multipoles 输出的正确打开方式Octopus 在 TDDFT 传播结束后会在 td.general 目录下生成多个文本文件名字类似 multipoles。这个文件里按列写时间、偶极矩的三个分量、以及更高阶多极矩。我们做吸收谱关心的只是偶极矩 d(t)它在 delta-kick 后会在零附近振荡振荡的振荡就是体系的激发模式。先别急着拿原始数据做傅里叶变换第一步是看时域信号有没有衰减——衰减是吸收边界存在时的预期行为但如果衰减过快说明盒子不够大、波函数被边界误吸收完全没有衰减并且出现拍频拍得很乱则可能是步长偏大或传播子阶数不够。提取吸收谱的频域做法是对 d(t) 做傅里叶变换实部对应吸收、虚部对应色散。单位换算上Octopus 输出时间是原子单位频率是原子频率单位将频率从 a.u. 转成 eV 要乘 27.2114转成 nm 要算 1239.84/E(eV)。这一步用脚本处理最稳手算容易错。需要注意的是Kick 沿 x 方向只激发 x 方向对称性允许的跃迁要拿完整的各向同性吸收谱例如溶液中的随机取向分子需要分别做 x/y/z 三个方向的 kick 再把谱平均或者用 Octopus 的 Kick 里指定各向同性选项。4.2 用 Python 从 multipoles 提取吸收曲线下面这段脚本是纯 Pythtb 处理 Octopus 输出标准流程的一部分依赖 numpy 和 matplotlib零额外依赖import numpy as np import matplotlib.pyplot as plt data np.loadtxt(td.general/multipoles) time data[:, 0] dipole_x data[:, 1] # 第一列是时间, 第二列开始是 x/y/z 偶极矩 # 去掉直流项被 kick 之后的平均位移 dipole_x dipole_x - np.mean(dipole_x) dt time[1] - time[0] n len(time) freq np.fft.rfftfreq(n, ddt) # 用 Hamming 窗抑制时域截断带来的谱泄漏 window np.hamming(n) spec np.fft.rfft(dipole_x * window) abs_spec np.abs(spec) # 换算成 eV 和相对强度 energy freq * 27.2114 intensity abs_spec / np.max(abs_spec) np.savetxt(absorption_spectrum.txt, np.c_[energy, intensity], headerenergy_eV intensity_arb) plt.plot(energy, intensity) plt.xlabel(Energy (eV)) plt.ylabel(Absorption (arb. units)) plt.xlim(0, 15) plt.savefig(absorption_spectrum.png, dpi300)这里的 FFT 处理的要点是直接对时域偶极矩做 FFT 会带出强烈的谱泄漏因为传播只做了有限步就截断了。加 Hamming 或 Blackman 窗可抑制泄漏但代价是峰形会变宽、弱峰会略微偏移本质上这是一个分辨率与泄漏的折中。我一般在确认峰位时加窗减泄漏在做振子强度积分时用没加窗的原始数据——因为窗函数会把峰面积改变。另外时间序列长度最好补零到 2 的幂做 FFT 加速但补零不增加真实分辨率只改插值密度。4.3 振子强度、光谱积分与实验对比的换算细节吸收谱的单位和实验对比时要小心理论谱输出的是任意单位强度要做绝对强度的振子强度计算需要从 TDDFT 的状态信息里取跃迁偶极矩。Octopus 的 casida 模式是一次完整的线性响应本征值求解等价于 Casida 方程会直接输出振子强度如果你跑的是时域方案则利用线性响应关系f 2 m_e ω |d(ω)|²在某个激发峰附近对 ω|d(ω)|² 积分再乘 2/3对随机取向体系做方向平均。实际操作时这个转换最容易踩坑的是单位d(ω) 是原子单位下的偶极矩谱密度积分的频率微元 dω 必须也是 a.u.最后乘的因子才正确。我见过不少人拿着 eV 刻度直接积分得到的振子强度偏大一个数量级。用输出的吸收谱和实验对比时还要注意一个物理细节理论上是零温、无展宽、无振动耦合的电子谱实验是室温、有振动精细结构和自然展宽的。所以不要指望峰形和实验一致看峰位和相对强度的趋势更有意义。峰位误差在 GGA 泛函下普遍红移 0.30.8 eV这个系统性偏差靠泛函升级PBE0、B3LYP、含 long-range corrected 的泛函改善但每个泛函引入的误差模式不同做文献对比时要标注泛函和赝势来源。5. TDDFT 模拟计算的五条血泪教训从盒子和时间步长到并行与重启5.1 盒子的坑轨道在边界被反射导致的“假吸收峰”现象吸收谱里出现了一个窄而高的尖峰能量位置不在任何文献报道的激发态附近而且换网格间距后这个峰的位置和强度乱跳。原因盒子太小高能轨道延伸到边界处没有掉到零被边界反射回来和真实波函数干涉产生虚假的激发模式。解决加大 Radius 或改用吸收边界验证方法是把 Radius 从 12 加到 16 a.u.再看这个峰是否稳定。如果峰位移动超过 0.2 eV就继续加直到不移动为止。5.2 时间步长的坑dt 太大导致密度振荡发散总能量漂移现象传播到几百步以后总能量开始持续上涨、偶极矩信号迅速出现锯齿状高频抖动。原因Dt 超过了体系所能承受的极限含时演化不稳定本质是 CFL 条件在实空间网格上的体现——电子振荡的特征周期进入了跨步长的数值失稳区间。解决把 Dt 从 2.0 降到 1.0 甚至 0.5 a.u.如果因为总传播时间太长导致 dt 降不下来考虑就是用更高阶传播子TdPropagator 改成 aetrs 或使用近似指数传播子公式它每步成本高一些但允许更大的步长。一个通用的检查指标是看总能量随时间的漂移量允许漂移在 1e-5 a.u. 内超过就降 Dt 或提阶数。5.3 单位换算的坑激光参数、能量条目、时间输出全混用现象算了一个激光驱动电子动力学输出的电子密度云在空间上看起来没问题但入射激光的频率和实验对不上也不知道该怎么换算。原因Octopus 全内部用原子单位但输入文件常写成 *angstrom、*eV 这类带单位的字面量。如果不写单位它就按 a.u. 解释——比如你把 3.1 eV 直接写成 3.1那就成了 3.1 Hartree ≈ 84 eV 的极端紫外光物理完全跑偏。解决在涉及能量、长度、时间的输入里养成每项都写清楚 units 的习惯例如 “ev”、“angstrom”、“fs”明确告诉 Octopus 执行换算输出解读时记住 1 Hartree 27.2114 eV1 a.u.时间 2.4189e-17 秒1 a.u.长度 0.52918 埃。5.4 并行与磁盘的坑小体系乱开核、大体系忘记清老文件现象一个只有十几个原子的小分子开了 64 核结果墙钟时间反而比 8 核更久另一个是算到一半报磁盘满Octopus 的分区格式中间文件把临时文件夹填爆了。原因Octopus 里的实空间网格按域分解并行每个核持有空间的一块子域如果体系格点总数太小域间通信开销超过了计算收益并行扩展不起来。另外 Octopus 的迭代过程要写大量临时文件包括 restart 用的 tmp 目录不清旧文件就会积累出几十 GB 的磁盘占用。解决小体系先测 1/2/4/8 核的加速比找到拐点就收手定期清理 restart/td.general 目录里不再用的中间结果或用 Octopus 自带的磁盘占用信息统计功能确认哪些文件可删。5.5 重启与续算的坑从中间状态接续 tddft 传播时偶极矩错乱现象用 restart 接着之前中断的传播往下跑偶极矩信号在重启点出现断崖式跳变后续演化信号明显偏掉。原因Octopus 的 TD 重启依赖所有态的信息完整写盘如果中途改了网格、盒子或传播参数重启读入的波函数和当前配置不一致密度演化的初态就错了。解决重启只用于完全一致的配置有任何参数改动就老老实实 FromScratch yes 重跑。另一个更隐蔽的情形是并行拓扑变了也能导致重启文件读不回来——重启时保持和中断时相同的并行分区数才最省事。6. 进阶技巧把 TDDFT 从线性响应推上强场与动力学TDDFT 的传播器天然支持非线性动力学delta-kick 只是拿它当黑匣子用拿强激光场直接驱动才是它的真本事。把 Kick 换成激光场定义时Octopus 提供 TDExternalFields 这个变量可以叠加多个激光脉冲每个脉冲能单独设置频率、强度、偏振和包络函数如 sin²、高斯、梯形上升下降沿。注意这套参数和外场相互作用项直接决定激发通道我习惯先用一个弱场跑线性区做基准再逐渐加场强寻找响应拐点——场强从 1e-4 到 1e-2 a.u. 之间能看到明显的非线性阈值这个阈值和电子能带结构的带隙有直接关联。另一个实用进阶是电子温度与激发态布居分析。做完传播后Octopus 可以把每个 Kohn-Sham 轨道投影到参考态上得到激发态的占据布居随时间的变化。实现上需要先把基态波函数存好作为参考投影基再在 TDDFT 传播过程里定时输出这些投影系数。这个功能是研究光致电荷转移、激子解离动力学的核心抓手你可以从布居变化的斜率读出电荷转移速率而不再停留在「看谱峰位置对不对」的表层。关于热电子动力学还有一个高频需求在金属纳米团簇上算光激发后的电子温度演化。这时需要把 Smearing 换成有限温度展开如 Fermi-DiracOctopus 里直接设置 SmearingFunction fermi-dirac温度单位设定成 Kelvin 即可。此时要注意的是传播算符的构建会依赖热平衡态密度如果直接用普通基态传播会看到温度相关的振动噪声。最后提一个我到现在都保持的习惯每套新体系的第一次 TDDFT 计算先给偶极矩信号做一次 10 步以内的短传播肉眼确认振荡平滑且能量守恒再放开步数跑长传播。这个习惯帮我至少挡掉了三成以上的无效长任务——步长、盒子、赝势的任何一项有问题短传播里都会有预兆。希望这套流程能帮你在 Octopus 的 TDDFT 模拟上少绕一些弯直接把算力花在该花的地方。本文还有配套的精品资源点击获取