Projected RDF Tutorial¶

The projected RDF decomposes the pair correlation $g(r)$ into orientational components using $l=2$ spherical harmonics. This notebook shows:

1. What the $Y_2^0$ (uniaxial) and $Y_2^2$ (shear) lobes look like in 3D
2. How to interpret the sign of $g_\text{uniaxial}(r)$ and $g_\text{shear}(r)$
3. A worked example: 375 000-atom SiO₂ glass under 15 % uniaxial compression along $x$

In [1]:

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.io as pio

pio.renderers.default = "notebook"

# ── shared plot theme ────────────────────────────────────────────────────────
FONT_FAMILY = "Inter, Arial, sans-serif"
COLORS = {
    "bg": "#FAFAFA",
    "grid": "#E8E8E8",
    "zero": "#AAAAAA",
    "SiO": "#2563EB",  # blue
    "SiSi": "#DC2626",  # red
    "OO": "#16A34A",  # green
    "pos": "#DC2626",
    "neg": "#2563EB",
    "iso": "#6B7280",
    "strained": "#DC2626",
}


def apply_theme(fig, title="", xaxis_title="", yaxis_title="", height=480, legend_pos=None):
    legend_pos = legend_pos or {"x": 0.02, "y": 0.98, "xanchor": "left", "yanchor": "top"}
    fig.update_layout(
        title={
            "text": title,
            "font": {"size": 16, "family": FONT_FAMILY, "color": "#111827"},
            "x": 0.5,
            "xanchor": "center",
        },
        font={"family": FONT_FAMILY, "size": 13, "color": "#374151"},
        paper_bgcolor=COLORS["bg"],
        plot_bgcolor=COLORS["bg"],
        xaxis={
            "title": xaxis_title,
            "gridcolor": COLORS["grid"],
            "linecolor": "#D1D5DB",
            "mirror": True,
            "ticks": "outside",
            "ticklen": 4,
        },
        yaxis={
            "title": yaxis_title,
            "gridcolor": COLORS["grid"],
            "linecolor": "#D1D5DB",
            "mirror": True,
            "ticks": "outside",
            "ticklen": 4,
            "zerolinecolor": COLORS["zero"],
            "zerolinewidth": 1.5,
        },
        legend={**legend_pos, "bgcolor": "rgba(255,255,255,0.85)", "bordercolor": "#D1D5DB", "borderwidth": 1},
        height=height,
        margin={"l": 70, "r": 30, "t": 60, "b": 60},
    )


def apply_3d_theme(fig, title="", height=520):
    scene_style = {
        "xaxis": {
            "backgroundcolor": COLORS["bg"],
            "gridcolor": COLORS["grid"],
            "showbackground": True,
            "title_font": {"family": FONT_FAMILY},
        },
        "yaxis": {
            "backgroundcolor": COLORS["bg"],
            "gridcolor": COLORS["grid"],
            "showbackground": True,
            "title_font": {"family": FONT_FAMILY},
        },
        "zaxis": {
            "backgroundcolor": COLORS["bg"],
            "gridcolor": COLORS["grid"],
            "showbackground": True,
            "title_font": {"family": FONT_FAMILY},
        },
        "aspectmode": "cube",
    }
    fig.update_layout(
        title={
            "text": title,
            "font": {"size": 15, "family": FONT_FAMILY, "color": "#111827"},
            "x": 0.5,
            "xanchor": "center",
        },
        font={"family": FONT_FAMILY, "size": 12},
        paper_bgcolor=COLORS["bg"],
        scene=scene_style,
        height=height,
        margin={"l": 0, "r": 0, "t": 50, "b": 0},
    )

import numpy as np import plotly.graph_objects as go from plotly.subplots import make_subplots import plotly.io as pio pio.renderers.default = "notebook" # ── shared plot theme ──────────────────────────────────────────────────────── FONT_FAMILY = "Inter, Arial, sans-serif" COLORS = { "bg": "#FAFAFA", "grid": "#E8E8E8", "zero": "#AAAAAA", "SiO": "#2563EB", # blue "SiSi": "#DC2626", # red "OO": "#16A34A", # green "pos": "#DC2626", "neg": "#2563EB", "iso": "#6B7280", "strained": "#DC2626", } def apply_theme(fig, title="", xaxis_title="", yaxis_title="", height=480, legend_pos=None): legend_pos = legend_pos or {"x": 0.02, "y": 0.98, "xanchor": "left", "yanchor": "top"} fig.update_layout( title={ "text": title, "font": {"size": 16, "family": FONT_FAMILY, "color": "#111827"}, "x": 0.5, "xanchor": "center", }, font={"family": FONT_FAMILY, "size": 13, "color": "#374151"}, paper_bgcolor=COLORS["bg"], plot_bgcolor=COLORS["bg"], xaxis={ "title": xaxis_title, "gridcolor": COLORS["grid"], "linecolor": "#D1D5DB", "mirror": True, "ticks": "outside", "ticklen": 4, }, yaxis={ "title": yaxis_title, "gridcolor": COLORS["grid"], "linecolor": "#D1D5DB", "mirror": True, "ticks": "outside", "ticklen": 4, "zerolinecolor": COLORS["zero"], "zerolinewidth": 1.5, }, legend={**legend_pos, "bgcolor": "rgba(255,255,255,0.85)", "bordercolor": "#D1D5DB", "borderwidth": 1}, height=height, margin={"l": 70, "r": 30, "t": 60, "b": 60}, ) def apply_3d_theme(fig, title="", height=520): scene_style = { "xaxis": { "backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True, "title_font": {"family": FONT_FAMILY}, }, "yaxis": { "backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True, "title_font": {"family": FONT_FAMILY}, }, "zaxis": { "backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True, "title_font": {"family": FONT_FAMILY}, }, "aspectmode": "cube", } fig.update_layout( title={ "text": title, "font": {"size": 15, "family": FONT_FAMILY, "color": "#111827"}, "x": 0.5, "xanchor": "center", }, font={"family": FONT_FAMILY, "size": 12}, paper_bgcolor=COLORS["bg"], scene=scene_style, height=height, margin={"l": 0, "r": 0, "t": 50, "b": 0}, )

1 — Spherical harmonic lobes¶

$Y_2^0$ — uniaxial deformation signal¶

$$Y_2^0(\theta) = \sqrt{\frac{5}{4\pi}} \left(\frac{3}{2}\cos^2\theta - \frac{1}{2}\right)$$

The lobe is a dumbbell along $z$: positive at the poles, negative at the equator. When pairs preferentially align along $z$ (compression), the first-shell peak in $g_\text{uniaxial}(r)$ becomes positive. When they tilt away (transverse expansion), it becomes negative.

In [2]:

theta = np.linspace(0, np.pi, 150)
phi = np.linspace(0, 2 * np.pi, 150)
THETA, PHI = np.meshgrid(theta, phi)

norm_Y20 = np.sqrt(5.0 / (4.0 * np.pi))
Y20 = norm_Y20 * (1.5 * np.cos(THETA) ** 2 - 0.5)
R = np.abs(Y20)
X = R * np.sin(THETA) * np.cos(PHI)
Y = R * np.sin(THETA) * np.sin(PHI)
Z = R * np.cos(THETA)


def box_edges(scale=1.0):
    s = scale
    corners = np.array(
        [
            [-s, -s, -s],
            [s, -s, -s],
            [s, s, -s],
            [-s, s, -s],
            [-s, -s, -s],
            [-s, -s, s],
            [s, -s, s],
            [s, s, s],
            [-s, s, s],
            [-s, -s, s],
            [None, None, None],
            [s, -s, -s],
            [s, -s, s],
            [None, None, None],
            [s, s, -s],
            [s, s, s],
            [None, None, None],
            [-s, s, -s],
            [-s, s, s],
        ]
    )
    return corners[:, 0], corners[:, 1], corners[:, 2]


bx, by, bz = box_edges()

fig_y20 = go.Figure()
fig_y20.add_trace(
    go.Surface(
        x=X,
        y=Y,
        z=Z,
        surfacecolor=Y20,
        colorscale="RdBu",
        cmid=0,
        cmin=-0.6,
        cmax=0.6,
        colorbar={"title": {"text": "Y₂⁰", "side": "right"}, "thickness": 14, "len": 0.6},
        opacity=0.92,
        name="Y₂⁰",
        hovertemplate="Y₂⁰ = %{surfacecolor:.3f}<extra></extra>",
    )
)
fig_y20.add_trace(
    go.Scatter3d(
        x=bx,
        y=by,
        z=bz,
        mode="lines",
        line={"color": "#6B7280", "width": 2},
        name="simulation box",
    )
)
fig_y20.add_trace(
    go.Cone(
        x=[0],
        y=[0],
        z=[0.65],
        u=[0],
        v=[0],
        w=[0.3],
        sizemode="absolute",
        sizeref=0.12,
        colorscale=[[0, "#FF0000"], [1, "#FF0000"]],
        showscale=False,
        name="deformation axis (z)",
    )
)
apply_3d_theme(fig_y20, title="Y₂⁰ — uniaxial signal (deformation axis = z)")
fig_y20.update_layout(
    scene={
        **fig_y20.layout.scene.to_plotly_json(),
        "xaxis_title": "x",
        "yaxis_title": "y",
        "zaxis_title": "z",
        "annotations": [
            {
                "x": 0,
                "y": 0,
                "z": 0.88,
                "text": "ε<sub>zz</sub>",
                "showarrow": False,
                "font": {"size": 20, "color": "#FF0000", "family": FONT_FAMILY},
            }
        ],
    }
)
fig_y20.show()

theta = np.linspace(0, np.pi, 150) phi = np.linspace(0, 2 * np.pi, 150) THETA, PHI = np.meshgrid(theta, phi) norm_Y20 = np.sqrt(5.0 / (4.0 * np.pi)) Y20 = norm_Y20 * (1.5 * np.cos(THETA) ** 2 - 0.5) R = np.abs(Y20) X = R * np.sin(THETA) * np.cos(PHI) Y = R * np.sin(THETA) * np.sin(PHI) Z = R * np.cos(THETA) def box_edges(scale=1.0): s = scale corners = np.array( [ [-s, -s, -s], [s, -s, -s], [s, s, -s], [-s, s, -s], [-s, -s, -s], [-s, -s, s], [s, -s, s], [s, s, s], [-s, s, s], [-s, -s, s], [None, None, None], [s, -s, -s], [s, -s, s], [None, None, None], [s, s, -s], [s, s, s], [None, None, None], [-s, s, -s], [-s, s, s], ] ) return corners[:, 0], corners[:, 1], corners[:, 2] bx, by, bz = box_edges() fig_y20 = go.Figure() fig_y20.add_trace( go.Surface( x=X, y=Y, z=Z, surfacecolor=Y20, colorscale="RdBu", cmid=0, cmin=-0.6, cmax=0.6, colorbar={"title": {"text": "Y₂⁰", "side": "right"}, "thickness": 14, "len": 0.6}, opacity=0.92, name="Y₂⁰", hovertemplate="Y₂⁰ = %{surfacecolor:.3f}", ) ) fig_y20.add_trace( go.Scatter3d( x=bx, y=by, z=bz, mode="lines", line={"color": "#6B7280", "width": 2}, name="simulation box", ) ) fig_y20.add_trace( go.Cone( x=[0], y=[0], z=[0.65], u=[0], v=[0], w=[0.3], sizemode="absolute", sizeref=0.12, colorscale=[[0, "#FF0000"], [1, "#FF0000"]], showscale=False, name="deformation axis (z)", ) ) apply_3d_theme(fig_y20, title="Y₂⁰ — uniaxial signal (deformation axis = z)") fig_y20.update_layout( scene={ **fig_y20.layout.scene.to_plotly_json(), "xaxis_title": "x", "yaxis_title": "y", "zaxis_title": "z", "annotations": [ { "x": 0, "y": 0, "z": 0.88, "text": "ε_zz", "showarrow": False, "font": {"size": 20, "color": "#FF0000", "family": FONT_FAMILY}, } ], } ) fig_y20.show()

$Y_2^2$ shear lobes — one per shear plane¶

Each shear plane maps to one $Y_2^m$ component whose Cartesian form is the off-diagonal product of that plane's axes:

Shear plane	Component	Cartesian form
$xy$	$\text{Im}\,Y_2^2$	$\propto \hat{u}_x \hat{u}_y$
$xz$	$\text{Re}\,Y_2^1$	$\propto \hat{u}_x \hat{u}_z$
$yz$	$\text{Im}\,Y_2^1$	$\propto \hat{u}_y \hat{u}_z$

In [3]:

norm_Y21 = np.sqrt(15.0 / (8.0 * np.pi))
sin_t = np.sin(THETA)
cos_t = np.cos(THETA)
ux = sin_t * np.cos(PHI)
uy = sin_t * np.sin(PHI)
uz = cos_t

components = {
    "xy  (Im Y₂²)": norm_Y21 * ux * uy,
    "xz  (Re Y₂¹)": norm_Y21 * ux * uz,
    "yz  (Im Y₂¹)": norm_Y21 * uy * uz,
}

fig_shear = make_subplots(
    rows=3,
    cols=1,
    specs=[[{"type": "surface"}], [{"type": "surface"}], [{"type": "surface"}]],
    subplot_titles=list(components.keys()),
)

SHEAR_CS = "RdBu"
for row, (_label, Ylm) in enumerate(components.items(), start=1):
    R_sh = np.abs(Ylm)
    Xs = R_sh * np.sin(THETA) * np.cos(PHI)
    Ys = R_sh * np.sin(THETA) * np.sin(PHI)
    Zs = R_sh * np.cos(THETA)
    fig_shear.add_trace(
        go.Surface(
            x=Xs,
            y=Ys,
            z=Zs,
            surfacecolor=Ylm,
            colorscale=SHEAR_CS,
            cmid=0,
            cmin=-0.4,
            cmax=0.4,
            showscale=(row == 3),
            colorbar={"title": {"text": "value", "side": "right"}, "thickness": 12, "x": 1.02},
            opacity=0.92,
            hovertemplate="Y_lm = %{surfacecolor:.3f}<extra></extra>",
        ),
        row=row,
        col=1,
    )
    fig_shear.add_trace(
        go.Scatter3d(
            x=bx,
            y=by,
            z=bz,
            mode="lines",
            line={"color": "#6B7280", "width": 2},
            showlegend=(row == 1),
            name="simulation box",
        ),
        row=row,
        col=1,
    )

fig_shear.update_layout(
    title={
        "text": "Shear Y_lm lobes — one per shear plane",
        "font": {"size": 15, "family": FONT_FAMILY, "color": "#111827"},
        "x": 0.5,
        "xanchor": "center",
    },
    font={"family": FONT_FAMILY, "size": 12},
    paper_bgcolor=COLORS["bg"],
    height=1400,
    margin={"l": 0, "r": 0, "t": 60, "b": 0},
)
for i in range(1, 4):
    fig_shear.update_scenes(
        {
            "xaxis_title": "x",
            "yaxis_title": "y",
            "zaxis_title": "z",
            "aspectmode": "cube",
            "xaxis": {"backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True},
            "yaxis": {"backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True},
            "zaxis": {"backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True},
        },
        row=i,
        col=1,
    )
fig_shear.show()

norm_Y21 = np.sqrt(15.0 / (8.0 * np.pi)) sin_t = np.sin(THETA) cos_t = np.cos(THETA) ux = sin_t * np.cos(PHI) uy = sin_t * np.sin(PHI) uz = cos_t components = { "xy (Im Y₂²)": norm_Y21 * ux * uy, "xz (Re Y₂¹)": norm_Y21 * ux * uz, "yz (Im Y₂¹)": norm_Y21 * uy * uz, } fig_shear = make_subplots( rows=3, cols=1, specs=[[{"type": "surface"}], [{"type": "surface"}], [{"type": "surface"}]], subplot_titles=list(components.keys()), ) SHEAR_CS = "RdBu" for row, (_label, Ylm) in enumerate(components.items(), start=1): R_sh = np.abs(Ylm) Xs = R_sh * np.sin(THETA) * np.cos(PHI) Ys = R_sh * np.sin(THETA) * np.sin(PHI) Zs = R_sh * np.cos(THETA) fig_shear.add_trace( go.Surface( x=Xs, y=Ys, z=Zs, surfacecolor=Ylm, colorscale=SHEAR_CS, cmid=0, cmin=-0.4, cmax=0.4, showscale=(row == 3), colorbar={"title": {"text": "value", "side": "right"}, "thickness": 12, "x": 1.02}, opacity=0.92, hovertemplate="Y_lm = %{surfacecolor:.3f}", ), row=row, col=1, ) fig_shear.add_trace( go.Scatter3d( x=bx, y=by, z=bz, mode="lines", line={"color": "#6B7280", "width": 2}, showlegend=(row == 1), name="simulation box", ), row=row, col=1, ) fig_shear.update_layout( title={ "text": "Shear Y_lm lobes — one per shear plane", "font": {"size": 15, "family": FONT_FAMILY, "color": "#111827"}, "x": 0.5, "xanchor": "center", }, font={"family": FONT_FAMILY, "size": 12}, paper_bgcolor=COLORS["bg"], height=1400, margin={"l": 0, "r": 0, "t": 60, "b": 0}, ) for i in range(1, 4): fig_shear.update_scenes( { "xaxis_title": "x", "yaxis_title": "y", "zaxis_title": "z", "aspectmode": "cube", "xaxis": {"backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True}, "yaxis": {"backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True}, "zaxis": {"backgroundcolor": COLORS["bg"], "gridcolor": COLORS["grid"], "showbackground": True}, }, row=i, col=1, ) fig_shear.show()

2 — Interpretation guide¶

The plot below shows a schematic $g_\text{uniaxial}(r)$:

• Blue (isotropic) — no deformation, signal flat at zero.
• Red (z-strained) — first-shell peak goes positive (bonds align along compression axis); subsequent shells may show oscillating sign as the structure reorganises.

In [4]:

r_mock = np.linspace(1.5, 8.0, 1200)

g_iso = np.zeros_like(r_mock)
g_strained = 0.20 * np.exp(-((r_mock - 2.35) ** 2) / (2 * 0.06**2)) - 0.10 * np.exp(
    -((r_mock - 2.52) ** 2) / (2 * 0.06**2)
)

fig_interp = go.Figure()
fig_interp.add_trace(
    go.Scatter(
        x=r_mock,
        y=g_iso,
        mode="lines",
        name="isotropic (reference)",
        line={"color": COLORS["iso"], "width": 2, "dash": "dot"},
    )
)
fig_interp.add_trace(
    go.Scatter(
        x=r_mock,
        y=g_strained,
        mode="lines",
        name="z-strained",
        line={"color": COLORS["strained"], "width": 2.5},
        hovertemplate="r = %{x:.2f} Å<br>g_uniaxial = %{y:.4f}<extra></extra>",
    )
)
fig_interp.add_annotation(
    x=2.35,
    y=0.20,
    text="bonds align <b>along</b> z",
    showarrow=True,
    arrowhead=2,
    arrowcolor=COLORS["strained"],
    ax=55,
    ay=-35,
    font={"color": COLORS["strained"], "size": 12, "family": FONT_FAMILY},
)
fig_interp.add_annotation(
    x=2.52,
    y=-0.10,
    text="bonds tilt <b>away</b> from z",
    showarrow=True,
    arrowhead=2,
    arrowcolor="#2563EB",
    ax=-60,
    ay=35,
    font={"color": "#2563EB", "size": 12, "family": FONT_FAMILY},
)
apply_theme(
    fig_interp,
    title="g_uniaxial(r) — sign convention (schematic)",
    xaxis_title="r (Å)",
    yaxis_title="g_uniaxial(r)",
    height=420,
    legend_pos={"x": 0.72, "y": 0.98, "xanchor": "left", "yanchor": "top"},
)
fig_interp.show()

r_mock = np.linspace(1.5, 8.0, 1200) g_iso = np.zeros_like(r_mock) g_strained = 0.20 * np.exp(-((r_mock - 2.35) ** 2) / (2 * 0.06**2)) - 0.10 * np.exp( -((r_mock - 2.52) ** 2) / (2 * 0.06**2) ) fig_interp = go.Figure() fig_interp.add_trace( go.Scatter( x=r_mock, y=g_iso, mode="lines", name="isotropic (reference)", line={"color": COLORS["iso"], "width": 2, "dash": "dot"}, ) ) fig_interp.add_trace( go.Scatter( x=r_mock, y=g_strained, mode="lines", name="z-strained", line={"color": COLORS["strained"], "width": 2.5}, hovertemplate="r = %{x:.2f} Å
g_uniaxial = %{y:.4f}", ) ) fig_interp.add_annotation( x=2.35, y=0.20, text="bonds align along z", showarrow=True, arrowhead=2, arrowcolor=COLORS["strained"], ax=55, ay=-35, font={"color": COLORS["strained"], "size": 12, "family": FONT_FAMILY}, ) fig_interp.add_annotation( x=2.52, y=-0.10, text="bonds tilt away from z", showarrow=True, arrowhead=2, arrowcolor="#2563EB", ax=-60, ay=35, font={"color": "#2563EB", "size": 12, "family": FONT_FAMILY}, ) apply_theme( fig_interp, title="g_uniaxial(r) — sign convention (schematic)", xaxis_title="r (Å)", yaxis_title="g_uniaxial(r)", height=420, legend_pos={"x": 0.72, "y": 0.98, "xanchor": "left", "yanchor": "top"}, ) fig_interp.show()

3 — Worked example: SiO₂ glass under uniaxial compression¶

Two LAMMPS snapshots from a 300 K compression simulation using a strain rate of $10^8$ 1/s:

Snapshot	Timestep	$L_x$ (Å)	$L_y$ (Å)	$L_z$ (Å)	$\varepsilon_{xx}$
Pristine	0	177.88	177.88	177.88	0 %
Compressed	1 500 000	151.20	182.10	182.13	-15 %

375 000 atoms per frame gives shot noise ≈ 0.007 per bin — small enough to resolve the ~0.1 signal from the deformed Si-O network unambiguously.

Deformation axis is x → we compute deformation_axis='x'.

In [5]:

from ase.io import read
from amorphouspy import compute_projected_rdf

pristine = read("data/SiO2_compression_X.0.extxyz")
compressed = read("data/SiO2_compression_X.1500000.extxyz")

print(f"Pristine   {len(pristine)} atoms | box: {pristine.get_cell().array.diagonal().round(2)} Å")
print(f"Compressed {len(compressed)} atoms | box: {compressed.get_cell().array.diagonal().round(2)} Å")

r, uni_0, _ = compute_projected_rdf(pristine, deformation_axis="x", r_max=10.0, n_bins=1000)
r, uni_2, _ = compute_projected_rdf(compressed, deformation_axis="x", r_max=10.0, n_bins=1000)

from ase.io import read from amorphouspy import compute_projected_rdf pristine = read("data/SiO2_compression_X.0.extxyz") compressed = read("data/SiO2_compression_X.1500000.extxyz") print(f"Pristine {len(pristine)} atoms | box: {pristine.get_cell().array.diagonal().round(2)} Å") print(f"Compressed {len(compressed)} atoms | box: {compressed.get_cell().array.diagonal().round(2)} Å") r, uni_0, _ = compute_projected_rdf(pristine, deformation_axis="x", r_max=10.0, n_bins=1000) r, uni_2, _ = compute_projected_rdf(compressed, deformation_axis="x", r_max=10.0, n_bins=1000)

Pristine   375000 atoms | box: [177.88 177.88 177.88] Å
Compressed 375000 atoms | box: [151.2  182.1  182.13] Å

In [6]:

import numpy as np

pair_styles = {
    (8, 14): {"name": "O-Si", "color": COLORS["SiO"]},
    (14, 14): {"name": "Si-Si", "color": COLORS["SiSi"]},
    (8, 8): {"name": "O-O", "color": COLORS["OO"]},
}
pairs_ordered = [(8, 14), (14, 14), (8, 8)]

# 2 rows × 3 cols: top = raw g_uni, bottom = Δg = compressed - pristine
fig_comp = make_subplots(
    rows=2,
    cols=3,
    shared_xaxes=True,
    row_heights=[2, 1],
    subplot_titles=[
        "O-Si — g_uniaxial",
        "Si-Si — g_uniaxial",
        "O-O — g_uniaxial",
        "O-Si — Δg (compressed - pristine)",
        "Si-Si — Δg",
        "O-O — Δg",
    ],
    vertical_spacing=0.08,
    horizontal_spacing=0.08,
)

for col, pair in enumerate(pairs_ordered, start=1):
    color = pair_styles[pair]["color"]
    g0 = uni_0[pair]
    g2 = uni_2[pair]
    dg = g2 - g0

    # --- top row: raw curves ---
    fig_comp.add_trace(
        go.Scatter(
            x=r,
            y=g0,
            mode="lines",
            name="pristine" if col == 1 else None,
            showlegend=(col == 1),
            legendgroup="pristine",
            line={"color": "#9CA3AF", "width": 1.2},
            hovertemplate="r=%{x:.3f} Å  g=%{y:.4f}<extra>pristine</extra>",
        ),
        row=1,
        col=col,
    )
    fig_comp.add_trace(
        go.Scatter(
            x=r,
            y=g2,
            mode="lines",
            name="compressed (ε_xx=-15%)" if col == 1 else None,
            showlegend=(col == 1),
            legendgroup="compressed",
            line={"color": color, "width": 1.5},
            hovertemplate="r=%{x:.3f} Å  g=%{y:.4f}<extra>compressed</extra>",
        ),
        row=1,
        col=col,
    )
    fig_comp.add_hline(y=0, line_dash="dot", line_color=COLORS["zero"], line_width=1, row=1, col=col)

    # --- bottom row: difference ---
    fig_comp.add_trace(
        go.Scatter(
            x=r,
            y=dg,
            mode="lines",
            name="Δg" if col == 1 else None,
            showlegend=(col == 1),
            legendgroup="dg",
            line={"color": color, "width": 1.5},
            hovertemplate="r=%{x:.3f} Å  Δg=%{y:.4f}<extra>Δg</extra>",
        ),
        row=2,
        col=col,
    )
    fig_comp.add_hline(y=0, line_dash="dot", line_color=COLORS["zero"], line_width=1, row=2, col=col)

fig_comp.update_layout(
    title={
        "text": "Uniaxial projected RDF — SiO₂ glass, ε_xx = -15 % (375 000 atoms)",
        "font": {"size": 15, "family": FONT_FAMILY, "color": "#111827"},
        "x": 0.5,
        "xanchor": "center",
    },
    font={"family": FONT_FAMILY, "size": 12, "color": "#374151"},
    paper_bgcolor=COLORS["bg"],
    plot_bgcolor=COLORS["bg"],
    height=700,
    legend={
        "x": 0.27,
        "y": 0.41,
        "xanchor": "right",
        "yanchor": "bottom",
        "bgcolor": "rgba(255,255,255,0.88)",
        "bordercolor": "#D1D5DB",
        "borderwidth": 1,
    },
    margin={"l": 60, "r": 30, "t": 80, "b": 60},
)
for row in range(1, 3):
    for col in range(1, 4):
        fig_comp.update_xaxes(
            gridcolor=COLORS["grid"],
            linecolor="#D1D5DB",
            mirror=True,
            ticks="outside",
            ticklen=4,
            row=row,
            col=col,
        )
        ylabel = "g_uni(r)" if row == 1 else "Δg(r)"
        fig_comp.update_yaxes(
            title_text=ylabel,
            gridcolor=COLORS["grid"],
            linecolor="#D1D5DB",
            zerolinecolor=COLORS["zero"],
            zerolinewidth=1.5,
            mirror=True,
            ticks="outside",
            ticklen=4,
            row=row,
            col=col,
        )
for col in range(1, 4):
    fig_comp.update_xaxes(title_text="r (Å)", row=2, col=col)
fig_comp.show()

import numpy as np pair_styles = { (8, 14): {"name": "O-Si", "color": COLORS["SiO"]}, (14, 14): {"name": "Si-Si", "color": COLORS["SiSi"]}, (8, 8): {"name": "O-O", "color": COLORS["OO"]}, } pairs_ordered = [(8, 14), (14, 14), (8, 8)] # 2 rows × 3 cols: top = raw g_uni, bottom = Δg = compressed - pristine fig_comp = make_subplots( rows=2, cols=3, shared_xaxes=True, row_heights=[2, 1], subplot_titles=[ "O-Si — g_uniaxial", "Si-Si — g_uniaxial", "O-O — g_uniaxial", "O-Si — Δg (compressed - pristine)", "Si-Si — Δg", "O-O — Δg", ], vertical_spacing=0.08, horizontal_spacing=0.08, ) for col, pair in enumerate(pairs_ordered, start=1): color = pair_styles[pair]["color"] g0 = uni_0[pair] g2 = uni_2[pair] dg = g2 - g0 # --- top row: raw curves --- fig_comp.add_trace( go.Scatter( x=r, y=g0, mode="lines", name="pristine" if col == 1 else None, showlegend=(col == 1), legendgroup="pristine", line={"color": "#9CA3AF", "width": 1.2}, hovertemplate="r=%{x:.3f} Å g=%{y:.4f}pristine", ), row=1, col=col, ) fig_comp.add_trace( go.Scatter( x=r, y=g2, mode="lines", name="compressed (ε_xx=-15%)" if col == 1 else None, showlegend=(col == 1), legendgroup="compressed", line={"color": color, "width": 1.5}, hovertemplate="r=%{x:.3f} Å g=%{y:.4f}compressed", ), row=1, col=col, ) fig_comp.add_hline(y=0, line_dash="dot", line_color=COLORS["zero"], line_width=1, row=1, col=col) # --- bottom row: difference --- fig_comp.add_trace( go.Scatter( x=r, y=dg, mode="lines", name="Δg" if col == 1 else None, showlegend=(col == 1), legendgroup="dg", line={"color": color, "width": 1.5}, hovertemplate="r=%{x:.3f} Å Δg=%{y:.4f}Δg", ), row=2, col=col, ) fig_comp.add_hline(y=0, line_dash="dot", line_color=COLORS["zero"], line_width=1, row=2, col=col) fig_comp.update_layout( title={ "text": "Uniaxial projected RDF — SiO₂ glass, ε_xx = -15 % (375 000 atoms)", "font": {"size": 15, "family": FONT_FAMILY, "color": "#111827"}, "x": 0.5, "xanchor": "center", }, font={"family": FONT_FAMILY, "size": 12, "color": "#374151"}, paper_bgcolor=COLORS["bg"], plot_bgcolor=COLORS["bg"], height=700, legend={ "x": 0.27, "y": 0.41, "xanchor": "right", "yanchor": "bottom", "bgcolor": "rgba(255,255,255,0.88)", "bordercolor": "#D1D5DB", "borderwidth": 1, }, margin={"l": 60, "r": 30, "t": 80, "b": 60}, ) for row in range(1, 3): for col in range(1, 4): fig_comp.update_xaxes( gridcolor=COLORS["grid"], linecolor="#D1D5DB", mirror=True, ticks="outside", ticklen=4, row=row, col=col, ) ylabel = "g_uni(r)" if row == 1 else "Δg(r)" fig_comp.update_yaxes( title_text=ylabel, gridcolor=COLORS["grid"], linecolor="#D1D5DB", zerolinecolor=COLORS["zero"], zerolinewidth=1.5, mirror=True, ticks="outside", ticklen=4, row=row, col=col, ) for col in range(1, 4): fig_comp.update_xaxes(title_text="r (Å)", row=2, col=col) fig_comp.show()

What the data tell us¶

The top sub-panel for each pair shows the raw $g_\text{uniaxial}(r)$ for pristine (grey) and compressed (colour). The pristine curve is not a clean zero baseline — it is shot noise from the finite snapshot.

The bottom sub-panel shows the difference $\Delta g = g_\text{uni}^\text{comp} - g_\text{uni}^\text{pristine}$, which cancels the shared structural noise and isolates the deformation signal.

O-Si — strongest signal. The difference troughs at −3.09 @ 1.63 Å, meaning Si-O bonds preferentially orient perpendicular to the compression axis $x$ under −15 % strain. This is the expected response: uniaxial compression along $x$ pushes atoms transversely apart, tilting Si-O bonds away from $x$, so $g_{20}(x) < 0$ at the first-shell distance (~1.61 Å). The signal rms is 14.4× the pristine noise floor.

Si-Si — trough at −1.05 @ 3.14 Å. The Si-Si vector connects centres of corner-sharing SiO₄ tetrahedra through the bridging oxygen. A negative signal means these Si-Si vectors also tilt away from $x$, consistent with the Si-O bond reorientation: if Si-O bonds tilt away from $x$, the Si-Si vectors along the Si-O-Si bridge follow. The difference rms is 10.8× the noise floor.

O-O — antisymmetric doublet: positive peak at +1.43 @ 2.46 Å and negative trough at −1.21 @ 2.68 Å, straddling the pristine O-O first-shell peak at ~2.61 Å. This shape is the fingerprint of intra-tetrahedral distortion. Under compression along $x$, SiO₄ tetrahedra flatten: O-O pairs whose vector is already roughly parallel to $x$ get squeezed to shorter distances (~2.46 Å), while those pointing transversely are stretched to longer distances (~2.68 Å). The positive lobe at short $r$ and negative lobe at long $r$ are exactly the antisymmetric doublet $\Delta g \approx -\delta r \cdot dg/dr$ expected from a peak shift $\delta r < 0$ (shift to shorter $r$).

This is physically distinct from the O-Si signal: Si-O bond lengths barely change (the O-Si doublet is a single trough, not a doublet), the bonds simply reorient. For O-O you are seeing the distortion of the O-Si-O bond angle inside the tetrahedron. The fact that the positive lobe is slightly larger than the negative indicates a net compression of O-O distances along $x$ on top of the reorientation — the tetrahedra are not only tilting but also slightly compressing in volume. The difference rms is 25.6× the noise floor, the largest SNR of the three pairs.

Pair	$\Delta g$ peak	at $r$ (Å)	rms $\Delta g$ / rms noise	Interpretation
O-Si	−3.09	1.63	14.4×	Si-O bonds tilt away from compression axis
Si-Si	−1.05	3.14	10.8×	Si-Si vectors follow Si-O reorientation
O-O	+1.43 / −1.21	2.46 / 2.68	25.6×	Intra-tetrahedral distortion: O-O peak shift along $x$

In [7]:

first_shell_ranges = {(8, 14): (1.4, 2.2), (14, 14): (2.5, 4.5), (8, 8): (2.0, 3.5)}
pair_labels = {(8, 14): "O-Si", (14, 14): "Si-Si", (8, 8): "O-O"}

print(f"{'Pair':<8} {'noise (pristine rms)':>21} {'Δg rms':>8} {'SNR':>6} {'Δg peak':>9} {'at r (Å)':>10}")
print("-" * 68)
for pair, (rlo, rhi) in first_shell_ranges.items():
    if pair not in uni_0:
        continue
    mask = (r > rlo) & (r < rhi)
    g0 = uni_0[pair]
    dg = uni_2[pair] - g0
    noise = float(np.sqrt(np.mean(g0**2)))
    rms_dg = float(np.sqrt(np.mean(dg**2)))
    peak_idx = np.argmax(np.abs(dg[mask]))
    peak_val = dg[mask][peak_idx]
    peak_r = r[mask][peak_idx]
    print(
        f"{pair_labels[pair]:<8} {noise:>21.5f} {rms_dg:>8.5f} {rms_dg / noise:>6.2f}x {peak_val:>9.4f} {peak_r:>10.3f}"
    )

first_shell_ranges = {(8, 14): (1.4, 2.2), (14, 14): (2.5, 4.5), (8, 8): (2.0, 3.5)} pair_labels = {(8, 14): "O-Si", (14, 14): "Si-Si", (8, 8): "O-O"} print(f"{'Pair':<8} {'noise (pristine rms)':>21} {'Δg rms':>8} {'SNR':>6} {'Δg peak':>9} {'at r (Å)':>10}") print("-" * 68) for pair, (rlo, rhi) in first_shell_ranges.items(): if pair not in uni_0: continue mask = (r > rlo) & (r < rhi) g0 = uni_0[pair] dg = uni_2[pair] - g0 noise = float(np.sqrt(np.mean(g0**2))) rms_dg = float(np.sqrt(np.mean(dg**2))) peak_idx = np.argmax(np.abs(dg[mask])) peak_val = dg[mask][peak_idx] peak_r = r[mask][peak_idx] print( f"{pair_labels[pair]:<8} {noise:>21.5f} {rms_dg:>8.5f} {rms_dg / noise:>6.2f}x {peak_val:>9.4f} {peak_r:>10.3f}" )

Pair      noise (pristine rms)   Δg rms    SNR   Δg peak   at r (Å)
--------------------------------------------------------------------
O-Si                   0.01621  0.25291  15.61x   -3.1046      1.625
Si-Si                  0.01287  0.15234  11.84x   -0.9823      3.135
O-O                    0.00721  0.21018  29.15x    1.3988      2.485

References¶

1. J. Zausch and J. Horbach, "The build-up and relaxation of stresses in a glass-forming soft-sphere mixture under shear: A computer simulation study," EPL 88, 60001 (2009). https://doi.org/10.1209/0295-5075/88/60001

2. S. Ganisetti, A. Atila, J. Guénolé, A. Prakash, J. Horbach, L. Wondraczek, and E. Bitzek, "The origin of deformation induced topological anisotropy in silica glass," Acta Mater. 257, 119108 (2023). https://doi.org/10.1016/j.actamat.2023.119108

3. A. Atila and E. Bitzek, "Atomistic origins of deformation-induced structural anisotropy in metaphosphate glasses and its influence on mechanical properties," J. Non-Cryst. Solids 627, 122822 (2024). https://doi.org/10.1016/j.jnoncrysol.2024.122822