GPUMD 替位合金建模与结构优化

本例使用鸢算科技集群的 GPUMD 5.9.1 P100 构建（GCC 13.3、CUDA 12.9），用 ASE 3.29、NumPy 和 Matplotlib 建模、提取和画图。每个原生任务只用一张 P100。所用数据是经典势实际计算结果。

下载本篇配套输入、真实结果与脚本，解压进入目录。先查看已有结果：

bash run.sh

analysis/RESULTS.json 为结果汇总，CSV 可继续作图，figures/ 为高清图片。这个步骤不需要 GPU。要重新计算，另解压一份，进入新目录：

sbatch compute.slurm

脚本按小体系静态计算、短热接触的顺序运行，申请一张 P100、一个 CPU 核心。inputs/ 是输入，results/ 是配套参考输出，work/ 是你新运行的结果。不要把新输出追加到原参考目录。完成后自动提取 CSV 和图片。

1. 确定原子数、成分和随机种子#

本例从金刚石 Si 常规胞扩展 4×4×4，共 512 原子。按固定随机种子 160031 选择 256 个位置替换为 Ge，得到 Si0.5Ge0.5 替位模型。使用官方 Si/Ge 双元素 Tersoff 1989 势，元素必须与势文件声明一致。

cat inputs/sige-alloy/CASE.json
head -n 6 inputs/sige-alloy/model.xyz
head -n 3 inputs/sige-alloy/potential.txt

建模脚本 prepare_alloy.py 的核心操作：

import numpy as np
from ase.build import bulk
a = bulk('Si', 'diamond', a=5.5, cubic=True).repeat((4, 4, 4))
ids = np.random.default_rng(160031).choice(len(a), 256, replace=False)
for i in ids:
    a[i].symbol = 'Ge'

完整脚本还写入零初速度和元素分组：方法 0 的组 0 为 Si、组 1 为 Ge。这里的组 0 只是元素标签，没有 fix 0 时不会自动固定。不同案例中标签意义由该案例自己定义，不要照搬热接触的固定组约定。

2. 同时优化原子位置与静水晶胞#

cat inputs/sige-alloy/run.in

potential potential.txt
minimize fire 0.0001 6000 1 1
ensemble nve
time_step 0
dump_thermo 1
dump_xyz 1 snapshot.xyz precision double mass velocity force potential virial group_labels
run 1

该设置让原子位置和晶胞体积一起变化，保持立方形状。若研究剪切应力或非立方晶胞，不要把这里的静水形变当作所有晶格自由度均已完全优化。

3. 查看结构、能量和力#

head -n 5 work/sige-alloy/snapshot.xyz
cat work/sige-alloy/thermo.out

本次实际得到：

结构图为原始与优化后的薄片投影，颜色区分 Si、Ge；完整三维结构保留在 XYZ 中：

若用 OVITO 查看，可打开 snapshot.xyz，按元素着色并添加切片；XYZ 中的分组标签同样保留。进一步的轨迹与绘图工具见GPUMDkit。

4. 怎样用于自己的合金#

修改元素比例、初始晶格和随机种子，先确认势函数确实覆盖这组元素及其交叉作用，再优化结构。一个随机小胞只是一份指定构型，不等于随机合金已统计收敛；本例没有给出混合自由能、相图或平衡组分。要比较构型能量，应保持相同原子数/成分、能量零点和优化约束，并增加独立构型或 SQS 等建模步骤。

与这个“原子混合的替位模型”不同，Si/Ge 热接触是两块材料沿 x 拼接的界面模型；双层层间势则演示层间作用组合。

依据：GPUMD 模型文件与分组、最小化。

用ASE继续分析化学近邻#

ASE RDF与化学分波沿用本篇512原子真实结构，比较SiSi/SiGe/GeGe分波、成分加权总RDF和第一近邻配位积分。