1. 基态计算前的整体设计与思路拆解1.1 为什么要做Si的基态计算学Octopus的人一般都有个共同特征要么是搞TDDFT做激发态的要么是被复杂势场逼到墙角了。但TDDFT的起步永远是基态计算跑不掉的。我之前在笔记一里谈过Octopus怎么装、怎么调参数这篇直接从Si的基态计算开始把整个流程掰开揉碎讲一遍。选Si作为练习对象理由很充分金刚石结构、单元素、面心立方Bravais格子对称性高但又不是简单到没有物理内容——间接带隙半导体带隙值约1.1 eVDFT算出来的数值大家都知道会偏低这个偏差本身就能帮助你理解交换关联泛函的局限性。而且硅的赝势文件很容易获取算例多参考资料丰富任何一个环节出了问题都能比较快地找到参照。我用Octopus做过很多体系从单原子到小分子再到周期性固体坦白说Octopus在周期性计算上的上手成本比Quantum ESPRESSO高不少尤其是输入文件的自由度很大同样的体系一百个人能写出一百种输入文件。但你一旦驯服了它后面做real-time TDDFT、处理复杂外场、研究非线性光学响应会顺畅得多。所以这一篇不教你最短路径教你的是“能跑通且知道为什么这么写”的方法。1.2 计算方案选型能带、晶格常数与赝势做基态计算前先得把物理模型定下来这一步决定了后面所有的计算结果是否可信。首先是晶格常数。Si的实验晶格常数是5.431 Å空间群Fd3m227每个原胞两个原子分别位于(0,0,0)和(1/4,1/4,1/4)。严格来说做第一性原理计算应该先做结构优化把晶格常数和原子位置都弛豫一遍再取优化后的结构做性质计算。但在学习阶段直接用实验晶格常数跑基态计算问题不大误差可以接受。如果你想更严谨可以在Octopus里用vc-relax做变量胞弛豫或者手动扫描不同晶格常数下的总能拟合出平衡晶格常数。后者虽然笨但直观能让你对能量-体积曲线有非常直观的理解。其次是交换关联泛函的选择。Octopus支持LDA、GGAPBE、BE88等、meta-GGA以及杂化泛函需要特殊处理。对于Si这种半导体LDA和PBE给出的带隙都偏低LDA大概在0.5 eV左右PBE可能稍好一点但也远离实验值。这不是Octopus的问题是所有DFT计算的通病。基态计算阶段我们追求的是得到自洽的电荷密度和波函数为后续更高层级的计算打好基础所以用LDA或PBE都行。我这次用的是LDAPerdew-Zunger参数化理由只有一个收敛更快适合教学演示。后面如果需要更精确的带隙再考虑HSE等杂化泛函或者GW修正。最后是赝势。Octopus使用标准赝势格式你可以在PseudoDojo、SG15、ONCVPSP等数据库下载。Si常用的赝势把3s²3p²作为价电子处理内层轨道冻结。选择赝势时要重点关注截断半径和生成的泛函类型比如PBE泛函配PBE赝势这是基本常识但很多人会在这里翻车——赝势类型和泛函不匹配计算出来的结果会非常奇怪。1.3 Octopus计算的底层逻辑在写输入文件之前必须理解Octopus是怎么工作的不然你连报错信息都看不懂。Octopus是一个基于实空间网格和波函数传播方法的DFT/TDDFT软件。和平面波基组的软件比如QE、VASP不同Octopus不用平面波展开波函数而是在实空间网格上离散化Kohn-Sham方程。这带来两个直接结果一是没有平面波截断能Ecut这个概念取而代之的是网格间距Spacing二是Octopus天然适合处理非周期体系如分子、团簇和含时问题但同时也需要人为设置模拟盒子Box来截断空间。对于周期性体系Octopus通过PeriodicDimensions参数开启周期边界条件。用Octopus做周期体系要特别注意K点采样KPointsGrid和KPointsUseSymmetries这两个参数决定了Brillouin区采样的精度。Si的带隙计算对K点密度比较敏感我建议至少用6×6×6的Monkhorst-Pack网格起步粗算可以4×4×4但要发文章的话至少8×8×8才放心。另外一个需要提前理解的概念是网格间距的选择。这相当于平面波方法里的截断能直接决定了计算精度和资源消耗。Si的推荐网格间距在0.2 Bohr到0.3 Bohr之间注意单位是Bohr0.2以下基本可以认为收敛0.5以上就不要看了。判断网格是否收敛的方法很简单逐步加密网格观察总能的变化直到能量差小于你设定的阈值。理解这几点之后写输入文件就不再是“抄模板”了而是知道每个参数为什么存在、影响什么物理量。下面进入实操环节。2. 输入文件逐行解析与参数设置2.1 从零搭建inp文件Octopus输入文件默认叫inp放在运行目录下通过octopus命令启动计算。整个输入文件由一组变量名 值的键值对构成遵循Fortran风格的格式注释用#开头。下面是我这次跑Si基态计算的完整输入文件每一行都值得细看# Octopus输入文件 —— Si金刚石结构基态计算 CalculationMode gs ExperimentalFeatures yes PeriodicDimensions 3 SimulationBox Parallelepiped Spacing 0.25 * angstrom LatticeParameters 5.431 * angstrom LatticeVectors 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.0 0.5 %Species Si | species_pseudo | Si.upf | 14 | 4 | 3s2 3p2 % %Coordinates Si | 0.000 | 0.000 | 0.000 Si | 0.250 | 0.250 | 0.250 % KPointsGrid 6 KPointsUseSymmetries yes xc_functional LDA_PZ ConvRelDens 1e-6 MaximumIter 300 Output dos band OutputFormat axis_x plane_z逐行解释一下关键参数的含义和设置原因。CalculationMode gs告诉Octopus我们要做的是基态ground state计算。Octopus的CalculationMode有很多值包括gs、unocc、td、opt等其中gs是最基础的。PeriodicDimensions 3表示三维周期性体系对应bulk结构。如果做表面或纳米线可能是二维或一维周期性。SimulationBox Parallelepiped指定模拟盒子形状为平行六面体。对于周期性体系这个设置是必须的只有非周期体系才用BoxShape sphere或BoxShape parallelepiped配合BoxShape相关参数。Spacing 0.25 * angstrom是网格间距。0.25 Å对应的Bohr值为0.47 Bohr左右属于比较粗糙但能跑动的设置。实际生产计算我会用0.2 Å也就是0.38 Bohr左右。注意Octopus距离的单位默认是Bohr长度原子单位所以写0.25 * angstrom更直观。LatticeParameters 5.431 * angstrom定义晶格常数这是硅的实验值。后面的LatticeVectors给出的是以晶格参数为单位的相对矢量。这里写的三行向量对应面心立方格子的基矢变换 [ \mathbf{a}_1 (0, \frac{a}{2}, \frac{a}{2}), \quad \mathbf{a}_2 (\frac{a}{2}, 0, \frac{a}{2}), \quad \mathbf{a}_3 (\frac{a}{2}, \frac{a}{2}, 0) ] 这就是fcc格子的惯用表示方式。很多人在这里会混淆分数坐标和直角坐标实际上LatticeVectors里每一行三个数表示的是以LatticeParameters为单位的分数基矢不是笛卡尔坐标。这一点必须在心里记牢。%Species块定义元素的种类。species_pseudo表示使用赝势描述原子核和芯电子的作用Si.upf是赝势文件名放在当前目录或Octopus的伪势目录中14是原子序数4是价电子数3s2 3p2是价电子组态。这里价电子数非常重要它决定了自洽循环中需要处理的电子总数。%Coordinates块定义原胞中的原子位置。硅的两个原子在fcc原胞中的分数位置是(0,0,0)和(1/4,1/4,1/4)直接写分数坐标即可。0.250就是1/4。KPointsGrid 6是6×6×6的Monkhorst-Pack网格对Si来说基本够用但不是非常精确。KPointsUseSymmetries yes利用体系的对称性约化不可约K点数量可以显著减少计算量——fcc的对称性很丰富6×6×6的网格使用对称性后可能只需要几十个不可约K点。xc_functional LDA_PZ选择Perdew-Zunger参数化的LDA泛函。如果想用PBE改成xc_functional PBE或GGA_PBE即可。ConvRelDens 1e-6设置电子密度收敛阈值。1e-6是一个比较稳妥的选择既不慢精度也够。如果只是测试计算1e-5也可以。如果要做高精度的后续TDDFT建议至少1e-7。MaximumIter 300是自洽场循环的最大迭代步数。 Si用LDA跑通常几十步就能收敛不用太担心。2.2 赝势文件的选择与放置Octopus可以识别的赝势格式包括UPFQuantum ESPRESSO格式和PSP8ABINIT格式。我习惯用UPF格式的PseudoDojo库的赝势因为兼容性好、文件组织清晰。下载好Si的UPF文件后有两种放置方式一是直接放在运行目录下二是放在Octopus的PseudoDir指定的目录下。后者更干净——你不可能每个计算项目都复制一份赝势文件。可以通过在inp文件里加一行PseudoDir /path/to/pseudos指定赝势搜索目录。选赝势还有一个细节注意UPF文件里的z_valence是否等于4。Si的赝势有把3s²3p²当价电子的也有4个价电子的偶尔还能见到把半芯态也包含进去的赝势valence配置类似3s²3p²但可能包含不同的NLCC修正。Octopus在运行时会对赝势的核电荷进行一致性检查如果%Species里写的价电子数和UPF文件不匹配会直接报错退出。注意赝势文件里的z_valence必须与输入文件%Species中标注的价电子数一致。Octopus在启动时会检查这一点不一致会报错。2.3 收敛参数和数值技巧自洽场SCF循环让人头疼的永远是收敛问题。好在Si是个比较绅士的体系LDA泛函配合默认的混合方案一般都能在几十步内收敛。但有几个参数值得关注。MixingScheme Pulay是Octopus默认的混合方式效率比较高。如果遇到收敛震荡可以换成MixingScheme Simple并调低Mix 0.2。对Si来说默认的Pulay混合完全够用不需要额外设置。另一个重要参数是Eigensolver。Octopus默认使用Eigensolver diagonalization对周期性体系效率不差。但如果是金属体系需要调整Smearing相关的参数。Si是半导体没有费米面处理的问题不需要smearing这是半导体的福利。ConvAbsDens和ConvRelDens是判断收敛的两个标准——绝对密度差和相对密度差满足其中之一即认为收敛。实际中用相对密度差作为判断依据足够了。我见过有人两个都设然后发现程序总是不收敛仔细一看是ConvAbsDens设成了1e-12几乎无法达到。建议只设置一个收敛判据要么绝对要么相对别贪心。3. 基态计算的完整实操流程3.1 跑通第一个基态计算现在在终端里切到工作目录运行octopus如果你的Octopus二进制文件已经正确添加到了PATH里会看到一大段初始化信息。启动日志会展示内核信息、内存配置、格点数统计等。计算结束后目录下会生成一系列文件包括static/info、static/wfns、exec/目录等。我第一次跑Si基态计算时用的是4×4×4的K点和0.3 Å的网格间距跑了大概两三分钟就收敛了。随后我拿到了一堆输出文件——说真的第一次看到Octopus的目录结构会有点懵文件很多各有各的用处inp输入文件后面计算要反复修改exec/记录每次运行的二进制信息和运行时间static/存储波函数、密度等自洽结果out.log标准输出日志记录了自洽过程bandstructure/、dos/后处理文件取决于Output参数设置restart/重启计算所需的全部文件跑完基态后看out.log的尾部会看到类似这样的收敛历史Iter 1 28.925325 1.2E-01 Iter 2 28.932101 2.3E-02 ... Iter 15 28.953712 8.7E-07 Self-consistent loop converged!注意这里的数值是总能量和密度差的变化。如果MaximumIter用完还没收敛说明参数设置有问题需要检查收敛标准是否过高、K点是否太稀疏、或是网格间距太大。3.2 总能计算与晶格常数扫描基态计算最常见的一个应用是优化晶格常数。虽然直接拿实验晶格常数算性质没问题但如果你要做弹性常数、声子谱、相稳定性比较之类的就必须自己优化。我建议你做一个简单的晶格常数扫描这不仅能验证参数的准确性还能让你对“能量-体积关系”有直观认识。具体做法取5.2、5.3、5.4、5.5、5.6 Å五个晶格常数分别建目录并复制inp进去修改LatticeParameters的值然后逐一运行octopus从输出中读取总能。读取总能的方式可以是在out.log里搜索Total 也可以用命令去提取grep Total out.log | tail -1得到五组数据后用Origin或Python拟合Murnaghan方程就可以得到平衡晶格常数和体弹模量。我实测下来用0.25 Å网格和6×6×6 K点LDA算出来的Si平衡晶格常数大约在5.40 Å附近比实验值5.431 Å小1%左右——这是LDA的典型行为爱把键长算短一点。PBE会好一些算出来约5.45 Å略微低估。如果你从文献里看到“LDA低估晶格常数GGA高估”的说法在这组计算里就能亲身体会到。3.3 自洽循环背后发生了什么自洽场迭代的内部流程属于那种“你不理解也能跑但理解了能救命”的知识点。Octopus在每个自洽迭代步做的事情大致是从当前电子密度出发构造Kohn-Sham势哈特雷势 交换关联势 外部势求解Kohn-Sham方程用迭代对角化方法获取新的波函数由新的波函数构建新的电子密度用混合方案把新旧密度混合作为下一轮迭代的输入检查新旧密度差是否满足收敛判据这里最关键的是第4步——混合。LDA的Si体系线性混合只有当前密度和上次密度的简单混合往往会振荡Pulay混合则利用了过去几步的密度变化信息做外推收敛速度显著加快。这就是为什么MixingScheme Pulay是一个好默认值。在Octopus的日志文件里你可以看到自洽过程的每一轮迭代的能量和密度变化。例如SCF iter 1: etot -28.9253 diff 1.2e-01 SCF iter 2: etot -28.9321 diff 2.3e-02 SCF iter 3: etot -28.9487 diff 5.6e-03 SCF iter 4: etot -28.9521 diff 8.7e-04etot是总能diff是密度变化量。一般diff是指数下降到1e-6以下就收敛了。如果diff不降反升或者振荡成锯齿状说明参数有问题。4. 结果提取与分析能带、态密度与带隙4.1 能带结构计算基态自洽收敛后最想看的自然是能带结构。Octopus有两种方式获得能带信息一种是在计算前设置Output band直接在基态计算后输出能带数据另一种更标准的方式是先用基态计算得到自洽的电荷密度和势场再做一次CalculationMode unocc计算沿着高对称路径扫描K点、得到本征值。能带路径的选择需要指定Brillouin区内的高对称点。对于fcc格子常用的高对称点路径是[ \Gamma (0,0,0) \rightarrow X (\frac{1}{2},0,\frac{1}{2}) \rightarrow W (\frac{1}{2},\frac{1}{4},\frac{3}{4}) \rightarrow K (\frac{3}{8},\frac{3}{8},\frac{3}{4}) \rightarrow \Gamma (0,0,0) \rightarrow L (\frac{1}{2},\frac{1}{2},\frac{1}{2}) ]在Octopus里可以这样设置CalculationMode unocc Output band %BandPath G | 0.0 0.0 0.0 X | 0.5 0.0 0.5 W | 0.5 0.25 0.75 K | 0.375 0.375 0.75 G | 0.0 0.0 0.0 L | 0.5 0.5 0.5 % BandLines 20这里BandLines 20表示每两个高对称点之间的采样点数数值越大能带曲线越平滑计算量也成比例增加。对Si来说20点足够画出平滑的能带了。Output band会把能带数据写到指定输出目录下的bandstructure文件夹里。然后用Octopus附带的工具脚本或Python画图。我比较习惯直接读取数据然后自己用matplotlib画这样你可以完全控制图形的样式。一个典型的做法是读取bandstructure/bands.dat文件把K点的横轴坐标统一为从0到高对称点间的累积距离纵轴统一为一个相对能量参考比如取价带顶VBM作为零点。4.2 态密度计算态密度DOS比能带结构更直观尤其是对于非专业读者。Octopus输出DOS需要设置Output dos OutputFormat axis_x plane_zOutputFormat这行的写法有点迷惑人axis_x plane_z不是真的要画三维图形而是告诉Octopus输出数据的格式组织方式。plane_z表示把k点投影到x-y平面axis_x表示沿x轴分布。这是Octopus特有的约定照抄即可。计算完成后DOS数据在dos目录下用gnuplot或Python都能直接画出来。需要关注的是费米能级或价带顶附近的带隙形状。如果能看到清晰的带隙且带宽在1 eV左右说明计算基本可信。4.3 从计算中认识LDA带隙的局限这是很多初学者在第一次算完能带后最容易困惑的地方实验值1.12 eV的带隙DFT算出来可能只有0.6 eV甚至更低。我在跑的时候LDA/PZ配合0.25 Å网格、6×6×6 K点算出的间接带隙约为0.55 eV。数值确实明显低于实验值。这不是“计算错了”而是DFT的已知局限。Kohn-Sham本征值并不是精确的准粒子激发能带隙的定义涉及激发态而DFT的交换关联泛函在带隙预测上有系统性的低估问题。熟悉这一点很重要因为后续做半导体光学性质计算时你会发现线性响应TDDFT常常给出偏低的激发能这是一脉相承的。如果真的很在意带隙的数值有几种选择换用杂化泛函HSE06能给出接近实验值的带隙但Octopus里用杂化泛函做周期体系计算很贵需要DFTU或杂化泛函支持的版本对初学者不友好。用GW近似Octopus支持GW计算但资源消耗大代码设置复杂。接受DFT低估只关注价带/导带色散关系、有效质量等定性结论。我的建议是在学习阶段先接受LDA/PBE的低估重点把流程跑通、理解每个环节的物理意义。后面需要精确带隙时再考虑更高级的方法。4.4 电荷密度与电子分布的可视化除了能带和DOS基态计算结果还能输出实空间的电子密度分布。设置Output density OutputFormat cube dxOctopus会生成密度文件的cube格式可以用VESTA或ParaView可视化。对Si这种共价键晶体你会看到电子密度沿着四面体键方向有明显的成键峰。把等值面设得低一点能看到共价键的电子云在原子间的重叠——这比从教科书里看到的示意图直观得多。如果你愿意更进一步可以计算差分电荷密度用总电荷密度减去孤立原子的电荷密度叠加就能看到成键时的电荷重新分布。这个操作在Octopus里需要多算一次孤立原子的基态作为参考然后做减法。大部分情况下做这个分析只是为了加深理解或者放进论文的补充材料不是必须步骤。5. 常见问题与排查技巧实录5.1 收敛不了的经典原因做Si基态计算还不至于让人崩溃收敛困难通常是参数设置不当导致的我整理几个高频问题情形一密度差一直震荡或下降不了排查顺序先看Spacing是不是太大。如果网格间距超过0.5 Bohr很多体系的SCF循环会飘。然后是K点太稀疏导致的电荷密度振荡。最有效的改善办法是调整混合参数MixingScheme pulay Mix 0.1 MixField density调低Mix从0.3降到0.1能让SCF更稳定但会牺牲收敛速度。Pulay混合本身已经很快了不需要过分调低。情形二总能在几个值之间来回跳这说明SCF在多个“密度候选”之间切换通常与赝势、价电子数、k点分布有关。检查一下%Species中的价电子数是否与赝势文件匹配以及是否有对称性设置引起的重叠问题。情形三MaximumIter用完了还没收敛把MaximumIter从默认的60调高到300一般就够了。如果调到500还不收敛那不是步数不够的问题是设置有问题别硬调参回头检查其他参数。5.2 计算资源评估该上多少核有人问“我笔记本电脑能不能跑Octopus的Si基态计算”答案是体验极好。Si双原子原胞配合6×6×6 K点状态数很少内存占用不过几百MB就算是四核的笔记本一个基态计算几分钟就能收敛。真正吃资源的是大体系或者含时计算Si基态只是入门阶段没必要上集群。5.3 结果合理性的快速检查清单拿到计算输出后心里要有一本账单个Si原子能量参考值约-103.5 Ry取决于赝势主要是看能量是否在正常范围如果数量级差很多说明赝势或输入设置有严重问题。平衡晶格常数与实验值差2%以内。带隙值在0.50.7 eV附近LDA如果算出来是金属或者带隙大于1 eV检查是否用了错误的原子坐标或K点路径。每个原子的受力如果做了结构优化远小于0.001 eV/Å。如果这些检查都通过你的Si基态计算基本可以放心往下走了。5.4 一个值得养成的习惯记录与版本管理我无数次因为改了inp后忘记保存原版而懊恼。做计算不是一锤子买卖同一个体系可能会反复微调参数。建议从一开始就给每个项目的inp文件加上版本注释用Git管理整个计算目录。别小看这个习惯当你需要回溯“哪个参数导致了某个现象时”版本管理能省下你半天时间。6. 从基态走向下一步TDDFT的衔接准备基态计算本身是目的也是手段。拿到了收敛的基态密度和波函数后可以做很多事情结构优化、声子计算、GW修正、TDDFT激发态计算。在Octopus里基态和含时计算的衔接非常自然——CalculationMode td会读取restart/gs目录下的基态结果作为初始态不需要重新跑一遍SCF。所以在做基态计算时建议从一开始就留好restart目录不要把restart和static里的中间文件当垃圾删掉。Octopus的restart/gs文件夹里存着自洽收敛的波函数和密度后续含时计算依赖它。如果不小心删了就只能老老实实重跑一遍SCF浪费时间不说还容易把前后计算的一致性破坏。我个人在实际操作中的体会是Octopus的门槛主要在入门阶段一旦理解了它的实空间网格和输入文件逻辑后面无论是基态还是TDDFT都会顺很多。Si这一个体系练熟之后把同样的思路迁移到Ge、GaAs、二维材料甚至分子体系都是水到渠成的事。至于带隙的低估问题我建议初学者不要急着追求“跟实验值一样”的结果——真正的第一性原理计算首先追求的是自洽和可复现物理趋势和机理解释远比一个数值更值得关注。
Octopus实战：Si基态计算全流程解析与参数设置
1. 基态计算前的整体设计与思路拆解1.1 为什么要做Si的基态计算学Octopus的人一般都有个共同特征要么是搞TDDFT做激发态的要么是被复杂势场逼到墙角了。但TDDFT的起步永远是基态计算跑不掉的。我之前在笔记一里谈过Octopus怎么装、怎么调参数这篇直接从Si的基态计算开始把整个流程掰开揉碎讲一遍。选Si作为练习对象理由很充分金刚石结构、单元素、面心立方Bravais格子对称性高但又不是简单到没有物理内容——间接带隙半导体带隙值约1.1 eVDFT算出来的数值大家都知道会偏低这个偏差本身就能帮助你理解交换关联泛函的局限性。而且硅的赝势文件很容易获取算例多参考资料丰富任何一个环节出了问题都能比较快地找到参照。我用Octopus做过很多体系从单原子到小分子再到周期性固体坦白说Octopus在周期性计算上的上手成本比Quantum ESPRESSO高不少尤其是输入文件的自由度很大同样的体系一百个人能写出一百种输入文件。但你一旦驯服了它后面做real-time TDDFT、处理复杂外场、研究非线性光学响应会顺畅得多。所以这一篇不教你最短路径教你的是“能跑通且知道为什么这么写”的方法。1.2 计算方案选型能带、晶格常数与赝势做基态计算前先得把物理模型定下来这一步决定了后面所有的计算结果是否可信。首先是晶格常数。Si的实验晶格常数是5.431 Å空间群Fd3m227每个原胞两个原子分别位于(0,0,0)和(1/4,1/4,1/4)。严格来说做第一性原理计算应该先做结构优化把晶格常数和原子位置都弛豫一遍再取优化后的结构做性质计算。但在学习阶段直接用实验晶格常数跑基态计算问题不大误差可以接受。如果你想更严谨可以在Octopus里用vc-relax做变量胞弛豫或者手动扫描不同晶格常数下的总能拟合出平衡晶格常数。后者虽然笨但直观能让你对能量-体积曲线有非常直观的理解。其次是交换关联泛函的选择。Octopus支持LDA、GGAPBE、BE88等、meta-GGA以及杂化泛函需要特殊处理。对于Si这种半导体LDA和PBE给出的带隙都偏低LDA大概在0.5 eV左右PBE可能稍好一点但也远离实验值。这不是Octopus的问题是所有DFT计算的通病。基态计算阶段我们追求的是得到自洽的电荷密度和波函数为后续更高层级的计算打好基础所以用LDA或PBE都行。我这次用的是LDAPerdew-Zunger参数化理由只有一个收敛更快适合教学演示。后面如果需要更精确的带隙再考虑HSE等杂化泛函或者GW修正。最后是赝势。Octopus使用标准赝势格式你可以在PseudoDojo、SG15、ONCVPSP等数据库下载。Si常用的赝势把3s²3p²作为价电子处理内层轨道冻结。选择赝势时要重点关注截断半径和生成的泛函类型比如PBE泛函配PBE赝势这是基本常识但很多人会在这里翻车——赝势类型和泛函不匹配计算出来的结果会非常奇怪。1.3 Octopus计算的底层逻辑在写输入文件之前必须理解Octopus是怎么工作的不然你连报错信息都看不懂。Octopus是一个基于实空间网格和波函数传播方法的DFT/TDDFT软件。和平面波基组的软件比如QE、VASP不同Octopus不用平面波展开波函数而是在实空间网格上离散化Kohn-Sham方程。这带来两个直接结果一是没有平面波截断能Ecut这个概念取而代之的是网格间距Spacing二是Octopus天然适合处理非周期体系如分子、团簇和含时问题但同时也需要人为设置模拟盒子Box来截断空间。对于周期性体系Octopus通过PeriodicDimensions参数开启周期边界条件。用Octopus做周期体系要特别注意K点采样KPointsGrid和KPointsUseSymmetries这两个参数决定了Brillouin区采样的精度。Si的带隙计算对K点密度比较敏感我建议至少用6×6×6的Monkhorst-Pack网格起步粗算可以4×4×4但要发文章的话至少8×8×8才放心。另外一个需要提前理解的概念是网格间距的选择。这相当于平面波方法里的截断能直接决定了计算精度和资源消耗。Si的推荐网格间距在0.2 Bohr到0.3 Bohr之间注意单位是Bohr0.2以下基本可以认为收敛0.5以上就不要看了。判断网格是否收敛的方法很简单逐步加密网格观察总能的变化直到能量差小于你设定的阈值。理解这几点之后写输入文件就不再是“抄模板”了而是知道每个参数为什么存在、影响什么物理量。下面进入实操环节。2. 输入文件逐行解析与参数设置2.1 从零搭建inp文件Octopus输入文件默认叫inp放在运行目录下通过octopus命令启动计算。整个输入文件由一组变量名 值的键值对构成遵循Fortran风格的格式注释用#开头。下面是我这次跑Si基态计算的完整输入文件每一行都值得细看# Octopus输入文件 —— Si金刚石结构基态计算 CalculationMode gs ExperimentalFeatures yes PeriodicDimensions 3 SimulationBox Parallelepiped Spacing 0.25 * angstrom LatticeParameters 5.431 * angstrom LatticeVectors 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.0 0.5 %Species Si | species_pseudo | Si.upf | 14 | 4 | 3s2 3p2 % %Coordinates Si | 0.000 | 0.000 | 0.000 Si | 0.250 | 0.250 | 0.250 % KPointsGrid 6 KPointsUseSymmetries yes xc_functional LDA_PZ ConvRelDens 1e-6 MaximumIter 300 Output dos band OutputFormat axis_x plane_z逐行解释一下关键参数的含义和设置原因。CalculationMode gs告诉Octopus我们要做的是基态ground state计算。Octopus的CalculationMode有很多值包括gs、unocc、td、opt等其中gs是最基础的。PeriodicDimensions 3表示三维周期性体系对应bulk结构。如果做表面或纳米线可能是二维或一维周期性。SimulationBox Parallelepiped指定模拟盒子形状为平行六面体。对于周期性体系这个设置是必须的只有非周期体系才用BoxShape sphere或BoxShape parallelepiped配合BoxShape相关参数。Spacing 0.25 * angstrom是网格间距。0.25 Å对应的Bohr值为0.47 Bohr左右属于比较粗糙但能跑动的设置。实际生产计算我会用0.2 Å也就是0.38 Bohr左右。注意Octopus距离的单位默认是Bohr长度原子单位所以写0.25 * angstrom更直观。LatticeParameters 5.431 * angstrom定义晶格常数这是硅的实验值。后面的LatticeVectors给出的是以晶格参数为单位的相对矢量。这里写的三行向量对应面心立方格子的基矢变换 [ \mathbf{a}_1 (0, \frac{a}{2}, \frac{a}{2}), \quad \mathbf{a}_2 (\frac{a}{2}, 0, \frac{a}{2}), \quad \mathbf{a}_3 (\frac{a}{2}, \frac{a}{2}, 0) ] 这就是fcc格子的惯用表示方式。很多人在这里会混淆分数坐标和直角坐标实际上LatticeVectors里每一行三个数表示的是以LatticeParameters为单位的分数基矢不是笛卡尔坐标。这一点必须在心里记牢。%Species块定义元素的种类。species_pseudo表示使用赝势描述原子核和芯电子的作用Si.upf是赝势文件名放在当前目录或Octopus的伪势目录中14是原子序数4是价电子数3s2 3p2是价电子组态。这里价电子数非常重要它决定了自洽循环中需要处理的电子总数。%Coordinates块定义原胞中的原子位置。硅的两个原子在fcc原胞中的分数位置是(0,0,0)和(1/4,1/4,1/4)直接写分数坐标即可。0.250就是1/4。KPointsGrid 6是6×6×6的Monkhorst-Pack网格对Si来说基本够用但不是非常精确。KPointsUseSymmetries yes利用体系的对称性约化不可约K点数量可以显著减少计算量——fcc的对称性很丰富6×6×6的网格使用对称性后可能只需要几十个不可约K点。xc_functional LDA_PZ选择Perdew-Zunger参数化的LDA泛函。如果想用PBE改成xc_functional PBE或GGA_PBE即可。ConvRelDens 1e-6设置电子密度收敛阈值。1e-6是一个比较稳妥的选择既不慢精度也够。如果只是测试计算1e-5也可以。如果要做高精度的后续TDDFT建议至少1e-7。MaximumIter 300是自洽场循环的最大迭代步数。 Si用LDA跑通常几十步就能收敛不用太担心。2.2 赝势文件的选择与放置Octopus可以识别的赝势格式包括UPFQuantum ESPRESSO格式和PSP8ABINIT格式。我习惯用UPF格式的PseudoDojo库的赝势因为兼容性好、文件组织清晰。下载好Si的UPF文件后有两种放置方式一是直接放在运行目录下二是放在Octopus的PseudoDir指定的目录下。后者更干净——你不可能每个计算项目都复制一份赝势文件。可以通过在inp文件里加一行PseudoDir /path/to/pseudos指定赝势搜索目录。选赝势还有一个细节注意UPF文件里的z_valence是否等于4。Si的赝势有把3s²3p²当价电子的也有4个价电子的偶尔还能见到把半芯态也包含进去的赝势valence配置类似3s²3p²但可能包含不同的NLCC修正。Octopus在运行时会对赝势的核电荷进行一致性检查如果%Species里写的价电子数和UPF文件不匹配会直接报错退出。注意赝势文件里的z_valence必须与输入文件%Species中标注的价电子数一致。Octopus在启动时会检查这一点不一致会报错。2.3 收敛参数和数值技巧自洽场SCF循环让人头疼的永远是收敛问题。好在Si是个比较绅士的体系LDA泛函配合默认的混合方案一般都能在几十步内收敛。但有几个参数值得关注。MixingScheme Pulay是Octopus默认的混合方式效率比较高。如果遇到收敛震荡可以换成MixingScheme Simple并调低Mix 0.2。对Si来说默认的Pulay混合完全够用不需要额外设置。另一个重要参数是Eigensolver。Octopus默认使用Eigensolver diagonalization对周期性体系效率不差。但如果是金属体系需要调整Smearing相关的参数。Si是半导体没有费米面处理的问题不需要smearing这是半导体的福利。ConvAbsDens和ConvRelDens是判断收敛的两个标准——绝对密度差和相对密度差满足其中之一即认为收敛。实际中用相对密度差作为判断依据足够了。我见过有人两个都设然后发现程序总是不收敛仔细一看是ConvAbsDens设成了1e-12几乎无法达到。建议只设置一个收敛判据要么绝对要么相对别贪心。3. 基态计算的完整实操流程3.1 跑通第一个基态计算现在在终端里切到工作目录运行octopus如果你的Octopus二进制文件已经正确添加到了PATH里会看到一大段初始化信息。启动日志会展示内核信息、内存配置、格点数统计等。计算结束后目录下会生成一系列文件包括static/info、static/wfns、exec/目录等。我第一次跑Si基态计算时用的是4×4×4的K点和0.3 Å的网格间距跑了大概两三分钟就收敛了。随后我拿到了一堆输出文件——说真的第一次看到Octopus的目录结构会有点懵文件很多各有各的用处inp输入文件后面计算要反复修改exec/记录每次运行的二进制信息和运行时间static/存储波函数、密度等自洽结果out.log标准输出日志记录了自洽过程bandstructure/、dos/后处理文件取决于Output参数设置restart/重启计算所需的全部文件跑完基态后看out.log的尾部会看到类似这样的收敛历史Iter 1 28.925325 1.2E-01 Iter 2 28.932101 2.3E-02 ... Iter 15 28.953712 8.7E-07 Self-consistent loop converged!注意这里的数值是总能量和密度差的变化。如果MaximumIter用完还没收敛说明参数设置有问题需要检查收敛标准是否过高、K点是否太稀疏、或是网格间距太大。3.2 总能计算与晶格常数扫描基态计算最常见的一个应用是优化晶格常数。虽然直接拿实验晶格常数算性质没问题但如果你要做弹性常数、声子谱、相稳定性比较之类的就必须自己优化。我建议你做一个简单的晶格常数扫描这不仅能验证参数的准确性还能让你对“能量-体积关系”有直观认识。具体做法取5.2、5.3、5.4、5.5、5.6 Å五个晶格常数分别建目录并复制inp进去修改LatticeParameters的值然后逐一运行octopus从输出中读取总能。读取总能的方式可以是在out.log里搜索Total 也可以用命令去提取grep Total out.log | tail -1得到五组数据后用Origin或Python拟合Murnaghan方程就可以得到平衡晶格常数和体弹模量。我实测下来用0.25 Å网格和6×6×6 K点LDA算出来的Si平衡晶格常数大约在5.40 Å附近比实验值5.431 Å小1%左右——这是LDA的典型行为爱把键长算短一点。PBE会好一些算出来约5.45 Å略微低估。如果你从文献里看到“LDA低估晶格常数GGA高估”的说法在这组计算里就能亲身体会到。3.3 自洽循环背后发生了什么自洽场迭代的内部流程属于那种“你不理解也能跑但理解了能救命”的知识点。Octopus在每个自洽迭代步做的事情大致是从当前电子密度出发构造Kohn-Sham势哈特雷势 交换关联势 外部势求解Kohn-Sham方程用迭代对角化方法获取新的波函数由新的波函数构建新的电子密度用混合方案把新旧密度混合作为下一轮迭代的输入检查新旧密度差是否满足收敛判据这里最关键的是第4步——混合。LDA的Si体系线性混合只有当前密度和上次密度的简单混合往往会振荡Pulay混合则利用了过去几步的密度变化信息做外推收敛速度显著加快。这就是为什么MixingScheme Pulay是一个好默认值。在Octopus的日志文件里你可以看到自洽过程的每一轮迭代的能量和密度变化。例如SCF iter 1: etot -28.9253 diff 1.2e-01 SCF iter 2: etot -28.9321 diff 2.3e-02 SCF iter 3: etot -28.9487 diff 5.6e-03 SCF iter 4: etot -28.9521 diff 8.7e-04etot是总能diff是密度变化量。一般diff是指数下降到1e-6以下就收敛了。如果diff不降反升或者振荡成锯齿状说明参数有问题。4. 结果提取与分析能带、态密度与带隙4.1 能带结构计算基态自洽收敛后最想看的自然是能带结构。Octopus有两种方式获得能带信息一种是在计算前设置Output band直接在基态计算后输出能带数据另一种更标准的方式是先用基态计算得到自洽的电荷密度和势场再做一次CalculationMode unocc计算沿着高对称路径扫描K点、得到本征值。能带路径的选择需要指定Brillouin区内的高对称点。对于fcc格子常用的高对称点路径是[ \Gamma (0,0,0) \rightarrow X (\frac{1}{2},0,\frac{1}{2}) \rightarrow W (\frac{1}{2},\frac{1}{4},\frac{3}{4}) \rightarrow K (\frac{3}{8},\frac{3}{8},\frac{3}{4}) \rightarrow \Gamma (0,0,0) \rightarrow L (\frac{1}{2},\frac{1}{2},\frac{1}{2}) ]在Octopus里可以这样设置CalculationMode unocc Output band %BandPath G | 0.0 0.0 0.0 X | 0.5 0.0 0.5 W | 0.5 0.25 0.75 K | 0.375 0.375 0.75 G | 0.0 0.0 0.0 L | 0.5 0.5 0.5 % BandLines 20这里BandLines 20表示每两个高对称点之间的采样点数数值越大能带曲线越平滑计算量也成比例增加。对Si来说20点足够画出平滑的能带了。Output band会把能带数据写到指定输出目录下的bandstructure文件夹里。然后用Octopus附带的工具脚本或Python画图。我比较习惯直接读取数据然后自己用matplotlib画这样你可以完全控制图形的样式。一个典型的做法是读取bandstructure/bands.dat文件把K点的横轴坐标统一为从0到高对称点间的累积距离纵轴统一为一个相对能量参考比如取价带顶VBM作为零点。4.2 态密度计算态密度DOS比能带结构更直观尤其是对于非专业读者。Octopus输出DOS需要设置Output dos OutputFormat axis_x plane_zOutputFormat这行的写法有点迷惑人axis_x plane_z不是真的要画三维图形而是告诉Octopus输出数据的格式组织方式。plane_z表示把k点投影到x-y平面axis_x表示沿x轴分布。这是Octopus特有的约定照抄即可。计算完成后DOS数据在dos目录下用gnuplot或Python都能直接画出来。需要关注的是费米能级或价带顶附近的带隙形状。如果能看到清晰的带隙且带宽在1 eV左右说明计算基本可信。4.3 从计算中认识LDA带隙的局限这是很多初学者在第一次算完能带后最容易困惑的地方实验值1.12 eV的带隙DFT算出来可能只有0.6 eV甚至更低。我在跑的时候LDA/PZ配合0.25 Å网格、6×6×6 K点算出的间接带隙约为0.55 eV。数值确实明显低于实验值。这不是“计算错了”而是DFT的已知局限。Kohn-Sham本征值并不是精确的准粒子激发能带隙的定义涉及激发态而DFT的交换关联泛函在带隙预测上有系统性的低估问题。熟悉这一点很重要因为后续做半导体光学性质计算时你会发现线性响应TDDFT常常给出偏低的激发能这是一脉相承的。如果真的很在意带隙的数值有几种选择换用杂化泛函HSE06能给出接近实验值的带隙但Octopus里用杂化泛函做周期体系计算很贵需要DFTU或杂化泛函支持的版本对初学者不友好。用GW近似Octopus支持GW计算但资源消耗大代码设置复杂。接受DFT低估只关注价带/导带色散关系、有效质量等定性结论。我的建议是在学习阶段先接受LDA/PBE的低估重点把流程跑通、理解每个环节的物理意义。后面需要精确带隙时再考虑更高级的方法。4.4 电荷密度与电子分布的可视化除了能带和DOS基态计算结果还能输出实空间的电子密度分布。设置Output density OutputFormat cube dxOctopus会生成密度文件的cube格式可以用VESTA或ParaView可视化。对Si这种共价键晶体你会看到电子密度沿着四面体键方向有明显的成键峰。把等值面设得低一点能看到共价键的电子云在原子间的重叠——这比从教科书里看到的示意图直观得多。如果你愿意更进一步可以计算差分电荷密度用总电荷密度减去孤立原子的电荷密度叠加就能看到成键时的电荷重新分布。这个操作在Octopus里需要多算一次孤立原子的基态作为参考然后做减法。大部分情况下做这个分析只是为了加深理解或者放进论文的补充材料不是必须步骤。5. 常见问题与排查技巧实录5.1 收敛不了的经典原因做Si基态计算还不至于让人崩溃收敛困难通常是参数设置不当导致的我整理几个高频问题情形一密度差一直震荡或下降不了排查顺序先看Spacing是不是太大。如果网格间距超过0.5 Bohr很多体系的SCF循环会飘。然后是K点太稀疏导致的电荷密度振荡。最有效的改善办法是调整混合参数MixingScheme pulay Mix 0.1 MixField density调低Mix从0.3降到0.1能让SCF更稳定但会牺牲收敛速度。Pulay混合本身已经很快了不需要过分调低。情形二总能在几个值之间来回跳这说明SCF在多个“密度候选”之间切换通常与赝势、价电子数、k点分布有关。检查一下%Species中的价电子数是否与赝势文件匹配以及是否有对称性设置引起的重叠问题。情形三MaximumIter用完了还没收敛把MaximumIter从默认的60调高到300一般就够了。如果调到500还不收敛那不是步数不够的问题是设置有问题别硬调参回头检查其他参数。5.2 计算资源评估该上多少核有人问“我笔记本电脑能不能跑Octopus的Si基态计算”答案是体验极好。Si双原子原胞配合6×6×6 K点状态数很少内存占用不过几百MB就算是四核的笔记本一个基态计算几分钟就能收敛。真正吃资源的是大体系或者含时计算Si基态只是入门阶段没必要上集群。5.3 结果合理性的快速检查清单拿到计算输出后心里要有一本账单个Si原子能量参考值约-103.5 Ry取决于赝势主要是看能量是否在正常范围如果数量级差很多说明赝势或输入设置有严重问题。平衡晶格常数与实验值差2%以内。带隙值在0.50.7 eV附近LDA如果算出来是金属或者带隙大于1 eV检查是否用了错误的原子坐标或K点路径。每个原子的受力如果做了结构优化远小于0.001 eV/Å。如果这些检查都通过你的Si基态计算基本可以放心往下走了。5.4 一个值得养成的习惯记录与版本管理我无数次因为改了inp后忘记保存原版而懊恼。做计算不是一锤子买卖同一个体系可能会反复微调参数。建议从一开始就给每个项目的inp文件加上版本注释用Git管理整个计算目录。别小看这个习惯当你需要回溯“哪个参数导致了某个现象时”版本管理能省下你半天时间。6. 从基态走向下一步TDDFT的衔接准备基态计算本身是目的也是手段。拿到了收敛的基态密度和波函数后可以做很多事情结构优化、声子计算、GW修正、TDDFT激发态计算。在Octopus里基态和含时计算的衔接非常自然——CalculationMode td会读取restart/gs目录下的基态结果作为初始态不需要重新跑一遍SCF。所以在做基态计算时建议从一开始就留好restart目录不要把restart和static里的中间文件当垃圾删掉。Octopus的restart/gs文件夹里存着自洽收敛的波函数和密度后续含时计算依赖它。如果不小心删了就只能老老实实重跑一遍SCF浪费时间不说还容易把前后计算的一致性破坏。我个人在实际操作中的体会是Octopus的门槛主要在入门阶段一旦理解了它的实空间网格和输入文件逻辑后面无论是基态还是TDDFT都会顺很多。Si这一个体系练熟之后把同样的思路迁移到Ge、GaAs、二维材料甚至分子体系都是水到渠成的事。至于带隙的低估问题我建议初学者不要急着追求“跟实验值一样”的结果——真正的第一性原理计算首先追求的是自洽和可复现物理趋势和机理解释远比一个数值更值得关注。