ABACUS 恒压与晶胞动力学实操

下载13个条件的完整输入、原生轨迹、续算状态与绘图脚本，解压进入附件根目录。先用已有数据生成图表：

bash analyze.sh

结果在reference/analysis、图片在reference/figures、ASE轨迹在reference/export。从头执行整套流程：

sbatch submit.slurm

脚本在新的work中依次运行13个条件，已有work时请新解压，避免覆盖已完成结果。想先跑一个条件，可在计算资源环境中加载下述模块、运行setup.py后用run_case.py指定目录名；普通用户不要在登录节点运行DFT。提交模板申请1张P100、1个MPI进程、1个线程，套餐用户按实际权限用yskit生成脚本。

本例ABACUS3.10.1模块为abacus/3.10.1-p100，GNU13.3、HPC-X2.50/OpenMPI5.0.10rc2、CUDA12.8环境，实际为P100双精度；后处理为ase/3.29.0。模型采用32原子fcc铝、3价电子LDA赝势Al.pz-vbc.UPF及PZ泛函，40Ry、Γ点、64条带，共96个价电子。不能用同名11价电子赝势替换后仍保留这些带数与精度设置。

这批轨迹只有40fs，用来指导输入、读取和续算，不用于扩散系数、熔点、平衡密度或收敛热力学性质。

1 先理解目标压力与瞬时压力#

本例仍从8.1Å的Al32晶胞出发，目标300K、0kbar。初始实际总压力约−80.916kbar（−8.092GPa），明显偏离目标，下面展示的是压浴响应过程，不能读成“0压平衡晶胞”。生产计算应先优化和预平衡，再采集统计。

完整等向NPT输入：

INPUT_PARAMETERS
suffix tutorial
pseudo_dir ../../PP_ORB
ntype 1
basis_type pw
calculation md
ecutwfc 40
nbands 64
ks_solver cg
symmetry 0
device gpu
precision double
scf_thr 1e-7
scf_nmax 100
smearing_method gaussian
smearing_sigma 0.01
cal_force 1
cal_stress 1
out_chg 0
md_type npt
md_nstep 80
md_dt 0.5
md_tfirst 300
md_seed 12345
md_dumpfreq 1
md_restartfreq 40
dft_functional pz
init_vel 1
md_tlast 300
dump_force 1
dump_vel 1
dump_virial 1
md_thermostat nhc
md_pmode iso
md_pcouple xyz
md_pfirst 0
md_plast 0
md_pfreq 0.002
md_tfreq 0.01
md_tchain 3
md_pchain 3

结构和初始速度与NVE完整STRU/KPT一致。cal_stress 1保证计算应力；md_pfirst/md_plast单位kbar，1kbar=0.1GPa。md_pfreq 0.002（fs⁻¹）控制本例压浴耦合，md_tfreq控制温浴，二者不是同一参数。

2 选择晶胞自由度#

• md_pmode iso、md_pcouple xyz：各向等比例变化。
• md_pmode aniso、md_pcouple none：允许不同方向长度独立变化。
• md_pmode tri、md_pcouple none：还允许剪切形状变化。

本例三个原始完整INPUT均在附件，对应另两个输入如下：

INPUT_PARAMETERS
suffix tutorial
pseudo_dir ../../PP_ORB
ntype 1
basis_type pw
calculation md
ecutwfc 40
nbands 64
ks_solver cg
symmetry 0
device gpu
precision double
scf_thr 1e-7
scf_nmax 100
smearing_method gaussian
smearing_sigma 0.01
cal_force 1
cal_stress 1
out_chg 0
md_type npt
md_nstep 80
md_dt 0.5
md_tfirst 300
md_seed 12345
md_dumpfreq 1
md_restartfreq 40
dft_functional pz
init_vel 1
md_tlast 300
dump_force 1
dump_vel 1
dump_virial 1
md_thermostat nhc
md_pmode aniso
md_pcouple none
md_pfirst 0
md_plast 0
md_pfreq 0.002
md_tfreq 0.01
md_tchain 3
md_pchain 3

INPUT_PARAMETERS
suffix tutorial
pseudo_dir ../../PP_ORB
ntype 1
basis_type pw
calculation md
ecutwfc 40
nbands 64
ks_solver cg
symmetry 0
device gpu
precision double
scf_thr 1e-7
scf_nmax 100
smearing_method gaussian
smearing_sigma 0.01
cal_force 1
cal_stress 1
out_chg 0
md_type npt
md_nstep 80
md_dt 0.5
md_tfirst 300
md_seed 12345
md_dumpfreq 1
md_restartfreq 40
dft_functional pz
init_vel 1
md_tlast 300
dump_force 1
dump_vel 1
dump_virial 1
md_thermostat nhc
md_pmode tri
md_pcouple none
md_pfirst 0
md_plast 0
md_pfreq 0.002
md_tfreq 0.01
md_tchain 3
md_pchain 3

3 同时看压力、体积和形状#

40fs时压力已经越过0，还在振荡；不能把这一帧直接称为目标压力下的平衡体积。温度、压力和体积都要观察足够长时间，并检查压浴质量/频率、体系尺寸及电子/应力精度。

python - <<'PY'
import csv
with open('reference/analysis/conditions.csv') as f:
    for r in csv.DictReader(f):
        if r['name'].startswith('al-npt-'):
            print(r['name'],'V/A3',r['final_V_A3'],
                  'P/GPa',float(r['final_P_kbar'])*0.1)
PY

4 别漏掉压力的动能部分#

MD_dump的VIRIAL (kbar)与热力学表中的总压力不是同一个量。本例总压力由virial张量对角平均再加2K/(3V)的离子动能贡献，附件实际重算并核对了二者。NPT物理势能加动能也不包含完整压浴/温浴扩展能量，不能用它的漂移替代NVE判据。

变胞会带来仿射位移；直接拿实验室坐标MSD拟合扩散还会混入晶胞伸缩。配套分析对NPT仅提供明确标注的参考晶胞映射位移，不将其当通用扩散系数。

参阅3.10.1压浴源码。