CP2K电荷密度与分子轨道

下载20份完整输入、真实输出、网格和绘图脚本。在已分配的CPU资源中运行 bash run.sh 示例名称；输出在 work/示例名称/，原结果在 reference/。不修改公共软件、不自动提交任务；正式提交流程参见CP2K入门。

1 在DFT中请求网格输出#

bash run.sh water-orbitals

完整输入中的DFT/PRINT段：

&PRINT
 &E_DENSITY_CUBE
  STRIDE 2 2 2
 &END E_DENSITY_CUBE
 &MO_CUBES
  NHOMO 1
  NLUMO 1
  STRIDE 2 2 2
 &END MO_CUBES
&END PRINT

本例价电子数8，闭壳层占据4个轨道，输出密度、HOMO第4轨道和LUMO第5轨道。文件名中的编号不是原子编号，也不是轨道能量。

ls work/water-orbitals/*.cube

2 先确认网格单位#

本例CUBE头中每方向75点，网格矢量约0.302356 Bohr；正的网格数对应Bohr坐标约定。电子密度单位e/Bohr³，轨道振幅单位Bohr⁻³ᐟ²。

STRIDE=2令输出网格比内部网格稀疏。直接把数据当e/ų积分会得到错误的电子数；若转换到Å坐标，密度和体积单位必须一起转换。

3 用实际积分核对输出#

包内 analyze.py 的 cube()处理CP2K原始数值字段，保留正负号。有些相邻负数没有空格，例如 0.29588E-018-0.26774E-018，不能简单按空格split后当作两列。

origin, grid_vectors, rho = cube(filename)
voxel_volume = abs(np.linalg.det(grid_vectors))
N = rho.sum() * voxel_volume

这里grid_vectors仍以Bohr计。对波函数应改为 sum(psi**2)*voxel_volume，不能直接积分有正负号的psi。

密度接近8个价电子、两个轨道接近归一化到1；偏差与稀疏输出网格、截断和文本精度有关，不声称数学上严格等于整数。需要更精细积分或等值面时应减小STRIDE并重新核对。

4 两种颜色表示什么#

密度图是log10显示。轨道图的红/蓝表示波函数的符号，整体乘−1不改变物理态；不是正负电荷，也不是上下自旋。不同轨道使用各自对称色阶，不能直接按颜色深浅比较两个轨道的总电子数。

HOMO为便于看节点，选取穿过轨道的y≈6.08 Å切面；密度和LUMO使用x≈6.08 Å切面。各图横轴和标题标出方向，图上原子为相应投影。单个切面不能替代完整三维等值面判断。

5 重绘与替换自己的数据#

完整包全部任务结束后，在解压目录：

bash analyze.sh

输出图在work/。只运行单个轨道示例时，可复用analyze.py中的CUBE读取和积分函数，指向自己的文件。先保持数据单位、晶胞、原子位置和选择的切面一致，再调整色阶与展示范围。

这里是GTH价电子密度，不是包含相同芯电子表示的全电子密度；需要核心区性质时应选择适合的GAPW/全电子设置，不能仅靠放大这张图获得缺失的芯电子信息。