原子计算示例

本节提供常见原子的详细计算示例和结果对比。

氢 (H, Z=1)

单电子体系，HF 与精确解相同。

配置

from atomscf.grid import radial_grid_linear
from atomscf.scf_hf import run_hf_minimal

r, w = radial_grid_linear(n=800, rmin=1e-6, rmax=50.0)
result = run_hf_minimal(Z=1, r=r, w=w)

结果

性质	数值	理论值
E_total	-0.500 Ha	-0.500 Ha（精确）
ε_1s	-0.500 Ha	-0.500 Ha
误差	< 0.001%

氦 (He, Z=2)

闭壳层，RHF 适用。

RHF 计算

from atomscf.scf_hf import run_hf_scf, HFSCFGeneralConfig

config = HFSCFGeneralConfig(
    Z=2, r=r, w=w,
    occ_by_l={0: [2.0]},
    eigs_per_l={0: 1},
    spin_mode='RHF',
    mix_alpha=0.5,
    tol=1e-6,
    maxiter=100
)
result = run_hf_scf(config)

结果对比

方法	E_total (Ha)	ε_1s (Ha)	相对误差
本代码 (RHF)	-2.787	-0.866	4.0%
Clementi (HF)	-2.862	-0.918	1.4%
实验	-2.904

备注

RHF 缺少关联能（~0.04 Ha）。使用 LSDA 可改善至 ~1% 误差。

锂 (Li, Z=3)

开壳层 (1s² 2s¹)，UHF 优于 RHF。

UHF 计算

config = HFSCFGeneralConfig(
    Z=3, r=r, w=w,
    occ_by_l={0: [2.0, 1.0]},
    occ_by_l_spin={
        0: {'up': [1.0, 1.0], 'down': [1.0, 0.0]}
    },
    eigs_per_l={0: 2},
    spin_mode='UHF',
    mix_alpha=0.3,
    tol=1e-6,
    maxiter=120
)
result = run_hf_scf(config)

结果对比

方法	E_total (Ha)	ε_2s (Ha)	误差 (mHa)
RHF	-7.159	-0.089	274
UHF	-7.213	-0.193 (↑)	219
Clementi (ROHF)	-7.433	-0.196

改善: UHF 比 RHF 降低误差 20%

碳 (C, Z=6)

开壳层 p² (³P 基态)，DFT 推荐。

LSDA-VWN 计算

from atomscf.scf import run_lsda_vwn, SCFConfig

config = SCFConfig(
    Z=6, r=r, w=w,
    lmax=2,
    eigs_per_l=2,
    eig_solver="fd5_aux",
    xc="VWN",
    mix_alpha=0.5,
    tol=5e-5,
    maxiter=140
)
result = run_lsda_vwn(config, verbose=True)

结果对比

性质	LSDA-VWN	Clementi (ROHF)	误差
E_total	-37.83 Ha	-37.69 Ha	+140 mHa
ε_2s	-1.02 Ha	-1.02 Ha	~0
ε_2p	-0.56 Ha	-0.43 Ha	+130 mHa

警告

碳原子 HF 计算需要 ROHF（非 UHF）以正确描述 ³P 态。 LSDA 包含关联，结果更合理。

氮 (N, Z=7)

开壳层 p³ (⁴S 基态)，高自旋。

LSDA-PZ81 计算

from atomscf.scf import run_lsda_pz81

config = SCFConfig(
    Z=7, r=r, w=w,
    lmax=2,
    eigs_per_l=2,
    xc="PZ81",  # PZ81 关联
    mix_alpha=0.4,
    tol=5e-5,
    maxiter=150
)
result = run_lsda_pz81(config)

结果

性质	LSDA-PZ81	文献值
E_total	-54.4 Ha	-54.6 Ha
ε_2p	-0.67 Ha	-0.54 Ha

氧 (O, Z=8)

开壳层 p⁴ (³P 基态)。

LSDA-VWN 计算

config = SCFConfig(
    Z=8, r=r, w=w,
    lmax=2,
    eigs_per_l=2,
    xc="VWN",
    mix_alpha=0.3,  # 慢收敛，小混合
    tol=5e-5,
    maxiter=180
)
result = run_lsda_vwn(config, progress_every=20)

收敛技巧

氧原子收敛较难，建议：

减小 mix_alpha 至 0.2-0.3
增加 maxiter 至 200+
放宽 tol 至 1e-4

氖 (Ne, Z=10)

闭壳层 (1s² 2s² 2p⁶)，RHF 适用。

RHF 计算

config = HFSCFGeneralConfig(
    Z=10, r=r, w=w,
    occ_by_l={0: [2.0, 2.0], 1: [6.0]},
    eigs_per_l={0: 2, 1: 1},
    spin_mode='RHF',
    mix_alpha=0.4,
    tol=1e-6,
    maxiter=150
)
result = run_hf_scf(config)

结果对比

性质	RHF	Clementi (HF)
E_total	-128.3 Ha	-128.5 Ha
ε_1s	-32.6 Ha	-32.8 Ha
ε_2s	-1.90 Ha	-1.93 Ha
ε_2p	-0.82 Ha	-0.85 Ha

批量计算脚本

计算多个原子：

from atomscf.grid import radial_grid_linear
from atomscf.scf import run_lsda_vwn, SCFConfig

# 原子列表
atoms = {
    'H': 1, 'He': 2, 'Li': 3, 'Be': 4,
    'B': 5, 'C': 6, 'N': 7, 'O': 8,
    'F': 9, 'Ne': 10
}

r, w = radial_grid_linear(n=1200, rmin=1e-6, rmax=70.0)
results = {}

for name, Z in atoms.items():
    print(f"\n计算 {name} (Z={Z})...")

    config = SCFConfig(
        Z=Z, r=r, w=w,
        lmax=2,
        eigs_per_l=2,
        xc="VWN",
        mix_alpha=0.3 if Z > 5 else 0.5,
        tol=5e-5,
        maxiter=200
    )

    try:
        result = run_lsda_vwn(config, verbose=False)
        results[name] = {
            'E': result.energies['E_total'],
            'converged': result.converged,
            'iters': result.iterations
        }
        print(f"  ✓ E = {result.energies['E_total']:.4f} Ha")
    except Exception as e:
        print(f"  ✗ 失败: {e}")

# 保存结果
import json
with open('results.json', 'w') as f:
    json.dump(results, f, indent=2)

网格收敛性测试

测试不同网格点数的影响：

import numpy as np
import matplotlib.pyplot as plt

grid_sizes = [400, 600, 800, 1000, 1200, 1500]
energies = []

for n in grid_sizes:
    r, w = radial_grid_linear(n=n, rmin=1e-6, rmax=50.0)
    result = run_hf_minimal(Z=2, r=r, w=w)  # He 原子
    energies.append(result.E_total)
    print(f"n={n:4d}: E = {result.E_total:.8f} Ha")

# 绘图
plt.figure(figsize=(8, 5))
plt.plot(grid_sizes, energies, 'o-')
plt.axhline(y=-2.8617, color='r', linestyle='--', label='HF 极限')
plt.xlabel('网格点数 n')
plt.ylabel('总能量 (Ha)')
plt.title('He 原子网格收敛性')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.savefig('grid_convergence.png', dpi=150)

预期趋势：能量单调收敛至 HF 极限。

参考数据来源

Clementi & Roetti (1974): HF 极限值（Z ≤ 54）
NIST Atomic Spectra Database: 实验光谱数据
Computational Chemistry Comparison and Benchmark Database: 多种方法对比

常见问题

为什么我的结果与文献不同？

可能原因：

方法差异: 本代码 RHF/UHF，文献可能 ROHF
基组差异: 网格离散 vs Gaussian 基组
关联缺失: HF 无关联，DFT 有近似关联
收敛不足: 检查 converged 标志

如何选择 HF vs DFT？

HF (RHF/UHF): - 优势：交换精确 - 适用：闭壳层、教学、基准 - 缺点：无关联
DFT (LSDA): - 优势：速度快、包含关联 - 适用：开壳层、过渡金属 - 缺点：交换近似、带隙低估

建议: 教学用 HF，实际计算用 DFT。

下一步

Hartree-Fock 方法: HF 理论深入
API 参考: 高级 API 使用