ASE晶胞与原子联合优化：外压和应力

只用普通原子优化器，通常不会自动调整晶胞。本篇使用ASE3.29.0的FrechetCellFilter，让原子位置、晶胞长度和形状一起变化，并分别在0与2 GPa外压下计算同一四原子Cu模型。

这里使用EMT经验势，不是DFT平衡晶格常数或实验高压铜数据。输入从相同轻微变形、剪切和原子扰动出发，目的是学习完整操作与单位。

1. 下载后先看全部实际曲线#

下载完整代码、102帧原生ASE轨迹、全应力与画图包：

python3 -m zipfile -e ase-dimensionality-cell-pressure-wave309-*.zip .
cd ase-dimensionality-cell-pressure-wave309
bash run.sh
cat results/cell-optimization.csv

这一步只做后处理，不重新优化。完整102个保存帧的4原子坐标、晶胞、能量、原子力与3×3应力均保留，包含零压分支的后续精化起点。

2. 准备相同初态#

import numpy as np
from ase.build import bulk
from ase.calculators.emt import EMT

initial = bulk('Cu', 'fcc', a=3.6, cubic=True)
initial.set_cell(initial.cell.array @ np.array([
    [1.05, 0.03, 0],
    [0, 0.97, 0.02],
    [0, 0, 1.02],
]), scale_atoms=True)
initial.positions += np.random.default_rng(20261109).normal(
    0, 0.02, initial.positions.shape
)

scale_atoms=True随晶胞变形缩放原子位置，然后加固定随机种子的扰动。两种外压都复制这个相同初态，不把零压已经优化的结构悄悄改成另一组起点。

3. 让原子与晶胞一起变化#

from ase.filters import FrechetCellFilter
from ase.optimize import FIRE
from ase.units import GPa

atoms = initial.copy()
atoms.calc = EMT()
pressure = 2.0 * GPa
cell_filter = FrechetCellFilter(atoms, scalar_pressure=pressure)
opt = FIRE(cell_filter, logfile='optimization.log',
           trajectory='optimization.traj')
opt.run(fmax=0.005, steps=1500)

这里scalar_pressure原生单位为eV/ų，因此用2.0 * GPa转换，不是直接写2。联合优化器还把晶胞变量的梯度与原子力一起处理；fmax是该优化问题的停止尺度，不能直接把它改名为GPa应力精度。ASE Filters说明

若只想改变原子而固定晶胞，直接对atoms优化；若只变晶胞而固定分数坐标，使用相应StrainFilter。要限制某个方向或体积，需明确选择自由度，不能仅删一列最终坐标。

4. 新计算通过Slurm提交#

资料包根目录包含完整calculate.py及模板：

cat run.slurm
sbatch run.slurm
squeue -u "$USER"

当前模板申请P100节点1张卡、1个CPU核；ASE/EMT实际在CPU运行，不据此声称GPU加速。没有外部DFT或POTCAR前置条件。已有native/时脚本停止；新数据写native/、my-fields/、my-results/和my-plots/。

一键新计算直接采用0.005优化器停止尺度。参考零压分支先采用0.02，再从已有末态精化，因此新轨迹的保存帧数可以不同；分析器处理所有实际帧，不强迫凑成102帧。

5. 原子力和应力分别检查#

import numpy as np
from ase.units import GPa

forces = atoms.get_forces()
sigma = atoms.get_stress(voigt=False)
print('Max atomic force:', np.linalg.norm(forces, axis=1).max())
print('Stress / GPa:', sigma / GPa)
print('Pressure residual / GPa:',
      np.max(np.abs(sigma + pressure * np.eye(3))) / GPa)

这里压缩对应负的应力，外部静水压P的目标为σ≈−P I，所以检查σ+P I。完整应力张量还包含剪切分量，不能只看平均压力而忽略形变方向。

本例零压分支在0.02停止尺度时，原子力已经小于0.02 eV/Å，但一个应力分量约为−0.215 GPa。希望最大应力残余小于0.2 GPa时，还应从已有结构继续优化；本次把停止尺度收紧到0.005，13步精化后残余约0.0592 GPa。这里没有重新从最初结构算一遍，也没有把原子力的阈值直接当成应力阈值。

6. 本次最终数值#

2 GPa下的体积小于零压，符合这个模型的压缩响应。这里是有限参数的小例，不把最后几位小数当作材料精度。

7. 有外压时看焓#

在每个固定外压条件下：

H = E + P V

图中每条焓曲线减去自己相同外压下的初始值。不能把0与2 GPa的焓绝对值放在一起，再据此判断哪个外压“更准确”；它们是不同外部条件。

python read_native_tools.py native-reference --output my-read
python analyze.py --data my-read --output my-analysis
python plot.py --data my-analysis --output my-plots
cat my-analysis/DATA.json

读取器从原生.traj重新提取每个保存帧的完整原子和晶胞信息。自己的新计算把native-reference换成native即可。原始几何的维度识别与连通分量介绍见0D到3D结构分析。