Calorine弹性常数与声子

下载本专题结构、势文件、真实参考数据和全部Python脚本。以下命令在解压目录执行；CPU计算使用已分配的CPU资源，GPU通过Slurm分配。原始结果放在 reference/，复跑结果写入当前目录。

1 先准备平衡结构#

弹性和声子使用包内先做零外压晶胞优化的Si结构；平衡Si8体积162.353548 ų，最大残余应力约3.52E-10 eV/ų。将预应力结构直接作为平衡结构，会改变弹性解释。

python cpu_workflows.py
ls equilibrium.xyz elastic.npz phonopy.yaml
python plot.py

2 夹持离子与松弛离子弹性#

from ase.io import read
from calorine.calculators import CPUNEP
from calorine.tools import get_elastic_stiffness_tensor
atoms = read('equilibrium.xyz')
atoms.calc = CPUNEP('si-nep.txt')
C_fixed = get_elastic_stiffness_tensor(atoms, clamped=True,
                                       epsilon=1e-3)
C_relax = get_elastic_stiffness_tensor(atoms, clamped=False,
                                       epsilon=1e-3)
print(C_fixed)
print(C_relax)

返回值已经是GPa，不能再除一次ASE的GPa。clamped=True指应变后内部原子保持分数坐标；False允许每个变形结构的内部坐标松弛，晶胞仍由该次应变确定。

剪切响应的差别来自内部原子弛豫。本例应变幅度1E-3；研究精度需至少比较几组幅度和更严格的内部松弛。

3 用Calorine生成Phonopy力常数#

import numpy as np
from ase.build import bulk
from calorine.calculators import CPUNEP
from calorine.tools import get_force_constants
atoms = bulk('Si', 'diamond', a=162.35354772209914**(1/3))
ph = get_force_constants(
    atoms, CPUNEP('si-nep.txt'), np.diag([3,3,3]),
    kwargs_generate_displacements={'distance':0.01})
ph.symmetrize_force_constants()
ph.save('phonopy.yaml', settings={'force_constants':True})
result = ph.run_qpoints([[0,0,0],[0.5,0,0.5]])
print(ph.qpoints.frequencies)

3×3×3原胞超胞共54个Si，晶体对称性将本例独立位移压缩到1个。实际Γ点光学频率14.326164 THz；声学零模约2E-7 THz数量级。图中沿Γ–X–L–Γ读取真实力常数，未用手画曲线替代计算。

4 替换结构时检查什么#

确认原胞、超胞矩阵、周期边界与模型一致；低对称结构通常需要更多位移。实质负频应检查结构是否稳定、超胞尺寸及位移幅度，不能统一删除或取绝对值。NEP声子反映模型能量面的二阶导数，必须评估模型在目标结构与形变附近的可靠性。