Melt-Quench Tutorial¶

This notebook walks through the full melt-quench workflow in amorphouspy: from composition input and interatomic potential selection through simulation setup, execution, and results inspection.

Sections:

Interatomic potentials — theory, supported elements, charges, and when to use each
Automatic potential selection
Structure generation
Potential generation
Melt-quench simulation — protocols and parameters
Post-processing and plotting
Quick reference

In [1]:

Copied!





from collections import Counter
import matplotlib.pyplot as plt
import numpy as np

import amorphouspy as am
from amorphouspy.lammps.potentials.potential import (
    compatible_potentials,
    generate_potential,
    get_supported_elements,
    select_potential,
)

import plotly.io as pio

pio.renderers.default = "notebook"
from collections import Counter
import matplotlib.pyplot as plt
import numpy as np

import amorphouspy as am
from amorphouspy.lammps.potentials.potential import (
    compatible_potentials,
    generate_potential,
    get_supported_elements,
    select_potential,
)

import plotly.io as pio

pio.renderers.default = "notebook"

1. Interatomic Potentials¶

amorphouspy ships four force-field models for oxide glass simulations. All use a Damped Shifted Force (DSF) Coulomb treatment for long-range electrostatics by default and differ in the short-range pair interaction.

Potential	Pair style	Elements	Default T_melt
PMMCS (Pedone)	Morse + Coulomb/DSF	29 elements (see below)	5000 K
BMP-harm (Bertani-Menziani-Pedone)	Morse + Buckingham + nb3b/harmonic + Coulomb/DSF	38 elements; B restricted	4000 K
BMP-shrm (Bertani-Menziani-Pedone)	Morse + Buckingham + nb3b/screened + Coulomb/DSF	38 elements; B restricted	4000 K
SHIK (Sundararaman)	Buckingham + r⁻²⁴ + Coulomb/DSF	9 elements	4000 K
BJP (Bouhadja)	Born-Mayer-Huggins + Coulomb/DSF	Ca, Al, Si, O	5000 K

Selection rule applied by select_potential: PMMCS → BMP-shrm → BMP-harm → SHIK → BJP (first that covers all elements). When B is present the order shifts to PMMCS → BMP-shrm → BMP-harm → SHIK → BJP.

1.1 PMMCS Potential (Pedone et al.)¶

References

Pedone et al., J. Phys. Chem. B 110, 11780 (2006) — original parameterisation

Functional form — Morse + repulsive term + Coulomb DSF:

$$V_{ij}(r) = D_{ij}\left[e^{-2a_{ij}(r - r_0)} - 2e^{-a_{ij}(r - r_0)}\right] + \frac{C_{ij}}{r^{12}} + \frac{q_i q_j}{r}$$

LAMMPS pair style: hybrid/overlay coul/dsf 0.25 8.0 pedone 5.5

Partial charges (formal charges scaled to 0.5×):

Element	Charge e
O	−1.2
Li, Na, K, Cu, Ag	+0.6
Be, Mg, Ca, Sr, Ba, Mn, Fe, Co, Ni, Zn	+1.2
Sc, Cr, Al, Fe3, Nd, Gd, Er	+1.8
Ti, Zr, Si, Ge, Sn	+2.4
P	+3.0

Supported elements (29):

Group	Elements
Alkali	Li, Na, K
Alkaline-earth	Be, Mg, Ca, Sr, Ba
Transition metals	Sc, Ti, Cr, Mn, Fe (Fe²⁺ / Fe³⁺), Co, Ni, Cu, Zn, Zr, Ag
Network formers	Al, Si, Ge, Sn, P
Rare-earth	Nd, Gd, Er
Anion	O

When to use: Broadest element coverage. The default choice for silicate, borate, phosphate, and mixed-modifier glasses. Suitable for most compositions without transition metals requiring charge-transfer.

Limitation: Does not include B.

1.2 SHIK Potential (Sundararaman et al.)¶

References

Sundararaman et al., J. Chem. Phys. 148, 194504 (2018) — SiO₂
Sundararaman et al., J. Chem. Phys. 150, 154505 (2019) — alkali and alkaline-earth aluminosilicate
Sundararaman et al., J. Chem. Phys. 152, 104501 (2020) — borate with mixed network formers
Shih et al., J. Non-Cryst. Sol. 565, 120853 (2021) — alkaline-earth silicate and borate

Functional form — Buckingham + r⁻²⁴ + Coulomb DSF:

$$V_{ij}(r) = A_{ij}\,e^{-B_{ij} r} - \frac{C_{ij}}{r^6} + \frac{D_{ij}}{r^{24}} + \frac{q_i q_j}{r}$$

The r⁻²⁴ term acts as a steep repulsive wall that improves stability at very short distances compared to a standard Buckingham potential.

LAMMPS pair style: hybrid/overlay coul/dsf 0.2 10.0 table spline 10000

Note: SHIK uses tabulated pair coefficients. generate_shik_potential() writes one .tbl file per pair to disk and embeds absolute paths in the LAMMPS input. The output_dir argument controls where these files land.

Partial charges (non-integer, composition-dependent):

Element	Charge e
Li	+0.5727
Na	+0.6018
K	+0.6849
Mg	+1.0850
Ca	+1.4977
B	+1.6126
Al	+1.6334
Si	+1.7755
O	computed for charge neutrality

Supported elements (9): Li, Na, K, Mg, Ca, B, Al, Si, O

When to use: Preferred for borosilicate, boroaluminosilicate, and alkali-silicate glasses, where B must be included. Highly validated for SiO₂. More limited element set than PMMCS.

Limitation: Only 9 elements; SHIK protocol initialises with a Langevin + NVE/limit pre-melt block at 5000 K (included automatically when melt=True, stripped in later protocol stages).

Special SHIK protocol note: Default temperature_high = 4000 K (lower than the others because the potential was parameterised with this ceiling in mind).

1.3 BJP Potential (Bouhadja et al.)¶

References

Bouhadja et al., J. Chem. Phys. 138, 224510 (2013) — CaO-Al₂O₃-SiO₂

Functional form — Born-Mayer-Huggins (BMH) + Coulomb DSF:

$$V_{ij}(r) = A_{ij}\,e^{(r_{ij} - r)/\rho_{ij}} - \frac{C_{ij}}{r^6} - \frac{D_{ij}}{r^8} + \frac{q_i q_j}{r}$$

LAMMPS pair style: born/coul/dsf 0.25 8.0

Charges (scaled formal charges):

Element	Charge e
O	−1.2
Ca	+1.2
Al	+1.8
Si	+2.4

Supported elements (4): Ca, Al, Si, O

When to use: Specifically developed and validated for calcium aluminosilicate (CAS) glasses. Best choice when studying CaO-Al₂O₃-SiO₂ with maximum fidelity to the original parameterisation.

Limitation: Very narrow element set — strictly Ca/Al/Si/O. Any other element will cause validation to fail.

1.4 BMP Potential (Bertani, Menziani & Pedone)¶

References

Bertani et al., Phys. Rev. Materials 5, 045602 (2021) — original multi-component parameterisation
Malavasi & Pedone, Acta Materialia 229, 117801 (2022) — catalase-mimetic cations in 45S5 Bioglass®
Bertani et al., J. Am. Ceram. Soc. 105, 7254 (2022) — borate and borosilicate extension
Bertani et al., J. Am. Ceram. Soc. 106, 5104 (2023) — erratum for the borosilicate parameters

Functional form — Pedone Morse + Buckingham cation–cation repulsion + three-body angle term + Coulomb DSF:

$$V_{ij}(r) = D_{ij}\left[e^{-2a_{ij}(r-r_0)} - 2e^{-a_{ij}(r-r_0)}\right] + \frac{C_{ij}}{r^{12}} + A_{ij}\,e^{-r/\rho_{ij}} + V_\text{3-body}$$

LAMMPS pair style: hybrid/overlay coul/dsf 0.25 8.0 pedone 7.0 buck 7.0 nb3b/harmonic (BMP-harm) or nb3b/screened (BMP-shrm)

The three-body term constrains cation–O–cation angles and is written to a separate parameter file (BMP.nb3b.harmonic or BMP.nb3b.shrm) at potential-generation time.

Partial charges (scaled formal charges, same convention as PMMCS):

Element	Charge e
O	−1.2
Li, Na, K, Cu, Ag, Cu2	+0.6
Be, Mg, Ca, Sr, Ba, Mn, Fe, Co, Ni, Zn, Cu2	+1.2
Sc, Al, Cr, Fe3, B, Ce3, Mn3, Nd, Gd, Er	+1.8
Ti, Zr, Si, Ge, Sn, Ce4, V4, Mn4	+2.4
P, V5	+3.0

Supported elements (38): Li, Na, K, Be, Mg, Ca, Sr, Ba, Sc, Ti, Zr, Cr, Mn, Fe, Fe3, Ce3, Ce4, Mn3, Mn4, Cu2, V5, V4, Co, Ni, Cu, Ag, Zn, Al, Si, Ge, Sn, P, Nd, Gd, Er, B, O — plus Ga/O Buckingham term

Three-body term — harm vs. shrm:

Variant	LAMMPS style	Active elements	Active triplets
`bmp-harmonic`	`nb3b/harmonic`	Si, O, P	Si–O–Si, Si–O–P, P–O–P
`bmp-screened-harmonic`	`nb3b/screened`	Si, O, P, B, V	above + B–O–B, B–O–Si, V–O–V, V–O–P, V–O–Si

Use BMP-shrm when the composition contains B or V — the screened three-body term explicitly models B–O–B and B–O–Si angles.
Use BMP-harm for B/V-free compositions (Si, P, and other oxides) when a simpler harmonic angle term is sufficient.

Boron restriction — Dell-Bray model:

When B is present, the B–O Morse depth D_ij is not a fixed parameter — it is computed at runtime from the glass composition using the Dell-Bray empirical model:

$$R = \frac{[\text{A}_2\text{O}] + [\text{AEO}]}{[\text{B}_2\text{O}_3]}, \quad K = \frac{[\text{SiO}_2]}{[\text{B}_2\text{O}_3]}$$

This model was parameterised only for systems containing B, Si, alkali oxides (Li₂O, Na₂O, K₂O), and alkaline-earth oxides (MgO, CaO, SrO, BaO). Consequently, when B is in the composition, generate_potential enforces:

Allowed with B: B, Si, O + alkali (Li, Na, K) + alkaline-earth (Mg, Ca, Sr, Ba)

Any other element co-existing with B (e.g. Al, P, Ti, transition metals) raises a ValueError because the Dell-Bray D_ij would be physically meaningless for those compositions. For B-free compositions this restriction does not apply and the full 38-element set is available.

In [2]:

Copied!





# --- BMP-harm: B-free soda-lime silicate glass ---
struct_dict_sp = am.get_structure_dict({"Na2O": 0.20, "CaO": 0.10, "SiO2": 0.70}, target_atoms=300)
pot_bmp_harmonic = generate_potential(struct_dict_sp, potential_type="bmp-harmonic", melt=True)
print(pot_bmp_harmonic[["Name", "Model", "Species"]])
print()
print("--- BMP-harm config lines ---")
for line in pot_bmp_harmonic.loc[0, "Config"]:
    print(line, end="")

print()
print("Three-body file written to:", pot_bmp_harmonic.loc[0, "Filename"])
# --- BMP-harm: B-free soda-lime silicate glass ---
struct_dict_sp = am.get_structure_dict({"Na2O": 0.20, "CaO": 0.10, "SiO2": 0.70}, target_atoms=300)
pot_bmp_harmonic = generate_potential(struct_dict_sp, potential_type="bmp-harmonic", melt=True)
print(pot_bmp_harmonic[["Name", "Model", "Species"]])
print()
print("--- BMP-harm config lines ---")
for line in pot_bmp_harmonic.loc[0, "Config"]:
    print(line, end="")

print()
print("Three-body file written to:", pot_bmp_harmonic.loc[0, "Filename"])

           Name         Model          Species
0  BMP-harmonic  BMP-harmonic  [Ca, Na, O, Si]

--- BMP-harm config lines ---
# BMP potential: Pedone (Morse) + Buckingham repulsion + three-body
units metal
dimension 3
atom_style charge

### Groups ###
group Ca type 1
group Na type 2
group O type 3
group Si type 4

### Charges ###
set type 1 charge 1.2
set type 2 charge 0.6
set type 3 charge -1.2
set type 4 charge 2.4

### BMP-harmonic Potential Parameters ###
pair_style hybrid/overlay coul/dsf 0.25 8.0 pedone 7.0 buck 7.0 nb3b/harmonic
pair_coeff * * coul/dsf
pair_coeff 1 3 pedone 0.030211 2.241334 2.923245 5.0
pair_coeff 2 3 pedone 0.023363 1.763867 3.006315 5.0
pair_coeff 3 3 pedone 0.042395 1.379316 3.618701 100.0
pair_coeff 4 3 pedone 0.340554 2.0067 2.1 1.0
pair_coeff 4 4 buck 7.093669 0.975598 0.0 7.0
pair_coeff * * nb3b/harmonic /Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/BMP.nb3b.harmonic NULL NULL O Si

pair_modify shift yes

fix langevinnve all langevin 4000 4000 0.01 48279

fix ensemblenve all nve/limit 0.5

run 10000

unfix langevinnve

unfix ensemblenve

Three-body file written to: ['/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/BMP.nb3b.harmonic']

In [3]:

Copied!





# --- BMP-shrm: alkali borosilicate glass (B present → shrm preferred) ---
struct_dict_bs = am.get_structure_dict({"Na2O": 0.15, "CaO": 0.05, "B2O3": 0.15, "SiO2": 0.65}, target_atoms=300)
pot_bmp_screened_harmonic = generate_potential(struct_dict_bs, potential_type="bmp-screened-harmonic", melt=True)
print(pot_bmp_screened_harmonic[["Name", "Model", "Species"]])
print()
print("--- BMP-shrm config lines ---")
for line in pot_bmp_screened_harmonic.loc[0, "Config"]:
    print(line, end="")

print()
print("Three-body file written to:", pot_bmp_screened_harmonic.loc[0, "Filename"])
# --- BMP-shrm: alkali borosilicate glass (B present → shrm preferred) ---
struct_dict_bs = am.get_structure_dict({"Na2O": 0.15, "CaO": 0.05, "B2O3": 0.15, "SiO2": 0.65}, target_atoms=300)
pot_bmp_screened_harmonic = generate_potential(struct_dict_bs, potential_type="bmp-screened-harmonic", melt=True)
print(pot_bmp_screened_harmonic[["Name", "Model", "Species"]])
print()
print("--- BMP-shrm config lines ---")
for line in pot_bmp_screened_harmonic.loc[0, "Config"]:
    print(line, end="")

print()
print("Three-body file written to:", pot_bmp_screened_harmonic.loc[0, "Filename"])

                    Name                  Model             Species
0  BMP-screened-harmonic  BMP-screened-harmonic  [B, Ca, Na, O, Si]

--- BMP-shrm config lines ---
# BMP potential: Pedone (Morse) + Buckingham repulsion + three-body
units metal
dimension 3
atom_style charge

### Groups ###
group B type 1
group Ca type 2
group Na type 3
group O type 4
group Si type 5

### Charges ###
set type 1 charge 1.8
set type 2 charge 1.2
set type 3 charge 0.6
set type 4 charge -1.2
set type 5 charge 2.4

### BMP-screened-harmonic Potential Parameters ###
pair_style hybrid/overlay coul/dsf 0.25 8.0 pedone 7.0 buck 7.0 nb3b/screened
pair_coeff * * coul/dsf
pair_coeff 1 4 pedone 2.8022423265597305 2.64333 1.43667 0.0
pair_coeff 2 4 pedone 0.030211 2.241334 2.923245 5.0
pair_coeff 3 4 pedone 0.023363 1.763867 3.006315 5.0
pair_coeff 4 4 pedone 0.042395 1.379316 3.618701 100.0
pair_coeff 5 4 pedone 0.340554 2.0067 2.1 1.0
pair_coeff 5 5 buck 7.093669 0.975598 0.0 7.0
pair_coeff 1 1 buck 8.9594 0.8012 0.0 7.0
pair_coeff 1 5 buck 8.9594 0.927 0.0 7.0
pair_coeff * * nb3b/screened /Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/BMP.nb3b.shrm B NULL NULL O Si

pair_modify shift yes

fix langevinnve all langevin 4000 4000 0.01 48279

fix ensemblenve all nve/limit 0.5

run 10000

unfix langevinnve

unfix ensemblenve

Three-body file written to: ['/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/BMP.nb3b.shrm']

In [4]:

Copied!





# --- Boron restriction: Al co-existing with B raises ValueError ---
struct_dict_bad = am.get_structure_dict({"Na2O": 0.10, "Al2O3": 0.10, "B2O3": 0.10, "SiO2": 0.70}, target_atoms=200)
try:
    generate_potential(struct_dict_bad, potential_type="bmp-screened-harmonic")
except ValueError as e:
    print("Expected error:", e)

# --- Without B, Al is perfectly fine with BMP ---
struct_dict_al = am.get_structure_dict({"Na2O": 0.10, "Al2O3": 0.10, "SiO2": 0.80}, target_atoms=200)
pot_al = generate_potential(struct_dict_al, potential_type="bmp-harmonic", melt=False)
print("\nAl without B — OK:", pot_al.loc[0, "Species"])
# --- Boron restriction: Al co-existing with B raises ValueError ---
struct_dict_bad = am.get_structure_dict({"Na2O": 0.10, "Al2O3": 0.10, "B2O3": 0.10, "SiO2": 0.70}, target_atoms=200)
try:
    generate_potential(struct_dict_bad, potential_type="bmp-screened-harmonic")
except ValueError as e:
    print("Expected error:", e)

# --- Without B, Al is perfectly fine with BMP ---
struct_dict_al = am.get_structure_dict({"Na2O": 0.10, "Al2O3": 0.10, "SiO2": 0.80}, target_atoms=200)
pot_al = generate_potential(struct_dict_al, potential_type="bmp-harmonic", melt=False)
print("\nAl without B — OK:", pot_al.loc[0, "Species"])

Expected error: BMP potential with boron is restricted to alkali/alkaline-earth borate and borosilicate glasses. Unsupported elements: ['Al']

Al without B — OK: ['Al', 'Na', 'O', 'Si']

2. Automatic Potential Selection¶

select_potential(elements) returns the best available potential name for a given element set (preference order: PMMCS → SHIK → BJP).
compatible_potentials(elements) returns all compatible names.

In [5]:

Copied!





# Inspect supported elements per potential
for name in ("pmmcs", "bmp-harmonic", "bmp-screened-harmonic", "shik", "bjp"):
    elems = sorted(get_supported_elements(name))
    print(f"{name.upper():10s} ({len(elems)} elements): {', '.join(elems)}")

print()

# Automatic selection for different glass families
examples = {
    "Na₂O-SiO₂ (soda-silicate)": {"Na", "Si", "O"},
    "Na₂O-B₂O₃-SiO₂ (borosilicate)": {"Na", "B", "Si", "O"},
    "Na₂O-CaO-B₂O₃-SiO₂ (alkali borosil.)": {"Na", "Ca", "B", "Si", "O"},
    "Na₂O-SiO₂-P₂O₅ (silicophosphate)": {"Na", "Si", "P", "O"},
    "CaO-Al₂O₃-SiO₂ (CAS)": {"Ca", "Al", "Si", "O"},
    "ZrO₂-SiO₂ (zircosilicate)": {"Zr", "Si", "O"},
}
print(f"{'Composition':<45} {'Best':>9}  All compatible")
print("-" * 75)
for label, elements in examples.items():
    best = select_potential(elements)
    all_compat = compatible_potentials(elements)
    print(f"{label:<45} {best!s:>9}  {all_compat}")
# Inspect supported elements per potential
for name in ("pmmcs", "bmp-harmonic", "bmp-screened-harmonic", "shik", "bjp"):
    elems = sorted(get_supported_elements(name))
    print(f"{name.upper():10s} ({len(elems)} elements): {', '.join(elems)}")

print()

# Automatic selection for different glass families
examples = {
    "Na₂O-SiO₂ (soda-silicate)": {"Na", "Si", "O"},
    "Na₂O-B₂O₃-SiO₂ (borosilicate)": {"Na", "B", "Si", "O"},
    "Na₂O-CaO-B₂O₃-SiO₂ (alkali borosil.)": {"Na", "Ca", "B", "Si", "O"},
    "Na₂O-SiO₂-P₂O₅ (silicophosphate)": {"Na", "Si", "P", "O"},
    "CaO-Al₂O₃-SiO₂ (CAS)": {"Ca", "Al", "Si", "O"},
    "ZrO₂-SiO₂ (zircosilicate)": {"Zr", "Si", "O"},
}
print(f"{'Composition':<45} {'Best':>9}  All compatible")
print("-" * 75)
for label, elements in examples.items():
    best = select_potential(elements)
    all_compat = compatible_potentials(elements)
    print(f"{label:<45} {best!s:>9}  {all_compat}")

PMMCS      (29 elements): Ag, Al, Ba, Be, Ca, Co, Cr, Cu, Er, Fe, Fe3, Gd, Ge, K, Li, Mg, Mn, Na, Nd, Ni, O, P, Sc, Si, Sn, Sr, Ti, Zn, Zr
BMP-HARMONIC (36 elements): Ag, Al, B, Ba, Be, Ca, Ce3, Ce4, Co, Cr, Cu, Cu2, Er, Fe, Fe3, Gd, Ge, K, Li, Mg, Mn, Mn3, Na, Nd, Ni, O, P, Sc, Si, Sn, Sr, Ti, V4, V5, Zn, Zr
BMP-SCREENED-HARMONIC (36 elements): Ag, Al, B, Ba, Be, Ca, Ce3, Ce4, Co, Cr, Cu, Cu2, Er, Fe, Fe3, Gd, Ge, K, Li, Mg, Mn, Mn3, Na, Nd, Ni, O, P, Sc, Si, Sn, Sr, Ti, V4, V5, Zn, Zr
SHIK       (9 elements): Al, B, Ca, K, Li, Mg, Na, O, Si
BJP        (4 elements): Al, Ca, O, Si

Composition                                        Best  All compatible
---------------------------------------------------------------------------
Na₂O-SiO₂ (soda-silicate)                         pmmcs  ['pmmcs', 'bmp-harmonic', 'bmp-screened-harmonic', 'shik', 'du_teter', 'yang2026']
Na₂O-B₂O₃-SiO₂ (borosilicate)                 bmp-screened-harmonic  ['bmp-harmonic', 'bmp-screened-harmonic', 'shik', 'du_teter', 'yang2026']
Na₂O-CaO-B₂O₃-SiO₂ (alkali borosil.)          bmp-screened-harmonic  ['bmp-harmonic', 'bmp-screened-harmonic', 'shik', 'du_teter', 'yang2026']
Na₂O-SiO₂-P₂O₅ (silicophosphate)                  pmmcs  ['pmmcs', 'bmp-harmonic', 'bmp-screened-harmonic', 'du_teter']
CaO-Al₂O₃-SiO₂ (CAS)                              pmmcs  ['pmmcs', 'bmp-harmonic', 'bmp-screened-harmonic', 'shik', 'bjp', 'du_teter']
ZrO₂-SiO₂ (zircosilicate)                         pmmcs  ['pmmcs', 'bmp-harmonic', 'bmp-screened-harmonic', 'du_teter']

3. Structure Generation¶

get_structure_dict builds the internal representation (a dict with "atoms" and "box" keys) from an oxide composition.
get_ase_structure converts that dict to an ASE Atoms object needed by the simulation functions.

The composition is specified as a dictionary of oxide formulas → molar fractions.

In [6]:

Copied!





# Example: 25 CaO - 25 Al2O3 - 50 SiO2 (mol%)
composition = {"CaO": 0.25, "Al2O3": 0.25, "SiO2": 0.50}

struct_dict = am.get_structure_dict(composition, target_atoms=500)
print("Box length:", struct_dict["box"], "Å")
print("Atom count:", len(struct_dict["atoms"]))

counts = Counter(a["element"] for a in struct_dict["atoms"])
print("Element counts:", dict(counts))
# Example: 25 CaO - 25 Al2O3 - 50 SiO2 (mol%)
composition = {"CaO": 0.25, "Al2O3": 0.25, "SiO2": 0.50}

struct_dict = am.get_structure_dict(composition, target_atoms=500)
print("Box length:", struct_dict["box"], "Å")
print("Atom count:", len(struct_dict["atoms"]))

counts = Counter(a["element"] for a in struct_dict["atoms"])
print("Element counts:", dict(counts))

Box length: 18.759772823007427 Å
Atom count: 499
Element counts: {'Ca': 39, 'O': 307, 'Al': 76, 'Si': 77}

In [7]:

Copied!





# Convert to ASE Atoms object
atoms = am.get_ase_structure(struct_dict)
print(atoms)
print("Cell:", atoms.cell.array.diagonal().tolist(), "Å")
print("PBC:", atoms.pbc)
# Convert to ASE Atoms object
atoms = am.get_ase_structure(struct_dict)
print(atoms)
print("Cell:", atoms.cell.array.diagonal().tolist(), "Å")
print("PBC:", atoms.pbc)

Atoms(symbols='Al76Ca39O307Si77', pbc=True, cell=[18.759772823007427, 18.759772823007427, 18.759772823007427], id=..., initial_charges=..., masses=..., mmcharges=..., type=...)
Cell: [18.759772823007427, 18.759772823007427, 18.759772823007427] Å
PBC: [ True  True  True]

4. Potential Generation¶

generate_potential(atoms_dict, potential_type) returns a pd.DataFrame with columns:

Column	Content
`Name`	Potential name (`"PMMCS"`, `"SHIK"`, `"BJP"`)
`Model`	Full model descriptor string
`Species`	List of element symbols in the structure
`Config`	List of LAMMPS input lines
`Filename`	Table file paths (SHIK only; empty list for PMMCS/BJP)

The Config lines are injected directly before the LAMMPS run commands by the simulation engine.

In [8]:

Copied!





# Generate BJP potential (CAS system)
potential_bjp = generate_potential(struct_dict, potential_type="bjp")
print(potential_bjp[["Name", "Model", "Species"]])
print()
print("--- LAMMPS config lines ---")
for line in potential_bjp.loc[0, "Config"]:
    print(line, end="")
# Generate BJP potential (CAS system)
potential_bjp = generate_potential(struct_dict, potential_type="bjp")
print(potential_bjp[["Name", "Model", "Species"]])
print()
print("--- LAMMPS config lines ---")
for line in potential_bjp.loc[0, "Config"]:
    print(line, end="")

  Name                           Model          Species
0  BJP  Born-Mayer-Huggins_coulomb_DSF  [Al, Ca, O, Si]

--- LAMMPS config lines ---
# Bouhadja et al., J. Chem. Phys. 138, 224510 (2013) 
units metal
dimension 3
atom_style charge

# create groups ###
group Al type 1
group Ca type 2
group O type 3
group Si type 4

### set charges ###
set type 1 charge 1.8
set type 2 charge 1.2
set type 3 charge -1.2
set type 4 charge 2.4

### Bouhadja Born-Mayer-Huggins + Coulomb Potential Parameters ###
pair_style born/coul/dsf 0.25 8.0
pair_coeff 1 1 0.002900 0.068000 1.570400 14.049800 0.000000
pair_coeff 1 2 0.003200 0.074000 1.957200 17.171000 0.000000
pair_coeff 1 3 0.007500 0.164000 2.606700 34.574700 0.000000
pair_coeff 1 4 0.002500 0.057000 1.505600 18.811600 0.000000
pair_coeff 2 2 0.003500 0.080000 2.344000 20.985600 0.000000
pair_coeff 2 3 0.007700 0.178000 2.993500 42.255600 0.000000
pair_coeff 2 4 0.002700 0.063000 1.892400 22.990700 0.000000
pair_coeff 3 3 0.012000 0.263000 3.643000 85.084000 0.000000
pair_coeff 3 4 0.007000 0.156000 2.541900 46.293000 0.000000
pair_coeff 4 4 0.001200 0.046000 1.440800 25.187300 0.000000

pair_modify shift yes

In [9]:

Copied!





# Generate PMMCS potential for the same system
potential_pmmcs = generate_potential(struct_dict, potential_type="pmmcs")
print("--- PMMCS config lines ---")
for line in potential_pmmcs.loc[0, "Config"]:
    print(line, end="")
# Generate PMMCS potential for the same system
potential_pmmcs = generate_potential(struct_dict, potential_type="pmmcs")
print("--- PMMCS config lines ---")
for line in potential_pmmcs.loc[0, "Config"]:
    print(line, end="")

--- PMMCS config lines ---
# A. Pedone et.al., JPCB (2006), https://doi.org/10.1021/jp0611018
units metal
dimension 3
atom_style charge

# create groups ###
group Al type 1
group Ca type 2
group O type 3
group Si type 4

### set charges ###
set type 1 charge 1.8
set type 2 charge 1.2
set type 3 charge -1.2
set type 4 charge 2.4

### Pmmcs Potential Parameters ###
pair_style hybrid/overlay coul/dsf 0.25 8.0 pedone 5.5
pair_coeff * * coul/dsf
pair_coeff 1 3 pedone 0.361581 1.900442 2.164818 0.9
pair_coeff 2 3 pedone 0.030211 2.241334 2.923245 5.0
pair_coeff 3 3 pedone 0.042395 1.379316 3.618701 22.0
pair_coeff 4 3 pedone 0.340554 2.0067 2.1 1.0

pair_modify shift yes

4.1 Custom electrostatics settings¶

By default all three potentials use DSF (Damped Shifted Force) Coulomb with their built-in cutoffs and damping parameter α. Pass a config object to override any of these:

Class	Parameters	Supported potentials
`DsfConfig`	`alpha`, `long_range_cutoff`	PMMCS, BJP, SHIK
`WolfConfig`	`alpha`, `long_range_cutoff`	PMMCS, BJP
`PppmConfig`	`long_range_cutoff`, `kspace_accuracy`	PMMCS, BJP
`EwaldConfig`	`long_range_cutoff`, `kspace_accuracy`	PMMCS, BJP

PppmConfig and EwaldConfig append a kspace_style line and do not use alpha.

In [10]:

Copied!





from amorphouspy.lammps.potentials.potential import DsfConfig, WolfConfig, PppmConfig

# --- BJP: stronger damping with DSF ---
cfg_bjp_custom = DsfConfig(alpha=0.20)
pot_bjp_custom = generate_potential(struct_dict, potential_type="bjp", melt=False, electrostatics=cfg_bjp_custom)
pair_style_line = next(line for line in pot_bjp_custom.loc[0, "Config"] if "pair_style" in line)
print("BJP custom DSF:", pair_style_line.strip())

# --- BJP: Wolf method instead of DSF ---
cfg_bjp_wolf = WolfConfig(alpha=0.20)
pot_bjp_wolf = generate_potential(struct_dict, potential_type="bjp", melt=False, electrostatics=cfg_bjp_wolf)
pair_style_line = next(line for line in pot_bjp_wolf.loc[0, "Config"] if "pair_style" in line)
print("BJP Wolf:      ", pair_style_line.strip())

# --- PMMCS: PPPM instead of DSF (appends a kspace_style line) ---
cfg_pmmcs_pppm = PppmConfig(long_range_cutoff=12.0, kspace_accuracy=1e-6)
pot_pmmcs_pppm = generate_potential(struct_dict, potential_type="pmmcs", melt=False, electrostatics=cfg_pmmcs_pppm)
config_lines = pot_pmmcs_pppm.loc[0, "Config"]
for line in config_lines:
    if "pair_style" in line or "kspace_style" in line:
        print("PMMCS PPPM:    ", line.strip())
from amorphouspy.lammps.potentials.potential import DsfConfig, WolfConfig, PppmConfig

# --- BJP: stronger damping with DSF ---
cfg_bjp_custom = DsfConfig(alpha=0.20)
pot_bjp_custom = generate_potential(struct_dict, potential_type="bjp", melt=False, electrostatics=cfg_bjp_custom)
pair_style_line = next(line for line in pot_bjp_custom.loc[0, "Config"] if "pair_style" in line)
print("BJP custom DSF:", pair_style_line.strip())

# --- BJP: Wolf method instead of DSF ---
cfg_bjp_wolf = WolfConfig(alpha=0.20)
pot_bjp_wolf = generate_potential(struct_dict, potential_type="bjp", melt=False, electrostatics=cfg_bjp_wolf)
pair_style_line = next(line for line in pot_bjp_wolf.loc[0, "Config"] if "pair_style" in line)
print("BJP Wolf:      ", pair_style_line.strip())

# --- PMMCS: PPPM instead of DSF (appends a kspace_style line) ---
cfg_pmmcs_pppm = PppmConfig(long_range_cutoff=12.0, kspace_accuracy=1e-6)
pot_pmmcs_pppm = generate_potential(struct_dict, potential_type="pmmcs", melt=False, electrostatics=cfg_pmmcs_pppm)
config_lines = pot_pmmcs_pppm.loc[0, "Config"]
for line in config_lines:
    if "pair_style" in line or "kspace_style" in line:
        print("PMMCS PPPM:    ", line.strip())

BJP custom DSF: pair_style born/coul/dsf 0.2 8.0
BJP Wolf:       pair_style born/coul/wolf 0.2 8.0
PMMCS PPPM:     pair_style hybrid/overlay coul/long 12.0 pedone 5.5
PMMCS PPPM:     kspace_style pppm 1e-06

In [11]:

Copied!





# SHIK potential — writes .tbl files to the "shik_tables" directory by default
# The table files encode V(r) and F(r) on a fine grid for each pair

# Use a borosilicate composition to exercise the B parameters
struct_dict_bs = am.get_structure_dict({"Na2O": 0.15, "B2O3": 0.15, "SiO2": 0.7}, target_atoms=500)
atoms = am.get_ase_structure(struct_dict_bs)
potential_shik = generate_potential(struct_dict_bs, potential_type="shik", melt=True)
print(potential_shik[["Name", "Model", "Species"]])
print()
print("--- SHIK config lines ---")
for line in potential_shik.loc[0, "Config"]:
    print(line, end="")
# SHIK potential — writes .tbl files to the "shik_tables" directory by default
# The table files encode V(r) and F(r) on a fine grid for each pair

# Use a borosilicate composition to exercise the B parameters
struct_dict_bs = am.get_structure_dict({"Na2O": 0.15, "B2O3": 0.15, "SiO2": 0.7}, target_atoms=500)
atoms = am.get_ase_structure(struct_dict_bs)
potential_shik = generate_potential(struct_dict_bs, potential_type="shik", melt=True)
print(potential_shik[["Name", "Model", "Species"]])
print()
print("--- SHIK config lines ---")
for line in potential_shik.loc[0, "Config"]:
    print(line, end="")

   Name Model         Species
0  SHIK  SHIK  [B, Na, O, Si]

--- SHIK config lines ---
# S. Sundararaman et al.,
# for silica glass, J. Chem. Phys. 2018, 148, 19, https://doi.org/10.1063/1.5023707
# for alkali and alkaline-earth aluminosilicate glasses, J. Chem. Phys. 2019, 150, 15, https://doi.org/10.1063/1.5079663
# for borate glasses with mixed network formers, J. Chem. Phys. 2020, 152, 10, https://doi.org/10.1063/1.5142605
# for alkaline earth silicate and borate glasses, Shih et al. J. Non-Cryst. Sol. 2021, 565, 120853, https://doi.org/10.1016/j.jnoncrysol.2021.120853

dimension 3
atom_style charge


### Group Definitions ###
group B type 1
group Na type 2
group O type 3
group Si type 4

### Charges ###
set type 1 charge 1.6126
set type 2 charge 0.6018
set type 3 charge -0.9541625
set type 4 charge 1.7755

### SHIK Potential ###
pair_style hybrid/overlay coul/dsf 0.2 10.0 table spline 10000
pair_coeff * * coul/dsf
pair_coeff 1 1 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_B_B.tbl" SHIK_Buck_r24 8.0
pair_coeff 1 2 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_B_Na.tbl" SHIK_Buck_r24 8.0
pair_coeff 1 3 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_B_O.tbl" SHIK_Buck_r24 8.0
pair_coeff 1 4 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_B_Si.tbl" SHIK_Buck_r24 8.0
pair_coeff 2 2 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_Na_Na.tbl" SHIK_Buck_r24 8.0
pair_coeff 2 3 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_Na_O.tbl" SHIK_Buck_r24 8.0
pair_coeff 2 4 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_Na_Si.tbl" SHIK_Buck_r24 8.0
pair_coeff 3 3 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_O_O.tbl" SHIK_Buck_r24 8.0
pair_coeff 3 4 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_O_Si.tbl" SHIK_Buck_r24 8.0
pair_coeff 4 4 table "/Users/achrafatila/Documents/Workflows/amorphouspy/notebooks/shik_tables/table_Si_Si.tbl" SHIK_Buck_r24 8.0

pair_modify shift yes


thermo_style custom step temp pe etotal pxx pxy pxz pyy pyz pzz vol

thermo_modify flush yes

thermo 100

fix langevinnve all langevin 4000 4000 0.01 48279

fix ensemblenve all nve/limit 0.5

run 10000

unfix langevinnve

unfix ensemblenve

5. Melt-Quench Simulation¶

melt_quench_simulation runs a multi-stage protocol that:

Heats the random structure to temperature_high
Equilibrates the melt
Cools to temperature_low at the specified cooling_rate
Releases pressure (NPT at 0 GPa)
Final NVT equilibration at room temperature

Each potential has a dedicated protocol with different ensemble choices and stage durations. These differ because the potentials have different stability characteristics.

Protocol comparison¶

Stage	PMMCS	BMP (harm/shrm)	SHIK	BJP
1 Heating	NVT ramp T_low → T_high (with Langevin+NVE/limit init)	NVT ramp T_low → T_high (with Langevin+NVE/limit init)	NVT from T_high (with Langevin+NVE/limit init)	NPT ramp T_low → T_high
2 Melt equil.	NVT, 1 M steps (default)	NVT, 1 M steps (default)	NVT, 100 ps	NPT, 100k steps
3 NPT equil.	—	—	NPT 0.1 GPa, 700 ps	—
4 Quench	NVT ramp T_high → T_low	NVT ramp T_high → T_low	NPT quench, pressure 0.1→0 GPa	NPT ramp T_high → T_low
5 Pressure release	NPT 0 GPa, 1 M steps	NPT 0 GPa, 1 M steps	—	NPT 0 GPa, 100k steps
6 Final equil.	NVT, 100k steps	NVT, 100k steps	NPT 0 GPa, 100 ps	NVT, 100k steps

BMP init block: Stage 1 runs with the full potential including the Langevin + NVE/limit pre-equilibration block embedded in the config (the melt=True block). Stages 2–5 use the same potential with that block stripped, so it fires exactly once.

Key parameters¶

Parameter	Default	Description
`temperature_high`	PMMCS/BJP: 5000 K, BMP/SHIK: 4000 K	Melt temperature
`temperature_low`	300 K	Final glass temperature
`timestep`	1.0 fs	MD integration timestep
`heating_rate`	1×10¹² K/s	Rate for heating ramp
`cooling_rate`	1×10¹² K/s	Rate for quenching ramp
`equilibration_steps`	None (protocol defaults)	Override all equilibration stages
`n_print`	1000	Output frequency (steps)
`langevin`	False	Use Langevin thermostat
`seed`	12345	Velocity initialisation seed

Step count calculation¶

The number of heating/cooling steps is derived from the rate:

steps = ΔT / (timestep [fs] × rate [K/s]) × 1e15 [fs/s]

For example, with ΔT = 4700 K, timestep = 1 fs, and rate = 1×10¹² K/s:

steps = 4700 / (1×10⁻¹⁵ × 1×10¹²) = 4 700 000 steps  (~4.7 ns)

A faster rate (e.g. 1×10¹³ K/s) gives 10× fewer steps and a shorter simulation.

In [12]:

Copied!





result = am.melt_quench_simulation(
    structure=atoms,
    potential=potential_shik,
    temperature_high=5000.0,
    temperature_low=300.0,
    timestep=1.0,
    heating_rate=1e15,  # fast for demonstration
    cooling_rate=1e15,
    n_print=100,
    equilibration_steps=10000,  # override protocol defaults — production: remove this line
    seed=42,
)

glass = result["structure"]
history = result["result"]  # list[dict | None], one entry per protocol stage
print("Final structure:", glass)
print(f"Protocol stages: {len(history)}")
for i, stage in enumerate(history):
    keys = list(stage.keys()) if stage else None
    print(f"  Stage {i + 1}: {keys}")
result = am.melt_quench_simulation(
    structure=atoms,
    potential=potential_shik,
    temperature_high=5000.0,
    temperature_low=300.0,
    timestep=1.0,
    heating_rate=1e15,  # fast for demonstration
    cooling_rate=1e15,
    n_print=100,
    equilibration_steps=10000,  # override protocol defaults — production: remove this line
    seed=42,
)

glass = result["structure"]
history = result["result"]  # list[dict | None], one entry per protocol stage
print("Final structure:", glass)
print(f"Protocol stages: {len(history)}")
for i, stage in enumerate(history):
    keys = list(stage.keys()) if stage else None
    print(f"  Stage {i + 1}: {keys}")

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity LogStep. Returning un-normalized quantity

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity

Final structure: Atoms(symbols='B46Na46O304Si106', pbc=True, cell=[[21.08330128860441, 3.872939615782298e-15, 3.872939615782299e-15], [-3.872939615782299e-15, 21.08330128860441, 3.872939615782299e-15], [-3.872939615782299e-15, -3.8729396157823e-15, 21.08330128860441]], id=..., indices=..., initial_charges=..., masses=..., mmcharges=..., momenta=..., type=...)
Protocol stages: 6
  Stage 1: ['steps', 'natoms', 'cells', 'indices', 'forces', 'velocities', 'unwrapped_positions', 'positions', 'temperature', 'energy_pot', 'energy_tot', 'volume', 'LogStep', 'pressures']
  Stage 2: ['steps', 'natoms', 'cells', 'indices', 'forces', 'velocities', 'unwrapped_positions', 'positions', 'temperature', 'energy_pot', 'energy_tot', 'volume', 'pressures']
  Stage 3: ['steps', 'natoms', 'cells', 'indices', 'forces', 'velocities', 'unwrapped_positions', 'positions', 'temperature', 'energy_pot', 'energy_tot', 'volume', 'pressures']
  Stage 4: ['steps', 'natoms', 'cells', 'indices', 'forces', 'velocities', 'unwrapped_positions', 'positions', 'temperature', 'energy_pot', 'energy_tot', 'volume', 'pressures']
  Stage 5: ['steps', 'natoms', 'cells', 'indices', 'forces', 'velocities', 'unwrapped_positions', 'positions', 'temperature', 'energy_pot', 'energy_tot', 'volume', 'pressures']
  Stage 6: ['steps', 'natoms', 'cells', 'indices', 'forces', 'velocities', 'unwrapped_positions', 'positions', 'temperature', 'energy_pot', 'energy_tot', 'volume', 'pressures']

/Users/achrafatila/miniforge3/envs/amorphouspy/lib/python3.13/site-packages/lammpsparser/units.py:242: UserWarning:

Warning: Couldn't determine the LAMMPS to pyiron unit conversion type of quantity Density. Returning un-normalized quantity

Parallel runs with executorlib¶

For multiple independent compositions or cooling rates, use executorlib.SingleNodeExecutor to dispatch jobs as futures.

from executorlib import SingleNodeExecutor

compositions = {
    "CAS_25_25_50": {"CaO": 0.25, "Al2O3": 0.25, "SiO2": 0.50},
    "CAS_30_10_60": {"CaO": 0.30, "Al2O3": 0.10, "SiO2": 0.60},
}


def run_one(label, comp):
    sd = am.get_structure_dict(comp, target_atoms=500)
    pot = generate_potential(sd, potential_type="bjp")
    ase = am.get_ase_structure(sd)
    res = am.melt_quench_simulation(
        structure=ase,
        potential=pot,
        heating_rate=1e15,
        cooling_rate=1e15,
        equilibration_steps=5000,
        seed=42,
    )
    return label, res["structure"]


with SingleNodeExecutor(max_workers=2) as exe:
    futures = {label: exe.submit(run_one, label, comp) for label, comp in compositions.items()}

results = {label: fut.result() for label, fut in futures.items()}
for label, (_name, struct) in results.items():
    print(f"{label}: {len(struct)} atoms")

6. Post-processing and Plotting¶

melt_quench_simulation now returns the full per-stage simulation history: result["result"] is a list[dict | None] with one thermo dict per protocol stage (in order). Each dict has the same keys as a single md_simulation result.

merge_thermo concatenates the list into one continuous trace; plot_schedule visualises temperature, potential energy, and volume per stage with colour coding.

In [13]:

Copied!





def merge_thermo(history: list[dict]) -> dict:
    """Concatenate per-stage thermo dicts into one continuous trace (steps monotonically increasing)."""
    if not history:
        return {}
    merged = {}
    step_offset = 0
    for stage in history:
        if stage is None:
            continue
        for key, values in stage.items():
            arr = np.asarray(values)
            if key in ("steps", "step"):
                arr = arr + step_offset
            merged.setdefault(key, [])
            merged[key].extend(arr.tolist())
        steps_key = "steps" if "steps" in stage else "step"
        if steps_key in stage:
            step_offset += int(np.asarray(stage[steps_key])[-1])
    return {k: np.array(v) for k, v in merged.items()}


def plot_schedule(history: list[dict], labels: list[str] | None = None, timestep_fs: float = 1.0):
    """Plot temperature, potential energy, and volume for a multi-stage schedule."""
    _fig, axes = plt.subplots(1, 3, figsize=(3.5 * 3, 3.5 / 1.3333), dpi=300)
    step_offset = 0
    colors = plt.cm.tab10.colors

    for i, (stage, color) in enumerate(zip(history, colors)):
        if stage is None:
            continue
        steps = np.asarray(stage.get("steps", stage.get("step", [])))
        time_ps = (steps + step_offset) * timestep_fs * 1e-3
        label = labels[i] if labels else f"Stage {i + 1}"

        for ax, key, ylabel in zip(
            axes,
            ["temperature", "energy_pot", "volume"],
            ["Temperature (K)", "Potential energy (eV)", "Volume (Å³)"],
        ):
            data = stage.get(key)
            if data is not None:
                ax.plot(time_ps, data, lw=0.8, color=color, label=label)
                ax.set_xlabel("Time (ps)")
                ax.set_ylabel(ylabel)

        if len(steps):
            step_offset += int(steps[-1])

    for ax in axes:
        ax.legend(fontsize=7)
    plt.tight_layout()
    plt.show()
def merge_thermo(history: list[dict]) -> dict:
    """Concatenate per-stage thermo dicts into one continuous trace (steps monotonically increasing)."""
    if not history:
        return {}
    merged = {}
    step_offset = 0
    for stage in history:
        if stage is None:
            continue
        for key, values in stage.items():
            arr = np.asarray(values)
            if key in ("steps", "step"):
                arr = arr + step_offset
            merged.setdefault(key, [])
            merged[key].extend(arr.tolist())
        steps_key = "steps" if "steps" in stage else "step"
        if steps_key in stage:
            step_offset += int(np.asarray(stage[steps_key])[-1])
    return {k: np.array(v) for k, v in merged.items()}


def plot_schedule(history: list[dict], labels: list[str] | None = None, timestep_fs: float = 1.0):
    """Plot temperature, potential energy, and volume for a multi-stage schedule."""
    _fig, axes = plt.subplots(1, 3, figsize=(3.5 * 3, 3.5 / 1.3333), dpi=300)
    step_offset = 0
    colors = plt.cm.tab10.colors

    for i, (stage, color) in enumerate(zip(history, colors)):
        if stage is None:
            continue
        steps = np.asarray(stage.get("steps", stage.get("step", [])))
        time_ps = (steps + step_offset) * timestep_fs * 1e-3
        label = labels[i] if labels else f"Stage {i + 1}"

        for ax, key, ylabel in zip(
            axes,
            ["temperature", "energy_pot", "volume"],
            ["Temperature (K)", "Potential energy (eV)", "Volume (Å³)"],
        ):
            data = stage.get(key)
            if data is not None:
                ax.plot(time_ps, data, lw=0.8, color=color, label=label)
                ax.set_xlabel("Time (ps)")
                ax.set_ylabel(ylabel)

        if len(steps):
            step_offset += int(steps[-1])

    for ax in axes:
        ax.legend(fontsize=7)
    plt.tight_layout()
    plt.show()

In [14]:

Copied!





# Per-stage view — each stage has its own colour
shik_stage_labels = [
    "Heat (NVT)",
    "Equil. melt (NVT)",
    "NPT equil.",
    "Quench (NPT)",
    "Anneal (NPT)",
    "Final equil. (NVT)",
]
plot_schedule(history, labels=shik_stage_labels)

# Continuous trace using merge_thermo
merged = merge_thermo(history)
time_ps = merged["steps"] * 1e-3  # 1 fs timestep → ps

_fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(time_ps, merged["temperature"], lw=0.8, color="tomato")
axes[0].set_xlabel("Time (ps)")
axes[0].set_ylabel("Temperature (K)")
axes[0].set_title("Full melt-quench — temperature trace")

axes[1].plot(time_ps, merged["volume"], lw=0.8, color="steelblue")
axes[1].set_xlabel("Time (ps)")
axes[1].set_ylabel("Volume (Å³)")
axes[1].set_title("Full melt-quench — volume trace")

plt.tight_layout()
plt.show()
# Per-stage view — each stage has its own colour
shik_stage_labels = [
    "Heat (NVT)",
    "Equil. melt (NVT)",
    "NPT equil.",
    "Quench (NPT)",
    "Anneal (NPT)",
    "Final equil. (NVT)",
]
plot_schedule(history, labels=shik_stage_labels)

# Continuous trace using merge_thermo
merged = merge_thermo(history)
time_ps = merged["steps"] * 1e-3  # 1 fs timestep → ps

_fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(time_ps, merged["temperature"], lw=0.8, color="tomato")
axes[0].set_xlabel("Time (ps)")
axes[0].set_ylabel("Temperature (K)")
axes[0].set_title("Full melt-quench — temperature trace")

axes[1].plot(time_ps, merged["volume"], lw=0.8, color="steelblue")
axes[1].set_xlabel("Time (ps)")
axes[1].set_ylabel("Volume (Å³)")
axes[1].set_title("Full melt-quench — volume trace")

plt.tight_layout()
plt.show()

No description has been provided for this image

Density from the equilibrated structure¶

The volume of the simulation box after NPT relaxation gives the glass density directly.

In [15]:

Copied!

from ase.data import atomic_masses

def compute_density(structure):
    """Density in g/cm³ from an ASE Atoms object."""
    mass_amu = sum(atomic_masses[n] for n in structure.get_atomic_numbers())
    volume_A3 = structure.get_volume()
    return (mass_amu * 1.66054e-24) / (volume_A3 * 1e-24)

rho = compute_density(glass)
print(f"Glass density: {rho:.3f} g/cm³")
from ase.data import atomic_masses

def compute_density(structure):
    """Density in g/cm³ from an ASE Atoms object."""
    mass_amu = sum(atomic_masses[n] for n in structure.get_atomic_numbers())
    volume_A3 = structure.get_volume()
    return (mass_amu * 1.66054e-24) / (volume_A3 * 1e-24)

rho = compute_density(glass)
print(f"Glass density: {rho:.3f} g/cm³")

Glass density: 1.665 g/cm³

Structural analysis¶

After the melt-quench the final Atoms structure can be passed directly to the analysis tools. See StructureTutorial.ipynb for a full walkthrough; here is a quick RDF check.

In [16]:

Copied!





# RDF for the quenched glass — pairs present in the SHIK borosilicate composition
r, rdfs, cumcn = am.compute_rdf(
    structure=glass,
    r_max=8.0,
    n_bins=500,
    type_pairs=[(14, 8), (5, 8), (11, 8)],  # Si-O, B-O, Na-O
)

fig, ax = plt.subplots(figsize=(3.5, 3.5 / 1.3333), dpi=300)
labels = {(8, 14): "Si-O", (5, 8): "B-O", (8, 11): "Na-O"}
for (zi, zj), g in rdfs.items():
    ax.plot(r, g, label=labels.get((zi, zj), f"{zi}-{zj}"))
ax.axhline(1, ls="--", color="gray", lw=0.7)
ax.set_xlabel("r (Å)")
ax.set_ylabel("g(r)")
ax.set_xlim(0, 5)
ax.set_ylim(0, None)
ax.set_title("Partial RDFs — quenched glass")
ax.legend()
plt.tight_layout()
plt.show()
# RDF for the quenched glass — pairs present in the SHIK borosilicate composition
r, rdfs, cumcn = am.compute_rdf(
    structure=glass,
    r_max=8.0,
    n_bins=500,
    type_pairs=[(14, 8), (5, 8), (11, 8)],  # Si-O, B-O, Na-O
)

fig, ax = plt.subplots(figsize=(3.5, 3.5 / 1.3333), dpi=300)
labels = {(8, 14): "Si-O", (5, 8): "B-O", (8, 11): "Na-O"}
for (zi, zj), g in rdfs.items():
    ax.plot(r, g, label=labels.get((zi, zj), f"{zi}-{zj}"))
ax.axhline(1, ls="--", color="gray", lw=0.7)
ax.set_xlabel("r (Å)")
ax.set_ylabel("g(r)")
ax.set_xlim(0, 5)
ax.set_ylim(0, None)
ax.set_title("Partial RDFs — quenched glass")
ax.legend()
plt.tight_layout()
plt.show()

OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

7. Quick Reference¶

Functions¶

Function	Signature	Returns
`get_structure_dict`	`(composition, target_atoms, mode="molar")`	`dict` with `"atoms"`, `"box"`
`get_ase_structure`	`(atoms_dict, replicate=(1,1,1))`	`ase.Atoms`
`generate_potential`	`(atoms_dict, potential_type="pmmcs", melt=False)`	`pd.DataFrame`
`melt_quench_simulation`	`(structure, potential, *, temperature_high, temperature_low, heating_rate, cooling_rate, ...)`	`{"structure": Atoms, "result": list[dict \\| None]}`
`frames_from_melt_quench_result`	`(result, initial_structure, *, stage=-1, stride=1)`	`list[Atoms]` — trajectory frames from one protocol stage

Potential quick-pick¶

B or V in system?          → bmp-screened-harmonic  (three-body covers B–O–B, B–O–Si, V–O–*)
B-free, needs P or broad?  → bmp-harmonic           (or pmmcs if transition metals needed)
Transition metals (no B)?  → pmmcs                  (check element coverage)
Strictly Ca/Al/Si/O?       → bjp                    (or pmmcs — both work)
Everything else?           → pmmcs

BMP boron restriction cheat-sheet¶

Composition type	Allowed?	Why
Na₂O-B₂O₃-SiO₂	✓	alkali borosilicate — Dell-Bray valid
CaO-B₂O₃-SiO₂	✓	alkaline-earth borosilicate — Dell-Bray valid
Na₂O-CaO-B₂O₃-SiO₂	✓	mixed modifier borosilicate — Dell-Bray valid
Al₂O₃-SiO₂ (no B)	✓	B absent — no restriction
Na₂O-Al₂O₃-B₂O₃-SiO₂	✗	Al + B — Dell-Bray invalid (Al not a modifier)
Na₂O-P₂O₅-B₂O₃	✗	P + B — Dell-Bray doesn't account for P as former

Output structure¶

result  = melt_quench_simulation(...)
glass   = result["structure"]   # ASE Atoms — pass to analysis functions
history = result["result"]      # list[dict | None] — one thermo dict per protocol stage

# Extract trajectory frames from the final stage for frame-averaged analysis
frames  = am.frames_from_melt_quench_result(result, atoms)  # list[Atoms]

For custom multi-stage cooling or annealing schedules, see MDTutorial.ipynb.