diff --git a/docs/source/user_guide/benchmarks/iron_properties.rst b/docs/source/user_guide/benchmarks/iron_properties.rst new file mode 100644 index 000000000..4cf54be93 --- /dev/null +++ b/docs/source/user_guide/benchmarks/iron_properties.rst @@ -0,0 +1,158 @@ +=============== +Iron Properties +=============== + +Summary +------- + +This benchmark evaluates MLIP performance on a comprehensive set of BCC iron +properties relevant to plasticity and fracture. The benchmark is based on +`Zhang et al. (2023) `_, which assessed the +efficiency, accuracy, and transferability of machine learning potentials for +dislocations and cracks in iron. + +Seven groups of properties are computed and compared against DFT (PBE) reference +values: equation of state, elastic constants, the Bain path, vacancy formation +energy, surface energies, generalised stacking fault energies, and +traction-separation curves. + +Metrics +------- + +EOS properties +^^^^^^^^^^^^^^ + +1. Lattice parameter error (%) + +The equilibrium BCC lattice parameter :math:`a_0` is obtained by fitting a +third-order Birch-Murnaghan equation of state to 30 energy-volume points +sampled around 2.834 Å. The percentage error relative to the DFT reference +value (2.831 Å) is reported. + +2. Bulk modulus error (%) + +The bulk modulus :math:`B_0` is extracted from the same EOS fit. The +percentage error relative to the DFT reference value (178.0 GPa) is reported. + + +Elastic constants +^^^^^^^^^^^^^^^^^ + +3. :math:`C_{11}` error (%) + +The elastic constant :math:`C_{11}` is computed using a stress-strain approach +on a 4x4x4 BCC supercell. Small positive and negative strains +(:math:`\pm 10^{-5}`) are applied along each Voigt direction, and the elastic +constants are extracted from the resulting stress differences. The percentage +error relative to the DFT reference value (296.7 GPa) is reported. + +4. :math:`C_{12}` error (%) + +Same as (3), for the elastic constant :math:`C_{12}`. Reference: 151.4 GPa. + +5. :math:`C_{44}` error (%) + +Same as (3), for the elastic constant :math:`C_{44}`. Reference: 104.7 GPa. + + +Vacancy formation energy +^^^^^^^^^^^^^^^^^^^^^^^^^ + +6. :math:`E_{\mathrm{vac}}` error (%) + +A single vacancy is created in a 4x4x4 BCC supercell by removing one atom. +Atomic positions are relaxed at fixed cell volume. The vacancy formation energy +is calculated as +:math:`E_{\mathrm{vac}} = E_{\mathrm{defect}} - E_{\mathrm{perfect}} + E_{\mathrm{coh}}`, +where :math:`E_{\mathrm{coh}}` is the cohesive energy per atom. The percentage +error relative to the DFT reference value (2.02 eV) is reported. + + +Bain path +^^^^^^^^^ + +7. BCC-FCC energy difference error (meV) + +The Bain path maps the continuous tetragonal distortion from BCC +(:math:`c/a = 1`) to FCC (:math:`c/a = \sqrt{2}`). For each of 65 target +:math:`c/a` ratios between 0.72 and 2.0, a tetragonally distorted cell is +created and its volume is relaxed isotropically (uniform scaling preserving the +:math:`c/a` ratio). The absolute error in the BCC-FCC energy difference +relative to the DFT reference (83.5 meV/atom) is reported. + + +Surface energies +^^^^^^^^^^^^^^^^ + +8. Surface energy MAE (J/m²) + +Surface energies are computed for the (100), (110), (111), and (112) cleavage +planes. For each surface, a slab is created with vacuum and the surface energy +is calculated as +:math:`\gamma = (E_{\mathrm{slab}} - E_{\mathrm{bulk}}) / 2A`. +Atomic positions are relaxed at fixed cell shape. The mean absolute error +across all four surfaces, relative to DFT reference values +(:math:`\gamma_{100}` = 2.41, :math:`\gamma_{110}` = 2.37, +:math:`\gamma_{111}` = 2.58, :math:`\gamma_{112}` = 2.48 J/m²), is reported. + + +Stacking fault energies +^^^^^^^^^^^^^^^^^^^^^^^ + +9. Max SFE :math:`\{110\}\langle111\rangle` error (%) + +The generalised stacking fault energy (GSFE) curve for the +:math:`\{110\}\langle111\rangle` slip system is computed by incrementally +displacing the upper half of a crystallographically oriented supercell along +the slip direction. The displacement covers one full Burgers vector +(:math:`b = a\sqrt{3}/2`) in 16 steps. Atoms are constrained to relax only +perpendicular to the fault plane. The percentage error in the maximum +(unstable) SFE relative to the DFT reference (0.75 J/m²) is reported. + +10. Max SFE :math:`\{112\}\langle111\rangle` error (%) + +Same as (9), for the :math:`\{112\}\langle111\rangle` slip system. +Reference: 1.12 J/m². + + +Traction-separation curves +^^^^^^^^^^^^^^^^^^^^^^^^^^ + +11. Max traction (100) error (%) + +A traction-separation curve is computed for the (100) cleavage plane by +incrementally separating crystal halves in 0.05 Å steps up to 5.0 Å, without +atomic relaxation, and measuring forces at each step. The traction (tensile +stress) is obtained from the sum of z-forces on the upper region divided by +the cross-sectional area. The percentage error in the maximum traction relative +to the DFT reference (35.0 GPa) is reported. + +12. Max traction (110) error (%) + +Same as (11), for the (110) cleavage plane. Reference: 30.0 GPa. + + +Computational cost +------------------ + +Medium: tests are likely to take minutes to run on GPU, or hours on CPU for each model. +The benchmark is marked as slow and excluded from default test runs. + + +Data availability +----------------- + +Input structures: + +* All structures are generated programmatically using ASE. BCC iron unit cells and + supercells are constructed from the equilibrium lattice parameter obtained via EOS + fitting. + +Reference data: + +* DFT (PBE) reference values from: + + * Zhang, L., Csányi, G., van der Giessen, E., & Maresca, F. (2023). + "Efficiency, Accuracy, and Transferability of Machine Learning Potentials: + Application to Dislocations and Cracks in Iron." + `arXiv:2307.10072 `_ diff --git a/docs/source/user_guide/benchmarks/physicality.rst b/docs/source/user_guide/benchmarks/physicality.rst index 0d5b55d25..7d392a6e7 100644 --- a/docs/source/user_guide/benchmarks/physicality.rst +++ b/docs/source/user_guide/benchmarks/physicality.rst @@ -135,3 +135,8 @@ Data availability ----------------- None required; diatomics are generated in ASE. + +.. toctree:: + :maxdepth: 2 + + iron_properties diff --git a/ml_peg/analysis/physicality/iron_properties/analyse_iron_properties.py b/ml_peg/analysis/physicality/iron_properties/analyse_iron_properties.py new file mode 100644 index 000000000..76ae8dc47 --- /dev/null +++ b/ml_peg/analysis/physicality/iron_properties/analyse_iron_properties.py @@ -0,0 +1,544 @@ +""" +Analyse BCC iron properties benchmark. + +This analysis combines EOS, elastic, Bain path, defect, surface, and stacking fault +properties. + +Reference +---------- +Zhang, L., Csányi, G., van der Giessen, E., & Maresca, F. (2023). +Efficiency, Accuracy, and Transferability of Machine Learning Potentials: +Application to Dislocations and Cracks in Iron. +""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any + +import numpy as np +import pandas as pd +import pytest + +from ml_peg.analysis.utils.decorators import build_table +from ml_peg.analysis.utils.utils import load_metrics_config +from ml_peg.app import APP_ROOT +from ml_peg.calcs import CALCS_ROOT +from ml_peg.models.get_models import get_model_names +from ml_peg.models.models import current_models + +MODELS = get_model_names(current_models) +CALC_PATH = CALCS_ROOT / "physicality" / "iron_properties" / "outputs" +OUT_PATH = APP_ROOT / "data" / "physicality" / "iron_properties" + +METRICS_CONFIG_PATH = Path(__file__).with_name("metrics.yml") +DEFAULT_THRESHOLDS, DEFAULT_TOOLTIPS, _ = load_metrics_config(METRICS_CONFIG_PATH) + +# DFT reference values +DFT_REFERENCE = { + # EOS properties + "a0": 2.831, # Lattice parameter (Å) + "B0": 178.0, # Bulk modulus (GPa) + "E_bcc_fcc": 83.5, # BCC-FCC energy difference (meV/atom) + # Elastic constants (GPa) + "C11": 296.7, + "C12": 151.4, + "C44": 104.7, + # Defect properties + "E_vac": 2.02, # Vacancy formation energy (eV) + "gamma_100": 2.41, # Surface energy (J/m²) + "gamma_110": 2.37, + "gamma_111": 2.58, + "gamma_112": 2.48, + "gamma_us_110": 0.75, # Unstable SFE (J/m²) + "gamma_us_112": 1.12, + # Traction-separation properties + "max_traction_100": 35.0, # Max traction for (100) cleavage (GPa) + "max_traction_110": 30.0, # Max traction for (110) cleavage (GPa) +} + +# DFT reference curves: (x_array, y_array) for plotting against MLIP results. +# fmt: off +DFT_CURVES = { + # Bain path: c/a ratio vs energy relative to BCC minimum (meV/atom) + "bain": ( + np.array([ + 0.698931, 0.746072, 0.788498, 0.815211, 0.878064, 0.923633, + 0.969202, 0.999057, 1.058768, 1.102766, 1.148334, 1.192332, + 1.237901, 1.283470, 1.327467, 1.371464, 1.413891, 1.462602, + 1.506600, 1.550597, 1.596166, 1.632307, 1.685732, 1.731301, + 1.776870, 1.820867, 1.864865, 1.910434, 1.956003, 2.000000, + ]), + np.array([ + 318.348, 225.434, 168.930, 132.529, 58.100, 21.163, + 3.248, 0.000, 9.265, 25.045, 45.717, 70.736, + 97.387, 123.493, 144.164, 156.140, 162.137, 157.267, + 144.242, 127.414, 116.020, 110.058, 111.168, 118.253, + 135.121, 157.423, 186.246, 221.048, 261.829, 308.044, + ]), + ), + # SFE {110}<111>: displacement (Å) vs SFE (J/m²) + "sfe_110": ( + np.array([0.24948, 0.49491, 0.73639, 0.97793, 1.22288, + 1.47228, 1.71376, 1.96316, 2.21261]), + np.array([0.13725, 0.44608, 0.74265, 0.91912, 0.97549, + 0.91422, 0.74510, 0.44608, 0.13480]), + ), + # SFE {112}<111>: displacement (Å) vs SFE (J/m²) + "sfe_112": ( + np.array([0.24948, 0.48691, 0.73639, 0.97793, 1.22343, + 1.46893, 1.71031, 1.95971, 2.20911]), + np.array([0.16137, 0.51100, 0.83130, 1.00733, 1.09780, + 1.11980, 0.93399, 0.55990, 0.17604]), + ), + # Traction-separation (100): separation (Å) vs traction (GPa) + "ts_100": ( + np.array([ + 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, + 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, + 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, + ]), + np.array([ + 0.022, 12.350, 21.600, 28.048, 32.251, 34.000, + 35.045, 34.788, 33.604, 31.294, 27.635, 23.440, + 20.571, 18.231, 16.387, 14.923, 13.783, 12.695, + 11.477, 10.608, 9.565, 8.949, 8.095, 7.273, + 6.499, 5.864, 5.176, 4.694, 4.101, 3.488, + 3.372, 2.903, 2.691, 2.228, 1.831, 1.291, + 1.270, 1.492, 0.971, 1.142, 0.629, + ]), + ), + # Traction-separation (110): separation (Å) vs traction (GPa) + "ts_110": ( + np.array([ + 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, + 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, + 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, + 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8, 3.9, 4.0, + ]), + np.array([ + 0.043, 13.291, 22.902, 29.408, 33.134, 34.666, + 34.930, 33.885, 32.144, 30.322, 28.024, 25.722, + 22.988, 20.383, 18.169, 16.304, 14.363, 12.924, + 11.271, 10.082, 8.839, 7.597, 6.480, 5.937, + 5.520, 4.540, 3.986, 3.310, 3.177, 2.730, + 2.518, 2.336, 1.620, 1.796, 1.625, 1.307, + 1.171, 1.283, 0.892, 0.916, 0.320, + ]), + ), +} +# fmt: on + + +# Curve file mapping for CSV loading +CURVE_FILES = { + "eos": "eos_curve.csv", + "bain": "bain_path.csv", + "sfe_110": "sfe_110_curve.csv", + "sfe_112": "sfe_112_curve.csv", + "ts_100": "ts_100_curve.csv", + "ts_110": "ts_110_curve.csv", +} + + +def load_model_results(model_name: str) -> dict[str, Any] | None: + """ + Load iron properties results for a model. + + Parameters + ---------- + model_name : str + Name of the model to load results for. + + Returns + ------- + dict[str, Any] | None + Dictionary of results, or None if file does not exist. + """ + json_path = CALC_PATH / model_name / "results.json" + if not json_path.exists(): + return None + return json.loads(json_path.read_text()) + + +def load_curve(model_name: str, curve_type: str) -> pd.DataFrame: + """ + Load curve data for a model. + + Parameters + ---------- + model_name : str + Name of the model to load curve for. + curve_type : str + Type of curve to load (e.g., 'eos', 'bain', 'sfe_110'). + + Returns + ------- + pd.DataFrame + Curve data, or empty DataFrame if file does not exist. + """ + filename = CURVE_FILES.get(curve_type) + if not filename: + return pd.DataFrame() + csv_path = CALC_PATH / model_name / filename + if not csv_path.exists(): + return pd.DataFrame() + return pd.read_csv(csv_path) + + +def compute_metrics(results: dict[str, Any]) -> dict[str, float]: + """ + Compute metrics from model results. + + Parameters + ---------- + results : dict[str, Any] + Dictionary containing model calculation results. + + Returns + ------- + dict[str, float] + Dictionary mapping metric names to computed values. + """ + metrics: dict[str, float] = {} + + # ========================================================================== + # EOS metrics + # ========================================================================== + eos = results.get("eos", {}) + if "a0" in eos: + a0_mlip = eos["a0"] + a0_error = abs(a0_mlip - DFT_REFERENCE["a0"]) / DFT_REFERENCE["a0"] * 100 + metrics["a0 error (%)"] = a0_error + + if "B0" in eos: + B0_mlip = eos["B0"] # noqa: N806 + B0_error = abs(B0_mlip - DFT_REFERENCE["B0"]) / DFT_REFERENCE["B0"] * 100 # noqa: N806 + metrics["B0 error (%)"] = B0_error + + # ========================================================================== + # Bain path metrics + # ========================================================================== + bain = results.get("bain_path", {}) + if "delta_E_meV" in bain: + E_bcc_fcc_mlip = bain["delta_E_meV"] # noqa: N806 + E_bcc_fcc_error = abs(E_bcc_fcc_mlip - DFT_REFERENCE["E_bcc_fcc"]) # noqa: N806 + metrics["BCC-FCC ΔE error (meV)"] = E_bcc_fcc_error + + # ========================================================================== + # Elastic constants metrics + # ========================================================================== + elastic = results.get("elastic", {}) + if "C11" in elastic: + C11_error = ( # noqa: N806 + abs(elastic["C11"] - DFT_REFERENCE["C11"]) / DFT_REFERENCE["C11"] * 100 + ) + metrics["C11 error (%)"] = C11_error + if "C12" in elastic: + C12_error = ( # noqa: N806 + abs(elastic["C12"] - DFT_REFERENCE["C12"]) / DFT_REFERENCE["C12"] * 100 + ) + metrics["C12 error (%)"] = C12_error + if "C44" in elastic: + C44_error = ( # noqa: N806 + abs(elastic["C44"] - DFT_REFERENCE["C44"]) / DFT_REFERENCE["C44"] * 100 + ) + metrics["C44 error (%)"] = C44_error + + # ========================================================================== + # Vacancy metrics + # ========================================================================== + vacancy = results.get("vacancy", {}) + if "E_vac" in vacancy: + E_vac_mlip = vacancy["E_vac"] # noqa: N806 + E_vac_error = ( # noqa: N806 + abs(E_vac_mlip - DFT_REFERENCE["E_vac"]) / DFT_REFERENCE["E_vac"] * 100 + ) + metrics["E_vac error (%)"] = E_vac_error + + # ========================================================================== + # Surface energy metrics + # ========================================================================== + surfaces = results.get("surfaces", {}) + surface_errors = [] + + for surface in ["100", "110", "111", "112"]: + key_mlip = f"gamma_{surface}" + if key_mlip in surfaces: + gamma_mlip = surfaces[key_mlip] + gamma_dft = DFT_REFERENCE[key_mlip] + error = abs(gamma_mlip - gamma_dft) + surface_errors.append(error) + + if surface_errors: + metrics["Surface MAE (J/m²)"] = np.mean(surface_errors) + + # ========================================================================== + # Stacking fault metrics + # ========================================================================== + sfe_110 = results.get("sfe_110", {}) + if "max_sfe" in sfe_110: + max_sfe_110_mlip = sfe_110["max_sfe"] + max_sfe_110_error = ( + abs(max_sfe_110_mlip - DFT_REFERENCE["gamma_us_110"]) + / DFT_REFERENCE["gamma_us_110"] + * 100 + ) + metrics["Max SFE 110 error (%)"] = max_sfe_110_error + + sfe_112 = results.get("sfe_112", {}) + if "max_sfe" in sfe_112: + max_sfe_112_mlip = sfe_112["max_sfe"] + max_sfe_112_error = ( + abs(max_sfe_112_mlip - DFT_REFERENCE["gamma_us_112"]) + / DFT_REFERENCE["gamma_us_112"] + * 100 + ) + metrics["Max SFE 112 error (%)"] = max_sfe_112_error + + # ========================================================================== + # Traction-separation metrics + # ========================================================================== + ts_100 = results.get("ts_100", {}) + if "max_traction" in ts_100: + traction_100_error = ( + abs(ts_100["max_traction"] - DFT_REFERENCE["max_traction_100"]) + / DFT_REFERENCE["max_traction_100"] + * 100 + ) + metrics["Max traction (100) error (%)"] = traction_100_error + + ts_110 = results.get("ts_110", {}) + if "max_traction" in ts_110: + traction_110_error = ( + abs(ts_110["max_traction"] - DFT_REFERENCE["max_traction_110"]) + / DFT_REFERENCE["max_traction_110"] + * 100 + ) + metrics["Max traction (110) error (%)"] = traction_110_error + + return metrics + + +def _load_all_results() -> dict[str, dict[str, Any]]: + """ + Load results for all models. + + Returns + ------- + dict[str, dict[str, Any]] + Dictionary mapping model names to their results. + """ + all_results: dict[str, dict[str, Any]] = {} + for model_name in MODELS: + results = load_model_results(model_name) + if results is not None: + all_results[model_name] = results + return all_results + + +def _load_curves_for_all_models(curve_type: str) -> dict[str, pd.DataFrame]: + """ + Load curves of given type for all models. + + Parameters + ---------- + curve_type : str + Type of curve to load (e.g., 'eos', 'bain', 'sfe_110'). + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their curve DataFrames. + """ + curves: dict[str, pd.DataFrame] = {} + for model_name in MODELS: + curve = load_curve(model_name, curve_type) + if not curve.empty: + curves[model_name] = curve + return curves + + +@pytest.fixture +def iron_eos_curves() -> dict[str, pd.DataFrame]: + """ + Load EOS curves for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their EOS curve DataFrames. + """ + return _load_curves_for_all_models("eos") + + +@pytest.fixture +def iron_bain_curves() -> dict[str, pd.DataFrame]: + """ + Load Bain path curves for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their Bain path curve DataFrames. + """ + return _load_curves_for_all_models("bain") + + +@pytest.fixture +def iron_sfe_110_curves() -> dict[str, pd.DataFrame]: + """ + Load SFE 110 curves for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their SFE 110 curve DataFrames. + """ + return _load_curves_for_all_models("sfe_110") + + +@pytest.fixture +def iron_sfe_112_curves() -> dict[str, pd.DataFrame]: + """ + Load SFE 112 curves for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their SFE 112 curve DataFrames. + """ + return _load_curves_for_all_models("sfe_112") + + +@pytest.fixture +def iron_ts_100_curves() -> dict[str, pd.DataFrame]: + """ + Load T-S (100) curves for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their T-S (100) curve DataFrames. + """ + return _load_curves_for_all_models("ts_100") + + +@pytest.fixture +def iron_ts_110_curves() -> dict[str, pd.DataFrame]: + """ + Load T-S (110) curves for all models. + + Returns + ------- + dict[str, pd.DataFrame] + Dictionary mapping model names to their T-S (110) curve DataFrames. + """ + return _load_curves_for_all_models("ts_110") + + +def collect_metrics() -> pd.DataFrame: + """ + Gather metrics for all models. + + Returns + ------- + pd.DataFrame + DataFrame containing metrics for all models. + """ + metrics_rows: list[dict[str, float | str]] = [] + + OUT_PATH.mkdir(parents=True, exist_ok=True) + + all_results = _load_all_results() + + for model_name, results in all_results.items(): + model_metrics = compute_metrics(results) + row = {"Model": model_name} | model_metrics + metrics_rows.append(row) + + columns = ["Model"] + list(DEFAULT_THRESHOLDS.keys()) + + return pd.DataFrame(metrics_rows).reindex(columns=columns) + + +@pytest.fixture +def iron_properties_collection() -> pd.DataFrame: + """ + Collect iron properties metrics across all models. + + Returns + ------- + pd.DataFrame + DataFrame containing iron properties metrics for all models. + """ + return collect_metrics() + + +@pytest.fixture +def iron_properties_metrics_dataframe( + iron_properties_collection: pd.DataFrame, +) -> pd.DataFrame: + """ + Provide the aggregated iron properties metrics dataframe. + + Parameters + ---------- + iron_properties_collection : pd.DataFrame + Collection of iron properties metrics. + + Returns + ------- + pd.DataFrame + The aggregated iron properties metrics DataFrame. + """ + return iron_properties_collection + + +@pytest.fixture +@build_table( + filename=OUT_PATH / "iron_properties_metrics_table.json", + metric_tooltips=DEFAULT_TOOLTIPS, + thresholds=DEFAULT_THRESHOLDS, + weights=None, +) +def metrics( + iron_properties_metrics_dataframe: pd.DataFrame, +) -> dict[str, dict]: + """ + Compute iron properties metrics for all models. + + Parameters + ---------- + iron_properties_metrics_dataframe + Aggregated per-model metrics. + + Returns + ------- + dict[str, dict] + Mapping of metric names to per-model results. + """ + metrics_df = iron_properties_metrics_dataframe + metrics_dict: dict[str, dict[str, float | None]] = {} + for column in metrics_df.columns: + if column == "Model": + continue + values = [ + value if pd.notna(value) else None for value in metrics_df[column].tolist() + ] + metrics_dict[column] = dict(zip(metrics_df["Model"], values, strict=False)) + return metrics_dict + + +def test_iron_properties(metrics: dict[str, dict]) -> None: + """ + Run iron properties analysis. + + Parameters + ---------- + metrics : dict[str, dict] + Dictionary of iron properties metrics from the metrics fixture. + """ + return diff --git a/ml_peg/analysis/physicality/iron_properties/metrics.yml b/ml_peg/analysis/physicality/iron_properties/metrics.yml new file mode 100644 index 000000000..011274642 --- /dev/null +++ b/ml_peg/analysis/physicality/iron_properties/metrics.yml @@ -0,0 +1,77 @@ +metrics: + # EOS properties + a0 error (%): + good: 0.0 + bad: 2.0 + unit: "%" + tooltip: "Lattice parameter error relative to DFT" + level_of_theory: PBE + B0 error (%): + good: 0.0 + bad: 20.0 + unit: "%" + tooltip: "Bulk modulus error relative to DFT" + level_of_theory: PBE + BCC-FCC ΔE error (meV): + good: 0.0 + bad: 50.0 + unit: "meV" + tooltip: "Error in BCC-FCC energy difference along Bain path" + level_of_theory: PBE + # Elastic constants + C11 error (%): + good: 0.0 + bad: 20.0 + unit: "%" + tooltip: "Elastic constant C11 error relative to DFT (ref: 243 GPa)" + level_of_theory: PBE + C12 error (%): + good: 0.0 + bad: 20.0 + unit: "%" + tooltip: "Elastic constant C12 error relative to DFT (ref: 145 GPa)" + level_of_theory: PBE + C44 error (%): + good: 0.0 + bad: 20.0 + unit: "%" + tooltip: "Elastic constant C44 error relative to DFT (ref: 116 GPa)" + level_of_theory: PBE + # Defect properties + E_vac error (%): + good: 0.0 + bad: 20.0 + unit: "%" + tooltip: "Vacancy formation energy error relative to DFT" + level_of_theory: PBE + Surface MAE (J/m²): + good: 0.0 + bad: 0.5 + unit: "J/m²" + tooltip: "Mean absolute error for surface energies (100, 110, 111, 112)" + level_of_theory: PBE + Max SFE 110 error (%): + good: 0.0 + bad: 30.0 + unit: "%" + tooltip: "Error in maximum stacking fault energy for {110}<111> slip system" + level_of_theory: PBE + Max SFE 112 error (%): + good: 0.0 + bad: 30.0 + unit: "%" + tooltip: "Error in maximum stacking fault energy for {112}<111> slip system" + level_of_theory: PBE + # Traction-separation properties + Max traction (100) error (%): + good: 0.0 + bad: 30.0 + unit: "%" + tooltip: "Maximum traction error relative to DFT (ref: 35 GPa)" + level_of_theory: PBE + Max traction (110) error (%): + good: 0.0 + bad: 30.0 + unit: "%" + tooltip: "Maximum traction error relative to DFT (ref: 30 GPa)" + level_of_theory: PBE diff --git a/ml_peg/app/physicality/iron_properties/app_iron_properties.py b/ml_peg/app/physicality/iron_properties/app_iron_properties.py new file mode 100644 index 000000000..0766f5462 --- /dev/null +++ b/ml_peg/app/physicality/iron_properties/app_iron_properties.py @@ -0,0 +1,319 @@ +"""Run iron properties app.""" + +from __future__ import annotations + +from dash import Dash, Input, Output, callback, dcc +from dash.dcc import Loading +from dash.exceptions import PreventUpdate +from dash.html import Div, Label +import numpy as np +import pandas as pd +import plotly.graph_objects as go + +from ml_peg.analysis.physicality.iron_properties.analyse_iron_properties import ( + DFT_CURVES, +) +from ml_peg.app import APP_ROOT +from ml_peg.app.base_app import BaseApp +from ml_peg.calcs import CALCS_ROOT +from ml_peg.models.get_models import get_model_names +from ml_peg.models.models import current_models + +# Get all models +MODELS = get_model_names(current_models) +BENCHMARK_NAME = "Iron Properties" +DATA_PATH = APP_ROOT / "data" / "physicality" / "iron_properties" +CALC_PATH = CALCS_ROOT / "physicality" / "iron_properties" / "outputs" +DOCS_URL = "https://ddmms.github.io/ml-peg/user_guide/benchmarks/physicality.html#iron-properties" + +# Curve configuration: file name, x column, y column, title, x label, y label +CURVE_CONFIG = { + "eos": { + "file": "eos_curve.csv", + "x": "volume", + "y": "energy", + "title": "Equation of State", + "x_label": "Volume (ų/atom)", + "y_label": "Energy (eV/atom)", + }, + "bain": { + "file": "bain_path.csv", + "x": "ca_ratio", + "y": "energy_meV", + "title": "Bain Path", + "x_label": "c/a ratio", + "y_label": "Energy (meV/atom)", + }, + "sfe_110": { + "file": "sfe_110_curve.csv", + "x": "displacement", + "y": "sfe_J_per_m2", + "title": "Stacking Fault Energy {110}<111>", + "x_label": "Displacement (Å)", + "y_label": "SFE (J/m²)", + }, + "sfe_112": { + "file": "sfe_112_curve.csv", + "x": "displacement", + "y": "sfe_J_per_m2", + "title": "Stacking Fault Energy {112}<111>", + "x_label": "Displacement (Å)", + "y_label": "SFE (J/m²)", + }, + "ts_100": { + "file": "ts_100_curve.csv", + "x": "separation", + "y": "traction", + "title": "Traction-Separation Curve (100)", + "x_label": "Separation (Å)", + "y_label": "Traction (GPa)", + }, + "ts_110": { + "file": "ts_110_curve.csv", + "x": "separation", + "y": "traction", + "title": "Traction-Separation Curve (110)", + "x_label": "Separation (Å)", + "y_label": "Traction (GPa)", + }, +} + + +def _load_curve_data(model_name: str, curve_type: str) -> pd.DataFrame | None: + """ + Load curve data for a model. + + Parameters + ---------- + model_name : str + Name of the model. + curve_type : str + Type of curve to load (e.g., 'eos', 'bain', 'sfe_110'). + + Returns + ------- + pd.DataFrame or None + DataFrame with curve data, or None if not found. + """ + model_dir = CALC_PATH / model_name + if not model_dir.exists(): + return None + + config = CURVE_CONFIG.get(curve_type) + if not config: + return None + + csv_path = model_dir / config["file"] + if not csv_path.exists(): + return None + + return pd.read_csv(csv_path) + + +def _create_figure(df: pd.DataFrame, curve_type: str, model_name: str) -> go.Figure: + """ + Create plotly figure for the given curve type. + + Parameters + ---------- + df : pd.DataFrame + DataFrame containing the curve data. + curve_type : str + Type of curve to plot (e.g., 'eos', 'bain', 'sfe_110'). + model_name : str + Name of the model for the title. + + Returns + ------- + go.Figure + Plotly figure object. + """ + config = CURVE_CONFIG.get(curve_type, {}) + + fig = go.Figure() + + # Get model data x-range for limiting DFT reference + model_x_max = df[config["x"]].max() if config.get("x") else None + + # Add DFT reference curve if available + dft_data = DFT_CURVES.get(curve_type) + if dft_data is not None: + x_dft, y_dft = dft_data[0].copy(), dft_data[1].copy() + + # For SFE curves, limit DFT data to model x-range (data is periodic) + if curve_type.startswith("sfe_") and model_x_max is not None: + mask = x_dft <= model_x_max + x_dft = np.array(x_dft)[mask] + y_dft = np.array(y_dft)[mask] + + fig.add_trace( + go.Scatter( + x=x_dft, + y=y_dft, + mode="lines", + name="DFT Reference", + line={"width": 2, "dash": "dash", "color": "gray"}, + ) + ) + + # Add model curve + fig.add_trace( + go.Scatter( + x=df[config["x"]], + y=df[config["y"]], + mode="lines+markers", + name=model_name, + line={"width": 2}, + marker={"size": 6}, + ) + ) + + # Special handling for Bain path (add BCC/FCC reference lines) + if curve_type == "bain": + fig.add_vline(x=1.0, line_dash="dash", line_color="gray", annotation_text="BCC") + fig.add_vline( + x=1.414, line_dash="dash", line_color="gray", annotation_text="FCC" + ) + + fig.update_layout( + title=f"{config['title']} - {model_name}", + xaxis_title=config["x_label"], + yaxis_title=config["y_label"], + template="plotly_white", + showlegend=True, + height=500, + ) + + return fig + + +class IronPropertiesApp(BaseApp): + """Iron properties benchmark app layout and callbacks.""" + + def register_callbacks(self) -> None: + """Register callbacks for curve visualization.""" + model_dropdown_id = f"{BENCHMARK_NAME}-model-dropdown" + curve_dropdown_id = f"{BENCHMARK_NAME}-curve-dropdown" + figure_id = f"{BENCHMARK_NAME}-figure" + + @callback( + Output(figure_id, "figure"), + Input(model_dropdown_id, "value"), + Input(curve_dropdown_id, "value"), + ) + def update_figure(model_name: str, curve_type: str) -> go.Figure: # noqa: F811 + # Invoked by Dash's callback system, not called directly. + """ + Update figure based on model and curve selection. + + Parameters + ---------- + model_name : str + Name of the selected model. + curve_type : str + Type of curve to display. + + Returns + ------- + go.Figure + Updated plotly figure. + """ + if not model_name or not curve_type: + raise PreventUpdate + + df = _load_curve_data(model_name, curve_type) + if df is None or df.empty: + # Return empty figure with message + fig = go.Figure() + fig.add_annotation( + text=f"No data available for {model_name} - {curve_type}", + xref="paper", + yref="paper", + x=0.5, + y=0.5, + showarrow=False, + font={"size": 16}, + ) + fig.update_layout( + template="plotly_white", + height=500, + ) + return fig + + return _create_figure(df, curve_type, model_name) + + +def get_app() -> IronPropertiesApp: + """ + Get iron properties benchmark app layout and callback registration. + + Returns + ------- + IronPropertiesApp + Benchmark layout and callback registration. + """ + model_options = [{"label": model, "value": model} for model in MODELS] + default_model = model_options[0]["value"] if model_options else None + + extra_components = [ + Div( + [ + Label("Select model for curve visualization:"), + dcc.Dropdown( + id=f"{BENCHMARK_NAME}-model-dropdown", + options=model_options, + value=default_model, + clearable=False, + style={"width": "300px", "marginBottom": "20px"}, + ), + dcc.Dropdown( + id=f"{BENCHMARK_NAME}-curve-dropdown", + options=[ + {"label": "EOS Curve", "value": "eos"}, + {"label": "Bain Path", "value": "bain"}, + {"label": "SFE {110}<111>", "value": "sfe_110"}, + {"label": "SFE {112}<111>", "value": "sfe_112"}, + {"label": "T-S Curve (100)", "value": "ts_100"}, + {"label": "T-S Curve (110)", "value": "ts_110"}, + ], + value="eos", + clearable=False, + style={"width": "300px"}, + ), + ], + style={"marginBottom": "20px"}, + ), + Loading( + dcc.Graph( + id=f"{BENCHMARK_NAME}-figure", + style={"height": "500px", "width": "100%", "marginTop": "20px"}, + ), + type="circle", + ), + ] + + return IronPropertiesApp( + name=BENCHMARK_NAME, + description=( + "Comprehensive BCC iron properties benchmark. " + "Includes equation of state (lattice parameter, bulk modulus), " + "elastic constants (C11, C12, C44), Bain path (BCC-FCC transformation), " + "vacancy formation energy, surface energies (100, 110, 111, 112), " + "generalized stacking fault energy curves for {110}<111> and " + "{112}<111> slip systems, and traction-separation curves for (100) " + "and (110) cleavage planes. " + "This benchmark is computationally expensive and marked with " + "@pytest.mark.slow." + ), + docs_url=DOCS_URL, + table_path=DATA_PATH / "iron_properties_metrics_table.json", + extra_components=extra_components, + ) + + +if __name__ == "__main__": + dash_app = Dash(__name__, assets_folder=DATA_PATH.parent.parent) + iron_properties_app = get_app() + dash_app.layout = iron_properties_app.layout + iron_properties_app.register_callbacks() + dash_app.run(port=8060, debug=True) diff --git a/ml_peg/calcs/physicality/iron_properties/calc_iron_properties.py b/ml_peg/calcs/physicality/iron_properties/calc_iron_properties.py new file mode 100644 index 000000000..0913a37df --- /dev/null +++ b/ml_peg/calcs/physicality/iron_properties/calc_iron_properties.py @@ -0,0 +1,817 @@ +""" +Run calculations for BCC iron properties benchmark. + +This benchmark computes fundamental properties of BCC iron including: +- Equation of state (lattice parameter, bulk modulus) +- Elastic constants (C11, C12, C44) +- Bain path energy curve +- Vacancy formation energy +- Surface energies (100, 110, 111, 112) +- Generalized stacking fault energy curves (110, 112) +- Traction-separation curves (100, 110) + +This benchmark is computationally expensive and marked with @pytest.mark.slow. +""" + +from __future__ import annotations + +import json +from pathlib import Path +from typing import Any + +from ase.build import bulk +from ase.constraints import FixedLine +from ase.filters import ExpCellFilter +from ase.optimize import BFGS +import numpy as np +import pandas as pd +import pytest + +from ml_peg.calcs.physicality.iron_properties.iron_utils import ( + EV_PER_A2_TO_J_PER_M2, + EV_PER_A3_TO_GPA, + apply_voigt_strain, + calculate_surface_energy, + create_bain_cell, + create_bcc_supercell, + create_sfe_110_structure, + create_sfe_112_structure, + create_surface_100, + create_surface_110, + create_surface_111, + create_surface_112, + fit_eos, + relax_volume_isotropic, +) +from ml_peg.models.get_models import load_models +from ml_peg.models.models import current_models + +MODELS = load_models(current_models) + +# Local directory for calculator outputs +OUT_PATH = Path(__file__).parent / "outputs" + +# ============================================================================= +# Test Parameters +# ============================================================================= + +# EOS calculation parameters +EOS_NUM_POINTS = 30 + +# BFGS optimization parameters +BFGS_FMAX = 1e-5 +BFGS_MAX_ITER = 100 + +# Elastic constants parameters +ELASTIC_STRAIN = 1.0e-5 +ELASTIC_SUPERCELL_SIZE = (4, 4, 4) +ELASTIC_ATOM_JIGGLE = 1.0e-5 # Random perturbation to prevent saddle points + +# Bain path parameters +BAIN_NUM_POINTS = 65 + +# Vacancy calculation parameters +VACANCY_SUPERCELL_SIZE = (4, 4, 4) + +# Surface calculation parameters +SURFACE_VACUUM = 10.0 # Angstroms + +# Stacking fault calculation parameters +SFE_STEPS = 16 # Number of steps to cover one Burgers vector + +# Traction-separation parameters +TS_MAX_SEPARATION = 5.0 # Angstroms +TS_STEP_SIZE = 0.05 # Angstroms + + +# ============================================================================= +# EOS Calculation +# ============================================================================= + + +def run_eos_calculation(calc: Any) -> dict[str, Any]: + """ + Run the energy-volume curve calculation. + + Parameters + ---------- + calc + ASE calculator object. + + Returns + ------- + dict + Dictionary with EOS results including a0, B0, V0, E0, volumes, energies. + """ + # Generate lattice parameters: 2.834 - 0.05 + (0.1/30)*i for i in 1..30 + lattice_params = np.array( + [2.834 - 0.05 + 0.1 / 30 * i for i in range(1, EOS_NUM_POINTS + 1)] + ) + + volumes = [] + energies = [] + + for lat in lattice_params: + atoms = bulk("Fe", "bcc", a=lat, cubic=True) + atoms.calc = calc + + # Relax atomic positions at fixed cell volume + # (matches LAMMPS minimize behavior) + opt = BFGS(atoms, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + + energy = atoms.get_potential_energy() + volume = atoms.get_volume() + + n_atoms = len(atoms) + volumes.append(volume / n_atoms) + energies.append(energy / n_atoms) + + volumes = np.array(volumes) + energies = np.array(energies) + + # Fit Birch-Murnaghan EOS + eos_results = fit_eos(volumes, energies) + + return { + "volumes": volumes.tolist(), + "energies": energies.tolist(), + "lattice_params": lattice_params.tolist(), + "a0": eos_results["a0"], + "E0": eos_results["E0"], + "B0": eos_results["B0"], + "Bp": eos_results["Bp"], + "V0": eos_results["V0"], + } + + +# ============================================================================= +# Elastic Constants Calculation +# ============================================================================= + + +def run_elastic_calculation(calc: Any, lattice_parameter: float) -> dict[str, Any]: + """ + Calculate elastic constants using the stress-strain method. + + Parameters + ---------- + calc + ASE calculator object. + lattice_parameter + Equilibrium lattice parameter from EOS fit. + + Returns + ------- + dict + Dictionary with elastic constants C11, C12, C44, bulk_modulus. + """ + # Create supercell + atoms_ref = create_bcc_supercell(lattice_parameter, ELASTIC_SUPERCELL_SIZE) + atoms_ref.calc = calc + + # Box relaxation + ecf = ExpCellFilter(atoms_ref) + opt = BFGS(ecf, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + + # Apply random jiggle to atoms to prevent staying on saddle points + rng = np.random.default_rng(seed=87287) + jiggle = rng.uniform( + -ELASTIC_ATOM_JIGGLE, ELASTIC_ATOM_JIGGLE, atoms_ref.positions.shape + ) + atoms_ref.positions += jiggle + + # Elastic constant matrix + C = np.zeros((6, 6)) # noqa: N806 + + for i in range(6): + direction = i + 1 + + # Positive strain with off-diagonal cell adjustment + atoms_pos = apply_voigt_strain(atoms_ref.copy(), direction, ELASTIC_STRAIN) + atoms_pos.calc = calc + + opt_pos = BFGS(atoms_pos, logfile=None) + opt_pos.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + stress_pos = atoms_pos.get_stress(voigt=True) + + # Negative strain with off-diagonal cell adjustment + atoms_neg = apply_voigt_strain(atoms_ref.copy(), direction, -ELASTIC_STRAIN) + atoms_neg.calc = calc + + opt_neg = BFGS(atoms_neg, logfile=None) + opt_neg.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + stress_neg = atoms_neg.get_stress(voigt=True) + + # Compute elastic constants using stress differences + # C_ij = dσ_i / dε_j = (σ_pos - σ_neg) / (2 * ε) + delta_stress = stress_pos - stress_neg + delta_strain = 2 * ELASTIC_STRAIN + + for j in range(6): + C[j, i] = delta_stress[j] / delta_strain * EV_PER_A3_TO_GPA + + # Symmetrize + C_sym = 0.5 * (C + C.T) # noqa: N806 + + # Extract cubic averages + C11 = (C_sym[0, 0] + C_sym[1, 1] + C_sym[2, 2]) / 3 # noqa: N806 + C12 = (C_sym[0, 1] + C_sym[0, 2] + C_sym[1, 2]) / 3 # noqa: N806 + C44 = (C_sym[3, 3] + C_sym[4, 4] + C_sym[5, 5]) / 3 # noqa: N806 + + bulk_modulus = (C11 + 2 * C12) / 3 + + return { + "C11": C11, + "C12": C12, + "C44": C44, + "bulk_modulus": bulk_modulus, + "C_matrix": C_sym.tolist(), + } + + +# ============================================================================= +# Bain Path Calculation +# ============================================================================= + + +def run_bain_path_calculation(calc: Any, lattice_parameter: float) -> dict[str, Any]: + """ + Calculate the Bain path energy curve. + + For each target c/a ratio, creates a tetragonally distorted cell and performs + isotropic volume relaxation (uniform scaling only) to find the minimum energy + while maintaining the c/a ratio. This matches the LAMMPS behavior where + 'fix box/relax aniso 0.0 couple xyz' is used. + + The 'couple xyz' constraint in LAMMPS couples all three diagonal stress + components together, meaning x, y, and z dimensions can only change by the + same fractional amount. This preserves the c/a ratio during relaxation. + + Parameters + ---------- + calc + ASE calculator object. + lattice_parameter + Equilibrium BCC lattice parameter. + + Returns + ------- + dict + Dictionary with ca_ratios, energies, E_bcc, E_fcc, delta_E. + """ + # Generate c/a ratios: 0.7 + 0.02*i for i in 1..65 + ca_ratios_target = np.array([0.7 + 0.02 * i for i in range(1, BAIN_NUM_POINTS + 1)]) + + ca_ratios = [] + energies = [] + + for ratio in ca_ratios_target: + # Create tetragonally distorted cell at target c/a ratio + atoms = create_bain_cell(lattice_parameter, ratio) + atoms.calc = calc + + # Step 1: Isotropic volume relaxation (maintains c/a ratio) + # This is equivalent to LAMMPS: fix box/relax aniso 0.0 couple xyz + # Only uniform scaling is allowed, preserving the cell shape + atoms_relaxed = relax_volume_isotropic(atoms, calc) + + # Step 2: Atomic position relaxation at fixed cell + opt = BFGS(atoms_relaxed, logfile=None) + try: + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception: + pass + + energy = atoms_relaxed.get_potential_energy() + cell = atoms_relaxed.get_cell() + ca_actual = cell[2, 2] / cell[1, 1] + + n_atoms = len(atoms_relaxed) + ca_ratios.append(ca_actual) + energies.append(energy / n_atoms) + + ca_ratios = np.array(ca_ratios) + energies = np.array(energies) + + # Normalize energies (subtract minimum, convert to meV) + E_min = np.min(energies) # noqa: N806 + energies_norm = (energies - E_min) * 1000 # meV/atom + + # Find BCC point (c/a ≈ 1.0) and FCC point (c/a ≈ 1.414) + idx_bcc = np.argmin(np.abs(ca_ratios - 1.0)) + idx_fcc = np.argmin(np.abs(ca_ratios - np.sqrt(2))) + + E_bcc = energies_norm[idx_bcc] # noqa: N806 + E_fcc = energies_norm[idx_fcc] # noqa: N806 + + return { + "ca_ratios": ca_ratios.tolist(), + "energies": energies.tolist(), + "energies_meV": energies_norm.tolist(), + "E_bcc_meV": E_bcc, + "E_fcc_meV": E_fcc, + "delta_E_meV": E_fcc - E_bcc, + } + + +# ============================================================================= +# Vacancy Calculation +# ============================================================================= + + +def run_vacancy_calculation(calc: Any, lattice_parameter: float) -> dict[str, Any]: + """ + Calculate the vacancy formation energy. + + Parameters + ---------- + calc : Any + ASE calculator object. + lattice_parameter : float + Equilibrium lattice parameter from EOS fit. + + Returns + ------- + dict + Dictionary with vacancy results including E_vac, E_coh, E_perfect, E_defect. + """ + atoms_perfect = create_bcc_supercell(lattice_parameter, VACANCY_SUPERCELL_SIZE) + atoms_perfect.calc = calc + + n_atoms = len(atoms_perfect) + E_perfect = atoms_perfect.get_potential_energy() # noqa: N806 + E_coh = E_perfect / n_atoms # noqa: N806 + + atoms_defect = atoms_perfect.copy() + del atoms_defect[0] + atoms_defect.calc = calc + + opt = BFGS(atoms_defect, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + E_defect = atoms_defect.get_potential_energy() # noqa: N806 + E_vac = (E_defect - E_perfect) + E_coh # noqa: N806 + + return { + "E_vac": E_vac, + "E_coh": E_coh, + "E_perfect": E_perfect, + "E_defect": E_defect, + } + + +# ============================================================================= +# Surface Calculations +# ============================================================================= + +# Surface configuration: create_fn, layers/size, area_axes, vacuum +SURFACE_CONFIG = { + "100": { + "create_fn": create_surface_100, + "layers": 10, + "area_axes": (0, 1), + "vacuum": SURFACE_VACUUM, + }, + "110": { + "create_fn": create_surface_110, + "layers": 10, + "area_axes": (0, 1), + "vacuum": SURFACE_VACUUM, + }, + "111": { + "create_fn": create_surface_111, + "size": (3, 15, 3), + "area_axes": (0, 2), + "vacuum": SURFACE_VACUUM, + }, + "112": { + "create_fn": create_surface_112, + "layers": 15, + "area_axes": (0, 1), + "vacuum": 5.0, + }, +} + + +def run_surface_calculations(calc: Any, lattice_parameter: float) -> dict[str, Any]: + """ + Calculate surface energies for (100), (110), (111), (112) surfaces. + + Parameters + ---------- + calc : Any + ASE calculator object. + lattice_parameter : float + Equilibrium lattice parameter from EOS fit. + + Returns + ------- + dict + Dictionary with surface energies gamma_100, gamma_110, gamma_111, gamma_112. + """ + surfaces = {} + + for name, cfg in SURFACE_CONFIG.items(): + print(f"Running surface calculation for {name}...") + create_fn = cfg["create_fn"] + area_axes = cfg["area_axes"] + vacuum = cfg["vacuum"] + + # Build kwargs for structure creation + if "size" in cfg: + bulk_kwargs = {"size": cfg["size"], "vacuum": 0.0} + slab_kwargs = {"size": cfg["size"], "vacuum": vacuum} + else: + bulk_kwargs = {"layers": cfg["layers"], "vacuum": 0.0} + slab_kwargs = {"layers": cfg["layers"], "vacuum": vacuum} + + # Bulk reference + atoms_bulk = create_fn(lattice_parameter, **bulk_kwargs) + atoms_bulk.calc = calc + e_bulk = atoms_bulk.get_potential_energy() + cell = atoms_bulk.get_cell() + area = np.linalg.norm(np.cross(cell[area_axes[0]], cell[area_axes[1]])) + + # Slab with vacuum + atoms_slab = create_fn(lattice_parameter, **slab_kwargs) + atoms_slab.calc = calc + opt = BFGS(atoms_slab, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + e_slab = atoms_slab.get_potential_energy() + + surfaces[name] = calculate_surface_energy(e_slab, e_bulk, area) + + return {f"gamma_{k}": v for k, v in surfaces.items()} + + +# ============================================================================= +# Stacking Fault Energy Calculations +# ============================================================================= + +# SFE configuration: create_fn, displacement axis +SFE_CONFIG = { + "110": {"create_fn": create_sfe_110_structure, "axis": 1}, + "112": {"create_fn": create_sfe_112_structure, "axis": 2}, +} + + +def run_sfe_calculation( + calc: Any, lattice_parameter: float, sfe_type: str +) -> dict[str, Any]: + """ + Calculate GSFE curve for specified slip system. + + Parameters + ---------- + calc : Any + ASE calculator object. + lattice_parameter : float + Equilibrium lattice parameter from EOS fit. + sfe_type : str + Type of SFE calculation ('110' or '112'). + + Returns + ------- + dict + Dictionary with displacements, sfe_J_per_m2, and max_sfe. + """ + config = SFE_CONFIG[sfe_type] + atoms = config["create_fn"](lattice_parameter) + atoms.calc = calc + + # Calculate Burgers vector magnitude: b = a * sqrt(3) / 2 + burgers_vector = lattice_parameter * np.sqrt(3) / 2 + step_size = burgers_vector / SFE_STEPS + + cell = atoms.get_cell() + ly = cell[1, 1] + lz = cell[2, 2] + area = ly * lz + + opt = BFGS(atoms, logfile=None) + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + e0 = atoms.get_potential_energy() + + positions = atoms.get_positions() + x_mid = (positions[:, 0].min() + positions[:, 0].max()) / 2 + 0.1 + upper_mask = positions[:, 0] < x_mid + upper_indices = np.where(upper_mask)[0] + + displacements = [0.0] + sfe_j_per_m2 = [0.0] + + constraints = [FixedLine(idx, direction=[1, 0, 0]) for idx in range(len(atoms))] + displacement_axis = config["axis"] + + for step in range(1, SFE_STEPS + 1): + positions = atoms.get_positions() + positions[upper_indices, displacement_axis] += step_size + atoms.set_positions(positions) + + atoms.set_constraint(constraints) + + opt = BFGS(atoms, logfile=None) + try: + opt.run(fmax=BFGS_FMAX, steps=BFGS_MAX_ITER) + except Exception: + pass + + atoms.set_constraint() + energy = atoms.get_potential_energy() + sfe = (energy - e0) / (2 * area) * EV_PER_A2_TO_J_PER_M2 + + displacements.append(step * step_size) + sfe_j_per_m2.append(sfe) + + return { + "displacements": displacements, + "sfe_J_per_m2": sfe_j_per_m2, + "max_sfe": max(sfe_j_per_m2), + } + + +# ============================================================================= +# Traction-Separation Calculations +# ============================================================================= + +# T-S configuration: structure creation function +TS_CONFIG = { + "100": lambda a: create_surface_100(a, layers=36, vacuum=0.0), + "110": lambda a: create_surface_110(a, layers=10, vacuum=0.0), +} + + +def run_ts_calculation( + calc: Any, lattice_parameter: float, direction: str +) -> dict[str, Any]: + """ + Calculate traction-separation curve for specified cleavage plane. + + The calculation incrementally separates crystal halves without relaxation + and measures energy and traction (stress from forces). + + Parameters + ---------- + calc : Any + ASE calculator object. + lattice_parameter : float + Equilibrium lattice parameter from EOS fit. + direction : str + Cleavage plane direction ('100' or '110'). + + Returns + ------- + dict + Dictionary with separations, energies, traction, and max_traction. + """ + create_fn = TS_CONFIG[direction] + num_steps = int(TS_MAX_SEPARATION / TS_STEP_SIZE) + 1 + + separations = [] + energies = [] + traction = [] + + for i in range(num_steps): + dd = TS_STEP_SIZE * i + + # Create fresh structure for each separation + atoms = create_fn(lattice_parameter) + atoms.calc = calc + + # Get cell dimensions + cell = atoms.get_cell() + lx = cell[0, 0] + ly = cell[1, 1] + lz = cell[2, 2] + area = lx * ly + + # Identify upper and lower atoms + positions = atoms.get_positions() + z_mid = lz / 2 - 0.1 + upper_mask = positions[:, 2] > z_mid + upper_indices = np.where(upper_mask)[0] + + # Expand cell in z direction + new_cell = cell.copy() + new_cell[2, 2] = lz + dd + atoms.set_cell(new_cell, scale_atoms=False) + + # Move upper atoms + positions = atoms.get_positions() + positions[upper_indices, 2] += dd + atoms.set_positions(positions) + + # Calculate energy (no relaxation!) + energy = atoms.get_potential_energy() + + # Calculate forces for stress + forces = atoms.get_forces() + + # Sum of z-forces on upper region + fz_upper = np.sum(forces[upper_indices, 2]) + + # Convert to stress (GPa): σ = F / A + # Negate because forces on upper atoms point downward (negative z) + # but traction (tensile stress) should be positive + sig_upper = -EV_PER_A3_TO_GPA * fz_upper / area + + separations.append(dd) + energies.append(energy) + traction.append(sig_upper) + + # Max traction from force-based calculation + max_traction = np.max(np.abs(traction)) + + return { + "separations": separations, + "energies": energies, + "traction": traction, + "max_traction": max_traction, + } + + +# ============================================================================= +# Helper Functions +# ============================================================================= + + +def _save_curve(write_dir: Path, name: str, data: dict[str, list]) -> None: + """ + Save curve data to CSV file. + + Parameters + ---------- + write_dir : Path + Directory to save the file. + name : str + Base name for the CSV file (without extension). + data : dict[str, list] + Column name to data mapping for the DataFrame. + """ + pd.DataFrame(data).to_csv(write_dir / f"{name}.csv", index=False) + + +# ============================================================================= +# Main Benchmark Function +# ============================================================================= + + +def run_iron_properties(model_name: str, model: Any) -> None: + """ + Run the full iron properties benchmark for a single model. + + This benchmark includes: + - Equation of state (lattice parameter, bulk modulus) + - Elastic constants (C11, C12, C44) + - Bain path energy curve + - Vacancy formation energy + - Surface energies (100, 110, 111, 112) + - Stacking fault energy curves (110, 112) + - Traction-separation curves (100, 110) + + Parameters + ---------- + model_name + Name of the model being evaluated. + model + Model wrapper providing ``get_calculator``. + """ + calc = model.get_calculator() + write_dir = OUT_PATH / model_name + write_dir.mkdir(parents=True, exist_ok=True) + + results: dict[str, Any] = {} + + # EOS calculation + print(f"[{model_name}] Running EOS calculation...") + eos_results = run_eos_calculation(calc) + results["eos"] = eos_results + a0 = eos_results["a0"] + print( + f"[{model_name}] Lattice parameter: {a0:.4f} Å, " + f"Bulk modulus: {eos_results['B0']:.1f} GPa" + ) + + # Save EOS curve data + _save_curve( + write_dir, + "eos_curve", + {"volume": eos_results["volumes"], "energy": eos_results["energies"]}, + ) + + # Elastic constants calculation + print(f"[{model_name}] Running elastic constants calculation...") + elastic_results = run_elastic_calculation(calc, a0) + results["elastic"] = elastic_results + print( + f"[{model_name}] C11={elastic_results['C11']:.1f}, " + f"C12={elastic_results['C12']:.1f}, C44={elastic_results['C44']:.1f} GPa" + ) + + # Bain path calculation + print(f"[{model_name}] Running Bain path calculation...") + bain_results = run_bain_path_calculation(calc, a0) + results["bain_path"] = bain_results + + # Save Bain path data + _save_curve( + write_dir, + "bain_path", + { + "ca_ratio": bain_results["ca_ratios"], + "energy": bain_results["energies"], + "energy_meV": bain_results["energies_meV"], + }, + ) + + # Vacancy calculation + print(f"[{model_name}] Running vacancy calculation...") + vacancy_results = run_vacancy_calculation(calc, a0) + results["vacancy"] = vacancy_results + print(f"[{model_name}] E_vac = {vacancy_results['E_vac']:.3f} eV") + + # Surface calculations + print(f"[{model_name}] Running surface calculations...") + surface_results = run_surface_calculations(calc, a0) + results["surfaces"] = surface_results + + # SFE calculations + sfe_results = {} + for sfe_type in ["110", "112"]: + print(f"[{model_name}] Running SFE {sfe_type} calculation...") + sfe_result = run_sfe_calculation(calc, a0, sfe_type) + sfe_results[sfe_type] = sfe_result + results[f"sfe_{sfe_type}"] = {"max_sfe": sfe_result["max_sfe"]} + _save_curve( + write_dir, + f"sfe_{sfe_type}_curve", + { + "displacement": sfe_result["displacements"], + "sfe_J_per_m2": sfe_result["sfe_J_per_m2"], + }, + ) + + # T-S calculations + ts_results = {} + for direction in ["100", "110"]: + print(f"[{model_name}] Running T-S ({direction}) calculation...") + ts_result = run_ts_calculation(calc, a0, direction) + ts_results[direction] = ts_result + results[f"ts_{direction}"] = {"max_traction": ts_result["max_traction"]} + _save_curve( + write_dir, + f"ts_{direction}_curve", + { + "separation": ts_result["separations"], + "energy": ts_result["energies"], + "traction": ts_result["traction"], + }, + ) + print( + f"[{model_name}] Max traction ({direction}): " + f"{ts_result['max_traction']:.2f} GPa" + ) + + # Save all results as JSON + (write_dir / "results.json").write_text(json.dumps(results, indent=2, default=str)) + + # Save summary metrics + summary: dict[str, Any] = { + "a0": a0, + "B0": eos_results["B0"], + "C11": elastic_results["C11"], + "C12": elastic_results["C12"], + "C44": elastic_results["C44"], + "E_bcc_fcc_meV": bain_results["delta_E_meV"], + "E_vac": vacancy_results["E_vac"], + "gamma_100": surface_results["gamma_100"], + "gamma_110": surface_results["gamma_110"], + "gamma_111": surface_results["gamma_111"], + "gamma_112": surface_results["gamma_112"], + "max_sfe_110": sfe_results["110"]["max_sfe"], + "max_sfe_112": sfe_results["112"]["max_sfe"], + "max_traction_100": ts_results["100"]["max_traction"], + "max_traction_110": ts_results["110"]["max_traction"], + } + + (write_dir / "summary.json").write_text(json.dumps(summary, indent=2)) + + print(f"[{model_name}] Done. Results saved to {write_dir}") + + +@pytest.mark.slow +@pytest.mark.parametrize("model_name", MODELS) +def test_iron_properties(model_name: str) -> None: + """ + Run iron properties benchmark for each registered model. + + This test is marked as slow and excluded from default test runs. + Run with ``pytest --run-slow`` to include. + + Parameters + ---------- + model_name + Name of the model to evaluate. + """ + run_iron_properties(model_name, MODELS[model_name]) diff --git a/ml_peg/calcs/physicality/iron_properties/iron_utils.py b/ml_peg/calcs/physicality/iron_properties/iron_utils.py new file mode 100644 index 000000000..fd1c15c58 --- /dev/null +++ b/ml_peg/calcs/physicality/iron_properties/iron_utils.py @@ -0,0 +1,697 @@ +""" +Utility functions for BCC iron property calculations. + +This module provides structure creation and EOS fitting functions for iron benchmarks. +""" + +from __future__ import annotations + +from typing import Any + +from ase import Atoms +from ase.build import bulk +import numpy as np +from scipy.optimize import leastsq, minimize_scalar + +# ============================================================================= +# Unit Conversion Constants +# ============================================================================= + +EV_TO_J = 1.60218e-19 +ANGSTROM_TO_M = 1.0e-10 +EV_PER_A2_TO_J_PER_M2 = 16.0217733 +EV_PER_A3_TO_GPA = 160.21765 + +# ============================================================================= +# Crystallographic Rotation Matrices +# ============================================================================= + +# Rotation matrix for [-110]/[111]/[11-2] crystallographic frame +# Used for (111) surface, (112) surface, and {110}<111> SFE +ROTATION_111_FRAME = np.array( + [ + [-1, 1, 0], # ex: [-110] + [1, 1, 1], # ey: [111] + [1, 1, -2], # ez: [11-2] + ], + dtype=float, +) / np.array([[np.sqrt(2)], [np.sqrt(3)], [np.sqrt(6)]]) + + +# ============================================================================= +# EOS Fitting Functions +# ============================================================================= + + +def eos_birch_murnaghan( + params: tuple[float, float, float, float], vol: np.ndarray +) -> np.ndarray: + """ + Birch-Murnaghan equation of state (3rd order). + + Parameters + ---------- + params + (E0, B0, Bp, V0). + vol + Volume array. + + Returns + ------- + np.ndarray + Energy array. + """ + E0, B0, Bp, V0 = params # noqa: N806 + eta = (vol / V0) ** (1.0 / 3.0) + return E0 + 9.0 * B0 * V0 / 16.0 * (eta**2 - 1.0) ** 2 * ( + 6.0 + Bp * (eta**2 - 1.0) - 4.0 * eta**2 + ) + + +def get_eos_initial_guess( + vol: np.ndarray, ene: np.ndarray +) -> tuple[float, float, float, float]: + """ + Get initial guess for EOS parameters using quadratic fit. + + Parameters + ---------- + vol + Volume array. + ene + Energy array. + + Returns + ------- + tuple + (E0, B0, Bp, V0) initial guess. + """ + a, b, c = np.polyfit(vol, ene, 2) + V0 = -b / (2 * a) # noqa: N806 + E0 = a * V0**2 + b * V0 + c # noqa: N806 + B0 = 2 * a * V0 # noqa: N806 + Bp = 4.0 # noqa: N806 + return E0, B0, Bp, V0 + + +def fit_eos( + vol: np.ndarray, + ene: np.ndarray, +) -> dict[str, Any]: + """ + Fit Birch-Murnaghan equation of state to energy-volume data. + + Parameters + ---------- + vol + Volume per atom array (Angstrom^3). + ene + Energy per atom array (eV). + + Returns + ------- + dict + Fitted parameters: + - E0: Equilibrium energy (eV) + - B0: Bulk modulus (GPa) + - Bp: Pressure derivative (dimensionless) + - V0: Equilibrium volume per atom (Angstrom^3) + - a0: Equilibrium lattice parameter (Angstrom) - for BCC + """ + x0 = get_eos_initial_guess(vol, ene) + + def residual(params, y, x): + """ + Compute residual for EOS fitting. + + Parameters + ---------- + params : tuple + EOS parameters (E0, B0, Bp, V0). + y : np.ndarray + Observed energies. + x : np.ndarray + Volumes. + + Returns + ------- + np.ndarray + Residual array (observed - predicted). + """ + return y - eos_birch_murnaghan(params, x) + + params, _ = leastsq(residual, x0, args=(ene, vol)) + E0, B0, Bp, V0 = params # noqa: N806 + + # Convert bulk modulus to GPa (from eV/Angstrom^3) + B0_GPa = B0 * EV_PER_A3_TO_GPA # noqa: N806 + + # Calculate lattice parameter for BCC (2 atoms per unit cell) + a0 = (V0 * 2) ** (1.0 / 3.0) + + return {"E0": E0, "B0": B0_GPa, "Bp": Bp, "V0": V0, "a0": a0} + + +# ============================================================================= +# Isotropic Volume Relaxation +# ============================================================================= + + +def relax_volume_isotropic( + atoms: Atoms, + calc: Any, + scale_bounds: tuple[float, float] = (0.9, 1.1), + xtol: float = 1e-8, +) -> Atoms: + """ + Relax cell volume isotropically (uniform scaling) to minimize energy. + + This maintains cell shape (all ratios between cell dimensions) while finding + the optimal volume. This is equivalent to LAMMPS 'fix box/relax aniso 0.0 + couple xyz' which couples all three diagonal stress components together, + allowing only uniform scaling during relaxation. + + For a tetragonal cell with c/a ratio, this preserves c/a while optimizing + the volume. + + Parameters + ---------- + atoms : Atoms + ASE Atoms object (will be copied, not modified). + calc : Any + ASE calculator. + scale_bounds : tuple[float, float], optional + Bounds for the scale factor search (min, max). Default: (0.9, 1.1). + xtol : float, optional + Tolerance for the scale factor optimization. Default: 1e-8. + + Returns + ------- + Atoms + New Atoms object at optimal volume with same cell shape. + + Notes + ----- + This function matches the LAMMPS behavior for Bain path calculations where + 'couple xyz' is used to maintain the c/a ratio during volume relaxation. + The optimization finds the uniform scale factor that minimizes the total + energy of the system. + """ + atoms = atoms.copy() + original_cell = atoms.cell.array.copy() + + def energy_at_scale(scale: float) -> float: + """ + Calculate energy at a given uniform scale factor. + + Parameters + ---------- + scale : float + Uniform scale factor to apply to the cell. + + Returns + ------- + float + Potential energy of the system at the given scale. + """ + test_atoms = atoms.copy() + test_atoms.set_cell(original_cell * scale, scale_atoms=True) + test_atoms.calc = calc + return test_atoms.get_potential_energy() + + # Find optimal scale factor that minimizes energy + result = minimize_scalar( + energy_at_scale, + bounds=scale_bounds, + method="bounded", + options={"xatol": xtol}, + ) + optimal_scale = result.x + + # Create relaxed structure at optimal volume + relaxed_atoms = atoms.copy() + relaxed_atoms.set_cell(original_cell * optimal_scale, scale_atoms=True) + relaxed_atoms.calc = calc + + return relaxed_atoms + + +# ============================================================================= +# Structure Creation Functions +# ============================================================================= + + +def _create_oriented_bcc_structure( + lattice_parameter: float, + rotation: np.ndarray, + cell_dims: tuple[float, float, float], + max_range: int, + symbol: str = "Fe", + wrap: bool = True, +) -> Atoms: + """ + Create BCC structure with given orientation using rotation matrix. + + This is a generic function used by several oriented structure creators. + + Parameters + ---------- + lattice_parameter : float + BCC lattice parameter in Angstroms. + rotation : np.ndarray + 3x3 rotation matrix (rows are the new basis vectors). + cell_dims : tuple[float, float, float] + Cell dimensions (lx, ly, lz) in Angstroms. + max_range : int + Range for scanning cubic positions. + symbol : str, optional + Chemical symbol (default: 'Fe'). + wrap : bool, optional + Whether to wrap positions into cell (default: True). + + Returns + ------- + Atoms + ASE Atoms object with oriented structure. + """ + a = lattice_parameter + lx, ly, lz = cell_dims + cell = np.array([[lx, 0, 0], [0, ly, 0], [0, 0, lz]]) + + positions = [] + eps = 1e-8 + + for i in range(-max_range, max_range + 1): + for j in range(-max_range, max_range + 1): + for k in range(-max_range, max_range + 1): + for basis in [(0, 0, 0), (0.5, 0.5, 0.5)]: + pos_cubic = a * np.array([i + basis[0], j + basis[1], k + basis[2]]) + pos_oriented = rotation @ pos_cubic + + frac_x = pos_oriented[0] / lx + frac_y = pos_oriented[1] / ly + frac_z = pos_oriented[2] / lz + + if ( + 0 - eps <= frac_x < 1 - eps + and 0 - eps <= frac_y < 1 - eps + and 0 - eps <= frac_z < 1 - eps + ): + positions.append(pos_oriented) + + if len(positions) == 0: + raise ValueError("No atoms found for oriented structure") + + positions = np.array(positions) + _, unique_idx = np.unique( + np.round(positions, decimals=6), axis=0, return_index=True + ) + positions = positions[unique_idx] + + atoms = Atoms( + symbols=[symbol] * len(positions), positions=positions, cell=cell, pbc=True + ) + + if wrap: + atoms.wrap() + + return atoms + + +def create_bcc_supercell( + lattice_parameter: float, size: tuple = (4, 4, 4), symbol: str = "Fe" +) -> Atoms: + """ + Create a BCC supercell. + + Parameters + ---------- + lattice_parameter + Lattice parameter in Angstroms. + size + Supercell size as (nx, ny, nz). + symbol + Chemical symbol (default: 'Fe'). + + Returns + ------- + Atoms + ASE Atoms object. + """ + unit_cell = bulk(symbol, "bcc", a=lattice_parameter, cubic=True) + return unit_cell * size + + +def create_bain_cell(lattice_parameter: float, ca_ratio: float) -> Atoms: + """ + Create a tetragonally distorted BCC cell for Bain path calculation. + + Parameters + ---------- + lattice_parameter + BCC lattice parameter. + ca_ratio + Target c/a ratio. + + Returns + ------- + Atoms + Tetragonally distorted cell. + """ + beta = (1.0 / ca_ratio) ** (1.0 / 3.0) + al = lattice_parameter * beta + alz = al * ca_ratio + + cell = np.array([[al, 0, 0], [0, al, 0], [0, 0, alz]]) + positions = np.array([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]]) @ cell + + return Atoms(symbols=["Fe", "Fe"], positions=positions, cell=cell, pbc=True) + + +def create_surface_100( + lattice_parameter: float, layers: int = 10, vacuum: float = 0.0, symbol: str = "Fe" +) -> Atoms: + """ + Create a (100) surface slab for BCC iron. + + Parameters + ---------- + lattice_parameter : float + Lattice parameter in Angstroms. + layers : int, optional + Number of atomic layers (default: 10). + vacuum : float, optional + Vacuum thickness in Angstroms (default: 0.0). + symbol : str, optional + Chemical symbol (default: 'Fe'). + + Returns + ------- + Atoms + ASE Atoms object with the (100) surface slab. + """ + a = lattice_parameter + cell = np.array([[a, 0, 0], [0, a, 0], [0, 0, a * layers]]) + + positions = [] + for k in range(layers): + positions.append([0, 0, k * a]) + positions.append([0.5 * a, 0.5 * a, (k + 0.5) * a]) + + atoms = Atoms( + symbols=[symbol] * len(positions), positions=positions, cell=cell, pbc=True + ) + if vacuum > 0: + atoms.center(vacuum=vacuum, axis=2) + return atoms + + +def create_surface_110( + lattice_parameter: float, layers: int = 10, vacuum: float = 0.0, symbol: str = "Fe" +) -> Atoms: + """ + Create a (110) surface slab for BCC iron. + + Parameters + ---------- + lattice_parameter : float + Lattice parameter in Angstroms. + layers : int, optional + Number of atomic layers (default: 10). + vacuum : float, optional + Vacuum thickness in Angstroms (default: 0.0). + symbol : str, optional + Chemical symbol (default: 'Fe'). + + Returns + ------- + Atoms + ASE Atoms object with the (110) surface slab. + """ + a = lattice_parameter + lx = a + ly = a * np.sqrt(2) + lz = a * np.sqrt(2) * layers + + cell = np.array([[lx, 0, 0], [0, ly, 0], [0, 0, lz]]) + positions = [] + d110 = a * np.sqrt(2) / 2 + + # Each (110) plane in a cell of a x a*sqrt(2) contains 2 atoms. + # The BCC (110) surface has a centered rectangular structure. + for k in range(layers * 2): + z = k * d110 + if k % 2 == 0: + # Even planes: atoms at (0, 0) and (a/2, ly/2) + positions.append([0, 0, z]) + positions.append([0.5 * a, 0.5 * ly, z]) + else: + # Odd planes: atoms at (0, ly/2) and (a/2, 0) + positions.append([0, 0.5 * ly, z]) + positions.append([0.5 * a, 0, z]) + + atoms = Atoms( + symbols=[symbol] * len(positions), positions=positions, cell=cell, pbc=True + ) + if vacuum > 0: + atoms.center(vacuum=vacuum, axis=2) + return atoms + + +def create_surface_111( + lattice_parameter: float, + size: tuple = (3, 15, 3), + vacuum: float = 0.0, + symbol: str = "Fe", +) -> Atoms: + """ + Create a (111) surface slab for BCC iron. + + Parameters + ---------- + lattice_parameter : float + Lattice parameter in Angstroms. + size : tuple, optional + Cell size as (nx, ny, nz) (default: (3, 15, 3)). + vacuum : float, optional + Vacuum thickness in Angstroms (default: 0.0). + symbol : str, optional + Chemical symbol (default: 'Fe'). + + Returns + ------- + Atoms + ASE Atoms object with the (111) surface slab. + """ + a = lattice_parameter + + cell_dims = ( + a * np.sqrt(2) * size[0], + a * np.sqrt(3) * size[1], + a * np.sqrt(6) * size[2], + ) + max_range = int(max(size) * 3 + 5) + + atoms = _create_oriented_bcc_structure( + lattice_parameter, ROTATION_111_FRAME, cell_dims, max_range, symbol + ) + + if vacuum > 0: + atoms.center(vacuum=vacuum, axis=1) + + return atoms + + +def create_surface_112( + lattice_parameter: float, layers: int = 15, vacuum: float = 0.0, symbol: str = "Fe" +) -> Atoms: + """ + Create a (112) surface slab for BCC iron. + + Parameters + ---------- + lattice_parameter : float + Lattice parameter in Angstroms. + layers : int, optional + Number of atomic layers (default: 15). + vacuum : float, optional + Vacuum thickness in Angstroms (default: 0.0). + symbol : str, optional + Chemical symbol (default: 'Fe'). + + Returns + ------- + Atoms + ASE Atoms object with the (112) surface slab. + """ + a = lattice_parameter + + cell_dims = (a * np.sqrt(2), a * np.sqrt(3), a * np.sqrt(6) * layers) + max_range = int(layers * 3 + 5) + + atoms = _create_oriented_bcc_structure( + lattice_parameter, ROTATION_111_FRAME, cell_dims, max_range, symbol + ) + + if vacuum > 0: + atoms.center(vacuum=vacuum, axis=2) + + return atoms + + +def create_sfe_110_structure(lattice_parameter: float) -> Atoms: + """ + Create structure for {110}<111> stacking fault calculation. + + Parameters + ---------- + lattice_parameter : float + Lattice parameter in Angstroms. + + Returns + ------- + Atoms + ASE Atoms object for SFE calculation. + """ + a = lattice_parameter + size = (20, 1, 3) + + cell_dims = ( + a * np.sqrt(2) * size[0], + a * np.sqrt(3) * size[1], + a * np.sqrt(6) * size[2], + ) + max_range = int(max(size) * 3 + 5) + + return _create_oriented_bcc_structure( + lattice_parameter, ROTATION_111_FRAME, cell_dims, max_range + ) + + +def create_sfe_112_structure(lattice_parameter: float) -> Atoms: + """ + Create structure for {112}<111> stacking fault calculation. + + Parameters + ---------- + lattice_parameter : float + Lattice parameter in Angstroms. + + Returns + ------- + Atoms + ASE Atoms object for SFE calculation. + """ + a = lattice_parameter + size = (15, 1, 1) + + # Rotation matrix for {112} orientation + ex = np.array([1, 1, -2]) / np.sqrt(6) + ey = np.array([-1, 1, 0]) / np.sqrt(2) + ez = np.array([1, 1, 1]) / np.sqrt(3) + rotation = np.array([ex, ey, ez]) + + cell_dims = ( + a * np.sqrt(6) * size[0], + a * np.sqrt(2) * size[1], + a * np.sqrt(3) * size[2], + ) + max_range = int(max(size) * 3 + 5) + + return _create_oriented_bcc_structure( + lattice_parameter, rotation, cell_dims, max_range + ) + + +# ============================================================================= +# Elastic Calculation Utilities +# ============================================================================= + + +def apply_voigt_strain(atoms: Atoms, direction: int, magnitude: float) -> Atoms: + """ + Apply Voigt strain with off-diagonal cell adjustment. + + For normal strains (directions 1-3), this scales the entire cell vector + rather than just the diagonal component. This maintains cell vector ratios + and is important for triclinic cells or pre-strained configurations. + + LAMMPS equivalent for direction 1 (xx): + change_box all x delta 0 ${delta} xy delta ${deltaxy} xz delta ${deltaxz} + where deltaxy = up * xy, deltaxz = up * xz + + For cubic/orthogonal cells (xy=xz=yz=0), this is equivalent to apply_strain(). + + Parameters + ---------- + atoms : Atoms + ASE Atoms object. + direction : int + Voigt direction (1-6): + 1=xx, 2=yy, 3=zz, 4=yz, 5=xz, 6=xy. + magnitude : float + Strain magnitude (e.g., 1e-5). + + Returns + ------- + Atoms + Strained ASE Atoms object. + """ + atoms_strained = atoms.copy() + cell = atoms_strained.cell.array.copy() + + if direction == 1: + # Scale entire x cell vector (a1): maintains xy/lx and xz/lx ratios + # LAMMPS: x -> x + delta, xy -> xy + up*xy, xz -> xz + up*xz + cell[0, :] *= 1 + magnitude + elif direction == 2: + # Scale entire y cell vector (a2): maintains yz/ly ratio + # LAMMPS: y -> y + delta, yz -> yz + up*yz + cell[1, :] *= 1 + magnitude + elif direction == 3: + # Scale entire z cell vector (a3) + # LAMMPS: z -> z + delta + cell[2, :] *= 1 + magnitude + elif direction == 4: + # yz shear: LAMMPS changes yz tilt only + # For LAMMPS compatibility: simple shear (not symmetric) + # cell[1, 2] is the yz tilt component + lz = cell[2, 2] + cell[1, 2] += magnitude * lz + elif direction == 5: + # xz shear: LAMMPS changes xz tilt only + lz = cell[2, 2] + cell[0, 2] += magnitude * lz + elif direction == 6: + # xy shear: LAMMPS changes xy tilt only + ly = cell[1, 1] + cell[0, 1] += magnitude * ly + + atoms_strained.set_cell(cell, scale_atoms=True) + return atoms_strained + + +def calculate_surface_energy( + E_slab: float, # noqa: N803 + E_bulk: float, # noqa: N803 + area: float, +) -> float: + """ + Calculate surface energy in J/m^2. + + Parameters + ---------- + E_slab : float + Total energy of the slab with vacuum (eV). + E_bulk : float + Total energy of the bulk reference (eV). + area : float + Surface area (Angstrom^2). + + Returns + ------- + float + Surface energy in J/m^2. + """ + delta_E = E_slab - E_bulk # noqa: N806 + return delta_E * EV_TO_J / (2 * area * ANGSTROM_TO_M**2)