进了 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 | 电动力学/光子学相关的本征模式计算 |
kubo | Kubo 线性响应,用于电导率等 |
对于这篇笔记,我只聚焦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 eV,6×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是收敛判据:前者要求电子密度的最大绝对变化量小于某个阈值,后者要求相对变化量满足条件。对硅这种半导体,ConvAbsDens=1e-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,扫Spacing:0.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 个 Bohr,Spacing = 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模式把硅的能带结构画出来,顺带把能带带隙为什么偏低的问题再演示一遍。