弹性常数计算模块 (elastic) 

备注

数学表述：

\[\min_{r,h} E(r,h) \quad \text{s.t.} \quad \sigma = 0\]

符号说明： - \(r\) 原子位置 - \(h\) 晶胞参数 - \(E\) 总能量 - \(\sigma\) 应力张量

示例

>>> relaxer = StructureRelaxer()
>>> converged = relaxer.full_relax(cell, potential)
>>> if converged:
...     print("成功制备无应力基态")

internal_relax(cell, potential)[源代码]

执行内部弛豫：仅优化原子位置，保持晶胞形状固定

内部弛豫用于形变后的结构优化，在保持宏观形变的前提下，寻找原子的最优位置配置。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{r} E(r,h_{\text{fixed}}) \quad \text{s.t.} \quad h = \text{const}\]

符号说明： - \(r\) 原子位置（优化变量） - \(h_{\text{fixed}}\) 固定的晶胞参数 - \(E\) 总能量

示例

>>> relaxer = StructureRelaxer()
>>> # 先施加形变
>>> cell.apply_deformation(deformation_matrix)
>>> # 再内部弛豫
>>> converged = relaxer.internal_relax(cell, potential)

uniform_lattice_relax(cell, potential)[源代码]

等比例晶格弛豫：只优化晶格常数，保持原子相对位置不变

这种弛豫方式特别适合基态制备，既避免了原子位置的数值噪声，又保持了晶体结构的完美对称性，同时计算效率极高。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{s} E(r_{\text{scaled}}, s \cdot L) \quad \text{其中} \quad r_{\text{scaled}} = s \cdot r_{\text{frac}} \cdot L\]

符号说明： - \(s\) 等比例缩放因子（唯一优化变量） - \(r_{\text{frac}}\) 固定的分数坐标 - \(L\) 原始晶格矢量矩阵 - \(E\) 总能量

优势： 1. 自由度仅为 1，计算极快 2. 保持晶体结构完美对称性 3. 避免原子位置数值噪声 4. 适合各种晶系（不限于 FCC）

示例

>>> relaxer = StructureRelaxer()
>>> converged = relaxer.uniform_lattice_relax(cell, potential)
>>> if converged:
...     print("成功优化晶格常数，保持结构对称性")

class thermoelasticsim.elastic.deformation_method.zero_temp.ZeroTempDeformationCalculator(cell, potential, delta=0.005, num_steps=5, relaxer_params=None, supercell_dims=None, base_relax_mode='uniform')[源代码]

基类：object

零温显式形变法计算器。

负责从基态制备到弹性常数求解的完整流程管理。

参数:

cell (Cell) -- 待计算的晶胞对象。
potential (Potential) -- 势能函数对象。
delta (float, optional) -- 应变幅度，默认 0.005 (0.5%)。
num_steps (int, optional) -- 每个应变分量的步数，默认 5（教学模式）。
relaxer_params (dict, optional) -- 结构弛豫器参数。
supercell_dims (tuple[int, int, int], optional) -- 超胞维度 (nx, ny, nz)。
base_relax_mode ({"uniform", "full"}, optional) -- 基态制备模式，等比例缩放或变胞+位置。

__init__(cell, potential, delta=0.005, num_steps=5, relaxer_params=None, supercell_dims=None, base_relax_mode='uniform')[源代码]

初始化零温形变计算器

参数:

cell (Cell) -- 晶胞对象
potential (Potential) -- 势能函数对象
delta (float, optional) -- 应变幅度，建议范围0.001-0.01，默认0.005
num_steps (int, optional) -- 每个应变分量的步数，1为生产模式，>1为教学模式，默认5
relaxer_params (dict, optional) -- 传递给StructureRelaxer的参数
supercell_dims (tuple, optional) -- 超胞维度(nx, ny, nz)，用于正确显示等效单胞参数
base_relax_mode (str)

抛出:

ValueError -- 如果delta或num_steps参数不合理

calculate()[源代码]

执行完整的零温弹性常数计算。

返回:: 弹性常数矩阵 (GPa) 和拟合优度 R^2。
返回类型:: tuple[numpy.ndarray, float]

class thermoelasticsim.elastic.deformation_method.zero_temp.ElasticConstantsSolver[源代码]

基类：object

弹性常数求解器

从应力应变数据通过线性回归求解弹性常数矩阵。支持最小二乘法和岭回归两种方法，并提供拟合质量评估。

solve(strains, stresses, method='least_squares')[源代码]

求解弹性常数矩阵

参数:

strains (ndarray)
stresses (ndarray)
method (str)
alpha (float)

返回类型:

tuple[ndarray, float]

备注

数学原理：

根据胡克定律：\(\sigma = C \cdot \varepsilon\)

其中：

\(\sigma\): 应力向量 (N×6)
\(C\): 弹性常数矩阵 (6×6)
\(\varepsilon\): 应变向量 (N×6)

最小二乘求解：

\[\min_C ||\sigma - C \cdot \varepsilon||^2\]

解析解：

\[C = (\varepsilon^T \varepsilon)^{-1} \varepsilon^T \sigma\]

示例

>>> solver = ElasticConstantsSolver()
>>> C_matrix, r2_score = solver.solve(strains, stresses)
>>> print(f"弹性常数矩阵:\\n{C_matrix}")
>>> print(f"拟合优度: {r2_score:.6f}")

__init__()[源代码]: 初始化弹性常数求解器

solve(strains, stresses, method='least_squares', alpha=1e-05)[源代码]

求解弹性常数矩阵

参数:

strains (numpy.ndarray) -- 应变数据，形状(N, 6)，N为数据点数
stresses (numpy.ndarray) -- 应力数据，形状(N, 6)
method (str, optional) -- 求解方法，支持'least_squares'和'ridge'，默认'least_squares'
alpha (float, optional) -- 岭回归正则化参数，仅在method='ridge'时使用，默认1e-5

返回:

弹性常数矩阵(6,6)和拟合优度R²

返回类型:

tuple[numpy.ndarray, float]

抛出:

ValueError -- 如果输入数据格式不正确或求解方法不支持

备注

支持的求解方法：

最小二乘法 ('least_squares')：标准线性回归，适用于大多数情况
岭回归 ('ridge')：加入L2正则化，适用于病态矩阵情况

拟合优度R²计算：

\[R^2 = 1 - \frac{SS_{res}}{SS_{tot}}\]

其中：

\(SS_{res} = \sum (y_i - \hat{y}_i)^2\): 残差平方和
\(SS_{tot} = \sum (y_i - \bar{y})^2\): 总平方和

示例

>>> solver = ElasticConstantsSolver()
>>> # 使用最小二乘法
>>> C, r2 = solver.solve(strains, stresses)
>>> # 使用岭回归（适用于病态情况）
>>> C, r2 = solver.solve(strains, stresses, method='ridge', alpha=1e-3)

class thermoelasticsim.elastic.deformation_method.zero_temp.ShearDeformationMethod[源代码]

基类：object

LAMMPS风格剪切形变方法

实现经过验证的LAMMPS风格剪切形变算法，用于计算剪切弹性常数(C44, C55, C66)。该方法通过晶胞剪切和原子位置调整来施加纯剪切应变，避免了传统形变矩阵方法在剪切计算中的数值问题。

apply_shear_deformation(cell, direction, strain_magnitude)[源代码]

对晶胞施加LAMMPS风格剪切形变

参数:

cell (Cell)
direction (int)
strain_magnitude (float)

返回类型:

calculate_c44_response(cell, potential, strain_amplitudes, relaxer)[源代码]

计算C44剪切响应

参数:

cell (Cell)
potential (Potential)
strain_amplitudes (ndarray)
relaxer (StructureRelaxer)
include_base_state (bool)

返回类型:

dict[str, Any]

备注

LAMMPS剪切形变的核心原理：

晶胞剪切：直接修改晶格矢量的非对角分量
原子位置调整：按照仿射变换调整所有原子位置
应变-应力关系：通过virial应力张量计算剪切应力

支持的剪切方向： 1. direction=4：yz 剪切 (σ23, C44) 2. direction=5：xz 剪切 (σ13, C55) 3. direction=6：xy 剪切 (σ12, C66)

剪切应变定义：

\[\gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\]

其中γ是工程剪切应变，与张量剪切应变的关系为 \(\gamma = 2\varepsilon\)

示例

>>> shear_method = ShearDeformationMethod()
>>>
>>> # 施加xy剪切形变
>>> deformed_cell = shear_method.apply_shear_deformation(
...     cell, direction=6, strain_magnitude=0.005
... )
>>>
>>> # 计算C44响应
>>> strains = np.linspace(-0.004, 0.004, 9)
>>> c44_data = shear_method.calculate_c44_response(
...     cell, potential, strains, relaxer
... )

__init__()[源代码]: 初始化剪切形变方法

apply_shear_deformation(cell, direction, strain_magnitude)[源代码]

对晶胞施加LAMMPS风格剪切形变

该方法实现了LAMMPS中的剪切形变算法，通过同时修改晶胞参数和原子位置来施加纯剪切应变。这种方法保持了晶胞的周期性边界条件，并确保剪切形变的一致性。

参数:

cell (Cell) -- 待形变的原始晶胞
direction (int) -- 剪切方向代码：4(yz), 5(xz), 6(xy)
strain_magnitude (float) -- 剪切应变幅度（工程剪切应变γ）

返回:

施加剪切形变后的新晶胞对象

返回类型:

抛出:

ValueError -- 如果指定了不支持的剪切方向

备注

LAMMPS剪切形变算法：

yz剪切 (direction=4)：

\[ \begin{align}\begin{aligned}h_{zy} \leftarrow h_{zy} + \gamma \cdot h_{zz}\\y_i \leftarrow y_i + \gamma \cdot z_i\end{aligned}\end{align} \]
xz剪切 (direction=5)：

\[ \begin{align}\begin{aligned}h_{zx} \leftarrow h_{zx} + \gamma \cdot h_{zz}\\x_i \leftarrow x_i + \gamma \cdot z_i\end{aligned}\end{align} \]
xy剪切 (direction=6)：

\[ \begin{align}\begin{aligned}h_{yx} \leftarrow h_{yx} + \gamma \cdot h_{yy}\\x_i \leftarrow x_i + \gamma \cdot y_i\end{aligned}\end{align} \]

其中：

h: 晶格矢量矩阵
γ: 工程剪切应变
(x_i, y_i, z_i): 第 i 个原子的坐标

示例

>>> method = ShearDeformationMethod()
>>>
>>> # 施加0.5%的xy剪切
>>> deformed = method.apply_shear_deformation(cell, 6, 0.005)
>>>
>>> # 验证剪切应变
>>> original_pos = cell.get_positions()
>>> deformed_pos = deformed.get_positions()
>>> shear_strain = np.mean(deformed_pos[:, 0] - original_pos[:, 0]) / np.mean(original_pos[:, 1])
>>> print(f"实际剪切应变: {shear_strain:.6f}")

calculate_c44_response(cell, potential, strain_amplitudes, relaxer, include_base_state=True)[源代码]

计算 C44 剪切响应。

参数:

cell (Cell) -- 基态晶胞对象。
potential (Potential) -- 势能函数实现。
strain_amplitudes (numpy.ndarray) -- 应变幅度数组，建议范围 [-0.005, 0.005]。
relaxer (StructureRelaxer) -- 结构弛豫器实例。
include_base_state (bool, optional) -- 是否包含零应变基态点，默认 True。

返回:

含方向列表、每方向详细数据与汇总拟合结果。

返回类型:

thermoelasticsim.elastic.deformation_method.zero_temp.calculate_zero_temp_elastic_constants(cell, potential, delta=0.005, num_steps=5, **kwargs)[源代码]

零温弹性常数计算的便捷函数

这是一个高级接口，封装了完整的零温弹性常数计算流程。适合快速计算和脚本使用。

参数:

cell (Cell) -- 晶胞对象
potential (Potential) -- 势能函数对象
delta (float, optional) -- 应变幅度，默认0.005 (0.5%)
num_steps (int, optional) -- 每个应变分量的步数，默认5
**kwargs -- 其他参数传递给ZeroTempDeformationCalculator

返回:

弹性常数矩阵(GPa)和拟合优度R²

返回类型:

tuple[numpy.ndarray, float]

示例

>>> from thermoelasticsim.elastic.deformation_method.zero_temp import calculate_zero_temp_elastic_constants
>>>
>>> # 快速计算
>>> C_matrix, r2 = calculate_zero_temp_elastic_constants(cell, potential)
>>> print(f"弹性常数矩阵(GPa):\\n{C_matrix}")
>>> print(f"拟合优度: {r2:.6f}")
>>>
>>> # 高精度计算
>>> C_matrix, r2 = calculate_zero_temp_elastic_constants(
...     cell, potential, delta=0.001, num_steps=1
... )

有限温计算

ThermoElasticSim - 有限温弹性常数计算模块

该模块实现有限温条件下，通过显式形变法计算弹性常数的工作流：

NPT 系综平衡
NVT 系综应力采样
弹性常数拟合

class thermoelasticsim.elastic.deformation_method.finite_temp.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.elastic.deformation_method.finite_temp.Deformer(delta=0.01, num_steps=10)[源代码]

基类：object

生成变形矩阵并应用于晶胞

参数:

delta (float)
num_steps (int)

delta

应变幅度，默认0.01

Type:: float

num_steps

每个应变分量的步数，默认10

Type:: int

logger

日志记录器

Type:: logging.Logger

STRAIN_COMPONENTS = [array([[1, 0, 0], [0, 0, 0], [0, 0, 0]]), array([[0, 0, 0], [0, 1, 0], [0, 0, 0]]), array([[0, 0, 0], [0, 0, 0], [0, 0, 1]]), array([[0, 0, 0], [0, 0, 1], [0, 1, 0]]), array([[0, 0, 1], [0, 0, 0], [1, 0, 0]]), array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])]

__init__(delta=0.01, num_steps=10)[源代码]

初始化Deformer

参数:

delta (float, optional) -- 应变幅度，建议范围0.001-0.05 (默认: 0.01)
num_steps (int, optional) -- 每个应变分量的步数，必须大于0 (默认: 10)

抛出:

ValueError -- 如果delta超出合理范围或num_steps非正

apply_deformation(cell, deformation_matrix)[源代码]

对晶胞施加形变矩阵 \(\mathbf{F}\)

参数:

cell (Cell) -- 包含原子的晶胞对象
deformation_matrix (numpy.ndarray) -- 3×3 形变矩阵

抛出:

ValueError -- 如果变形矩阵不是3x3矩阵或行列式接近零

返回类型:

None

generate_deformation_matrices()[源代码]

生成形变矩阵列表

为 6 个 Voigt 分量依次生成均匀分布的形变序列，采用：

\[\mathbf{F}(\delta) = \mathbf{I} + \delta\,\mathbf{E}_k\]

其中 \(\mathbf{E}_k\) 为第 \(k\) 个分量的基矩阵（剪切为对称填充）。

返回:: 形变矩阵列表，每个为 (3, 3)
返回类型:: list of numpy.ndarray

备注

预分配结果列表以减少内存分配开销。

class thermoelasticsim.elastic.deformation_method.finite_temp.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)[源代码]

已移除：有限差分应力路径不再提供，避免误用。

返回类型:: ndarray

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/Å³
返回类型:: numpy.ndarray

calculate_stress_basic(cell, potential)[源代码]

向后兼容别名 - 指向总应力方法

返回类型:: ndarray

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

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

总应力张量 (3, 3)，单位 eV/Å³

返回类型:

np.ndarray

calculate_virial_kinetic_stress(cell, potential)[源代码]

向后兼容别名 - 指向总应力方法

返回类型:: ndarray

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\) 的力

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

维里应力张量 (3, 3)，单位 eV/Å³

返回类型:

np.ndarray

compute_stress(cell, potential)[源代码]

计算应力张量（使用总应力作为主要方法）

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

3x3 应力张量矩阵（总应力 = 动能 + 维里）

返回类型:

np.ndarray

get_all_stress_components(cell, potential)[源代码]

计算应力张量的所有分量

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

包含应力张量分量的字典，键包括： - "kinetic": 动能应力张量 - "virial": 维里应力张量 - "total": 总应力张量（kinetic + virial）

返回类型:

Dict[str, np.ndarray]

validate_tensor_symmetry(tensor, tolerance=1e-10)[源代码]

验证应力张量是否对称

参数:

tensor (np.ndarray) -- 应力张量 (3x3)
tolerance (float, optional) -- 对称性的容差, 默认值为1e-10

返回:

如果对称则为True, 否则为False

返回类型:

class thermoelasticsim.elastic.deformation_method.finite_temp.AndersenBarostat(target_pressure, mass, temperature)[源代码]

基类：Barostat

Andersen 恒压器（已弃用）

自 4.0 版本弃用: 该旧架构恒压器已弃用；请使用新架构中的 thermoelasticsim.md.propagators.MTKBarostatPropagator。

参数:

target_pressure (float) -- 目标压力（eV/Å³）
mass (float) -- 活塞质量参数（内部单位）
temperature (float) -- 系统温度 (K)

__init__(target_pressure, mass, temperature)[源代码]

参数:

target_pressure (float)
mass (float)
temperature (float)

apply(cell, dt, potential)[源代码]

应用 Andersen 恒压器

参数:: dt (float)

class thermoelasticsim.elastic.deformation_method.finite_temp.MDSimulator(cell, potential, integrator, ensemble='NVE', thermostat=None, barostat=None)[源代码]

基类：object

分子动力学模拟器，支持NVE、NVT和NPT系综

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象
integrator (Integrator) -- 积分器对象
ensemble (str) -- 系综类型: 'NVE', 'NVT', 或 'NPT'
thermostat (Optional[Dict] or Thermostat) -- 恒温器配置或对象
barostat (Optional[Barostat]) -- 恒压器对象

__init__(cell, potential, integrator, ensemble='NVE', thermostat=None, barostat=None)[源代码]

calculate_total_energy()[源代码]

计算系统总能量

返回类型:: float

plot_pressure()[源代码]: 绘制压力演化图

plot_temperature()[源代码]: 绘制温度演化图

run(steps, dt)[源代码]

运行分子动力学模拟

参数:

steps (int) -- 模拟步数
dt (float) -- 时间步长

run_npt(steps, dt)[源代码]: 运行NPT系综模拟

run_nve(steps, dt)[源代码]: 运行NVE系综模拟

run_nvt(steps, dt)[源代码]: 运行NVT系综模拟

class thermoelasticsim.elastic.deformation_method.finite_temp.AndersenThermostat(target_temperature, collision_frequency)[源代码]

基类：Thermostat

Andersen 恒温器

通过随机碰撞来控制温度，实现符合正则系综的采样。该方法适合平衡态采样，但由于速度随机化，动力学轨迹会被打断。

参数:

target_temperature (float) -- 目标温度，单位为开尔文 (K)。
collision_frequency (float) -- 碰撞频率，单位为 fs⁻¹。较高的碰撞频率会更频繁地随机分配速度，导致温度波动较大；较低的碰撞频率则速度随机化较少，温度控制更稳定。

__init__(target_temperature, collision_frequency)[源代码]

参数:

target_temperature (float)
collision_frequency (float)

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

参数:

cell (Cell) -- 晶胞对象，包含原子集合
dt (float) -- 时间步长 (fs)
potential (Potential) -- 势函数对象（此实现中未直接使用）

返回类型:

None

class thermoelasticsim.elastic.deformation_method.finite_temp.TensorConverter[源代码]

基类：object

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

static from_voigt(voigt, tensor_type)[源代码]

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

参数:

voigt (np.ndarray) -- 形状为 (6,) 的 Voigt 向量。
tensor_type (str) -- 张量类型，必须是 'stress' 或 'strain'。

返回:

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

返回类型:

np.ndarray

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

class thermoelasticsim.elastic.deformation_method.finite_temp.FiniteTempElasticityWorkflow(cell, potential, temperature, integrator, deformation_delta=0.1, num_steps=10, npt_equilibration_steps=10000, nvt_sampling_steps=500000, time_step=1, save_path=None)[源代码]

基类：object

有限温弹性常数计算工作流

实现完整的有限温度下弹性常数计算流程，包括： - NPT系综平衡 - NVT系综应力采样 - 弹性常数拟合

参数:

cell (Cell) -- 初始晶胞对象
potential (Potential) -- 势能函数对象
temperature (float) -- 目标温度(K)
integrator (Integrator) -- 积分器对象
deformation_delta (float, optional) -- 最大变形量，默认为0.1
num_steps (int, optional) -- 变形步数，默认为10
npt_equilibration_steps (int, optional) -- NPT平衡步数，默认为10000
nvt_sampling_steps (int, optional) -- NVT采样步数，默认为500000
time_step (int, optional) -- 时间步长(fs)，默认为1
save_path (str, optional) -- 结果保存路径，默认为当前时间戳命名的目录

__init__(cell, potential, temperature, integrator, deformation_delta=0.1, num_steps=10, npt_equilibration_steps=10000, nvt_sampling_steps=500000, time_step=1, save_path=None)[源代码]

run_npt_equilibration(cell)[源代码]

在NPT系综下进行平衡

参数:: cell (Cell) -- 待平衡的晶胞对象
返回:: 平衡后的晶胞对象
返回类型:: Cell

run_nvt_sampling(cell)[源代码]

在 NVT 系综下采样应力

参数:

cell (Cell) -- 待采样的晶胞对象

返回:

avg_stress (numpy.ndarray) -- 平均应力张量 (3, 3)，单位 eV/Å³
std_stress (numpy.ndarray) -- 应力标准差 (3, 3)，单位 eV/Å³

返回类型:

tuple[ndarray, ndarray]

calculate_elastic_constants()[源代码]

计算弹性常数的主流程

执行完整工作流： 1. 生成变形矩阵 2. 对每个变形进行 NPT 平衡 3. 在 NVT 系综下采样应力 4. 使用最小二乘法拟合弹性常数

返回:: 弹性常数矩阵 (6, 6)，单位 GPa
返回类型:: numpy.ndarray

基础工具与材料参数

变形模块

用于生成并对晶胞施加小形变矩阵的工具。

理论基础

小形变近似下，形变梯度可写为：

\[\mathbf{F} = \mathbf{I} + \boldsymbol{\varepsilon}\]

其中 \(\boldsymbol{\varepsilon}\) 为对称的小应变张量。为便于逐分量扫描，本模块构造 6 个基方向（Voigt 记号）：\(\varepsilon_{11}, \varepsilon_{22}, \varepsilon_{33}, \varepsilon_{23}, \varepsilon_{13}, \varepsilon_{12}\)。

class thermoelasticsim.elastic.deformation.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.elastic.deformation.Deformer(delta=0.01, num_steps=10)[源代码]

基类：object

生成变形矩阵并应用于晶胞

参数:

delta (float)
num_steps (int)

delta

应变幅度，默认0.01

Type:: float

num_steps

每个应变分量的步数，默认10

Type:: int

logger

日志记录器

Type:: logging.Logger

__init__(delta=0.01, num_steps=10)[源代码]

初始化Deformer

参数:

delta (float, optional) -- 应变幅度，建议范围0.001-0.05 (默认: 0.01)
num_steps (int, optional) -- 每个应变分量的步数，必须大于0 (默认: 10)

抛出:

ValueError -- 如果delta超出合理范围或num_steps非正

STRAIN_COMPONENTS = [array([[1, 0, 0], [0, 0, 0], [0, 0, 0]]), array([[0, 0, 0], [0, 1, 0], [0, 0, 0]]), array([[0, 0, 0], [0, 0, 0], [0, 0, 1]]), array([[0, 0, 0], [0, 0, 1], [0, 1, 0]]), array([[0, 0, 1], [0, 0, 0], [1, 0, 0]]), array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])]

generate_deformation_matrices()[源代码]

生成形变矩阵列表

为 6 个 Voigt 分量依次生成均匀分布的形变序列，采用：

\[\mathbf{F}(\delta) = \mathbf{I} + \delta\,\mathbf{E}_k\]

其中 \(\mathbf{E}_k\) 为第 \(k\) 个分量的基矩阵（剪切为对称填充）。

返回:: 形变矩阵列表，每个为 (3, 3)
返回类型:: list of numpy.ndarray

备注

预分配结果列表以减少内存分配开销。

apply_deformation(cell, deformation_matrix)[源代码]

对晶胞施加形变矩阵 \(\mathbf{F}\)

参数:

cell (Cell) -- 包含原子的晶胞对象
deformation_matrix (numpy.ndarray) -- 3×3 形变矩阵

抛出:

ValueError -- 如果变形矩阵不是3x3矩阵或行列式接近零

返回类型:

None

材料参数配置模块

该模块定义了用于弹性常数计算的标准材料参数，包括晶格常数、原子质量、文献弹性常数等。支持多种晶体结构和材料元素。

主要组件：

MaterialParameters

材料参数数据类

预定义材料常量

ALUMINUM_FCC, COPPER_FCC 等

材料查询和验证工具

材料检索、枚举与常数比较

基本使用：

>>> from thermoelasticsim.elastic.materials import ALUMINUM_FCC, COPPER_FCC
>>> print(f"Al晶格常数: {ALUMINUM_FCC.lattice_constant:.3f} Å")
Al晶格常数: 4.050 Å

thermoelasticsim.elastic.materials.dataclass(cls=None, /, *, init=True, repr=True, eq=True, order=False, unsafe_hash=False, frozen=False, match_args=True, kw_only=False, slots=False, weakref_slot=False)[源代码]

Add dunder methods based on the fields defined in the class.

Examines PEP 526 __annotations__ to determine fields.

If init is true, an __init__() method is added to the class. If repr is true, a __repr__() method is added. If order is true, rich comparison dunder methods are added. If unsafe_hash is true, a __hash__() method is added. If frozen is true, fields may not be assigned to after instance creation. If match_args is true, the __match_args__ tuple is added. If kw_only is true, then by default all fields are keyword-only. If slots is true, a new class with a __slots__ attribute is returned.

class thermoelasticsim.elastic.materials.MaterialParameters(name, symbol, mass_amu, lattice_constant, structure, literature_elastic_constants, temperature=0.0, description='')[源代码]

基类：object

材料参数数据类

封装材料的所有基本物理参数，包括结构、热力学和弹性性质。使用frozen=True确保参数不可变，避免意外修改。

参数:

name (str)
symbol (str)
mass_amu (float)
lattice_constant (float)
structure (str)
literature_elastic_constants (dict[str, float])
temperature (float)
description (str)

name

材料全名，如"Aluminum"

Type:: str

symbol

化学元素符号，如"Al"

Type:: str

mass_amu

原子质量 (原子质量单位)

Type:: float

lattice_constant

晶格常数 (Å)

Type:: float

structure

晶体结构类型，如"fcc", "bcc", "hcp"

Type:: str

literature_elastic_constants

文献弹性常数参考值 (GPa) 包含键: "C11", "C12", "C44"等

Type:: dict

temperature

参考温度 (K)，默认为0K (零温)

Type:: float, optional

description

材料描述信息

Type:: str, optional

示例

创建自定义材料参数:

>>> gold_fcc = MaterialParameters(
...     name="Gold",
...     symbol="Au",
...     mass_amu=196.966569,
...     lattice_constant=4.08,
...     structure="fcc",
...     literature_elastic_constants={
...         "C11": 192.0,
...         "C12": 163.0,
...         "C44": 41.5
...     }
... )
>>> print(f"{gold_fcc.name}: C11 = {gold_fcc.literature_elastic_constants['C11']} GPa")
Gold: C11 = 192.0 GPa

备注

弹性常数命名约定：

C11, C22, C33: 正应力-正应变常数
C12, C13, C23: 泊松效应常数
C44, C55, C66: 剪切常数
立方晶系: C11=C22=C33, C44=C55=C66, C12=C13=C23

name: str

symbol: str

mass_amu: float

lattice_constant: float

structure: str

literature_elastic_constants: dict[str, float]

temperature: float = 0.0

description: str = ''

__post_init__()[源代码]: 参数验证

property bulk_modulus: float

计算体积模量

对于立方晶系: K = (C11 + 2*C12) / 3

返回:: 体积模量 (GPa)
返回类型:: float

property shear_modulus: float

计算剪切模量

对于立方晶系: G = C44 (Voigt平均的一种近似)

返回:: 剪切模量 (GPa)
返回类型:: float

property young_modulus: float

计算杨氏模量

使用关系: E = 9*K*G / (3*K + G) 其中K为体积模量，G为剪切模量

返回:: 杨氏模量 (GPa)
返回类型:: float

property poisson_ratio: float

计算泊松比

使用关系: ν = (3*K - 2*G) / (6*K + 2*G)

返回:: 泊松比 (无量纲)
返回类型:: float

property theoretical_density: float

基于晶体学参数估算材料密度（g/cm³）。

对立方晶系/金刚石结构，使用常规晶胞的原子数与晶格常数估算：

ρ = (N_atoms × m_atom) / V_cell

其中 m_atom = mass_amu × 1.66053906660e-24 g， V_cell = (a[Å] × 1e-8 cm/Å)^3。

返回:: 理论密度（g/cm³）。
返回类型:: float

备注

仅在结构为 "fcc" 或 "diamond" 时使用该公式；其他结构返回 NotImplementedError。
结果为近似理论值，忽略温度导致的体膨胀与缺陷等因素。

__init__(name, symbol, mass_amu, lattice_constant, structure, literature_elastic_constants, temperature=0.0, description='')

参数:

name (str)
symbol (str)
mass_amu (float)
lattice_constant (float)
structure (str)
literature_elastic_constants (dict[str, float])
temperature (float)
description (str)

返回类型:

None

thermoelasticsim.elastic.materials.get_material_by_symbol(symbol)[源代码]

根据元素符号获取材料参数

参数:: symbol (str) -- 化学元素符号，如"Al", "Cu"
返回:: 匹配的材料参数，未找到返回None
返回类型:: MaterialParameters or None

示例

>>> al_params = get_material_by_symbol("Al")
>>> if al_params:
...     print(f"找到材料: {al_params.name}")
找到材料: Aluminum

thermoelasticsim.elastic.materials.get_all_materials()[源代码]

获取所有预定义材料参数

返回:: 键为材料名称，值为MaterialParameters对象的字典
返回类型:: dict

示例

>>> materials = get_all_materials()
>>> for name, params in materials.items():
...     print(f"{name}: {params.structure.upper()}")
Aluminum: FCC
Copper: FCC
Gold: FCC

thermoelasticsim.elastic.materials.compare_elastic_constants(material1, material2)[源代码]

比较两种材料的弹性常数

参数:

material1 (MaterialParameters) -- 要比较的材料参数
material2 (MaterialParameters) -- 要比较的材料参数

返回:

包含比较结果的嵌套字典

返回类型:

示例

>>> al = ALUMINUM_FCC
>>> cu = COPPER_FCC
>>> comparison = compare_elastic_constants(al, cu)
>>> print(f"C11比值: {comparison['ratios']['C11']:.2f}")
C11比值: 0.65

材料常量

预定义示例：ALUMINUM_FCC, COPPER_FCC 等（详见模块源码）。为便于交叉引用，以下显式导出常量：

thermoelasticsim.elastic.materials.ALUMINUM_FCC

材料参数数据类

封装材料的所有基本物理参数，包括结构、热力学和弹性性质。使用frozen=True确保参数不可变，避免意外修改。

thermoelasticsim.elastic.materials.name

材料全名，如"Aluminum"

Type:: str

thermoelasticsim.elastic.materials.symbol

化学元素符号，如"Al"

Type:: str

thermoelasticsim.elastic.materials.mass_amu

原子质量 (原子质量单位)

Type:: float

thermoelasticsim.elastic.materials.lattice_constant

晶格常数 (Å)

Type:: float

thermoelasticsim.elastic.materials.structure

晶体结构类型，如"fcc", "bcc", "hcp"

Type:: str

thermoelasticsim.elastic.materials.literature_elastic_constants

文献弹性常数参考值 (GPa) 包含键: "C11", "C12", "C44"等

Type:: dict

thermoelasticsim.elastic.materials.temperature

参考温度 (K)，默认为0K (零温)

Type:: float, optional

thermoelasticsim.elastic.materials.description

材料描述信息

Type:: str, optional

示例

创建自定义材料参数:

>>> gold_fcc = MaterialParameters(
...     name="Gold",
...     symbol="Au",
...     mass_amu=196.966569,
...     lattice_constant=4.08,
...     structure="fcc",
...     literature_elastic_constants={
...         "C11": 192.0,
...         "C12": 163.0,
...         "C44": 41.5
...     }
... )
>>> print(f"{gold_fcc.name}: C11 = {gold_fcc.literature_elastic_constants['C11']} GPa")
Gold: C11 = 192.0 GPa

备注

弹性常数命名约定：

C11, C22, C33: 正应力-正应变常数
C12, C13, C23: 泊松效应常数
C44, C55, C66: 剪切常数
立方晶系: C11=C22=C33, C44=C55=C66, C12=C13=C23

thermoelasticsim.elastic.materials.COPPER_FCC

材料参数数据类

封装材料的所有基本物理参数，包括结构、热力学和弹性性质。使用frozen=True确保参数不可变，避免意外修改。

thermoelasticsim.elastic.materials.name

材料全名，如"Aluminum"

Type:: str

thermoelasticsim.elastic.materials.symbol

化学元素符号，如"Al"

Type:: str

thermoelasticsim.elastic.materials.mass_amu

原子质量 (原子质量单位)

Type:: float

thermoelasticsim.elastic.materials.lattice_constant

晶格常数 (Å)

Type:: float

thermoelasticsim.elastic.materials.structure

晶体结构类型，如"fcc", "bcc", "hcp"

Type:: str

thermoelasticsim.elastic.materials.literature_elastic_constants

文献弹性常数参考值 (GPa) 包含键: "C11", "C12", "C44"等

Type:: dict

thermoelasticsim.elastic.materials.temperature

参考温度 (K)，默认为0K (零温)

Type:: float, optional

thermoelasticsim.elastic.materials.description

材料描述信息

Type:: str, optional

示例

创建自定义材料参数:

>>> gold_fcc = MaterialParameters(
...     name="Gold",
...     symbol="Au",
...     mass_amu=196.966569,
...     lattice_constant=4.08,
...     structure="fcc",
...     literature_elastic_constants={
...         "C11": 192.0,
...         "C12": 163.0,
...         "C44": 41.5
...     }
... )
>>> print(f"{gold_fcc.name}: C11 = {gold_fcc.literature_elastic_constants['C11']} GPa")
Gold: C11 = 192.0 GPa

备注

弹性常数命名约定：

C11, C22, C33: 正应力-正应变常数
C12, C13, C23: 泊松效应常数
C44, C55, C66: 剪切常数
立方晶系: C11=C22=C33, C44=C55=C66, C12=C13=C23

thermoelasticsim.elastic.materials.GOLD_FCC

材料参数数据类

封装材料的所有基本物理参数，包括结构、热力学和弹性性质。使用frozen=True确保参数不可变，避免意外修改。

thermoelasticsim.elastic.materials.name

材料全名，如"Aluminum"

Type:: str

thermoelasticsim.elastic.materials.symbol

化学元素符号，如"Al"

Type:: str

thermoelasticsim.elastic.materials.mass_amu

原子质量 (原子质量单位)

Type:: float

thermoelasticsim.elastic.materials.lattice_constant

晶格常数 (Å)

Type:: float

thermoelasticsim.elastic.materials.structure

晶体结构类型，如"fcc", "bcc", "hcp"

Type:: str

thermoelasticsim.elastic.materials.literature_elastic_constants

文献弹性常数参考值 (GPa) 包含键: "C11", "C12", "C44"等

Type:: dict

thermoelasticsim.elastic.materials.temperature

参考温度 (K)，默认为0K (零温)

Type:: float, optional

thermoelasticsim.elastic.materials.description

材料描述信息

Type:: str, optional

示例

创建自定义材料参数:

>>> gold_fcc = MaterialParameters(
...     name="Gold",
...     symbol="Au",
...     mass_amu=196.966569,
...     lattice_constant=4.08,
...     structure="fcc",
...     literature_elastic_constants={
...         "C11": 192.0,
...         "C12": 163.0,
...         "C44": 41.5
...     }
... )
>>> print(f"{gold_fcc.name}: C11 = {gold_fcc.literature_elastic_constants['C11']} GPa")
Gold: C11 = 192.0 GPa

备注

弹性常数命名约定：

C11, C22, C33: 正应力-正应变常数
C12, C13, C23: 泊松效应常数
C44, C55, C66: 剪切常数
立方晶系: C11=C22=C33, C44=C55=C66, C12=C13=C23

基准与工作流

零温弹性常数基准工作流

目标： 1. 将示例脚本中经过验证的计算流程沉淀为可复用 API 2. 统一材料常量来源，减少硬编码 3. 尽可能保持输出与绘图一致（可视化与 CSV 导出可选）

主要接口： 1. calculate_c11_c12_traditional(cell, potential, relaxer, mat) 2. calculate_c44_lammps_shear(cell, potential, relaxer, mat)

注意：本模块不负责日志目录管理，调用方可按需配置日志与输出目录。

class thermoelasticsim.elastic.benchmark.Callable: 基类：object

thermoelasticsim.elastic.benchmark.dataclass(cls=None, /, *, init=True, repr=True, eq=True, order=False, unsafe_hash=False, frozen=False, match_args=True, kw_only=False, slots=False, weakref_slot=False)[源代码]

Add dunder methods based on the fields defined in the class.

Examines PEP 526 __annotations__ to determine fields.

If init is true, an __init__() method is added to the class. If repr is true, a __repr__() method is added. If order is true, rich comparison dunder methods are added. If unsafe_hash is true, a __hash__() method is added. If frozen is true, fields may not be assigned to after instance creation. If match_args is true, the __match_args__ tuple is added. If kw_only is true, then by default all fields are keyword-only. If slots is true, a new class with a __slots__ attribute is returned.

class thermoelasticsim.elastic.benchmark.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.elastic.benchmark.CrystallineStructureBuilder[源代码]

基类：object

统一的晶体结构生成器

提供多种晶格类型的标准化生成方法，支持任意超胞尺寸和材料参数。所有生成的结构都经过周期性边界条件优化，适用于分子动力学计算。

create_fcc(element, lattice_constant, supercell)[源代码]

创建面心立方(FCC)结构

参数:

element (str)
lattice_constant (float)
supercell (tuple[int, int, int])

返回类型:

create_bcc(element, lattice_constant, supercell)[源代码]

创建体心立方(BCC)结构 (预留)

参数:

element (str)
lattice_constant (float)
supercell (tuple[int, int, int])

返回类型:

create_hcp(element, lattice_constant, supercell)[源代码]

创建密排六方(HCP)结构 (预留)

参数:

element (str)
lattice_constant (float)
supercell (tuple[int, int, int])

返回类型:

示例

创建3×3×3的铝FCC超胞:

>>> builder = CrystallineStructureBuilder()
>>> al_cell = builder.create_fcc("Al", 4.05, (3, 3, 3))
>>> print(f"原子数: {al_cell.num_atoms}")
原子数: 108

创建不同尺寸的铜结构:

>>> cu_cell = builder.create_fcc("Cu", 3.615, (4, 4, 4))
>>> print(f"晶胞体积: {cu_cell.volume:.1f} Å³")
晶胞体积: 1885.4 Å³

ATOMIC_MASSES = {'Ag': 107.8682, 'Al': 26.9815, 'Ar': 39.95, 'Au': 196.966569, 'C': 12.0107, 'Cr': 51.9961, 'Cu': 63.546, 'Fe': 55.845, 'Mg': 24.305, 'Ni': 58.6934, 'Pb': 207.2, 'Ti': 47.867, 'V': 50.9415, 'Zn': 65.38}

__init__()[源代码]: 初始化晶体结构生成器

create_bcc(element, lattice_constant, supercell)[源代码]

创建体心立方(BCC)结构

BCC结构特点： - 原胞包含2个原子 - 基矢：(0,0,0), (0.5,0.5,0.5) - 配位数：8 - 致密度：68%

参数:

element (str) -- 元素符号，如"Fe"、"Cr"等
lattice_constant (float) -- 晶格常数 (Å)
supercell (tuple of int) -- 超胞尺寸 (nx, ny, nz)

返回:

包含BCC结构的晶胞对象

返回类型:

备注

此方法为预留接口，完整实现待后续版本。当前主要支持FCC结构的弹性常数计算。

create_diamond(element, lattice_constant, supercell)[源代码]

创建金刚石（diamond cubic）结构的常规立方胞（8原子/胞）并复制成超胞。

参数:

element (str) -- 元素符号，例如 "C"、"Si"。
lattice_constant (float) -- 常规立方胞边长 a（Å）。
supercell (tuple[int,int,int]) -- 超胞复制 (nx, ny, nz)。

返回:

周期性金刚石结构晶胞。

返回类型:

create_fcc(element, lattice_constant, supercell)[源代码]

创建面心立方(FCC)结构

FCC结构特点： - 原胞包含4个原子 - 基矢：(0,0,0), (0.5,0.5,0), (0.5,0,0.5), (0,0.5,0.5) - 配位数：12 - 致密度：74%

参数:

element (str) -- 元素符号，如"Al"、"Cu"等
lattice_constant (float) -- 晶格常数 (Å)
supercell (tuple of int) -- 超胞尺寸 (nx, ny, nz)

返回:

包含FCC结构的晶胞对象

返回类型:

抛出:

ValueError -- 如果元素不在支持列表中
TypeError -- 如果参数类型不正确

示例

创建标准Al FCC结构:

>>> builder = CrystallineStructureBuilder()
>>> cell = builder.create_fcc("Al", 4.05, (2, 2, 2))
>>> cell.num_atoms
32

备注

生成的结构特性：

原子数计算：总原子数 = 4 × nx × ny × nz
晶胞矢量：正交晶胞，a = b = c = lattice_constant
周期边界：默认启用，适用于无限系统模拟
原子编号：从0开始连续编号

create_hcp(element, lattice_constant, supercell)[源代码]

创建密排六方(HCP)结构

HCP结构特点： - 原胞包含2个原子 - 六方晶系，c/a ≈ 1.633 - 配位数：12 - 致密度：74%

参数:

element (str) -- 元素符号，如"Zn"、"Mg"等
lattice_constant (float) -- a轴晶格常数 (Å)
supercell (tuple of int) -- 超胞尺寸 (nx, ny, nz)

返回:

包含HCP结构的晶胞对象

返回类型:

备注

此方法为预留接口，完整实现待后续版本。当前主要支持FCC结构的弹性常数计算。

static get_atomic_mass(element)[源代码]

获取指定元素的原子质量

参数:: element (str) -- 元素符号
返回:: 原子质量 (amu)
返回类型:: float
抛出:: ValueError -- 如果元素不在支持列表中

static get_supported_elements()[源代码]

获取支持的元素列表

返回:: 支持的元素符号列表
返回类型:: list of str

class thermoelasticsim.elastic.benchmark.ShearDeformationMethod[源代码]

基类：object

LAMMPS风格剪切形变方法

实现经过验证的LAMMPS风格剪切形变算法，用于计算剪切弹性常数(C44, C55, C66)。该方法通过晶胞剪切和原子位置调整来施加纯剪切应变，避免了传统形变矩阵方法在剪切计算中的数值问题。

apply_shear_deformation(cell, direction, strain_magnitude)[源代码]

对晶胞施加LAMMPS风格剪切形变

参数:

cell (Cell)
direction (int)
strain_magnitude (float)

返回类型:

calculate_c44_response(cell, potential, strain_amplitudes, relaxer)[源代码]

计算C44剪切响应

参数:

cell (Cell)
potential (Potential)
strain_amplitudes (ndarray)
relaxer (StructureRelaxer)
include_base_state (bool)

返回类型:

dict[str, Any]

备注

LAMMPS剪切形变的核心原理：

晶胞剪切：直接修改晶格矢量的非对角分量
原子位置调整：按照仿射变换调整所有原子位置
应变-应力关系：通过virial应力张量计算剪切应力

支持的剪切方向： 1. direction=4：yz 剪切 (σ23, C44) 2. direction=5：xz 剪切 (σ13, C55) 3. direction=6：xy 剪切 (σ12, C66)

剪切应变定义：

\[\gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\]

其中γ是工程剪切应变，与张量剪切应变的关系为 \(\gamma = 2\varepsilon\)

示例

>>> shear_method = ShearDeformationMethod()
>>>
>>> # 施加xy剪切形变
>>> deformed_cell = shear_method.apply_shear_deformation(
...     cell, direction=6, strain_magnitude=0.005
... )
>>>
>>> # 计算C44响应
>>> strains = np.linspace(-0.004, 0.004, 9)
>>> c44_data = shear_method.calculate_c44_response(
...     cell, potential, strains, relaxer
... )

__init__()[源代码]: 初始化剪切形变方法

apply_shear_deformation(cell, direction, strain_magnitude)[源代码]

对晶胞施加LAMMPS风格剪切形变

该方法实现了LAMMPS中的剪切形变算法，通过同时修改晶胞参数和原子位置来施加纯剪切应变。这种方法保持了晶胞的周期性边界条件，并确保剪切形变的一致性。

参数:

cell (Cell) -- 待形变的原始晶胞
direction (int) -- 剪切方向代码：4(yz), 5(xz), 6(xy)
strain_magnitude (float) -- 剪切应变幅度（工程剪切应变γ）

返回:

施加剪切形变后的新晶胞对象

返回类型:

抛出:

ValueError -- 如果指定了不支持的剪切方向

备注

LAMMPS剪切形变算法：

yz剪切 (direction=4)：

\[ \begin{align}\begin{aligned}h_{zy} \leftarrow h_{zy} + \gamma \cdot h_{zz}\\y_i \leftarrow y_i + \gamma \cdot z_i\end{aligned}\end{align} \]
xz剪切 (direction=5)：

\[ \begin{align}\begin{aligned}h_{zx} \leftarrow h_{zx} + \gamma \cdot h_{zz}\\x_i \leftarrow x_i + \gamma \cdot z_i\end{aligned}\end{align} \]
xy剪切 (direction=6)：

\[ \begin{align}\begin{aligned}h_{yx} \leftarrow h_{yx} + \gamma \cdot h_{yy}\\x_i \leftarrow x_i + \gamma \cdot y_i\end{aligned}\end{align} \]

其中：

h: 晶格矢量矩阵
γ: 工程剪切应变
(x_i, y_i, z_i): 第 i 个原子的坐标

示例

>>> method = ShearDeformationMethod()
>>>
>>> # 施加0.5%的xy剪切
>>> deformed = method.apply_shear_deformation(cell, 6, 0.005)
>>>
>>> # 验证剪切应变
>>> original_pos = cell.get_positions()
>>> deformed_pos = deformed.get_positions()
>>> shear_strain = np.mean(deformed_pos[:, 0] - original_pos[:, 0]) / np.mean(original_pos[:, 1])
>>> print(f"实际剪切应变: {shear_strain:.6f}")

calculate_c44_response(cell, potential, strain_amplitudes, relaxer, include_base_state=True)[源代码]

计算 C44 剪切响应。

参数:

cell (Cell) -- 基态晶胞对象。
potential (Potential) -- 势能函数实现。
strain_amplitudes (numpy.ndarray) -- 应变幅度数组，建议范围 [-0.005, 0.005]。
relaxer (StructureRelaxer) -- 结构弛豫器实例。
include_base_state (bool, optional) -- 是否包含零应变基态点，默认 True。

返回:

含方向列表、每方向详细数据与汇总拟合结果。

返回类型:

class thermoelasticsim.elastic.benchmark.StructureRelaxer(optimizer_type='L-BFGS', optimizer_params=None, supercell_dims=None, trajectory_recorder=None)[源代码]

基类：object

结构弛豫计算引擎

负责执行结构优化计算，包括完全弛豫和内部弛豫两种模式。

参数:

optimizer_type (str, optional) -- 优化器类型，默认 'L-BFGS'
optimizer_params (dict, optional) -- 优化器参数字典
supercell_dims (tuple[int, int, int] | None)

optimizer_type

使用的优化器类型

Type:: str

optimizer_params

优化器配置参数

Type:: dict

备注

两种弛豫模式的区别：

完全弛豫: 同时优化原子位置和晶胞参数，用于制备无应力基态。
内部弛豫: 仅优化原子位置，保持晶胞形状固定，用于形变后的平衡。

示例

>>> relaxer = StructureRelaxer()
>>> # 完全弛豫制备基态
>>> converged = relaxer.full_relax(cell, potential)
>>> # 形变后内部弛豫
>>> converged = relaxer.internal_relax(deformed_cell, potential)

__init__(optimizer_type='L-BFGS', optimizer_params=None, supercell_dims=None, trajectory_recorder=None)[源代码]

初始化结构弛豫器

参数:

optimizer_type (str, optional) -- 优化器类型，支持 'L-BFGS', 'BFGS', 'GD'，默认 'L-BFGS'
optimizer_params (dict, optional) -- 优化器参数，如果为None则使用默认参数
supercell_dims (tuple, optional) -- 超胞维度(nx, ny, nz)，用于正确显示等效单胞参数
trajectory_recorder (ElasticTrajectoryRecorder, optional) -- 轨迹记录器，用于记录优化过程中的系统状态

抛出:

ValueError -- 如果指定了不支持的优化器类型

full_relax(cell, potential)[源代码]

执行完全弛豫：同时优化原子位置和晶胞参数

完全弛豫用于制备无应力基态，是零温弹性常数计算的第一步。通过同时优化原子位置和晶胞参数，消除系统内应力。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{r,h} E(r,h) \quad \text{s.t.} \quad \sigma = 0\]

符号说明： - \(r\) 原子位置 - \(h\) 晶胞参数 - \(E\) 总能量 - \(\sigma\) 应力张量

示例

>>> relaxer = StructureRelaxer()
>>> converged = relaxer.full_relax(cell, potential)
>>> if converged:
...     print("成功制备无应力基态")

internal_relax(cell, potential)[源代码]

执行内部弛豫：仅优化原子位置，保持晶胞形状固定

内部弛豫用于形变后的结构优化，在保持宏观形变的前提下，寻找原子的最优位置配置。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{r} E(r,h_{\text{fixed}}) \quad \text{s.t.} \quad h = \text{const}\]

符号说明： - \(r\) 原子位置（优化变量） - \(h_{\text{fixed}}\) 固定的晶胞参数 - \(E\) 总能量

示例

>>> relaxer = StructureRelaxer()
>>> # 先施加形变
>>> cell.apply_deformation(deformation_matrix)
>>> # 再内部弛豫
>>> converged = relaxer.internal_relax(cell, potential)

uniform_lattice_relax(cell, potential)[源代码]

等比例晶格弛豫：只优化晶格常数，保持原子相对位置不变

这种弛豫方式特别适合基态制备，既避免了原子位置的数值噪声，又保持了晶体结构的完美对称性，同时计算效率极高。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{s} E(r_{\text{scaled}}, s \cdot L) \quad \text{其中} \quad r_{\text{scaled}} = s \cdot r_{\text{frac}} \cdot L\]

符号说明： - \(s\) 等比例缩放因子（唯一优化变量） - \(r_{\text{frac}}\) 固定的分数坐标 - \(L\) 原始晶格矢量矩阵 - \(E\) 总能量

优势： 1. 自由度仅为 1，计算极快 2. 保持晶体结构完美对称性 3. 避免原子位置数值噪声 4. 适合各种晶系（不限于 FCC）

示例

>>> relaxer = StructureRelaxer()
>>> converged = relaxer.uniform_lattice_relax(cell, potential)
>>> if converged:
...     print("成功优化晶格常数，保持结构对称性")

class thermoelasticsim.elastic.benchmark.MaterialParameters(name, symbol, mass_amu, lattice_constant, structure, literature_elastic_constants, temperature=0.0, description='')[源代码]

基类：object

材料参数数据类

封装材料的所有基本物理参数，包括结构、热力学和弹性性质。使用frozen=True确保参数不可变，避免意外修改。

参数:

name (str)
symbol (str)
mass_amu (float)
lattice_constant (float)
structure (str)
literature_elastic_constants (dict[str, float])
temperature (float)
description (str)

name

材料全名，如"Aluminum"

Type:: str

symbol

化学元素符号，如"Al"

Type:: str

mass_amu

原子质量 (原子质量单位)

Type:: float

lattice_constant

晶格常数 (Å)

Type:: float

structure

晶体结构类型，如"fcc", "bcc", "hcp"

Type:: str

literature_elastic_constants

文献弹性常数参考值 (GPa) 包含键: "C11", "C12", "C44"等

Type:: dict

temperature

参考温度 (K)，默认为0K (零温)

Type:: float, optional

description

材料描述信息

Type:: str, optional

示例

创建自定义材料参数:

>>> gold_fcc = MaterialParameters(
...     name="Gold",
...     symbol="Au",
...     mass_amu=196.966569,
...     lattice_constant=4.08,
...     structure="fcc",
...     literature_elastic_constants={
...         "C11": 192.0,
...         "C12": 163.0,
...         "C44": 41.5
...     }
... )
>>> print(f"{gold_fcc.name}: C11 = {gold_fcc.literature_elastic_constants['C11']} GPa")
Gold: C11 = 192.0 GPa

备注

弹性常数命名约定：

C11, C22, C33: 正应力-正应变常数
C12, C13, C23: 泊松效应常数
C44, C55, C66: 剪切常数
立方晶系: C11=C22=C33, C44=C55=C66, C12=C13=C23

__init__(name, symbol, mass_amu, lattice_constant, structure, literature_elastic_constants, temperature=0.0, description='')

参数:

name (str)
symbol (str)
mass_amu (float)
lattice_constant (float)
structure (str)
literature_elastic_constants (dict[str, float])
temperature (float)
description (str)

返回类型:

None

__post_init__()[源代码]: 参数验证

property bulk_modulus: float

计算体积模量

对于立方晶系: K = (C11 + 2*C12) / 3

返回:: 体积模量 (GPa)
返回类型:: float

description: str = ''

property poisson_ratio: float

计算泊松比

使用关系: ν = (3*K - 2*G) / (6*K + 2*G)

返回:: 泊松比 (无量纲)
返回类型:: float

property shear_modulus: float

计算剪切模量

对于立方晶系: G = C44 (Voigt平均的一种近似)

返回:: 剪切模量 (GPa)
返回类型:: float

temperature: float = 0.0

property theoretical_density: float

基于晶体学参数估算材料密度（g/cm³）。

对立方晶系/金刚石结构，使用常规晶胞的原子数与晶格常数估算：

ρ = (N_atoms × m_atom) / V_cell

其中 m_atom = mass_amu × 1.66053906660e-24 g， V_cell = (a[Å] × 1e-8 cm/Å)^3。

返回:: 理论密度（g/cm³）。
返回类型:: float

备注

仅在结构为 "fcc" 或 "diamond" 时使用该公式；其他结构返回 NotImplementedError。
结果为近似理论值，忽略温度导致的体膨胀与缺陷等因素。

property young_modulus: float

计算杨氏模量

使用关系: E = 9*K*G / (3*K + G) 其中K为体积模量，G为剪切模量

返回:: 杨氏模量 (GPa)
返回类型:: float

name: str

symbol: str

mass_amu: float

lattice_constant: float

structure: str

literature_elastic_constants: dict[str, float]

class thermoelasticsim.elastic.benchmark.EAMAl1Potential(cutoff=6.5)[源代码]

基类：Potential

铝的嵌入式原子方法 (EAM) 势实现。

基于 Mendelev et al. (2008) Phil. Mag. 88(12) 的Al参数化。

参数:: cutoff (float, optional) -- 截断距离（Å），默认 6.5。

备注

绑定 C++ 后端 eam_al1；单位：能量 eV，长度 Å，力 eV/Å。
理论背景与方法见模块“References”。

__init__(cutoff=6.5)[源代码]

参数:: cutoff (float)

calculate_energy(cell, neighbor_list=None)[源代码]

使用EAM势计算系统的总势能。

参数:

cell (Cell) -- 包含原子信息的晶胞对象。
neighbor_list (NeighborList, optional) -- 在此实现中未使用，但为保持接口一致性而保留。

返回:

系统的总势能，单位 eV。

返回类型:

float

calculate_forces(cell, neighbor_list=None)[源代码]

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

参数:

cell (Cell) -- 包含原子信息的晶胞对象。
neighbor_list (NeighborList, optional) -- 在此实现中未使用，但为保持接口一致性而保留。

返回类型:

None

class thermoelasticsim.elastic.benchmark.ResponsePlotter(dpi=300, figsize_scale=1.0, literature_values=None)[源代码]

基类：object

弹性常数响应曲线绘制器

提供统一的绘图接口，支持： - C11/C12联合响应图 - C44/C55/C66剪切响应图 - 弹性常数对比图 - 收敛性分析图

示例

>>> plotter = ResponsePlotter()
>>> plotter.plot_c11_c12_response(c11_data, c12_data, 'output.png')
>>> plotter.plot_shear_response(c44_data, 'shear_output.png')

参数:

dpi (int)
figsize_scale (float)
literature_values (dict[str, float] | None)

__init__(dpi=300, figsize_scale=1.0, literature_values=None)[源代码]

初始化绘图器

参数:

dpi (int) -- 图像分辨率
figsize_scale (float) -- 图像尺寸缩放因子
literature_values (dict[str, float] | None)

plot_c11_c12_combined_response(c11_data, c12_data, supercell_size, output_path, slope_override_c11=None, slope_override_c12=None, subtitle_c11=None, subtitle_c12=None)[源代码]

生成C11/C12联合应力-应变响应关系图

从v7的plot_c11_c12_combined_response提取和重构

参数:

c11_data (List[Dict]) -- C11数据点列表
c12_data (List[Dict]) -- C12数据点列表
supercell_size (Tuple[int, int, int]) -- 系统尺寸
output_path (str) -- 输出文件路径
slope_override_c11 (float | None)
slope_override_c12 (float | None)
subtitle_c11 (str | None)
subtitle_c12 (str | None)

返回:

生成的图像文件名

返回类型:

plot_shear_response(detailed_results, supercell_size, output_path)[源代码]

生成C44/C55/C66剪切应力-应变响应关系图

从v7的plot_stress_strain_response提取和重构

参数:

detailed_results (List[Dict]) -- 详细的剪切响应结果
supercell_size (Tuple[int, int, int]) -- 系统尺寸
output_path (str) -- 输出文件路径

返回:

生成的图像文件名

返回类型:

class thermoelasticsim.elastic.benchmark.BenchmarkConfig(supercell_size=(3, 3, 3), delta=0.003, num_steps=1, shear_strains=(-0.004, -0.003, -0.002, -0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015, 0.002, 0.003, 0.004), optimizer_type='L-BFGS', optimizer_params=None, precision_mode=False, small_linear_strains=(-2e-05, -1e-05, 0.0, 1e-05, 2e-05), base_stress_tol_gpa=1e-05, base_relax_max_passes=3, small_shear_strains=(-0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015))[源代码]

基类：object

基准参数配置

参数:

supercell_size (tuple[int, int, int])
delta (float)
num_steps (int)
shear_strains (tuple[float, ...])
optimizer_type (str)
optimizer_params (dict[str, Any])
precision_mode (bool)
small_linear_strains (tuple[float, ...])
base_stress_tol_gpa (float)
base_relax_max_passes (int)
small_shear_strains (tuple[float, ...])

supercell_size: tuple[int, int, int] = (3, 3, 3)

delta: float = 0.003

num_steps: int = 1

shear_strains: tuple[float, ...] = (-0.004, -0.003, -0.002, -0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015, 0.002, 0.003, 0.004)

optimizer_type: str = 'L-BFGS'

optimizer_params: dict[str, Any] = None

precision_mode: bool = False

small_linear_strains: tuple[float, ...] = (-2e-05, -1e-05, 0.0, 1e-05, 2e-05)

base_stress_tol_gpa: float = 1e-05

base_relax_max_passes: int = 3

small_shear_strains: tuple[float, ...] = (-0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015)

build_relaxer()[源代码]

构建结构弛豫器

返回类型:: StructureRelaxer

__init__(supercell_size=(3, 3, 3), delta=0.003, num_steps=1, shear_strains=(-0.004, -0.003, -0.002, -0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015, 0.002, 0.003, 0.004), optimizer_type='L-BFGS', optimizer_params=None, precision_mode=False, small_linear_strains=(-2e-05, -1e-05, 0.0, 1e-05, 2e-05), base_stress_tol_gpa=1e-05, base_relax_max_passes=3, small_shear_strains=(-0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015))

参数:

supercell_size (tuple[int, int, int])
delta (float)
num_steps (int)
shear_strains (tuple[float, ...])
optimizer_type (str)
optimizer_params (dict[str, Any])
precision_mode (bool)
small_linear_strains (tuple[float, ...])
base_stress_tol_gpa (float)
base_relax_max_passes (int)
small_shear_strains (tuple[float, ...])

返回类型:

None

thermoelasticsim.elastic.benchmark.calculate_c11_c12_traditional(cell, potential, relaxer, material_params, strain_points=None, do_internal_relax=True)[源代码]

使用单轴应变法计算 C11/C12，测试点与示例保持一致。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (list[float] | ndarray | None)
do_internal_relax (bool)

返回类型:

thermoelasticsim.elastic.benchmark.calculate_c44_lammps_shear(cell, potential, relaxer, material_params, strain_points)[源代码]

使用LAMMPS风格剪切方法计算C44/C55/C66。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (ndarray | list[float])

返回类型:

thermoelasticsim.elastic.benchmark.calculate_c11_c12_robust(cell, potential, relaxer, material_params, strain_points=None, do_internal_relax=True)[源代码]

稳健法计算 C11/C12：正交应变(δ,-δ,0) + 等压应变(η,η,η)。

正交：σxx = (C11 - C12)·δ, σyy = -(C11 - C12)·δ
等压：σxx = (C11 + 2C12)·η = 3K·η

通过两次一维拟合解出 C11、C12，通常比单轴法对势细节更鲁棒。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (list[float] | ndarray | None)
do_internal_relax (bool)

返回类型:

thermoelasticsim.elastic.benchmark.calculate_c11_c12_biaxial_orthorhombic(cell, potential, relaxer, material_params, strain_points=None, do_internal_relax=True)[源代码]

通过两种无体积拟合获取 C11/C12：

正交: F = diag(1+δ, 1-δ, 1) → 斜率 Δ = C11 - C12
双轴: F = diag(1+ε, 1+ε, 1) → 斜率 Σ = C11 + C12

最终：C11 = (Σ+Δ)/2, C12 = (Σ-Δ)/2 避免直接使用等压路径对体积敏感的数值漂移。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (list[float] | ndarray | None)
do_internal_relax (bool)

返回类型:

thermoelasticsim.elastic.benchmark.run_zero_temp_benchmark(material_params, potential, supercell_size=(3, 3, 3), output_dir=None, save_json=True, precision=False, log_level=None, *, optimizer_type=None, uniaxial_strains=None, shear_strains=None, optimizer_params=None)[源代码]

运行通用零温弹性基准（结构自适应）。

根据 material_params.structure 自动生成晶体结构（当前支持 fcc、diamond）。
其它流程保持与铝基准一致：弛豫、C11/C12 单轴、C44 剪切、绘图与CSV/JSON。

参数:

material_params (MaterialParameters)
supercell_size (tuple[int, int, int])
output_dir (str | None)
save_json (bool)
precision (bool)
log_level (int | None)
optimizer_type (str | None)
uniaxial_strains (list[float] | tuple[float, ...] | None)
shear_strains (list[float] | tuple[float, ...] | None)
optimizer_params (dict[str, Any] | None)

返回类型:

thermoelasticsim.elastic.benchmark.run_size_sweep(sizes=None, output_root=None, material_params=None, potential_factory=None, precision=False)[源代码]

依次运行 2x2x2、3x3x3、4x4x4 超胞的有限尺寸效应基准。

默认材料为铝（ALUMINUM_FCC），默认势为 EAMAl1Potential。
可通过传入 material_params 与 potential_factory 进行泛化。

参数:

sizes (list[tuple[int, int, int]] | None)
output_root (str | None)
material_params (MaterialParameters | None)
potential_factory (Callable[[], Any] | None)
precision (bool)

返回类型:

list[dict]

顶层导出（便于交叉引用）

弹性常数计算模块

class thermoelasticsim.elastic.Deformer(delta=0.01, num_steps=10)[源代码]

基类：object

生成变形矩阵并应用于晶胞

参数:

delta (float)
num_steps (int)

delta

应变幅度，默认0.01

Type:: float

num_steps

每个应变分量的步数，默认10

Type:: int

logger

日志记录器

Type:: logging.Logger

STRAIN_COMPONENTS = [array([[1, 0, 0], [0, 0, 0], [0, 0, 0]]), array([[0, 0, 0], [0, 1, 0], [0, 0, 0]]), array([[0, 0, 0], [0, 0, 0], [0, 0, 1]]), array([[0, 0, 0], [0, 0, 1], [0, 1, 0]]), array([[0, 0, 1], [0, 0, 0], [1, 0, 0]]), array([[0, 1, 0], [1, 0, 0], [0, 0, 0]])]

__init__(delta=0.01, num_steps=10)[源代码]

初始化Deformer

参数:

delta (float, optional) -- 应变幅度，建议范围0.001-0.05 (默认: 0.01)
num_steps (int, optional) -- 每个应变分量的步数，必须大于0 (默认: 10)

抛出:

ValueError -- 如果delta超出合理范围或num_steps非正

apply_deformation(cell, deformation_matrix)[源代码]

对晶胞施加形变矩阵 \(\mathbf{F}\)

参数:

cell (Cell) -- 包含原子的晶胞对象
deformation_matrix (numpy.ndarray) -- 3×3 形变矩阵

抛出:

ValueError -- 如果变形矩阵不是3x3矩阵或行列式接近零

返回类型:

None

generate_deformation_matrices()[源代码]

生成形变矩阵列表

为 6 个 Voigt 分量依次生成均匀分布的形变序列，采用：

\[\mathbf{F}(\delta) = \mathbf{I} + \delta\,\mathbf{E}_k\]

其中 \(\mathbf{E}_k\) 为第 \(k\) 个分量的基矩阵（剪切为对称填充）。

返回:: 形变矩阵列表，每个为 (3, 3)
返回类型:: list of numpy.ndarray

备注

预分配结果列表以减少内存分配开销。

class thermoelasticsim.elastic.StrainCalculator[源代码]

基类：object

应变计算器类

参数:: F (numpy.ndarray) -- 3x3 变形矩阵

compute_strain(F)[源代码]

计算应变张量并返回 Voigt 表示法

参数:: F (numpy.ndarray) -- 3x3 变形矩阵
返回:: 应变向量，形状为 (6,)
返回类型:: numpy.ndarray

class thermoelasticsim.elastic.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)[源代码]

已移除：有限差分应力路径不再提供，避免误用。

返回类型:: ndarray

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/Å³
返回类型:: numpy.ndarray

calculate_stress_basic(cell, potential)[源代码]

向后兼容别名 - 指向总应力方法

返回类型:: ndarray

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

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

总应力张量 (3, 3)，单位 eV/Å³

返回类型:

np.ndarray

calculate_virial_kinetic_stress(cell, potential)[源代码]

向后兼容别名 - 指向总应力方法

返回类型:: ndarray

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\) 的力

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

维里应力张量 (3, 3)，单位 eV/Å³

返回类型:

np.ndarray

compute_stress(cell, potential)[源代码]

计算应力张量（使用总应力作为主要方法）

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

3x3 应力张量矩阵（总应力 = 动能 + 维里）

返回类型:

np.ndarray

get_all_stress_components(cell, potential)[源代码]

计算应力张量的所有分量

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

包含应力张量分量的字典，键包括： - "kinetic": 动能应力张量 - "virial": 维里应力张量 - "total": 总应力张量（kinetic + virial）

返回类型:

Dict[str, np.ndarray]

validate_tensor_symmetry(tensor, tolerance=1e-10)[源代码]

验证应力张量是否对称

参数:

tensor (np.ndarray) -- 应力张量 (3x3)
tolerance (float, optional) -- 对称性的容差, 默认值为1e-10

返回:

如果对称则为True, 否则为False

返回类型:

thermoelasticsim.elastic.ElasticConstantsWorkflow: ZeroTempDeformationCalculator 的别名

class thermoelasticsim.elastic.ElasticConstantsSolver[源代码]

基类：object

弹性常数求解器

从应力应变数据通过线性回归求解弹性常数矩阵。支持最小二乘法和岭回归两种方法，并提供拟合质量评估。

solve(strains, stresses, method='least_squares')[源代码]

求解弹性常数矩阵

参数:

strains (ndarray)
stresses (ndarray)
method (str)
alpha (float)

返回类型:

tuple[ndarray, float]

备注

数学原理：

根据胡克定律：\(\sigma = C \cdot \varepsilon\)

其中：

\(\sigma\): 应力向量 (N×6)
\(C\): 弹性常数矩阵 (6×6)
\(\varepsilon\): 应变向量 (N×6)

最小二乘求解：

\[\min_C ||\sigma - C \cdot \varepsilon||^2\]

解析解：

\[C = (\varepsilon^T \varepsilon)^{-1} \varepsilon^T \sigma\]

示例

>>> solver = ElasticConstantsSolver()
>>> C_matrix, r2_score = solver.solve(strains, stresses)
>>> print(f"弹性常数矩阵:\\n{C_matrix}")
>>> print(f"拟合优度: {r2_score:.6f}")

__init__()[源代码]: 初始化弹性常数求解器

solve(strains, stresses, method='least_squares', alpha=1e-05)[源代码]

求解弹性常数矩阵

参数:

strains (numpy.ndarray) -- 应变数据，形状(N, 6)，N为数据点数
stresses (numpy.ndarray) -- 应力数据，形状(N, 6)
method (str, optional) -- 求解方法，支持'least_squares'和'ridge'，默认'least_squares'
alpha (float, optional) -- 岭回归正则化参数，仅在method='ridge'时使用，默认1e-5

返回:

弹性常数矩阵(6,6)和拟合优度R²

返回类型:

tuple[numpy.ndarray, float]

抛出:

ValueError -- 如果输入数据格式不正确或求解方法不支持

备注

支持的求解方法：

最小二乘法 ('least_squares')：标准线性回归，适用于大多数情况
岭回归 ('ridge')：加入L2正则化，适用于病态矩阵情况

拟合优度R²计算：

\[R^2 = 1 - \frac{SS_{res}}{SS_{tot}}\]

其中：

\(SS_{res} = \sum (y_i - \hat{y}_i)^2\): 残差平方和
\(SS_{tot} = \sum (y_i - \bar{y})^2\): 总平方和

示例

>>> solver = ElasticConstantsSolver()
>>> # 使用最小二乘法
>>> C, r2 = solver.solve(strains, stresses)
>>> # 使用岭回归（适用于病态情况）
>>> C, r2 = solver.solve(strains, stresses, method='ridge', alpha=1e-3)

class thermoelasticsim.elastic.StructureRelaxer(optimizer_type='L-BFGS', optimizer_params=None, supercell_dims=None, trajectory_recorder=None)[源代码]

基类：object

结构弛豫计算引擎

负责执行结构优化计算，包括完全弛豫和内部弛豫两种模式。

参数:

optimizer_type (str, optional) -- 优化器类型，默认 'L-BFGS'
optimizer_params (dict, optional) -- 优化器参数字典
supercell_dims (tuple[int, int, int] | None)

optimizer_type

使用的优化器类型

Type:: str

optimizer_params

优化器配置参数

Type:: dict

备注

两种弛豫模式的区别：

完全弛豫: 同时优化原子位置和晶胞参数，用于制备无应力基态。
内部弛豫: 仅优化原子位置，保持晶胞形状固定，用于形变后的平衡。

示例

>>> relaxer = StructureRelaxer()
>>> # 完全弛豫制备基态
>>> converged = relaxer.full_relax(cell, potential)
>>> # 形变后内部弛豫
>>> converged = relaxer.internal_relax(deformed_cell, potential)

__init__(optimizer_type='L-BFGS', optimizer_params=None, supercell_dims=None, trajectory_recorder=None)[源代码]

初始化结构弛豫器

参数:

optimizer_type (str, optional) -- 优化器类型，支持 'L-BFGS', 'BFGS', 'GD'，默认 'L-BFGS'
optimizer_params (dict, optional) -- 优化器参数，如果为None则使用默认参数
supercell_dims (tuple, optional) -- 超胞维度(nx, ny, nz)，用于正确显示等效单胞参数
trajectory_recorder (ElasticTrajectoryRecorder, optional) -- 轨迹记录器，用于记录优化过程中的系统状态

抛出:

ValueError -- 如果指定了不支持的优化器类型

full_relax(cell, potential)[源代码]

执行完全弛豫：同时优化原子位置和晶胞参数

完全弛豫用于制备无应力基态，是零温弹性常数计算的第一步。通过同时优化原子位置和晶胞参数，消除系统内应力。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{r,h} E(r,h) \quad \text{s.t.} \quad \sigma = 0\]

符号说明： - \(r\) 原子位置 - \(h\) 晶胞参数 - \(E\) 总能量 - \(\sigma\) 应力张量

示例

>>> relaxer = StructureRelaxer()
>>> converged = relaxer.full_relax(cell, potential)
>>> if converged:
...     print("成功制备无应力基态")

internal_relax(cell, potential)[源代码]

执行内部弛豫：仅优化原子位置，保持晶胞形状固定

内部弛豫用于形变后的结构优化，在保持宏观形变的前提下，寻找原子的最优位置配置。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{r} E(r,h_{\text{fixed}}) \quad \text{s.t.} \quad h = \text{const}\]

符号说明： - \(r\) 原子位置（优化变量） - \(h_{\text{fixed}}\) 固定的晶胞参数 - \(E\) 总能量

示例

>>> relaxer = StructureRelaxer()
>>> # 先施加形变
>>> cell.apply_deformation(deformation_matrix)
>>> # 再内部弛豫
>>> converged = relaxer.internal_relax(cell, potential)

uniform_lattice_relax(cell, potential)[源代码]

等比例晶格弛豫：只优化晶格常数，保持原子相对位置不变

这种弛豫方式特别适合基态制备，既避免了原子位置的数值噪声，又保持了晶体结构的完美对称性，同时计算效率极高。

参数:

cell (Cell) -- 待优化的晶胞对象，会被直接修改
potential (Potential) -- 势能函数对象

返回:

优化是否成功收敛

返回类型:

备注

数学表述：

\[\min_{s} E(r_{\text{scaled}}, s \cdot L) \quad \text{其中} \quad r_{\text{scaled}} = s \cdot r_{\text{frac}} \cdot L\]

符号说明： - \(s\) 等比例缩放因子（唯一优化变量） - \(r_{\text{frac}}\) 固定的分数坐标 - \(L\) 原始晶格矢量矩阵 - \(E\) 总能量

优势： 1. 自由度仅为 1，计算极快 2. 保持晶体结构完美对称性 3. 避免原子位置数值噪声 4. 适合各种晶系（不限于 FCC）

示例

>>> relaxer = StructureRelaxer()
>>> converged = relaxer.uniform_lattice_relax(cell, potential)
>>> if converged:
...     print("成功优化晶格常数，保持结构对称性")

class thermoelasticsim.elastic.ZeroTempDeformationCalculator(cell, potential, delta=0.005, num_steps=5, relaxer_params=None, supercell_dims=None, base_relax_mode='uniform')[源代码]

基类：object

零温显式形变法计算器。

负责从基态制备到弹性常数求解的完整流程管理。

参数:

cell (Cell) -- 待计算的晶胞对象。
potential (Potential) -- 势能函数对象。
delta (float, optional) -- 应变幅度，默认 0.005 (0.5%)。
num_steps (int, optional) -- 每个应变分量的步数，默认 5（教学模式）。
relaxer_params (dict, optional) -- 结构弛豫器参数。
supercell_dims (tuple[int, int, int], optional) -- 超胞维度 (nx, ny, nz)。
base_relax_mode ({"uniform", "full"}, optional) -- 基态制备模式，等比例缩放或变胞+位置。

__init__(cell, potential, delta=0.005, num_steps=5, relaxer_params=None, supercell_dims=None, base_relax_mode='uniform')[源代码]

初始化零温形变计算器

参数:

cell (Cell) -- 晶胞对象
potential (Potential) -- 势能函数对象
delta (float, optional) -- 应变幅度，建议范围0.001-0.01，默认0.005
num_steps (int, optional) -- 每个应变分量的步数，1为生产模式，>1为教学模式，默认5
relaxer_params (dict, optional) -- 传递给StructureRelaxer的参数
supercell_dims (tuple, optional) -- 超胞维度(nx, ny, nz)，用于正确显示等效单胞参数
base_relax_mode (str)

抛出:

ValueError -- 如果delta或num_steps参数不合理

calculate()[源代码]

执行完整的零温弹性常数计算。

返回:: 弹性常数矩阵 (GPa) 和拟合优度 R^2。
返回类型:: tuple[numpy.ndarray, float]

class thermoelasticsim.elastic.FiniteTempElasticityWorkflow(cell, potential, temperature, integrator, deformation_delta=0.1, num_steps=10, npt_equilibration_steps=10000, nvt_sampling_steps=500000, time_step=1, save_path=None)[源代码]

基类：object

有限温弹性常数计算工作流

实现完整的有限温度下弹性常数计算流程，包括： - NPT系综平衡 - NVT系综应力采样 - 弹性常数拟合

参数:

cell (Cell) -- 初始晶胞对象
potential (Potential) -- 势能函数对象
temperature (float) -- 目标温度(K)
integrator (Integrator) -- 积分器对象
deformation_delta (float, optional) -- 最大变形量，默认为0.1
num_steps (int, optional) -- 变形步数，默认为10
npt_equilibration_steps (int, optional) -- NPT平衡步数，默认为10000
nvt_sampling_steps (int, optional) -- NVT采样步数，默认为500000
time_step (int, optional) -- 时间步长(fs)，默认为1
save_path (str, optional) -- 结果保存路径，默认为当前时间戳命名的目录

__init__(cell, potential, temperature, integrator, deformation_delta=0.1, num_steps=10, npt_equilibration_steps=10000, nvt_sampling_steps=500000, time_step=1, save_path=None)[源代码]

calculate_elastic_constants()[源代码]

计算弹性常数的主流程

执行完整工作流： 1. 生成变形矩阵 2. 对每个变形进行 NPT 平衡 3. 在 NVT 系综下采样应力 4. 使用最小二乘法拟合弹性常数

返回:: 弹性常数矩阵 (6, 6)，单位 GPa
返回类型:: numpy.ndarray

run_npt_equilibration(cell)[源代码]

在NPT系综下进行平衡

参数:: cell (Cell) -- 待平衡的晶胞对象
返回:: 平衡后的晶胞对象
返回类型:: Cell

run_nvt_sampling(cell)[源代码]

在 NVT 系综下采样应力

参数:

cell (Cell) -- 待采样的晶胞对象

返回:

avg_stress (numpy.ndarray) -- 平均应力张量 (3, 3)，单位 eV/Å³
std_stress (numpy.ndarray) -- 应力标准差 (3, 3)，单位 eV/Å³

返回类型:

tuple[ndarray, ndarray]

class thermoelasticsim.elastic.ShearDeformationMethod[源代码]

基类：object

LAMMPS风格剪切形变方法

实现经过验证的LAMMPS风格剪切形变算法，用于计算剪切弹性常数(C44, C55, C66)。该方法通过晶胞剪切和原子位置调整来施加纯剪切应变，避免了传统形变矩阵方法在剪切计算中的数值问题。

apply_shear_deformation(cell, direction, strain_magnitude)[源代码]

对晶胞施加LAMMPS风格剪切形变

参数:

cell (Cell)
direction (int)
strain_magnitude (float)

返回类型:

calculate_c44_response(cell, potential, strain_amplitudes, relaxer)[源代码]

计算C44剪切响应

参数:

cell (Cell)
potential (Potential)
strain_amplitudes (ndarray)
relaxer (StructureRelaxer)
include_base_state (bool)

返回类型:

dict[str, Any]

备注

LAMMPS剪切形变的核心原理：

晶胞剪切：直接修改晶格矢量的非对角分量
原子位置调整：按照仿射变换调整所有原子位置
应变-应力关系：通过virial应力张量计算剪切应力

支持的剪切方向： 1. direction=4：yz 剪切 (σ23, C44) 2. direction=5：xz 剪切 (σ13, C55) 3. direction=6：xy 剪切 (σ12, C66)

剪切应变定义：

\[\gamma_{ij} = \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\]

其中γ是工程剪切应变，与张量剪切应变的关系为 \(\gamma = 2\varepsilon\)

示例

>>> shear_method = ShearDeformationMethod()
>>>
>>> # 施加xy剪切形变
>>> deformed_cell = shear_method.apply_shear_deformation(
...     cell, direction=6, strain_magnitude=0.005
... )
>>>
>>> # 计算C44响应
>>> strains = np.linspace(-0.004, 0.004, 9)
>>> c44_data = shear_method.calculate_c44_response(
...     cell, potential, strains, relaxer
... )

__init__()[源代码]: 初始化剪切形变方法

apply_shear_deformation(cell, direction, strain_magnitude)[源代码]

对晶胞施加LAMMPS风格剪切形变

该方法实现了LAMMPS中的剪切形变算法，通过同时修改晶胞参数和原子位置来施加纯剪切应变。这种方法保持了晶胞的周期性边界条件，并确保剪切形变的一致性。

参数:

cell (Cell) -- 待形变的原始晶胞
direction (int) -- 剪切方向代码：4(yz), 5(xz), 6(xy)
strain_magnitude (float) -- 剪切应变幅度（工程剪切应变γ）

返回:

施加剪切形变后的新晶胞对象

返回类型:

抛出:

ValueError -- 如果指定了不支持的剪切方向

备注

LAMMPS剪切形变算法：

yz剪切 (direction=4)：

\[ \begin{align}\begin{aligned}h_{zy} \leftarrow h_{zy} + \gamma \cdot h_{zz}\\y_i \leftarrow y_i + \gamma \cdot z_i\end{aligned}\end{align} \]
xz剪切 (direction=5)：

\[ \begin{align}\begin{aligned}h_{zx} \leftarrow h_{zx} + \gamma \cdot h_{zz}\\x_i \leftarrow x_i + \gamma \cdot z_i\end{aligned}\end{align} \]
xy剪切 (direction=6)：

\[ \begin{align}\begin{aligned}h_{yx} \leftarrow h_{yx} + \gamma \cdot h_{yy}\\x_i \leftarrow x_i + \gamma \cdot y_i\end{aligned}\end{align} \]

其中：

h: 晶格矢量矩阵
γ: 工程剪切应变
(x_i, y_i, z_i): 第 i 个原子的坐标

示例

>>> method = ShearDeformationMethod()
>>>
>>> # 施加0.5%的xy剪切
>>> deformed = method.apply_shear_deformation(cell, 6, 0.005)
>>>
>>> # 验证剪切应变
>>> original_pos = cell.get_positions()
>>> deformed_pos = deformed.get_positions()
>>> shear_strain = np.mean(deformed_pos[:, 0] - original_pos[:, 0]) / np.mean(original_pos[:, 1])
>>> print(f"实际剪切应变: {shear_strain:.6f}")

calculate_c44_response(cell, potential, strain_amplitudes, relaxer, include_base_state=True)[源代码]

计算 C44 剪切响应。

参数:

cell (Cell) -- 基态晶胞对象。
potential (Potential) -- 势能函数实现。
strain_amplitudes (numpy.ndarray) -- 应变幅度数组，建议范围 [-0.005, 0.005]。
relaxer (StructureRelaxer) -- 结构弛豫器实例。
include_base_state (bool, optional) -- 是否包含零应变基态点，默认 True。

返回:

含方向列表、每方向详细数据与汇总拟合结果。

返回类型:

class thermoelasticsim.elastic.MaterialParameters(name, symbol, mass_amu, lattice_constant, structure, literature_elastic_constants, temperature=0.0, description='')[源代码]

基类：object

材料参数数据类

封装材料的所有基本物理参数，包括结构、热力学和弹性性质。使用frozen=True确保参数不可变，避免意外修改。

参数:

name (str)
symbol (str)
mass_amu (float)
lattice_constant (float)
structure (str)
literature_elastic_constants (dict[str, float])
temperature (float)
description (str)

name

材料全名，如"Aluminum"

Type:: str

symbol

化学元素符号，如"Al"

Type:: str

mass_amu

原子质量 (原子质量单位)

Type:: float

lattice_constant

晶格常数 (Å)

Type:: float

structure

晶体结构类型，如"fcc", "bcc", "hcp"

Type:: str

literature_elastic_constants

文献弹性常数参考值 (GPa) 包含键: "C11", "C12", "C44"等

Type:: dict

temperature

参考温度 (K)，默认为0K (零温)

Type:: float, optional

description

材料描述信息

Type:: str, optional

示例

创建自定义材料参数:

>>> gold_fcc = MaterialParameters(
...     name="Gold",
...     symbol="Au",
...     mass_amu=196.966569,
...     lattice_constant=4.08,
...     structure="fcc",
...     literature_elastic_constants={
...         "C11": 192.0,
...         "C12": 163.0,
...         "C44": 41.5
...     }
... )
>>> print(f"{gold_fcc.name}: C11 = {gold_fcc.literature_elastic_constants['C11']} GPa")
Gold: C11 = 192.0 GPa

备注

弹性常数命名约定：

C11, C22, C33: 正应力-正应变常数
C12, C13, C23: 泊松效应常数
C44, C55, C66: 剪切常数
立方晶系: C11=C22=C33, C44=C55=C66, C12=C13=C23

__init__(name, symbol, mass_amu, lattice_constant, structure, literature_elastic_constants, temperature=0.0, description='')

参数:

name (str)
symbol (str)
mass_amu (float)
lattice_constant (float)
structure (str)
literature_elastic_constants (dict[str, float])
temperature (float)
description (str)

返回类型:

None

__post_init__()[源代码]: 参数验证

property bulk_modulus: float

计算体积模量

对于立方晶系: K = (C11 + 2*C12) / 3

返回:: 体积模量 (GPa)
返回类型:: float

description: str = ''

property poisson_ratio: float

计算泊松比

使用关系: ν = (3*K - 2*G) / (6*K + 2*G)

返回:: 泊松比 (无量纲)
返回类型:: float

property shear_modulus: float

计算剪切模量

对于立方晶系: G = C44 (Voigt平均的一种近似)

返回:: 剪切模量 (GPa)
返回类型:: float

temperature: float = 0.0

property theoretical_density: float

基于晶体学参数估算材料密度（g/cm³）。

对立方晶系/金刚石结构，使用常规晶胞的原子数与晶格常数估算：

ρ = (N_atoms × m_atom) / V_cell

其中 m_atom = mass_amu × 1.66053906660e-24 g， V_cell = (a[Å] × 1e-8 cm/Å)^3。

返回:: 理论密度（g/cm³）。
返回类型:: float

备注

仅在结构为 "fcc" 或 "diamond" 时使用该公式；其他结构返回 NotImplementedError。
结果为近似理论值，忽略温度导致的体膨胀与缺陷等因素。

property young_modulus: float

计算杨氏模量

使用关系: E = 9*K*G / (3*K + G) 其中K为体积模量，G为剪切模量

返回:: 杨氏模量 (GPa)
返回类型:: float

name: str

symbol: str

mass_amu: float

lattice_constant: float

structure: str

literature_elastic_constants: dict[str, float]

thermoelasticsim.elastic.get_material_by_symbol(symbol)[源代码]

根据元素符号获取材料参数

参数:: symbol (str) -- 化学元素符号，如"Al", "Cu"
返回:: 匹配的材料参数，未找到返回None
返回类型:: MaterialParameters or None

示例

>>> al_params = get_material_by_symbol("Al")
>>> if al_params:
...     print(f"找到材料: {al_params.name}")
找到材料: Aluminum

thermoelasticsim.elastic.get_all_materials()[源代码]

获取所有预定义材料参数

返回:: 键为材料名称，值为MaterialParameters对象的字典
返回类型:: dict

示例

>>> materials = get_all_materials()
>>> for name, params in materials.items():
...     print(f"{name}: {params.structure.upper()}")
Aluminum: FCC
Copper: FCC
Gold: FCC

thermoelasticsim.elastic.compare_elastic_constants(material1, material2)[源代码]

比较两种材料的弹性常数

参数:

material1 (MaterialParameters) -- 要比较的材料参数
material2 (MaterialParameters) -- 要比较的材料参数

返回:

包含比较结果的嵌套字典

返回类型:

示例

>>> al = ALUMINUM_FCC
>>> cu = COPPER_FCC
>>> comparison = compare_elastic_constants(al, cu)
>>> print(f"C11比值: {comparison['ratios']['C11']:.2f}")
C11比值: 0.65

class thermoelasticsim.elastic.BenchmarkConfig(supercell_size=(3, 3, 3), delta=0.003, num_steps=1, shear_strains=(-0.004, -0.003, -0.002, -0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015, 0.002, 0.003, 0.004), optimizer_type='L-BFGS', optimizer_params=None, precision_mode=False, small_linear_strains=(-2e-05, -1e-05, 0.0, 1e-05, 2e-05), base_stress_tol_gpa=1e-05, base_relax_max_passes=3, small_shear_strains=(-0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015))[源代码]

基类：object

基准参数配置

参数:

supercell_size (tuple[int, int, int])
delta (float)
num_steps (int)
shear_strains (tuple[float, ...])
optimizer_type (str)
optimizer_params (dict[str, Any])
precision_mode (bool)
small_linear_strains (tuple[float, ...])
base_stress_tol_gpa (float)
base_relax_max_passes (int)
small_shear_strains (tuple[float, ...])

__init__(supercell_size=(3, 3, 3), delta=0.003, num_steps=1, shear_strains=(-0.004, -0.003, -0.002, -0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015, 0.002, 0.003, 0.004), optimizer_type='L-BFGS', optimizer_params=None, precision_mode=False, small_linear_strains=(-2e-05, -1e-05, 0.0, 1e-05, 2e-05), base_stress_tol_gpa=1e-05, base_relax_max_passes=3, small_shear_strains=(-0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015))

参数:

supercell_size (tuple[int, int, int])
delta (float)
num_steps (int)
shear_strains (tuple[float, ...])
optimizer_type (str)
optimizer_params (dict[str, Any])
precision_mode (bool)
small_linear_strains (tuple[float, ...])
base_stress_tol_gpa (float)
base_relax_max_passes (int)
small_shear_strains (tuple[float, ...])

返回类型:

None

base_relax_max_passes: int = 3

base_stress_tol_gpa: float = 1e-05

build_relaxer()[源代码]

构建结构弛豫器

返回类型:: StructureRelaxer

delta: float = 0.003

num_steps: int = 1

optimizer_params: dict[str, Any] = None

optimizer_type: str = 'L-BFGS'

precision_mode: bool = False

shear_strains: tuple[float, ...] = (-0.004, -0.003, -0.002, -0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015, 0.002, 0.003, 0.004)

small_linear_strains: tuple[float, ...] = (-2e-05, -1e-05, 0.0, 1e-05, 2e-05)

small_shear_strains: tuple[float, ...] = (-0.0015, -0.001, -0.0005, 0.0, 0.0005, 0.001, 0.0015)

supercell_size: tuple[int, int, int] = (3, 3, 3)

thermoelasticsim.elastic.calculate_c11_c12_traditional(cell, potential, relaxer, material_params, strain_points=None, do_internal_relax=True)[源代码]

使用单轴应变法计算 C11/C12，测试点与示例保持一致。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (list[float] | ndarray | None)
do_internal_relax (bool)

返回类型:

thermoelasticsim.elastic.calculate_c11_c12_robust(cell, potential, relaxer, material_params, strain_points=None, do_internal_relax=True)[源代码]

稳健法计算 C11/C12：正交应变(δ,-δ,0) + 等压应变(η,η,η)。

正交：σxx = (C11 - C12)·δ, σyy = -(C11 - C12)·δ
等压：σxx = (C11 + 2C12)·η = 3K·η

通过两次一维拟合解出 C11、C12，通常比单轴法对势细节更鲁棒。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (list[float] | ndarray | None)
do_internal_relax (bool)

返回类型:

thermoelasticsim.elastic.calculate_c44_lammps_shear(cell, potential, relaxer, material_params, strain_points)[源代码]

使用LAMMPS风格剪切方法计算C44/C55/C66。

参数:

relaxer (StructureRelaxer)
material_params (MaterialParameters)
strain_points (ndarray | list[float])

返回类型:

thermoelasticsim.elastic.run_zero_temp_benchmark(material_params, potential, supercell_size=(3, 3, 3), output_dir=None, save_json=True, precision=False, log_level=None, *, optimizer_type=None, uniaxial_strains=None, shear_strains=None, optimizer_params=None)[源代码]

运行通用零温弹性基准（结构自适应）。

根据 material_params.structure 自动生成晶体结构（当前支持 fcc、diamond）。
其它流程保持与铝基准一致：弛豫、C11/C12 单轴、C44 剪切、绘图与CSV/JSON。

参数:

material_params (MaterialParameters)
supercell_size (tuple[int, int, int])
output_dir (str | None)
save_json (bool)
precision (bool)
log_level (int | None)
optimizer_type (str | None)
uniaxial_strains (list[float] | tuple[float, ...] | None)
shear_strains (list[float] | tuple[float, ...] | None)
optimizer_params (dict[str, Any] | None)

返回类型:

thermoelasticsim.elastic.run_size_sweep(sizes=None, output_root=None, material_params=None, potential_factory=None, precision=False)[源代码]

依次运行 2x2x2、3x3x3、4x4x4 超胞的有限尺寸效应基准。

默认材料为铝（ALUMINUM_FCC），默认势为 EAMAl1Potential。
可通过传入 material_params 与 potential_factory 进行泛化。

参数:

sizes (list[tuple[int, int, int]] | None)
output_root (str | None)
material_params (MaterialParameters | None)
potential_factory (Callable[[], Any] | None)
precision (bool)

返回类型:

list[dict]

class thermoelasticsim.elastic.ElasticWaveAnalyzer(C11, C12, C44, density)[源代码]

基类：object

基于弹性常数的波速解析计算器（立方晶系）。

参数:

C11 (float) -- 立方晶系弹性常数，单位 GPa。
C12 (float) -- 立方晶系弹性常数，单位 GPa。
C44 (float) -- 立方晶系弹性常数，单位 GPa。
density (float) -- 材料密度，单位 g/cm^3。

备注

使用Christoffel方程 Γ·u = v^2 u 求解本征值 v^2 与本征向量 u。
在当前单位系统下，Γ 以 (km/s)^2 表示，故速度可直接取平方根。
纵波的判定依据：偏振方向与传播方向夹角最小（|n·u| 最大）。

C11: float

C12: float

C44: float

density: float

__init__(C11, C12, C44, density)

参数:

C11 (float)
C12 (float)
C44 (float)
density (float)

返回类型:

None

calculate_wave_velocities(direction)[源代码]

计算指定方向的纵横波速度与偏振方向。

参数:

direction (Iterable[float]) -- 传播方向向量或米勒指数，例如 [1,0,0], [1,1,0], [1,1,1]。

返回:

包含下列键：

'longitudinal': float，纵波速度 (km/s)
'transverse1': float，横波速度1 (km/s)
'transverse2': float，横波速度2 (km/s)
'polarizations': dict，包含三个模式的单位偏振向量

返回类型:

generate_report()[源代码]

生成标准方向 [100]、[110]、[111] 的波速报告。

返回:: 形如 {"[100]": {...}, "[110]": {...}, "[111]": {...}} 的字典。
返回类型:: dict

class thermoelasticsim.elastic.WaveExcitation(direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, phase=0.0, mode='standing', phase_speed_km_s=None)[源代码]

基类：object

平面波激发器（教学版）

参数:

direction (Axis) -- 传播方向，当前仅支持 "x"（等价于[100]）。
polarization ({"L", "Ty", "Tz"}) -- 极化：纵波 L（平行 x），或横波 Ty/Tz（沿 y/z）。
n_waves (int) -- 超胞长度内包含的整波数（满足 PBC 的允许波矢）。
amplitude_velocity (float) -- 初始速度幅值（Å/fs），建议 1e-4–1e-3。
amplitude_displacement (float | None) -- 初始位移幅值（Å）；若提供则叠加位移型激发（可为空）。
phase (float) -- 初始相位（弧度）。
mode (Literal['standing', 'traveling'])
phase_speed_km_s (float | None)

备注

平面波形式： s(r) = A cos(k·r + φ)
k = 2π n / Lx · ex（n 为整数，满足 PBC）
纵波：极化 e = ex；横波：e = ey 或 ez。
施加后移除质心平动，避免总动量漂移。

direction: Literal['x', 'y', 'z']

polarization: Literal['L', 'Ty', 'Tz']

n_waves: int

amplitude_velocity: float

amplitude_displacement: float | None

phase: float

mode: Literal['standing', 'traveling']

phase_speed_km_s: float | None

__init__(direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, phase=0.0, mode='standing', phase_speed_km_s=None)

参数:

direction (Literal['x', 'y', 'z'])
polarization (Literal['L', 'Ty', 'Tz'])
n_waves (int)
amplitude_velocity (float)
amplitude_displacement (float | None)
phase (float)
mode (Literal['standing', 'traveling'])
phase_speed_km_s (float | None)

返回类型:

None

apply(cell)[源代码]

对晶胞施加平面波激发（原地修改 cell）。

参数:: cell (thermoelasticsim.core.structure.Cell)
返回类型:: None

class thermoelasticsim.elastic.DynamicsConfig(supercell=(64, 12, 12), dt_fs=0.5, steps=6000, sample_every=50, direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, use_source=True, source_slab_fraction=0.06, source_amplitude_velocity=0.0005, source_t0_fs=200.0, source_sigma_fs=80.0, source_type='gaussian', source_cycles=4, source_freq_thz=1.0, record_trajectory=False, trajectory_file=None, detector_frac_a=0.2, detector_frac_b=0.7, measure_method='auto', v_max_km_s=None, absorber_enabled=False, absorber_slab_fraction=0.1, absorber_tau_fs=250.0, absorber_profile='cos2')[源代码]

基类：object

MD 波传播最小配置。

参数:

supercell (tuple[int, int, int])
dt_fs (float)
steps (int)
sample_every (int)
direction (Literal['x', 'y', 'z'])
polarization (Literal['L', 'Ty', 'Tz'])
n_waves (int)
amplitude_velocity (float)
amplitude_displacement (float | None)
use_source (bool)
source_slab_fraction (float)
source_amplitude_velocity (float)
source_t0_fs (float)
source_sigma_fs (float)
source_type (Literal['gaussian', 'tone_burst'])
source_cycles (int)
source_freq_thz (float)
record_trajectory (bool)
trajectory_file (str | None)
detector_frac_a (float)
detector_frac_b (float)
measure_method (str)
v_max_km_s (float | None)
absorber_enabled (bool)
absorber_slab_fraction (float)
absorber_tau_fs (float)
absorber_profile (Literal['cos2', 'linear'])

supercell: tuple[int, int, int]

dt_fs: float

steps: int

sample_every: int

direction: Literal['x', 'y', 'z']

polarization: Literal['L', 'Ty', 'Tz']

n_waves: int

amplitude_velocity: float

amplitude_displacement: float | None

use_source: bool

source_slab_fraction: float

source_amplitude_velocity: float

source_t0_fs: float

source_sigma_fs: float

source_type: Literal['gaussian', 'tone_burst']

source_cycles: int

source_freq_thz: float

record_trajectory: bool

trajectory_file: str | None

detector_frac_a: float

detector_frac_b: float

measure_method: str

v_max_km_s: float | None

absorber_enabled: bool

absorber_slab_fraction: float

absorber_tau_fs: float

absorber_profile: Literal['cos2', 'linear']

__init__(supercell=(64, 12, 12), dt_fs=0.5, steps=6000, sample_every=50, direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, use_source=True, source_slab_fraction=0.06, source_amplitude_velocity=0.0005, source_t0_fs=200.0, source_sigma_fs=80.0, source_type='gaussian', source_cycles=4, source_freq_thz=1.0, record_trajectory=False, trajectory_file=None, detector_frac_a=0.2, detector_frac_b=0.7, measure_method='auto', v_max_km_s=None, absorber_enabled=False, absorber_slab_fraction=0.1, absorber_tau_fs=250.0, absorber_profile='cos2')

参数:

supercell (tuple[int, int, int])
dt_fs (float)
steps (int)
sample_every (int)
direction (Literal['x', 'y', 'z'])
polarization (Literal['L', 'Ty', 'Tz'])
n_waves (int)
amplitude_velocity (float)
amplitude_displacement (float | None)
use_source (bool)
source_slab_fraction (float)
source_amplitude_velocity (float)
source_t0_fs (float)
source_sigma_fs (float)
source_type (Literal['gaussian', 'tone_burst'])
source_cycles (int)
source_freq_thz (float)
record_trajectory (bool)
trajectory_file (str | None)
detector_frac_a (float)
detector_frac_b (float)
measure_method (str)
v_max_km_s (float | None)
absorber_enabled (bool)
absorber_slab_fraction (float)
absorber_tau_fs (float)
absorber_profile (Literal['cos2', 'linear'])

返回类型:

None

thermoelasticsim.elastic.simulate_plane_wave_mvp(material_symbol='Al', dynamics=None, excitation=None, out_xt_path=None)[源代码]

运行最小版 MD 平面波传播并测速（返回结果字典）。

为教学简化，当前固定使用 EAM Al1 势与 FCC 结构。后续可抽象至 CLI。

参数:

material_symbol (str)
dynamics (DynamicsConfig | None)
excitation (WaveExcitation | None)
out_xt_path (str | None)

返回类型:

thermoelasticsim.elastic.plot_polar_plane(analyzer, plane='001', n_angles=360, outpath=None, dpi=300, annotate_hkls=None)[源代码]

绘制给定晶面内的声速各向异性极图。

参数:

analyzer (ElasticWaveAnalyzer) -- 解析计算器实例。
plane (str, optional) -- 平面标识（"001"/"110"/"111"），默认"001"。
n_angles (int, optional) -- 角度采样点数，默认360。
outpath (str, optional) -- 输出文件路径；若为None则使用当前目录下的 'anisotropy_polar.png'。
dpi (int, optional) -- 保存DPI，默认300。
annotate_hkls (Sequence[Sequence[int]] | None)

返回:

生成的图像路径。

返回类型:

thermoelasticsim.elastic.plot_velocity_surface_3d(analyzer, mode='L', n_theta=60, n_phi=120, output_html=None, output_png=None)[源代码]

绘制三维速度面（plotly，可交互）。

参数:

analyzer (ElasticWaveAnalyzer) -- 解析计算器实例。
mode ({"L", "Tmin", "Tmax"}) -- 绘制纵波面或两支横波中的较小/较大。
n_theta (int) -- 角度网格分辨率（θ: [0,π], φ: [0,2π)）。
n_phi (int) -- 角度网格分辨率（θ: [0,π], φ: [0,2π)）。
output_html (str, optional) -- 交互版HTML文件路径（推荐）。
output_png (str, optional) -- 静态PNG文件路径（需要安装kaleido）。

返回:

实际生成的文件路径；若依赖缺失则返回 (None, None)。

返回类型:

(html_path, png_path)

力学计算

class thermoelasticsim.elastic.mechanics.TensorConverter[源代码]

基类：object

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

static from_voigt(voigt, tensor_type)[源代码]

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

参数:

voigt (np.ndarray) -- 形状为 (6,) 的 Voigt 向量。
tensor_type (str) -- 张量类型，必须是 'stress' 或 'strain'。

返回:

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

返回类型:

np.ndarray

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

class thermoelasticsim.elastic.mechanics.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_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/Å³
返回类型:: numpy.ndarray

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\) 的力

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

维里应力张量 (3, 3)，单位 eV/Å³

返回类型:

np.ndarray

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

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

总应力张量 (3, 3)，单位 eV/Å³

返回类型:

np.ndarray

calculate_finite_difference_stress(cell, potential, dr=1e-06)[源代码]

已移除：有限差分应力路径不再提供，避免误用。

返回类型:: ndarray

get_all_stress_components(cell, potential)[源代码]

计算应力张量的所有分量

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

包含应力张量分量的字典，键包括： - "kinetic": 动能应力张量 - "virial": 维里应力张量 - "total": 总应力张量（kinetic + virial）

返回类型:

Dict[str, np.ndarray]

compute_stress(cell, potential)[源代码]

计算应力张量（使用总应力作为主要方法）

参数:

cell (Cell) -- 包含原子的晶胞对象
potential (Potential) -- 势能对象

返回:

3x3 应力张量矩阵（总应力 = 动能 + 维里）

返回类型:

np.ndarray

calculate_virial_kinetic_stress(cell, potential)[源代码]

向后兼容别名 - 指向总应力方法

返回类型:: ndarray

calculate_stress_basic(cell, potential)[源代码]

向后兼容别名 - 指向总应力方法

返回类型:: ndarray

validate_tensor_symmetry(tensor, tolerance=1e-10)[源代码]

验证应力张量是否对称

参数:

tensor (np.ndarray) -- 应力张量 (3x3)
tolerance (float, optional) -- 对称性的容差, 默认值为1e-10

返回:

如果对称则为True, 否则为False

返回类型:

class thermoelasticsim.elastic.mechanics.StrainCalculator[源代码]

基类：object

应变计算器类

参数:: F (numpy.ndarray) -- 3x3 变形矩阵

compute_strain(F)[源代码]

计算应变张量并返回 Voigt 表示法

参数:: F (numpy.ndarray) -- 3x3 变形矩阵
返回:: 应变向量，形状为 (6,)
返回类型:: numpy.ndarray

涨落法（开发中）

涨落法子包

弹性波传播模拟

解析计算

弹性波解析计算模块

本模块实现立方晶系材料中，基于弹性常数 C11, C12, C44 与材料密度 ρ 的弹性波速解析计算（阶段A）。支持任意传播方向的Christoffel方程求解，自动区分纵波与横波，并给出相应偏振向量。

单位与约定

弹性常数输入单位：GPa
密度输入单位：g/cm^3
输出速度单位：km/s

在上述单位下，波速满足 v[km/s] = sqrt(C[GPa] / ρ[g/cm^3])。

示例

>>> ana = ElasticWaveAnalyzer(C11=110.0, C12=61.0, C44=33.0, density=2.70)
>>> ana.calculate_wave_velocities([1, 0, 0])["longitudinal"]
6.38...
>>> report = ana.generate_report()
>>> sorted(report.keys())
['[100]', '[110]', '[111]']

class thermoelasticsim.elastic.wave.analytical.Iterable: 基类：object

thermoelasticsim.elastic.wave.analytical.dataclass(cls=None, /, *, init=True, repr=True, eq=True, order=False, unsafe_hash=False, frozen=False, match_args=True, kw_only=False, slots=False, weakref_slot=False)[源代码]

Add dunder methods based on the fields defined in the class.

Examines PEP 526 __annotations__ to determine fields.

If init is true, an __init__() method is added to the class. If repr is true, a __repr__() method is added. If order is true, rich comparison dunder methods are added. If unsafe_hash is true, a __hash__() method is added. If frozen is true, fields may not be assigned to after instance creation. If match_args is true, the __match_args__ tuple is added. If kw_only is true, then by default all fields are keyword-only. If slots is true, a new class with a __slots__ attribute is returned.

class thermoelasticsim.elastic.wave.analytical.ElasticWaveAnalyzer(C11, C12, C44, density)[源代码]

基类：object

基于弹性常数的波速解析计算器（立方晶系）。

参数:

C11 (float) -- 立方晶系弹性常数，单位 GPa。
C12 (float) -- 立方晶系弹性常数，单位 GPa。
C44 (float) -- 立方晶系弹性常数，单位 GPa。
density (float) -- 材料密度，单位 g/cm^3。

备注

使用Christoffel方程 Γ·u = v^2 u 求解本征值 v^2 与本征向量 u。
在当前单位系统下，Γ 以 (km/s)^2 表示，故速度可直接取平方根。
纵波的判定依据：偏振方向与传播方向夹角最小（|n·u| 最大）。

C11: float

C12: float

C44: float

density: float

calculate_wave_velocities(direction)[源代码]

计算指定方向的纵横波速度与偏振方向。

参数:

direction (Iterable[float]) -- 传播方向向量或米勒指数，例如 [1,0,0], [1,1,0], [1,1,1]。

返回:

包含下列键：

'longitudinal': float，纵波速度 (km/s)
'transverse1': float，横波速度1 (km/s)
'transverse2': float，横波速度2 (km/s)
'polarizations': dict，包含三个模式的单位偏振向量

返回类型:

generate_report()[源代码]

生成标准方向 [100]、[110]、[111] 的波速报告。

返回:: 形如 {"[100]": {...}, "[110]": {...}, "[111]": {...}} 的字典。
返回类型:: dict

__init__(C11, C12, C44, density)

参数:

C11 (float)
C12 (float)
C44 (float)
density (float)

返回类型:

None

动力学模拟

弹性波动力学模拟模块

本模块提供弹性波在晶格中传播的分子动力学模拟功能。

主要功能

平面波激发：支持位移型和速度型激发，以及时间域源注入
NVE传播模拟：零温条件下的短时间线性响应模拟
速度测量：多种算法（互相关、到达时间拟合、k-ω谱分析）
可视化输出：时空图（x-t图）显示位移场和包络

物理模型

传播方向：当前仅支持[100]方向（x轴）
边界条件：y/z方向周期性边界，x方向准无限介质
激发模式：小振幅线性响应区，避免非线性效应
吸收边界：可选的海绵层减少边界反射

单位约定

长度：Å（埃）
时间：fs（飞秒）
速度：Å/fs（内部），km/s（输出，换算关系：1 Å/fs = 100 km/s）

备注

本模块为教学演示设计，追求物理清晰性和代码可理解性。后续可扩展支持更多传播方向和材料类型。

class thermoelasticsim.elastic.wave.dynamics.WaveExcitation(direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, phase=0.0, mode='standing', phase_speed_km_s=None)[源代码]

基类：object

平面波激发器（教学版）

参数:

direction (Axis) -- 传播方向，当前仅支持 "x"（等价于[100]）。
polarization ({"L", "Ty", "Tz"}) -- 极化：纵波 L（平行 x），或横波 Ty/Tz（沿 y/z）。
n_waves (int) -- 超胞长度内包含的整波数（满足 PBC 的允许波矢）。
amplitude_velocity (float) -- 初始速度幅值（Å/fs），建议 1e-4–1e-3。
amplitude_displacement (float | None) -- 初始位移幅值（Å）；若提供则叠加位移型激发（可为空）。
phase (float) -- 初始相位（弧度）。
mode (Literal['standing', 'traveling'])
phase_speed_km_s (float | None)

备注

平面波形式： s(r) = A cos(k·r + φ)
k = 2π n / Lx · ex（n 为整数，满足 PBC）
纵波：极化 e = ex；横波：e = ey 或 ez。
施加后移除质心平动，避免总动量漂移。

direction: Literal['x', 'y', 'z']

polarization: Literal['L', 'Ty', 'Tz']

n_waves: int

amplitude_velocity: float

amplitude_displacement: float | None

phase: float

mode: Literal['standing', 'traveling']

phase_speed_km_s: float | None

apply(cell)[源代码]

对晶胞施加平面波激发（原地修改 cell）。

参数:: cell (thermoelasticsim.core.structure.Cell)
返回类型:: None

__init__(direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, phase=0.0, mode='standing', phase_speed_km_s=None)

参数:

direction (Literal['x', 'y', 'z'])
polarization (Literal['L', 'Ty', 'Tz'])
n_waves (int)
amplitude_velocity (float)
amplitude_displacement (float | None)
phase (float)
mode (Literal['standing', 'traveling'])
phase_speed_km_s (float | None)

返回类型:

None

class thermoelasticsim.elastic.wave.dynamics.DynamicsConfig(supercell=(64, 12, 12), dt_fs=0.5, steps=6000, sample_every=50, direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, use_source=True, source_slab_fraction=0.06, source_amplitude_velocity=0.0005, source_t0_fs=200.0, source_sigma_fs=80.0, source_type='gaussian', source_cycles=4, source_freq_thz=1.0, record_trajectory=False, trajectory_file=None, detector_frac_a=0.2, detector_frac_b=0.7, measure_method='auto', v_max_km_s=None, absorber_enabled=False, absorber_slab_fraction=0.1, absorber_tau_fs=250.0, absorber_profile='cos2')[源代码]

基类：object

MD 波传播最小配置。

参数:

supercell (tuple[int, int, int])
dt_fs (float)
steps (int)
sample_every (int)
direction (Literal['x', 'y', 'z'])
polarization (Literal['L', 'Ty', 'Tz'])
n_waves (int)
amplitude_velocity (float)
amplitude_displacement (float | None)
use_source (bool)
source_slab_fraction (float)
source_amplitude_velocity (float)
source_t0_fs (float)
source_sigma_fs (float)
source_type (Literal['gaussian', 'tone_burst'])
source_cycles (int)
source_freq_thz (float)
record_trajectory (bool)
trajectory_file (str | None)
detector_frac_a (float)
detector_frac_b (float)
measure_method (str)
v_max_km_s (float | None)
absorber_enabled (bool)
absorber_slab_fraction (float)
absorber_tau_fs (float)
absorber_profile (Literal['cos2', 'linear'])

supercell: tuple[int, int, int]

dt_fs: float

steps: int

sample_every: int

direction: Literal['x', 'y', 'z']

polarization: Literal['L', 'Ty', 'Tz']

n_waves: int

amplitude_velocity: float

amplitude_displacement: float | None

use_source: bool

source_slab_fraction: float

source_amplitude_velocity: float

source_t0_fs: float

source_sigma_fs: float

source_type: Literal['gaussian', 'tone_burst']

source_cycles: int

source_freq_thz: float

record_trajectory: bool

trajectory_file: str | None

detector_frac_a: float

detector_frac_b: float

measure_method: str

v_max_km_s: float | None

absorber_enabled: bool

absorber_slab_fraction: float

absorber_tau_fs: float

absorber_profile: Literal['cos2', 'linear']

__init__(supercell=(64, 12, 12), dt_fs=0.5, steps=6000, sample_every=50, direction='x', polarization='L', n_waves=2, amplitude_velocity=0.0001, amplitude_displacement=None, use_source=True, source_slab_fraction=0.06, source_amplitude_velocity=0.0005, source_t0_fs=200.0, source_sigma_fs=80.0, source_type='gaussian', source_cycles=4, source_freq_thz=1.0, record_trajectory=False, trajectory_file=None, detector_frac_a=0.2, detector_frac_b=0.7, measure_method='auto', v_max_km_s=None, absorber_enabled=False, absorber_slab_fraction=0.1, absorber_tau_fs=250.0, absorber_profile='cos2')

参数:

supercell (tuple[int, int, int])
dt_fs (float)
steps (int)
sample_every (int)
direction (Literal['x', 'y', 'z'])
polarization (Literal['L', 'Ty', 'Tz'])
n_waves (int)
amplitude_velocity (float)
amplitude_displacement (float | None)
use_source (bool)
source_slab_fraction (float)
source_amplitude_velocity (float)
source_t0_fs (float)
source_sigma_fs (float)
source_type (Literal['gaussian', 'tone_burst'])
source_cycles (int)
source_freq_thz (float)
record_trajectory (bool)
trajectory_file (str | None)
detector_frac_a (float)
detector_frac_b (float)
measure_method (str)
v_max_km_s (float | None)
absorber_enabled (bool)
absorber_slab_fraction (float)
absorber_tau_fs (float)
absorber_profile (Literal['cos2', 'linear'])

返回类型:

None

thermoelasticsim.elastic.wave.dynamics.simulate_plane_wave_mvp(material_symbol='Al', dynamics=None, excitation=None, out_xt_path=None)[源代码]

运行最小版 MD 平面波传播并测速（返回结果字典）。

为教学简化，当前固定使用 EAM Al1 势与 FCC 结构。后续可抽象至 CLI。

参数:

material_symbol (str)
dynamics (DynamicsConfig | None)
excitation (WaveExcitation | None)
out_xt_path (str | None)

返回类型:

可视化工具

弹性波各向异性可视化

提供在指定晶面内采样传播方向，计算纵/横波速并绘制极坐标图的工具。

功能

sample_plane_directions：在给定面（"001"/"110"/"111"）上等角采样方向
compute_velocities_over_directions：批量计算 v_L, v_T1, v_T2
plot_polar_plane：生成各向异性极图（matplotlib）

说明

仅依赖 matplotlib，默认使用项目统一的中文字体与Agg后端。
角度 θ ∈ [0, 2π)，方向向量 n(θ) = u cosθ + v sinθ，其中 {u, v} 为该平面的正交基，法向量为平面的米勒指数方向。

thermoelasticsim.elastic.wave.visualization.sample_plane_directions(plane, n_angles=360)[源代码]

在指定晶面上等角采样传播方向。

参数:

plane (str) -- 平面标识（"001"/"110"/"111"）。
n_angles (int, optional) -- 角度采样数（默认360）。

返回:

(theta, directions)，其中 theta 形状为 (n_angles,)， directions 形状为 (n_angles, 3)，均为单位长度方向向量。

返回类型:

tuple[ndarray, ndarray]

thermoelasticsim.elastic.wave.visualization.compute_velocities_over_directions(analyzer, directions, plane=None)[源代码]

批量计算多个方向上的纵/横波速。

参数:

analyzer (ElasticWaveAnalyzer) -- 解析计算器实例。
directions (Iterable[Iterable[float]]) -- 多个单位方向向量。
plane (str | None)

返回:

包含 'vL', 'vT1', 'vT2' 三个键，对应等长的速度数组（km/s）。

返回类型:

thermoelasticsim.elastic.wave.visualization.plot_polar_plane(analyzer, plane='001', n_angles=360, outpath=None, dpi=300, annotate_hkls=None)[源代码]

绘制给定晶面内的声速各向异性极图。

参数:

analyzer (ElasticWaveAnalyzer) -- 解析计算器实例。
plane (str, optional) -- 平面标识（"001"/"110"/"111"），默认"001"。
n_angles (int, optional) -- 角度采样点数，默认360。
outpath (str, optional) -- 输出文件路径；若为None则使用当前目录下的 'anisotropy_polar.png'。
dpi (int, optional) -- 保存DPI，默认300。
annotate_hkls (Sequence[Sequence[int]] | None)

返回:

生成的图像路径。

返回类型: