核心模块 (core)

本模块提供分子动力学模拟的基础数据结构。

结构模块

分子动力学结构模块

该模块提供分子动力学模拟中的基础数据结构，包括原子和晶胞的表示。实现了周期性边界条件、最小镜像约定等关键算法。

理论基础

周期性边界条件 (Periodic Boundary Conditions, PBC)：

\[\mathbf{r}_{frac} = \mathbf{L}^{-T} \cdot \mathbf{r}_{cart}\]

其中： - \(\mathbf{r}_{cart}\) - 笛卡尔坐标 - \(\mathbf{r}_{frac}\) - 分数坐标 - \(\mathbf{L}\) - 晶格矢量矩阵

最小镜像约定 (Minimum Image Convention)：

\[\mathbf{r}_{min} = \mathbf{r} - \mathbf{L} \cdot \text{round}(\mathbf{L}^{-1} \cdot \mathbf{r})\]

Classes

Atom: 表示单个原子，包含位置、速度、力等属性
Cell: 表示晶胞结构，管理原子集合和晶格矢量

Functions

_apply_pbc_numba: JIT优化的周期性边界条件应用函数

备注

本模块是ThermoElasticSim的核心组件，为分子动力学模拟提供基础数据结构。所有长度单位为埃(Å)，时间单位为飞秒(fs)，能量单位为电子伏特(eV)。

示例

创建简单的铝晶胞：

>>> from thermoelasticsim.core.structure import Atom, Cell
>>> import numpy as np
>>> # 创建FCC铝的原胞
>>> a = 4.05  # 晶格常数
>>> lattice = np.array([[a, 0, 0], [0, a, 0], [0, 0, a]])
>>> atoms = [
...     Atom(0, "Al", 26.98, [0, 0, 0]),
...     Atom(1, "Al", 26.98, [a/2, a/2, 0]),
...     Atom(2, "Al", 26.98, [a/2, 0, a/2]),
...     Atom(3, "Al", 26.98, [0, a/2, a/2])
... ]
>>> cell = Cell(lattice, atoms)

Added in version 4.0.0.

class thermoelasticsim.core.structure.Atom(id, symbol, mass_amu, position, velocity=None)[源代码]

基类：object

原子对象，分子动力学模拟的基本单元

表示一个原子的完整状态，包括位置、速度、受力等物理属性。在MD模拟中，原子遵循牛顿运动方程：

\[m \frac{d^2\mathbf{r}}{dt^2} = \mathbf{F}\]

参数:

id (int) -- 原子的唯一标识符
symbol (str) -- 元素符号 (如 'Al', 'Cu', 'Si')
mass_amu (float) -- 原子质量，原子质量单位 (amu)
position (array_like) -- 初始位置，3D笛卡尔坐标 (Å)
velocity (array_like, optional) -- 初始速度，3D笛卡尔坐标 (Å/fs)，默认为零

id

原子唯一标识符

Type:: int

symbol

元素符号

Type:: str

mass_amu

原始质量 (amu)

Type:: float

mass

转换后的质量 (eV·fs²/Å²)

Type:: float

position

当前位置向量 (3,)

Type:: numpy.ndarray

velocity

当前速度向量 (3,)

Type:: numpy.ndarray

force

当前受力向量 (3,)

Type:: numpy.ndarray

备注

质量单位转换关系： 1 amu = 1.66054e-27 kg = 103.62 eV·fs²/Å²

示例

创建铝原子：

>>> atom = Atom(id=1, symbol="Al", mass_amu=26.98, position=[0, 0, 0])
>>> atom.velocity = np.array([1.0, 0.0, 0.0])  # 设置速度
>>> print(f"动能: {0.5 * atom.mass * np.dot(atom.velocity, atom.velocity):.4f} eV")

Methods

`accelerate_by`(velocity_change)	通过速度增量改变原子速度
`copy`()	创建 Atom 的深拷贝
`move_by`(displacement)	通过位置增量移动原子

__init__(id, symbol, mass_amu, position, velocity=None)[源代码]

参数:

id (int)
symbol (str)
mass_amu (float)
position (ndarray)
velocity (ndarray | None)

返回类型:

None

move_by(displacement)[源代码]

通过位置增量移动原子

Args:: displacement: 位置增量向量，形状为(3,)

抛出:: ValueError -- 如果位置增量不是3D向量:
参数:: displacement (ndarray)
返回类型:: None

示例

>>> atom = Atom(1, 'H', 1.0, [0, 0, 0])
>>> atom.move_by([0.1, 0.2, 0.3])
>>> print(atom.position)
[0.1 0.2 0.3]

accelerate_by(velocity_change)[源代码]

通过速度增量改变原子速度

Args:: velocity_change: 速度增量向量，形状为(3,)

抛出:: ValueError -- 如果velocity_change不是3D向量:
参数:: velocity_change (ndarray)
返回类型:: None

示例

>>> atom = Atom(1, 'H', 1.0, [0, 0, 0])
>>> atom.accelerate_by([0.1, 0.2, 0.3])
>>> print(atom.velocity)
[0.1 0.2 0.3]

copy()[源代码]

创建 Atom 的深拷贝

返回类型:: 新的Atom对象，包含当前原子的所有属性的副本

示例

>>> atom1 = Atom(1, 'H', 1.0, [0, 0, 0])
>>> atom2 = atom1.copy()
>>> atom2.id == atom1.id
True
>>> atom2.position is atom1.position
False

class thermoelasticsim.core.structure.Cell(lattice_vectors, atoms, pbc_enabled=True)[源代码]

基类：object

晶胞对象，管理原子集合和晶格结构

晶胞是分子动力学模拟的基本容器，定义了模拟盒子的几何形状和其中包含的原子。支持周期性边界条件和最小镜像约定。

晶格矢量定义为行向量：

\[\begin{split}\mathbf{L} = \begin{pmatrix} \mathbf{a}_1 \\ \mathbf{a}_2 \\ \mathbf{a}_3 \end{pmatrix}\end{split}\]

晶胞体积计算：

\[V = |\mathbf{a}_1 \cdot (\mathbf{a}_2 \times \mathbf{a}_3)|\]

参数:

lattice_vectors (array_like) -- 3×3晶格矢量矩阵，每行为一个晶格矢量 (Å)
atoms (list of Atom) -- 晶胞中的原子列表
pbc_enabled (bool, optional) -- 是否启用周期性边界条件，默认True

lattice_vectors

晶格矢量矩阵 (3, 3)

Type:: numpy.ndarray

atoms

原子对象列表

Type:: list of Atom

volume

晶胞体积 (Å³)

Type:: float

pbc_enabled

周期性边界条件标志

Type:: bool

lattice_locked

晶格锁定标志（用于内部弛豫）

Type:: bool

lattice_inv

晶格逆矩阵，用于坐标转换

Type:: numpy.ndarray

num_atoms

原子数量（属性）

Type:: int

apply_periodic_boundary(positions)[源代码]: 应用周期性边界条件

minimum_image(displacement)[源代码]: 计算最小镜像位移

calculate_temperature()[源代码]

计算瞬时温度

返回类型:: float

build_supercell(repetition)[源代码]

构建超胞

参数:: repetition (tuple)
返回类型:: Cell

备注

分数坐标与笛卡尔坐标转换： - 笛卡尔→分数: \(\mathbf{s} = \mathbf{L}^{-T} \cdot \mathbf{r}\) - 分数→笛卡尔: \(\mathbf{r} = \mathbf{L}^T \cdot \mathbf{s}\)

示例

创建FCC铝晶胞：

>>> a = 4.05  # 晶格常数
>>> lattice = a * np.eye(3)  # 立方晶格
>>> atoms = [
...     Atom(0, "Al", 26.98, [0, 0, 0]),
...     Atom(1, "Al", 26.98, [a/2, a/2, 0]),
...     Atom(2, "Al", 26.98, [a/2, 0, a/2]),
...     Atom(3, "Al", 26.98, [0, a/2, a/2])
... ]
>>> cell = Cell(lattice, atoms)
>>> print(f"晶胞体积: {cell.volume:.2f} Å³")

Attributes:

num_atoms: 返回原子数量

参数:

lattice_vectors (ndarray)
atoms (list[Atom])
pbc_enabled (bool)

Methods

`apply_deformation`(deformation_matrix)	对晶胞和原子坐标施加变形矩阵F
`apply_periodic_boundary`(positions)	应用周期性边界条件
`build_supercell`(repetition)	构建超胞，返回一个新的 Cell 对象。
`calculate_kinetic_energy`()	计算当前系统的总动能
`calculate_stress_tensor`(potential)	计算应力张量
`calculate_temperature`()	计算系统的瞬时温度
`calculate_volume`()	计算晶胞的体积
`copy`()	创建 Cell 的深拷贝
`get_box_lengths`()	返回模拟盒子在 x、y、z 方向的长度
`get_com_position`()	计算系统质心位置。
`get_com_velocity`()	计算系统质心速度。
`get_forces`()	获取所有原子的力信息
`get_fractional_coordinates`()	获取所有原子的分数坐标
`get_positions`()	获取所有原子的位置信息
`get_velocities`()	获取所有原子的速度信息
`get_volume`()	计算晶胞体积
`lock_lattice_vectors`()	锁定晶格向量，防止在优化过程中被修改
`minimum_image`(displacement)	计算最小镜像位移向量
`remove_com_motion`()	移除系统的质心运动。
`scale_lattice`(scale_factor)	按给定因子等比例缩放晶格
`set_fractional_coordinates`(fractional_coords)	根据分数坐标设置所有原子的笛卡尔坐标
`set_lattice_vectors`(lattice_vectors)	设置新的晶格矢量并更新相关属性。
`set_positions`(positions)	设置所有原子的笛卡尔坐标
`unlock_lattice_vectors`()	解锁晶格向量，允许在需要时修改

__init__(lattice_vectors, atoms, pbc_enabled=True)[源代码]

参数:

lattice_vectors (ndarray)
atoms (list[Atom])
pbc_enabled (bool)

返回类型:

None

calculate_volume()[源代码]

计算晶胞的体积

返回:: 晶胞体积，单位为Å^3
返回类型:: float

get_box_lengths()[源代码]

返回模拟盒子在 x、y、z 方向的长度

返回:: 包含三个方向长度的数组，单位为Å
返回类型:: numpy.ndarray

calculate_stress_tensor(potential)[源代码]

计算应力张量

参数:: potential (Potential) -- 用于计算应力的势能对象
返回:: 3x3应力张量矩阵，单位为eV/Å^3
返回类型:: numpy.ndarray

lock_lattice_vectors()[源代码]

锁定晶格向量，防止在优化过程中被修改

备注

锁定后，晶格向量将不能被修改，直到调用unlock_lattice_vectors()

返回类型:: None

unlock_lattice_vectors()[源代码]

解锁晶格向量，允许在需要时修改

备注

解锁后，晶格向量可以被修改，直到再次调用lock_lattice_vectors()

返回类型:: None

apply_deformation(deformation_matrix)[源代码]

对晶胞和原子坐标施加变形矩阵F

Args:: deformation_matrix: 3x3变形矩阵F

备注

当晶格未锁定时： 1. 更新晶格矢量： L_new = F * L_old 2. 更新原子坐标： r_new = r_old * F^T

当晶格已锁定时：仅对原子坐标施加F变形： r_new = r_old * F^T

示例

>>> import numpy as np
>>> # 小幅变形矩阵
>>> F = np.array([[1.01, 0, 0], [0, 1, 0], [0, 0, 1]])
>>> cell.apply_deformation(F)

参数:: deformation_matrix (ndarray)
返回类型:: None

apply_periodic_boundary(positions)[源代码]

应用周期性边界条件

将原子位置映射回主晶胞内，使用标准的最小镜像约定。

算法：

\[ \begin{align}\begin{aligned}\mathbf{s} = \mathbf{L}^{-T} \cdot \mathbf{r}\\\mathbf{s}' = \mathbf{s} - \lfloor \mathbf{s} + 0.5 \rfloor\\\mathbf{r}' = \mathbf{L}^T \cdot \mathbf{s}'\end{aligned}\end{align} \]

参数:: positions (numpy.ndarray) -- 原子位置坐标，形状为 (3,) 或 (N, 3)
返回:: 应用PBC后的位置，保持输入形状
返回类型:: numpy.ndarray
抛出:: ValueError -- 如果位置数组形状不正确或包含非有限值

备注

该实现使用分数坐标进行周期性映射，确保数值稳定性。对于三斜晶系也能正确处理。

示例

>>> # 单个原子位置
>>> pos = np.array([5.0, 6.0, 7.0])
>>> new_pos = cell.apply_periodic_boundary(pos)
>>> # 批量处理
>>> positions = np.random.randn(100, 3) * 10
>>> new_positions = cell.apply_periodic_boundary(positions)

copy()[源代码]: 创建 Cell 的深拷贝

calculate_temperature()[源代码]

计算系统的瞬时温度

使用均分定理计算温度，扣除质心运动贡献：

\[T = \frac{2 E_{kin}}{k_B \cdot N_{dof}}\]

其中动能（扣除质心）：

\[E_{kin} = \sum_i \frac{1}{2} m_i |\mathbf{v}_i - \mathbf{v}_{cm}|^2\]

自由度数： - 多原子系统: \(N_{dof} = 3N - 3\) (扣除质心平动) - 单原子系统: \(N_{dof} = 3\) (不扣除)

返回:: 系统温度 (K)
返回类型:: float

备注

瞬时温度会有涨落，通常需要时间平均来获得稳定值。对于NVE系综，温度是守恒量的函数。

示例

>>> temp = cell.calculate_temperature()
>>> print(f"瞬时温度: {temp:.1f} K")
>>> # 计算时间平均温度
>>> temps = [cell.calculate_temperature() for _ in range(1000)]
>>> avg_temp = np.mean(temps)

参见

calculate_kinetic_energy: 计算总动能（含质心运动）

calculate_kinetic_energy()[源代码]

计算当前系统的总动能

返回类型:: 系统总动能，单位为eV

备注

计算所有原子的动能总和，包含质心运动

示例

>>> cell = Cell(lattice_vectors, atoms)
>>> kinetic = cell.calculate_kinetic_energy()
>>> print(f"系统动能: {kinetic:.6f} eV")

get_positions()[源代码]

获取所有原子的位置信息

返回类型:: 原子位置数组，形状为(num_atoms, 3)

示例

>>> positions = cell.get_positions()
>>> print(f"原子数量: {positions.shape[0]}")

get_velocities()[源代码]

获取所有原子的速度信息

返回:: 原子速度数组，形状为 (num_atoms, 3)
返回类型:: numpy.ndarray

get_forces()[源代码]

获取所有原子的力信息

返回:: 原子力数组，形状为 (num_atoms, 3)
返回类型:: numpy.ndarray

minimum_image(displacement)[源代码]

计算最小镜像位移向量

根据最小镜像约定，找到最近的周期镜像之间的位移。这对于正确计算周期性系统中的距离和力至关重要。

数学原理：

\[\mathbf{d}_{min} = \mathbf{d} - \mathbf{L} \cdot \text{round}(\mathbf{L}^{-1} \cdot \mathbf{d})\]

参数:: displacement (numpy.ndarray) -- 原始位移向量 (3,)
返回:: 最小镜像位移向量 (3,)
返回类型:: numpy.ndarray
抛出:: ValueError -- 如果位移向量不是3D或包含非有限值

备注

该方法确保返回的位移向量长度最小，对应于最近的周期镜像。在计算原子间相互作用时必须使用此方法。

示例

>>> # 计算两原子间的最小距离
>>> r12 = atom2.position - atom1.position
>>> r12_min = cell.minimum_image(r12)
>>> distance = np.linalg.norm(r12_min)

参见

apply_periodic_boundary: 应用周期性边界条件

build_supercell(repetition)[源代码]

构建超胞，返回一个新的 Cell 对象。

采用更简单、更健壮的算法，直接在笛卡尔坐标系下操作。

参数:: repetition (tuple of int) -- 在 x, y, z 方向上的重复次数，例如 (2, 2, 2)
返回:: 新的超胞对象
返回类型:: Cell

get_com_velocity()[源代码]

计算系统质心速度。

返回:: 质心速度向量
返回类型:: numpy.ndarray

remove_com_motion()[源代码]

移除系统的质心运动。

该方法计算并移除系统的净平动，保证总动量为零。在应用恒温器和恒压器时应该调用此方法。

get_com_position()[源代码]

计算系统质心位置。

返回:: 质心位置向量
返回类型:: numpy.ndarray

property num_atoms: 返回原子数量

set_lattice_vectors(lattice_vectors)[源代码]

设置新的晶格矢量并更新相关属性。

参数:: lattice_vectors (np.ndarray) -- 新的 3x3 晶格矢量矩阵。

get_fractional_coordinates()[源代码]

获取所有原子的分数坐标

返回:: 原子分数坐标数组，形状为(num_atoms, 3)
返回类型:: numpy.ndarray

备注

分数坐标：r_frac = (L^T)^-1 * r_cart 其中L是晶格矢量矩阵，r_cart是笛卡尔坐标

set_fractional_coordinates(fractional_coords)[源代码]

根据分数坐标设置所有原子的笛卡尔坐标

参数:: fractional_coords (numpy.ndarray) -- 分数坐标数组，形状为(num_atoms, 3)
返回类型:: None

备注

笛卡尔坐标：r_cart = L^T * r_frac 其中L是晶格矢量矩阵，r_frac是分数坐标

set_positions(positions)[源代码]

设置所有原子的笛卡尔坐标

参数:: positions (numpy.ndarray) -- 笛卡尔坐标数组，形状为(num_atoms, 3)
返回类型:: None

get_volume()[源代码]

计算晶胞体积

返回:: 晶胞体积 (Å³)
返回类型:: float

备注

体积通过晶格矢量的标量三重积计算： V = |a · (b × c)|

scale_lattice(scale_factor)[源代码]

按给定因子等比例缩放晶格

参数:: scale_factor (float) -- 缩放因子，>1放大，<1缩小
返回类型:: None

备注

这个方法只缩放晶格矢量，不改变原子坐标。通常与原子坐标的相应缩放一起使用。

示例

>>> cell.scale_lattice(1.1)  # 放大10%
>>> # 同时需要缩放原子坐标以保持相对位置
>>> for atom in cell.atoms:
...     atom.position *= 1.1

配置模块

配置模块

包含 ConfigManager 类，用于加载和管理 YAML 格式的配置文件

Classes:: ConfigManager: 加载和访问配置文件的配置管理器

class thermoelasticsim.core.config.ConfigManager(config_file)[源代码]

基类：object

加载和管理 YAML 配置文件的配置管理器

参数:: config_file (str)

config

存储加载的配置数据

Type:: dict

Methods

`get`(key[, default])	获取配置参数的值
`load_config`(config_file)	从 YAML 文件加载配置数据

__init__(config_file)[源代码]

初始化配置管理器并加载配置文件

参数:

config_file (str) -- 配置文件的路径

抛出:

FileNotFoundError -- 如果配置文件不存在
yaml.YAMLError -- 如果配置文件不是有效的YAML格式

返回类型:

None

static load_config(config_file)[源代码]

从 YAML 文件加载配置数据

参数:

config_file (str) -- 配置文件的路径

返回:

包含配置数据的字典

返回类型:

dict

抛出:

FileNotFoundError -- 如果配置文件不存在
yaml.YAMLError -- 如果配置文件不是有效的YAML格式

备注

使用yaml.safe_load()来避免潜在的安全问题

get(key, default=None)[源代码]

获取配置参数的值

参数:

key (str) -- 要获取的配置键名
default (Any, optional) -- 如果键不存在时返回的默认值 (默认: None)

返回:

配置值或默认值

返回类型:

Any