diff --git a/amorphouspy/src/amorphouspy/lammps/potentials/du_teter_potential.py b/amorphouspy/src/amorphouspy/lammps/potentials/du_teter_potential.py index 8a6d6abe..e393585c 100644 --- a/amorphouspy/src/amorphouspy/lammps/potentials/du_teter_potential.py +++ b/amorphouspy/src/amorphouspy/lammps/potentials/du_teter_potential.py @@ -28,6 +28,7 @@ import itertools from pathlib import Path +from typing import Literal import numpy as np import pandas as pd @@ -601,13 +602,19 @@ def fit_BO_params(K: float, R: float, N4: float | None = None) -> dict: } -def get_all_BO_params(structure_dict: dict, *, original_dbx_approach: bool = True) -> dict: +def get_all_BO_params( + structure_dict: dict, + *, + n4_model: Literal["dbx", "dbx_generalized"] = "dbx", +) -> dict: """Compute all B-O interaction parameters for the given composition. Args: structure_dict: Dictionary containing the structure information, - original_dbx_approach: If True, use the original DBX approach with R = cNa2O / cB2O3. - If False, use the modified approach with R = (sum of all modifiers - Al2O3) / cB2O3. + n4_model: Model used to estimate the fraction of 4-coordinated boron (N4). + ``"dbx"`` — original Dell-Bray-Xiao model with R = cNa2O / cB2O3. + ``"dbx_generalized"`` — generalized approach with + R = (sum of all modifiers - Al2O3) / cB2O3. Returns: Dictionary with Du/Teter parameters for the B-O pair. @@ -632,7 +639,7 @@ def get_all_BO_params(structure_dict: dict, *, original_dbx_approach: bool = Tru cBaO = mol_fraction.get("BaO", 0) cAl2O3 = mol_fraction.get("Al2O3", 0) - if original_dbx_approach: + if n4_model == "dbx": R = max(cNa2O / cB2O3, 0) else: modifier_sum = cLi2O + cNa2O + cK2O + cBeO + cMgO + cCaO + cSrO + cBaO - cAl2O3 @@ -708,12 +715,12 @@ def _build_all_pair_params( species: list[str], structure_dict: dict, *, - original_dbx_approach: bool = True, + n4_model: Literal["dbx", "dbx_generalized"] = "dbx", ) -> dict[str, dict]: """Build per-pair Du/Teter parameter dicts for all X-O interactions.""" pair_params: dict[str, dict] = {} if "B" in species: - pair_params["B-O"] = get_all_BO_params(structure_dict, original_dbx_approach=original_dbx_approach) + pair_params["B-O"] = get_all_BO_params(structure_dict, n4_model=n4_model) for elem in species: pair_name = "O-O" if elem == "O" else f"{elem}-O" if pair_name not in pair_params: @@ -740,7 +747,7 @@ def generate_du_teter_potential( *, melt: bool = True, use_three_body: bool = False, - original_dbx_approach: bool = True, + n4_model: Literal["dbx", "dbx_generalized"] = "dbx", ) -> pd.DataFrame: """Generate a LAMMPS potential file for the Du/Teter potential. @@ -756,11 +763,11 @@ def generate_du_teter_potential( use_three_body: If True, write a Stillinger-Weber ``.sw`` file and add the ``sw`` pair style to the LAMMPS config for O-P-O / P-O-P three-body interactions. Requires P to be present in the structure. - original_dbx_approach: If True (default), use the original Dell-Bray-Xiao - model with R = cNa2O / cB2O3. If False, use a custom, generalized approach - (unpublished) with: - R = (sum of all alkaline and earth-alkaline oxides - Al2O3) / cB2O3 - is used. + n4_model: Model used to estimate the fraction of 4-coordinated boron (N4). + ``"dbx"`` (default) — original Dell-Bray-Xiao model with + R = cNa2O / cB2O3. + ``"dbx_generalized"`` — generalized approach (unpublished) with + R = (sum of all alkaline and earth-alkaline oxides - Al2O3) / cB2O3. Returns: A DataFrame containing the potential configuration. @@ -773,7 +780,7 @@ def generate_du_teter_potential( _validate_du_teter_inputs(species, use_three_body=use_three_body) - pair_params = _build_all_pair_params(species, structure_dict, original_dbx_approach=original_dbx_approach) + pair_params = _build_all_pair_params(species, structure_dict, n4_model=n4_model) # ------------------------------------------------------------------ # Write table files diff --git a/amorphouspy/src/amorphouspy/lammps/potentials/potential.py b/amorphouspy/src/amorphouspy/lammps/potentials/potential.py index 1998ba00..36004b28 100644 --- a/amorphouspy/src/amorphouspy/lammps/potentials/potential.py +++ b/amorphouspy/src/amorphouspy/lammps/potentials/potential.py @@ -17,13 +17,23 @@ __all__ = ["DsfConfig", "EwaldConfig", "InteractionConfig", "PppmConfig", "WolfConfig"] # Preference order: pmmcs covers the most elements, shik adds B, bjp is most limited. -POTENTIAL_PREFERENCE = ("pmmcs", "bmp-harmonic", "bmp-screened-harmonic", "shik", "bjp", "du_teter", "yang2026") +POTENTIAL_PREFERENCE = ( + "pmmcs", + "bmp-harmonic", + "bmp-screened-harmonic", + "shik", + "bjp", + "du_teter", + "du_teter_dbx_generalized", + "yang2026", +) _POTENTIAL_MODULES = { "pmmcs": pmmcs, "bjp": bjp, "shik": shik, "du_teter": du_teter, + "du_teter_dbx_generalized": du_teter, "bmp-harmonic": bmp, "bmp-screened-harmonic": bmp, "yang2026": yang2026, @@ -122,8 +132,9 @@ def generate_potential( return bjp.generate_bjp_potential(atoms_dict, melt=melt, electrostatics=electrostatics) if potential_type.lower() == "shik": return shik.generate_shik_potential(atoms_dict, melt=melt, electrostatics=electrostatics) - if potential_type.lower() == "du_teter": - return du_teter.generate_du_teter_potential(atoms_dict, melt=melt) + if potential_type.lower() in ("du_teter", "du_teter_dbx_generalized"): + n4_model = "dbx_generalized" if potential_type.lower() == "du_teter_dbx_generalized" else "dbx" + return du_teter.generate_du_teter_potential(atoms_dict, melt=melt, n4_model=n4_model) if potential_type.lower() in ("bmp-harmonic", "bmp-screened-harmonic"): variant = "harmonic" if potential_type.lower() == "bmp-harmonic" else "screened-harmonic" return bmp.generate_bmp_potential(atoms_dict, variant=variant, melt=melt, electrostatics=electrostatics) diff --git a/amorphouspy/src/tests/test_potentials.py b/amorphouspy/src/tests/test_potentials.py index a7a72596..7bf6a192 100644 --- a/amorphouspy/src/tests/test_potentials.py +++ b/amorphouspy/src/tests/test_potentials.py @@ -146,13 +146,6 @@ def test_compatible_potentials_preserves_preference_order(): assert indices == sorted(indices) -def test_compatible_potentials_all_subset_of_known(): - """Results contain only recognised potential names.""" - result = compatible_potentials({"Si", "O"}) - available_potentials = ["pmmcs", "bmp-harmonic", "bmp-screened-harmonic", "shik", "bjp", "du_teter", "yang2026"] - assert all(p in available_potentials for p in result) - - # --------------------------------------------------------------------------- # Du/Teter SW three-body # --------------------------------------------------------------------------- @@ -1010,7 +1003,7 @@ def test_get_all_BO_params_default_mol_fractions(): def test_get_all_BO_params_original_approach_uses_only_Na2O(): - """original_dbx_approach=True uses R = cNa2O / cB2O3 (ignores other modifiers).""" + """n4_model='dbx' uses R = cNa2O / cB2O3 (ignores other modifiers).""" # Same composition but different modifier: Ca instead of Na struct_na = { "atoms": [{"element": "B"}, {"element": "O"}, {"element": "Na"}], @@ -1020,15 +1013,15 @@ def test_get_all_BO_params_original_approach_uses_only_Na2O(): "atoms": [{"element": "B"}, {"element": "O"}, {"element": "Ca"}], "mol_fraction": {"B2O3": 0.3, "Na2O": 0.0, "CaO": 0.2}, } - result_na = get_all_BO_params(struct_na, original_dbx_approach=True) - result_ca = get_all_BO_params(struct_ca, original_dbx_approach=True) + result_na = get_all_BO_params(struct_na, n4_model="dbx") + result_ca = get_all_BO_params(struct_ca, n4_model="dbx") # Ca contributes nothing to R in the original approach (only Na2O counted) # → result_ca should have R=0, result_na has R=0.2/0.3 > 0, giving different A assert result_na["A"] != result_ca["A"] def test_get_all_BO_params_modified_approach_includes_all_modifiers(): - """original_dbx_approach=False sums all modifiers (Li2O, Na2O, K2O, MgO, CaO, SrO, BaO, BeO minus Al2O3).""" + """n4_model='dbx_generalized' sums all modifiers (Li2O, Na2O, K2O, MgO, CaO, SrO, BaO, BeO minus Al2O3).""" struct_na = { "atoms": [{"element": "B"}, {"element": "O"}, {"element": "Na"}], "mol_fraction": {"B2O3": 0.3, "Na2O": 0.2, "CaO": 0.0}, @@ -1037,17 +1030,17 @@ def test_get_all_BO_params_modified_approach_includes_all_modifiers(): "atoms": [{"element": "B"}, {"element": "O"}, {"element": "Ca"}], "mol_fraction": {"B2O3": 0.3, "Na2O": 0.0, "CaO": 0.2}, } - result_na = get_all_BO_params(struct_na, original_dbx_approach=False) - result_ca = get_all_BO_params(struct_ca, original_dbx_approach=False) + result_na = get_all_BO_params(struct_na, n4_model="dbx_generalized") + result_ca = get_all_BO_params(struct_ca, n4_model="dbx_generalized") # Both have the same total modifier / B2O3 ratio → identical parameters assert result_na["A"] == pytest.approx(result_ca["A"], rel=1e-6) assert result_na["rc"] == pytest.approx(result_ca["rc"], rel=1e-4) def test_get_all_BO_params_original_is_default(): - """original_dbx_approach defaults to True.""" + """n4_model defaults to 'dbx'.""" result_default = get_all_BO_params(_NABS_STRUCTURE) - result_explicit = get_all_BO_params(_NABS_STRUCTURE, original_dbx_approach=True) + result_explicit = get_all_BO_params(_NABS_STRUCTURE, n4_model="dbx") assert result_default["A"] == pytest.approx(result_explicit["A"]) assert result_default["rc"] == pytest.approx(result_explicit["rc"]) @@ -1150,20 +1143,20 @@ def test_generate_du_teter_with_boron(tmp_path): def test_generate_du_teter_original_dbx_approach_default(tmp_path): - """original_dbx_approach defaults to True (Na2O-only R).""" + """n4_model defaults to 'dbx' (Na2O-only R).""" atoms = { "atoms": [{"element": "B"}, {"element": "O"}, {"element": "Na"}], "mol_fraction": {"B2O3": 0.3, "Na2O": 0.2}, } generate_du_teter_potential(atoms, output_dir=str(tmp_path / "default"), melt=False) - generate_du_teter_potential(atoms, output_dir=str(tmp_path / "explicit"), melt=False, original_dbx_approach=True) + generate_du_teter_potential(atoms, output_dir=str(tmp_path / "explicit"), melt=False, n4_model="dbx") # Both should produce a table file for B-O assert list((tmp_path / "default").glob("table_B_O*")) assert list((tmp_path / "explicit").glob("table_B_O*")) def test_generate_du_teter_original_vs_modified_dbx_approach_differ(tmp_path): - """original_dbx_approach=False (all modifiers) changes B-O parameters vs True (Na only).""" + """n4_model='dbx_generalized' (all modifiers) changes B-O parameters vs 'dbx' (Na only).""" # Na-B-O: original approach counts Na2O; modified approach also counts Na2O — same here # Use CaO only so original (Na-only) gives R=0, modified gives R=cCaO/cB2O3 > 0 → different A atoms = { @@ -1171,12 +1164,12 @@ def test_generate_du_teter_original_vs_modified_dbx_approach_differ(tmp_path): "mol_fraction": {"B2O3": 0.3, "CaO": 0.2}, } - result_orig = get_all_BO_params(atoms, original_dbx_approach=True) # R = cNa2O / cB2O3 = 0 - result_mod = get_all_BO_params(atoms, original_dbx_approach=False) # R = cCaO / cB2O3 > 0 + result_orig = get_all_BO_params(atoms, n4_model="dbx") # R = cNa2O / cB2O3 = 0 + result_mod = get_all_BO_params(atoms, n4_model="dbx_generalized") # R = cCaO / cB2O3 > 0 assert result_orig["A"] != result_mod["A"] # generate_du_teter_potential should accept and pass through the flag - df = generate_du_teter_potential(atoms, output_dir=str(tmp_path), melt=False, original_dbx_approach=False) + df = generate_du_teter_potential(atoms, output_dir=str(tmp_path), melt=False, n4_model="dbx_generalized") assert "pair_coeff" in "".join(df["Config"].iloc[0])