★ FEATURED ARTICLE 简介这份资源围绕第一性原理计算软件Octopus展开聚焦时间依赖密度泛函理论TDDFT模拟计算面向具备一定量子化学与凝聚态物理基础、希望上手电子结构模拟的科研人员与研究生。内容涵盖基态DFT预计算、TDDFT参数设置、激发态响应计算、吸收与荧光光谱分析以及电子密度与轨道可视化等完整流程可帮助读者理解分子与材料光学性质的模拟思路。压缩包为octopus-4.1.0版本整体约4.58MB通常包含源代码、编译脚本、文档与示例输入文件等便于自行编译配置计算环境。目前已有409人学习下载适合作为TDDFT入门与实操的参考材料配合官方文档逐步掌握参数设置与输出解读为太阳能电池材料、光催化剂等方向的研究提供计算支持。1. octopus 与 TDDFT 模拟计算从一束激光到电子云的真实响应很多做第一性原理的同行第一次接触含时密度泛函理论TDDFT时都会卡在同一个地方基态算得挺顺一到“加光场、看激发、追电子动力学”就翻车。要么是体系太大跑不动要么是结果对不上实验吸收谱要么是根本不知道时间步长该取多少。Octopus 这个开源实空间 TDDFT 代码恰好是解决这类问题的常用工具之一——它不依赖平面波基组而是在实空间网格上直接传播含时 Kohn-Sham 方程天然适合做强场电离、高次谐波、光吸收谱和实时电子动力学。这篇笔记面向的是已经会跑基态 DFT、想把手伸向含时模拟的从业者也适合做光学响应、分子激发态、纳米结构等方向的人。我会按“概念立住 → 最小可跑 → 参数怎么调 → 坑在哪 → 进阶技巧”的顺序把 Octopus 做 TDDFT 的落地路径拆开讲每一步都尽量给到能直接抄的命令和输入片段。2. Octopus 做 TDDFT 的底层逻辑与最小可跑案例2.1 为什么是实空间网格而不是平面波基组TDDFT 的核心是求解含时 Kohn-Sham 方程i ∂ψ_j(r,t)/∂t [ -∇²/2 V_KS(r,t) ] ψ_j(r,t)平面波基组在处理周期性体系时很优雅但做含时传播时每个时间步都要做大量 FFT而且对真空层和孤立分子的描述需要很大的超胞。Octopus 的选择是实空间离散化把波函数直接存在三维网格上用有限差分近似动能算符。这样做的好处是真空区域不需要额外基函数边界条件用吸收边界或零边界直接处理时间传播可以用显式算法每个时间步只做局域乘法强场下的电离、电子云外溢实空间网格天然能描述。代价是网格步长必须足够小通常要保证每个方向至少 0.2–0.3 Bohr 的分辨率否则动能算符误差会直接毁掉激发能。常见做法是先用收敛测试确定网格间距再固定下来做含时计算。2.2 最小可跑案例一氧化碳分子的光吸收谱下面这个例子是我一般用来验证环境是否配好的最小体系CO 分子用 LDA 交换关联先做基态再做含时传播最后通过极化率得到吸收谱。输入文件分成两块基态 inp 和含时 inp。基态输入inp# 基态计算CO 分子LDA实空间网格 CalculationMode gs UnitsOutput eV_angstrom %Coordinates C | 0.0 | 0.0 | 0.0 O | 0.0 | 0.0 | 1.13 % BoxShape sphere Radius 6.0 Spacing 0.25 %Species C | species_pseudo | set | standard | lda | 1.0 O | species_pseudo | set | standard | lda | 1.0 % TheoryLevel dft XCFunctional lda_x lda_c_pz ConvRelDens 1e-6 MaximumIter 200 %Output wfs density %含时输入inp覆盖基态文件中的对应块# 含时计算施加弱电场计算极化率 CalculationMode td UnitsOutput eV_angstrom %Coordinates C | 0.0 | 0.0 | 0.0 O | 0.0 | 0.0 | 1.13 % BoxShape sphere Radius 6.0 Spacing 0.25 %Species C | species_pseudo | set | standard | lda | 1.0 O | species_pseudo | set | standard | lda | 1.0 % TheoryLevel dft XCFunctional lda_x lda_c_pz TDPropagator aetrs TDTimeStep 0.002 TDMaxSteps 5000 %TDExternalFields electric_field | 0.01 | 0.0 | 0.0 | envelope_constant | 0.0 | 0.0 % %TDFrequencies 0.0 | 0.0 | 0.0 % %TDOutput multipoles %运行命令# 先跑基态 octopus gs.log # 再跑含时注意会读取基态波函数 octopus td.log逻辑说明基态计算生成restart/目录下的波函数和密度含时计算默认从restart/读取初态。TDExternalFields里给了一个沿 x 方向的常数电场幅度 0.01 a.u.属于弱场线性响应区间。TDPropagator aetrs是近似 enforced time-reversal symmetry 传播子对大多数体系稳定。TDTimeStep 0.002对应约 0.05 as对价电子动力学足够。TDMaxSteps 5000总时长约 10 fs频率分辨率约 0.4 eV对宽谱吸收峰够用但想看精细振动结构要加长。参数说明Spacing是最关键参数0.25 Å 对 C、O 的赝势够用但如果用全电子或含 3d 元素要降到 0.15–0.18。Radius要保证波函数在边界衰减到可忽略否则吸收谱会出现虚假振荡。TDTimeStep不能太大一般取 0.001–0.005太大直接发散。2.3 从极化率到吸收谱后处理怎么做含时跑完后td.log里会有multipoles输出记录每个时间步的偶极矩。吸收谱正比于频率乘以极化率虚部。Octopus 自带oct-propagation_spectrum工具# 提取吸收谱 oct-propagation_spectrum它会读取td.general/multipoles输出cross_section文件。用 gnuplot 或 Python 画出来import numpy as np import matplotlib.pyplot as plt data np.loadtxt(cross_section) plt.plot(data[:,0], data[:,1]) plt.xlabel(Energy (eV)) plt.ylabel(Cross section (arb. units)) plt.xlim(0, 30) plt.show()逻辑说明cross_section第一列是能量第二列是截面。CO 的典型吸收峰在 8–14 eV 之间如果算出来峰位偏移超过 0.5 eV优先检查网格间距和交换关联泛函。LDA 通常会低估激发能这是泛函本身的问题不是代码 bug。3. 参数怎么设网格、时间步、传播子与交换关联的取舍3.1 网格间距与盒子半径的收敛测试实空间方法最怕的就是“网格不够细”。我一般会做一组收敛测试固定盒子半径逐步减小 Spacing看总能量和第一激发能的变化。下面是一个批量脚本的骨架#!/bin/bash # 收敛测试不同 Spacing 下的基态总能量 for sp in 0.30 0.25 0.20 0.18 0.15; do mkdir -p sp_$sp cp inp sp_$sp/ sed -i s/Spacing .*/Spacing $sp/ sp_$sp/inp cd sp_$sp octopus gs.log grep Total energy gs.log | tail -1 cd .. done逻辑说明每次改 Spacing 后重新跑基态记录总能量。当相邻两次总能量差小于 1 meV/atom 时认为收敛。对含时计算还要看第一激发峰位置是否稳定。经验值C、H、O、N 用 0.20–0.25 Å过渡金属用 0.15–0.18 Å含 f 电子用 0.12–0.15 Å。盒子半径的测试类似但更简单从 4 Å 开始每次加 1 Å看总能量和偶极矩是否变化。对弱场吸收谱半径不够会导致低能端出现虚假峰因为波函数被边界反射。3.2 时间步长与传播子的选择Octopus 支持多种传播子常用的有传播子关键词特点适用场景近似时间反演对称aetrs稳定默认大多数 TDDFT四阶 Runge-Kuttaetrs精度高耗时强场、短脉冲Crank-Nicolsoncn无条件稳定大时间步指数传播子exp精度高小体系时间步长的选择有个经验公式TDTimeStep ≤ 0.1 / E_max其中 E_max 是体系最高占据态能量的两倍。对价电子体系0.002–0.005 通常安全。如果跑着跑着能量不守恒先查时间步长再查传播子。# 检查能量守恒td.log 中 Total energy 应基本恒定 grep Total energy td.log | head -20如果能量漂移超过 1e-4 a.u.说明时间步长太大或传播子不合适。强场下建议用 etrs虽然慢一点但不会翻车。3.3 交换关联泛函与赝势的匹配TDDFT 的激发能对泛函很敏感。LDA 通常低估 0.3–0.8 eVGGA如 PBE略好但仍有系统误差。如果做吸收谱常见做法是用 LDA 或 PBE 先跑再用剪刀算符平移。如果做强场电离泛函影响相对小因为电离主要由外场决定。赝势方面Octopus 自带standard和hgh等系列。注意赝势的价电子数要和 XCFunctional 匹配。比如用lda赝势就配lda_x lda_c_pz不要混用。混用会导致总能量异常含时传播直接发散。# 检查赝势文件是否被正确读取 grep Species gs.log grep valence gs.log4. 避坑与排查TDDFT 跑不起来时先看这五条4.1 含时计算第一步就发散能量变成 NaN现象td.log里 Total energy 在几步内变成 NaN 或极大值。原因最常见的是时间步长太大或者网格间距太粗导致动能算符误差爆炸。其次可能是赝势与泛函不匹配或者初始波函数没读进去。解决先把 TDTimeStep 降到 0.001Spacing 降到 0.18重跑。如果还发散检查restart/目录是否存在基态是否收敛。基态没收敛就做含时等于在错误初态上传播必炸。4.2 吸收谱出现负值或剧烈振荡现象cross_section在低能端出现负截面或者整条谱线毛刺很多。原因总传播时间太短频率分辨率不够或者盒子半径太小波函数被边界反射或者电场幅度太大进入非线性区。解决把 TDMaxSteps 加到 10000 以上Radius 加 2 Å电场幅度降到 0.001。如果还振荡检查multipoles输出是否等间隔时间步长是否稳定。4.3 激发能比实验值低 1 eV 以上现象主吸收峰位置明显偏低。原因LDA/PBE 的固有误差不是代码问题。TDDFT 在绝热近似下对电荷转移激发和双激发描述很差。解决换用杂化泛函如 B3LYP或做剪刀修正。Octopus 支持XCFunctional b3lyp但计算量会大很多。如果只是看趋势LDA 加 0.5 eV 平移也能用。4.4 并行效率低跑几天没进展现象MPI 并行后加速比远低于线性。原因实空间网格的并行域分解在体系较小时通信开销大或者 FFT 在含时传播中成为瓶颈。解决小体系50 原子用 OpenMP 线程并行更划算大体系用 MPI但每个域至少要有 10^4 个网格点。检查octopus编译时是否链接了并行 FFTW。4.5 重启含时计算后结果对不上现象中断后从restart/继续跑谱线形状变了。原因含时传播的相位信息存在波函数里重启时如果时间步编号或外场相位没对齐会引入跳变。解决尽量一次跑完。如果必须重启用TDMaxSteps分段但每段都要从同一初态重新激发不要接着上一段的末态跑。或者用FromScratch no并确保td.general/multipoles是追加模式。5. 进阶技巧用 delta 脉冲和频谱分析快速定位激发态5.1 delta 脉冲一次传播拿到全谱弱场吸收谱的常规做法是加一个常数电场但更高效的是加 delta 脉冲在 t0 时刻给一个瞬时电场然后自由传播。这样一次计算就能得到整个频率范围的响应。%TDExternalFields electric_field | 0.001 | 0.0 | 0.0 | envelope_delta | 0.0 | 0.0 %逻辑说明envelope_delta在第一个时间步施加电场之后为零。幅度要小保证线性响应。传播时间要足够长频率分辨率 ΔE 2π / T_max。比如 T_max 20 fsΔE ≈ 0.2 eV。5.2 用 Python 做峰值拟合与振子强度提取跑完 delta 脉冲后multipoles文件里是偶极矩随时间的变化。做傅里叶变换得到极化率再提取峰位和振子强度import numpy as np from scipy.signal import find_peaks # 读取偶极矩 data np.loadtxt(td.general/multipoles) t data[:,0] dip data[:,1] # x 方向偶极矩 # 去趋势加窗 dip dip - np.mean(dip) window np.hanning(len(dip)) dip_fft np.fft.rfft(dip * window) freq np.fft.rfftfreq(len(t), dt[1]-t[0]) # 转成能量取虚部 energy freq * 27.2114 # Hartree - eV strength -np.imag(dip_fft) # 找峰 peaks, _ find_peaks(strength, height0.1*np.max(strength)) for p in peaks: print(fPeak at {energy[p]:.2f} eV, strength {strength[p]:.3f})逻辑说明加 Hanning 窗减少截断效应rfft输出复数虚部对应吸收。峰位就是激发能峰高正比于振子强度。如果峰位和实验差 0.5 eV 以内说明泛函和网格都合理。5.3 一个我常犯的错误忘了检查盒子边界早期我做纳米团簇的 TDDFT吸收谱低能端总有一个莫名其妙的峰。查了两天才发现是盒子半径不够电子云在边界被反射形成了驻波。后来养成习惯每次含时计算前先看基态密度在边界处的值如果大于 1e-6就把 Radius 加 2 Å。这个检查用octopus输出的density文件就能做# 用 octopus 自带工具看密度边界值 oct-center-geom inp或者直接在gs.log里搜Density at boundary。这个习惯帮我省了很多后悔药。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
