Phonopy 跨软件结构转换

从VASP结构转到QE、ABACUS或LAMMPS时，先转换几何，再补相应软件的计算参数，最后求力。本文用Si两原子原胞和一个双元素三斜结构实际完成14条格式路线；这是结构文件操作，尚未运行这些软件的电子计算或MD。

1. 下载与准备#

下载结构转换和读回实例包，在已分配的CPU资源内操作：

unzip phonopy-interfaces-wave149-*.zip -d phonopy-interfaces
cd phonopy-interfaces
bash run.sh list

模块为 phonopy/4.8.3：Python3.11.16、NumPy2.3.5、SciPy1.16.3、spglib2.7.0；Phonopy C/OpenMP部分由GCC13.3编译，另有phonors0.6.0 Rust后端。CP2K结构解析采用cp2k-input-tools0.9.1的固定393e87bb版本。每次转换、分析使用已分配范围内1个CPU，不需要GPU。

一份新解压目录运行 bash run-all.sh，依次转换28个例子、读回为VASP格式、比较几何并画图；结果在 analysis/roundtrip-errors.csv、analysis/CHECKS.json 和 figures/。另一份新解压目录只提取包内真实结果：

bash analyze.sh
ls analysis-reference figures-reference

同一输出目录不会重复覆盖，继续试另一条路线或重新解压到新目录即可。包内 run.slurm 是CPU分区56核申请模板，实际工具只使用其中1核；按可用资源提交。大型超胞的计算另外申请资源。

2. 原始结构是什么#

Si使用晶格常数5.43Å的金刚石原胞，完整输入为：

generated by phonopy
   1.0
     0.0000000000000000    2.7149999999999999    2.7149999999999999
     2.7149999999999999    0.0000000000000000    2.7149999999999999
     2.7149999999999999    2.7149999999999999    0.0000000000000000
Si
   2
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000
  0.2500000000000000  0.2500000000000000  0.2500000000000000

第二个NaCl两原子例子刻意选用一般三斜晶胞和非特殊坐标，用于检查两种元素及坐标是否保持；它没有结构优化，不能当作岩盐相的平衡晶胞。

Two-species skew-cell geometry example; not a relaxed material
1.0
5.2 0.0 0.0
0.7 5.7 0.0
0.2 0.4 6.1
Na Cl
1 1
Direct
0.08 0.17 0.29
0.59 0.63 0.74

3. 先选择自己要的格式#

下表中的命令在包的根目录运行，每条生成自己的 work/名称/structure 和 roundtrip.vasp。同一格式另有对应 nacl-名称 例子。

这些是格式接口示例，不表示集群已安装表中每一种电子结构程序或提供其授权。exciting读回本例时按标准元素文件名识别元素；下载包不含species电子数据。

4. QE：原生结构卡片之后，还要补namelist#

需要看原生命令时，可在新目录转换：

module load phonopy/4.8.3
phonopy-calc-convert --calcin vasp --calcout qe \
  -i inputs/Si.vasp -o Si-cards.in \
  --additional-info Si Si.upf
cat Si-cards.in

--additional-info填写元素和赝势文件名，不会替你下载或验证赝势。此版本写出的是结构卡片，第一行 ! ibrav=... 是注释；不能把它当作已经填写好的&SYSTEM。

配套 bash run.sh si-qe 会明确补最小namelist及Γ点卡片，原始结构卡片另外保存在 structure-cards-only。得到的几何骨架为：

&CONTROL
 calculation="scf"
/
&SYSTEM
 ibrav=0, nat=2, ntyp=1
/
&ELECTRONS
/
!    ibrav = 0, nat = 2, ntyp = 1
CELL_PARAMETERS bohr
    0.0000000000000000    5.1306064620899088    5.1306064620899088
    5.1306064620899088    0.0000000000000000    5.1306064620899088
    5.1306064620899088    5.1306064620899088    0.0000000000000000
ATOMIC_SPECIES
 Si   28.08550   Si.upf
ATOMIC_POSITIONS crystal
 Si   0.0000000000000000  0.0000000000000000  0.0000000000000000
 Si   0.2500000000000000  0.2500000000000000  0.2500000000000000
K_POINTS gamma

这里仍缺 ecutwfc/ecutrho、实际 pseudo_dir、收敛参数和适合体系的k采样。K_POINTS gamma 仅为结构读回准备，不能代替收敛选择。先按QE操作教程补齐；需要真实“位移→SCF求力→力常数”流程，可接着看实际求力教程。

5. CP2K：明确提供模板，再写入结构#

bash run.sh si-cp2k
cat inputs/cp2k-structure-template.inp
cat work/si-cp2k/structure

CP2K路线通过Phonopy原生Python API读取结构模板的配置树，再替换晶胞和坐标。它是配套脚本明确提供的辅助步骤；此版本的 phonopy-calc-convert --calcout cp2k 本身不会生成配置树。核心操作为：

from phonopy.interface.calculator import (
    read_crystal_structure, write_crystal_structure, Cp2kStructureInfo,
)
from phonopy.interface.cp2k import read_cp2k
cell, _ = read_crystal_structure("inputs/Si.vasp", "vasp")
_, tree = read_cp2k("inputs/cp2k-structure-template.inp")
info = Cp2kStructureInfo(unitcell_filename="Si.inp", config_tree=tree)
write_crystal_structure("Si.inp", cell, "cp2k", info)

生成文件含 &CELL、&COORD、ENERGY_FORCE 以及力输出设置。还要按照CP2K教程添加DFT、KIND、基组/赝势及SCF控制，才能实际求力。

6. 最后得到什么#

bash run-all.sh
head -5 analysis/roundtrip-errors.csv
ls figures

全部28条路线保留了原子数量、按元素排列的顺序、分数坐标、晶胞几何和周期最短距离。本例读回晶胞边长最大差约2.20×10⁻⁹Å，周期原子间距最大差约5.84×10⁻¹⁰Å；这反映文本写入精度，不代表电子计算精度。

长度、原子顺序和LAMMPS晶胞旋转的检查方法见单位与晶胞核对。格式及接口约定参考Phonopy官方接口说明；本文命令对应集群安装的4.8.3。