LAMMPS代算声子谱：分子动力学声子谱计算方法

声子谱的物理意义

声子谱（Phonon Spectrum）是描述晶体中晶格振动频率与波矢关系的物理量，是分析材料热力学稳定性、热膨胀、热导率和超导性质的重要基础数据。声子谱的计算方法主要分为两类：基于力常数的直接法（DFPT/有限位移法）和基于分子动力学轨迹的间接法。

LAMMPS作为经典的分子动力学软件，可以通过动力学矩阵方法或速度自相关函数方法计算声子谱。与VASP/CASTEP的密度泛函微扰理论（DFPT）方法相比，LAMMPS方法的优势在于可以处理大体系、高温和含缺陷的声子性质计算。

LAMMPS计算声子谱的方法

方法一：动力学矩阵法（PHONON包）

LAMMPS的USER-PHONON包提供了基于有限位移法的声子谱计算功能。其原理是对每个原子施加微小位移，计算由此产生的力响应，构建动力学矩阵，然后对角化得到声子频率。

计算步骤：

1. 在平衡结构上对每个原子施加微小位移
2. 计算位移后体系中所有原子上的力
3. 构建动力学矩阵D(q)： D(q) = (1/√(m_i * m_j)) * Σ_R Φ(R) * exp(iq·R)
4. 对角化D(q)得到声子频率

安装USER-PHONON包：

bash

复制

cd lammps/src
make yes-USER-PHONON
make mpi

方法二：速度自相关函数法（VACF）

通过分子动力学模拟轨迹计算速度自相关函数，然后通过傅里叶变换得到声子态密度：

1. 运行NVE系综MD模拟，记录原子速度
2. 计算速度自相关函数： C(t) = <v(0)·v(t)>
3. 傅里叶变换得到振动态密度（VDOS）： g(ω) = ∫C(t)exp(-iωt)dt

这种方法可以得到声子态密度，但无法直接获得色散关系（ω vs q）。

方法三：结合Phonopy使用

最常用的方案是用LAMMPS生成力数据，配合Phonopy后处理：

1. 用LAMMPS计算平衡结构和力常数
2. 导出力常数矩阵
3. 用Phonopy计算声子谱和热力学性质

LAMMPS + Phonopy完整计算流程

1. 准备力场和结构

data文件（以Si为例）：

8 atoms
2 atom types

0.0 5.43 xlo xhi
0.0 5.43 ylo yhi
0.0 5.43 zlo zhi

Masses

1 28.0855
2 28.0855

Atoms

1 1 0.0 0.0 0.0
2 1 0.0 2.715 2.715
...

2. LAMMPS输入脚本

平衡结构优化（in.minimize）：

units metal
dimension 3
boundary p p p
atom_style atomic

read_data Si.data

# Tersoff势
pair_style tersoff
pair_coeff * * Si.tersoff Si

# 能量最小化
thermo 100
thermo_style custom step pe pxx pyy pzz
min_style cg
minimize 1e-10 1e-10 1000 10000

# 写出优化后的结构
write_data Si_optimized.data

力常数计算（in.phonon）：

units metal
dimension 3
boundary p p p
atom_style atomic

read_data Si_optimized.data

pair_style tersoff
pair_coeff * * Si.tersoff Si

# PHONON计算设置
fix 1 all phonon 20 5000 0.4 Si_map.in Si_binary
run 100000

# 输出声子数据
fix 1 all phonon 20 5000 0.4 Si_map.in Si_binary

参数说明：

• 20：每20步输出一次
• 5000：总共输出5000个数据点
• 0.4：q点间距（Å⁻¹）
• Si_map.in：原子映射文件
• Si_binary：输出文件前缀

3. q点路径设置

Si_map.in文件格式：

# 原子映射文件
# ngroup nchain
2 1
# 组1：原子1
1
0 0 0
# 组2：原子2
1
1 0 0

高对称路径： 对于面心立方（FCC）结构，标准高对称路径为： Γ(0,0,0) → X(0.5,0,0) → W(0.5,0.25,0) → L(0.5,0.5,0.5) → Γ(0,0,0) → K(0.375,0.375,0.75)

4. 后处理和可视化

LAMMPS PHONON包输出的二进制文件需要后处理转换为可读格式：

bash

复制

# 使用phana后处理程序
phana -r Si_binary
# 生成声子谱数据
phana -p Si_map.in 20 0.4 0 0 0 0.5 0 0 100

输出数据可用Gnuplot、Origin或Python Matplotlib绘制声子谱图。

力场选择对声子谱的影响

力场的质量直接决定声子谱计算的准确性：

材料类型	推荐力场	说明
共价晶体(Si,Ge)	Tersoff	精确描述共价键和角度
金属体系	EAM/MEAM	描述金属键和电子效应
离子晶体	Buckingham/Born-Mayer	描述长程库仑相互作用
分子晶体	OPLS/CHARMM	描述分子间范德华力
二维材料	SW/Rebo/LCBOBII	需要验证面内振动
含氢体系	ReaxFF/AIREBO	描述化学反应和键断裂

力场验证建议：

• 先计算平衡晶格常数，与实验对比
• 检查弹性常数是否合理
• 与DFT声子谱（如果有）对比

声子谱分析的物理量

1. 动力学稳定性判断

声子谱无虚频 → 材料动力学稳定 声子谱有虚频（Γ点附近） → 结构不稳定，可能发生相变

虚频的原因：

• 结构优化不充分
• 力场参数不当
• 真实的结构不稳定性（相变前兆）

2. 热力学性质计算

从声子谱可以计算：

• 零点振动能：E_ZP = Σ (1/2)ℏω
• 自由能：F = E₀ + Σ [1/2ℏω + kT·ln(1-exp(-ℏω/kT))]
• 热容：Cv = Σ k(ℏω/kT)² exp(ℏω/kT)/(exp(ℏω/kT)-1)²
• 熵：S = Σ [ℏω/(kT(e^(ℏω/kT)-1)) – ln(1-e^(-ℏω/kT))]

3. 声子态密度分析

声子态密度（PDOS）反映不同频率的声子模式数量：

• 低频区（<5 THz）：声学声子，对应整体平移和旋转
• 中频区（5-15 THz）：光学声子，对应原子相对振动
• 高频区（>15 THz）：含氢体系的C-H/N-H伸缩振动

分波声子态密度（Projected PDOS）： 可以分析每个原子对声子模式的贡献，识别特定原子的振动特征。

常见问题与解决方案

问题1：声子谱出现大量虚频

原因和解决方案：

1. 结构优化不充分：使用更严格的收敛标准重新优化

   minimize 1e-12 1e-12 5000 50000
   

2. 力场不适合该体系：更换力场或重新拟合参数
3. 超胞太小：增大超胞到至少2×2×2
4. 对称性破缺：检查结构对称性是否正确处理

问题2：声子频率与实验偏差大

• 检查力场参数是否适用于该材料
• 考虑使用DFT方法（VASP+Phonopy）替代
• 检查温度效应（实验通常在室温，MD可模拟有限温度）
• 检查非谐效应（高阶力常数）

问题3：计算效率低

LAMMPS声子谱计算的计算量与原子数的平方成正比：

优化策略：

• 使用对称性减少需要计算的位移数量
• 并行计算（MPI加速）
• 先用小超胞测试，确认参数合理后再用大超胞
• 考虑使用MLIP（机器学习势）替代经验力场

问题4：有限温度声子谱

LAMMPS可以计算有限温度下的声子谱，反映温度效应：

# 先在NPT下热平衡
fix 1 all npt temp 300 300 0.1 iso 0 0 1
run 100000
unfix 1

# NVE下采集数据
fix 1 all nve
fix 2 all phonon 20 50000 0.4 map.in binary_300K
run 200000

温度效应包括：

• 热膨胀：晶格常数增大，声子频率红移
• 非谐效应：声子线宽展宽
• 声子软化：某些模式频率随温度升高而降低

力场参数标定方法

如果现有力场不适合声子谱计算，可以进行参数标定：

1. 收集DFT参考数据：声子频率、弹性常数、晶格常数
2. 选择力场函数形式：Tersoff、SW、MEAM等
3. 参数优化：使用GULP或potfit程序拟合DFT数据
4. 验证：对比拟合后力场计算的声子谱与DFT结果

实操案例：石墨烯声子谱计算

以石墨烯为例：

1. 结构：4×4×1超胞，32个C原子，真空层15 Å
2. 力场：Tersoff势（ optimized for graphene）
3. 计算：
   • 先做结构优化（最小化到1e-10）
   • PHONON计算：每10步输出，共10000步
4. 结果分析：
   • Γ点：3个声学声子（ZA, TA, LA）+ 3个光学声子（ZO, TO, LO）
   • K点：G带频率约1350 cm⁻¹（拉曼活性）
   • Γ点LO声子频率约1580 cm⁻¹（G带）
   • ZA声子在K附近有Kohn反常

总结

LAMMPS计算声子谱是分子动力学方法在晶格振动分析中的重要应用。通过合理选择力场、设置计算参数和正确分析结果，可以获得材料的声子色散关系、态密度和热力学性质。对于需要考虑温度效应、缺陷或大体系声子性质的计算，LAMMPS方法具有独特优势。

我们提供专业的LAMMPS声子谱计算服务，涵盖各类晶体材料、二维材料和含缺陷体系的声子性质计算，支持从力场选择到结果分析的全流程外包。