弹性常数计算模块 (elastic)
本模块提供材料弹性常数的计算功能。
显式形变法
零温计算
零温显式形变法弹性常数计算模块
该模块实现零温条件下通过显式形变法计算弹性常数。采用三层架构： 底层结构弛豫 → 中层工作流管理 → 高层数据求解
基本原理： 1. 制备无应力基态 2. 施加微小形变 3. 优化原子位置 4. 测量应力响应 5. 线性拟合求解弹性常数
数学基础： 胡克定律的张量形式：\(\\sigma = C : \\varepsilon\)
在Voigt记号下：\(\\sigma = C \\cdot \\varepsilon\)
其中： - \(\\sigma\)：应力向量 (6×1) - \(C\)：弹性常数矩阵 (6×6) - \(\\varepsilon\)：应变向量 (6×1)
- class thermoelasticsim.elastic.deformation_method.zero_temp.DeformationResult(strain_voigt, stress_voigt, converged, deformation_matrix)[源代码]
基类：
object单次形变计算结果的数据容器
- strain_voigt
Voigt记号应变向量，形状 (6,)
- Type:
- stress_voigt
Voigt记号应力向量，形状 (6,)
- Type:
- deformation_matrix
形变矩阵F，形状 (3,3)
- Type:
- class thermoelasticsim.elastic.deformation_method.zero_temp.StructureRelaxer(optimizer_type='L-BFGS', optimizer_params=None, supercell_dims=None, trajectory_recorder=None)[源代码]
基类：
object结构弛豫计算引擎
负责执行结构优化计算，包括完全弛豫和内部弛豫两种模式。
- 参数:
备注
两种弛豫模式的区别：
完全弛豫：同时优化原子位置和晶胞参数，用于制备无应力基态
内部弛豫：仅优化原子位置，保持晶胞形状固定，用于形变后的平衡
示例
>>> relaxer = StructureRelaxer() >>> # 完全弛豫制备基态 >>> converged = relaxer.full_relax(cell, potential) >>> # 形变后内部弛豫 >>> converged = relaxer.internal_relax(deformed_cell, potential)
Methods
full_relax(cell, potential)执行完全弛豫：同时优化原子位置和晶胞参数
internal_relax(cell, potential)执行内部弛豫：仅优化原子位置，保持晶胞形状固定
uniform_lattice_relax(cell, potential)等比例晶格弛豫：只优化晶格常数，保持原子相对位置不变
- __init__(optimizer_type='L-BFGS', optimizer_params=None, supercell_dims=None, trajectory_recorder=None)[源代码]
初始化结构弛豫器
- 参数:
- 抛出:
ValueError -- 如果指定了不支持的优化器类型
- full_relax(cell, potential)[源代码]
执行完全弛豫：同时优化原子位置和晶胞参数
完全弛豫用于制备无应力基态，是零温弹性常数计算的第一步。 通过同时优化原子位置和晶胞参数，消除系统内应力。
备注
数学表述：
\[\begin{split}\\min_{r,h} E(r,h) \\quad \\text{s.t.} \\quad \\sigma = 0\end{split}\]其中： - \(r\)：原子位置 - \(h\)：晶胞参数 - \(E\)：总能量 - \(\\sigma\)：应力张量
示例
>>> relaxer = StructureRelaxer() >>> converged = relaxer.full_relax(cell, potential) >>> if converged: ... print("成功制备无应力基态")
- internal_relax(cell, potential)[源代码]
执行内部弛豫：仅优化原子位置，保持晶胞形状固定
内部弛豫用于形变后的结构优化，在保持宏观形变的前提下， 寻找原子的最优位置配置。
备注
数学表述：
\[\begin{split}\\min_{r} E(r,h_{\\text{fixed}}) \\quad \\text{s.t.} \\quad h = \\text{const}\end{split}\]其中： - \(r\)：原子位置（优化变量） - \(h_{\\text{fixed}}\)：固定的晶胞参数 - \(E\)：总能量
示例
>>> relaxer = StructureRelaxer() >>> # 先施加形变 >>> cell.apply_deformation(deformation_matrix) >>> # 再内部弛豫 >>> converged = relaxer.internal_relax(cell, potential)
- uniform_lattice_relax(cell, potential)[源代码]
等比例晶格弛豫：只优化晶格常数，保持原子相对位置不变
这种弛豫方式特别适合基态制备，既避免了原子位置的数值噪声， 又保持了晶体结构的完美对称性，同时计算效率极高。
备注
数学表述：
\[\begin{split}\\min_{s} E(r_{\\text{scaled}}, s \\cdot L) \\quad \\text{其中} \\quad r_{\\text{scaled}} = s \\cdot r_{\\text{frac}} \\cdot L\end{split}\]其中： - \(s\)：等比例缩放因子（唯一优化变量） - \(r_{\\text{frac}}\)：固定的分数坐标 - \(L\)：原始晶格矢量矩阵 - \(E\)：总能量
优势： - 自由度仅为1，计算极快 - 保持晶体结构完美对称性 - 避免原子位置数值噪声 - 适合各种晶系（不限于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零温显式形变法计算器
管理从基态制备到弹性常数求解的完整计算流程。 实现标准的显式形变法：制备基态→施加形变→内部弛豫→测量应力→线性拟合。
- 参数:
- relaxer
结构弛豫器实例
- Type:
- reference_stress
基态参考应力张量
- Type:
备注
计算流程包含5个步骤：
基态制备：完全弛豫获得无应力基态
形变生成：生成6个Voigt分量的形变矩阵序列
应力计算：对每个形变施加内部弛豫并测量应力
数据收集：收集所有应力-应变数据对
线性拟合：通过最小二乘法求解弹性常数矩阵
应变生成策略： - 对称应变：\(\\varepsilon_{11}, \\varepsilon_{22}, \\varepsilon_{33}\) - 剪切应变：\(\\varepsilon_{23}, \\varepsilon_{13}, \\varepsilon_{12}\) - 应变范围：\([-\\delta, +\\delta]\)，均匀分布
示例
>>> calculator = ZeroTempDeformationCalculator(cell, potential, delta=0.005) >>> elastic_matrix, r2_score = calculator.calculate() >>> print(f"弹性常数矩阵 (GPa):\\n{elastic_matrix}") >>> print(f"拟合优度 R²: {r2_score:.6f}")
Methods
执行完整的零温弹性常数计算
- __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²的元组
- 返回类型:
备注
该方法执行5步计算流程：
制备无应力基态
生成形变矩阵序列
计算每个形变的应力响应
收集应力-应变数据
线性拟合求解弹性常数
拟合质量评估： - R² > 0.999：优秀 - R² > 0.99：良好 - R² < 0.99：需改进
示例
>>> calculator = ZeroTempDeformationCalculator(cell, potential) >>> C_matrix, r2 = calculator.calculate() >>> >>> # 检查拟合质量 >>> if r2 > 0.999: ... print("拟合质量优秀") >>> else: ... print(f"拟合质量R²={r2:.6f}，可能需要调整参数")
- class thermoelasticsim.elastic.deformation_method.zero_temp.ElasticConstantsSolver[源代码]
基类：
object弹性常数求解器
从应力应变数据通过线性回归求解弹性常数矩阵。 支持最小二乘法和岭回归两种方法，并提供拟合质量评估。
备注
数学原理：
根据胡克定律：\(\\sigma = C \\cdot \\varepsilon\)
其中： - \(\\sigma\)：应力向量 (N×6) - \(C\)：弹性常数矩阵 (6×6) - \(\\varepsilon\)：应变向量 (N×6)
最小二乘求解：
\[\begin{split}\\min_C ||\\sigma - C \\cdot \\varepsilon||^2\end{split}\]解析解：
\[\begin{split}C = (\\varepsilon^T \\varepsilon)^{-1} \\varepsilon^T \\sigma\end{split}\]示例
>>> solver = ElasticConstantsSolver() >>> C_matrix, r2_score = solver.solve(strains, stresses) >>> print(f"弹性常数矩阵:\\n{C_matrix}") >>> print(f"拟合优度: {r2_score:.6f}")
Methods
solve(strains, stresses[, method, alpha])求解弹性常数矩阵
- 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²
- 返回类型:
- 抛出:
ValueError -- 如果输入数据格式不正确或求解方法不支持
备注
支持的求解方法：
最小二乘法 ('least_squares')： 标准线性回归，适用于大多数情况
岭回归 ('ridge')： 加入L2正则化，适用于病态矩阵情况
拟合优度R²计算：
\[\begin{split}R^2 = 1 - \\frac{SS_{res}}{SS_{tot}}\end{split}\]其中： - \(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[源代码]
基类：
objectLAMMPS风格剪切形变方法
实现经过验证的LAMMPS风格剪切形变算法，用于计算剪切弹性常数(C44, C55, C66)。 该方法通过晶胞剪切和原子位置调整来施加纯剪切应变，避免了传统形变矩阵方法 在剪切计算中的数值问题。
备注
LAMMPS剪切形变的核心原理：
晶胞剪切：直接修改晶格矢量的非对角分量
原子位置调整：按照仿射变换调整所有原子位置
应变-应力关系：通过virial应力张量计算剪切应力
支持的剪切方向： - direction=4: yz剪切 (σ23, C44) - direction=5: xz剪切 (σ13, C55) - direction=6: xy剪切 (σ12, C66)
剪切应变定义：
\[\begin{split}\\gamma_{ij} = \\frac{\\partial u_i}{\\partial x_j} + \\frac{\\partial u_j}{\\partial x_i}\end{split}\]其中γ是工程剪切应变，与张量剪切应变的关系为 \(\\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 ... )
Methods
apply_shear_deformation(cell, direction, ...)对晶胞施加LAMMPS风格剪切形变
calculate_c44_response(cell, potential, ...)计算C44剪切响应
- apply_shear_deformation(cell, direction, strain_magnitude)[源代码]
对晶胞施加LAMMPS风格剪切形变
该方法实现了LAMMPS中的剪切形变算法，通过同时修改晶胞参数和原子位置 来施加纯剪切应变。这种方法保持了晶胞的周期性边界条件，并确保 剪切形变的一致性。
- 参数:
- 返回:
施加剪切形变后的新晶胞对象
- 返回类型:
- 抛出:
ValueError -- 如果指定了不支持的剪切方向
备注
LAMMPS剪切形变算法：
yz剪切 (direction=4)：
\[ \begin{align}\begin{aligned}\begin{split}h_{zy} \\leftarrow h_{zy} + \\gamma \\cdot h_{zz}\end{split}\\\begin{split}y_i \\leftarrow y_i + \\gamma \\cdot z_i\end{split}\end{aligned}\end{align} \]xz剪切 (direction=5)：
\[ \begin{align}\begin{aligned}\begin{split}h_{zx} \\leftarrow h_{zx} + \\gamma \\cdot h_{zz}\end{split}\\\begin{split}x_i \\leftarrow x_i + \\gamma \\cdot z_i\end{split}\end{aligned}\end{align} \]xy剪切 (direction=6)：
\[ \begin{align}\begin{aligned}\begin{split}h_{yx} \\leftarrow h_{yx} + \\gamma \\cdot h_{yy}\end{split}\\\begin{split}x_i \\leftarrow x_i + \\gamma \\cdot y_i\end{split}\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剪切响应
使用LAMMPS风格剪切方法计算三个剪切弹性常数(C44, C55, C66)的应力响应。 该方法对每个剪切方向施加一系列应变，通过内部弛豫消除原子间力， 然后测量相应的剪切应力。
- 参数:
cell (Cell) -- 基态晶胞对象
potential (Potential) -- 势能函数
strain_amplitudes (numpy.ndarray) -- 应变幅度数组，建议范围[-0.005, 0.005]
relaxer (StructureRelaxer) -- 结构弛豫器实例
include_base_state (bool, optional) -- 是否包含零应变基态点，默认True
- 返回:
包含以下键的计算结果： - 'directions': 剪切方向信息 - 'detailed_results': 每个方向的详细数据 - 'summary': 拟合结果汇总
- 返回类型:
备注
计算流程：
基态制备：确保从无应力基态开始
多方向扫描：对三个剪切方向分别计算
应变序列：对每个方向施加应变序列
内部弛豫：固定晶胞形状，优化原子位置
应力测量：计算剪切分量应力响应
线性拟合：通过最小二乘法求解弹性常数
立方晶系的剪切弹性常数：
\[\begin{split}C_{44} = C_{55} = C_{66} = \\frac{\\partial\\sigma_{ij}}{\\partial\\gamma_{ij}}\end{split}\]其中σ是剪切应力，γ是工程剪切应变。
示例
>>> shear_method = ShearDeformationMethod() >>> strain_points = np.linspace(-0.004, 0.004, 9) >>> >>> results = shear_method.calculate_c44_response( ... base_cell, potential, strain_points, relaxer ... ) >>> >>> # 查看结果 >>> print(f"C44 = {results['summary']['C44']:.2f} GPa") >>> print(f"拟合优度 R² = {results['summary']['r2_score']:.6f}")
- thermoelasticsim.elastic.deformation_method.zero_temp.calculate_zero_temp_elastic_constants(cell, potential, delta=0.005, num_steps=5, **kwargs)[源代码]
零温弹性常数计算的便捷函数
这是一个高级接口，封装了完整的零温弹性常数计算流程。 适合快速计算和脚本使用。
- 参数:
- 返回:
弹性常数矩阵(GPa)和拟合优度R²
- 返回类型:
示例
>>> 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 - 有限温弹性常数计算模块
该模块实现了有限温条件下，通过显式形变法计算弹性常数的工作流。
- 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) -- 结果保存路径，默认为当前时间戳命名的目录
Methods
计算弹性常数的主流程
run_npt_equilibration(cell)在NPT系综下进行平衡
run_nvt_sampling(cell)在NVT系综下采样应力
- __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)[源代码]
力学计算
- class thermoelasticsim.elastic.mechanics.StressCalculator[源代码]
基类：
object应力张量计算器
处理三个张量贡献： 1. 动能张量: p_i⊗p_i/m_i 2. 维里张量: r_i⊗f_i 3. 晶格应变张量: ∂U/∂h·h^T
Methods
calculate_finite_difference_stress(cell, ...)计算有限差分应力张量（基于能量导数）
calculate_kinetic_stress(cell)计算纯动能应力张量
calculate_lattice_stress(cell, potential[, dr])为保持向后兼容性的别名方法 - 指向有限差分应力方法
calculate_stress_basic(cell, potential)向后兼容别名 - 指向总应力方法
calculate_total_stress(cell, potential)计算总应力张量（动能+维里项）
calculate_virial_kinetic_stress(cell, potential)向后兼容别名 - 指向总应力方法
calculate_virial_stress(cell, potential)计算纯维里应力张量（原子间相互作用力项）
compute_stress(cell, potential)计算应力张量（使用总应力作为主要方法）
get_all_stress_components(cell, potential)计算应力张量的所有分量
validate_tensor_symmetry(tensor[, tolerance])验证应力张量是否对称
- calculate_kinetic_stress(cell)[源代码]
计算纯动能应力张量
使用公式：σ_αβ^kinetic = -1/V * Σᵢ m_i v_iα v_iβ 在零温下此项为零
- 参数:
cell (Cell) -- 包含原子的晶胞对象
- 返回:
动能应力张量 (3x3)
- 返回类型:
np.ndarray
- calculate_virial_stress(cell, potential)[源代码]
计算纯维里应力张量（原子间相互作用力项）
使用公式：σ_αβ^virial = -1/V * Σ(i<j) r_ij^α * F_ij^β 其中： - r_ij = r_i - r_j (从原子j指向原子i的向量) - F_ij是原子j对原子i的力
- calculate_total_stress(cell, potential)[源代码]
计算总应力张量（动能+维里项）
使用公式：σ_αβ = -1/V * (Σᵢ m_i v_iα v_iβ + Σ(i<j) r_ij^α * F_ij^β)
- class thermoelasticsim.elastic.mechanics.StrainCalculator[源代码]
基类：
object应变计算器类
- 参数:
F (numpy.ndarray) -- 3x3 变形矩阵
Methods
计算应变张量并返回 Voigt 表示法
- compute_strain(F)[源代码]
计算应变张量并返回 Voigt 表示法
- 参数:
F (numpy.ndarray) -- 3x3 变形矩阵
- 返回:
应变向量，形状为 (6,)
- 返回类型:
涨落法（开发中）
涨落法子包