Formal in Silico

§1

为什么要这样做

科学计算软件通常用三种办法建立信任：与解析解对比，与实验对比，与别的程序对比。这些办法都有效，但只覆盖被测到的情形，而且当结果出现偏差时，很难说清偏差来自哪一层。

1991 · FEM

Sleipner A 平台沉没

北海 Sleipner A 平台的混凝土基座在注水试验中破坏并沉没。事后调查指出，有限元线弹性分析中的单元形状不佳，使三格室壁的剪应力被低估约 47%，配筋因此不足。程序正常运行，没有报错。

教训：网格质量和离散误差应当是可检查的前提条件，而不是工程师的经验判断。

2016 · DFT

Δ 量规与可复现性

Lejaeghere 等人比较了 15 个固体计算代码、40 种赝势或基组方案对 71 种单质晶体的 PBE 状态方程。较新的方法之间吻合良好，较早的方案存在明显偏差，而差异的来源只能靠逐项对比去定位。

教训：结果依赖大量隐式的数值设置。这些设置应当写进规范，并与误差估计绑定。

并行化让问题更难。MD 程序在数千个 MPI 进程上做区域分解，原子跨越子域边界时如果丢失或重复，能量曲线上可能只出现一个很难察觉的跳变。浮点归约的顺序随进程数改变，同一个输入在不同规模的机器上会给出逐位不同的轨迹。这些错误与物理无关，是并发协议的错误，而并发协议正是 TLA+ 擅长处理的对象。

§2

核心主张：可信度取决于链条上最弱的一环

我们把一次模拟从方程到浮点数的过程拆成六层，并为每两层之间写下明确的证明义务。前三段义务属于 PDE 与数值分析，第四段属于软件规范，第五段属于浮点分析。

L1
物理模型
方程本身：Kohn–Sham 方程、牛顿或 Langevin 运动方程、线弹性方程。
数学适定性：解存在、唯一，并连续依赖于数据
L2
连续问题
函数空间中的变分形式，例如在 H¹ 中求 u 使 a(u, v) = ℓ(v) 对所有 v 成立。
数学离散误差估计：u_h → u，并给出收敛阶
L3
离散问题
有限维代数系统：刚度矩阵、平面波系数、粒子坐标。
数学算法收敛：迭代在第 k 步停止时，代数误差受控
L4
抽象算法
SCF 迭代、Velocity Verlet、共轭梯度法的数学描述，按精确算术理解。
TLA+ 精化：并行程序的每一步都对应抽象算法的一步，或不改变抽象状态
L5
并行程序
MPI 进程、消息、共享内存线程、GPU 流，以及检查点与重启。
浮点舍入误差界：浮点结果与精确算术结果的距离
L6
浮点执行
在 IEEE 754 硬件上实际算出来的那个数。

‖u − ũ_h^(k)‖ ≤ ‖u − u_h‖ + ‖u_h − u_h^(k)‖ + ‖u_h^(k) − ũ_h^(k)‖

数学离散误差与代数误差：由 PDE 理论和算法理论给出
浮点舍入误差：由浮点分析给出
TLA+ 前提：程序实际算出的就是 ũ_h^(k)，由精化保证

这个不等式是项目的骨架。它的每一项都要有定理支撑，而它能够成立的前提（并行程序忠实地执行了抽象算法）要由软件层的规范来保证。模型本身与真实物理之间的差距不在这个不等式里，那是确认（validation）问题，见 §7。

§3

两层形式化，一座桥，一个地基

数学层

PDE 理论与数值分析

适定性：Lax–Milgram、inf-sup 条件、Kohn–Sham 能量泛函极小元的存在性
离散收敛：先验误差估计（Céa 引理、插值估计、平面波截断误差）与后验误差估计
结构保持：辛性、时间可逆性、动量守恒、电荷守恒、晶体对称性
迭代收敛：SCF 不动点、Newton 法、Krylov 方法的收敛条件与停止准则

工具　Lean 4 + Mathlib。Mathlib 已有 Hilbert 空间、测度论和 Lax–Milgram，可以直接作为起点。策略　先把定理写成带文献出处的公理，让上层可以立刻依赖它；再逐步把公理替换为机器证明。每一条尚未证明的公理都登记在可信账本里。

软件层

TLA+ 规范与模型检验

安全性：原子不丢不重，每个自由度只有一个所有者，并行组装无写冲突，无死锁
活性：消息最终送达；SCF 循环要么收敛，要么在有限步内报告失败
协议：区域分解与 halo 交换、分布式 FFT 转置、检查点与重启、自适应网格加密与负载再平衡

工具　PlusCal 写算法；TLC 在小规模实例上穷举；Apalache 做符号检验；TLAPS 写对任意规模成立的证明。与代码的连接　实现运行时输出事件日志，用 trace validation 检查日志是否是规范允许的行为。

桥

精化：把数学定理传递给并行程序

抽象算法写成 TLA+ 规范 A：状态是全局密度或全局相空间坐标，每一步是一次数学上的迭代映射。并行实现写成规范 P：状态分散在各个进程上。精化映射把 P 的状态投影成 A 的状态，例如把各进程上的局部密度拼成全局密度。证明 P ⇒ A 之后，数学层关于 A 的收敛定理就对并行程序成立。

性能优化也按同样的方式处理：每一次优化都是一次新的精化，必须证明它仍然实现上一层规范。

地基

浮点是一等公民

IEEE 754 加法不满足结合律。浮点结果不直接进入结论：大型求解器的输出交给经过证明的检查器，检查器要么全程用精确有理数，要么用向外舍入的区间算术，所以需要信任的浮点性质缩小到“就近舍入的结果是最近的 double”。

全局归约采用确定性的求和顺序，使结果不随进程数或线程数变化。舍入模型写进规范，与离散误差、代数误差放在同一个误差预算里。

§4

证明义务矩阵

三个领域在六层上各自需要证明的内容。这张表是项目的工作清单：每个单元格最终都要对应一份规范、一个定理或一条登记过的假设。

层	DFT · 密度泛函理论 · 电子结构	MD · 分子动力学 · 原子轨迹	FEM · 有限元 · 连续介质
模型	Kohn–Sham 方程，一个非线性特征值问题：H[ρ]ψ_i = ε_iψ_i，ρ = Σf_i\|ψ_i\|²	哈密顿系统；恒温模拟为 Langevin 动力学及其 Fokker–Planck 方程	椭圆、抛物、双曲型 PDE，例如线弹性 −∇·σ(u) = f 及其弱形式
连续理论	Lieb 变分原理；LDA 型 Kohn–Sham 模型能量极小元的存在性；Hamiltonian 的自伴性	力场 Lipschitz 时解存在唯一；能量、动量、角动量守恒；相空间体积守恒（Liouville）；Langevin 动力学的不变测度与遍历性	Lax–Milgram 给出适定性；混合问题（Stokes、近不可压弹性）需要 inf-sup 条件；解的正则性决定可达的收敛阶
离散化	平面波截断的先验误差估计；Brillouin 区 k 点积分误差；赝势带来的模型误差单独登记	Velocity Verlet 的辛性与时间可逆性；向后误差分析：存在修正哈密顿量，能量误差 O(Δt²) 在指数长的时间内有界	Céa 引理加插值估计：‖u − u_h‖_H¹ ≤ C h^k\|u\|_H^k+1；网格形状正则性作为显式前提；后验估计量的可靠性与有效性
算法	SCF 不动点迭代与 Anderson / Pulay 混合的局部收敛条件；LOBPCG、Davidson 保持轨道正交；SCF 状态机：收敛、振荡、重启	近邻表：只要上次重建后最大位移小于 skin 厚度的一半，就不会漏掉相互作用对；SHAKE / RATTLE 约束迭代的收敛	CG / GMRES 与预条件子的收敛和停止准则；自适应循环 SOLVE → ESTIMATE → MARK → REFINE 的收敛性
并行执行	k 点、能带、平面波三级并行的数据划分；分布式三维 FFT 转置不丢不重；全局归约后各进程看到同一个 ρ	区域分解下原子迁移守恒；ghost 原子交换一致；无死锁；检查点重启前后状态等价	并行组装无写冲突；每个自由度恰有一个所有者；分布式加密后网格保持协调，负载再平衡不丢单元
浮点	正交化过程中舍入误差的累积；电荷守恒 ∫ρ = N 的数值偏差界	长时间积分中舍入误差的累积；归约顺序导致的逐位不可复现	刚度矩阵条件数 κ ∼ h⁻² 对可达精度的限制；数值求积误差

§5

可信等级

不是所有东西都能马上被证明。每个组件都要如实标出自己所处的等级，并在可信账本中列出它依赖的假设。等级只升不虚标。

测试

回归测试、解析解对比、与既有代码对标。只对测过的输入成立。

文献证明

有纸面证明和文献出处，假设已显式列出，但尚未机器检查。

模型检验

TLC 或 Apalache 在有限规模实例上穷举了全部行为。

机器证明

Lean、Rocq 或 TLAPS 证明，对任意规模成立，但停留在抽象模型、精确算术或代码的逐行模型层面。

端到端

机器证明一直覆盖到可执行代码和浮点舍入。先例：Boldo 等人对一维波动方程 C 程序的完整证明。

§6

原则

规范先于代码

每个模块先有数学陈述和 TLA+ 规范，再有实现。

假设必须可见

正则性、网格条件、步长限制、收敛阈值都写进规范，不藏在输入文件的默认值里。

等级如实

宁可标 T0，也不虚标 T3。可信账本是项目里最重要的文件。

小的可信内核

只有少数内核需要达到 T4；上层通过组合定理继承它们的保证。

浮点进入规范

舍入模型和误差界与离散误差一起进入同一个误差预算。

默认可复现

相同输入在任意进程数下给出逐位相同的结果，除非使用者显式放弃这一点。

先正确，后性能

先得到可证明的版本，再在不破坏证明的前提下优化。每次优化都是一次新的精化。

验证检查器，不验证求解器

LAPACK、FFT 库这类大型数值库不做证明，而是为它们的输出写经过证明的检查器。求解器可以替换，保证来自检查器。

每次运行附带报告

报告的每一行都标出可信账本中的命题编号，由编号查到它的等级与所依赖的假设；报告末尾列出本次运行用到的假设。

每个输出都有定理

程序输出的每个数都要对应一条机器证明的定理。结果是近似值时（迭代、离散化、k 点采样、有限温度、浮点），要证明的是误差界，测试、模型检验与 sanitizer 只作辅助证据。还没有这样一条定理的输出算未完成，列为缺口。

开放

规范、证明、基准和可信账本都将公开，任何人都可以重新检查。

§7

不做什么

不判断物理模型是否正确

交换关联泛函、力场参数、本构关系与真实世界的差距属于确认（validation）。本项目只做验证（verification）：正确地求解给定的方程。

不在一开始就替代现有软件

Quantum ESPRESSO、VASP、LAMMPS、GROMACS、deal.II、FEniCS 是对标对象和参照。

不要求所有代码都被证明

输入解析、I/O、可视化停留在 T0 即可，只要它们不进入误差预算。

§8

路线图

当前目标 G1：DFTK 功能对齐 2026-09-29 设定

让 mini-DFT 覆盖 DFTK.jl 的主要功能。只有实现、与 DFTK 在同一离散问题上一致、输出有 Lean 证书三条都满足，才算完成。顺序是先在 1D 做全，再把 Lean 模型推广到 d 维做 3D，最后是基础设施与外部接口。FEM 改由并行目标 G2 推进；阶段 2 的 MD 在 2026-10-01 做出了基础版本 mini-MD v0。逐项状态 →

并行目标 G2：带形式化证明的有限元，对照 scikit-fem 2026-09-30 设定

用自己的 C 实现覆盖 scikit-fem 的主要功能，scikit-fem 在同一离散问题上作独立参照。完成标准与 G1 相同；能做到的地方，证书的结论针对连续问题的真解，而不只是离散解。先在 1D 做全，再做 2D 三角形网格。逐项状态 →

并行目标 G3：带形式化证明的 CALPHAD，对照 pycalphad 2026-10-01 设定

用自己的 C 实现读 TDB 热力学数据库、求相平衡，pycalphad 在同一数据库、同一组相上作独立参照。完成标准与 G1、G2 相同；证书的结论针对数据库模型真正的全局平衡，而不只是驱动找到的局部解。先把二元体系做全（点平衡、相图、不变反应、磁性），再做多元与亚点阵模型。逐项状态 →

阶段 0最小 DFT进行中

一维周期约化 Hartree 模型，平面波离散，串行执行。用它把 L1 到 L6 完整走通一遍，同时建立整个项目会复用的基础设施：可信账本、证书检查器、运行报告格式、trace validation。

退出标准义务清单中的命题全部达到目标等级；四组基准通过；每次运行都输出可信报告。
阶段 1DFT 扩展部分提前完成

一维 LDA 交换关联（凸性丢失，唯一性改为局部结论）；k 点采样，以及 k 点的 MPI 并行——项目第一个真正的并行协议；三维平面波与局域赝势，与 DFTK.jl 交叉校验。

退出标准k 点并行在 1 到 64 个进程上结果逐位一致，运行日志全部通过 trace validation。
阶段 2推广到 MD 与 FEM已达退出标准

复用阶段 0 建立的框架，做两个最小版本：一维 Lennard-Jones 原子链加 Velocity Verlet；一维 Poisson 方程加 P1 线性元。TLA+ 模块覆盖 halo 交换、原子迁移、并行组装、检查点。

退出标准两个最小版本都输出可信报告；Céa 引理与 Verlet 辛性达到 T3。
阶段 3可信内核未开始

选定少数内核达到 T4：平面波动能算子的作用、FEM 单元刚度矩阵组装、Verlet 单步；浮点误差界的机器证明与代码相连。

退出标准这些内核的浮点误差界经过机器检查，并进入全局误差预算。
阶段 4对标未开始

在标准基准上与既有代码比较：DFT 用 Δ 量规，MD 用 NVE 能量漂移和径向分布函数，FEM 用制造解方法检验收敛阶。

退出标准公开基准结果，以及每个结果对应的可信账本。

逐项现状见进度页，以及每个领域各自的页面：DFT、FEM、MD、CALPHAD。

§9

第一个里程碑：mini-DFT 与 mini-FEM

在扩展到 MD 和 FEM 之前，先用一个尽量小的 DFT 程序把 L1 到 L6 全部走通。它要小到每一层的证明义务都能写清楚，同时保留 Kohn–Sham 计算的完整骨架：非线性特征值问题、自洽迭代、平面波离散。

模型	一维周期约化 Hartree（rHF）模型，不含交换关联项。区域为周期区间 [0, L)，N 个电子，自旋简并。能量为动能、外势能与 ½ Hartree 能之和，Hartree 核在非零频率上为 4π/q²。
外势	光滑周期函数：常数、余弦、周期化的高斯双势阱（可重复多份）。一维 Coulomb 势在原点不可积，因此不用它。
离散	平面波，\|k\| ≤ K；只取 Γ 点，或 N_k 个等权 k 点（Γ 中心或 Monkhorst–Pack 网格）。实空间网格点数不少于 4K + 1，保证密度和 Hamiltonian 作用没有混叠。
占据	零温时假设能隙存在，每次迭代都用区间算术检查；有限温度时用 Fermi–Dirac、Gaussian、Methfessel–Paxton 或 Marzari–Vanderbilt 占据，所有 k 点共用一个化学势。
算法	SCF 迭代，线性或 Anderson 混合，可加 Kerker 预条件；稠密矩阵直接对角化；k 点用共享内存线程并行，按固定顺序累加。

为什么选 rHF

它的能量关于密度矩阵是凸的，而 Hartree 项在非零频率上严格凸，所以基态密度唯一。在 DFT 家族里，这是数学性质最干净的模型，适合作为第一个机器证明的对象。加入 LDA 交换关联后凸性会丢失，这一步放在阶段 1。

LAPACK 求解，检查器担保

LAPACK 本身停在 T0。驱动收敛后写出证书，用 Lean 写成并编译的检查器以精确有理数重新检查，输出能量、自由能、截断极限、密度、本征值、能带与态密度的严格区间。检查器的可靠性是 Lean 定理；它不信任求解器给出的任何数，求解器算错只会让检查失败或区间变宽。

mini-FEM 在目标 G2 下沿用同一模式：一维 Poisson 方程的 P1 线性元（Dirichlet、Neumann、Robin 边界），梁方程的三次 Hermite 元（固支、简支、自由、滑动），右端项由多项式、指数、正弦、余弦组合而成。Lean 检查器给出连续问题真解在每个节点上的严格区间，以及整个区间上最大误差的严格上界。自适应加密以这些经过证明的单元误差界作为标记量。

mini-MD（2026-10-01）对一维 Lennard-Jones 原子链加 Velocity Verlet 做了同样的事：检查器给出每一步与 Newton 方程真解的距离，以及从初值出发约一个振动周期内的整体误差界。mini-CALPHAD（目标 G3）为读自 TDB 数据库的二元相平衡出具证书（包括磁性相与整张相图），给出化学势、驱动力，以及任意平衡中可能出现哪些相。详见 MD 与 CALPHAD 页。

§10

决定与待决问题

D1 · 已决定 · 2026-09-27 · 修订 · 2026-09-29

证明器分层：数学层与浮点模型用 Lean 4，代码层连接待选工具

分工

Lean 4 + Mathlib 负责 L1–L4（适定性、离散误差估计、结构保持、迭代收敛），以及 L6 的两部分：证书检查器（精确有理数，编译成可执行程序）与区间算术内核的逐行 binary64 模型。
Frama-C/WP 负责 C 可信内核的运行时安全。
把 Lean 中的浮点结论落到实际的 C 代码上，候选工具是 Rocq 的 Flocq、VST、CompCert，或 Frama-C 的浮点模型。尚未开始。

理由

数学层是项目的主体；Mathlib 的泛函分析与测度论是一个统一的库，AI 辅助证明工具也大多优先支持 Lean。“验证检查器，不验证求解器”把浮点义务缩小了，Lean 中一个很小的 binary64 模型就够用。代码级验证方面 Rocq 的工具链和 T4 先例仍然更成熟，所以这一步的工具尚未选定。

代价

Lean 模型与 C 代码的对应目前只有逐行对照与逐位差分测试，没有机器证明，区间算术的正确性停在 T3。

D2 · 已决定 · 2026-09-27

实现语言：C11，分成经过证明的内核和不受信任的驱动