Calorine结构优化与状态方程

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

1 原子位置优化与晶胞优化分开#

固定晶胞只允许原子移动，适合比较同一体积或在给定晶格下松弛；晶胞优化还改变晶格和体积，适合寻找指定外压下的平衡结构。二维模型要保留真空约束，不能照搬本例三维各向晶胞优化。

from ase.build import bulk
from calorine.calculators import CPUNEP
from calorine.tools import relax_structure
from ase.io import write
atoms = bulk('Si', 'diamond', a=5.431, cubic=True)
atoms.rattle(0.03, seed=28)
atoms.calc = CPUNEP('si-nep.txt')
relax_structure(atoms, constant_cell=True, fmax=0.02,
                steps=200, logfile='relax.log')
write('relaxed.xyz', atoms)
# 再求零外压晶胞；力和晶胞收敛由过滤器共同控制
relax_structure(atoms, constant_cell=False, scalar_pressure=0.0,
                fmax=1e-5, steps=200, logfile='cell.log')
print(atoms.get_volume(), atoms.get_stress())

fmax的原子力单位为eV/Å；本例固定晶胞优化使用0.02阈值，实际最大力从1.23564降到0.007101 eV/Å。外压scalar_pressure的单位为GPa。有限步数退出后仍要检查实际力和应力，不能仅看脚本返回。

2 七个体积拟合状态方程#

python cpu_workflows.py
cat eos.csv
cat CPU-RESULTS.json
python plot.py

完整脚本使用0.94、0.96、0.98、1.00、1.02、1.04、1.06倍参考体积。线性晶胞因子是体积倍率的立方根：

atoms = bulk('Si', 'diamond', a=5.431, cubic=True)
volumes, energies = [], []
for ratio in [0.94, 0.96, 0.98, 1, 1.02, 1.04, 1.06]:
    sample = atoms.copy()
    sample.set_cell(atoms.cell * ratio**(1/3), scale_atoms=True)
    sample.calc = CPUNEP('si-nep.txt')
    volumes.append(sample.get_volume())
    energies.append(sample.get_potential_energy())
from ase.eos import EquationOfState
from ase.units import GPa
v0, e0, b0 = EquationOfState(volumes, energies,
                            eos='birchmurnaghan').fit()
print(v0, e0, b0 / GPa)

本次Si8拟合平衡体积 162.358699 ų，体模量 92.1846 GPa。这是该势与该采样范围的结果；不是新的DFT，也不自动等于实验室温体模量。直接晶胞优化体积162.353548 ų，与七点拟合存在有限采样差异。

3 换成自己的体系#

先用ASE读取自己的POSCAR或extxyz，确认周期方向、原子类型和模型元素完全一致。优化前后分别保存结构；EOS必须让最低点落在采样区间中间，若最低能量在端点，应扩大或平移体积范围。合金、低对称结构还需在每个固定体积松弛内部坐标，不能直接把本例理想Si流程当作所有体系的最终做法。