Oxidation phase diagram of a metallic surface

This tutorial reproduces the canonical mcpy workflow: sweep the oxygen chemical potential \(\mu_{\mathrm{O}}\) across a range of environmental conditions, run an independent GCMC trajectory at each value, then assemble the resulting configurations into a surface oxidation phase diagram.

It is intentionally longer than Your first simulation — by the end you will have a publication-style figure, not just a log file.

Goal 

For an Ag(111) slab in contact with an oxygen reservoir, determine which oxide coverages are thermodynamically stable as a function of \(\mu_{\mathrm{O}}\), and locate the transitions between them.

The deliverable is a plot of the surface Gibbs energy \(\gamma\) vs. \(\mu_{\mathrm{O}}\), with the lower envelope (convex hull) marking the stable phases.

Prerequisites 

A working mcpy installation with a MACE checkpoint (see Installation).
A successful run of Your first simulation — confirms that the calculator and cell are wired correctly.
Calibrated species_radii for Ag and O on your potential (see Calibrating species_radii).

Inputs 

Setting	Value
Surface	fcc(111), 4×4×3 Ag slab, 8 Å vacuum
Reservoir species	O
Temperature	500 K
\(\mu_{\mathrm{O}}\)	sweep: −6.0, −5.5, −5.0, −4.5, −4.0 eV
GCMC steps per condition	2 000 (production), 500 (equilibration)
Move set	Insertion + Deletion + Displacement of O
Cell	`CustomCell` 5 Å above the top Ag layer

Workflow 

Per-condition GCMC. Loop over the target \(\mu_{\mathrm{O}}\) values; for each one start from the same clean slab, run equilibration followed by production, and store the trajectory and log under a condition-specific filename.
Reference energies. Compute the energy of the clean Ag slab \(E_{\mathrm{ref}}\) once, separately, and record the bulk Ag chemical potential \(\mu_{\mathrm{Ag}}\).
Phase-diagram assembly. For every frame written during production, compute the surface Gibbs energy

\[\gamma(\mu_{\mathrm{O}}) = \frac{E - n_{\mathrm{Ag}}\,\mu_{\mathrm{Ag}} - n_{\mathrm{O}}\,\mu_{\mathrm{O}}}{A} - \gamma_{\mathrm{ref}},\]

then take the lower envelope across all conditions.

Code 

Per-condition driver script:

for delta_mu in (-1.0, -0.5, 0.0, 0.5, 1.0):
    run_gcmc(
        delta_mu_gas=delta_mu,
        gcmc_steps=2000,
        outdir=f'sweep/dmu_{delta_mu:+.1f}',
    )

The body of run_gcmc is the same as in Your first simulation, parameterised on delta_mu_gas so the only difference between runs is \(\mu_{\mathrm{O}} = \mu_{\mathrm{O}}^{\mathrm{ref}} + \Delta\mu\).

Phase-diagram post-processing (sketch):

from mcpy.utils.phase_diagram import surface_gibbs

gammas = surface_gibbs(
    trajectories=glob('sweep/dmu_*/*.xyz'),
    mu_metal=-3.27,
    mu_gas_ref=-4.91,
    e_ref=E_clean_slab,
    area=slab_area_A2,
)
plot_lower_envelope(gammas)

Expected output 

You should obtain:

One gcmc_*.xyz and gcmc_*.out pair per \(\mu_{\mathrm{O}}\) value, under sweep/dmu_*/.
A phase-diagram figure phase_diagram.pdf, showing one curve per condition and a highlighted convex hull marking the stable phases.

A representative log line at the end of a converged run looks like:

2000       62          -174.55          ins: 22.3%, del: 10.1%, disp: 41.7%

Interpretation 

Each \(\mu_{\mathrm{O}}\) value should converge to a roughly constant \(N_{\mathrm{O}}\) in the second half of its run. If it does not, extend the production stage or add displacement moves.
Intersections between \(\gamma(\mu)\) curves locate the \(\mu_{\mathrm{O}}\) values at which the surface switches between oxide phases. Convert those values to \((T, p_{\mathrm{O}_2})\) with the standard relation given in Examples.
A single curve sitting above all others across the whole range usually means the structure is metastable and was never visited at equilibrium — discard it.

Common pitfalls 

Mis-calibrated radii — insertions all overlap; acceptance is near zero. Re-run the calibration in Calibrating species_radii.
Too-short equilibration — early frames bias \(\gamma\). Drop the first 25–50 % of every trajectory before assembling the phase diagram.
Reference energy mismatch — if \(E_{\mathrm{ref}}\) is computed with a different relaxation tolerance than the production runs, the curves all shift by a constant. Use the same calculator settings everywhere.

Next steps 

Run the sweep in parallel with Replica Exchange: see the ReplicaExchange section of Ensembles.
Apply the same recipe to a supported nanoparticle: see examples/gcmc_nano_supported.py.
Extend to a bimetallic reservoir by adding a second species to mu and registering the corresponding insertion/deletion moves.