进了 Octopus 这个实空间 DFT 的坑之后我一直想找一套完整的“能跑通 明白每个参数在干嘛”的流程而不是对着官方手册复制粘贴。上一篇笔记写到了安装和入门今天直接拿硅(Si)开刀把基态计算完整过一遍。这篇的内容不只是给你一份可以运行的输入文件更重要的是把输入文件里每个关键参数背后的物理意义讲清楚包括为什么硅使用金刚石结构、为什么格子要这样写、k点网格和实空间网格间距分别在控制什么精度、SCF 收敛参数到底怎么调。如果你也已经装好 Octopus 但还在纠结第一步怎么跑或者以前只用过平面波程序QE、VASP想转到实空间方法这篇笔记应该能帮你少走很多弯路。1. 为什么选硅以及基态计算到底在算什么1.1 硅的晶体结构一个由两个原子撑起的“空洞”体系很多人刚开始算硅的时候习惯直接搜一个 cif 文件丢进去但 Octopus 不是结构文件导向的程序它更希望你用晶格矢量加原子坐标把体系描述清楚。所以你得先对硅的晶体结构有个底。硅在常温下是金刚石结构空间群是 Fd-3m。这个结构可以理解成“面心立方FCC格子 一个由两个硅原子组成的基元”每个硅原子周围有四个最近邻形成四面体配位。用晶体学术语说惯用晶胞里一共有 8 个原子但惯用晶胞不是最小单元。Octopus 做周期性基态计算时建议用它的初基胞描述这样只需要两个原子Si(0) : 分数坐标 (0.00, 0.00, 0.00) Si(1) : 分数坐标 (0.25, 0.25, 0.25)再加一套 FCC 晶格基矢。实验上硅的晶格常数在室温附近大约是 5.43 Å换算成原子单位Bohr就是 5.43 / 0.529177 ≈ 10.26 Bohr。这个数字先记着后面写输入文件有三种不同单位制写法我建议直接用Units eV_Angstrom把长度用 Å 写能量用 eV 写这样最符合我们做材料计算的人的习惯。我对硅这个体系情有独钟是因为它特别适合当“第一个 Octopus 周期性计算”它是半导体所以基态 SCF 不需要加 smearing也不会有费米面附近的收敛困难它又是标准的四配位共价键体系可以很好地检验赝势和网格精度更关键的是硅的光学性质、能带结构、激子效应全都是 TDDFT 里的经典算例这次把基态练扎实了后面可以直接接激发态计算。1.2 基态计算是后面一切计算的“地基”在 TDDFT 的框架里Octopus 的很多计算模式都要求你先有一个收敛的基态波函数和电子密度作为起点。比如接下来我要跑的td含时演化、unocc非占据态、lr线性响应都需要从基态的restart目录读取波函数和密度信息。换句话说基态自洽计算做得越稳后面所有模拟就越不容易出幺蛾子。可以打一个不严谨但很贴切的比方基态计算就是盖房子时打地基。如果你只是在地基没打好的情况下直接用kubo或者spectrum跑光学谱结果往往是密度矩阵有虚部分量、含时演化总能量飘掉、吸收谱出现一堆莫名其妙的负峰。我在跑 Octopus 之前也踩过这种坑——基态密度没有完全自洽就强行开始 TDDFT最后能量一直不守恒。这次我们就把地基打扎实。先算总能量再看电子本征值最后确认电子密度收敛到指定的阈值这样后面无论接能带计算还是光学响应心里都有底。1.3 Octopus 和平面波程序的本质区别这节不聊性能纯粹讲思路。很多从 VASP、Quantum ESPRESSO 转过来的人第一次看到 Octopus 的输入文件会懵它没有ENCUT截断能没有KPOINTS那种 l 坐标系一列列写也没有POSCAR。因为 Octopus 用的是实空间网格 有限差分方法而不是平面波基组。这意味着两件事第一波函数是直接定义在三维实空间网格点上的动能算符用高阶有限差分近似表示而不是在倒空间里对角化。它的好处是处理局域势、外场、含时演化非常自然做 TDDFT 时不需要处理平面波基组遇到的“对偶空间变换”和含时势的局域性问题。第二收敛测试的“旋钮”完全不同。平面波程序中你调ENCUT即平面波动能截断值而在 Octopus 里你需要调的是实空间网格间距Spacing单位通常为 Bohr以及模拟盒子的尺寸。对于周期性体系Spacing决定实空间精度k点网格决定布里渊区取样密度两者各自收敛后才算真正收敛。从客观角度说实空间网格方法对内存的消耗一般比平面波方法更大因为你需要在整个实空间盒子里存储波函数实部和虚部。但好处是对 MPI 并行更友好而且对于纳米结构、非周期体系、含外场问题省去了许多超胞技巧。所以 Octopus 不是要替代 VASP而是它更适合做“强场、超快、复杂环境”这一类问题。2. 输入文件拆解核心参数逐个讲2.1 从 CalculationMode 开始gs 到底是什么模式Octopus 的输入文件叫inp所有计算都通过开头的CalculationMode告诉程序你要干什么。基态计算就是CalculationMode gsgs是 Ground State基态的缩写。程序会做 Kohn-Sham 方程的自洽场SCF求解直到电子密度和总能量满足你设定的收敛判据。除了gs之外Octopus 中你还可能用到计算模式用途gs基态自洽计算获得电子密度、波函数、总能量unocc在固定密度下计算附加的未占据态常用于能带和光学矩阵元td含时演化比如飞秒激光激发下的电子动力学lr线性响应计算吸收谱、激发能em电动力学/光子学相关的本征模式计算kuboKubo 线性响应用于电导率等对于这篇笔记我只聚焦gs。但你要知道基态计算结束后 Octopus 会自动生成一个restart目录里面存了收敛后的波函数、密度和 Kohn-Sham 哈密顿量信息。下一次如果要跑unocc或td程序会优先从restart读取所以每个算例尽量单独建一个目录避免不同体系之间互相干扰。2.2 周期性设置如何告诉 Octopus 这是一个 3D 晶体Octopus 默认是孤立体系非周期所以做晶体计算必须显式打开周期性维度PeriodicDimensions 3这行告诉程序“我在处理三维周期体系”。接下来的关键是定义晶格矢量。对于面心立方硅通常有两种写法。第一种是直接给出三个晶格基矢单位依赖Units设置。FCC 的初基胞基矢可以写成矩阵形式a1 (0, a/2, a/2) a2 (a/2, 0, a/2) a3 (a/2, a/2, 0)如果 a 5.43 Å那么在Units eV_Angstrom下你可以这样写%LatticeVectors 0.000 | 2.715 | 2.715 2.715 | 0.000 | 2.715 2.715 | 2.715 | 0.000 %第二种更省事的写法是用LatticeParameters告诉程序晶格常数和角度让它从这些参数自动生成基矢。FCC 初基胞的三个基矢长度相等两两夹角都是 60 度于是可以写%LatticeParameters 5.43 | 5.43 | 5.43 | 60.0 | 60.0 | 60.0 %我推荐第二种写法因为它和结构描述的习惯一致。请注意这里的角度单位是度不是弧度。如果你使用的是默认原子单位Hartree Bohr那要写成 10.26 Bohr 和 60 度。为了少踩单位换算的坑直接在文件开头写Units eV_Angstrom后面所有坐标、晶格长度都按 Å 写能量相关量按 eV 写非常直观。我个人强烈建议新入门者用这种方式尤其是当你从 VASP 转过来时许多参考数据本身就是用 Å 和 eV 描述的直接对照更方便。定义完晶格之后就轮到原子坐标。对于周期性体系我建议直接用分数坐标这样和晶格常数解耦FractionalCoordinates yes %Coordinates Si | 0.00 | 0.00 | 0.00 Si | 0.25 | 0.25 | 0.25 %注意Coordinates块的排列是第一列是元素标签字符串带引号后面三列是坐标分量列与列之间用|分隔。块指令以%开头、%结尾。2.3 网格精度三件套Spacing、k点、盒子在平面波程序里你要担心的是截断能在 Octopus 里要担心的是实空间网格间距Spacing。Spacing指的是相邻网格点之间的距离单位是 Bohr即 a.u.1 Bohr 0.529 Å。这个值是 Octopus 最核心的收敛参数网格越密波函数描述越精确但内存和耗时呈立方增长。为什么网格间距会直接决定精度因为 Octopus 用有限差分表示动能算符如果网格太疏二阶导数的离散误差会非常大直接导致总能量偏低或电子密度振荡。考虑到硅体系比较温和建议先用Spacing 0.45跑通整个流程然后依次测试 0.40、0.35、0.30直到总能量基本不变变化小于 0.001 Hartree。关于这个收敛测试的细节我在第 3 节会展开。k 点网格控制布里渊区采样。这个变量在 Octopus 里写作%KPointsGrid 6 | 6 | 6 %它表示在倒空间三个方向上各取 6 个 Monkhorst-Pack 网格。硅是半导体带隙约 1.1 eV6×6×6 已经能给出非常合理的总能和能带。如果你后面想和实验带隙对比或做精确的光学计算再提高到 10×10×10 或 12×12×12。但基态阶段不建议一上来就开大 k 点网格因为 Octopus 对内存和计算量很敏感先小后大才是正道。盒子尺寸对于周期体系由晶格矢量决定通常不需要额外设置。但如果你后面想在晶体里加入一个局部缺陷、吸附分子或者外场探针就需要考虑扩大盒子或者用超胞。这篇笔记只约束在单胞基态所以盒子这一项基本是“自动”的。2.4 赝势与自洽收敛控制Octopus 支持多种赝势格式包括 UPF、PSML、假原子HGH等。我的建议是直接用 PseudoDojo 或者 SG15 系列这些都是经过系统测试的标准赝势可靠性高。具体流程是去 PseudoDojo 官网下载 Si 的标准 UPF 文件例如Si.UPF放到一个专门目录里然后在输入文件中指定PseudopotentialSet custom PseudoPotentialDir ./pseudos或者干脆在原子坐标里直接写%Coordinates Si | 0.00 | 0.00 | 0.00 | species_pseudo_dir Si.UPF ...但我更推荐先用PseudoPotentialDir统一管理。和 VASP 里的POTCAR一样赝势必须和元素一一对应不能张冠李戴。自洽场收敛控制是另一个重点。先看四个核心变量MaximumIter 200 ConvAbsDens 1e-7 ConvRelDens 1e-6 MixingScheme diis Mixing 0.3MaximumIter是 SCF 最大迭代步数200 步通常是够的。如果跑到 200 步还不收敛问题大概率不在参数上而在结构或者赝势。ConvAbsDens和ConvRelDens是收敛判据前者要求电子密度的最大绝对变化量小于某个阈值后者要求相对变化量满足条件。对硅这种半导体ConvAbsDens1e-7是合理的如果你只是临时测试流程可以先放宽到1e-5但最终计算一定要收紧。MixingScheme控制 SCF 迭代中密度或者势的混合方式diis是比较稳健的选择。如果你一开始发现 SCF 非常容易震荡可以把Mixing调小到 0.2 甚至 0.1。对于半导体体系不需要开Smearing默认关闭即可。如果算金属或有缺陷体系才要考虑KPointsSmearing之类的参数。3. 实操记录从结构建模到 SCF 收敛3.1 准备赝势文件在 Octopus 中赝势文件的摆放位置非常重要。我是这样组织的~/octopus-work/ ├── Si-gs/ │ ├── inp │ └── pseudos/ │ └── Si.UPFinp放在算例目录下pseudos子目录里放赝势。输入文件里写上PseudoPotentialDir ./pseudos跑之前看一眼octopus --version确认版本。Octopus 不同版本对 UPF 格式的兼容性有一些细微差别我用的是 12.x 及以上版本PseudoDojo 标准 UPF 基本没有兼容问题。如果你的版本比较老遇到解析赝势报错可以考虑换成 PSML 格式或者更新程序。3.2 一份能跑通的完整 inp 示例下面是我在硅基态计算中最开始用的输入文件结构不复杂但跑得非常稳CalculationMode gs Units eV_Angstrom PeriodicDimensions 3 ExperimentalFeatures yes %LatticeParameters 5.43 | 5.43 | 5.43 | 60.0 | 60.0 | 60.0 % FractionalCoordinates yes %Coordinates Si | 0.00 | 0.00 | 0.00 Si | 0.25 | 0.25 | 0.25 % Spacing 0.40 KPointsGrid 6 | 6 | 6 PseudopotentialSet custom PseudoPotentialDir ./pseudos MaximumIter 200 ConvAbsDens 1e-7 ConvRelDens 1e-6 MixingScheme diis Mixing 0.3关于ExperimentalFeatures yes这行我在 10.x 和 12.x 版本中见过不同表现。有些版本默认开启了一些高级特性有些则需要显式打开。实际测试时如果程序报“block structure needs ExperimentalFeatures yes”就加上这一行如果没报错不加也行。这是很正常的 Octopus 版本差异别被它吓住。3.3 运行 octopus 并观察 log进入Si-gs目录执行mpirun -np 4 octopus如果只是小体系测试单核跑也行。运行时注意看终端输出的 SCF 块通常长这样SCF CYCLE ITER # 1 : etot -15.78448211 H abs_den 0.1234E-01 ... SCF CYCLE ITER # 2 : etot -15.82103974 H abs_den 0.8732E-02 ... ...etot是总能Hartree 单位abs_den是密度变化。正常情况下你会看到etot快速下降然后慢慢趋于平稳最终abs_den收敛到目标阈值。程序的最终输出会写到static/info文件里这是最关键的几个结果文件之一。有一点值得提醒Octopus 的log文件是程序运行日志包含大量调试信息和时间戳static/info才是整理好的结果汇总。一开始跑完不要盯着终端看直接看static/info和log的末尾几行。3.4 收敛测试怎么做收敛测试不能只做一次更不能只看最终“收敛了没”。我的习惯是用三步走第一步固定Spacing 0.45扫 k 点网格4×4×4、6×6×6、8×8×8记录总能。如果 6×6×6 到 8×8×8 的总能差小于 1 meV/atom就说明 k 点网格够了。第二步固定 k 点为 6×6×6扫Spacing0.45、0.40、0.35、0.30。每一步都从头跑一遍 SCF记录总能和耗时。把总能和Spacing画成曲线你会发现它逐渐趋于平稳。对硅这种 sp 杂化体系一般Spacing 0.30 ~ 0.35左右可达到大约 meV/原子精度。如果算力和内存足够可以再跑到 0.25 做最终确认。第三步用测试出来的最终参数再跑一次完整计算记录static/info、输出文件路径和运行时长。这一步是为了将来复现时有个“基准”。我见过很多新手一上来就Spacing 0.18加12×12×12k 点直接把内存吃满然后抱怨 Octopus 跑不动。其实基态计算最重要的是流程跑通精度测试放在后面慢慢加一点都不迟。4. 看懂输出与检查结果4.1 static/info 里到底有什么计算成功结束后static/info里的信息很有用。我第一次跑硅的时候最关心的三个量分别是总能量Total -xxxx.xxxxxxxx Ha最高占据态HOMO和最低未占据态LUMO的本征能量Brillouin zone 积分后的电子态数N electrons。总能量本身在凝聚态物理里并不是一个可直接和实验对照的量因为 DFT 总能包含电子间的交换关联、离子间的经典相互作用等复杂项所以测试精度的主要手段是“看它收敛到了什么值”。但如果你想要和文献或者 VASP、QE 对照可以参考同一赝势族下的结果。PBE 泛函 标准赝势下硅单胞总能大多落在几十个 Hartree 的区间内各程序之间的差异主要来自赝势和截断方式。另一个值得关注的是本征值。static/info里会列出所有被占据的 Kohn-Sham 能级取决于 k 点可见最高占据态和最低空态间的带隙会明显低于实验值。比如 Si 实验带隙约 1.17 eV但 PBE 算出来通常在 0.5~0.7 eV 之间。这不是 Octopus 算错了而是 LDA/PBE 系统性低估带隙是 DFT 的常识。如果你要做带隙修正得靠unocc 混合泛函、GW 或者《TDDFT 学习笔记》里后续会写的含时方法。4.2 输出文件总览与后续用途一个成功的 Octopus 基态计算会在目录下产生一堆文件和子目录log运行日志保存所有中间迭代信息static/info最终结果汇总static/eigenvalues各 k 点的本征值表restart/gs基态波函数、密度、Kohn-Sham 势的二进制存储目录charge_density或density-restart电子密度的二进制或文本文件potential-restart有效势存储后续做含时演化时会用到。特别提醒restart/gs是后续所有计算的“命根子”。如果你准备做能带结构或者吸收谱不要轻易删除restart/gs。而且如果你换了赝势或者改了晶格参数最好重新建一个新目录把旧的 restart 留在原地别让程序读到旧数据造成混乱。4.3 和实验/文献数据对比基态计算最常做的对比是晶格常数的结构优化。你可以扫描多个LatticeParameters的a值比如 5.35、5.40、5.43、5.45、5.50 Å分别算基态总能找到总能最低点对应的晶格常数。通常 PBE 算出的硅晶格常数会比实验值偏大一点大致在 5.45~5.48 Å 范围。这个数量级偏差完全正常因为它和电子交换关联泛函的系统误差有关。如果你只想验证“我的 Octopus 基态流程没问题”最直接的参照是看绝对能量的自洽性同一输入文件连续跑两次结果应该完全一致。如果两次结果不一致说明重启读取出了问题或者并行环境下有随机性被引入需要排查。这个细节可能很多人没意识到但它真的很能测试程序状态。5. 常见问题与排查技巧实录5.1 SCF 振荡不收敛如果你看到 SCF 的能量在某个值附近来回震荡或者abs_den降到 1e-3 之后就卡住最常见原因是混合参数太大或者初始波函数质量太差。我的做法是把Mixing从默认值降到 0.2 或 0.1增加密度混合的稳定性把MaximumIter提高到 300 或 400给程序更多迭代步数如果依然不行检查是否在inp里错误设置了Smearing。硅是半导体别乱开。还有一个容易忽略的点如果你修改了Spacing或原子坐标后没有删除restart/gs里的旧波函数程序可能会从旧波函数重新开始导致初始波函数和新网格不匹配SCF 很容易跑飞。这种情况下删除restart/gs重新跑往往一次就收敛。5.2 k 点太少导致总能误差明显这个坑我踩过。有段时间我用 2×2×2 的 k 点算硅总能比 8×8×8 的结果高了差不多 0.05 Hartree。一开始我还以为是赝势选错了后来把 k 点加密到 8×8×8 才回归正常。这就是典型的布里渊区采样不足。排查方法很简单把 k 点网格翻倍比如从 4×4×4 到 8×8×8观察总能变化量。如果变化大于 1 meV/atom说明原来采样太稀。对于绝缘体/半导体k 点网格对总能的误差下降通常比较平滑如果你的体系是金属则需要配合Smearing使用否则 SCF 收敛会非常痛苦。这也是我建议新手拿硅当第一个算例的原因之一——不需要处理 smearing 的烦恼。5.3 Spacing 与内存爆掉实空间网格的一个现实问题是内存消耗大。假设你的盒子边长是 a 个 BohrSpacing 0.35那么每个方向大约有 a/0.35 个网格点总网格点数为三个方向的乘积。每个波函数通常需要存两个甚至更多复数数组每个复数双精度占 16 字节。算到几百个波段、几百万网格点的时候内存轻松上几 GB/核。经验法则跑 Octopus 的基态计算内存通常比 CPU 核数更先成为瓶颈。我自己的 32 GB 工作站算硅单胞开 8×8×8 k 点、Spacing 0.30就已经比较吃紧。建议先用Spacing 0.45跑通再逐步加密同时观察log里打印的“number of grid points”和“memory estimated”提前判断是否会爆内存。5.4 赝势报错Octopus 对 UPF 文件的兼容性整体不错但偶尔你会在运行初期看到类似“cannot read pseudopotential”或者“file format not recognized”的错误。这种情况一般有三个原因伪势路径写错检查PseudoPotentialDir是否指向了正确的目录UPF 文件本身损坏重新下载注意别用文本编辑器“智能”改掉换行符Octopus 版本太老对新版 UPF 格式支持不完整尝试换 PSML 格式或者换成 Octopus 官方测试过的赝势库。5.5 并行与重启技巧Octopus 是基于 MPI 的并行程序。我的经验是在 4~8 核规模下跑硅单胞基态比较合适核心数继续增加不一定带来线性加速反而可能因为通信开销导致总时间变长。另外它不像某些程序那样“随时断点续跑”但基于restart/gs实现了“自然再启动”只要你在同一目录再执行一次 octopus程序会读取restart/gs并继续迭代。所以如果中途因为超时被掐断直接重新运行即可不需要重头再来。需要警惕的是如果你改了输入文件里的关键参数比如晶格常数、k 点、赝势最好把旧的 restart 彻底删掉再重跑。我见过有人改了网格间距却没删 restart结果 Octopus 报了一个“网格尺寸不匹配”的诡异错误排查了很久才意识到是重启文件在捣乱。就我个人的使用体会来说Octopus 的基态计算流程一旦跑顺后续能带、态密度、含时演化这些计算就只是“换一个 CalculationMode 指定输出”的问题了。硅这个体系又是少有的“结果可靠、收敛友好、物理图像清楚”的试金石SCF 如果不收敛大概率是参数设置问题结果如果和实验对不上也能很自然地想到是泛函近似而不是程序 bug。希望你也能和我一样第一次在 Octopus 里完成硅基态计算后对这套“实空间网格 周期性边界 赝势框架”建立起手感。下一篇我打算直接接着写如何用unocc模式把硅的能带结构画出来顺带把能带带隙为什么偏低的问题再演示一遍。