pyscf

Name: pyscf
Author: tondevrel

github.com/tondevrel其他2458 次安装38 周前更新

代码仓库

关于 pyscf

PySCF 是用于从头算电子结构与量子化学的 Python 框架，适合研究分子及周期性体系的科研人员。它涵盖 Hartree-Fock、DFT、后 Hartree-Fock 方法、几何优化、激发态和分子性质计算。

主要功能

01Hartree-Fock 与 DFT 计算
02MP2、CCSD 与 CI 方法
03分子几何优化
04激发态计算
05分子性质分析
06周期性体系计算

适用场景

01计算分子能量与梯度
02优化分子几何结构
03分析激发态与分子轨道
04模拟固体与表面

安装

复制下面的命令，在终端中运行。

$npx skills add https://github.com/tondevrel/scientific-agent-skills --skill pyscf

安装后在 Claude Code 中输入/pyscf确认已加载。

详细介绍

PySCF - 量子化学框架

用于从头算电子结构计算和量子化学的 Python 库。

适用场景

运行 Hartree-Fock（HF）和密度泛函理论（DFT）计算
计算分子能量和梯度
优化分子几何结构
计算后 HF 方法（MP2、CCSD、CI、耦合簇）
分析分子轨道和电子密度
计算激发态（TD-DFT、CI、EOM-CCSD）
计算分子性质（偶极矩、电荷、NMR、IR）
执行周期性计算（固体、表面）
对量子化学方法进行基准测试
将量子计算整合到 ML/AI 工作流中

参考文档

官方文档：https://pyscf.org/
搜索模式：gto.M、scf.RHF、dft.RKS、mp.MP2、optimize、tddft.TDDFT

核心原则

使用 PySCF 执行

任务	模块	示例
构建分子	`gto`	`gto.M(atom='H 0 0 0; H 0 0 1')`
Hartree-Fock	`scf`	`scf.RHF(mol).run()`
DFT 计算	`dft`	`dft.RKS(mol, xc='B3LYP')`
MP2 相关	`mp`	`mp.MP2(mf).run()`
耦合簇	`cc`	`cc.CCSD(mf).run()`
几何优化	`geomopt`	`optimize(mf)`
激发态	`tddft`	`tddft.TDDFT(mf).run()`
周期性体系	`pbc`	`pbc.gto.Cell()`

不要用于

分子动力学模拟（使用 GROMACS、OpenMM、ASE）
超大体系（>1000 个原子）——使用半经验方法或力场
交互式可视化（使用 PyMOL、VMD）
高通量虚拟筛选（速度太慢）
实时量子模拟

快速参考

安装

# pip (recommended)
pip install pyscf

# With extensions
pip install pyscf[geomopt,dftd3,dmrgscf]

# From conda
conda install -c pyscf pyscf

# Development version
pip install git+https://github.com/pyscf/pyscf

标准导入

# Core modules
from pyscf import gto, scf, dft
from pyscf import mp, cc, ci
from pyscf import grad, geomopt
from pyscf import tddft, tdscf
from pyscf import lo, ao2mo

# Tools
from pyscf.tools import molden, cubegen
import numpy as np

基本模式 - 单点能量

from pyscf import gto, scf

# 1. Build molecule
mol = gto.M(
    atom='O 0 0 0; H 0 1 0; H 0 0 1',
    basis='6-31g',
    charge=0,
    spin=0  # 2S, 0 = singlet
)

# 2. Run SCF
mf = scf.RHF(mol)
energy = mf.kernel()

# 3. Check convergence
if not mf.converged:
    raise RuntimeError("SCF did not converge")

print(f"E(HF) = {energy:.8f} Hartree")

基本模式 - DFT 计算

from pyscf import gto, dft

mol = gto.M(atom='C 0 0 0; O 0 0 1.2', basis='def2-tzvp')

# DFT with B3LYP
mf = dft.RKS(mol)
mf.xc = 'B3LYP'
energy = mf.kernel()

print(f"E(B3LYP) = {energy:.8f} Hartree")

基本模式 - 几何优化

from pyscf import gto, scf
from pyscf.geomopt.geometric_solver import optimize

mol = gto.M(atom='H 0 0 0; H 0 0 1.5', basis='6-31g')
mf = scf.RHF(mol)

# Optimize geometry
mol_eq = optimize(mf)

print(f"Optimized geometry:\n{mol_eq.atom}")

关键规则

✅ 应该做

检查收敛性 - 始终确认 mf.converged == True
从小基组开始 - 先用 STO-3G 或 3-21G 测试
指定电荷和自旋 - 明确电子态
保存结果 - 使用 chkfile 以支持重启
利用对称性 - 适用时启用以提升速度
提供良好的初始猜测 - 使用 init_guess='minao' 或先前的结果
监控内存 - 适当设置 max_memory
确认单位 - PySCF 内部使用 Bohr 和 Hartree
使用密度拟合 - 在大型体系中启用
检查多重度 - 确保自旋与预期态一致

❌ 不要做

跳过收敛性检查 - 绝不要使用未收敛计算的结果
一开始就使用超大基组 - 从小基组开始，逐步扩大
忽略电荷/自旋 - 默认值可能与体系不匹配
运行时不使用 chkfile - 崩溃时会丢失数据
不必要地禁用对称性 - 对称性可以加快计算
混用坐标系 - 统一使用埃或 Bohr
使用错误的多重度 - 应为 2S+1，而不是 2S
忘记考虑内存 - 大型计算可能导致 OOM
相信未收敛的结果 - 这些结果没有意义
比较使用不同基组的方法 - 使用一致的基组

反模式（绝不要）

# ❌ BAD: Not checking convergence
mf = scf.RHF(mol)
energy = mf.kernel()
# Use energy without checking convergence!

# ✅ GOOD: Always check convergence
mf = scf.RHF(mol)
energy = mf.kernel()
if not mf.converged:
    raise RuntimeError("SCF not converged, results unreliable")

# ❌ BAD: Using huge basis immediately
mol = gto.M(atom='protein.xyz', basis='def2-qzvppd')
# This will take forever or crash!

# ✅ GOOD: Start small, expand if needed
mol = gto.M(atom='protein.xyz', basis='sto-3g')
# Test first, then use def2-svp, def2-tzvp, etc.

# ❌ BAD: Ignoring electronic state
mol = gto.M(atom='O 0 0 0; O 0 0 1.2')
mf = scf.RHF(mol)  # What charge? What spin?

# ✅ GOOD: Explicit state specification
mol = gto.M(
    atom='O 0 0 0; O 0 0 1.2',
    charge=0,
    spin=2  # Triplet O2
)
mf = scf.ROHF(mol)  # Restricted open-shell for triplet

# ❌ BAD: No checkpoint file
mf = scf.RHF(mol)
mf.kernel()  # If crash, lose everything!

# ✅ GOOD: Use checkpoint files
mf = scf.RHF(mol)
mf.chkfile = 'calculation.chk'
mf.kernel()
# Can restart if needed

# ❌ BAD: Wrong multiplicity
mol = gto.M(atom='O 0 0 0; O 0 0 1.2', spin=2)  # Wrong!
# spin should be 2S, not multiplicity!

# ✅ GOOD: Correct spin specification
mol = gto.M(
    atom='O 0 0 0; O 0 0 1.2',
    spin=2  # 2S = 2 for triplet (multiplicity = 2S+1 = 3)
)

构建分子

分子定义基础

from pyscf import gto

# Cartesian coordinates (Angstrom by default)
mol = gto.M(
    atom='''
        O  0.0  0.0  0.0
        H  0.0  1.0  0.0
        H  0.0  0.0  1.0
    ''',
    basis='6-31g',
    charge=0,
    spin=0
)

# Compact format
mol = gto.M(
    atom='O 0 0 0; H 0 1 0; H 0 0 1',
    basis='6-31g'
)

# Using tuple format
mol = gto.M(
    atom=[
        ('O', (0.0, 0.0, 0.0)),
        ('H', (0.0, 1.0, 0.0)),
        ('H', (0.0, 0.0, 1.0))
    ],
    basis='6-31g'
)

从文件加载

from pyscf import gto

# From XYZ file
mol = gto.M(atom='molecule.xyz', basis='6-31g')

# From PDB file
mol = gto.M(atom='protein.pdb', basis='sto-3g')

# From Z-matrix
mol = gto.M(
    atom='''
    H
    H  1  0.74
    ''',
    basis='6-31g'
)

分子性质和信息

from pyscf import gto

mol = gto.M(atom='H2O', basis='6-31g')

# Basic properties
print(f"Number of atoms: {mol.natm}")
print(f"Number of electrons: {mol.nelectron}")
print(f"Number of basis functions: {mol.nao_nr()}")
print(f"Charge: {mol.charge}")
print(f"Spin (2S): {mol.spin}")
print(f"Nuclear repulsion: {mol.energy_nuc():.8f}")

# Atom information
for i in range(mol.natm):
    atom_symbol = mol.atom_symbol(i)
    atom_charge = mol.atom_charge(i)
    coord = mol.atom_coord(i)
    print(f"Atom {i}: {atom_symbol} (Z={atom_charge}) at {coord}")

基组指定

from pyscf import gto

# Single basis for all atoms
mol = gto.M(atom='H2O', basis='6-31g')

# Different basis for different atoms
mol = gto.M(
    atom='H2O',
    basis={
        'O': '6-311g**',
        'H': '6-31g'
    }
)

# Custom basis sets
mol = gto.M(
    atom='H 0 0 0; H 0 0 1',
    basis={
        'H': gto.basis.parse('''
            H    S
                 13.0107010              0.19682158E-01
                  1.9622572              0.13796524
                  0.44453796             0.47831935
            H    S
                  0.12194962             1.0000000
        ''')
    }
)

分子对称性

from pyscf import gto

# Auto-detect symmetry
mol = gto.M(atom='H2O', basis='6-31g', symmetry=True)
print(f"Point group: {mol.groupname}")
print(f"Irreps: {mol.irrep_name}")

# Specify symmetry
mol = gto.M(atom='H2O', basis='6-31g', symmetry='C2v')

# Disable symmetry
mol = gto.M(atom='H2O', basis='6-31g', symmetry=False)

Hartree-Fock 方法

限制性 Hartree-Fock（RHF）

from pyscf import gto, scf

# Closed-shell molecule
mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)

# Run calculation
energy = mf.kernel()

if mf.converged:
    print(f"E(RHF) = {energy:.8f} Hartree")
    print(f"Orbital energies: {mf.mo_energy}")
else:
    print("SCF not converged")

非限制性 Hartree-Fock（UHF）

from pyscf import gto, scf

# Open-shell molecule (triplet O2)
mol = gto.M(
    atom='O 0 0 0; O 0 0 1.2',
    basis='6-31g',
    spin=2  # Triplet
)

mf = scf.UHF(mol)
energy = mf.kernel()

# Analyze spin
print(f"<S^2> = {mf.spin_square()[0]:.4f}")  # Should be ~2.0 for triplet

限制性开壳层 HF（ROHF）

from pyscf import gto, scf

# Open-shell with restricted spatial orbitals
mol = gto.M(
    atom='O 0 0 0; O 0 0 1.2',
    basis='6-31g',
    spin=2
)

mf = scf.ROHF(mol)
energy = mf.kernel()

print(f"E(ROHF) = {energy:.8f} Hartree")

SCF 收敛控制

from pyscf import gto, scf

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)

# Convergence parameters
mf.conv_tol = 1e-8  # Energy convergence threshold
mf.conv_tol_grad = 1e-5  # Gradient convergence
mf.max_cycle = 100  # Maximum iterations
mf.diis_space = 12  # DIIS history size

# Initial guess
mf.init_guess = 'minao'  # Options: 'minao', 'atom', '1e', 'huckel'

# Level shift (helps difficult cases)
mf.level_shift = 0.5

energy = mf.kernel()

困难情况的收敛策略

from pyscf import gto, scf

mol = gto.M(atom='difficult_molecule.xyz', basis='6-31g')

# Strategy 1: Smaller basis first
mol_small = gto.M(atom='difficult_molecule.xyz', basis='sto-3g')
mf_small = scf.RHF(mol_small)
mf_small.kernel()

# Use small basis result as initial guess
mol_large = gto.M(atom='difficult_molecule.xyz', basis='6-31g')
mf = scf.RHF(mol_large)
mf.init_guess = mf_small.make_rdm1()  # Use small basis density

# Strategy 2: Use level shift
mf.level_shift = 0.3
mf.kernel()

# Strategy 3: Try different initial guess
if not mf.converged:
    mf.init_guess = 'atom'
    mf.kernel()

密度泛函理论

DFT 计算基础

from pyscf import gto, dft

mol = gto.M(atom='H2O', basis='def2-tzvp')

# Setup DFT
mf = dft.RKS(mol)
mf.xc = 'B3LYP'  # Hybrid functional

energy = mf.kernel()
print(f"E(B3LYP) = {energy:.8f} Hartree")

常用泛函

from pyscf import gto, dft

mol = gto.M(atom='CH4', basis='6-31g')

# LDA
mf = dft.RKS(mol)
mf.xc = 'LDA,VWN'
e_lda = mf.kernel()

# GGA
mf.xc = 'PBE'
e_pbe = mf.kernel()

# Hybrid GGA
mf.xc = 'B3LYP'
e_b3lyp = mf.kernel()

# Meta-GGA
mf.xc = 'M06'
e_m06 = mf.kernel()

# Range-separated
mf.xc = 'wB97X-D'
e_wb97xd = mf.kernel()

print(f"LDA:    {e_lda:.6f}")
print(f"PBE:    {e_pbe:.6f}")
print(f"B3LYP:  {e_b3lyp:.6f}")
print(f"M06:    {e_m06:.6f}")
print(f"wB97X-D: {e_wb97xd:.6f}")

非限制性 DFT

from pyscf import gto, dft

# Radical species
mol = gto.M(
    atom='C 0 0 0; H 0 0 1; H 0 1 0; H 1 0 0',  # CH3 radical
    basis='6-31g*',
    spin=1  # Doublet
)

mf = dft.UKS(mol)
mf.xc = 'B3LYP'
energy = mf.kernel()

print(f"<S^2> = {mf.spin_square()[0]:.4f}")  # Should be ~0.75 for doublet

网格和积分

from pyscf import gto, dft

mol = gto.M(atom='H2O', basis='6-31g')
mf = dft.RKS(mol)
mf.xc = 'B3LYP'

# Set grid quality
mf.grids.level = 3  # 0=coarse, 3=fine, 5=very fine
mf.grids.prune = True  # Prune grid points

# Custom grid
from pyscf.dft import gen_grid
mf.grids = gen_grid.Grids(mol)
mf.grids.atom_grid = {"H": (50, 194), "O": (75, 302)}  # (radial, angular)
mf.grids.build()

energy = mf.kernel()

色散修正

from pyscf import gto, dft

mol = gto.M(atom='benzene.xyz', basis='def2-tzvp')
mf = dft.RKS(mol)

# DFT-D3 dispersion
mf.xc = 'B3LYP'
mf._numint.libxc = dft.xcfun  # Use XCFun for better D3

# Or use PySCF's D3 interface
from pyscf import dftd3
mf = dftd3.dftd3(mf)

energy = mf.kernel()

后 Hartree-Fock 方法

MP2 - 二阶 Møller-Plesset 方法

from pyscf import gto, scf, mp

mol = gto.M(atom='H2O', basis='cc-pvdz')

# First run HF
mf = scf.RHF(mol)
mf.kernel()

# MP2 calculation
mp2 = mp.MP2(mf)
e_corr, t2 = mp2.kernel()

print(f"E(HF)  = {mf.e_tot:.8f}")
print(f"E(MP2) = {mp2.e_tot:.8f}")
print(f"E_corr = {e_corr:.8f}")

CCSD - 耦合簇单激发与双激发

from pyscf import gto, scf, cc

mol = gto.M(atom='H2O', basis='cc-pvdz')

# HF reference
mf = scf.RHF(mol)
mf.kernel()

# CCSD
ccsd = cc.CCSD(mf)
e_corr, t1, t2 = ccsd.kernel()

print(f"E(HF)   = {mf.e_tot:.8f}")
print(f"E(CCSD) = {ccsd.e_tot:.8f}")
print(f"E_corr  = {e_corr:.8f}")

CCSD(T) - 含微扰三激发

from pyscf import gto, scf, cc

mol = gto.M(atom='N2', basis='cc-pvdz')
mf = scf.RHF(mol)
mf.kernel()

# CCSD(T)
ccsd = cc.CCSD(mf)
ccsd.kernel()

# Perturbative triples correction
et = ccsd.ccsd_t()

print(f"E(HF)     = {mf.e_tot:.8f}")
print(f"E(CCSD)   = {ccsd.e_tot:.8f}")
print(f"E(T)      = {et:.8f}")
print(f"E(CCSD(T))= {ccsd.e_tot + et:.8f}")

Configuration Interaction（CI）

from pyscf import gto, scf, ci

mol = gto.M(atom='LiH', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# CISD
cisd = ci.CISD(mf)
e_cisd = cisd.kernel()[0]

# Full CI (only for small systems!)
fci = ci.FCI(mf)
e_fci = fci.kernel()[0]

print(f"E(HF)   = {mf.e_tot:.8f}")
print(f"E(CISD) = {e_cisd:.8f}")
print(f"E(FCI)  = {e_fci:.8f}")

CASSCF - 完全活性空间

from pyscf import gto, scf, mcscf

mol = gto.M(atom='O 0 0 0; O 0 0 1.2', basis='6-31g', spin=2)
mf = scf.ROHF(mol)
mf.kernel()

# CASSCF(6,6) - 6 electrons in 6 orbitals
mc = mcscf.CASSCF(mf, 6, 6)
mc.kernel()

print(f"E(CASSCF) = {mc.e_tot:.8f}")

# Analyze active space
from pyscf import tools
tools.molden.from_mo(mol, 'cas_orbitals.molden', mc.mo_coeff)

几何优化

基本优化

from pyscf import gto, scf
from pyscf.geomopt.geometric_solver import optimize

# Start with non-equilibrium geometry
mol = gto.M(
    atom='H 0 0 0; H 0 0 1.5',  # Too far apart
    basis='6-31g'
)

mf = scf.RHF(mol)

# Optimize
mol_eq = optimize(mf)

print("Optimized geometry:")
print(mol_eq.atom)
print(f"Final energy: {mf.e_tot:.8f}")

带约束的优化

from pyscf import gto, scf
from pyscf.geomopt.geometric_solver import optimize

mol = gto.M(atom='H2O2', basis='6-31g')
mf = scf.RHF(mol)

# Constrained optimization (fix O-O distance)
def callback(mol):
    """Keep O-O distance fixed."""
    coords = mol.atom_coords()
    distance = np.linalg.norm(coords[0] - coords[1])
    return distance - 1.5  # Target distance

mol_eq = optimize(mf, constraints=callback)

过渡态优化

from pyscf import gto, scf
from pyscf.geomopt.berny_solver import optimize

mol = gto.M(atom='transition_state_guess.xyz', basis='6-31g')
mf = scf.RHF(mol)

# TS optimization (maximize along one mode)
mol_ts = optimize(mf, transition_state=True)

print(f"TS energy: {mf.e_tot:.8f}")

反应路径跟踪

from pyscf import gto, scf
from pyscf.geomopt.geometric_solver import optimize

# Optimize reactant
mol_reactant = gto.M(atom='reactant.xyz', basis='6-31g')
mf_reactant = scf.RHF(mol_reactant)
mol_reactant_opt = optimize(mf_reactant)

# Optimize product
mol_product = gto.M(atom='product.xyz', basis='6-31g')
mf_product = scf.RHF(mol_product)
mol_product_opt = optimize(mf_product)

# NEB or IRC for path
# (requires additional libraries)

分子性质

偶极矩

from pyscf import gto, scf

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Dipole moment
dip = mf.dip_moment(unit='Debye')
print(f"Dipole: {dip} Debye")
print(f"|μ| = {np.linalg.norm(dip):.4f} Debye")

Mulliken 布居分析

from pyscf import gto, scf
from pyscf.tools import dump_mat

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Mulliken charges
mulliken = mf.mulliken_pop()
print("\nMulliken charges:")
for i, charge in enumerate(mulliken[1]):
    print(f"Atom {mol.atom_symbol(i)}: {charge:.4f}")

Löwdin 布居分析

from pyscf import gto, scf
from pyscf import lo

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Löwdin analysis
lowdin_charges = lo.orth.lowdin(mol, mf.make_rdm1())
print("Löwdin charges:", lowdin_charges)

自然键轨道（NBO）

from pyscf import gto, scf
from pyscf import nao

mol = gto.M(atom='CH3OH', basis='6-31g*')
mf = scf.RHF(mol)
mf.kernel()

# NBO analysis (requires pyscf-nao)
# nao.analyze(mf)

分子轨道分析

from pyscf import gto, scf

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Orbital energies
print("Orbital energies (eV):")
for i, e in enumerate(mf.mo_energy * 27.2114):  # Convert to eV
    occ = "occ" if i < mol.nelectron//2 else "vir"
    print(f"  MO {i}: {e:8.4f} eV ({occ})")

# HOMO-LUMO gap
homo_idx = mol.nelectron // 2 - 1
lumo_idx = homo_idx + 1
gap = (mf.mo_energy[lumo_idx] - mf.mo_energy[homo_idx]) * 27.2114
print(f"\nHOMO-LUMO gap: {gap:.4f} eV")

电子密度分析

from pyscf import gto, scf
from pyscf.tools import cubegen

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Generate cube file for electron density
cubegen.density(mol, 'h2o_density.cube', mf.make_rdm1())

# Generate cube for specific orbital
cubegen.orbital(mol, 'h2o_homo.cube', mf.mo_coeff[:, homo_idx])

激发态

激发态 TD-DFT

from pyscf import gto, dft, tddft

mol = gto.M(atom='H2O', basis='6-31g')

# Ground state DFT
mf = dft.RKS(mol)
mf.xc = 'B3LYP'
mf.kernel()

# TD-DFT for excited states
td = tddft.TDDFT(mf)
td.nroots = 5  # Number of excited states

excitations = td.kernel()[0]

print("Excitation energies:")
for i, e in enumerate(excitations):
    print(f"State {i+1}: {e * 27.2114:.4f} eV")

含振子强度的 TD-DFT

from pyscf import gto, dft, tddft

mol = gto.M(atom='formaldehyde', basis='6-31g')
mf = dft.RKS(mol)
mf.xc = 'B3LYP'
mf.kernel()

td = tddft.TDDFT(mf)
td.nroots = 10

energies, xy = td.kernel()

# Calculate oscillator strengths
f = td.oscillator_strength()

print("\nExcitation   Energy (eV)   f")
print("-" * 40)
for i in range(len(energies)):
    print(f"  {i+1:2d}       {energies[i]*27.2114:8.4f}    {f[i]:.6f}")

激发态 EOM-CCSD

from pyscf import gto, scf, cc

mol = gto.M(atom='H2O', basis='cc-pvdz')
mf = scf.RHF(mol)
mf.kernel()

# EOM-CCSD
eom = cc.eom_rccsd.EOMIP(cc.CCSD(mf))
e_ip, v_ip = eom.kernel(nroots=3)

print("Ionization potentials:")
for i, e in enumerate(e_ip):
    print(f"IP {i+1}: {e * 27.2114:.4f} eV")

多参考态 CASSCF/CASPT2

from pyscf import gto, scf, mcscf

mol = gto.M(atom='ethylene', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# CASSCF for π → π* states
mc = mcscf.CASSCF(mf, 2, 2)  # 2e in 2 orbitals (π bonding/antibonding)

# State-averaged CASSCF (for excited states)
mc = mc.state_average_([0.5, 0.5])  # Equal weight S0 and S1
mc.kernel()

print("State-averaged energies:")
for i, e in enumerate(mc.e_states):
    print(f"State {i}: {e:.8f}")

高级工作流

势能面扫描

from pyscf import gto, scf
import numpy as np

def calculate_pes_scan(atom_template, distances, basis='6-31g'):
    """Scan potential energy surface along distance."""
    energies = []
    
    for d in distances:
        # Build molecule with current distance
        atom = atom_template.format(d=d)
        mol = gto.M(atom=atom, basis=basis)
        
        mf = scf.RHF(mol)
        energy = mf.kernel()
        
        if mf.converged:
            energies.append(energy)
        else:
            energies.append(np.nan)
    
    return np.array(energies)

# Example: H2 dissociation curve
distances = np.linspace(0.5, 3.0, 20)
atom_template = 'H 0 0 0; H 0 0 {d}'

energies = calculate_pes_scan(atom_template, distances)

# Find minimum
min_idx = np.nanargmin(energies)
print(f"Equilibrium distance: {distances[min_idx]:.3f} Angstrom")
print(f"Energy at minimum: {energies[min_idx]:.8f} Hartree")

热化学

from pyscf import gto, scf, hessian
import numpy as np

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Calculate Hessian (force constants)
h = hessian.RHF(mf)
hess = h.kernel()

# Vibrational frequencies
from pyscf.hessian import thermo
freq_info = thermo.harmonic_analysis(mol, hess)

print("Vibrational frequencies (cm^-1):")
print(freq_info['freq_wavenumber'])

# Thermochemistry at 298.15 K
results = thermo.thermo(mf, freq_info['freq_wavenumber'], 298.15)
print(f"\nZero-point energy: {results['ZPE']:.6f} Hartree")
print(f"Enthalpy: {results['H']:.6f} Hartree")
print(f"Entropy: {results['S']:.6f} Hartree/K")
print(f"Gibbs free energy: {results['G']:.6f} Hartree")

NMR 化学位移

from pyscf import gto, scf
from pyscf.prop import nmr

mol = gto.M(atom='CH4', basis='6-31g*')
mf = scf.RHF(mol)
mf.kernel()

# NMR shielding tensors
nmr_calc = nmr.RHF(mf)
shielding = nmr_calc.kernel()

print("NMR shielding tensors:")
for i in range(mol.natm):
    print(f"{mol.atom_symbol(i)}: {shielding[i]} ppm")

IR 光谱计算

from pyscf import gto, scf, hessian
from pyscf.hessian import thermo

mol = gto.M(atom='CO2', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Calculate Hessian
h = hessian.RHF(mf)
hess = h.kernel()

# Get frequencies and intensities
freq_info = thermo.harmonic_analysis(mol, hess)

print("IR Frequencies and Intensities:")
print("Freq (cm^-1)    Intensity")
for freq, intensity in zip(freq_info['freq_wavenumber'], 
                           freq_info['IR_intensity']):
    if freq > 0:  # Skip imaginary frequencies
        print(f"{freq:8.2f}      {intensity:.4f}")

QMMM 计算

from pyscf import gto, scf, qmmm

# QM region
mol = gto.M(
    atom='C 0 0 0; H 0 1 0; H 1 0 0; H 0 0 1; H -1 0 0',
    basis='6-31g'
)

# MM point charges (environment)
mm_coords = np.array([
    [3.0, 0.0, 0.0],
    [0.0, 3.0, 0.0],
    [0.0, 0.0, 3.0]
])
mm_charges = np.array([-0.5, -0.5, 1.0])

# Add MM charges to QM calculation
mf = scf.RHF(mol)
mf = qmmm.mm_charge(mf, mm_coords, mm_charges)

energy = mf.kernel()
print(f"QM/MM energy: {energy:.8f}")

溶剂化模型（PCM）

from pyscf import gto, scf, solvent

mol = gto.M(atom='H2O', basis='6-31g')

# Polarizable Continuum Model
mf = scf.RHF(mol)
mf = solvent.ddCOSMO(mf)  # ddCOSMO solvation
mf.with_solvent.eps = 78.3553  # Water dielectric

energy = mf.kernel()
print(f"Energy in solution: {energy:.8f}")

周期性体系

基本周期性计算

from pyscf.pbc import gto, scf

# Define unit cell
cell = gto.Cell()
cell.atom = '''
    C 0 0 0
    C 1.68 1.68 1.68
'''
cell.basis = 'gth-dzvp'
cell.pseudo = 'gth-pade'
cell.a = '''
    0.0 3.37 3.37
    3.37 0.0 3.37
    3.37 3.37 0.0
'''
cell.unit = 'angstrom'
cell.build()

# k-points
kpts = cell.make_kpts([2, 2, 2])  # 2x2x2 k-point mesh

# Run calculation
kmf = scf.KRHF(cell, kpts)
energy = kmf.kernel()

print(f"Energy per unit cell: {energy:.8f}")

能带结构计算

from pyscf.pbc import gto, scf, tools
import numpy as np

cell = gto.Cell()
# ... (build cell as above)

# Self-consistent calculation
kmf = scf.KRHF(cell)
kmf.kernel()

# Band structure along high-symmetry path
kpath = tools.pbc.get_kpath(cell, [10, 10, 10])
bands = []

for k in kpath:
    kmf_k = scf.KRHF(cell, k)
    kmf_k.kernel(dm0=kmf.make_rdm1())
    bands.append(kmf_k.mo_energy)

# Plot bands (requires matplotlib)

性能优化

密度拟合

from pyscf import gto, scf

mol = gto.M(atom='protein_fragment.xyz', basis='def2-svp')

# Use density fitting for faster calculation
mf = scf.RHF(mol).density_fit()
mf.kernel()

# Can also specify auxiliary basis
mf = scf.RHF(mol).density_fit(auxbasis='def2-svp-jkfit')
mf.kernel()

并行计算

from pyscf import gto, scf, lib

# Set number of threads
lib.num_threads(8)

mol = gto.M(atom='large_molecule.xyz', basis='6-31g')
mf = scf.RHF(mol)

# Enable parallel computation
mf.max_memory = 4000  # MB per thread
energy = mf.kernel()

内存管理

from pyscf import gto, scf

mol = gto.M(atom='large_system.xyz', basis='6-31g')
mf = scf.RHF(mol)

# Set maximum memory (MB)
mf.max_memory = 8000  # 8 GB

# Use less memory (slower)
mf.direct_scf = True  # Direct SCF, no integrals stored

energy = mf.kernel()

检查点与重启

from pyscf import gto, scf, lib

mol = gto.M(atom='molecule.xyz', basis='6-31g')
mf = scf.RHF(mol)

# Save checkpoint
mf.chkfile = 'calculation.chk'
energy = mf.kernel()

# Restart from checkpoint
mf2 = scf.RHF(mol)
mf2.chkfile = 'calculation.chk'
mf2.init_guess = 'chkfile'
energy2 = mf2.kernel()

与其他工具集成

导出为 Molden 格式

from pyscf import gto, scf
from pyscf.tools import molden

mol = gto.M(atom='H2O', basis='6-31g')
mf = scf.RHF(mol)
mf.kernel()

# Write Molden file for visualization
with open('molecule.molden', 'w') as f:
    molden.header(mol, f)
    molden.orbital_coeff(mol, f, mf.mo_coeff, ene=mf.mo_energy, occ=mf.mo_occ)

与 ASE 接口

from pyscf import gto, scf
from ase import Atoms
from ase.optimize import BFGS

# Create ASE calculator wrapper
class PySCFCalculator:
    def __init__(self, basis='6-31g', method='RHF'):
        self.basis = basis
        self.method = method
    
    def get_potential_energy(self, atoms):
        mol = gto.M(
            atom=[(atoms.get_chemical_symbols()[i], atoms.positions[i]) 
                  for i in range(len(atoms))],
            basis=self.basis,
            unit='angstrom'
        )
        
        if self.method == 'RHF':
            mf = scf.RHF(mol)
        elif self.method == 'B3LYP':
            mf = dft.RKS(mol)
            mf.xc = 'B3LYP'
        
        return mf.kernel()

# Use with ASE
atoms = Atoms('H2O', positions=[[0,0,0], [0,1,0], [0,0,1]])
atoms.calc = PySCFCalculator()

# Optimize with ASE
opt = BFGS(atoms)
opt.run(fmax=0.01)

常见问题与解决方案

SCF 收敛失败

from pyscf import gto, scf

mol = gto.M(atom='difficult.xyz', basis='6-31g')
mf = scf.RHF(mol)

# Problem: SCF not converging
energy = mf.kernel()
if not mf.converged:
    # Solution 1: Use level shift
    mf.level_shift = 0.5
    mf.kernel()
    
    # Solution 2: Increase DIIS space
    if not mf.converged:
        mf.diis_space = 15
        mf.kernel()
    
    # Solution 3: Try different initial guess
    if not mf.converged:
        mf.init_guess = 'atom'
        mf.kernel()
    
    # Solution 4: Use smaller basis first
    if not mf.converged:
        mol_small = gto.M(atom='difficult.xyz', basis='sto-3g')
        mf_small = scf.RHF(mol_small).run()
        mf.init_guess = mf_small.make_rdm1()
        mf.kernel()

内存问题

from pyscf import gto, scf

# Problem: Running out of memory
mol = gto.M(atom='huge.xyz', basis='def2-tzvp')
mf = scf.RHF(mol)

# Solution 1: Use density fitting
mf = mf.density_fit()

# Solution 2: Set memory limit
mf.max_memory = 4000  # MB

# Solution 3: Use direct SCF
mf.direct_scf = True

energy = mf.kernel()

自旋污染

from pyscf import gto, scf

mol = gto.M(atom='radical.xyz', basis='6-31g', spin=1)

# Problem: UHF has spin contamination
mf = scf.UHF(mol)
mf.kernel()
s2 = mf.spin_square()[0]
print(f"<S^2> = {s2:.4f}  (expected 0.75)")

# Solution: Use ROHF if appropriate
if abs(s2 - 0.75) > 0.1:
    mf = scf.ROHF(mol)
    mf.kernel()
    print(f"Using ROHF instead")

基组线性相关性

from pyscf import gto, scf

# Problem: Large basis sets can have linear dependencies
mol = gto.M(atom='molecule.xyz', basis='aug-cc-pvqz')
mf = scf.RHF(mol)

# Solution: Remove linear dependencies
from pyscf import lo
mf.kernel()

if not mf.converged:
    # Check condition number
    s = mf.get_ovlp()
    cond = np.linalg.cond(s)
    if cond > 1e10:
        print("Linear dependency detected")
        # Use canonical orthogonalization
        mf = scf.addons.remove_linear_dep_(mf)
        mf.kernel()

这份全面的 PySCF 指南涵盖了所有主要量子化学工作流中的 50 多个示例！

支持 Markdown

登录后评论

基本信息

作者: tondevrel
仓库: tondevrel/scientific-agent-skills
分类: 其他
安装次数: 58
GitHub Stars: 24

报告问题