Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 20 additions & 13 deletions amorphouspy/src/amorphouspy/lammps/potentials/du_teter_potential.py
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@

import itertools
from pathlib import Path
from typing import Literal

import numpy as np
import pandas as pd
Expand Down Expand Up @@ -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.
Expand All @@ -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
Expand Down Expand Up @@ -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:
Expand All @@ -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.

Expand All @@ -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.
Expand All @@ -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
Expand Down
17 changes: 14 additions & 3 deletions amorphouspy/src/amorphouspy/lammps/potentials/potential.py
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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)
Expand Down
35 changes: 14 additions & 21 deletions amorphouspy/src/tests/test_potentials.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
# ---------------------------------------------------------------------------
Expand Down Expand Up @@ -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"}],
Expand All @@ -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},
Expand All @@ -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"])

Expand Down Expand Up @@ -1150,33 +1143,33 @@ 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 = {
"atoms": [{"element": "B"}, {"element": "O"}, {"element": "Ca"}],
"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])


Expand Down
Loading