OVITO电荷积分与三维表面

下载输入、完整脚本和实际结果，在新目录解压后运行：

bash run.sh

分析使用OVITO3.16.1的Python接口，作图使用ASE3.29.0环境中的Matplotlib。共享环境提供预编译程序，用户无需自己编译。inputs/保存结构和轨迹，reference/保存参考数据，新结果写入work/。下面的单项命令在完整运行之后执行。

1. 计算电荷与参考密度的区别#

本例复用集群已有NaCl实算的64×64×64网格：价电子密度来自CHGCAR，全电子重构参考密度来自AECCAR0+AECCAR2。包内NaCl-charge-grid.npz保留读入后的数值、原胞、位置与元素类型，来源文件哈希在SOURCES.json；不需要重新运行VASP。

module purge
module load ovito/3.16.1
cd work
python volume_surface.py
cat results/NaCl-charge.csv

volume_surface.py中的source()恢复体素网格。随后执行：

from ovito.modifiers import BaderChargeIntegrationModifier
p = source()
p.modifiers.append(BaderChargeIntegrationModifier(
    operate_on='charge-density',
    density='Reference',
    integrand='Density'))
d = p.compute()
print(d.particles['Integrated Density'])
print(d.particles['Integrated Volume'])

density用来确定分区，integrand才是逐区积分的量。若直接导入自己的CHGCAR，可先打印p.compute().grids及各网格的属性名称；不要凭本例的重命名字段猜自己的输入字段。

两原子共8.000000个价电子，分区总体积46.298250 ų。这里赝势价电子数是Cl 7、Na 1，所以净电荷近似为Cl −0.833328e、Na +0.833328e。若换了Na_pv等赝势，必须换成对应ZVAL，不能继续减1。

图左是分数坐标c=0.5的密度截面，以a、b方向分数坐标展示；原胞为非正交晶胞，图的方框不代表实际正方形截面。两个程序采用的分区算法不同，本例最大差约0.00226个电子。结果展示了操作流程，并未证明64网格已经收敛；正式比较转移电荷应增加网格并保持赝势、结构和分区参考一致。

2. 创建电荷等值面#

from ovito.modifiers import CreateIsosurfaceModifier
from ovito.io import export_file
p = source()
p.modifiers.append(CreateIsosurfaceModifier(
    operate_on='charge-density', property='Density',
    isolevel=0.2, identify_regions=True))
d = p.compute()
export_file(d, 'density.vtk', 'vtk/trimesh', key='isosurface')

本例阈值0.2 e/ų，生成11052个顶点的表面。输出results/NaCl-density.vtk可在支持VTK三角网格的可视化软件中打开；等值面大小会随阈值变化，不能把一个阈值下的表面当成原子物理半径。

3. 原子团簇的几何表面#

同一脚本还对非周期Al团簇使用ConstructSurfaceModifier(radius=3, smoothing_level=0, identify_regions=True)。本例108原子，几何表面积604.694982 Ų、包围体积1032.434859 ų；结果保存在results/Al-cluster.vtk。

这是alpha-shape几何表面，半径参数不同会改变结果。非周期系统外部空间无界，因此“整个空间的填充分数”等量没有有限数值，本包将这类非有限统计记为null，不把它们画成有限材料性质。

参考：OVITO Python处理器说明。