py4vasp电荷密度与MD轨迹提取

下载本篇输入、完整脚本、数据与高清图。解压后在新目录运行，reference保存鸢算科技实际运行的参考结果。

1 电荷密度保存在哪#

先按HDF5读取教程运行下载包。NaCl例子在VASP中设置LH5=.TRUE.、LCHARG=.TRUE.，密度保存在vaspwave.h5，结构和其他结果在vaspout.h5。

import py4vasp
import numpy as np
c = py4vasp.Calculation.from_path('input/nacl-charge-hdf5')
d = c.density.read()
q = d['charge']
volume = abs(np.linalg.det(d['structure']['lattice_vectors']))
rho = q / volume
print(q.shape, q.mean())

本次网格为64×64×64。该读取数组沿用“密度乘晶胞体积”的归一化，数组平均值为 7.9999998231，即8个价电子；除以46.29825 ų得到e/ų密度。不要将数组元素直接标成e/ų，也不要再除一次网格点总数。

这是PAW价电子密度，不含全部芯电子。该粗网格下局部最小值约−0.00975 e/ų，原始数值保留，不应为得到好看的图而直接截为零。定量密度拓扑研究需检查密度类型和网格；Bader的全电子参考方法见Bader教程。

2 读取全部MD帧#

默认read()通常返回最后一步；读取轨迹需要使用切片：

c = py4vasp.Calculation.from_path('input/si-md-short-hdf5')
e = c.energy[:].read()
f = c.force[:].read()['forces']
s = c.stress[:].read()['stress']
positions = c.structure[:].read()['positions']
print(f.shape, s.shape)

本例为2原子Si、20步、步长1 fs的短NVE流程，使用3×3×3 k点，电子收敛1E-6。力数组为(20,2,3)，应力数组为(20,3,3)。

关键INCAR设置为：

IBRION = 0
NSW = 20
POTIM = 1
TEBEG = 300
TEEND = 300
SMASS = -3
ISYM = 0
LH5 = .TRUE.

从energy字典读取total energy ETOTAL、temperature TEIN等键，脚本输出md.csv。

总能跨度约 0.000101 eV/胞。瞬时温度在约1.34–300 K之间显著变化，这是极小体系动能与势能交换的短演示，不代表已平衡的300 K材料性质。正式AIMD需要更大体系、合适统计系综、足够平衡和采样，且应重新收敛k点及步长。

参考：py4vasp密度接口、py4vasp接口说明。