分子动力学引擎 (md)
本模块实现了基于算符分离的分子动力学引擎。
接口定义
MD 算符分离架构接口定义
本模块定义了基于刘维尔算符 Trotter 分解的 MD 架构接口； 核心思想是将完整积分步骤拆解为可组合的基础算符。
主要类
- MDComponent
MD 组件基类
- Propagator
算符传播子基类，对应 \(\exp(iL\,\Delta t)\)
- IntegrationScheme
积分方案基类，通过 Trotter 分解组合多个传播子
参考文献
Martyna et al., Mol. Phys. 87, 1117 (1996); Tuckerman et al., J. Chem. Phys. 97, 1990 (1992); Martyna et al., J. Chem. Phys. 97, 2635 (1992); Yoshida, Phys. Lett. A 150, 262 (1990).
- class thermoelasticsim.md.interfaces.Propagator[源代码]
基类：
MDComponent算符传播子基类。
在特定物理过程下更新系统状态（如位置、速度、热浴变量等）。
- class thermoelasticsim.md.interfaces.IntegrationScheme[源代码]
基类：
MDComponent积分方案基类。
通过 Trotter 分解组合多个传播子，实现完整积分步骤。
- class thermoelasticsim.md.interfaces.ThermostatInterface[源代码]
基类：
MDComponent恒温器接口。
传播子
基础算符传播子实现
每个传播子对应一个基本物理过程，可通过对称 Trotter 分解组合为完整积分方案。
已实现的传播子
- PositionPropagator
位置演化：\(\mathbf{r} \leftarrow \mathbf{r} + \mathbf{v}\,\Delta t\)
- VelocityPropagator
速度演化：\(\mathbf{v} \leftarrow \mathbf{v} + \mathbf{F}\,\Delta t/m\)
- ForcePropagator
力计算：\(\mathbf{F} = -\nabla U(\mathbf{r})\)
设计要点
单一职责：每个传播子仅负责一种过程
与 Velocity-Verlet 分解一致，减少重复力计算
数值稳定、支持时间可逆的对称分解
- class thermoelasticsim.md.propagators.StressCalculator[源代码]
基类：
object应力张量计算器
总应力由两部分组成（体积 \(V\)）：
- 动能项（动量通量）
- \[\sigma^{\text{kin}}_{\alpha\beta} = -\,\frac{1}{V} \sum_i m_i\, v_{i\alpha} v_{i\beta}\]
- 维里项（力-位矢项）
- \[\sigma^{\text{vir}}_{\alpha\beta} = -\,\frac{1}{V} \sum_{i<j} r_{ij,\alpha}\, F_{ij,\beta}\]
备注
有限差分/晶格导数法（如 \(-\partial U/\partial \varepsilon\) 或 \(-V^{-1}\,\partial U/\partial \mathbf{h}\, \mathbf{h}^T\)）是另一种等价 的应力计算方式，用于数值校验；它并非额外“第三项”，不与维里项相加。
本实现采用 \(\sigma = \sigma^{\text{kin}} + \sigma^{\text{vir}}\)。
对 EAM 势，使用经特殊处理的多体维里解析形式；如解析不可用，可用有限差分做校验。
正负号及指标约定与项目其它模块保持一致。
- calculate_kinetic_stress(cell)[源代码]
计算动能应力张量
使用的约定：
\[\sigma^{\text{kin}}_{\alpha\beta} = -\,\frac{1}{V} \sum_i m_i\, v_{i\alpha} v_{i\beta}.\]说明：在静止构型或 \(T\to 0\) 极限，速度趋近于零，因此该项趋近于 0。
- 参数:
cell (Cell) -- 包含原子的晶胞对象
- 返回:
动能应力张量 (3, 3)，单位 eV/ų
- 返回类型:
- calculate_total_stress(cell, potential)[源代码]
计算总应力张量（动能项 + 维里项）
\[\sigma_{\alpha\beta} = -\,\frac{1}{V} \left( \sum_i m_i v_{i\alpha} v_{i\beta} + \sum_{i<j} r_{ij,\alpha} F_{ij,\beta} \right)\]
- calculate_virial_stress(cell, potential)[源代码]
计算维里应力张量（相互作用力贡献）
公式：
\[\sigma^{\text{vir}}_{\alpha\beta} = -\,\frac{1}{V} \sum_{i<j} r_{ij,\alpha}\, F_{ij,\beta}\]其中：
- \(\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j\)
从原子 \(j\) 指向原子 \(i\) 的位移向量
- \(\mathbf{F}_{ij}\)
原子 \(j\) 作用于 \(i\) 的力
- class thermoelasticsim.md.propagators.Propagator[源代码]
基类：
MDComponent算符传播子基类。
在特定物理过程下更新系统状态（如位置、速度、热浴变量等）。
- class thermoelasticsim.md.propagators.PositionPropagator[源代码]
基类：
Propagator位置传播子：\(\exp(iL_r\,\Delta t)\)
实现位置的时间演化：\(\mathbf{r}(t+\Delta t)=\mathbf{r}(t)+\mathbf{v}(t)\,\Delta t\)。
备注
处理周期性边界；不改变速度与力；与势能无关。
- class thermoelasticsim.md.propagators.VelocityPropagator[源代码]
基类：
Propagator速度传播子：\(\exp(iL_v\,\Delta t)\)
实现速度的时间演化：\(\mathbf{v}(t+\Delta t)=\mathbf{v}(t)+\mathbf{F}(t)\,\Delta t/m\)。
备注
假设
atom.force已计算；不在此处重新计算力。
- class thermoelasticsim.md.propagators.ForcePropagator(potential)[源代码]
基类：
Propagator力计算传播子：\(\mathbf{F} = -\nabla U(\mathbf{r})\)
封装势函数的力计算调用，控制力计算时机，避免重复计算。
- class thermoelasticsim.md.propagators.BerendsenThermostatPropagator(target_temperature, tau=None, mode='equilibration')[源代码]
基类：
PropagatorBerendsen 恒温器传播子（速度缩放）
缩放因子 \(\lambda\)：
\[\lambda = \sqrt{1 + \frac{\Delta t}{\tau}\left(\frac{T_0}{T}-1\right)}\]备注
后处理式恒温；不严格产生正则分布，常用于快速平衡。
引用
Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. (1984). Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics, 81(8), 3684–3690. https://doi.org/10.1063/1.448118
- propagate(cell, dt, **kwargs)[源代码]
执行 Berendsen 温度调节
备注
标准流程：
计算瞬时温度 \(T\)
计算缩放因子 \(\lambda = \sqrt{1 + \tfrac{\Delta t}{\tau_T}(\tfrac{T_0}{T}-1)}\)
对所有原子速度缩放 \(\mathbf{v} \leftarrow \lambda\,\mathbf{v}\)
可选：数值稳定性限制与统计监控
- get_statistics()[源代码]
获取 Berendsen 恒温器统计信息
- 返回:
统计信息字典：
- total_steps
总调节步数
- average_scaling
平均缩放因子
- max_scaling
最大缩放因子
- min_scaling
最小缩放因子
- target_temperature
目标温度 (K)
- tau
时间常数 \(\tau_T\) (fs)
- mode
运行模式（equilibration/production）
- limited_steps
被限制的步数
- upper_limited / lower_limited
上/下限限制次数
- limitation_rate
限制比例 (%)
- effective_coupling
有效耦合强度评估
- 返回类型:
- class thermoelasticsim.md.propagators.AndersenThermostatPropagator(target_temperature, collision_frequency=None)[源代码]
基类：
PropagatorAndersen 随机碰撞恒温器传播子
每个原子在每个时间步以概率 \(p=\nu\,\Delta t\) 与热浴碰撞， 并从 Maxwell–Boltzmann 分布重新采样速度。
备注
产生正确的正则 (NVT) 分布，但会打断动力学连续性。
引用
Andersen, H. C. (1980). Molecular dynamics simulations at constant pressure and/or temperature. The Journal of Chemical Physics, 72(4), 2384–2393. https://doi.org/10.1063/1.439486
- __init__(target_temperature, collision_frequency=None)[源代码]
初始化 Andersen 恒温器
- 参数:
- 抛出:
ValueError -- 如果target_temperature或collision_frequency为非正数
- propagate(cell, dt, **kwargs)[源代码]
执行 Andersen 随机碰撞
备注
计算系统级碰撞概率 \(P=\nu\,\Delta t\)
逐原子以概率 \(P/N\) 判定是否碰撞
若碰撞，从 Maxwell–Boltzmann 分布采样新速度
- get_statistics()[源代码]
获取 Andersen 恒温器统计信息
- 返回:
统计信息字典：
- total_steps
总步数
- total_collisions
总碰撞次数
- collision_rate
平均碰撞率 (次/步)
- expected_rate
理论碰撞率 (fs⁻¹)
- target_temperature
目标温度 (K)
- effective_collision_frequency
有效碰撞频率 (fs⁻¹)
- base_collision_frequency
基础频率 (fs⁻¹)
- recent_collisions
最近 100 步的碰撞数
- system_size_initialized
是否已根据系统尺寸初始化
- 返回类型:
- class thermoelasticsim.md.propagators.NoseHooverChainPropagator(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
基类：
PropagatorNose–Hoover 链恒温器传播子（四阶 Suzuki–Yoshida 分解）
备注
链式热浴变量自洽演化，时间可逆对称分解，保证正则系综采样质量。
- 参考文献：
Martyna et al., J. Chem. Phys. 97, 2635 (1992); Yoshida, Phys. Lett. A 150, 262 (1990).
- __init__(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
初始化Nose-Hoover链恒温器
- 参数:
- 抛出:
ValueError -- 如果参数设置不合理
- propagate(cell, dt, **kwargs)[源代码]
执行 Nose–Hoover 链传播（四阶 Suzuki–Yoshida）
实现完整的 NHC 算符：
\[\exp\big(iL_{\mathrm{NHC}}\,\Delta t\big)\]采用三步四阶 Suzuki–Yoshida 分解保证高精度和长期稳定性。
算法流程：
遍历每个 SY 系数 \(w_k\)
对每个系数执行 \(t_{loop}\) 次内循环
调用
_nhc_integration_loop(w_k cdot Delta t / t_{loop})
- 参数:
- 返回类型:
None
备注
与 NVE 的对称包裹关系：
\[\exp\big(iL_{\mathrm{NHC}}\tfrac{\Delta t}{2}\big)\,\exp\big(iL_{\mathrm{NVE}}\,\Delta t\big)\,\exp\big(iL_{\mathrm{NHC}}\tfrac{\Delta t}{2}\big)\]
- thermoelasticsim.md.propagators.create_nve_propagators(potential)[源代码]
创建 NVE 积分所需的基础传播子
- 参数:
potential (Potential) -- 势函数对象
- 返回:
传播子字典：
'position''velocity''force'
- 返回类型:
示例
>>> propagators = create_nve_propagators(eam_potential) >>> pos_prop = propagators['position'] >>> vel_prop = propagators['velocity'] >>> force_prop = propagators['force']
- class thermoelasticsim.md.propagators.LangevinThermostatPropagator(target_temperature, friction=1.0)[源代码]
基类：
PropagatorLangevin 动力学恒温器传播子（BBK 积分）
通过阻尼力和随机力模拟热浴作用，实现温度控制。
备注
Langevin 运动方程：
\[m_i \frac{d^2\mathbf{r}_i}{dt^2} = \mathbf{F}_i(\mathbf{r}) - \gamma m_i \frac{d\mathbf{r}_i}{dt} + \mathbf{R}_i(t)\]其中：
\(\mathbf{F}_i(\mathbf{r})\) - 保守力
\(\gamma\) - 摩擦系数（friction coefficient）
\(\mathbf{R}_i(t)\) - 随机力，满足涨落–耗散定理
BBK 积分常数：
\(c_1 = \exp(-\gamma\,\Delta t)\) - 速度阻尼因子
\(\sigma = \sqrt{k_B T\,(1-c_1^2)/m}\) - 随机力标准差
- 参考文献：
Brünger, Brooks & Karplus, Chem. Phys. Lett. 105, 495 (1984).
- __init__(target_temperature, friction=1.0)[源代码]
初始化 Langevin 恒温器。
- 参数:
- 抛出:
ValueError -- 如果 target_temperature 或 friction 为非正数。
- propagate(cell, dt, **kwargs)[源代码]
执行 Langevin 恒温器传播（BBK 速度更新）。
- 参数:
- 抛出:
ValueError -- 如果 dt <= 0。
RuntimeError -- 如果温度计算或随机数生成失败。
- 返回类型:
None
- get_statistics()[源代码]
获取 Langevin 恒温器统计信息
- 返回:
统计信息字典：
- total_steps
总传播步数
- target_temperature
目标温度 (K)
- friction
摩擦系数 (ps⁻¹)
- damping_time
阻尼时间 (ps)
- coupling_strength
耦合强度描述
- average_temperature
平均温度 (K)
- temperature_std
温度标准差 (K)
- temperature_error
平均温度误差 (%)
- average_friction_work
平均摩擦做功 (eV)
- average_random_work
平均随机做功 (eV)
- energy_balance
摩擦做功与随机做功比值
- recent_temperatures
最近 100 步的温度历史
- 返回类型:
- set_friction(new_friction)[源代码]
动态调整摩擦系数
- 参数:
new_friction (float) -- 新的摩擦系数 (ps⁻¹)
- 抛出:
ValueError -- 如果new_friction为非正数
- 返回类型:
None
备注
这允许在模拟过程中调整恒温器的耦合强度， 例如在平衡阶段使用大摩擦系数，在生产阶段使用小摩擦系数。
- thermoelasticsim.md.propagators.exprel(x)[源代码]
计算 (exp(x) - 1) / x，避免x接近0时的数值误差
- 参数:
x (np.ndarray) -- 输入数组
- 返回:
计算结果
- 返回类型:
np.ndarray
- class thermoelasticsim.md.propagators.MTKBarostatPropagator(target_temperature, target_pressure, pdamp, pchain=3, ploop=1)[源代码]
基类：
PropagatorMTK（Martyna–Tobias–Klein）恒压器传播子
备注
对称分解 + 矩阵指数实现：
晶格共轭动量 \(\mathbf{P}_g\) 的演化由应力与外压驱动
在 \(\mathbf{P}_g\) 的特征基中进行精确积分（矩阵指数）
核心关系：
\[d\mathbf{P}_g/dt = V\,(\boldsymbol{\sigma}_{\mathrm{int}} - P_{\mathrm{ext}}\,\mathbf{I}) + \text{kinetic correction}\]引用
Martyna, G. J.; Tobias, D. J.; Klein, M. L. (1994). Constant pressure molecular dynamics algorithms. The Journal of Chemical Physics, 101(5), 4177–4189. https://doi.org/10.1063/1.467468
- 参数:
- propagate(cell, dt, potential=None)[源代码]
执行恒压器传播一个时间步dt
基于对称Trotter分解的恒压器链积分，包含： 1. 恒压器链传播 (dt/2) 2. 晶格动量更新 3. 恒压器链传播 (dt/2)
- integrate_momenta(cell, potential, dt)[源代码]
更新粒子动量（与恒压器耦合的半步，MTK 公式）
在 \(\mathbf{P}_g\) 的特征基下：
\[\mathbf{p}' = \mathbf{p}\,\exp\!\Big(-\frac{\boldsymbol{\kappa}\,\Delta t}{W}\Big) + \Delta t\,\mathbf{F}\,\operatorname{exprel}\!\Big(-\frac{\boldsymbol{\kappa}\,\Delta t}{W}\Big)\]其中 \(\boldsymbol{\kappa} = \boldsymbol{\lambda} + \operatorname{Tr}(\mathbf{P}_g)/(3N)\)。
备注
此方法已包含力贡献，无需额外的 NVE 速度半步
调用前后需配合
_update_cell_momenta(dt/2)与位置/晶格传播
- propagate_positions_and_cell(cell, dt)[源代码]
使用矩阵指数精确传播粒子位置与晶格（MTK 核心步骤）
流程：
特征分解 \(\mathbf{P}_g = \mathbf{U}\,\operatorname{diag}(\boldsymbol{\lambda})\,\mathbf{U}^\top\)
在特征坐标系中用 \(\exp(\boldsymbol{\lambda}\,\Delta t/W)\) 精确积分
变换回原坐标系
积分方案
积分方案实现
基于刘维尔算符的对称 Trotter 分解，组合基础传播子形成不同系综的积分方案。
已实现的方案
- NVEScheme
微正则系综的 Velocity-Verlet 算法实现。
设计要点
基于对称 Trotter 分解保持时间可逆性与稳定性
复用基础传播子，避免重复力计算
每个方案保证其对应系综的物理正确性
支持延迟初始化以处理势函数依赖
- class thermoelasticsim.md.schemes.IntegrationScheme[源代码]
基类：
MDComponent积分方案基类。
通过 Trotter 分解组合多个传播子，实现完整积分步骤。
- class thermoelasticsim.md.schemes.AndersenThermostatPropagator(target_temperature, collision_frequency=None)[源代码]
基类：
PropagatorAndersen 随机碰撞恒温器传播子
每个原子在每个时间步以概率 \(p=\nu\,\Delta t\) 与热浴碰撞， 并从 Maxwell–Boltzmann 分布重新采样速度。
备注
产生正确的正则 (NVT) 分布，但会打断动力学连续性。
引用
Andersen, H. C. (1980). Molecular dynamics simulations at constant pressure and/or temperature. The Journal of Chemical Physics, 72(4), 2384–2393. https://doi.org/10.1063/1.439486
- __init__(target_temperature, collision_frequency=None)[源代码]
初始化 Andersen 恒温器
- 参数:
- 抛出:
ValueError -- 如果target_temperature或collision_frequency为非正数
- get_statistics()[源代码]
获取 Andersen 恒温器统计信息
- 返回:
统计信息字典：
- total_steps
总步数
- total_collisions
总碰撞次数
- collision_rate
平均碰撞率 (次/步)
- expected_rate
理论碰撞率 (fs⁻¹)
- target_temperature
目标温度 (K)
- effective_collision_frequency
有效碰撞频率 (fs⁻¹)
- base_collision_frequency
基础频率 (fs⁻¹)
- recent_collisions
最近 100 步的碰撞数
- system_size_initialized
是否已根据系统尺寸初始化
- 返回类型:
- class thermoelasticsim.md.schemes.BerendsenThermostatPropagator(target_temperature, tau=None, mode='equilibration')[源代码]
基类：
PropagatorBerendsen 恒温器传播子（速度缩放）
缩放因子 \(\lambda\)：
\[\lambda = \sqrt{1 + \frac{\Delta t}{\tau}\left(\frac{T_0}{T}-1\right)}\]备注
后处理式恒温；不严格产生正则分布，常用于快速平衡。
引用
Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. (1984). Molecular dynamics with coupling to an external bath. The Journal of Chemical Physics, 81(8), 3684–3690. https://doi.org/10.1063/1.448118
- get_statistics()[源代码]
获取 Berendsen 恒温器统计信息
- 返回:
统计信息字典：
- total_steps
总调节步数
- average_scaling
平均缩放因子
- max_scaling
最大缩放因子
- min_scaling
最小缩放因子
- target_temperature
目标温度 (K)
- tau
时间常数 \(\tau_T\) (fs)
- mode
运行模式（equilibration/production）
- limited_steps
被限制的步数
- upper_limited / lower_limited
上/下限限制次数
- limitation_rate
限制比例 (%)
- effective_coupling
有效耦合强度评估
- 返回类型:
- class thermoelasticsim.md.schemes.ForcePropagator(potential)[源代码]
基类：
Propagator力计算传播子：\(\mathbf{F} = -\nabla U(\mathbf{r})\)
封装势函数的力计算调用，控制力计算时机，避免重复计算。
- class thermoelasticsim.md.schemes.LangevinThermostatPropagator(target_temperature, friction=1.0)[源代码]
基类：
PropagatorLangevin 动力学恒温器传播子（BBK 积分）
通过阻尼力和随机力模拟热浴作用，实现温度控制。
备注
Langevin 运动方程：
\[m_i \frac{d^2\mathbf{r}_i}{dt^2} = \mathbf{F}_i(\mathbf{r}) - \gamma m_i \frac{d\mathbf{r}_i}{dt} + \mathbf{R}_i(t)\]其中：
\(\mathbf{F}_i(\mathbf{r})\) - 保守力
\(\gamma\) - 摩擦系数（friction coefficient）
\(\mathbf{R}_i(t)\) - 随机力，满足涨落–耗散定理
BBK 积分常数：
\(c_1 = \exp(-\gamma\,\Delta t)\) - 速度阻尼因子
\(\sigma = \sqrt{k_B T\,(1-c_1^2)/m}\) - 随机力标准差
- 参考文献：
Brünger, Brooks & Karplus, Chem. Phys. Lett. 105, 495 (1984).
- __init__(target_temperature, friction=1.0)[源代码]
初始化 Langevin 恒温器。
- 参数:
- 抛出:
ValueError -- 如果 target_temperature 或 friction 为非正数。
- get_effective_parameters()[源代码]
获取当前有效参数
- 返回:
参数字典：
- friction
摩擦系数 (ps⁻¹)
- damping_time
阻尼时间 (ps)，\(\tau_{\mathrm{damp}} = 1/\gamma\)
- target_temperature
目标温度 (K)
- coupling_description
耦合强度描述（weak/moderate/strong）
- 返回类型:
- get_statistics()[源代码]
获取 Langevin 恒温器统计信息
- 返回:
统计信息字典：
- total_steps
总传播步数
- target_temperature
目标温度 (K)
- friction
摩擦系数 (ps⁻¹)
- damping_time
阻尼时间 (ps)
- coupling_strength
耦合强度描述
- average_temperature
平均温度 (K)
- temperature_std
温度标准差 (K)
- temperature_error
平均温度误差 (%)
- average_friction_work
平均摩擦做功 (eV)
- average_random_work
平均随机做功 (eV)
- energy_balance
摩擦做功与随机做功比值
- recent_temperatures
最近 100 步的温度历史
- 返回类型:
- propagate(cell, dt, **kwargs)[源代码]
执行 Langevin 恒温器传播（BBK 速度更新）。
- 参数:
- 抛出:
ValueError -- 如果 dt <= 0。
RuntimeError -- 如果温度计算或随机数生成失败。
- 返回类型:
None
- set_friction(new_friction)[源代码]
动态调整摩擦系数
- 参数:
new_friction (float) -- 新的摩擦系数 (ps⁻¹)
- 抛出:
ValueError -- 如果new_friction为非正数
- 返回类型:
None
备注
这允许在模拟过程中调整恒温器的耦合强度， 例如在平衡阶段使用大摩擦系数，在生产阶段使用小摩擦系数。
- class thermoelasticsim.md.schemes.MTKBarostatPropagator(target_temperature, target_pressure, pdamp, pchain=3, ploop=1)[源代码]
基类：
PropagatorMTK（Martyna–Tobias–Klein）恒压器传播子
备注
对称分解 + 矩阵指数实现：
晶格共轭动量 \(\mathbf{P}_g\) 的演化由应力与外压驱动
在 \(\mathbf{P}_g\) 的特征基中进行精确积分（矩阵指数）
核心关系：
\[d\mathbf{P}_g/dt = V\,(\boldsymbol{\sigma}_{\mathrm{int}} - P_{\mathrm{ext}}\,\mathbf{I}) + \text{kinetic correction}\]引用
Martyna, G. J.; Tobias, D. J.; Klein, M. L. (1994). Constant pressure molecular dynamics algorithms. The Journal of Chemical Physics, 101(5), 4177–4189. https://doi.org/10.1063/1.467468
- 参数:
- get_barostat_energy()[源代码]
计算恒压器贡献的能量
\[E_{\mathrm{baro}} = \sum_j \frac{p_{\xi_j}^2}{2 R_j} + \frac{\operatorname{Tr}(\mathbf{P}_g^{\top}\mathbf{P}_g)}{2W} + \text{external work}\]- 返回:
恒压器能量 (eV)
- 返回类型:
- integrate_momenta(cell, potential, dt)[源代码]
更新粒子动量（与恒压器耦合的半步，MTK 公式）
在 \(\mathbf{P}_g\) 的特征基下：
\[\mathbf{p}' = \mathbf{p}\,\exp\!\Big(-\frac{\boldsymbol{\kappa}\,\Delta t}{W}\Big) + \Delta t\,\mathbf{F}\,\operatorname{exprel}\!\Big(-\frac{\boldsymbol{\kappa}\,\Delta t}{W}\Big)\]其中 \(\boldsymbol{\kappa} = \boldsymbol{\lambda} + \operatorname{Tr}(\mathbf{P}_g)/(3N)\)。
备注
此方法已包含力贡献，无需额外的 NVE 速度半步
调用前后需配合
_update_cell_momenta(dt/2)与位置/晶格传播
- class thermoelasticsim.md.schemes.NoseHooverChainPropagator(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
基类：
PropagatorNose–Hoover 链恒温器传播子（四阶 Suzuki–Yoshida 分解）
备注
链式热浴变量自洽演化，时间可逆对称分解，保证正则系综采样质量。
- 参考文献：
Martyna et al., J. Chem. Phys. 97, 2635 (1992); Yoshida, Phys. Lett. A 150, 262 (1990).
- __init__(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
初始化Nose-Hoover链恒温器
- 参数:
- 抛出:
ValueError -- 如果参数设置不合理
- get_conserved_energy(cell, potential=None)[源代码]
计算 NHC 扩展哈密顿量（守恒量）
\[H' = E_{\mathrm{kin}} + E_{\mathrm{pot}} + E_{\mathrm{thermo}}\]其中热浴能量 \(E_{\mathrm{thermo}}\) 为：
\[E_{\mathrm{thermo}} = \sum_{j=0}^{M-1} \frac{p_{\zeta_j}^2}{2 Q_j} + N_f k_B T_0\, \zeta_0 + k_B T_0 \sum_{j=1}^{M-1} \zeta_j\]
- propagate(cell, dt, **kwargs)[源代码]
执行 Nose–Hoover 链传播（四阶 Suzuki–Yoshida）
实现完整的 NHC 算符：
\[\exp\big(iL_{\mathrm{NHC}}\,\Delta t\big)\]采用三步四阶 Suzuki–Yoshida 分解保证高精度和长期稳定性。
算法流程：
遍历每个 SY 系数 \(w_k\)
对每个系数执行 \(t_{loop}\) 次内循环
调用
_nhc_integration_loop(w_k cdot Delta t / t_{loop})
- 参数:
- 返回类型:
None
备注
与 NVE 的对称包裹关系：
\[\exp\big(iL_{\mathrm{NHC}}\tfrac{\Delta t}{2}\big)\,\exp\big(iL_{\mathrm{NVE}}\,\Delta t\big)\,\exp\big(iL_{\mathrm{NHC}}\tfrac{\Delta t}{2}\big)\]
- class thermoelasticsim.md.schemes.PositionPropagator[源代码]
基类：
Propagator位置传播子：\(\exp(iL_r\,\Delta t)\)
实现位置的时间演化：\(\mathbf{r}(t+\Delta t)=\mathbf{r}(t)+\mathbf{v}(t)\,\Delta t\)。
备注
处理周期性边界；不改变速度与力；与势能无关。
- class thermoelasticsim.md.schemes.VelocityPropagator[源代码]
基类：
Propagator速度传播子：\(\exp(iL_v\,\Delta t)\)
实现速度的时间演化：\(\mathbf{v}(t+\Delta t)=\mathbf{v}(t)+\mathbf{F}(t)\,\Delta t/m\)。
备注
假设
atom.force已计算；不在此处重新计算力。
- class thermoelasticsim.md.schemes.NVEScheme[源代码]
-
微正则系综 (NVE) 的 Velocity-Verlet 积分方案
对称算符分离形式：
\[\exp(iL\,\Delta t) \approx \exp(iL_v\,\tfrac{\Delta t}{2})\,\exp(iL_r\,\Delta t)\,\exp(iL_v\,\tfrac{\Delta t}{2})\]其中
- \(iL_v\)
速度传播算符（力对动量的作用）
- \(iL_r\)
位置传播算符（动量对位置的作用）
备注
基本步骤（v-半步，r-整步，再 v-半步）：
\(\mathbf{v}(t+\tfrac{\Delta t}{2}) = \mathbf{v}(t) + \mathbf{F}(t)\,\Delta t/(2m)\)
\(\mathbf{r}(t+\Delta t) = \mathbf{r}(t) + \mathbf{v}(t+\tfrac{\Delta t}{2})\,\Delta t\)
计算 \(\mathbf{F}(t+\Delta t)\)
\(\mathbf{v}(t+\Delta t) = \mathbf{v}(t+\tfrac{\Delta t}{2}) + \mathbf{F}(t+\Delta t)\,\Delta t/(2m)\)
- step(cell, potential, dt)[源代码]
执行一步Velocity-Verlet积分
- 参数:
- 返回类型:
None
备注
该方法严格按照Velocity-Verlet算法执行：
第一次半步速度更新：使用当前位置的力
整步位置更新：使用更新后的速度
力重新计算：在新位置处计算力
第二次半步速度更新：使用新位置的力
这确保了算法的对称性和时间可逆性。
- 抛出:
ValueError -- 如果dt <= 0
RuntimeError -- 如果势函数计算失败
- 参数:
- 返回类型:
None
- class thermoelasticsim.md.schemes.BerendsenNVTScheme(target_temperature, tau=100.0)[源代码]
-
Berendsen NVT 积分方案（NVE + 速度缩放恒温器）
备注
两阶段流程：
NVE：标准 Velocity-Verlet 积分
NVT：应用 Berendsen 恒温器进行速度缩放
- __init__(target_temperature, tau=100.0)[源代码]
初始化Berendsen NVT积分方案
- 参数:
- 抛出:
ValueError -- 如果target_temperature或tau为非正数
- step(cell, potential, dt)[源代码]
执行一步Berendsen NVT积分
- 参数:
- 返回类型:
None
备注
该方法执行两阶段积分：
NVE阶段：标准Velocity-Verlet积分 - 速度半步更新 - 位置整步更新 - 力重新计算 - 速度再次半步更新
NVT阶段：Berendsen温度调节 - 计算当前温度 - 计算速度缩放因子 - 调节所有原子速度
- 抛出:
ValueError -- 如果dt <= 0
RuntimeError -- 如果积分过程失败
- 参数:
- 返回类型:
None
- class thermoelasticsim.md.schemes.AndersenNVTScheme(target_temperature, collision_frequency=0.01)[源代码]
-
Andersen NVT 积分方案（NVE + 随机碰撞恒温器）
备注
两阶段流程：
NVE：标准 Velocity-Verlet 积分
NVT：按碰撞概率重新采样 Maxwell 速度
- __init__(target_temperature, collision_frequency=0.01)[源代码]
初始化Andersen NVT积分方案
- 参数:
- 抛出:
ValueError -- 如果target_temperature或collision_frequency为非正数
- step(cell, potential, dt)[源代码]
执行一步Andersen NVT积分
- 参数:
- 返回类型:
None
备注
该方法执行两阶段积分：
NVE阶段：标准Velocity-Verlet积分 - 速度半步更新 - 位置整步更新 - 力重新计算 - 速度再次半步更新
NVT阶段：Andersen随机碰撞 - 为每个原子计算碰撞概率 - 随机选择碰撞原子 - 重新采样碰撞原子的Maxwell速度
- 抛出:
ValueError -- 如果dt <= 0
RuntimeError -- 如果积分过程失败
- 参数:
- 返回类型:
None
- class thermoelasticsim.md.schemes.NoseHooverNVTScheme(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
-
Nose–Hoover NVT 积分方案（算符分离）
基于对称分解：
\[\exp(iL\,\Delta t) \approx \exp(iL_{NH}\,\tfrac{\Delta t}{2})\,\exp(iL_v\,\tfrac{\Delta t}{2})\,\exp(iL_r\,\Delta t)\,\exp(iL_v\,\tfrac{\Delta t}{2})\,\exp(iL_{NH}\,\tfrac{\Delta t}{2})\]- __init__(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
初始化Nose-Hoover NVT积分方案
- 参数:
- 抛出:
ValueError -- 如果target_temperature为非正数或tdamp为非正数
- step(cell, potential, dt)[源代码]
执行一步Nose-Hoover NVT积分
- 参数:
- 返回类型:
None
备注
该方法执行严格的算符分离积分：
Nose-Hoover算符半步：处理热浴变量和速度摩擦
标准Velocity-Verlet积分： - 速度半步更新 - 位置整步更新 - 力重新计算 - 速度再次半步更新
Nose-Hoover算符再次半步：完成热浴变量演化
这种分解保证算法的时间可逆性和相空间体积保持。
- 抛出:
ValueError -- 如果dt == 0
RuntimeError -- 如果积分过程失败
- 参数:
- 返回类型:
None
- class thermoelasticsim.md.schemes.LangevinNVTScheme(target_temperature, friction=1.0)[源代码]
-
Langevin NVT 积分方案（Velocity-Verlet + BBK 速度修正）
- __init__(target_temperature, friction=1.0)[源代码]
初始化Langevin NVT积分方案
- 参数:
- 抛出:
ValueError -- 如果target_temperature或friction为非正数
- step(cell, potential, dt)[源代码]
执行一步Langevin NVT积分（BBK算法）
- 参数:
- 返回类型:
None
备注
该方法执行完整的BBK积分算法：
NVE阶段：标准Velocity-Verlet积分 - 半步速度更新（仅保守力） - 整步位置更新 - 力重新计算 - 半步速度更新（仅保守力）
NVT阶段：Langevin恒温器修正 - 计算BBK参数：c1, σ - 应用随机-摩擦速度修正 - 满足涨落-耗散定理
与其他恒温器的区别： - Berendsen: 速度缩放，不产生正确涨落 - Andersen: 随机重置，破坏动力学连续性更强 - Nose-Hoover: 确定性，但可能有长周期振荡 - Langevin: 连续随机修正，涨落-耗散平衡
- 抛出:
ValueError -- 如果dt <= 0
RuntimeError -- 如果积分过程失败
- 参数:
- 返回类型:
None
- get_statistics()[源代码]
获取积分与恒温器统计信息
- 返回:
统计信息字典：
- step_count
已执行的积分步数
- total_time
总积分时间 (fs)
- average_dt
平均时间步长 (fs)
- target_temperature
目标温度 (K)
- friction
摩擦系数 (ps⁻¹)
- damping_time
阻尼时间 (ps)
- thermostat_stats
恒温器详细统计
- 返回类型:
- set_friction(new_friction)[源代码]
动态调整摩擦系数
- 参数:
new_friction (float) -- 新的摩擦系数 (ps⁻¹)
- 抛出:
ValueError -- 如果new_friction为非正数
- 返回类型:
None
备注
这允许在模拟过程中调整恒温器的耦合强度。 典型用法： - 平衡阶段：使用大摩擦系数快速达到目标温度 - 生产阶段：使用小摩擦系数保持动力学真实性
- thermoelasticsim.md.schemes.create_nose_hoover_nvt_scheme(target_temperature, tdamp=100.0, tchain=3, tloop=1)[源代码]
创建标准的Nose-Hoover NVT积分方案
- 参数:
- 返回:
配置好的Nose-Hoover NVT积分方案实例
- 返回类型:
示例
>>> scheme = create_nose_hoover_nvt_scheme(300.0, tdamp=50.0) >>> for step in range(10000): ... scheme.step(cell, potential, dt=0.5) >>> stats = scheme.get_statistics() >>> nhc_vars = scheme.get_thermostat_variables() >>> print(f"热浴动量: {nhc_vars['p_zeta']}")
- thermoelasticsim.md.schemes.create_andersen_nvt_scheme(target_temperature, collision_frequency=0.01)[源代码]
创建标准的Andersen NVT积分方案
- 参数:
- 返回:
配置好的Andersen NVT积分方案实例
- 返回类型:
示例
>>> scheme = create_andersen_nvt_scheme(300.0, collision_frequency=0.02) >>> for step in range(10000): ... scheme.step(cell, potential, dt=0.5) >>> stats = scheme.get_statistics() >>> collision_stats = scheme.get_collision_statistics() >>> print(f"总碰撞次数: {collision_stats['total_collisions']}")
- thermoelasticsim.md.schemes.create_berendsen_nvt_scheme(target_temperature, tau=100.0)[源代码]
创建标准的Berendsen NVT积分方案
- 参数:
- 返回:
配置好的Berendsen NVT积分方案实例
- 返回类型:
示例
>>> scheme = create_berendsen_nvt_scheme(300.0, tau=50.0) >>> for step in range(10000): ... scheme.step(cell, potential, dt=0.5) >>> stats = scheme.get_statistics() >>> print(f"平均温度调节: {stats['thermostat_stats']['average_scaling']:.3f}")
- thermoelasticsim.md.schemes.create_langevin_nvt_scheme(target_temperature, friction=1.0)[源代码]
创建标准的Langevin NVT积分方案
基于BBK积分算法，提供基于物理模型的随机恒温。 结合摩擦阻力和随机力，通过涨落-耗散定理确保严格的正则系综采样。
- 参数:
- 返回:
配置好的Langevin NVT积分方案实例
- 返回类型:
示例
平衡阶段使用强耦合快速达到目标温度:
>>> # 强耦合快速平衡 >>> scheme = create_langevin_nvt_scheme(300.0, friction=5.0) # 强耦合 >>> for step in range(5000): # 快速平衡 ... scheme.step(cell, potential, dt=1.0)
生产阶段使用弱耦合保持动力学真实性:
>>> # 弱耦合生产模拟 >>> scheme.set_friction(0.5) # 切换到弱耦合 >>> for step in range(100000): # 长时间生产 ... scheme.step(cell, potential, dt=1.0) >>> stats = scheme.get_statistics() >>> print(f"平均温度: {stats['thermostat_stats']['average_temperature']:.1f} K") >>> print(f"能量平衡: {stats['thermostat_stats']['energy_balance']:.3f}")
备注
Langevin恒温器的优势：
严格正则系综: 基于涨落-耗散定理，确保正确的统计力学采样
数值稳定: 无Nose-Hoover的遍历性问题和长周期振荡
物理清晰: 明确的摩擦和随机力物理模型
参数灵活: 可动态调节耦合强度，适应不同模拟阶段
广泛适用: 特别适合生物分子和隐式溶剂系统
与其他恒温器的比较：
vs Berendsen: 产生正确涨落，不只是平均温度控制
vs Andersen: 连续随机修正，不是突然重置
vs Nose-Hoover: 随机而非确定性，避免长周期问题
参数选择建议：
快速平衡: γ = 5-10 ps⁻¹，快速达到热平衡
动力学研究: γ = 0.1-1 ps⁻¹，最小化对真实动力学的扰动
构象采样: γ = 1-5 ps⁻¹，平衡采样效率和动力学保真度
布朗动力学: γ > 10 ps⁻¹，主要由扩散控制
- thermoelasticsim.md.schemes.create_nve_scheme()[源代码]
创建标准的NVE积分方案
- 返回:
配置好的NVE积分方案实例
- 返回类型:
示例
>>> scheme = create_nve_scheme() >>> for step in range(1000): ... scheme.step(cell, potential, dt=0.5)
- class thermoelasticsim.md.schemes.MTKNPTScheme(target_temperature, target_pressure, tdamp, pdamp, tchain=3, pchain=3, tloop=1, ploop=1)[源代码]
-
MTK–NPT 积分方案（Nose–Hoover 链 + MTK 恒压器）
备注
对称分解示意：
\[\exp(iL\,\Delta t) \approx \exp(iL_{baro}\tfrac{\Delta t}{2})\,\exp(iL_{thermo}\tfrac{\Delta t}{2})\,\exp(iL_{p\_cell}\tfrac{\Delta t}{2})\,\exp(iL_p\tfrac{\Delta t}{2})\,\exp(iL_q\,\Delta t)\,\exp(iL_p\tfrac{\Delta t}{2})\,\exp(iL_{p\_cell}\tfrac{\Delta t}{2})\,\exp(iL_{thermo}\tfrac{\Delta t}{2})\,\exp(iL_{baro}\tfrac{\Delta t}{2})\]功能要点：
注意事项: - 比NVE和NVT积分需要更多计算时间 - 压力控制参数需要仔细调节 - 适中的时间步长确保数值稳定性
- 参数:
target_temperature (float) -- 目标温度 (K)
target_pressure (float) -- 目标压力 (GPa) - 自动转换为内部单位
tdamp (float) -- 温度阻尼时间 (fs)，典型值50-100*dt
pdamp (float) -- 压力阻尼时间 (fs)，典型值100-1000*dt
tchain (int, optional) -- 恒温器链长度，默认为3
pchain (int, optional) -- 恒压器链长度，默认为3
tloop (int, optional) -- 恒温器积分循环次数，默认为1
ploop (int, optional) -- 恒压器积分循环次数，默认为1
示例
>>> # 300K, 0.1GPa的NPT模拟 >>> scheme = MTKNPTScheme( ... target_temperature=300.0, ... target_pressure=0.1, # GPa ... tdamp=50.0, # fs ... pdamp=500.0 # fs ... ) >>> >>> # 运行10ps的NPT弛豫 >>> dt = 1.0 # fs >>> for step in range(10000): ... scheme.step(cell, potential, dt) ... ... if step % 1000 == 0: ... stats = scheme.get_statistics() ... print(f"T={stats['temperature']:.1f}K, " ... f"P={stats['pressure']:.3f}GPa")
- __init__(target_temperature, target_pressure, tdamp, pdamp, tchain=3, pchain=3, tloop=1, ploop=1)[源代码]
- thermoelasticsim.md.schemes.create_mtk_npt_scheme(target_temperature, target_pressure, tdamp, pdamp, tchain=3, pchain=3)[源代码]
创建MTK-NPT积分方案
- 参数:
- 返回:
配置好的MTK-NPT积分方案
- 返回类型:
示例
>>> # 300K, 大气压的NPT系综 >>> scheme = create_mtk_npt_scheme( ... target_temperature=300.0, ... target_pressure=0.0001, # ~大气压 ... tdamp=50.0, ... pdamp=500.0 ... ) >>> >>> # 高压高温条件 >>> scheme_hp = create_mtk_npt_scheme( ... target_temperature=1000.0, ... target_pressure=10.0, # 10 GPa ... tdamp=100.0, ... pdamp=1000.0 ... )
恒温器
- class thermoelasticsim.md.thermostats.CppInterface(lib_name)[源代码]
基类：
object用于调用 C++ 实现的函数的接口类
@class CppInterface @brief 用于调用 C++ 实现的函数的接口类
- 参数:
lib_name (str) -- 库的名称，用于确定可用的函数集合。
- __init__(lib_name)[源代码]
- calculate_eam_al1_energy(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Al1势的总能量。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
- 返回:
系统的总能量，单位为eV
- 返回类型:
- calculate_eam_al1_forces(num_atoms, positions, lattice_vectors, forces)[源代码]
计算EAM Al1势的原子力。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
forces (numpy.ndarray) -- 输出的力数组，形状为(num_atoms, 3)，将被更新
- 返回类型:
None
- calculate_eam_al1_virial(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Al1的维里张量（未除以体积）。返回形状(3,3)。
- calculate_eam_cu1_energy(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Cu1势的总能量。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
- 返回:
系统的总能量，单位为eV
- 返回类型:
- calculate_eam_cu1_forces(num_atoms, positions, lattice_vectors, forces)[源代码]
计算EAM Cu1势的原子力。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
forces (numpy.ndarray) -- 输出的力数组，形状为(num_atoms, 3)，将被更新
- 返回类型:
None
- calculate_eam_cu1_virial(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Cu1的维里张量（未除以体积）。返回形状(3,3)。
- calculate_lj_energy(num_atoms, positions, epsilon, sigma, cutoff, box_lengths, neighbor_pairs, num_pairs)[源代码]
调用 C++ 接口计算能量。
- 参数:
num_atoms (int) -- 原子数。
positions (numpy.ndarray) -- 原子位置数组，形状为 (num_atoms, 3)。
epsilon (float) -- Lennard-Jones 势参数 ε，单位 eV。
sigma (float) -- Lennard-Jones 势参数 σ，单位 Å。
cutoff (float) -- 截断距离，单位 Å。
box_lengths (numpy.ndarray) -- 盒子长度数组，形状为 (3,)。
neighbor_pairs (numpy.ndarray) -- 邻居对数组，形状为 (2*num_pairs,)。
num_pairs (int) -- 邻居对的数量。
- 返回:
总势能，单位 eV。
- 返回类型:
- calculate_lj_forces(num_atoms, positions, forces, epsilon, sigma, cutoff, box_lengths, neighbor_pairs, num_pairs)[源代码]
调用 C++ 接口计算作用力。
- 参数:
num_atoms (int) -- 原子数。
positions (numpy.ndarray) -- 原子位置数组，形状为 (num_atoms, 3)。
forces (numpy.ndarray) -- 力数组，形状为 (num_atoms, 3)，将被更新。
epsilon (float) -- Lennard-Jones 势参数 ε，单位 eV。
sigma (float) -- Lennard-Jones 势参数 σ，单位 Å。
cutoff (float) -- 截断距离，单位 Å。
box_lengths (numpy.ndarray) -- 盒子长度数组，形状为 (3,)。
neighbor_pairs (numpy.ndarray) -- 邻居对数组，形状为 (2*num_pairs,)。
num_pairs (int) -- 邻居对的数量。
- 返回类型:
None
- calculate_tersoff_c1988_energy(num_atoms, positions, lattice_vectors, *, shift_flag=False, delta=0.0)[源代码]
计算Tersoff C(1988)势能量。
- calculate_tersoff_c1988_forces(num_atoms, positions, lattice_vectors, forces, *, shift_flag=False, delta=0.0)[源代码]
计算Tersoff C(1988)势力。
- calculate_tersoff_c1988_virial(num_atoms, positions, lattice_vectors, *, shift_flag=False, delta=0.0)[源代码]
计算Tersoff C(1988)势维里张量。
- calculate_tersoff_energy(num_atoms, positions, lattice_vectors, *, A, B, lambda1, lambda2, lambda3, beta, n, c, d, h, R, D, m=3, shift_flag=False, delta=0.0)[源代码]
计算Tersoff势能量。
- calculate_tersoff_forces(num_atoms, positions, lattice_vectors, forces, *, A, B, lambda1, lambda2, lambda3, beta, n, c, d, h, R, D, m=3, shift_flag=False, delta=0.0)[源代码]
计算Tersoff势力。
- calculate_tersoff_virial(num_atoms, positions, lattice_vectors, *, A, B, lambda1, lambda2, lambda3, beta, n, c, d, h, R, D, m=3, shift_flag=False, delta=0.0)[源代码]
计算Tersoff势维里张量。
- compute_stress(num_atoms, positions, velocities, forces, masses, volume, box_lengths, stress_tensor)[源代码]
计算应力张量。
本函数允许输入的 positions, velocities, forces 既可以是 (num_atoms, 3) 也可以是 (3*num_atoms,) 的形状。 同理，对于 stress_tensor，既可以是 (3,3) 也可以是 (9,)。
- 参数:
num_atoms (int)
positions (np.ndarray) -- 原子位置数组，可为 (num_atoms, 3) 或 (3*num_atoms, )
velocities (np.ndarray) -- 原子速度数组，可为 (num_atoms, 3) 或 (3*num_atoms, )
forces (np.ndarray) -- 原子力数组，可为 (num_atoms, 3) 或 (3*num_atoms, )
masses (np.ndarray) -- 原子质量数组 (num_atoms,)
volume (float) -- 晶胞体积
box_lengths (np.ndarray) -- 晶胞长度数组 (3,)
stress_tensor (np.ndarray) -- 输出应力张量，可为 (3,3) 或 (9,)
- 返回类型:
None
- nose_hoover(dt, num_atoms, masses, velocities, forces, xi_array, Q, target_temperature)[源代码]
实现 Nose-Hoover 恒温器算法。!!!没有去除质心运动
- 参数:
dt (float) -- 时间步长。
num_atoms (int) -- 原子数。
masses (numpy.ndarray) -- 原子质量数组。
velocities (numpy.ndarray) -- 原子速度数组。
forces (numpy.ndarray) -- 原子受力数组。
xi_array (numpy.ndarray) -- Nose-Hoover 热浴变量数组。
Q (float) -- 热浴质量参数。
target_temperature (float) -- 目标温度。
- 返回:
更新后的 Nose-Hoover 热浴变量。
- 返回类型:
- nose_hoover_chain(dt, num_atoms, masses, velocities, forces, xi_chain, Q, chain_length, target_temperature)[源代码]
实现 Nose-Hoover 链恒温器算法。
- 参数:
dt (float) -- 时间步长。
num_atoms (int) -- 原子数。
masses (numpy.ndarray) -- 原子质量数组。
velocities (numpy.ndarray) -- 原子速度数组。
forces (numpy.ndarray) -- 原子受力数组。
xi_chain (numpy.ndarray) -- Nose-Hoover 链的热浴变量数组。
Q (numpy.ndarray) -- 热浴质量参数数组。
chain_length (int) -- Nose-Hoover 链的长度。
target_temperature (float) -- 目标温度。
- parrinello_rahman_hoover(dt, num_atoms, masses, velocities, forces, lattice_vectors, xi, Q, total_stress, target_pressure, W)[源代码]
执行Parrinello-Rahman-Hoover恒压器积分步骤
- class thermoelasticsim.md.thermostats.Thermostat(target_temperature)[源代码]
基类：
object恒温器基类，定义恒温器的接口和通用属性
- 参数:
target_temperature (float) -- 目标温度，单位K
- class thermoelasticsim.md.thermostats.BerendsenThermostat(target_temperature, tau)[源代码]
基类：
ThermostatBerendsen 恒温器
通过速度缩放实现温度控制，使系统快速达到目标温度。 该方法不产生严格的正则系综，但具有良好的数值稳定性。
- 参数:
- apply(cell, dt, potential)[源代码]
应用 Berendsen 恒温器以控制系统温度。
该恒温器通过以下缩放因子 \(\lambda\) 调整原子速度：
\[\lambda = \sqrt{1 + \frac{\Delta t}{\tau} \left(\frac{T_{\text{target}}}{T_{\text{current}}} - 1\right)}\]其中
- \(\Delta t\)
时间步长 dt
- \(T_{\text{target}}\)
目标温度 target_temperature
- \(T_{\text{current}}\)
当前系统温度
- class thermoelasticsim.md.thermostats.AndersenThermostat(target_temperature, collision_frequency)[源代码]
基类：
ThermostatAndersen 恒温器
通过随机碰撞来控制温度，实现符合正则系综的采样。 该方法适合平衡态采样，但由于速度随机化，动力学轨迹会被打断。
- 参数:
- __init__(target_temperature, collision_frequency)[源代码]
- apply(cell, dt, potential)[源代码]
应用 Andersen 恒温器以控制系统温度。
在每个时间步，针对每个原子，发生碰撞的概率为 \(p\)：
\[p = \nu \times \Delta t\]若发生碰撞，则根据目标温度 \(T_{\text{target}}\) 的 Maxwell-Boltzmann 分布重新分配原子速度：
\[v_i \sim \mathcal{N}(0, \sigma)\]其中 \(\sigma\) 为速度分布的标准差：
\[\sigma = \sqrt{\frac{k_B\, T_{\text{target}}}{m}}\]
- class thermoelasticsim.md.thermostats.NoseHooverThermostat(target_temperature, time_constant, Q=None)[源代码]
基类：
ThermostatNosé–Hoover 恒温器（已弃用）
自 4.0 版本弃用: 该旧架构实现已弃用；请使用新架构中的
thermoelasticsim.md.propagators.NoseHooverChainPropagator。备注
历史说明：本类为早期 Python 直接实现版本，用对称分解包裹 NVE。 新版实现提供更稳健的 NHC 链与严格的算符分离。
- __init__(target_temperature, time_constant, Q=None)[源代码]
- class thermoelasticsim.md.thermostats.NoseHooverChainThermostat(target_temperature, time_constant, chain_length=3, Q=None)[源代码]
基类：
ThermostatNosé–Hoover 链恒温器（已弃用）
自 4.0 版本弃用: 该旧架构实现已弃用；请使用新架构中的
thermoelasticsim.md.propagators.NoseHooverChainPropagator。备注
历史说明：本类依赖 C++ 接口；新版以纯 Python 传播子架构重写， 文档与接口更一致，建议迁移。
- __init__(target_temperature, time_constant, chain_length=3, Q=None)[源代码]
恒压器
- class thermoelasticsim.md.barostats.StressCalculator[源代码]
基类：
object应力张量计算器
总应力由两部分组成（体积 \(V\)）：
- 动能项（动量通量）
- \[\sigma^{\text{kin}}_{\alpha\beta} = -\,\frac{1}{V} \sum_i m_i\, v_{i\alpha} v_{i\beta}\]
- 维里项（力-位矢项）
- \[\sigma^{\text{vir}}_{\alpha\beta} = -\,\frac{1}{V} \sum_{i<j} r_{ij,\alpha}\, F_{ij,\beta}\]
备注
有限差分/晶格导数法（如 \(-\partial U/\partial \varepsilon\) 或 \(-V^{-1}\,\partial U/\partial \mathbf{h}\, \mathbf{h}^T\)）是另一种等价 的应力计算方式，用于数值校验；它并非额外“第三项”，不与维里项相加。
本实现采用 \(\sigma = \sigma^{\text{kin}} + \sigma^{\text{vir}}\)。
对 EAM 势，使用经特殊处理的多体维里解析形式；如解析不可用，可用有限差分做校验。
正负号及指标约定与项目其它模块保持一致。
- __init__()[源代码]
- calculate_finite_difference_stress(cell, potential, dr=1e-06)[源代码]
已移除：有限差分应力路径不再提供，避免误用。
- 返回类型:
- calculate_kinetic_stress(cell)[源代码]
计算动能应力张量
使用的约定：
\[\sigma^{\text{kin}}_{\alpha\beta} = -\,\frac{1}{V} \sum_i m_i\, v_{i\alpha} v_{i\beta}.\]说明：在静止构型或 \(T\to 0\) 极限，速度趋近于零，因此该项趋近于 0。
- 参数:
cell (Cell) -- 包含原子的晶胞对象
- 返回:
动能应力张量 (3, 3)，单位 eV/ų
- 返回类型:
- calculate_total_stress(cell, potential)[源代码]
计算总应力张量（动能项 + 维里项）
\[\sigma_{\alpha\beta} = -\,\frac{1}{V} \left( \sum_i m_i v_{i\alpha} v_{i\beta} + \sum_{i<j} r_{ij,\alpha} F_{ij,\beta} \right)\]
- calculate_virial_stress(cell, potential)[源代码]
计算维里应力张量（相互作用力贡献）
公式：
\[\sigma^{\text{vir}}_{\alpha\beta} = -\,\frac{1}{V} \sum_{i<j} r_{ij,\alpha}\, F_{ij,\beta}\]其中：
- \(\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j\)
从原子 \(j\) 指向原子 \(i\) 的位移向量
- \(\mathbf{F}_{ij}\)
原子 \(j\) 作用于 \(i\) 的力
- compute_stress(cell, potential)[源代码]
计算应力张量（使用总应力作为主要方法）
- get_all_stress_components(cell, potential)[源代码]
计算应力张量的所有分量
- class thermoelasticsim.md.barostats.CppInterface(lib_name)[源代码]
基类：
object用于调用 C++ 实现的函数的接口类
@class CppInterface @brief 用于调用 C++ 实现的函数的接口类
- 参数:
lib_name (str) -- 库的名称，用于确定可用的函数集合。
- __init__(lib_name)[源代码]
- calculate_eam_al1_energy(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Al1势的总能量。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
- 返回:
系统的总能量，单位为eV
- 返回类型:
- calculate_eam_al1_forces(num_atoms, positions, lattice_vectors, forces)[源代码]
计算EAM Al1势的原子力。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
forces (numpy.ndarray) -- 输出的力数组，形状为(num_atoms, 3)，将被更新
- 返回类型:
None
- calculate_eam_al1_virial(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Al1的维里张量（未除以体积）。返回形状(3,3)。
- calculate_eam_cu1_energy(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Cu1势的总能量。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
- 返回:
系统的总能量，单位为eV
- 返回类型:
- calculate_eam_cu1_forces(num_atoms, positions, lattice_vectors, forces)[源代码]
计算EAM Cu1势的原子力。
- 参数:
num_atoms (int) -- 原子数量
positions (numpy.ndarray) -- 原子位置数组，形状为(num_atoms, 3)
lattice_vectors (numpy.ndarray) -- 晶格向量数组，形状为(3, 3)或(9,)
forces (numpy.ndarray) -- 输出的力数组，形状为(num_atoms, 3)，将被更新
- 返回类型:
None
- calculate_eam_cu1_virial(num_atoms, positions, lattice_vectors)[源代码]
计算EAM Cu1的维里张量（未除以体积）。返回形状(3,3)。
- calculate_lj_energy(num_atoms, positions, epsilon, sigma, cutoff, box_lengths, neighbor_pairs, num_pairs)[源代码]
调用 C++ 接口计算能量。
- 参数:
num_atoms (int) -- 原子数。
positions (numpy.ndarray) -- 原子位置数组，形状为 (num_atoms, 3)。
epsilon (float) -- Lennard-Jones 势参数 ε，单位 eV。
sigma (float) -- Lennard-Jones 势参数 σ，单位 Å。
cutoff (float) -- 截断距离，单位 Å。
box_lengths (numpy.ndarray) -- 盒子长度数组，形状为 (3,)。
neighbor_pairs (numpy.ndarray) -- 邻居对数组，形状为 (2*num_pairs,)。
num_pairs (int) -- 邻居对的数量。
- 返回:
总势能，单位 eV。
- 返回类型:
- calculate_lj_forces(num_atoms, positions, forces, epsilon, sigma, cutoff, box_lengths, neighbor_pairs, num_pairs)[源代码]
调用 C++ 接口计算作用力。
- 参数:
num_atoms (int) -- 原子数。
positions (numpy.ndarray) -- 原子位置数组，形状为 (num_atoms, 3)。
forces (numpy.ndarray) -- 力数组，形状为 (num_atoms, 3)，将被更新。
epsilon (float) -- Lennard-Jones 势参数 ε，单位 eV。
sigma (float) -- Lennard-Jones 势参数 σ，单位 Å。
cutoff (float) -- 截断距离，单位 Å。
box_lengths (numpy.ndarray) -- 盒子长度数组，形状为 (3,)。
neighbor_pairs (numpy.ndarray) -- 邻居对数组，形状为 (2*num_pairs,)。
num_pairs (int) -- 邻居对的数量。
- 返回类型:
None
- calculate_tersoff_c1988_energy(num_atoms, positions, lattice_vectors, *, shift_flag=False, delta=0.0)[源代码]
计算Tersoff C(1988)势能量。
- calculate_tersoff_c1988_forces(num_atoms, positions, lattice_vectors, forces, *, shift_flag=False, delta=0.0)[源代码]
计算Tersoff C(1988)势力。
- calculate_tersoff_c1988_virial(num_atoms, positions, lattice_vectors, *, shift_flag=False, delta=0.0)[源代码]
计算Tersoff C(1988)势维里张量。
- calculate_tersoff_energy(num_atoms, positions, lattice_vectors, *, A, B, lambda1, lambda2, lambda3, beta, n, c, d, h, R, D, m=3, shift_flag=False, delta=0.0)[源代码]
计算Tersoff势能量。
- calculate_tersoff_forces(num_atoms, positions, lattice_vectors, forces, *, A, B, lambda1, lambda2, lambda3, beta, n, c, d, h, R, D, m=3, shift_flag=False, delta=0.0)[源代码]
计算Tersoff势力。
- calculate_tersoff_virial(num_atoms, positions, lattice_vectors, *, A, B, lambda1, lambda2, lambda3, beta, n, c, d, h, R, D, m=3, shift_flag=False, delta=0.0)[源代码]
计算Tersoff势维里张量。
- compute_stress(num_atoms, positions, velocities, forces, masses, volume, box_lengths, stress_tensor)[源代码]
计算应力张量。
本函数允许输入的 positions, velocities, forces 既可以是 (num_atoms, 3) 也可以是 (3*num_atoms,) 的形状。 同理，对于 stress_tensor，既可以是 (3,3) 也可以是 (9,)。
- 参数:
num_atoms (int)
positions (np.ndarray) -- 原子位置数组，可为 (num_atoms, 3) 或 (3*num_atoms, )
velocities (np.ndarray) -- 原子速度数组，可为 (num_atoms, 3) 或 (3*num_atoms, )
forces (np.ndarray) -- 原子力数组，可为 (num_atoms, 3) 或 (3*num_atoms, )
masses (np.ndarray) -- 原子质量数组 (num_atoms,)
volume (float) -- 晶胞体积
box_lengths (np.ndarray) -- 晶胞长度数组 (3,)
stress_tensor (np.ndarray) -- 输出应力张量，可为 (3,3) 或 (9,)
- 返回类型:
None
- nose_hoover(dt, num_atoms, masses, velocities, forces, xi_array, Q, target_temperature)[源代码]
实现 Nose-Hoover 恒温器算法。!!!没有去除质心运动
- 参数:
dt (float) -- 时间步长。
num_atoms (int) -- 原子数。
masses (numpy.ndarray) -- 原子质量数组。
velocities (numpy.ndarray) -- 原子速度数组。
forces (numpy.ndarray) -- 原子受力数组。
xi_array (numpy.ndarray) -- Nose-Hoover 热浴变量数组。
Q (float) -- 热浴质量参数。
target_temperature (float) -- 目标温度。
- 返回:
更新后的 Nose-Hoover 热浴变量。
- 返回类型:
- nose_hoover_chain(dt, num_atoms, masses, velocities, forces, xi_chain, Q, chain_length, target_temperature)[源代码]
实现 Nose-Hoover 链恒温器算法。
- 参数:
dt (float) -- 时间步长。
num_atoms (int) -- 原子数。
masses (numpy.ndarray) -- 原子质量数组。
velocities (numpy.ndarray) -- 原子速度数组。
forces (numpy.ndarray) -- 原子受力数组。
xi_chain (numpy.ndarray) -- Nose-Hoover 链的热浴变量数组。
Q (numpy.ndarray) -- 热浴质量参数数组。
chain_length (int) -- Nose-Hoover 链的长度。
target_temperature (float) -- 目标温度。
- parrinello_rahman_hoover(dt, num_atoms, masses, velocities, forces, lattice_vectors, xi, Q, total_stress, target_pressure, W)[源代码]
执行Parrinello-Rahman-Hoover恒压器积分步骤
- class thermoelasticsim.md.barostats.Barostat(target_pressure)[源代码]
基类：
object恒压器基类，定义恒压器的接口
- 参数:
target_pressure (ndarray)
- apply(cell, potential, dt)[源代码]
应用恒压器，更新晶胞参数
- calculate_internal_pressure(stress_tensor)[源代码]
根据应力张量计算内部压力
- 参数:
stress_tensor (np.ndarray) -- 应力张量展平后的数组，形状 (9,)
- 返回:
内部压力（eV/ų）
- 返回类型:
- class thermoelasticsim.md.barostats.ParrinelloRahmanHooverBarostat(target_pressure, time_constant, compressibility=4.57e-05, W=None, Q=None, stress_calculator=None)[源代码]
基类：
BarostatParrinello–Rahman–Hoover 恒压器（已弃用）
自 4.0 版本弃用: 该旧架构恒压器已弃用；请使用新架构中的
thermoelasticsim.md.propagators.MTKBarostatPropagator。- 参数:
target_pressure (np.ndarray) -- 目标压力张量 (3×3，单位 eV/ų)
time_constant (float) -- 压力耦合时间常数 (fs)
compressibility (float, optional) -- 等温压缩系数（无量纲或与内部单位一致），默认 4.57e-5
W (float, optional) -- 晶胞“质量”参数（内部单位）
Q (np.ndarray, optional) -- 热浴质量参数数组，形状 (9,)
stress_calculator (StressCalculator) -- 应力张量计算器实例
- __init__(target_pressure, time_constant, compressibility=4.57e-05, W=None, Q=None, stress_calculator=None)[源代码]
- apply(cell, potential, dt)[源代码]
应用 Parrinello-Rahman-Hoover 恒压器，更新晶胞参数和原子速度
- class thermoelasticsim.md.barostats.BerendsenBarostat(target_pressure, tau_p, compressibility=4.57e-05)[源代码]
基类：
BarostatBerendsen 恒压器（已弃用）
自 4.0 版本弃用: 该旧架构恒压器已弃用；请使用新架构中的
thermoelasticsim.md.propagators.MTKBarostatPropagator。- 参数:
- __init__(target_pressure, tau_p, compressibility=4.57e-05)[源代码]
- class thermoelasticsim.md.barostats.AndersenBarostat(target_pressure, mass, temperature)[源代码]
基类：
BarostatAndersen 恒压器（已弃用）
自 4.0 版本弃用: 该旧架构恒压器已弃用；请使用新架构中的
thermoelasticsim.md.propagators.MTKBarostatPropagator。- 参数:
- __init__(target_pressure, mass, temperature)[源代码]
MD模拟器
- class thermoelasticsim.md.md_simulator.AndersenThermostat(target_temperature, collision_frequency)[源代码]
基类：
ThermostatAndersen 恒温器
通过随机碰撞来控制温度，实现符合正则系综的采样。 该方法适合平衡态采样，但由于速度随机化，动力学轨迹会被打断。
- 参数:
- class thermoelasticsim.md.md_simulator.BerendsenThermostat(target_temperature, tau)[源代码]
基类：
ThermostatBerendsen 恒温器
通过速度缩放实现温度控制，使系统快速达到目标温度。 该方法不产生严格的正则系综，但具有良好的数值稳定性。
- 参数:
- class thermoelasticsim.md.md_simulator.NoseHooverChainThermostat(target_temperature, time_constant, chain_length=3, Q=None)[源代码]
基类：
ThermostatNosé–Hoover 链恒温器（已弃用）
自 4.0 版本弃用: 该旧架构实现已弃用；请使用新架构中的
thermoelasticsim.md.propagators.NoseHooverChainPropagator。备注
历史说明：本类依赖 C++ 接口；新版以纯 Python 传播子架构重写， 文档与接口更一致，建议迁移。
- class thermoelasticsim.md.md_simulator.NoseHooverThermostat(target_temperature, time_constant, Q=None)[源代码]
基类：
ThermostatNosé–Hoover 恒温器（已弃用）
自 4.0 版本弃用: 该旧架构实现已弃用；请使用新架构中的
thermoelasticsim.md.propagators.NoseHooverChainPropagator。备注
历史说明：本类为早期 Python 直接实现版本，用对称分解包裹 NVE。 新版实现提供更稳健的 NHC 链与严格的算符分离。