工具模块 (utils)

本模块提供各种辅助工具和实用函数。

通用工具

工具模块

包含 TensorConverter 类用于张量与 Voigt 表示之间的转换，DataCollector 类用于收集模拟过程中的数据， 及一些常用的单位转换常量

thermoelasticsim.utils.utils.AMU_TO_GRAM: float = 1.6605390402e-24

原子质量单位 amu → 克 g 的换算系数。

thermoelasticsim.utils.utils.ANGSTROM_TO_CM: float = 1e-08

长度单位换算系数：Å → cm。

class thermoelasticsim.utils.utils.TensorConverter

基类：object

张量转换工具类，支持应力和应变张量的 Voigt 与 3x3 矩阵表示之间的相互转换。

static to_voigt(tensor, tensor_type, tol=1e-08)

将一个对称的 3x3 张量转换为 Voigt 表示的 6 元素向量。

参数:

• tensor (np.ndarray) -- 形状为 (3, 3) 的张量。

• tensor_type (str) -- 张量类型，必须是 'stress' 或 'strain'。

• tol (float, optional) -- 用于检查张量对称性的容差。如果非对称分量差异大于此值，将记录警告。

返回:

形状为 (6,) 的 Voigt 表示向量。

返回类型:

np.ndarray

static from_voigt(voigt, tensor_type)

将 Voigt 表示的 6 元素向量转换为 3x3 的对称张量。

参数:

• voigt (np.ndarray) -- 形状为 (6,) 的 Voigt 向量。

• tensor_type (str) -- 张量类型，必须是 'stress' 或 'strain'。

返回:

形状为 (3, 3) 的对称张量。

返回类型:

np.ndarray

class thermoelasticsim.utils.utils.DataCollector(capacity=None, save_interval=1000, use_threading=False, output_file='trajectory.h5')

基类：object

数据收集器，用于记录和分析模拟轨迹

提供了系统状态的收集、统计和保存功能，支持异步数据处理和自动保存。

参数:

• capacity (int, optional) -- 存储容量限制，None 表示无限制

• save_interval (int, optional) -- 数据自动保存间隔，默认为 1000 步

• use_threading (bool, optional) -- 是否使用线程进行数据处理，默认为 False

• output_file (str, optional) -- 数据保存的文件名，默认为 'trajectory.h5'

__init__(capacity=None, save_interval=1000, use_threading=False, output_file='trajectory.h5')

参数:

• capacity (int | None)

• save_interval (int)

• use_threading (bool)

• output_file (str)

__del__()

确保线程池正确关闭

collect(cell, potential=None)

收集完整的系统状态数据

参数:

• cell (Cell) -- 当前的晶胞对象

• potential (Potential, optional) -- 势能对象，用于计算系统能量

save_trajectory()

将轨迹数据保存到 HDF5 文件

class thermoelasticsim.utils.utils.NeighborList(cutoff, skin=0.3)

基类：object

邻居列表类，用于生成和维护原子的邻居列表

提供高效的邻居搜索和更新功能，支持周期性边界条件。

参数:

• cutoff (float) -- 截断半径，单位为 Å

• skin (float, optional) -- 皮肤厚度，默认为 0.3 Å

__init__(cutoff, skin=0.3)

参数:

• cutoff (float)

• skin (float)

get_neighbor_stats()

返回邻居列表的统计信息

返回:

包含最小、最大、平均邻居数等统计信息的字典

返回类型:

dict

build(cell)

构建邻居列表

need_refresh()

判断是否需要更新邻居列表。

返回:

如果需要更新，返回 True；否则返回 False。

返回类型:

bool

update()

更新邻居列表，如果需要的话。

get_neighbors(atom_index)

获取指定原子的邻居列表。

参数:

atom_index (int) -- 原子的索引。

返回:

邻居原子的索引列表。

返回类型:

list of int

debug_neighbor_distribution()

分析和打印邻居分布情况

thermoelasticsim.utils.utils.get_atomic_mass(symbol)

根据原子符号获取原子质量（amu）

参数:

symbol (str) -- 原子符号，例如 'Al', 'Si'

返回:

原子质量（amu）

返回类型:

float

抛出:

KeyError -- 如果原子符号不存在

常量

thermoelasticsim.utils.utils.EV_TO_GPA: float = 160.2176634

应力单位换算系数：eV/ų → GPa。

thermoelasticsim.utils.utils.KB_IN_EV: float = 8.617332385e-05

玻尔兹曼常数 kB（单位 eV/K）。

thermoelasticsim.utils.utils.AMU_TO_EVFSA2: float = 104.3968445

原子质量单位 (amu) 转换为 eV·fs²/Ų 的系数。

优化器

分子动力学模拟优化器模块

提供多种优化算法用于原子结构优化，包括： - GradientDescentOptimizer: 带动量项的梯度下降 - BFGSOptimizer: 基于scipy.optimize.minimize的BFGS - LBFGSOptimizer: 改进的L-BFGS(有限内存BFGS)

class thermoelasticsim.utils.optimizers.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

备注

质量与常数采用全局统一定义（参见 thermoelasticsim.utils.utils 模块）：

AMU_TO_EVFSA2

质量单位转换常量（amu → eV·fs²/Ų）。

KB_IN_EV

玻尔兹曼常数（eV/K）。

示例

创建铝原子：

>>> 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")

__init__(id, symbol, mass_amu, position, velocity=None)

参数:

• id (int)

• symbol (str)

• mass_amu (float)

• position (ndarray)

• velocity (ndarray | None)

返回类型:

None

accelerate_by(velocity_change)

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

参数:

velocity_change (numpy.ndarray | array_like) -- 速度增量向量，形状为 (3,)

抛出:

ValueError -- 如果速度增量不是3D向量

返回类型:

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 对象，包含当前原子的所有属性的副本

返回类型:

Atom

示例

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

move_by(displacement)

通过位置增量移动原子

参数:

displacement (numpy.ndarray | array_like) -- 位置增量向量，形状为 (3,)

抛出:

ValueError -- 如果位置增量不是3D向量

返回类型:

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]

class thermoelasticsim.utils.optimizers.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)

返回类型:

thermoelasticsim.core.structure.Cell

备注

分数/笛卡尔坐标转换（列向量记号）：

\[\mathbf{r} = \mathbf{L}\,\mathbf{s},\qquad \mathbf{s} = \mathbf{L}^{-1}\,\mathbf{r}\]

实现采用“行向量右乘”的等价式：

\[\mathbf{s}^{\top} = \mathbf{r}^{\top}\,\mathbf{L}^{-1},\qquad \mathbf{r}^{\top} = \mathbf{s}^{\top}\,\mathbf{L}^{\top}\]

示例

创建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} ų")

__init__(lattice_vectors, atoms, pbc_enabled=True)

参数:

• lattice_vectors (ndarray)

• atoms (list[Atom])

• pbc_enabled (bool)

返回类型:

None

apply_deformation(deformation_matrix)

对晶格与原子坐标施加形变矩阵 \(\mathbf{F}\)

参数:

deformation_matrix (numpy.ndarray) -- 形变矩阵 \(\mathbf{F}\)，形状为 (3, 3)

返回类型:

None

备注

采用列向量约定的连续介质力学记号：

• 晶格未锁定（允许变胞）时：

  \[\mathbf{L}_{\text{new}} = \mathbf{F}\, \mathbf{L}_{\text{old}},\qquad \mathbf{r}_{\text{new}} = \mathbf{F}\, \mathbf{r}_{\text{old}}\]

• 晶格已锁定（固定胞形）时：

  \[\mathbf{L}_{\text{new}} = \mathbf{L}_{\text{old}},\qquad \mathbf{r}_{\text{new}} = \mathbf{F}\, \mathbf{r}_{\text{old}}\]

实现细节：本实现使用“行向量右乘”形式，等价写法为

\[\mathbf{r}_{\text{new}}^{\top} = \mathbf{r}_{\text{old}}^{\top} \mathbf{F}^{\top}.\]

示例

>>> import numpy as np
>>> F = np.array([[1.01, 0, 0], [0, 1, 0], [0, 0, 1]])
>>> cell.apply_deformation(F)

apply_periodic_boundary(positions)

应用周期性边界条件

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

算法（列向量记号）：

\[\mathbf{s} = \mathbf{L}^{-1}\,\mathbf{r},\qquad \mathbf{s}' = \mathbf{s} - \lfloor \mathbf{s} + 0.5 \rfloor,\qquad \mathbf{r}' = \mathbf{L}\,\mathbf{s}'\]

参数:

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)

build_supercell(repetition)

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

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

参数:

repetition (tuple of int) -- 在 x, y, z 方向上的重复次数，例如 (2, 2, 2)

返回:

新的超胞对象

返回类型:

Cell

calculate_kinetic_energy()

计算当前系统的总动能

返回:

系统总动能，单位为 \(\mathrm{eV}\)

返回类型:

float

备注

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

示例

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

calculate_stress_tensor(potential)

计算应力张量

参数:

potential (Potential) -- 用于计算应力的势能对象

返回:

3×3 应力张量矩阵 (eV/ų)

返回类型:

numpy.ndarray

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\]

自由度数定义：

\[\begin{split}N_{\mathrm{dof}} = \begin{cases} 3N - 3, & N > 1, \\ 3, & N = 1. \end{cases}\end{split}\]

返回:

系统温度，单位为 \(\mathrm{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_volume()

计算晶胞的体积

返回:

晶胞体积 (ų)

返回类型:

float

copy()

创建 Cell 的深拷贝

get_box_lengths()

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

返回:

包含三个方向长度的数组 (Å)

返回类型:

numpy.ndarray

get_com_position()

计算系统质心位置。

返回:

质心位置向量

返回类型:

numpy.ndarray

get_com_velocity()

计算系统质心速度。

返回:

质心速度向量

返回类型:

numpy.ndarray

get_forces()

获取所有原子的力信息

返回:

原子力数组，形状为 (num_atoms, 3)

返回类型:

numpy.ndarray

get_fractional_coordinates()

获取所有原子的分数坐标

返回:

原子分数坐标数组，形状为(num_atoms, 3)

返回类型:

numpy.ndarray

备注

坐标关系（列向量记号）：

\[\mathbf{s} = \mathbf{L}^{-1}\,\mathbf{r}\]

代码实现采用行向量右乘：\(\mathbf{s}^{\top} = \mathbf{r}^{\top}\,\mathbf{L}^{-1}\)。

get_positions()

获取所有原子的位置信息

返回:

原子位置数组，形状为 (num_atoms, 3)

返回类型:

numpy.ndarray

示例

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

get_velocities()

获取所有原子的速度信息

返回:

原子速度数组，形状为 (num_atoms, 3)

返回类型:

numpy.ndarray

get_volume()

计算晶胞体积

返回:

晶胞体积 (ų)

返回类型:

float

备注

体积通过晶格矢量的标量三重积计算：

\[V = \left|\, \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) \,\right|\]

lock_lattice_vectors()

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

备注

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

返回类型:

None

minimum_image(displacement)

计算最小镜像位移向量

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

数学原理（列向量记号）：

\[\mathbf{d}_{\min} = \mathbf{d} - \mathbf{L}\, \operatorname{round}(\mathbf{L}^{-1} \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

应用周期性边界条件

property num_atoms

返回原子数量

remove_com_motion()

移除系统的质心运动。

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

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

set_fractional_coordinates(fractional_coords)

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

参数:

fractional_coords (numpy.ndarray) -- 分数坐标数组，形状为(num_atoms, 3)

返回类型:

None

备注

坐标关系（列向量记号）：

\[\mathbf{r} = \mathbf{L}\,\mathbf{s}\]

代码实现采用行向量右乘：\(\mathbf{r}^{\top} = \mathbf{s}^{\top}\,\mathbf{L}^{\top}\)。

set_lattice_vectors(lattice_vectors)

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

参数:

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

set_positions(positions)

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

参数:

positions (numpy.ndarray) -- 笛卡尔坐标数组，形状为(num_atoms, 3)

返回类型:

None

set_velocities(velocities)

设置所有原子的速度

参数:

velocities (numpy.ndarray) -- 速度数组，形状为 (num_atoms, 3)

抛出:

ValueError -- 当输入数组形状与原子数不匹配时抛出

返回类型:

None

unlock_lattice_vectors()

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

备注

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

返回类型:

None

class thermoelasticsim.utils.optimizers.Potential(parameters, cutoff)

基类：ABC

势能计算的抽象基类。

提供原子间势能模型的统一接口，约定最小方法集以支持分子动力学 (MD) 中的力与能量评估。

参数:

• parameters (dict) -- 势能参数字典（按具体模型定义）。

• cutoff (float) -- 势能截断距离，单位 Å。

备注

• 所有具体势能模型至少需实现 calculate_forces 与 calculate_energy。

• 单位约定：能量 eV，长度 Å，力 eV/Å。

__init__(parameters, cutoff)

abstractmethod calculate_energy(cell, neighbor_list)

计算系统的总势能。

参数:

• cell (Cell) -- 晶胞与原子集合。

• neighbor_list (NeighborList) -- 邻居列表（按模型需要使用）。

返回:

系统总势能，单位 eV。

返回类型:

float

备注

抽象方法，必须由子类实现。

abstractmethod calculate_forces(cell, neighbor_list)

计算系统中所有原子的作用力。

参数:

• cell (Cell) -- 晶胞与原子集合。

• neighbor_list (NeighborList) -- 邻居列表（按模型需要使用）。

返回类型:

None

备注

• 需就地写入 force （单位 eV/Å）。

• 抽象方法，必须由子类实现。

class thermoelasticsim.utils.optimizers.Optimizer

基类：object

优化器抽象基类，定义优化方法的接口

子类需要实现optimize()方法来执行具体的优化算法。

optimize(cell, potential)

执行优化算法，需子类实现

参数:

• cell (Cell)

• potential (Potential)

返回类型:

tuple[bool, list]

optimize(cell, potential)

执行优化算法

参数:

• cell (Cell) -- 晶胞对象

• potential (Potential) -- 势能函数

返回:

返回(是否收敛, 轨迹数据)

返回类型:

tuple[bool, list]

备注

轨迹数据包含每步的: - 原子位置 - 晶胞体积 - 晶格矢量

class thermoelasticsim.utils.optimizers.GradientDescentOptimizer(maxiter=10000, tol=0.001, step_size=0.001, energy_tol=0.0001, beta=0.9)

基类：Optimizer

带动量项的梯度下降优化器

参数:

• maxiter (int, optional) -- 最大迭代步数，默认10000

• tol (float, optional) -- 力收敛阈值，默认1e-3 eV/Å

• step_size (float, optional) -- 步长，默认1e-3 Å

• energy_tol (float, optional) -- 能量变化阈值，默认1e-4 eV

• beta (float, optional) -- 动量系数[0,1)，默认0.9

converged

优化是否收敛的标志

Type:

bool

trajectory

记录优化轨迹的列表

Type:

list

__init__(maxiter=10000, tol=0.001, step_size=0.001, energy_tol=0.0001, beta=0.9)

初始化带动量项的梯度下降优化器

optimize(cell, potential)

执行带动量项的梯度下降优化

参数:

• cell (Cell) -- 晶胞对象

• potential (Potential) -- 势能函数

返回:

返回(是否收敛, 轨迹数据字典列表)

返回类型:

tuple[bool, list[dict]]

备注

收敛条件: - 最大原子力 < tol - 能量变化 < energy_tol

class thermoelasticsim.utils.optimizers.BFGSOptimizer(tol=1e-06, maxiter=10000)

基类：Optimizer

BFGS准牛顿优化器

参数:

• tol (float, optional) -- 收敛阈值，默认1e-6

• maxiter (int, optional) -- 最大迭代步数，默认10000

converged

优化是否收敛的标志

Type:

bool

trajectory

记录优化轨迹的列表

Type:

list

备注

使用scipy.optimize.minimize的BFGS实现 适用于中等规模系统优化

__init__(tol=1e-06, maxiter=10000)

初始化 BFGS 优化器

optimize(cell, potential)

执行 BFGS 优化

参数:

• cell (Cell) -- 晶胞对象

• potential (Potential) -- 势能函数

返回:

返回(是否收敛, 轨迹数据字典列表)

返回类型:

tuple[bool, list[dict]]

class thermoelasticsim.utils.optimizers.CGOptimizer(tol=1e-06, maxiter=10000)

基类：Optimizer

共轭梯度优化器（分数坐标），对剪切后的浅谷更稳健。

参数:

• tol (float, optional) -- 梯度收敛阈值（映射到 scipy 的 gtol），默认 1e-6

• maxiter (int, optional) -- 最大迭代步数，默认 10000

__init__(tol=1e-06, maxiter=10000)

参数:

• tol (float)

• maxiter (int)

optimize(cell, potential)

执行 CG 优化（变量为分数坐标）。

参数:

• cell (Cell)

• potential (Potential)

返回类型:

tuple[bool, list[dict]]

class thermoelasticsim.utils.optimizers.LBFGSOptimizer(ftol=1e-08, gtol=1e-06, maxiter=1000, supercell_dims=None, **kwargs)

基类：Optimizer

L-BFGS优化器(有限内存BFGS)

该优化器能够处理两种模式： 1. 仅优化原子位置 (默认)。 2. 同时优化原子位置和晶胞形状 (完全弛豫)。

参数:

• ftol (float, optional) -- 函数收敛阈值，默认1e-8

• gtol (float, optional) -- 梯度收敛阈值，默认1e-6

• maxiter (int, optional) -- 最大迭代步数，默认1000

• **kwargs (dict, optional) -- 其他要传递给 scipy.optimize.minimize 的选项, 例如 {'disp': True}

__init__(ftol=1e-08, gtol=1e-06, maxiter=1000, supercell_dims=None, **kwargs)

optimize(cell, potential, relax_cell=False)

执行 L-BFGS 优化。

参数:

• cell (Cell)

• potential (Potential)

• relax_cell (bool)

返回类型:

tuple[bool, list[dict]]

轨迹记录

轨迹存储和管理系统

基于H5MD (HDF5 for Molecular Dynamics) 标准的轨迹存储系统。 提供高效的数据压缩、并行I/O支持和丰富的元数据管理功能。

主要特性： - HDF5格式的高效存储 - 支持大规模轨迹数据 - 灵活的压缩选项 - 增量式写入支持 - 丰富的元数据记录

Author: Gilbert Young Created: 2025-08-15

class thermoelasticsim.utils.trajectory.datetime(year, month, day[, hour[, minute[, second[, microsecond[, tzinfo]]]]])

基类：date

The year, month and day arguments are required. tzinfo may be None, or an instance of a tzinfo subclass. The remaining arguments may be ints.

astimezone()

tz -> convert to local time in new timezone tz

classmethod combine()

date, time -> datetime with same date and time fields

ctime()

Return ctime() style string.

date()

Return date object with same year, month and day.

dst()

Return self.tzinfo.dst(self).

fold

classmethod fromisoformat()

string -> datetime from a string in most ISO 8601 formats

classmethod fromtimestamp()

timestamp[, tz] -> tz's local time from POSIX timestamp.

hour

isoformat()

[sep] -> string in ISO 8601 format, YYYY-MM-DDT[HH[:MM[:SS[.mmm[uuu]]]]][+HH:MM]. sep is used to separate the year from the time, and defaults to 'T'. The optional argument timespec specifies the number of additional terms of the time to include. Valid options are 'auto', 'hours', 'minutes', 'seconds', 'milliseconds' and 'microseconds'.

max = datetime.datetime(9999, 12, 31, 23, 59, 59, 999999)

microsecond

min = datetime.datetime(1, 1, 1, 0, 0)

minute

classmethod now(tz=None)

Returns new datetime object representing current time local to tz.

tz

Timezone object.

If no tz is specified, uses local timezone.

replace()

Return datetime with new specified fields.

resolution = datetime.timedelta(microseconds=1)

second

classmethod strptime()

string, format -> new datetime parsed from a string (like time.strptime()).

time()

Return time object with same time but with tzinfo=None.

timestamp()

Return POSIX timestamp as float.

timetuple()

Return time tuple, compatible with time.localtime().

timetz()

Return time object with same time and tzinfo.

tzinfo

tzname()

Return self.tzinfo.tzname(self).

classmethod utcfromtimestamp()

Construct a naive UTC datetime from a POSIX timestamp.

classmethod utcnow()

Return a new datetime representing UTC day and time.

utcoffset()

Return self.tzinfo.utcoffset(self).

utctimetuple()

Return UTC time tuple, compatible with time.localtime().

class thermoelasticsim.utils.trajectory.Path(*args, **kwargs)

基类：PurePath

PurePath subclass that can make system calls.

Path represents a filesystem path but unlike PurePath, also offers methods to do system calls on path objects. Depending on your system, instantiating a Path will return either a PosixPath or a WindowsPath object. You can also instantiate a PosixPath or WindowsPath directly, but cannot instantiate a WindowsPath on a POSIX system or vice versa.

absolute()

Return an absolute version of this path by prepending the current working directory. No normalization or symlink resolution is performed.

Use resolve() to get the canonical path to a file.

chmod(mode, *, follow_symlinks=True)

Change the permissions of the path, like os.chmod().

classmethod cwd()

Return a new path pointing to the current working directory (as returned by os.getcwd()).

exists()

Whether this path exists.

expanduser()

Return a new path with expanded ~ and ~user constructs (as returned by os.path.expanduser)

glob(pattern)

Iterate over this subtree and yield all existing files (of any kind, including directories) matching the given relative pattern.

group()

Return the group name of the file gid.

hardlink_to(target)

Make this path a hard link pointing to the same file as target.

Note the order of arguments (self, target) is the reverse of os.link's.

classmethod home()

Return a new path pointing to the user's home directory (as returned by os.path.expanduser('~')).

is_block_device()

Whether this path is a block device.

is_char_device()

Whether this path is a character device.

is_dir()

Whether this path is a directory.

is_fifo()

Whether this path is a FIFO.

is_file()

Whether this path is a regular file (also True for symlinks pointing to regular files).

is_mount()

Check if this path is a POSIX mount point

is_socket()

Whether this path is a socket.

is_symlink()

Whether this path is a symbolic link.

iterdir()

Iterate over the files in this directory. Does not yield any result for the special paths '.' and '..'.

lchmod(mode)

Like chmod(), except if the path points to a symlink, the symlink's permissions are changed, rather than its target's.

link_to(target)

Make the target path a hard link pointing to this path.

Note this function does not make this path a hard link to target, despite the implication of the function and argument names. The order of arguments (target, link) is the reverse of Path.symlink_to, but matches that of os.link.

Deprecated since Python 3.10 and scheduled for removal in Python 3.12. Use hardlink_to() instead.

lstat()

Like stat(), except if the path points to a symlink, the symlink's status information is returned, rather than its target's.

mkdir(mode=511, parents=False, exist_ok=False)

Create a new directory at this given path.

open(mode='r', buffering=-1, encoding=None, errors=None, newline=None)

Open the file pointed by this path and return a file object, as the built-in open() function does.

owner()

Return the login name of the file owner.

read_bytes()

Open the file in bytes mode, read it, and close the file.

read_text(encoding=None, errors=None)

Open the file in text mode, read it, and close the file.

readlink()

Return the path to which the symbolic link points.

rename(target)

Rename this path to the target path.

The target path may be absolute or relative. Relative paths are interpreted relative to the current working directory, not the directory of the Path object.

Returns the new Path instance pointing to the target path.

replace(target)

Rename this path to the target path, overwriting if that path exists.

The target path may be absolute or relative. Relative paths are interpreted relative to the current working directory, not the directory of the Path object.

Returns the new Path instance pointing to the target path.

resolve(strict=False)

Make the path absolute, resolving all symlinks on the way and also normalizing it.

rglob(pattern)

Recursively yield all existing files (of any kind, including directories) matching the given relative pattern, anywhere in this subtree.

rmdir()

Remove this directory. The directory must be empty.

samefile(other_path)

Return whether other_path is the same or not as this file (as returned by os.path.samefile()).

stat(*, follow_symlinks=True)

Return the result of the stat() system call on this path, like os.stat() does.

symlink_to(target, target_is_directory=False)

Make this path a symlink pointing to the target path. Note the order of arguments (link, target) is the reverse of os.symlink.

touch(mode=438, exist_ok=True)

Create this file with the given access mode, if it doesn't exist.

unlink(missing_ok=False)

Remove this file or link. If the path is a directory, use rmdir() instead.

write_bytes(data)

Open the file in bytes mode, write to it, and close the file.

write_text(data, encoding=None, errors=None, newline=None)

Open the file in text mode, write to it, and close the file.

class thermoelasticsim.utils.trajectory.Any(*args, **kwargs)

基类：object

Special type indicating an unconstrained type.

• Any is compatible with every type.

• Any assumed to have all methods.

• All values assumed to be instances of Any.

Note that all the above statements are true from the point of view of static type checkers. At runtime, Any should not be used with instance checks.

class thermoelasticsim.utils.trajectory.TrajectoryWriter(filename, mode='w', compression='gzip', compression_opts=4, chunk_size=100)

基类：object

H5MD格式轨迹文件写入器

遵循H5MD v1.1规范，提供高效的分子动力学轨迹存储。

参数:

• filename (str) -- 输出文件名（.h5或.hdf5扩展名）

• mode (str) -- 文件打开模式：'w'(覆盖), 'a'(追加), 'x'(创建新文件)

• compression (str, optional) -- 压缩算法：'gzip', 'lzf', 'szip'等

• compression_opts (int, optional) -- 压缩级别（gzip: 1-9）

• chunk_size (int, optional) -- 块大小，影响I/O性能

示例

>>> writer = TrajectoryWriter('trajectory.h5', compression='gzip')
>>> writer.initialize(n_atoms=100, n_frames_estimate=1000)
>>> writer.write_frame(positions, box, time=0.0, step=0)
>>> writer.close()

__init__(filename, mode='w', compression='gzip', compression_opts=4, chunk_size=100)

参数:

• filename (str)

• mode (str)

• compression (str | None)

• compression_opts (int | None)

• chunk_size (int)

open()

打开HDF5文件

initialize(n_atoms, n_frames_estimate=1000, atom_types=None, metadata=None)

初始化H5MD文件结构

参数:

• n_atoms (int) -- 原子数量

• n_frames_estimate (int) -- 预估帧数（用于优化存储）

• atom_types (list, optional) -- 原子类型列表

• metadata (dict, optional) -- 额外的元数据

write_frame(positions, box=None, time=None, step=None, velocities=None, forces=None, stress=None, energy=None, temperature=None, **kwargs)

写入单帧数据

参数:

• positions (np.ndarray) -- 原子位置 (n_atoms, 3)

• box (np.ndarray, optional) -- 晶格矢量 (3, 3)

• time (float, optional) -- 时间（ps）

• step (int, optional) -- MD步数

• velocities (np.ndarray, optional) -- 速度 (n_atoms, 3)

• forces (np.ndarray, optional) -- 力 (n_atoms, 3)

• stress (np.ndarray, optional) -- 应力张量 (3, 3)

• energy (float, optional) -- 总能量

• temperature (float, optional) -- 温度

• **kwargs -- 其他观测量

write_metadata(metadata)

写入额外的元数据

参数:

metadata (dict) -- 元数据字典

close()

关闭文件

class thermoelasticsim.utils.trajectory.TrajectoryReader(filename, mode='r', cache_size=100)

基类：object

H5MD格式轨迹文件读取器

提供高效的轨迹数据读取和分析功能。

参数:

• filename (str) -- 输入文件名

• mode (str) -- 读取模式：'r'(只读), 'r+'(读写)

• cache_size (int) -- 内存缓存帧数

示例

>>> reader = TrajectoryReader('trajectory.h5')
>>> reader.open()
>>> for frame_idx in range(reader.n_frames):
...     frame = reader.read_frame(frame_idx)
...     process_frame(frame)
>>> reader.close()

__init__(filename, mode='r', cache_size=100)

参数:

• filename (str)

• mode (str)

• cache_size (int)

open()

打开HDF5文件

read_frame(frame_idx)

读取指定帧

参数:

frame_idx (int) -- 帧索引

返回:

帧数据字典

返回类型:

dict

read_frames(start=0, stop=None, stride=1)

批量读取多帧

参数:

• start (int) -- 起始帧

• stop (int, optional) -- 结束帧（不包含）

• stride (int) -- 步长

返回:

帧数据列表

返回类型:

list

get_metadata()

获取元数据

返回类型:

dict[str, Any]

get_trajectory_info()

获取轨迹概要信息

返回类型:

dict[str, Any]

close()

关闭文件

thermoelasticsim.utils.trajectory.save_trajectory(filename, positions_list, boxes_list=None, times=None, **kwargs)

便捷函数：保存轨迹数据

参数:

• filename (str) -- 输出文件名

• positions_list (list) -- 位置列表

• boxes_list (list, optional) -- 晶格矢量列表

• times (list, optional) -- 时间列表

• **kwargs -- 其他数据

返回类型:

None

thermoelasticsim.utils.trajectory.load_trajectory(filename, start=0, stop=None, stride=1)

便捷函数：加载轨迹数据

参数:

• filename (str) -- 输入文件名

• start (int) -- 起始帧

• stop (int, optional) -- 结束帧

• stride (int) -- 步长

返回:

轨迹数据字典

返回类型:

dict

带轨迹记录的优化器扩展

为现有优化器添加轨迹记录功能，在优化过程中自动捕获原子位置、 晶格矢量、能量、应力等信息并保存到HDF5文件。

Author: Gilbert Young Created: 2025-08-15

class thermoelasticsim.utils.trajectory_recorder.Path(*args, **kwargs)

基类：PurePath

PurePath subclass that can make system calls.

Path represents a filesystem path but unlike PurePath, also offers methods to do system calls on path objects. Depending on your system, instantiating a Path will return either a PosixPath or a WindowsPath object. You can also instantiate a PosixPath or WindowsPath directly, but cannot instantiate a WindowsPath on a POSIX system or vice versa.

absolute()

Return an absolute version of this path by prepending the current working directory. No normalization or symlink resolution is performed.

Use resolve() to get the canonical path to a file.

chmod(mode, *, follow_symlinks=True)

Change the permissions of the path, like os.chmod().

classmethod cwd()

Return a new path pointing to the current working directory (as returned by os.getcwd()).

exists()

Whether this path exists.

expanduser()

Return a new path with expanded ~ and ~user constructs (as returned by os.path.expanduser)

glob(pattern)

Iterate over this subtree and yield all existing files (of any kind, including directories) matching the given relative pattern.

group()

Return the group name of the file gid.

hardlink_to(target)

Make this path a hard link pointing to the same file as target.

Note the order of arguments (self, target) is the reverse of os.link's.

classmethod home()

Return a new path pointing to the user's home directory (as returned by os.path.expanduser('~')).

is_block_device()

Whether this path is a block device.

is_char_device()

Whether this path is a character device.

is_dir()

Whether this path is a directory.

is_fifo()

Whether this path is a FIFO.

is_file()

Whether this path is a regular file (also True for symlinks pointing to regular files).

is_mount()

Check if this path is a POSIX mount point

is_socket()

Whether this path is a socket.

is_symlink()

Whether this path is a symbolic link.

iterdir()

Iterate over the files in this directory. Does not yield any result for the special paths '.' and '..'.

lchmod(mode)

Like chmod(), except if the path points to a symlink, the symlink's permissions are changed, rather than its target's.

link_to(target)

Make the target path a hard link pointing to this path.

Note this function does not make this path a hard link to target, despite the implication of the function and argument names. The order of arguments (target, link) is the reverse of Path.symlink_to, but matches that of os.link.

Deprecated since Python 3.10 and scheduled for removal in Python 3.12. Use hardlink_to() instead.

lstat()

Like stat(), except if the path points to a symlink, the symlink's status information is returned, rather than its target's.

mkdir(mode=511, parents=False, exist_ok=False)

Create a new directory at this given path.

open(mode='r', buffering=-1, encoding=None, errors=None, newline=None)

Open the file pointed by this path and return a file object, as the built-in open() function does.

owner()

Return the login name of the file owner.

read_bytes()

Open the file in bytes mode, read it, and close the file.

read_text(encoding=None, errors=None)

Open the file in text mode, read it, and close the file.

readlink()

Return the path to which the symbolic link points.

rename(target)

Rename this path to the target path.

The target path may be absolute or relative. Relative paths are interpreted relative to the current working directory, not the directory of the Path object.

Returns the new Path instance pointing to the target path.

replace(target)

Rename this path to the target path, overwriting if that path exists.

The target path may be absolute or relative. Relative paths are interpreted relative to the current working directory, not the directory of the Path object.

Returns the new Path instance pointing to the target path.

resolve(strict=False)

Make the path absolute, resolving all symlinks on the way and also normalizing it.

rglob(pattern)

Recursively yield all existing files (of any kind, including directories) matching the given relative pattern, anywhere in this subtree.

rmdir()

Remove this directory. The directory must be empty.

samefile(other_path)

Return whether other_path is the same or not as this file (as returned by os.path.samefile()).

stat(*, follow_symlinks=True)

Return the result of the stat() system call on this path, like os.stat() does.

symlink_to(target, target_is_directory=False)

Make this path a symlink pointing to the target path. Note the order of arguments (link, target) is the reverse of os.symlink.

touch(mode=438, exist_ok=True)

Create this file with the given access mode, if it doesn't exist.

unlink(missing_ok=False)

Remove this file or link. If the path is a directory, use rmdir() instead.

write_bytes(data)

Open the file in bytes mode, write to it, and close the file.

write_text(data, encoding=None, errors=None, newline=None)

Open the file in text mode, write to it, and close the file.

class thermoelasticsim.utils.trajectory_recorder.Any(*args, **kwargs)

基类：object

Special type indicating an unconstrained type.

• Any is compatible with every type.

• Any assumed to have all methods.

• All values assumed to be instances of Any.

Note that all the above statements are true from the point of view of static type checkers. At runtime, Any should not be used with instance checks.

class thermoelasticsim.utils.trajectory_recorder.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)

返回类型:

thermoelasticsim.core.structure.Cell

备注

分数/笛卡尔坐标转换（列向量记号）：

\[\mathbf{r} = \mathbf{L}\,\mathbf{s},\qquad \mathbf{s} = \mathbf{L}^{-1}\,\mathbf{r}\]

实现采用“行向量右乘”的等价式：

\[\mathbf{s}^{\top} = \mathbf{r}^{\top}\,\mathbf{L}^{-1},\qquad \mathbf{r}^{\top} = \mathbf{s}^{\top}\,\mathbf{L}^{\top}\]

示例

创建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} ų")

__init__(lattice_vectors, atoms, pbc_enabled=True)

参数:

• lattice_vectors (ndarray)

• atoms (list[Atom])

• pbc_enabled (bool)

返回类型:

None

apply_deformation(deformation_matrix)

对晶格与原子坐标施加形变矩阵 \(\mathbf{F}\)

参数:

deformation_matrix (numpy.ndarray) -- 形变矩阵 \(\mathbf{F}\)，形状为 (3, 3)

返回类型:

None

备注

采用列向量约定的连续介质力学记号：

• 晶格未锁定（允许变胞）时：

  \[\mathbf{L}_{\text{new}} = \mathbf{F}\, \mathbf{L}_{\text{old}},\qquad \mathbf{r}_{\text{new}} = \mathbf{F}\, \mathbf{r}_{\text{old}}\]

• 晶格已锁定（固定胞形）时：

  \[\mathbf{L}_{\text{new}} = \mathbf{L}_{\text{old}},\qquad \mathbf{r}_{\text{new}} = \mathbf{F}\, \mathbf{r}_{\text{old}}\]

实现细节：本实现使用“行向量右乘”形式，等价写法为

\[\mathbf{r}_{\text{new}}^{\top} = \mathbf{r}_{\text{old}}^{\top} \mathbf{F}^{\top}.\]

示例

>>> import numpy as np
>>> F = np.array([[1.01, 0, 0], [0, 1, 0], [0, 0, 1]])
>>> cell.apply_deformation(F)

apply_periodic_boundary(positions)

应用周期性边界条件

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

算法（列向量记号）：

\[\mathbf{s} = \mathbf{L}^{-1}\,\mathbf{r},\qquad \mathbf{s}' = \mathbf{s} - \lfloor \mathbf{s} + 0.5 \rfloor,\qquad \mathbf{r}' = \mathbf{L}\,\mathbf{s}'\]

参数:

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)

build_supercell(repetition)

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

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

参数:

repetition (tuple of int) -- 在 x, y, z 方向上的重复次数，例如 (2, 2, 2)

返回:

新的超胞对象

返回类型:

Cell

calculate_kinetic_energy()

计算当前系统的总动能

返回:

系统总动能，单位为 \(\mathrm{eV}\)

返回类型:

float

备注

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

示例

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

calculate_stress_tensor(potential)

计算应力张量

参数:

potential (Potential) -- 用于计算应力的势能对象

返回:

3×3 应力张量矩阵 (eV/ų)

返回类型:

numpy.ndarray

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\]

自由度数定义：

\[\begin{split}N_{\mathrm{dof}} = \begin{cases} 3N - 3, & N > 1, \\ 3, & N = 1. \end{cases}\end{split}\]

返回:

系统温度，单位为 \(\mathrm{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_volume()

计算晶胞的体积

返回:

晶胞体积 (ų)

返回类型:

float

copy()

创建 Cell 的深拷贝

get_box_lengths()

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

返回:

包含三个方向长度的数组 (Å)

返回类型:

numpy.ndarray

get_com_position()

计算系统质心位置。

返回:

质心位置向量

返回类型:

numpy.ndarray

get_com_velocity()

计算系统质心速度。

返回:

质心速度向量

返回类型:

numpy.ndarray

get_forces()

获取所有原子的力信息

返回:

原子力数组，形状为 (num_atoms, 3)

返回类型:

numpy.ndarray

get_fractional_coordinates()

获取所有原子的分数坐标

返回:

原子分数坐标数组，形状为(num_atoms, 3)

返回类型:

numpy.ndarray

备注

坐标关系（列向量记号）：

\[\mathbf{s} = \mathbf{L}^{-1}\,\mathbf{r}\]

代码实现采用行向量右乘：\(\mathbf{s}^{\top} = \mathbf{r}^{\top}\,\mathbf{L}^{-1}\)。

get_positions()

获取所有原子的位置信息

返回:

原子位置数组，形状为 (num_atoms, 3)

返回类型:

numpy.ndarray

示例

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

get_velocities()

获取所有原子的速度信息

返回:

原子速度数组，形状为 (num_atoms, 3)

返回类型:

numpy.ndarray

get_volume()

计算晶胞体积

返回:

晶胞体积 (ų)

返回类型:

float

备注

体积通过晶格矢量的标量三重积计算：

\[V = \left|\, \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) \,\right|\]

lock_lattice_vectors()

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

备注

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

返回类型:

None

minimum_image(displacement)

计算最小镜像位移向量

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

数学原理（列向量记号）：

\[\mathbf{d}_{\min} = \mathbf{d} - \mathbf{L}\, \operatorname{round}(\mathbf{L}^{-1} \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

应用周期性边界条件

property num_atoms

返回原子数量

remove_com_motion()

移除系统的质心运动。

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

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

set_fractional_coordinates(fractional_coords)

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

参数:

fractional_coords (numpy.ndarray) -- 分数坐标数组，形状为(num_atoms, 3)

返回类型:

None

备注

坐标关系（列向量记号）：

\[\mathbf{r} = \mathbf{L}\,\mathbf{s}\]

代码实现采用行向量右乘：\(\mathbf{r}^{\top} = \mathbf{s}^{\top}\,\mathbf{L}^{\top}\)。

set_lattice_vectors(lattice_vectors)

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

参数:

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

set_positions(positions)

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

参数:

positions (numpy.ndarray) -- 笛卡尔坐标数组，形状为(num_atoms, 3)

返回类型:

None

set_velocities(velocities)

设置所有原子的速度

参数:

velocities (numpy.ndarray) -- 速度数组，形状为 (num_atoms, 3)

抛出:

ValueError -- 当输入数组形状与原子数不匹配时抛出

返回类型:

None

unlock_lattice_vectors()

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

备注

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

返回类型:

None

class thermoelasticsim.utils.trajectory_recorder.Potential(parameters, cutoff)

基类：ABC

势能计算的抽象基类。

提供原子间势能模型的统一接口，约定最小方法集以支持分子动力学 (MD) 中的力与能量评估。

参数:

• parameters (dict) -- 势能参数字典（按具体模型定义）。

• cutoff (float) -- 势能截断距离，单位 Å。

备注

• 所有具体势能模型至少需实现 calculate_forces 与 calculate_energy。

• 单位约定：能量 eV，长度 Å，力 eV/Å。

__init__(parameters, cutoff)

abstractmethod calculate_energy(cell, neighbor_list)

计算系统的总势能。

参数:

• cell (Cell) -- 晶胞与原子集合。

• neighbor_list (NeighborList) -- 邻居列表（按模型需要使用）。

返回:

系统总势能，单位 eV。

返回类型:

float

备注

抽象方法，必须由子类实现。

abstractmethod calculate_forces(cell, neighbor_list)

计算系统中所有原子的作用力。

参数:

• cell (Cell) -- 晶胞与原子集合。

• neighbor_list (NeighborList) -- 邻居列表（按模型需要使用）。

返回类型:

None

备注

• 需就地写入 force （单位 eV/Å）。

• 抽象方法，必须由子类实现。

class thermoelasticsim.utils.trajectory_recorder.LBFGSOptimizer(ftol=1e-08, gtol=1e-06, maxiter=1000, supercell_dims=None, **kwargs)

基类：Optimizer

L-BFGS优化器(有限内存BFGS)

该优化器能够处理两种模式： 1. 仅优化原子位置 (默认)。 2. 同时优化原子位置和晶胞形状 (完全弛豫)。

参数:

• ftol (float, optional) -- 函数收敛阈值，默认1e-8

• gtol (float, optional) -- 梯度收敛阈值，默认1e-6

• maxiter (int, optional) -- 最大迭代步数，默认1000

• **kwargs (dict, optional) -- 其他要传递给 scipy.optimize.minimize 的选项, 例如 {'disp': True}

__init__(ftol=1e-08, gtol=1e-06, maxiter=1000, supercell_dims=None, **kwargs)

optimize(cell, potential, relax_cell=False)

执行 L-BFGS 优化。

参数:

• cell (Cell)

• potential (Potential)

• relax_cell (bool)

返回类型:

tuple[bool, list[dict]]

class thermoelasticsim.utils.trajectory_recorder.TrajectoryWriter(filename, mode='w', compression='gzip', compression_opts=4, chunk_size=100)

基类：object

H5MD格式轨迹文件写入器

遵循H5MD v1.1规范，提供高效的分子动力学轨迹存储。

参数:

• filename (str) -- 输出文件名（.h5或.hdf5扩展名）

• mode (str) -- 文件打开模式：'w'(覆盖), 'a'(追加), 'x'(创建新文件)

• compression (str, optional) -- 压缩算法：'gzip', 'lzf', 'szip'等

• compression_opts (int, optional) -- 压缩级别（gzip: 1-9）

• chunk_size (int, optional) -- 块大小，影响I/O性能

示例

>>> writer = TrajectoryWriter('trajectory.h5', compression='gzip')
>>> writer.initialize(n_atoms=100, n_frames_estimate=1000)
>>> writer.write_frame(positions, box, time=0.0, step=0)
>>> writer.close()

__init__(filename, mode='w', compression='gzip', compression_opts=4, chunk_size=100)

参数:

• filename (str)

• mode (str)

• compression (str | None)

• compression_opts (int | None)

• chunk_size (int)

close()

关闭文件

initialize(n_atoms, n_frames_estimate=1000, atom_types=None, metadata=None)

初始化H5MD文件结构

参数:

• n_atoms (int) -- 原子数量

• n_frames_estimate (int) -- 预估帧数（用于优化存储）

• atom_types (list, optional) -- 原子类型列表

• metadata (dict, optional) -- 额外的元数据

open()

打开HDF5文件

write_frame(positions, box=None, time=None, step=None, velocities=None, forces=None, stress=None, energy=None, temperature=None, **kwargs)

写入单帧数据

参数:

• positions (np.ndarray) -- 原子位置 (n_atoms, 3)

• box (np.ndarray, optional) -- 晶格矢量 (3, 3)

• time (float, optional) -- 时间（ps）

• step (int, optional) -- MD步数

• velocities (np.ndarray, optional) -- 速度 (n_atoms, 3)

• forces (np.ndarray, optional) -- 力 (n_atoms, 3)

• stress (np.ndarray, optional) -- 应力张量 (3, 3)

• energy (float, optional) -- 总能量

• temperature (float, optional) -- 温度

• **kwargs -- 其他观测量

write_metadata(metadata)

写入额外的元数据

参数:

metadata (dict) -- 元数据字典

class thermoelasticsim.utils.trajectory_recorder.TrajectoryRecorder(output_file, record_interval=1, record_forces=True, record_stress=True, record_energy=True, compression='gzip', compression_opts=4)

基类：object

轨迹记录器

可以作为回调函数集成到优化器中，自动记录优化轨迹。

参数:

• output_file (str) -- 输出轨迹文件名

• record_interval (int) -- 记录间隔（每N步记录一次）

• record_forces (bool) -- 是否记录力

• record_stress (bool) -- 是否记录应力

• record_energy (bool) -- 是否记录能量

• compression (str) -- 压缩算法

• compression_opts (int)

示例

>>> recorder = TrajectoryRecorder('optimization.h5')
>>> optimizer = LBFGSOptimizerWithTrajectory(recorder=recorder)
>>> optimizer.optimize(cell, potential)

__init__(output_file, record_interval=1, record_forces=True, record_stress=True, record_energy=True, compression='gzip', compression_opts=4)

参数:

• output_file (str)

• record_interval (int)

• record_forces (bool)

• record_stress (bool)

• record_energy (bool)

• compression (str)

• compression_opts (int)

initialize(cell, metadata=None)

初始化记录器

参数:

• cell (Cell) -- 初始晶胞

• metadata (dict, optional) -- 元数据

record(cell, potential, step=None)

记录当前状态

参数:

• cell (Cell) -- 当前晶胞

• potential (Potential) -- 势能函数

• step (int, optional) -- 优化步数

finalize()

完成记录并关闭文件

class thermoelasticsim.utils.trajectory_recorder.LBFGSOptimizerWithTrajectory(recorder=None, trajectory_file=None, record_interval=1, **kwargs)

基类：LBFGSOptimizer

带轨迹记录功能的L-BFGS优化器

在原有L-BFGS优化器基础上添加轨迹记录功能。

参数:

• recorder (TrajectoryRecorder, optional) -- 轨迹记录器

• trajectory_file (str, optional) -- 轨迹文件名（如果不提供recorder）

• record_interval (int) -- 记录间隔

• **kwargs -- 其他L-BFGS参数

__init__(recorder=None, trajectory_file=None, record_interval=1, **kwargs)

参数:

• recorder (TrajectoryRecorder | None)

• trajectory_file (str | None)

• record_interval (int)

optimize(cell, potential, relax_cell=False)

执行优化并记录轨迹

参数:

• cell (Cell) -- 晶胞对象

• potential (Potential) -- 势能函数

• relax_cell (bool) -- 是否优化晶格

返回:

(是否收敛, 优化信息字典)

返回类型:

tuple

thermoelasticsim.utils.trajectory_recorder.create_optimizer_with_trajectory(optimizer_type='L-BFGS', trajectory_file=None, record_interval=1, **optimizer_params)

创建带轨迹记录的优化器

参数:

• optimizer_type (str) -- 优化器类型

• trajectory_file (str, optional) -- 轨迹文件名

• record_interval (int) -- 记录间隔

• **optimizer_params -- 优化器参数

返回:

带轨迹记录的优化器实例

返回类型:

optimizer

class thermoelasticsim.utils.trajectory_recorder.DeformationTrajectoryRecorder(output_file)

基类：object

专门用于形变计算的轨迹记录器

记录形变过程中的完整信息，包括应变、应力、能量等。

参数:

output_file (str)

__init__(output_file)

参数:

output_file (str)

initialize(base_cell, metadata=None)

初始化

参数:

• base_cell (Cell)

• metadata (dict[str, Any] | None)

record_deformation(cell, potential, strain, stress, deformation_matrix, mode, converged)

记录形变状态

参数:

• cell (Cell) -- 形变后的晶胞

• potential (Potential) -- 势能函数

• strain (np.ndarray) -- 应变张量

• stress (np.ndarray) -- 应力张量

• deformation_matrix (np.ndarray) -- 形变矩阵

• mode (str) -- 形变模式

• converged (bool) -- 是否收敛

finalize()

完成记录

绘图与可视化（可选）

备注

可视化相关模块依赖较多第三方库，且与文档构建关系不大。为降低构建复杂度，本节默认不自动展开 API。需要详细说明时，可在本地环境中查看源码注释或启用相应模块文档。