Python 科研编程教程

1. 定位说明：本教程在科研编程中的角色

计算化学的日常代码，绝大多数是三类脚本，而非大型软件工程：

  1. 数据搬运与清洗：把 VASP/Gaussian/CP2K/LAMMPS 的输出文本、CSV、JSON 变成结构化的数字。
  2. 批量流程控制：对几十上百个结构文件做重命名、改参数、批量提交作业、汇总结果。
  3. 分析与可视化：算完的能量/力/吸附能，整理成表格、画图（配合 matplotlib、numpy、ase、pymatgen）。

本教程的目标是让你写出的脚本具备三个工程品质：

  • 可复现：参数从命令行/配置文件传入，而不是写死在代码里。
  • 可维护：一个月后回来看还看得懂（命名、类型注解、函数拆分）。
  • 健壮：文件缺失、编码错误、计算没收敛时能给出清晰报错，而不是默默产出错误数据。

原则：科研代码的正确性永远排在“优雅”前面。 一个能量取错符号的漂亮脚本，比一个能跑的丑脚本危险一百倍。

2. 环境与基础

2.1 虚拟环境（conda vs venv）

不同项目对依赖版本要求不同（例如 ase 的某版本改了 VASP 接口），永远不要在 base 环境里混装所有包。

# conda 方式（推荐，配合 Miniconda）
conda create -n dft python=3.14 -y
conda activate dft
conda install -c conda-forge numpy scipy pandas matplotlib ase -y

# pip 方式（在已激活的 conda 环境里也可用 pip）
python -m venv .venv
# Windows 激活：
.venv\Scripts\activate
pip install -r requirements.txt

# 导出 / 复现环境（科研复现的命根子）
conda env export > environment.yml      # conda 方式
pip freeze > requirements.txt           # pip 方式

2.2 运行方式

python script.py            # 直接跑脚本
python -i script.py         # 跑完进入交互式，变量还留着（调试神器）
python -m pdb script.py     # 用内置调试器
python -m script            # 把脚本当模块跑（需要脚本在可导入路径下）

在 VS Code 里，用「Run and Debug」打断点比满屏 print 高效得多；print 留给快速探查即可。

2.3 PEP8 风格要点（只讲最常犯的）

  • 4 空格缩进，不用 Tab。
  • 每行 ≤ 79 字符（VS Code 装 Ruff 扩展可自动提示/格式化）。
  • 变量/函数用 snake_case，类用 PascalCase，常量用 UPPER_CASE。
  • 导入顺序：标准库 → 第三方 → 本地，组间空一行。
import os
import re
from pathlib import Path

import numpy as np

from my_pkg.vasp import parse_outcar

2.4 Python 3.14 值得知道的新特性

  • 模板字符串（PEP 750）：t"..." 返回 string.templatelib.Template 对象而非 str，把「模板」和「插值值」分开存储，适合做安全 HTML/SQL 拼接或结构化日志。日常科研用得少，了解即可。
name = "CO2"
tpl = t"吸附物：{name}"     # 返回 Template 对象，不是 str
print(tpl.strings)         # ('吸附物：', '')
print(tpl.values)          # ('CO2',)
  • 注解惰性求值（PEP 649）：3.14 起函数/类的类型注解默认延迟求值，定义时不会真的去执行注解表达式，前向引用（类还没定义就拿来当注解）不再需要 from __future__ import annotations。注意访问 fn.__annotations__ 时才真正求值。

版本敏感速查：

  • 3.10 引入 match 结构化匹配、int | str 联合类型写法。
  • 3.11 引入异常组（ExceptionGroup / except*）、BaseException.add_note()、更清晰的 traceback。
  • 3.12 增强 f-string（任意表达式、嵌套引号、多行、表达式内注释）。
  • 3.13 引入 copy.replace()、自由线程（实验性，默认关闭）。
  • 3.14 模板字符串、注解惰性求值。

3. 内置数据结构

3.1 list

energies = [-3.1, -3.2, -3.0, -3.15]
energies.append(-3.25)          # 追加
energies.extend([-3.3, -3.35])  # 扩展多个
energies.insert(0, -2.9)        # 指定位置插入
last = energies.pop()           # 弹出并返回末尾
energies.remove(-3.0)           # 按值删除（只删第一个）

# 排序（计算化学常用：按能量从低到高 = 数值从小到大）
energies.sort()                  # 原地排序
ranked = sorted(energies)        # 返回新列表
ranked_rev = sorted(energies, reverse=True)

# 切片：取能量最低的三个
lowest3 = sorted(energies)[:3]

3.2 dict（科研里的“万能容器”）

results = {"formula": "Cu(111)", "E_total": -432.1, "converged": True}

# 取值
e = results["E_total"]               # 键不存在会 KeyError
e = results.get("E_total")           # 不存在返回 None
e = results.get("E_total", 0.0)      # 不存在返回默认值（推荐）

# setdefault：键不存在时写入默认值并返回，存在则返回已有值（计数器/分组高频）
counts = {}
for site in ["bridge", "top", "bridge", "hollow", "top"]:
    counts[site] = counts.setdefault(site, 0) + 1

# 合并字典（3.9+ 的 | 运算符）
defaults = {"encut": 400, "kpts": [3, 3, 1]}
user = {"encut": 520}
merged = defaults | user           # {'encut': 520, 'kpts': [3, 3, 1]}

# 遍历
for k, v in results.items():
    print(k, v)

3.3 set（去重、集合运算）

ads_sites_A = {"top", "bridge", "fcc"}
ads_sites_B = {"bridge", "hcp"}

ads_sites_A | ads_sites_B   # 并集
ads_sites_A & ads_sites_B   # 交集
ads_sites_A - ads_sites_B   # 差集

# 去重（元素类型必须可哈希；dict/list 不可哈希）
unique = set([1, 2, 2, 3, 3, 3])    # {1, 2, 3}

3.4 tuple（不可变，常用于坐标、常数打包）

lattice = (3.63, 3.63, 3.63)       # 立方晶胞参数
a, b, c = lattice                   # 解包
# tuple 可作为 dict 键（因为不可变、可哈希）
coords = {(0.0, 0.0, 0.0): "Cu", (0.5, 0.5, 0.5): "O"}

4. 字符串与文件

4.1 格式化（f-string）

formula = "Cu"
E_ads = -1.2345

# 基本
s = f"{formula} 上的吸附能为 {E_ads} eV"

# 精度 / 宽度 / 对齐（写表格输出必备）
s = f"{formula:>6} {E_ads:12.4f}"     # 右对齐宽度6；总宽12、4位小数
s = f"{0.001234:.3e}"                 # 科学计数法

# 3.12+：表达式内可用嵌套引号、注释、多行
s = f"{E_ads:.{2}f}"                  # 精度也可是表达式

4.2 pathlib（推荐用它取代 os.path）

from pathlib import Path

p = Path("D:/计算/CO2/Cu_111")        # 正斜杠在 Windows 下也合法
p = Path("data") / "POSCAR" / "relax"  # 用 / 拼路径

p.exists(); p.is_file(); p.is_dir()
p.name           # 'relax'
p.stem           # 'relax'（去后缀）
p.suffix         # ''（无后缀则为空）
p.parent         # PosixPath('data/POSCAR')
p.with_suffix(".vasp")     # 换后缀
p.resolve()                 # 绝对路径

# 遍历
for f in Path(".").glob("**/*.cif"):        # 递归匹配
    print(f)
for f in Path(".").rglob("OUTCAR"):         # rglob 等价于 **/OUTCAR
    print(f)

# 读写
text = Path("INCAR").read_text(encoding="utf-8")
Path("out.txt").write_text(text, encoding="utf-8")

4.3 读写文本

# 推荐：with 自动关闭；显式指定 encoding（Windows 默认 GBK，见第 13 节）
with open("OUTCAR", encoding="utf-8", errors="replace") as f:
    content = f.read()          # 一次读完

with open("POSCAR", encoding="utf-8") as f:
    lines = f.readlines()       # 读成 list

# 逐行读（大文件，避免一次吃满内存）
with open("big.log", encoding="utf-8") as f:
    for line in f:
        if "TOTEN" in line:
            ...

4.4 读写 CSV

import csv

# 写：用 csv.writer 处理分隔符/引号转义
with open("results.csv", "w", newline="", encoding="utf-8") as f:
    w = csv.writer(f)
    w.writerow(["formula", "E_ads_eV", "site"])
    w.writerow(["Cu(111)", -1.2345, "bridge"])

# 读
with open("results.csv", encoding="utf-8") as f:
    for row in csv.DictReader(f):      # 首行作列名，返回 dict
        print(row["formula"], row["E_ads_eV"])

4.5 读写 JSON（配置、结果交换通用格式）

import json

data = {"encut": 520, "kpts": [3, 3, 1], "functional": "PBE"}

with open("config.json", "w", encoding="utf-8") as f:
    json.dump(data, f, indent=2, ensure_ascii=False)

with open("config.json", encoding="utf-8") as f:
    loaded = json.load(f)

# 注意：JSON 里没有 tuple，会变成 list；没有 NaN（可用 allow_nan=True 写出 NaN/Infinity）

5. 函数

5.1 参数

def calc_E_ads(E_slab_ads, E_slab, E_gas, *, n=1.0, unit="eV"):
    """吸附能 = (E_slab+ads − E_slab − n·E_gas) / n，吸附能越负越稳。"""
    return (E_slab_ads - E_slab - n * E_gas) / n

calc_E_ads(-435.1, -430.0, -3.1)             # 位置参数
calc_E_ads(-435.1, -430.0, -3.1, n=1.0)      # 关键字参数
# * 之后的参数只能以关键字形式传，防止 n 被误当位置参数

def submit(name, *args, **kwargs):
    """*args 收多余位置参数为 tuple；**kwargs 收多余关键字为 dict。"""
    print(f"提交任务 {name}")
    for a in args:
        print("位置参数:", a)
    for k, v in kwargs.items():
        print(f"参数 {k} = {v}")

submit("relax", 400, 520, kpts=[3, 3, 1], xc="PBE")
# 解包反向操作：submit(*tup, **d)

5.2 默认参数陷阱（经典坑，见第 13 节）

# 错误：默认值在函数定义时只创建一次，list 会被多次调用共享！
def bad(paths=[]):
    paths.append("x")
    return paths

# 正确：用 None + 在函数体内创建
def good(paths=None):
    paths = paths or []
    paths.append("x")
    return paths

5.3 lambda（短匿名函数）

# 常用于 sorted 的 key
species = [("Cu", 63.5), ("O", 16.0), ("C", 12.0)]
sorted(species, key=lambda x: x[1])     # 按原子量排序
# key= 其实可以直接给函数引用
sorted(species, key=lambda x: x[0])

lambda 只能写单条表达式；复杂逻辑请写 def，别硬塞进 lambda。

5.4 作用域（LEGB）

查找顺序：Local → Enclosing（闭包）→ Global → Builtin。

x = "global"
def outer():
    x = "enclosing"
    def inner():
        x = "local"     # 本地赋值，不影响外层
        print(x)
    inner()
    print(x)            # enclosing
outer()

# 在函数里修改全局/外层变量需显式声明（科研脚本少用，多传参更清晰）
count = 0
def inc():
    global count
    count += 1

5.5 类型注解

from typing import Iterable, Optional, Callable

def parse_energies(lines: Iterable[str], scale: float = 1.0) -> list[float]:
    """从若干行里提取能量并缩放。"""
    ...

def get_ref(system: str) -> Optional[dict]:
    """可能返回 None 的引用数据。"""
    ...

# 3.10+ 可用 | 代替 Optional
def f(x: int | None) -> int | float:
    ...

注解只是文档/提示，运行时不强制检查；要真正校验需用 pydantic 或手写 assert。3.14 起注解惰性求值（见 2.4）。

6. 推导式与迭代器 / 生成器

6.1 推导式

# 列表推导
energies = [float(line.split("=")[1]) for line in lines if "TOTEN" in line]

# 字典推导
label = {"C": "carbon", "O": "oxygen"}
reverse = {v: k for k, v in label.items()}

# 集合推导
unique = {round(e, 3) for e in energies}

# 生成器表达式（圆括号，惰性，不一次性占内存）
gen = (e * 27.2114 for e in energies)   # Hartree -> eV
total = sum(gen)                         # 消费时才计算

何时用生成器表达式而非列表推导？数据量大、且只需要遍历一次时（例如逐条求和）。

6.2 迭代器与生成器（yield）

# 生成器函数：读大文件时按“逻辑块”产出，而不是一次性载入
def iter_convergence(path):
    """逐块产出 OUTCAR 里每次电子步迭代的能量。"""
    with open(path, encoding="utf-8") as f:
        for line in f:
            if "free  energy   TOTEN" in line:
                yield float(line.split("=")[-1].split("eV")[0].strip())

for e in iter_convergence("OUTCAR"):
    print(e)

# 判断收敛：最后两个电子步能量差 < 阈值
import itertools
es = iter_convergence("OUTCAR")
last2 = list(itertools.islice(es, -2, None))   # 取最后两个（需先耗尽生成器）
# 更直接：记录所有值再比较

yield 让函数变成一个惰性序列：调用时函数体不执行，直到 next()/for 时才跑到下一个 yield 并暂停，现场（局部变量）保留。

def countdown(n):
    while n > 0:
        yield n
        n -= 1

c = countdown(3)
print(next(c))   # 3
print(next(c))   # 2
print(list(c))   # [1]

7. 面向对象（科研里封装结构 / 参数很常用）

7.1 类与继承

class Molecule:
    """一个分子/团簇的最小封装。"""
    def __init__(self, name, symbols, coords):
        self.name = name
        self.symbols = symbols
        self.coords = coords

    def formula(self):
        from collections import Counter
        c = Counter(self.symbols)
        # C1 H2 O1 -> "CH2O" 的顺序按自定义规则
        order = ["C", "H", "O", "N", "Cu"]
        parts = [el if c[el] == 1 else f"{el}{c[el]}" for el in order if c[el]]
        return "".join(parts)

    def n_atoms(self):
        return len(self.symbols)

mol = Molecule("formaldehyde", ["C", "H", "H", "O"], [[0,0,0],[0,0,1],[1,0,0],[0,1,0]])
print(mol.formula(), mol.n_atoms())

class Slab(Molecule):
    """继承 Molecule，扩展表面模型特有属性。"""
    def __init__(self, name, symbols, coords, vacuum=15.0):
        super().__init__(name, symbols, coords)
        self.vacuum = vacuum

    def surface_area(self, a, b):
        return a * b

7.2 @property（把方法当属性读，且可加校验）

class Metal:
    def __init__(self, name, atomic_mass):
        self.name = name
        self._mass = atomic_mass

    @property
    def mass(self):
        return self._mass

    @mass.setter
    def mass(self, value):
        if value <= 0:
            raise ValueError("原子量必须为正")
        self._mass = value

cu = Metal("Cu", 63.5)
print(cu.mass)      # 像属性一样读
cu.mass = 64.0      # 走 setter 校验

7.3 dataclass（科研首选：自动生成 __init__/__repr__/__eq__）

from dataclasses import dataclass, field

@dataclass
class DFTJob:
    """一次 DFT 计算任务的全部参数。"""
    directory: str
    encut: float = 400.0
    xc: str = "PBE"
    kpts: list = field(default_factory=lambda: [3, 3, 1])  # 可变默认值必须用 default_factory
    converged: bool = False

    @property
    def label(self):
        return f"{self.directory}_{self.encut:.0f}"

job = DFTJob("Cu_111_relax", encut=520)
print(job)                      # 可读的 repr
job2 = DFTJob("Cu_111_relax", encut=520)
print(job == job2)              # 自动按字段比较

# 不可变版本：@dataclass(frozen=True) —— 适合当“常量参数包”传遍全场

用 field(default_factory=list) 而不是 kpts: list = []——这正是第 5.2/13 节可变默认参数陷阱在类层面的翻版。

8. 异常处理与日志

8.1 try / except / finally

try:
    energy = parse_energy("OUTCAR")
except FileNotFoundError:
    print("OUTCAR 不存在，跳过")
except ValueError as e:
    print(f"能量格式错误：{e}")
else:
    print("成功解析，无异常")     # 仅当 try 无异常时执行
finally:
    print("无论是否异常都执行（收尾：关文件、删临时文件）")

捕获要具体（FileNotFoundError、ValueError），别裸 except: 吞掉所有错误——那会把 bug 也藏起来。

# 主动抛出 + 附加说明（3.11+）
def require_converged(e_hist, tol=1e-4):
    if abs(e_hist[-1] - e_hist[-2]) > tol:
        raise RuntimeError(f"未收敛：最后两步能量差 {abs(e_hist[-1]-e_hist[-2]):.6f}")

8.2 自定义异常

class NotConvergedError(RuntimeError):
    """计算未收敛。"""
    def __init__(self, system, e_diff):
        self.system = system
        self.e_diff = e_diff
        super().__init__(f"{system} 未收敛，ΔE = {e_diff:.6f} eV")

raise NotConvergedError("Cu_111", 0.0012)

8.3 logging（替代满屏 print 的正规做法）

import logging

logging.basicConfig(
    level=logging.INFO,
    format="%(asctime)s [%(levelname)s] %(message)s",
    datefmt="%H:%M:%S",
)

logger = logging.getLogger("dft")

logger.debug("调试信息（默认不显示）")
logger.info("开始解析 OUTCAR")
logger.warning("k 点网格较稀疏")
logger.error("OUTCAR 缺失，任务失败")

科研脚本里至少做到：错误用 logger.error + 抛出异常，进度用 logger.info，把「机器可解析的输出」和「给人看的日志」分开（人看的走 stderr/日志文件，别污染 stdout 的 CSV）。

9. 装饰器与上下文管理器

9.1 装饰器（给函数包一层“外套”）

import time
import functools

def timed(func):
    @functools.wraps(func)          # 保留原函数的 __name__/__doc__
    def wrapper(*args, **kwargs):
        t0 = time.perf_counter()
        result = func(*args, **kwargs)
        print(f"{func.__name__} 用时 {time.perf_counter() - t0:.3f}s")
        return result
    return wrapper

@timed
def heavy_parse(path):
    return sum(1 for _ in open(path, encoding="utf-8"))

heavy_parse("OUTCAR")

# 带参数的装饰器（外面再套一层）
def retry(n=3):
    def deco(func):
        @functools.wraps(func)
        def wrapper(*a, **kw):
            for i in range(n):
                try:
                    return func(*a, **kw)
                except Exception:
                    if i == n - 1:
                        raise
        return wrapper
    return deco

@retry(3)
def flaky_read(path):
    ...

9.2 上下文管理器（with）

# 类方式：__enter__ 进入，__exit__ 离开（保证收尾）
class WorkingDir:
    """进入某目录执行，结束自动切回（批量提交常用）。"""
    import os
    def __init__(self, path):
        self.path = path
        self._prev = None
    def __enter__(self):
        self._prev = os.getcwd()
        os.chdir(self.path)
        return self
    def __exit__(self, *exc):
        os.chdir(self._prev)

with WorkingDir("Cu_111_relax"):
    ...   # 在此目录下操作

# 生成器方式（更简洁）：@contextlib.contextmanager
from contextlib import contextmanager

@contextmanager
def timer(name):
    t0 = time.perf_counter()
    yield
    print(f"{name}: {time.perf_counter() - t0:.3f}s")

with timer("解析"):
    heavy_parse("OUTCAR")

with open(...) 就是最常见的内置上下文管理器。contextlib.suppress(FileNotFoundError) 可优雅地忽略指定异常。

10. 常用标准库速查

import os            # os.getcwd/chdir/listdir、os.path 已由 pathlib 取代大部分
import sys           # sys.argv、sys.exit(code)、sys.stdin/stdout/stderr
import json          # 结构化数据交换
import csv           # 表格数据

10.1 collections

from collections import defaultdict, Counter, deque, OrderedDict, namedtuple

# defaultdict：访问不存在的键时用工厂函数生成默认值
d = defaultdict(list)
for site in ["bridge", "top", "bridge"]:
    d[site].append(1)             # 不用先判断键是否存在

# Counter：计数
c = Counter("formaldehyde")
c.most_common(3)                   # [('a', 2), ('d', 1), ('e', 1)]（示意）

# deque：两端快速增删（当队列用）
q = deque(maxlen=100)              # 只保留最近 100 个（滑动窗口，收敛监测常用）
for e in energies:
    q.append(e)

# namedtuple：轻量只读“结构体”
Point = namedtuple("Point", ["x", "y", "z"])
p = Point(1.0, 2.0, 3.0)
print(p.x, p.y, p.z)

10.2 itertools

import itertools as it

it.product([0, 1], repeat=3)                 # 笛卡尔积：超胞组合、网格点
it.combinations([1,2,3], 2)                  # 组合（选吸附位配对）
it.permutations([1,2,3], 2)                  # 排列
it.chain(a, b)                               # 把多个可迭代对象串成一个
it.groupby(data, key=...)                    # 相邻元素分组（需先按 key 排序！）
it.islice(gen, start, stop)                  # 切片惰性序列
it.zip_longest(a, b, fillvalue=None)         # 不等长 zip
it.count(), it.cycle(), it.repeat()          # 无限迭代器（配合 islice 截断）

# 组合扫描示例：吸附位 × 覆盖率 × 泛函
for site, cov, xc in it.product(["bridge", "top"], [0.25, 0.5], ["PBE", "RPBE"]):
    print(f"{site} {cov} {xc}")

10.3 functools

import functools as ft

@ft.lru_cache(maxsize=128)          # 缓存纯函数结果（重复算能量读取时加速）
def read_energy(path: str) -> float:
    return parse_energy(path)

ft.reduce(lambda a, b: a + b, [1,2,3,4])     # 累加
ft.partial(sorted, reverse=True)             # 固定部分参数

# functools.cached_property：第一次访问后缓存（配合类使用）
from functools import cached_property
class Job:
    @cached_property
    def energy(self):
        return parse_energy("OUTCAR")        # 只解析一次

10.4 argparse（命令行工具必备，见第 12.3）

10.5 glob / re / time

from pathlib import Path
Path(".").glob("*.cif")            # pathlib 版本（推荐）

import time
time.time()          # 时间戳（秒）
time.perf_counter()  # 高精度计时（测性能用这个）

11. 正则表达式 re（解析输出文件必备）

11.1 核心语法

模式含义模式含义
.任意字符（除换行）\d数字
\w字母数字下划线\s空白
*0+ 次+1+ 次
?0/1 次{m,n}m~n 次
[abc]字符集[^abc]取反
^行首$行尾
( )捕获组(?: )非捕获组
|或(?P<name>)命名组
import re

line = "  free  energy   TOTEN  =     -432.12345678 eV"

# 匹配并捕获
m = re.search(r"TOTEN\s*=\s*([-\d.]+)\s*eV", line)
print(m.group(1))            # '-432.12345678'（第一个括号）
print(float(m.group(1)))     # -432.12345678

# 命名组更可读
m = re.search(r"TOTEN\s*=\s*(?P<energy>[-\d.]+)\s*eV", line)
print(m.group("energy"))

# findall / finditer：找所有匹配
text = open("OUTCAR").read()
all_e = re.findall(r"TOTEN\s*=\s*([-\d.]+)", text)
for m in re.finditer(r"TOTEN\s*=\s*([-\d.]+)", text):
    print(m.start(), m.group(1))

11.2 关键注意点

# 1) 用 r"..." 原始字符串，避免反斜杠被转义
# 2) 贪婪 vs 非贪婪：.* 贪婪，.*? 非贪婪
re.search(r"<tag>.*</tag>", s)    # 贪到最后一个 </tag>
re.search(r"<tag>.*?</tag>", s)   # 到第一个 </tag>

# 3) 跨行匹配：re.DOTALL 让 . 匹配换行；re.MULTILINE 让 ^$ 匹配每行
re.search(r"TOTAL-FORCE.*?total drift", text, re.DOTALL)

# 4) 预编译（同一正则反复用时更快）
pat = re.compile(r"TOTEN\s*=\s*([-\d.]+)")
pat.findall(text)

# 5) 数字科学计数法：VASP 常用 E+03 形式
re.findall(r"[-\d.Ee+]+", "-.12345678E+03")   # 别漏 E/e/+ 符号

正则最大的坑：能匹配对，但格式稍变就漏。写完后拿真实输出文件里 3~5 个不同样本验证；能用 split/strip/in 解决就别上正则。

12. 科研实战示例

12.1 批量重命名 / 处理文件

from pathlib import Path
import shutil

root = Path("D:/计算/CO2")
# 把所有子目录里的 POSCAR 复制成带目录名的备份，避免覆盖
for pos in root.rglob("POSCAR"):
    tag = "_".join(pos.parent.parts[-2:])       # 用上层目录名做标签
    shutil.copy(pos, pos.parent / f"POSCAR_{tag}.bak")
    print("备份 ->", pos.parent / f"POSCAR_{tag}.bak")

# 批量改后缀 / 清理
for f in root.rglob("*.cif"):
    f.rename(f.with_suffix(".CIF"))

12.2 解析 VASP OUTCAR / OSZICAR 提取能量

import re
from pathlib import Path

def extract_total_energy(outcar: Path) -> float:
    """从 OUTCAR 提取最后一次 'free energy TOTEN'（单位 eV）。"""
    text = outcar.read_text(encoding="utf-8", errors="replace")
    matches = re.findall(r"free\s+energy\s+TOTEN\s*=\s*([-\d.]+)", text)
    if not matches:
        raise ValueError(f"{outcar} 中找不到 TOTEN")
    return float(matches[-1])          # 取最后一次（最终结构）

def extract_energy_history(outcar: Path) -> list[float]:
    """提取所有电子步能量，用于画收敛曲线 / 判断收敛。"""
    text = outcar.read_text(encoding="utf-8", errors="replace")
    return [float(x) for x in re.findall(r"free\s+energy\s+TOTEN\s*=\s*([-\d.]+)", text)]

# 用 OSZICAR 的 E0= 作为补充（每电子步一行，更快）
def parse_oszicar(oszicar: Path) -> list[float]:
    out = []
    for line in oszicar.read_text(encoding="utf-8", errors="replace").splitlines():
        m = re.search(r"E0=\s*([-\d.Ee+]+)", line)
        if m:
            out.append(float(m.group(1)))
    return out

if __name__ == "__main__":
    e = extract_total_energy(Path("OUTCAR"))
    hist = extract_energy_history(Path("OUTCAR"))
    print(f"总能 = {e:.6f} eV")
    print(f"电子步数 = {len(hist)}")
    if len(hist) >= 2:
        dE = abs(hist[-1] - hist[-2])
        print(f"最后两步能量差 = {dE:.6f} eV", "（已收敛）" if dE < 1e-4 else "（未收敛）")

VASP 能量单位：OUTCAR 默认 eV。转 Hartree 乘 1/27.211386，转 kJ/mol 乘 96.485。处理吸附能时务必确认正负号约定（E_ads = E_slab+ads − E_slab − E_gas，负值表示放热吸附）。

12.3 写一个带 argparse 的命令行工具

"""extract_energy.py — 从 OUTCAR 提取能量并可选写 CSV。

用法：
    python extract_energy.py OUTCAR
    python extract_energy.py OUTCAR -o energy.csv --unit hartree
    python extract_energy.py OUTCAR --history
"""
import argparse
import csv
import re
from pathlib import Path

EV2HARTREE = 1 / 27.211386

def parse_args():
    p = argparse.ArgumentParser(description="从 VASP OUTCAR 提取总能量")
    p.add_argument("outcar", type=Path, help="OUTCAR 文件路径")
    p.add_argument("-o", "--output", type=Path, help="写 CSV 到该路径")
    p.add_argument("--unit", choices=["eV", "hartree"], default="eV",
                   help="能量单位（默认 eV）")
    p.add_argument("--history", action="store_true", help="输出所有电子步能量")
    p.add_argument("--tol", type=float, default=1e-4,
                   help="收敛判据（最后两步能量差）")
    return p.parse_args()

def extract(outcar: Path) -> list[float]:
    text = outcar.read_text(encoding="utf-8", errors="replace")
    vals = [float(x) for x in re.findall(r"free\s+energy\s+TOTEN\s*=\s*([-\d.]+)", text)]
    if not vals:
        raise SystemExit(f"错误：{outcar} 里没有 TOTEN 能量")
    return vals

def main():
    args = parse_args()
    if not args.outcar.exists():
        raise SystemExit(f"错误：文件不存在 {args.outcar}")

    hist = extract(args.outcar)
    factor = EV2HARTREE if args.unit == "hartree" else 1.0

    if args.history:
        for i, e in enumerate(hist, 1):
            print(f"{i:4d}  {e * factor:.8f} {args.unit}")
    else:
        total = hist[-1] * factor
        print(f"总能量 = {total:.8f} {args.unit}")
        if len(hist) >= 2:
            dE = abs(hist[-1] - hist[-2])
            status = "收敛" if dE < args.tol else "未收敛"
            print(f"最后两步 ΔE = {dE:.8f} {args.unit}（{status}）")

    if args.output:
        with open(args.output, "w", newline="", encoding="utf-8") as f:
            w = csv.writer(f)
            w.writerow(["step", f"energy_{args.unit}"])
            for i, e in enumerate(hist, 1):
                w.writerow([i, e * factor])
        print(f"已写入 {args.output}")

if __name__ == "__main__":
    main()

13. 常见坑与调试

13.1 可变默认参数

def add(x, memo=[]):      # 错！memo 在定义时创建一次，多次调用共享
    memo.append(x)
    return memo

print(add(1))   # [1]
print(add(2))   # [1, 2]  —— 预期 [2]，结果被污染了

def add_fixed(x, memo=None):
    memo = memo if memo is not None else []
    memo.append(x)
    return memo

13.2 浅拷贝 vs 深拷贝

import copy
a = [[1, 2], [3, 4]]
b = a.copy()          # 浅拷贝：外层新列表，内层仍是同一批 list
b[0][0] = 999
print(a[0][0])        # 999 —— 内层被“共享”改掉了

c = copy.deepcopy(a)  # 深拷贝：完全独立
c[0][0] = 0
print(a[0][0])        # 999 —— 不受影响

numpy 数组同理：b = a 只是别名，b = a.copy() 才是新副本。改动共享数组是科研脚本里最隐蔽的 bug 来源。

13.3 编码 GBK / UTF-8（Windows 重灾区）

# Windows 中文环境下 open() 默认用 GBK，读 UTF-8 文件会报
# UnicodeDecodeError。显式指定 encoding 是唯一稳妥做法。
with open("data.txt", encoding="utf-8") as f:      # 正确
    ...

# 写中文到文件同理；写 JSON 时加 ensure_ascii=False 让中文直接可读
import json
json.dump({"说明": "吸附能"}, f, ensure_ascii=False)

# 若确实不确定来源编码，用 errors="replace" 防止崩溃，再手动核查
text = open("legacy.dat", encoding="gbk", errors="replace").read()

提醒：Windows 的 bash（如 Git Bash）里用 curl 发微信推送中文会因 GBK 出问题，需改用 Python 做 UTF-8 编码（这是本项目已踩过的坑）。

13.4 Windows 路径

# 反斜杠会被当转义：C:\Users 里 \U、\n 都是转义序列
p = "C:\\Users\\33451\\data"     # 老式：双反斜杠
p = r"C:\Users\33451\data"       # 原始字符串
p = "C:/Users/33451/data"        # 正斜杠，Windows 也认，最省心

# 首选 pathlib，自动处理分隔符
from pathlib import Path
p = Path("D:/计算") / "CO2" / "Cu_111"

13.5 调试技巧清单

  • repr() 看字符串的真实内容（能看出隐藏的空格、\r、编码）。
  • 怀疑文件没读进来：先 print(len(text))，再 print(text[:200])。
  • 二分定位：在函数入口 print(type(x), x)，确认类型对（尤其 str vs bytes、str vs list）。
  • 用 python -i 保留现场，逐行检查中间变量。
  • breakpoint() 进入 pdb：p 变量 打印、n 下一行、c 继续。
  • 断言前置条件，把错误尽早暴露：
def parse_energy(text):
    assert isinstance(text, str), "输入必须是 str"
    ...

14. 延伸资源

  • 官方文档：docs.python.org/3/ 的 Tutorial + Library Reference（标准库的最终权威）。
  • 科研库：
    • ase（Atomic Simulation Environment）—— 读写 POSCAR/CIF、建结构、跑 VASP/LAMMPS 接口。
    • pymatgen —— Materials Project 的 Python 库，结构/相图/高通量分析。
    • numpy / scipy / pandas / matplotlib —— 数值、数据处理与绘图基础四件套。
    • ase.io 或 pymatgen.io.vasp 读 VASP 输出，比自己写正则更稳（但自己写一遍能真正理解格式）。
  • 风格与工具：ruff（lint+format）、black（格式化）、mypy（类型检查）、pytest（测试）。
  • 版本变化：docs.python.org/3/whatsnew/ 里逐版本 What's New，关注 3.10~3.14 的变化。
  • 正则调试：regex101.com（可视化匹配过程）。
  • 建议：把你反复用的解析/提交脚本沉淀成自己的 utils 模块，写一次、处处复用；给关键脚本配一个最小 pytest 测试，防止改一处坏一片。

完。写代码时记住：科研脚本的正确性 > 可读性 > 性能。先让它对，再让它清楚，最后再让它快。