ALAMODE 准谐结构优化：ZnO 晶胞与原子位置

QHA 可以回答一个常见问题：温度改变以后，晶胞的尺寸与形状、晶胞内部的原子位置怎样变化？这与固定晶胞的声子热力学不同，也与自洽声子计算的频率重整化不同。本例使用 ZnO，实际同时优化晶胞与原子位置。

先看流程#

下载资料并解压 → 准备完整力常数和应变耦合 → 修改温度及优化方式 → 运行 ANPHON → 检查每个温度的结构收敛 → 提取晶格、体积、原子位移 → 画图。

本例使用 ALAMODE 官方 ZnO 示例中的 DFT 力常数及配套弹性/应变数据；在鸢算科技集群实际重新运行 QHA。本次没有重新做 ZnO 的 DFT。谐波力常数来自 4×4×2 超胞，高阶力常数来自 3×3×2 超胞，原胞包含四个原子。

资料包：ZnO QHA 输入与绘图脚本。解压后进入目录：

unzip alamode-qha-wave191.zip -d zno-qha
cd zno-qha
module load alamode/1.5.0
python prepare.py

prepare.py 会解压官方 XML，建立三种方案的独立目录。完整计算和画图也可在已申请的 CPU 计算资源中直接运行：

bash run.sh

这里采用 ALAMODE 1.5.0，oneAPI 2024.2 编译器、Intel MPI 与 MKL 顺序数学库。正式计算应将命令放入本站 ALAMODE CPU 作业模板，经 Slurm 提交；本例单独限制为一个线程，便于快速操作演示。

1. 先检查需要哪些文件#

ls data/*.xml
ls data/strain_IFC

需要谐波和高阶力常数，以及二/三阶弹性常数、应变引起的力、应变下的谐波力常数。strain_IFC/ 不是可有可无的附属文件；换自己的材料时，不能只改 KD = Zn O 就直接使用本例的 XML 和耦合数据。

elastic_constants.in 使用上游约定的 Ry 单位及晶胞体积因子；不要直接把以 GPa 输出的弹性矩阵原样填入。应变 XML、参考晶胞、原子顺序和力常数必须属于同一套数据。

2. 输入文件里改什么#

cd cases/scheme-0
cat input.in

本例实际输入如下：

&general
  PREFIX = ZnO_qha
  MODE = QHA
  EMIN = 0
  EMAX = 700
  DELTA_E = 0.5
  NKD = 2
  FCSXML = ZnO332_500K_cutoff1208_nbody233.xml
  FC2XML = ZnO442_harmonic.xml
  KD = Zn O

  TMIN = 0
  TMAX = 600
  DT = 300
/

&qha
  KMESH_INTERPOLATE = 4 4 2
  KMESH_QHA = 4 4 2

  RELAX_STR = 2
  QHA_SCHEME = 0
/ 
&relax
  RELAX_ALGO = 2
  MAX_STR_ITER = 1000
  COORD_CONV_TOL = 1.0e-5
  MIXBETA_COORD = 0.5
  CELL_CONV_TOL = 1.0e-5
  MIXBETA_CELL = 0.7

  SET_INIT_STR = 2
  ADD_HESS_DIAG = 0.0

  RENORM_3TO2ND = 2
  RENORM_2TO1ST = 2
  RENORM_34TO1ST = 0

  STRAIN_IFC_DIR = strain_IFC
/

&cell
  1.88972612462577
  3.235859326375770 0.000000000000000 0.000000000000000
  -1.617929663187880 2.802336379714220 0.000000000000000
  0.000000000000000 0.000000000000000 5.224712025937350
/

&displace
1
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
/

&strain
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
/

&kpoint
  2 
  8 8 8
/

MODE = QHA 选择准谐结构优化；RELAX_STR = 2 同时优化内部原子位置与晶胞，RELAX_STR = 1 只优化内部坐标、固定晶胞。QHA 不使用 RELAX_STR = 0 来代替普通固定结构声子计算。

KMESH_QHA 是结构优化用的声子网格，KMESH_INTERPOLATE 是插值网格，每个方向前者应为后者的整数倍。这里都用 4×4×2，&kpoint 中的 8×8×8 用于后处理；三者作用不同。它们是教学参数，正式材料计算要检查网格收敛。

TMIN = 0、TMAX = 600、DT = 300 给出 0、300、600 K。结构优化同时检查 COORD_CONV_TOL 与 CELL_CONV_TOL；其中原子位移更新量使用 Bohr，晶胞更新量是无量纲应变。它们不是 VASP 的电子或原子力收敛阈值。

DELTA_E = 0.5 是输出态密度的频率采样间隔，单位 cm⁻¹。它不改变温度步长，也不是能量截断。输入里还给出 0–700 cm⁻¹；1.5.0 的该后处理流程会按频率自动更新输出范围，读图时以实际文件第一列为准。

3. 运行并检查每个温度#

在作业申请到的 CPU 资源中运行：

anphon input.in > anphon.log
grep 'Structural optimization converged' anphon.log
tail -20 anphon.log

三个温度应各有一条结构优化收敛信息，再检查 .umn_tensor、.atom_disp 与 .qha_thermo 是否包含三个温度。只看退出码、只看到最后一行 Job finished，不足以判断所有温度都已优化好。

4. 提取晶格、体积与内部位移#

cat ZnO_qha.umn_tensor
cat ZnO_qha.atom_disp
cat ZnO_qha.qha_thermo
cd ../..
/public/toolkit/miniconda3/envs/ase-3.29.0/bin/python analyze.py

最后一条命令在三种方案都算完后执行；一键 run.sh 会顺序完成这些步骤。提取表为 analysis/qha-summary.csv，三张图片在 figures/。

.umn_tensor 第一列为温度，后九列依次是 3×3 位移梯度张量。以输入参考晶胞为 (A_{ref})，晶格矩阵按 (A(T)=A_{ref}(I+u)^T) 计算，体积为 (V(T)=V_{ref}\det(I+u))。本例 u_xx ≈ u_yy，说明六方基面方向保持等价。

特别注意：输出应变相对于输入的参考结构，不是自动相对于 0 K。要看热变化，需计算 a(T)/a(0 K)-1 或 V(T)/V(0 K)-1。0 K 仍有零点振动贡献，不能强行把文件里 0 K 的应变改成零。

.atom_disp 中位移为 Bohr，乘 0.529177210903 转为 Å。它与晶胞应变是两组不同结果：晶胞增大并不意味着内部原子分数坐标完全不变。

本次 Full QHA 的实际结果：

5. 自由能和热容怎样读#

.qha_thermo 给出温度、声子热容 Cv、振动自由能、重整静态势能 Phi0 和总自由能。这里的声子 Cv 对应每个温度下的优化结构，单位为每四原子原胞的 kB；它不是自动求得的定压热容 Cp。总自由能包含 Fvib + Phi0，也不能只拿 Fvib 一列比较结构的总稳定性。

本例只有三个温度点，适合学会运行、读表和画图。要计算热膨胀系数，应采用更密的温度点、更严格的结构阈值，并检查声子网格；不能用这三个点直接宣称得到了收敛的热膨胀曲线。

更多方法对照见QHA、ZSISA 与 v-ZSISA。数据来源与方法说明见 ALAMODE 官方 ZnO 教程。