ASE 计算化学建模与工作流教程

适用版本：ASE 3.28.0 · Python 3.14.4 · Windows + Miniconda · VS Code
说明：文中所有 import 路径与函数名均基于 ASE 3.28.0 真实 API。
凡不完全确定的细节，都会明确标注"以官方文档为准"。
EMT 等纯 Python 计算器示例可直接运行；VASP 相关示例为配置模板，需装有 VASP + POTCAR。

1. 简介与定位：ASE 在 DFT 工作流中的角色

ASE（Atomic Simulation Environment）不是做 DFT 计算的"引擎"，而是一套围绕计算引擎的外围工具：它帮你完成建模、提交、结果解析、后处理、可视化的全流程，底层真正的量子/力场计算由 VASP、GPAW、LAMMPS、Quantum Espresso 等计算器完成。

在你的日常 DFT 工作流（尤其是 CO₂ 还原、单原子催化、MOF、slab 模型）里，ASE 的角色通常是这样的三段式：

[建模]                [计算]                  [后处理]
ase.build 构建 slab    →  VASP / GPAW 算单点/优化  →  ase.io 读能量受力
ase.io 读写 POSCAR         ase.calculators 封装参数    批量提取、算吸附能
ase.constraints 固定原子    ase.optimize 驱动弛豫       可视化、画 DOS/势能面

ASE 能让你做三件很值的事：

用 Python 批量建模：循环生成不同晶面、不同吸附位、不同覆盖度的 slab，而不是手改 POSCAR。
统一的计算器接口：get_potential_energy() / get_forces() 这套方法对 EMT、VASP、GPAW、LAMMPS 都一样，换引擎几乎不用改脚本。
自动化后处理：read('vasprun.xml') 一次把能量、受力、轨迹都读回来，批量算吸附能/形成能。

一个贯穿本教程的最小示例（可真实运行）：

from ase.build import bulk
from ase.calculators.emt import EMT

atoms = bulk("Cu", "fcc", a=3.61)   # 构建 Cu 面心立方晶胞
atoms.calc = EMT()                  # 挂上 EMT 纯 Python 计算器（仅支持部分金属）
print("能量 =", atoms.get_potential_energy(), "eV")   # 单点能
print("受力 =", atoms.get_forces())                    # 受力，单位 eV/Å

这个 4 行的框架，就是后面所有章节的骨架。

2. 安装与版本检查、可视化依赖

2.1 安装

推荐在 Miniconda 里建一个独立环境，避免污染基础环境：

conda create -n ase python=3.14 -y
conda activate ase
pip install ase matplotlib

注意：ASE 3.28.0 是纯 Python 包，pip install ase 即可。matplotlib 用于画 PNG，pillow（通常随 matplotlib 装好）也是可视化依赖之一。

2.2 版本检查

import ase
print(ase.__version__)          # 应输出 3.28.0
print(ase.__file__)             # 安装路径，出问题时排查用

命令行工具检查：

ase info          # 输出版本、Python 版本、已装的可选依赖
ase --version

ase info 会列出哪些可选依赖（matplotlib、spglib、flask 等）可用，缺什么一目了然。

2.3 可视化依赖

ase.visualize.view（弹出 ase-gui 图形窗口）：需要 ase-gui 命令，随 ASE 自带。Windows 上应能直接弹窗。
write('x.png', atoms)（导出静态图）：需要 matplotlib。
可选高级渲染：povray（POV-Ray 渲染）、blender（Blender 渲染），不是必需，先跳过。

pip install matplotlib   # 至少要这个，否则画不了 PNG

3. Atoms 对象：核心数据结构

ASE 里一切结构都是 Atoms 对象，它内部封装了原子种类、坐标、晶胞、周期性边界条件四类信息。

3.1 手工创建一个 Atoms

from ase import Atoms

# 一个 H2 分子，放在原点附近
h2 = Atoms(
    symbols="H2",                       # 化学式；也可用 ["H", "H"]
    positions=[[0.0, 0.0, 0.0],
               [0.0, 0.0, 0.74]],       # 笛卡尔坐标，单位 Å
)
print(h2)                               # Atoms 对象可直接打印

带晶胞和周期性边界条件（PBC）的晶体：

from ase import Atoms

# 一个简单立方晶胞，边长 3 Å，PBC 三个方向全开
atoms = Atoms(
    symbols="Cu",
    positions=[[0.0, 0.0, 0.0]],
    cell=[3.0, 3.0, 3.0],               # 晶胞，单位 Å
    pbc=[True, True, True],             # 周期性边界条件
)

3.2 常用属性与方法

atoms.positions            # numpy 数组 (N, 3)，原子笛卡尔坐标（Å）
atoms.get_positions()      # 同上，等价写法
atoms.cell                 # Cell 对象，3x3 矩阵
atoms.cell.lengths()       # 晶胞三个边长 a,b,c（Å）
atoms.cell.angles()        # 晶胞三个夹角 α,β,γ（度）
atoms.pbc                  # 布尔列表，如 [True, True, True]
atoms.get_chemical_symbols()   # ["Cu", "O", ...]
atoms.get_atomic_numbers()     # [29, 8, ...]
atoms.get_masses()             # 原子质量（原子质量单位 amu）
atoms.get_volume()             # 晶胞体积（Å³）
atoms.get_number_of_atoms()    # 原子数 N
len(atoms)                     # 原子数，等价写法
atoms.get_center_of_mass()     # 质心坐标（默认按质量加权）
atoms.get_positions()

3.3 单位系统（重要，第 9 节还会细讲）

ASE 内部约定一套统一的原子单位，所有计算器返回的结果都遵循它：

物理量	ASE 单位
长度	Å（埃，1 Å = 10⁻¹⁰ m）
能量	eV
力	eV/Å
应力	eV/Å³
质量	原子质量单位（amu）
温度	K（开尔文）
时间	ASE 时间单位（约 10.18 fs）

VASP 的输出正是这套单位，所以 POSCAR（Å）、OUTCAR 里的能量（eV）几乎可以无缝对接，这是 ASE 对 VASP 用户特别友好的原因。

4. 构建结构：ase.build 的常用函数

ase.build 提供了一批一键建结构的函数。最常用的几类：

4.1 bulk —— 块体晶体

from ase.build import bulk

cu  = bulk("Cu", "fcc", a=3.61)              # 面心立方
fe  = bulk("Fe", "bcc", a=2.87)              # 体心立方
mg  = bulk("Mg", "hcp", a=3.21, c=5.21)      # 六方密堆积
si  = bulk("Si", "diamond", a=5.43)          # 金刚石结构
nacl= bulk("NaCl", "rocksalt", a=5.64)       # 岩盐结构

a 是晶格常数；hcp 等还需 c。
加 cubic=True 可把 fcc/bcc 的元胞展开成简单立方超胞（方便看晶格）：
```
cu_cubic = bulk("Cu", "fcc", a=3.61, cubic=True)
```
加 orthorhombic=True 得到正交胞。

注意：ase.build 里没有独立的 fcc() / bcc() 函数。构建块体 fcc/bcc 一律用 bulk(name, 'fcc'/'bcc', a=...)。名为 fcc111 / bcc100 的函数是表面（slab）构建器，见 4.3。

4.2 molecule —— 分子（来自内置 G2 数据库）

from ase.build import molecule

co2 = molecule("CO2")    # CO2 线性分子
h2o = molecule("H2O")
ch4 = molecule("CH4")
c2h4= molecule("C2H4")   # 乙烯，单原子催化里常见的吸附探针分子

molecule 返回一个孤立的、无 PBC 的分子（默认一个方向留了真空），坐标来自 G2 参考几何。

坑：molecule() 只负责搭结构，不给能量。你无法用 EMT 算 CO2 的能量，因为 EMT 只支持有限几种金属（详见第 7 节）。

4.3 表面（slab）——fcc111 / bcc100 / surface

构建表面最常用的是这些便捷函数，它们自动按 Miller 指数切面并加真空：

from ase.build import fcc111, fcc100, bcc100, hcp0001

# Pt(111) 表面：2x2 面内超胞，3 层原子，10 Å 真空
pt111 = fcc111("Pt", size=(2, 2, 3), a=3.92, vacuum=10.0)

# Fe(100)：2x2，4 层
fe100 = bcc100("Fe", size=(2, 2, 4), a=2.87, vacuum=10.0)

# 其他常见：fcc110 / fcc211 / bcc110 / bcc111 / hcp0001 等

size=(nx, ny, nz)：前两个是面内重复次数，第三个是层数。
a：晶格常数（一般用你 DFT 优化后的实验值/文献值）。
vacuum：在 z 方向加的真空厚度（Å）。

通用切面函数 surface()：当便捷函数覆盖不了你的晶格（比如某些 MOF / 氧化物 / 自定义晶体）时，用更通用的 surface：

from ase.build import surface

slab = surface("Pt", (1, 1, 1), layers=3, vacuum=10.0)   # 1x1 元胞
# 想做大超胞就手动 repeat
slab = slab.repeat((2, 2, 1))                            # 2x2

注意 surface() 生成的默认是 1×1 元胞；便捷函数 fcc111 通过 size 一步到位，通常更方便。

4.4 吸附：add_adsorbate

from ase.build import fcc111, molecule, add_adsorbate

slab = fcc111("Pt", size=(2, 2, 3), a=3.92, vacuum=10.0)
co2  = molecule("CO2")

# 把 CO2 放在 slab 顶位（ontop），C 原子距表面 2.0 Å
add_adsorbate(slab, co2, height=2.0, position="ontop")

position 可取 "ontop" / "bridge" / "fcc" / "hcp" 等，针对面内高对称位。
height：吸附原子（默认 mol_index=0，即分子第一个原子）到表面的垂直距离（Å）。
offset=(dx, dy)：在高对称位上再做面内微调。

4.5 直接手写 Atoms 构建（最灵活）

当你需要精确控制每个原子时，手写：

from ase import Atoms

# 手写一个 3 层 2x2 Pt(111) 的最小替代：一个 fcc 位吸附的 CO（示意）
atoms = Atoms(
    symbols=["Pt", "Pt", "C", "O"],
    positions=[
        [0.0, 0.0, 0.0],        # Pt
        [1.385, 0.0, 0.0],      # Pt（示意坐标）
        [0.693, 0.4, 2.0],      # C
        [0.693, 0.4, 3.15],     # O
    ],
    cell=[11.0, 11.0, 20.0],
    pbc=[True, True, True],
)

4.6 读写后构建（从已有结构出发）

from ase.io import read

atoms = read("POSCAR")              # 读入 VASP 结构
supercell = atoms.repeat((2, 2, 1)) # 做超胞
supercell.center(vacuum=10.0, axis=2)  # 沿 z 居中并加真空

center(vacuum=..., axis=...) 是"读进来后补真空"最常用的手段。

5. 读写格式：POSCAR / CONTCAR / CIF / XYZ / traj

ase.io 的 read 和 write 根据文件扩展名自动识别格式，这是它最好用的地方。

from ase.io import read, write

atoms = read("POSCAR")            # 自动识别为 VASP POSCAR 格式
write("CONTCAR", atoms)           # 写成 VASP 格式
write("out.cif", atoms)           # 写成 CIF
write("out.xyz", atoms)           # 写成 XYZ
write("out.traj", atoms)          # 写成 ASE 轨迹文件

5.1 VASP 的 POSCAR 读写（重点）

读：ASE 会自动处理 VASP4 / VASP5 两种格式，以及"笛卡尔 / 分数坐标（Direct）"两种坐标，无需手动指定：

atoms = read("POSCAR")
print(atoms.cell)          # 晶格，Å
print(atoms.positions)     # 始终是笛卡尔坐标（Å），ASE 已帮你转好
print(atoms.pbc)           # 通常 [True, True, True]

常用参数：

atoms = read("POSCAR", sort=True)      # 按元素排序（和 VASP 内部原子顺序对齐）
atoms = read("POSCAR", index=-1)       # 多帧文件里取最后一帧（POSCAR 本身是单帧，此处仅作通用写法演示）

写：写入 POSCAR 时有两个关键参数：

write("POSCAR", atoms)                    # 默认：笛卡尔坐标，VASP4 风格
write("POSCAR", atoms, direct=True)       # 用分数坐标（Direct）写
write("POSCAR", atoms, vasp5=True)        # VASP5 风格（带元素名那行）
write("POSCAR", atoms, vasp5=True, direct=True)  # 两者组合
write("POSCAR", atoms, sort=True)         # 写之前按元素排序

具体默认行为（direct 默认 True 还是 False）不同版本可能略有差异，建议显式写出 direct= 和 vasp5=，避免歧义。此点以官方文档 ase/io/vasp.py 为准。

5.2 CONTCAR 与"优化结果回读"

优化完的 CONTCAR 就是优化后结构，直接 read("CONTCAR") 即可。批量处理时（第 12 节）会用 glob 遍历目录读 CONTCAR / vasprun.xml。

5.3 CIF

atoms = read("mof.cif")
write("mof_from_ase.cif", atoms)

CIF 适合从晶体学数据库（如 CSD、Materials Project 导出的 CIF）导入结构。注意 CIF 里的对称性信息 ASE 默认会展开成 P1（即列出所有原子），这通常是你要的。

5.4 XYZ

atoms = read("cluster.xyz")
write("cluster.xyz", atoms)

XYZ 格式本身不含晶胞信息，所以周期结构写 XYZ 会丢失晶胞（PBC 也丢）。XYZ 适合分子、团簇、以及给其他软件（如 LAMMPS 前处理、VMD、Ovito）交换坐标。

5.5 traj（轨迹）与多帧读写

traj 是 ASE 原生二进制格式，能存多帧（每帧一个 Atoms，带能量受力）。常用于优化/MD 轨迹：

from ase.io import read, write
from ase.io.trajectory import Trajectory

# 读整个轨迹（返回 Atoms 列表）
images = read("opt.traj", index=":")     # 所有帧
last   = read("opt.traj", index=-1)      # 最后一帧

# 逐帧写入
traj = Trajectory("my.traj", "w")
for atoms in some_structures:
    traj.write(atoms)     # 每帧可带 energy/forces
traj.close()

优化器（第 8 节）会自动写 traj，你只需读它。

5.6 vasprun.xml（最推荐的能量读取入口）

这是 ASE 读 VASP 结果的首选，一次拿到能量 + 受力 + 优化轨迹：

atoms = read("vasprun.xml", index=-1)   # -1 取最后一步
E = atoms.get_potential_energy()        # 能量，eV
F = atoms.get_forces()                  # 受力，eV/Å

原因：vasprun.xml 是 VASP 的结构化输出，ASE 会附上一个 SinglePointCalculator，所以能量受力直接可用，比手动 grep OUTCAR 稳得多。

6. 可视化

6.1 弹出交互窗口（ase-gui）

from ase.visualize import view

view(atoms)                    # 打开 ase-gui，可旋转缩放
view([atoms1, atoms2])         # 多结构对比

view 会启动 ase-gui 进程。在无显示的环境（远程服务器、某些 WSL）会失败，此时改用下面的 PNG 导出。

6.2 导出 PNG

from ase.io import write

write("slab.png", atoms)                              # 默认视角
write("slab.png", atoms, rotation="10z,-80x")         # 先绕 z 转 10°，再绕 x 转 -80°
write("slab.png", atoms, show_unit_cell=2)            # 重复绘制晶胞边框（2 个周期）
write("slab.png", atoms, radii=0.6, colors=None)      # 自定义原子半径/颜色

需要 matplotlib。
rotation 字符串形如 "10z,-80x"：绕指定轴旋转指定角度，逗号分隔依次叠加。

6.3 命令行

ase gui POSCAR        # 用 GUI 打开
ase gui opt.traj      # 打开轨迹，可播放

7. 计算器接口

7.1 统一接口（最重要的抽象）

所有 ASE 计算器都实现同一套方法，挂到 Atoms 上即可：

atoms.calc = some_calculator
E = atoms.get_potential_energy()   # 总能量（eV）
F = atoms.get_forces()             # 受力 (N,3)（eV/Å）
S = atoms.get_stress()             # 应力，Voigt 6 分量（eV/Å³）

受力是原子受力（VASP 里 TOTEN 对应的力的负梯度方向约定以 ASE 为准）。能量是"总电子能量"，不包含零点能、熵修正——吸附能等热力学量需要你自己额外处理，别拿裸的 get_potential_energy() 当自由能。

7.2 EMT —— 可真实运行的快速示例

EMT（有效介质理论）是纯 Python 的半经验势，只支持有限的几种金属（Al, Cu, Ag, Au, Ni, Pd, Pt）。适合本教程里"真正跑起来"的演示：

from ase.build import bulk
from ase.calculators.emt import EMT

atoms = bulk("Cu", "fcc", a=3.61)
atoms.calc = EMT()
print("Cu 单原子能量 =", atoms.get_potential_energy(), "eV")

# 受力演示：手动挪一个原子制造不平衡
atoms[0].position += (0.1, 0.0, 0.0)
print("受力 =", atoms.get_forces())

关键坑：EMT 不能算 C、O、H、N 等元素，molecule("CO2") 配 EMT 会直接报"参数缺失"。所以涉及 CO2 的能量/受力计算必须换真实 DFT 计算器（VASP/GPAW），EMT 只用于金属体系演示。

7.3 VASP 计算器配置（模板，需装 VASP）

from ase.calculators.vasp import Vasp
from ase.io import read

atoms = read("POSCAR")

calc = Vasp(
    xc="PBE",              # 泛函
    encut=450,             # 平面波截断（eV）
    kpts=(4, 4, 1),        # k 点网格（slab 常用 z 方向 1）
    gamma=True,            # Gamma 中心 k 点网格
    ismear=0,              # 展宽：半导体/分子用 0，金属用 1
    sigma=0.05,            # 展宽宽度（eV）
    ibrion=2,              # 离子弛豫算法：2 = 共轭梯度 CG
    nsw=100,               # 最大离子步数
    isif=2,                # 2 = 只弛豫离子（保持晶胞），3 = 离子+晶胞
    prec="Normal",         # 计算精度
    ediff=1e-5,            # 电子步收敛
    nelm=100,              # 最大电子步数
    ncore=8,               # 并行核数
    lwave=False,           # 不写 WAVECAR
    lcharg=False,          # 不写 CHGCAR
    setups="recommended",  # 推荐 PAW 赝势
)

atoms.calc = calc
e = atoms.get_potential_energy()   # 触发 VASP 运行，返回能量

几个参数含义补充（与你 VASP 基础对应）：

参数	含义
`xc`	交换关联泛函，`PBE` / `RPBE` / `SCAN` 等
`encut`	ENCUT，平面波截断
`kpts`	KPOINTS；可给 tuple `(4,4,1)`，也可给密度字典 `{"density": 2.0, "gamma": True}`
`ibrion`	`-1` 不动 / `0` MD / `1` quasi-Newton / `2` CG
`isif`	`0` 不弛豫 / `2` 弛豫离子 / `3` 弛豫离子+晶胞 / `4` 只弛豫胞形
`ismear`	`0` Gaussian（半导体/分子）/ `1` Methfessel-Paxton（金属）

运行前环境要求：

系统 PATH 里要有 vasp_std（或用 VASP_COMMAND / ASE_VASP_COMMAND 环境变量指定可执行文件路径）。
VASP_PP_PATH 环境变量指向你的 POTCAR 目录（或用 setups 指定）。

# Windows (PowerShell) 示例
$env:VASP_PP_PATH = "D:\vasp\potpaw_PBE"
$env:ASE_VASP_COMMAND = "D:\vasp\bin\vasp_std.exe"

ASE 的 Vasp 计算器只是"帮你写 INCAR/POSCAR/KPOINTS/POTCAR 并调 vasp、解析结果"，它不含 VASP 本体。本教程 VASP 部分均为配置示例，实际运行需你本机装好 VASP + 赝势。

7.4 计算器通用参数与目录管理

calc = Vasp(xc="PBE", encut=400, kpts=(2, 2, 2))
calc.directory = "run_relax"     # 在该子目录里跑，结果也写进去
atoms.calc = calc

# 换个参数复用同一结构
atoms.calc.set(encut=500)        # 覆盖已有参数

calc.set(**kwargs) 是通用方法，所有计算器都支持，用来在循环里批量改参数。

8. 几何优化

ASE 内置多种优化器，统一用 run(fmax=..., steps=...) 驱动：

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.optimize import BFGS, LBFGS, FIRE

atoms = bulk("Cu", "fcc", a=3.61)
atoms.calc = EMT()

dyn = BFGS(atoms, trajectory="opt.traj", logfile="opt.log")
dyn.run(fmax=0.05, steps=200)    # 收敛判据：最大受力 < 0.05 eV/Å

8.1 常用优化器

优化器	特点	适用
`BFGS`	拟牛顿，稳健，收敛快	默认首选，小到中型体系
`LBFGS`	有限内存 BFGS，省内存	大体系、长轨迹
`FIRE`	快速惯性弛豫，无需 Hessian	粗优化、MD 后接力的快速下降
`GPMin` / `MDMin`	较简单	NEB 里常用 MDMin

导入路径：

from ase.optimize import BFGS, LBFGS, FIRE, GPMin, MDMin, QuasiNewton

8.2 收敛判据

fmax：最大原子受力阈值（eV/Å）。默认 0.1，做吸附能建议收紧到 0.03–0.05 eV/Å，甚至 0.02。
steps：最大迭代步数上限，防止死循环。

8.3 轨迹 traj 的回读与续算

# 优化完读取最终结构
from ase.io import read
opt = read("opt.traj", index=-1)    # 最后一帧 = 优化结果
print("优化后能量 =", opt.get_potential_energy(), "eV")

续算（从一个已有轨迹继续）：

dyn = BFGS(atoms, trajectory="opt.traj")   # 若文件已存在，默认 append 模式续写

8.4 带约束的优化

slab 弛豫通常要固定底层原子（详见第 10 节），约束会贯穿优化器自动生效。

9. 单点能与受力、能量/受力单位换算（ase.units）

9.1 单点计算

from ase.build import bulk
from ase.calculators.emt import EMT

atoms = bulk("Cu", "fcc", a=3.61)
atoms.calc = EMT()

E = atoms.get_potential_energy()      # eV
F = atoms.get_forces()                # eV/Å，形状 (N, 3)
S = atoms.get_stress()                # eV/Å³，Voigt 六分量 (xx,yy,zz,yz,xz,xy)

受力方向：F[i] 是第 i 个原子受的合力矢量；平衡结构里应趋近 0。
应力符号约定：ASE 里正值表示拉伸（张应力）。VASP 转过来时注意符号一致性。

9.2 ase.units 换算

from ase.units import eV, kJ, mol, kcal, Hartree, Bohr, Angstrom, nm

# 能量换算：1 eV = ?
E = 1.0                                   # eV
print(E / (kJ / mol), "kJ/mol")           # 96.49 kJ/mol
print(E / (kcal / mol), "kcal/mol")       # 23.06 kcal/mol
print(E / Hartree, "Hartree")             # 0.03675 Hartree

# 长度换算：1 Å = ?
L = 1.0                                   # Å
print(L / Bohr, "Bohr")                   # 1.890 Bohr
print(L / nm, "nm")                       # 0.1 nm

常用换算数值（背下来省事）：

关系	数值
1 eV	96.485 kJ/mol
1 eV	23.06 kcal/mol
1 Hartree	27.211 eV
1 Bohr	0.5292 Å
1 eV/Å	1.602×10⁻⁹ N（力）

吸附能常用 kJ/mol 汇报，但 ASE 内部能量是 eV，所以脚本里统一用 eV 算，最后一步再乘 96.485 转成 kJ/mol。

10. 约束：FixAtoms（吸附结构弛豫底层原子）

slab 模型弛豫的黄金法则：固定底层 1–2 层，只弛豫表面几层和吸附物，模拟体相"刚性衬底"。

10.1 用 layer tag 固定底层

fcc111 等便捷函数会自动给每层打 tag（底 = 1，往上递增）：

from ase.build import fcc111, molecule, add_adsorbate
from ase.constraints import FixAtoms

slab = fcc111("Pt", size=(2, 2, 4), a=3.92, vacuum=10.0)  # 4 层
co2  = molecule("CO2")
add_adsorbate(slab, co2, height=2.0, position="ontop")    # 吸附物 tag=0

# 固定 tag <= 2 的原子（即底部两层）；顶层(tag 3)和吸附物(tag 0)自由
mask = [atom.tag <= 2 for atom in slab]
slab.set_constraint(FixAtoms(mask=mask))
print("固定原子数:", sum(mask), " 自由原子数:", len(slab) - sum(mask))

10.2 用索引固定

当你从 POSCAR 读入、没有 tag 信息时，用原子索引（按 z 坐标排序取底层）：

import numpy as np
from ase.io import read
from ase.constraints import FixAtoms

atoms = read("POSCAR")
z = atoms.positions[:, 2]                       # 取 z 坐标
bottom_indices = np.argsort(z)[:N_bottom]       # z 最小的 N_bottom 个原子
atoms.set_constraint(FixAtoms(indices=bottom_indices.tolist()))

10.3 约束效果

挂上 FixAtoms 后，优化器 / 受力计算会自动把这些原子的受力置零，位置冻结。检查：

print(atoms.constraints)               # 查看已挂约束

配合第 8 节的 BFGS，就是标准的"固定衬底 + 弛豫表面/吸附物"流程。

10.4 其他约束（简介）

from ase.constraints import FixAtoms, FixScaled, FixCartesian, FixedPlane

FixScaled：固定分数坐标（VASP 的 Selective Dynamics + T 类似）。
FixCartesian：只锁某些笛卡尔分量（如锁 z 方向）。

VASP 用户常想直接写 Selective Dynamics，ASE 的 write("POSCAR", atoms, direct=True) 会自动把 FixAtoms/FixScaled 约束翻译成 VASP 的 Selective Dynamics 标记。

11. 进阶：NEB 过渡态、分子动力学、声子

11.1 NEB 过渡态搜索（简要）

from ase.neb import NEB
from ase.optimize import MDMin
from ase.build import fcc111
from ase.calculators.emt import EMT

# 初态 / 末态（示意：同一 slab 上两个吸附位）
initial = fcc111("Cu", size=(2, 2, 3), a=3.61, vacuum=10.0)
final   = initial.copy()
# ... 这里把吸附原子从初态位搬到末态位（实际需手动移动某原子）...

# 插值生成中间像
images = [initial]
for _ in range(3):
    images.append(initial.copy())
images.append(final)

for image in images:
    image.calc = EMT()

neb = NEB(images)
neb.interpolate()                      # 线性/IDPP 插值生成初始链
dyn = MDMin(neb, trajectory="neb.traj")
dyn.run(fmax=0.05)

# 能量势垒 = 最高像能量 - 初态能量
E_barrier = max(img.get_potential_energy() for img in images) - images[0].get_potential_energy()
print("势垒 =", E_barrier, "eV")

更精确的过渡态常用 CI-NEB（climbing image）：NEB(images, climb=True)。NEB 示例里的 interpolate() 默认方法在新版本已用 IDPP，具体参数以官方文档为准。EMT 只适合金属，CO2 解离之类的反应请用真实 DFT 计算器。

11.2 分子动力学（简要）

Langevin（恒温，近似正则系综）：

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.md.langevin import Langevin
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase import units

atoms = bulk("Cu", "fcc", a=3.61).repeat((2, 2, 2))
atoms.calc = EMT()

# 按 300 K 麦克斯韦分布初始化速度
MaxwellBoltzmannDistribution(atoms, temperature_K=300)

dyn = Langevin(
    atoms,
    timestep=1.0 * units.fs,     # 时间步长（ASE 时间单位）
    temperature_K=300,           # 目标温度 K
    friction=0.01,               # 摩擦系数
)
dyn.run(1000)                    # 跑 1000 步

VelocityVerlet（微正则 NVE，无恒温）：

from ase.md.verlet import VelocityVerlet

dyn = VelocityVerlet(atoms, timestep=1.0 * units.fs)
dyn.run(1000)

ASE 时间单位约 10.18 fs，所以 1.0 * units.fs 才是"1 飞秒"（units.fs ≈ 0.098）。这是新手最容易错的点：直接写 timestep=1.0 等于 10.18 fs，太大了。

11.3 声子（简要）

from ase.build import bulk
from ase.calculators.emt import EMT
from ase.phonons import Phonons

atoms = bulk("Al", "fcc", a=4.05)

ph = Phonons(atoms, EMT(), supercell=(2, 2, 2), delta=0.05)
ph.run()                        # 有限位移法算力常数
ph.read(acoustic=True)          # 声学支置零
ph.clean()                      # 清理临时文件
ph.write_dos()                  # 写声子 DOS

# 声子能带
band = ph.get_band_structure()  # 具体返回对象以官方文档为准

delta 是有限位移步长（Å）。声子对超胞尺寸敏感，(2,2,2) 通常不够，实际科研用 (4,4,4) 甚至更大。

12. 科研实战示例

12.1 构建 CO2 在金属表面的吸附 slab 模型

"""
CO2 吸附在 Pt(111) 上的 slab 建模脚本
说明：本脚本只负责"搭结构 + 固定衬底 + 写 POSCAR"，
      能量计算需配合真实 DFT 计算器（VASP/GPAW），EMT 不支持 C/O。
"""
from ase.build import fcc111, molecule, add_adsorbate
from ase.constraints import FixAtoms
from ase.io import write

# 1) 构建 Pt(111) slab：2x2，4 层，10 Å 真空
slab = fcc111("Pt", size=(2, 2, 4), a=3.92, vacuum=10.0)

# 2) 构建 CO2 分子，吸附到 top 位
co2 = molecule("CO2")
add_adsorbate(slab, co2, height=2.0, position="ontop")

# 3) 固定底部两层
mask = [atom.tag <= 2 for atom in slab]
slab.set_constraint(FixAtoms(mask=mask))

# 4) 写出 POSCAR（VASP5 风格 + 分数坐标）
write("POSCAR_co2_pt111", slab, vasp5=True, direct=True)
print("原子总数:", len(slab), " 固定:", sum(mask))

12.2 读 VASP POSCAR 提取能量并批量后处理

"""
批量读多个 VASP 计算的 vasprun.xml，提取最后一步能量。
目录结构假设：
    ./calc_01/vasprun.xml
    ./calc_02/vasprun.xml
    ...
"""
import glob, os
from ase.io import read

for d in sorted(glob.glob("calc_*")):
    vrun = os.path.join(d, "vasprun.xml")
    if not os.path.exists(vrun):
        print(f"[跳过] {d} 无 vasprun.xml")
        continue
    atoms = read(vrun, index=-1)              # 最后一步
    E = atoms.get_potential_energy()          # eV
    Fmax = (atoms.get_forces() ** 2).sum(axis=1).max() ** 0.5  # 最大受力
    print(f"{d}: E = {E:.4f} eV, Fmax = {Fmax:.4f} eV/Å")

12.3 吸附能计算（标准流程）

吸附能定义：E_ads = E(slab+ads) − E(slab) − E(ads)，负值表示吸附放热（稳定）。

"""
吸附能计算模板：需要三个独立 VASP 单点/优化计算的能量。
三个结构：吸附体系、干净 slab、孤立吸附分子（放足够大盒子）。
"""
from ase.io import read
from ase.units import kJ, mol

def energy_from_vrun(path):
    """从 vasprun.xml 读取最后一步总能量（eV）"""
    atoms = read(path, index=-1)
    return atoms.get_potential_energy()

E_slab_ads = energy_from_vrun("slab+co2/vasprun.xml")
E_slab     = energy_from_vrun("slab/vasprun.xml")
E_ads      = energy_from_vrun("co2_box/vasprun.xml")   # 孤立 CO2，放 15~20 Å 盒子

E_adsorb = E_slab_ads - E_slab - E_ads

print(f"吸附能 = {E_adsorb:.4f} eV = {E_adsorb / (kJ / mol):.1f} kJ/mol")

实际科研注意：裸的电子吸附能还要做零点能（ZPE）、vdW 修正（CO2 弱吸附时很重要，如用 DFT-D3 / vdW-DF）等校正，且干净 slab 和吸附体系必须用完全相同的 slab 层数、超胞大小、k 点、encut，否则误差相减会放大。

13. 常见坑

POSCAR 晶格 vs 笛卡尔坐标：read("POSCAR") 后 atoms.positions 永远是笛卡尔坐标（Å），ASE 已把 Direct 转好，别自己再手动乘晶格。写回时若要 Direct 用 write(..., direct=True)。
PBC 设置遗漏：读 POSCAR 通常自带 PBC，但手写 Atoms 或读 XYZ 后 PBC 默认是 False。做 slab/晶体计算前务必检查 atoms.pbc，必要时 atoms.set_pbc([True, True, True])。
units 混淆：ASE 一律 eV/Å。LAMMPS 用户转过来时尤其注意（LAMMPS 常用 metal units，但力单位是 eV/Å 或不同）；换算用 ase.units，别手写错 96.485 这个因子。
EMT 不支持 C/O/H/N：molecule("CO2") 配 EMT() 会报错。EMT 只用于 Al/Cu/Ag/Au/Ni/Pd/Pt 等金属。
timestep 单位：MD 里 timestep=1.0 是 10.18 fs，要写 1.0 * units.fs 才是 1 fs。
吸附能相减误差：三个计算必须同一 slab、同一参数，否则相减放大误差。见 12.3 提醒。
ase.build 没有 fcc()/bcc()：块体用 bulk(name, 'fcc'/'bcc', a=...)；fcc111/bcc100 是表面函数。
get_potential_energy() 是电子能：不是自由能，缺 ZPE、熵、vdW 校正，汇报热力学量时别直接用。
vasprun.xml 优于 OUTCAR：读能量用 read("vasprun.xml", index=-1)，比 grep OUTCAR 稳；OUTCAR 作为 fallback。
受力方向/应力符号：ASE 应力正值=拉伸；和 VASP 对接时确认符号约定一致。

14. 延伸资源

ASE 官方文档（最重要）：https://wiki.fysik.dtu.dk/ase/
- 结构构建：ase.build 模块页
- 文件读写：ase.io 模块页
- VASP 计算器：ase.calculators.vasp 页（参数列表以这里为准）
- 单位：ase.units 页
ASE 官方教程：https://wiki.fysik.dtu.dk/ase/tutorials/tutorials.html （sulfur、NEB、MD 等经典入门）
命令行工具：ase --help、ase gui --help、ase info
Materials Project / ase.io：MP 导出的 CIF/结构可直接 read()，配合 pymatgen 也可（但本教程纯 ASE 已够用）
GPAW：https://wiki.fysik.dtu.dk/gpaw/ —— ASE 原生的 DFT 代码，无需 VASP 即可跑真实 DFT（对没有 VASP 的环境很有用）

最后一句经验：遇到不确定的函数/参数，第一反应应该是 python -c "import ase.build as b; help(b.surface)" 或查官方文档，而不是凭记忆猜 API。ASE 版本迭代快，本教程以 3.28.0 为准。

ZTPchemistry

Python ASE库教程