💡 Tip: In Colab, run the first cell once to install dependencies, then restart the runtime before running the rest of the notebook. For local Jupyter, activate your synth-pdb environment first.

In [ ]:

Copied!





import os
import sys

print(f"🐍 Python: {sys.executable}")
try:
    import google.colab
    IN_COLAB = True
except ImportError:
    IN_COLAB = False

if IN_COLAB:
    if not os.path.exists("installed.marker"):
        import subprocess
        subprocess.run([sys.executable, "-m", "pip", "install", "-q",
                        "synth-pdb", "synth-nmr", "biotite",
                        "matplotlib", "ipywidgets", "numpy"])
        with open("installed.marker", "w") as f: f.write("done")
        print("✅ Installed. Please run this cell once more if imports fail.")
    else:
        print("✅ Already installed.")
else:
    print("✅ Running locally — no install needed.")
import os
import sys

print(f"🐍 Python: {sys.executable}")
try:
    import google.colab
    IN_COLAB = True
except ImportError:
    IN_COLAB = False

if IN_COLAB:
    if not os.path.exists("installed.marker"):
        import subprocess
        subprocess.run([sys.executable, "-m", "pip", "install", "-q",
                        "synth-pdb", "synth-nmr", "biotite",
                        "matplotlib", "ipywidgets", "numpy"])
        with open("installed.marker", "w") as f: f.write("done")
        print("✅ Installed. Please run this cell once more if imports fail.")
    else:
        print("✅ Already installed.")
else:
    print("✅ Running locally — no install needed.")

🧲 The RDC Alignment Tensor Explorer¶

Residual Dipolar Couplings (RDCs) are one of the most powerful — and most misunderstood — NMR observables. This notebook makes the key concepts visual and interactive.

🎯 What You'll Learn¶

Section	Concept
1	What is an alignment tensor? (interactive 3D polar map)
2	How Da and R change the RDC pattern for a real protein
3	What alignment media experimentalists actually use
4	Why RDCs and NOEs are complementary restraints

🔬 The Physics in One Sentence¶

In a weak alignment medium, a protein's N–H bond vectors are partially oriented relative to the magnetic field. The residual coupling you measure depends on the angle of each bond vector:

$$D(\theta, \phi) = D_a \left[ (3\cos^2\theta - 1) + \tfrac{3}{2}R\,\sin^2\theta\,\cos 2\phi \right]$$

Where:

$\theta$, $\phi$ = polar / azimuthal angle of the N–H bond in the alignment tensor frame
$D_a$ = axial component of the alignment tensor (Hz); controls magnitude
$R$ = rhombicity (0 ≤ R ≤ 2/3); controls asymmetry

Reference: Tjandra, N. & Bax, A. (1997). Science, 278, 1111–1114. DOI: 10.1126/science.278.5340.1111

In [ ]:

Copied!





# ============================================================
# Core imports and RDC formula
# ============================================================
import io

import ipywidgets as widgets
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display

plt.rcParams.update({
    'figure.dpi': 110,
    'font.size': 11,
    'axes.spines.top': False,
    'axes.spines.right': False,
})

def rdc_formula(theta_deg, phi_deg, Da, R):
    """Compute D(θ, φ) = Da * [(3cos²θ − 1) + (3/2)·R·sin²θ·cos(2φ)].

    Tjandra & Bax (1997), Science 278:1111.

    Parameters
    ----------
    theta_deg : float or ndarray  — polar angle in degrees (0° = Z-axis)
    phi_deg   : float or ndarray  — azimuthal angle in degrees
    Da        : float             — axial component (Hz)
    R         : float             — rhombicity (0 ≤ R ≤ 2/3)

    """
    theta = np.deg2rad(theta_deg)
    phi   = np.deg2rad(phi_deg)
    cos2  = np.cos(theta) ** 2
    sin2  = np.sin(theta) ** 2
    return Da * ((3 * cos2 - 1) + 1.5 * R * sin2 * np.cos(2 * phi))

print("✅ RDC formula loaded.")
print(f"\nSanity check (Z-axis, Da=10, R=0):  {rdc_formula(0, 0, 10, 0):.1f} Hz  (expected: 20.0)")
print(f"Sanity check (X-axis, Da=10, R=0):  {rdc_formula(90, 0, 10, 0):.1f} Hz  (expected: -10.0)")
# ============================================================
# Core imports and RDC formula
# ============================================================
import io

import ipywidgets as widgets
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display

plt.rcParams.update({
    'figure.dpi': 110,
    'font.size': 11,
    'axes.spines.top': False,
    'axes.spines.right': False,
})

def rdc_formula(theta_deg, phi_deg, Da, R):
    """Compute D(θ, φ) = Da * [(3cos²θ − 1) + (3/2)·R·sin²θ·cos(2φ)].

    Tjandra & Bax (1997), Science 278:1111.

    Parameters
    ----------
    theta_deg : float or ndarray  — polar angle in degrees (0° = Z-axis)
    phi_deg   : float or ndarray  — azimuthal angle in degrees
    Da        : float             — axial component (Hz)
    R         : float             — rhombicity (0 ≤ R ≤ 2/3)

    """
    theta = np.deg2rad(theta_deg)
    phi   = np.deg2rad(phi_deg)
    cos2  = np.cos(theta) ** 2
    sin2  = np.sin(theta) ** 2
    return Da * ((3 * cos2 - 1) + 1.5 * R * sin2 * np.cos(2 * phi))

print("✅ RDC formula loaded.")
print(f"\nSanity check (Z-axis, Da=10, R=0):  {rdc_formula(0, 0, 10, 0):.1f} Hz  (expected: 20.0)")
print(f"Sanity check (X-axis, Da=10, R=0):  {rdc_formula(90, 0, 10, 0):.1f} Hz  (expected: -10.0)")

🌍 Section 1: The Alignment Tensor as a 3D Color Map¶

The big idea: Every point on this sphere represents a possible N–H bond vector direction. The color of that point shows the RDC value you would measure if your N–H bond pointed that way.

🔴 Red / warm = large positive RDC (bond nearly parallel to the Z-axis)
🔵 Blue / cool = large negative RDC (bond nearly perpendicular to Z)
⚪ White = near-zero RDC (the "magic angle", θ ≈ 54.7°)

Key insight: $D_a$ stretches the color scale (bigger Da = hotter colors); $R$ breaks the circular symmetry around the equator (try R = 0 vs R = 0.5).

💡 Try it: Set R = 0 (axially symmetric). The sphere looks like a simple gradient from top to equator. Now increase R to 0.5 and watch the equator split into two distinct regions.

In [ ]:

Copied!





import os

IN_CI = bool(os.getenv("CI"))

if not IN_CI:
    # ============================================================
    # Section 1: Interactive 3D Polar Map of the Alignment Tensor
    # ============================================================

    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

    # Output widget so the plot updates cleanly in place
    polar_output = widgets.Output()

    # Build sliders
    da_slider_1 = widgets.FloatSlider(
        value=10.0, min=1.0, max=25.0, step=0.5,
        description='Dₐ (Hz)',
        style={'description_width': '80px'},
        layout=widgets.Layout(width='420px'),
        continuous_update=False,
        tooltip="Axial component of the alignment tensor. Controls magnitude of RDC values."
    )
    r_slider_1 = widgets.FloatSlider(
        value=0.1, min=0.0, max=0.665, step=0.025,
        description='R (rhombicity)',
        style={'description_width': '100px'},
        layout=widgets.Layout(width='420px'),
        continuous_update=False,
        tooltip="Rhombicity breaks axial symmetry. R=0 = axially symmetric; R=2/3 = maximum."
    )

    def draw_tensor_sphere(Da, R):
        """Render a colored unit sphere showing D(theta, phi) at each surface point."""
        # Clamp R defensively — floating-point slider steps can accumulate
        # to a value infinitesimally above 2/3, which the synth-nmr validator rejects.
        R = float(np.clip(R, 0.0, 2.0/3.0 - 1e-9))
        # Grid over the sphere
        theta_arr = np.linspace(0, 180, 120)     # polar angle (degrees)
        phi_arr   = np.linspace(0, 360, 240)     # azimuthal angle (degrees)
        THETA, PHI = np.meshgrid(theta_arr, phi_arr)

        # RDC at every (theta, phi)
        D = rdc_formula(THETA, PHI, Da, R)

        # Convert to Cartesian for plotting
        T = np.deg2rad(THETA)
        P = np.deg2rad(PHI)
        X = np.sin(T) * np.cos(P)
        Y = np.sin(T) * np.sin(P)
        Z = np.cos(T)

        # Symmetric color scale centred on zero
        limit = 2 * abs(Da)

        plt.ioff()
        fig = plt.figure(figsize=(9, 7))
        ax  = fig.add_subplot(111, projection='3d')

        # Color the surface by RDC value
        norm = mcolors.TwoSlopeNorm(vmin=-limit, vcenter=0, vmax=limit)
        fc   = plt.cm.RdBu_r(norm(D))

        ax.plot_surface(
            X, Y, Z,
            facecolors=fc,
            rstride=2, cstride=2,
            alpha=0.92, linewidth=0,
            antialiased=True
        )

        # Color bar
        mappable = plt.cm.ScalarMappable(cmap='RdBu_r', norm=norm)
        mappable.set_array(D)
        cb = fig.colorbar(mappable, ax=ax, shrink=0.5, pad=0.12)
        cb.set_label('D_{NH} (Hz)', fontsize=11)

        # Axis labels
        ax.set_xlabel('X', fontsize=10, labelpad=2)
        ax.set_ylabel('Y', fontsize=10, labelpad=2)
        ax.set_zlabel('Z (principal axis)', fontsize=10, labelpad=2)
        ax.set_title(
            f'Alignment Tensor: Dₐ = {Da:.1f} Hz,  R = {R:.3f}\n'
            f'Range: [{-limit:.1f}, {2*Da:.1f}] Hz  |  Magic angle: cos⁻¹(1/√3) ≈ 54.7°',
            fontsize=11, pad=10
        )

        # Draw the magic-angle circle (theta ≈ 54.7° → RDC = 0 when R=0)
        theta_magic = np.arccos(1.0 / np.sqrt(3.0))
        phi_ring    = np.linspace(0, 2 * np.pi, 300)
        ax.plot(
            np.sin(theta_magic) * np.cos(phi_ring),
            np.sin(theta_magic) * np.sin(phi_ring),
            np.full_like(phi_ring, np.cos(theta_magic)),
            'k--', lw=1.5, alpha=0.55, label='Magic angle θ ≈ 54.7° (D=0 when R=0)'
        )
        ax.legend(loc='upper left', fontsize=8, framealpha=0.5)

        ax.set_box_aspect([1, 1, 1])
        plt.tight_layout()

        buf = io.BytesIO()
        fig.savefig(buf, format='png', bbox_inches='tight')
        plt.close(fig)
        buf.seek(0)
        return buf.read()

    sphere_image = widgets.Image(format='png')

    def update_sphere(change=None):
        sphere_image.value = draw_tensor_sphere(da_slider_1.value, r_slider_1.value)

    da_slider_1.observe(update_sphere, names='value')
    r_slider_1.observe(update_sphere, names='value')

    update_sphere()   # Initial render

    display(
        widgets.VBox([
            widgets.HTML("<b>Alignment Tensor — 3D Surface Map</b><br>"
                         "<i>Drag the sliders to see how Dₐ and R reshape the tensor landscape.</i>"),
            da_slider_1,
            r_slider_1,
            sphere_image,
        ])
    )

# (Widget output skipped in CI)
import os

IN_CI = bool(os.getenv("CI"))

if not IN_CI:
    # ============================================================
    # Section 1: Interactive 3D Polar Map of the Alignment Tensor
    # ============================================================

    from mpl_toolkits.mplot3d import Axes3D  # noqa: F401

    # Output widget so the plot updates cleanly in place
    polar_output = widgets.Output()

    # Build sliders
    da_slider_1 = widgets.FloatSlider(
        value=10.0, min=1.0, max=25.0, step=0.5,
        description='Dₐ (Hz)',
        style={'description_width': '80px'},
        layout=widgets.Layout(width='420px'),
        continuous_update=False,
        tooltip="Axial component of the alignment tensor. Controls magnitude of RDC values."
    )
    r_slider_1 = widgets.FloatSlider(
        value=0.1, min=0.0, max=0.665, step=0.025,
        description='R (rhombicity)',
        style={'description_width': '100px'},
        layout=widgets.Layout(width='420px'),
        continuous_update=False,
        tooltip="Rhombicity breaks axial symmetry. R=0 = axially symmetric; R=2/3 = maximum."
    )

    def draw_tensor_sphere(Da, R):
        """Render a colored unit sphere showing D(theta, phi) at each surface point."""
        # Clamp R defensively — floating-point slider steps can accumulate
        # to a value infinitesimally above 2/3, which the synth-nmr validator rejects.
        R = float(np.clip(R, 0.0, 2.0/3.0 - 1e-9))
        # Grid over the sphere
        theta_arr = np.linspace(0, 180, 120)     # polar angle (degrees)
        phi_arr   = np.linspace(0, 360, 240)     # azimuthal angle (degrees)
        THETA, PHI = np.meshgrid(theta_arr, phi_arr)

        # RDC at every (theta, phi)
        D = rdc_formula(THETA, PHI, Da, R)

        # Convert to Cartesian for plotting
        T = np.deg2rad(THETA)
        P = np.deg2rad(PHI)
        X = np.sin(T) * np.cos(P)
        Y = np.sin(T) * np.sin(P)
        Z = np.cos(T)

        # Symmetric color scale centred on zero
        limit = 2 * abs(Da)

        plt.ioff()
        fig = plt.figure(figsize=(9, 7))
        ax  = fig.add_subplot(111, projection='3d')

        # Color the surface by RDC value
        norm = mcolors.TwoSlopeNorm(vmin=-limit, vcenter=0, vmax=limit)
        fc   = plt.cm.RdBu_r(norm(D))

        ax.plot_surface(
            X, Y, Z,
            facecolors=fc,
            rstride=2, cstride=2,
            alpha=0.92, linewidth=0,
            antialiased=True
        )

        # Color bar
        mappable = plt.cm.ScalarMappable(cmap='RdBu_r', norm=norm)
        mappable.set_array(D)
        cb = fig.colorbar(mappable, ax=ax, shrink=0.5, pad=0.12)
        cb.set_label('D_{NH} (Hz)', fontsize=11)

        # Axis labels
        ax.set_xlabel('X', fontsize=10, labelpad=2)
        ax.set_ylabel('Y', fontsize=10, labelpad=2)
        ax.set_zlabel('Z (principal axis)', fontsize=10, labelpad=2)
        ax.set_title(
            f'Alignment Tensor: Dₐ = {Da:.1f} Hz,  R = {R:.3f}\n'
            f'Range: [{-limit:.1f}, {2*Da:.1f}] Hz  |  Magic angle: cos⁻¹(1/√3) ≈ 54.7°',
            fontsize=11, pad=10
        )

        # Draw the magic-angle circle (theta ≈ 54.7° → RDC = 0 when R=0)
        theta_magic = np.arccos(1.0 / np.sqrt(3.0))
        phi_ring    = np.linspace(0, 2 * np.pi, 300)
        ax.plot(
            np.sin(theta_magic) * np.cos(phi_ring),
            np.sin(theta_magic) * np.sin(phi_ring),
            np.full_like(phi_ring, np.cos(theta_magic)),
            'k--', lw=1.5, alpha=0.55, label='Magic angle θ ≈ 54.7° (D=0 when R=0)'
        )
        ax.legend(loc='upper left', fontsize=8, framealpha=0.5)

        ax.set_box_aspect([1, 1, 1])
        plt.tight_layout()

        buf = io.BytesIO()
        fig.savefig(buf, format='png', bbox_inches='tight')
        plt.close(fig)
        buf.seek(0)
        return buf.read()

    sphere_image = widgets.Image(format='png')

    def update_sphere(change=None):
        sphere_image.value = draw_tensor_sphere(da_slider_1.value, r_slider_1.value)

    da_slider_1.observe(update_sphere, names='value')
    r_slider_1.observe(update_sphere, names='value')

    update_sphere()   # Initial render

    display(
        widgets.VBox([
            widgets.HTML("Alignment Tensor — 3D Surface Map
"
                         "Drag the sliders to see how Dₐ and R reshape the tensor landscape."),
            da_slider_1,
            r_slider_1,
            sphere_image,
        ])
    )

# (Widget output skipped in CI)

💡 Key Observations from the Sphere¶

Configuration	What You See	Why
R = 0	Perfectly rotationally symmetric around the Z-axis	R = 0 → cos(2φ) term vanishes → D depends only on θ
R = 0.33	Equator shows mild east–west asymmetry	The cos(2φ) term breaks ϕ-symmetry
R = 0.67	Strong four-lobed pattern at equator	Maximum rhombicity — tensor is squashed into a brick shape
Large Dₐ	Stronger colors (larger range)	Dₐ simply multiplies the whole function

The dashed ring marks the magic angle (θ ≈ 54.7°), where $3\cos^2\theta - 1 = 0$. For an axially symmetric tensor (R=0), bond vectors at this angle give zero RDC — they can't distinguish tensor orientation at all!

🧬 Section 2: Live RDC Profile for a Synthetic Protein¶

Now connect the abstract tensor to a real protein. Below we generate a synthetic peptide and compute its backbone N–H RDC values using synth_pdb.rdc.calculate_rdcs.

Each bar is one residue. The color tells you the sign:

🔴 Red = positive RDC (N–H bond closer to the tensor Z-axis)
🔵 Blue = negative RDC (N–H bond perpendicular to Z-axis)

Dashed horizontal lines show the theoretical limits — no measured value can lie outside these bounds.

Structural implication: Residues in a regular α-helix all have similar N–H orientations → consecutive residues show similar, smoothly varying RDC values. Loops and disordered regions show more irregular patterns.

In [ ]:

Copied!





# ============================================================
# Section 2: Generate structure + live RDC bar chart
# ============================================================
import biotite.structure.io.pdb as pdbio

from synth_pdb.generator import generate_pdb_content
from synth_pdb.rdc import calculate_rdcs

# Use a sequence with clear secondary structure mix so RDC variation is visible
SEQUENCE = "VKITVGGTLTVALGGALALALALA"
STRUCT_DEF = "1-5:beta,6-8:random,9-13:beta,14-16:random,17-24:alpha"

print("🧬 Generating synthetic protein (beta + loop + beta + loop + helix)...")
print("   This may take ~30 seconds without energy minimization.")

pdb_content = generate_pdb_content(
    sequence_str=SEQUENCE,
    structure=STRUCT_DEF,
    optimize_sidechains=True,
    minimize_energy=False,   # Fast for tutorial; flip to True for better H positions
)

pdb_file  = pdbio.PDBFile.read(io.StringIO(pdb_content))
structure = pdb_file.get_structure(model=1)

print(f"✅ Structure ready — {len(structure)} atoms, sequence: {SEQUENCE}")
print("\nNote: RDC calculation requires backbone N and H atoms.")
print("The generator adds idealised H positions even without --minimize.")
# ============================================================
# Section 2: Generate structure + live RDC bar chart
# ============================================================
import biotite.structure.io.pdb as pdbio

from synth_pdb.generator import generate_pdb_content
from synth_pdb.rdc import calculate_rdcs

# Use a sequence with clear secondary structure mix so RDC variation is visible
SEQUENCE = "VKITVGGTLTVALGGALALALALA"
STRUCT_DEF = "1-5:beta,6-8:random,9-13:beta,14-16:random,17-24:alpha"

print("🧬 Generating synthetic protein (beta + loop + beta + loop + helix)...")
print("   This may take ~30 seconds without energy minimization.")

pdb_content = generate_pdb_content(
    sequence_str=SEQUENCE,
    structure=STRUCT_DEF,
    optimize_sidechains=True,
    minimize_energy=False,   # Fast for tutorial; flip to True for better H positions
)

pdb_file  = pdbio.PDBFile.read(io.StringIO(pdb_content))
structure = pdb_file.get_structure(model=1)

print(f"✅ Structure ready — {len(structure)} atoms, sequence: {SEQUENCE}")
print("\nNote: RDC calculation requires backbone N and H atoms.")
print("The generator adds idealised H positions even without --minimize.")

In [ ]:

Copied!





import os

IN_CI = bool(os.getenv("CI"))

# Secondary-structure color lookup for region shading (defined here so later cells
# can use SS_REGIONS even when the interactive widget is skipped in CI)
SS_REGIONS = [
    (1,  5,  'β-strand 1', '#2980b9', 0.12),
    (6,  8,  'Loop',       '#7f8c8d', 0.06),
    (9,  13, 'β-strand 2', '#2980b9', 0.12),
    (14, 16, 'Loop',       '#7f8c8d', 0.06),
    (17, 24, 'α-helix',    '#e74c3c', 0.12),
]

if not IN_CI:

    # Create persistent image widget
    rdc_bar_image = widgets.Image(format='png')

    da_slider_2 = widgets.FloatSlider(
        value=10.0, min=1.0, max=25.0, step=0.5,
        description='Dₐ (Hz)',
        style={'description_width': '80px'},
        layout=widgets.Layout(width='450px'),
        continuous_update=False,
    )
    r_slider_2 = widgets.FloatSlider(
        value=0.1, min=0.0, max=0.665, step=0.025,
        description='R (rhombicity)',
        style={'description_width': '100px'},
        layout=widgets.Layout(width='450px'),
        continuous_update=False,
    )

    def draw_rdc_bar(Da, R):
        # Clamp R defensively — floating-point slider steps can accumulate
        # to a value infinitesimally above 2/3, which the synth-nmr validator rejects.
        R = float(np.clip(R, 0.0, 2.0/3.0 - 1e-9))
        rdcs = calculate_rdcs(structure, Da=Da, R=R)
        if not rdcs:
            print("⚠️ No RDCs computed — structure may lack backbone H atoms.")
            print("   Try re-running with minimize_energy=True above.")
            return None

        residues  = sorted(rdcs.keys())
        rdc_vals  = [rdcs[r] for r in residues]
        bar_colors = ['#e74c3c' if v >= 0 else '#2980b9' for v in rdc_vals]

        # Physical limits for this tensor
        rdc_max  =  2.0 * Da
        rdc_min  =  Da * (-1.0 - 1.5 * R)

        plt.ioff()
        fig, ax = plt.subplots(figsize=(12, 4.5))

        # Secondary structure shading
        for r_start, r_end, label, color, alpha in SS_REGIONS:
            ax.axvspan(r_start - 0.5, r_end + 0.5, color=color, alpha=alpha,
                       label=label if r_start in (1, 14, 17) else "")

        ax.bar(residues, rdc_vals, color=bar_colors, alpha=0.85,
               edgecolor='white', linewidth=0.5, zorder=3)

        # Physical limit lines
        ax.axhline(rdc_max, color='#c0392b', ls='--', lw=1.4, alpha=0.7,
                   label=f'Max: +2·Dₐ = {rdc_max:.1f} Hz  (θ=0°)')
        ax.axhline(rdc_min, color='#1a5276', ls='--', lw=1.4, alpha=0.7,
                   label=f'Min: {rdc_min:.1f} Hz  (θ=90°, φ=90°)')
        ax.axhline(0, color='black', lw=0.8, alpha=0.4)

        ax.set_xlabel('Residue Number', fontsize=12)
        ax.set_ylabel('¹D_{NH} (Hz)', fontsize=12)
        ax.set_title(
            f'Synthetic Backbone N–H RDCs   |   Dₐ = {Da:.1f} Hz,  R = {R:.3f}\n'
            f'Red = positive (N–H near Z-axis)   Blue = negative (N–H near XY-plane)',
            fontsize=11
        )
        ax.set_xlim(0.5, max(residues) + 0.5)
        ax.set_ylim(rdc_min * 1.15, rdc_max * 1.15)
        ax.grid(axis='y', ls='--', alpha=0.35, zorder=0)
        ax.legend(fontsize=9, loc='upper right', framealpha=0.8)

        plt.tight_layout()
        buf = io.BytesIO()
        fig.savefig(buf, format='png', bbox_inches='tight')
        plt.close(fig)
        buf.seek(0)
        return buf.read()

    def update_rdc_bar(change=None):
        data = draw_rdc_bar(da_slider_2.value, r_slider_2.value)
        if data: rdc_bar_image.value = data

    da_slider_2.observe(update_rdc_bar, names='value')
    r_slider_2.observe(update_rdc_bar, names='value')

    update_rdc_bar()  # Initial render

    display(
        widgets.VBox([
            widgets.HTML(
                "<b>Live RDC Profile — drag sliders to explore</b><br>"
                "<i>Shaded regions: 🔵 β-strands, 🔴 α-helix, ⬜ loops</i>"
            ),
            da_slider_2,
            r_slider_2,
            rdc_bar_image,
        ])
    )

# (Widget output skipped in CI)
import os

IN_CI = bool(os.getenv("CI"))

# Secondary-structure color lookup for region shading (defined here so later cells
# can use SS_REGIONS even when the interactive widget is skipped in CI)
SS_REGIONS = [
    (1,  5,  'β-strand 1', '#2980b9', 0.12),
    (6,  8,  'Loop',       '#7f8c8d', 0.06),
    (9,  13, 'β-strand 2', '#2980b9', 0.12),
    (14, 16, 'Loop',       '#7f8c8d', 0.06),
    (17, 24, 'α-helix',    '#e74c3c', 0.12),
]

if not IN_CI:

    # Create persistent image widget
    rdc_bar_image = widgets.Image(format='png')

    da_slider_2 = widgets.FloatSlider(
        value=10.0, min=1.0, max=25.0, step=0.5,
        description='Dₐ (Hz)',
        style={'description_width': '80px'},
        layout=widgets.Layout(width='450px'),
        continuous_update=False,
    )
    r_slider_2 = widgets.FloatSlider(
        value=0.1, min=0.0, max=0.665, step=0.025,
        description='R (rhombicity)',
        style={'description_width': '100px'},
        layout=widgets.Layout(width='450px'),
        continuous_update=False,
    )

    def draw_rdc_bar(Da, R):
        # Clamp R defensively — floating-point slider steps can accumulate
        # to a value infinitesimally above 2/3, which the synth-nmr validator rejects.
        R = float(np.clip(R, 0.0, 2.0/3.0 - 1e-9))
        rdcs = calculate_rdcs(structure, Da=Da, R=R)
        if not rdcs:
            print("⚠️ No RDCs computed — structure may lack backbone H atoms.")
            print("   Try re-running with minimize_energy=True above.")
            return None

        residues  = sorted(rdcs.keys())
        rdc_vals  = [rdcs[r] for r in residues]
        bar_colors = ['#e74c3c' if v >= 0 else '#2980b9' for v in rdc_vals]

        # Physical limits for this tensor
        rdc_max  =  2.0 * Da
        rdc_min  =  Da * (-1.0 - 1.5 * R)

        plt.ioff()
        fig, ax = plt.subplots(figsize=(12, 4.5))

        # Secondary structure shading
        for r_start, r_end, label, color, alpha in SS_REGIONS:
            ax.axvspan(r_start - 0.5, r_end + 0.5, color=color, alpha=alpha,
                       label=label if r_start in (1, 14, 17) else "")

        ax.bar(residues, rdc_vals, color=bar_colors, alpha=0.85,
               edgecolor='white', linewidth=0.5, zorder=3)

        # Physical limit lines
        ax.axhline(rdc_max, color='#c0392b', ls='--', lw=1.4, alpha=0.7,
                   label=f'Max: +2·Dₐ = {rdc_max:.1f} Hz  (θ=0°)')
        ax.axhline(rdc_min, color='#1a5276', ls='--', lw=1.4, alpha=0.7,
                   label=f'Min: {rdc_min:.1f} Hz  (θ=90°, φ=90°)')
        ax.axhline(0, color='black', lw=0.8, alpha=0.4)

        ax.set_xlabel('Residue Number', fontsize=12)
        ax.set_ylabel('¹D_{NH} (Hz)', fontsize=12)
        ax.set_title(
            f'Synthetic Backbone N–H RDCs   |   Dₐ = {Da:.1f} Hz,  R = {R:.3f}\n'
            f'Red = positive (N–H near Z-axis)   Blue = negative (N–H near XY-plane)',
            fontsize=11
        )
        ax.set_xlim(0.5, max(residues) + 0.5)
        ax.set_ylim(rdc_min * 1.15, rdc_max * 1.15)
        ax.grid(axis='y', ls='--', alpha=0.35, zorder=0)
        ax.legend(fontsize=9, loc='upper right', framealpha=0.8)

        plt.tight_layout()
        buf = io.BytesIO()
        fig.savefig(buf, format='png', bbox_inches='tight')
        plt.close(fig)
        buf.seek(0)
        return buf.read()

    def update_rdc_bar(change=None):
        data = draw_rdc_bar(da_slider_2.value, r_slider_2.value)
        if data: rdc_bar_image.value = data

    da_slider_2.observe(update_rdc_bar, names='value')
    r_slider_2.observe(update_rdc_bar, names='value')

    update_rdc_bar()  # Initial render

    display(
        widgets.VBox([
            widgets.HTML(
                "Live RDC Profile — drag sliders to explore
"
                "Shaded regions: 🔵 β-strands, 🔴 α-helix, ⬜ loops"
            ),
            da_slider_2,
            r_slider_2,
            rdc_bar_image,
        ])
    )

# (Widget output skipped in CI)

🔍 Experiments to Try¶

Helix vs strand: Which secondary structure shows more uniform RDC values? Can you see the periodic variation in the β-strand regions?
Dₐ only → uniform scaling: Set R=0 and change Dₐ. All bars scale proportionally — same pattern, different magnitude.
R reveals asymmetry: Now set Dₐ=10, R=0.667 (maximum). Some residues that were mildly positive become strongly negative. The rhombic term cos(2φ) discriminates between X and Y orientations — information you simply can't get from an axially symmetric tensor.
The prolines: Notice Proline residues are absent from the plot. Proline has no backbone amide H (its nitrogen is a tertiary amine in a pyrrolidine ring), so it contributes no ¹D_{NH} RDC.

🧪 Section 3: Alignment Media Comparison¶

Different experimental alignment media give different Da and R values. This panel compares five common media used in NMR labs:

Medium	Dₐ (Hz)	R	Reference
DMPC/DHPC bicelles	8–15	0.0–0.1	Tjandra & Bax (1997)
Pf1 filamentous phage	10–20	0.05–0.15	Hansen et al. (2000)
Stretched polyacrylamide	4–12	0.1–0.3	Tycko et al. (2000)
Cellulose crystallites	5–10	0.0–0.05	Bax & Grishaev (2005)
Purple membrane	15–25	0.2–0.4	Sanders & Schwonek (1992)

The next cell shows the same protein's RDC profile under each alignment medium (using representative Da and R values). This illustrates why experimentalists sometimes measure in multiple media — each one provides a different view of the molecular orientation.

In [ ]:

Copied!





# ============================================================
# Section 3: Alignment Media Comparison
# ============================================================

MEDIA = [
    # (label, Da, R, color)
    # Representative central values from the table above
    ('DMPC/DHPC\nbicelles',     11.0, 0.05, '#1abc9c'),
    ('Pf1 phage',               15.0, 0.10, '#9b59b6'),
    ('Stretched\npolyacrylamide', 8.0, 0.20, '#e67e22'),
    ('Cellulose\ncrystallites',   7.0, 0.02, '#2ecc71'),
    ('Purple\nmembrane',         20.0, 0.30, '#e74c3c'),
]

fig, axes = plt.subplots(len(MEDIA), 1, figsize=(12, 2.8 * len(MEDIA)), sharex=True)

for ax, (label, Da, R, color) in zip(axes, MEDIA):
    rdcs = calculate_rdcs(structure, Da=Da, R=R)
    if not rdcs:
        ax.text(0.5, 0.5, 'No RDCs (add H atoms with minimize_energy=True)',
                ha='center', va='center', transform=ax.transAxes)
        continue
    residues = sorted(rdcs.keys())
    vals     = [rdcs[r] for r in residues]
    rdc_max  =  2.0 * Da
    rdc_min  =  Da * (-1.0 - 1.5 * R)

    # SS shading
    for r_start, r_end, _, ss_col, ss_alpha in SS_REGIONS:
        ax.axvspan(r_start - 0.5, r_end + 0.5, color=ss_col, alpha=ss_alpha)

    bar_col = [color if v >= 0 else '#7f8c8d' for v in vals]
    ax.bar(residues, vals, color=bar_col, alpha=0.80,
           edgecolor='white', linewidth=0.4, zorder=3)
    ax.axhline(rdc_max, color='#c0392b', ls='--', lw=1.0, alpha=0.6)
    ax.axhline(rdc_min, color='#1a5276', ls='--', lw=1.0, alpha=0.6)
    ax.axhline(0, color='black', lw=0.6, alpha=0.35)
    ax.set_ylabel('D_{NH} (Hz)', fontsize=9)
    ax.set_title(f'{label}   (Dₐ={Da} Hz, R={R})', fontsize=10, loc='left', pad=4)
    ax.grid(axis='y', ls='--', alpha=0.3)
    ax.set_xlim(0.5, max(residues) + 0.5)

axes[-1].set_xlabel('Residue Number', fontsize=12)
fig.suptitle(
    'RDC Profile of the Same Protein in Five Different Alignment Media\n'
    'Each medium provides a different "view" of the molecular orientation',
    fontsize=12, fontweight='bold', y=1.01
)
plt.tight_layout()
plt.show()

print("\n📚 Why use multiple media?")
print("Each medium orients the protein differently (different Da, R, and Euler angles).")
print("Combining data from two or more media over-determines the alignment tensor,")
print("dramatically improving the precision of orientation constraints.")
print("Reference: Bax & Grishaev (2005), Curr Opin Struct Biol 15:563.")
# ============================================================
# Section 3: Alignment Media Comparison
# ============================================================

MEDIA = [
    # (label, Da, R, color)
    # Representative central values from the table above
    ('DMPC/DHPC\nbicelles',     11.0, 0.05, '#1abc9c'),
    ('Pf1 phage',               15.0, 0.10, '#9b59b6'),
    ('Stretched\npolyacrylamide', 8.0, 0.20, '#e67e22'),
    ('Cellulose\ncrystallites',   7.0, 0.02, '#2ecc71'),
    ('Purple\nmembrane',         20.0, 0.30, '#e74c3c'),
]

fig, axes = plt.subplots(len(MEDIA), 1, figsize=(12, 2.8 * len(MEDIA)), sharex=True)

for ax, (label, Da, R, color) in zip(axes, MEDIA):
    rdcs = calculate_rdcs(structure, Da=Da, R=R)
    if not rdcs:
        ax.text(0.5, 0.5, 'No RDCs (add H atoms with minimize_energy=True)',
                ha='center', va='center', transform=ax.transAxes)
        continue
    residues = sorted(rdcs.keys())
    vals     = [rdcs[r] for r in residues]
    rdc_max  =  2.0 * Da
    rdc_min  =  Da * (-1.0 - 1.5 * R)

    # SS shading
    for r_start, r_end, _, ss_col, ss_alpha in SS_REGIONS:
        ax.axvspan(r_start - 0.5, r_end + 0.5, color=ss_col, alpha=ss_alpha)

    bar_col = [color if v >= 0 else '#7f8c8d' for v in vals]
    ax.bar(residues, vals, color=bar_col, alpha=0.80,
           edgecolor='white', linewidth=0.4, zorder=3)
    ax.axhline(rdc_max, color='#c0392b', ls='--', lw=1.0, alpha=0.6)
    ax.axhline(rdc_min, color='#1a5276', ls='--', lw=1.0, alpha=0.6)
    ax.axhline(0, color='black', lw=0.6, alpha=0.35)
    ax.set_ylabel('D_{NH} (Hz)', fontsize=9)
    ax.set_title(f'{label}   (Dₐ={Da} Hz, R={R})', fontsize=10, loc='left', pad=4)
    ax.grid(axis='y', ls='--', alpha=0.3)
    ax.set_xlim(0.5, max(residues) + 0.5)

axes[-1].set_xlabel('Residue Number', fontsize=12)
fig.suptitle(
    'RDC Profile of the Same Protein in Five Different Alignment Media\n'
    'Each medium provides a different "view" of the molecular orientation',
    fontsize=12, fontweight='bold', y=1.01
)
plt.tight_layout()
plt.show()

print("\n📚 Why use multiple media?")
print("Each medium orients the protein differently (different Da, R, and Euler angles).")
print("Combining data from two or more media over-determines the alignment tensor,")
print("dramatically improving the precision of orientation constraints.")
print("Reference: Bax & Grishaev (2005), Curr Opin Struct Biol 15:563.")

🗺️ Section 4: RDCs vs NOEs — Local Distance vs Global Orientation¶

This is the most important conceptual diagram in RDC science.

NOEs (Nuclear Overhauser Effects):

Through space, signal ∝ r⁻⁶
Measure distances between nearby atoms
Sensitive only to contacts < 5–6 Å → local restraints
Build up secondary structure, but secondary elements can be mis-docked

RDCs (Residual Dipolar Couplings):

Measure bond vector orientations vs the alignment tensor
No distance cut-off — every N–H contributes
Constrain global fold topology — how secondary structure elements dock together
Completely insensitive to local packing but highly sensitive to long-range orientation

"RDCs are like a GPS giving global coordinates; NOEs are like bumping into walls to map a room." — paraphrased from Prestegard et al. (2000)

In [ ]:

Copied!





# ============================================================
# Section 4: RDC vs NOE — Direct Complementarity Visualisation
# ============================================================
from synth_pdb.contact import compute_contact_map
from synth_pdb.nmr import calculate_synthetic_noes

# Compute RDCs (current slider values) and NOEs
Da_final, R_final = 10.0, 0.1
rdcs_final = calculate_rdcs(structure, Da=Da_final, R=R_final)
noes       = calculate_synthetic_noes(structure, cutoff=5.0)
ca_dists   = compute_contact_map(structure, method='ca', power=None)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# ── Left: NOE Contact Map (distance matrix) ──────────────────────────────
ax0 = axes[0]
im0 = ax0.imshow(ca_dists, cmap='viridis_r', interpolation='nearest',
                 vmin=0, vmax=30)
plt.colorbar(im0, ax=ax0, label='Cα–Cα distance (Å)', shrink=0.80)
ax0.set_title('NOE / Contact Map\n(local: Cα–Cα distances)', fontsize=11)
ax0.set_xlabel('Residue Index'); ax0.set_ylabel('Residue Index')

# Overlay NOE pairs as scatter points
noe_i = [n['index_1'] - 1 for n in noes if n['index_1'] != n['index_2']]
noe_j = [n['index_2'] - 1 for n in noes if n['index_1'] != n['index_2']]
ax0.scatter(noe_j, noe_i, c='white', s=2, alpha=0.4, label=f'{len(noes)} NOE pairs')
ax0.legend(fontsize=8, loc='upper right')

# ── Middle: RDC Bar Chart (global information) ────────────────────────────
ax1 = axes[1]
if rdcs_final:
    r_ids  = sorted(rdcs_final.keys())
    r_vals = [rdcs_final[r] for r in r_ids]
    c_vals = ['#e74c3c' if v >= 0 else '#2980b9' for v in r_vals]
    ax1.barh(r_ids, r_vals, color=c_vals, alpha=0.85,
             edgecolor='white', linewidth=0.4)
    ax1.axvline(0, color='black', lw=0.8)
    ax1.axvline( 2*Da_final, color='#c0392b', ls='--', lw=1.2,
                label=f'+2Dₐ = {2*Da_final:.0f} Hz')
    ax1.axvline(-Da_final,   color='#1a5276', ls='--', lw=1.2,
                label=f'−Dₐ = {-Da_final:.0f} Hz')
    ax1.set_xlabel('¹D_{NH} (Hz)', fontsize=11)
    ax1.set_ylabel('Residue Number', fontsize=11)
    ax1.set_title(f'RDC Per Residue\n(global: bond orientations, Dₐ={Da_final} Hz)', fontsize=11)
    ax1.legend(fontsize=8)
    ax1.grid(axis='x', ls='--', alpha=0.35)
    # Shade SS regions on the y-axis
    for r_start, r_end, _, ss_col, ss_alpha in SS_REGIONS:
        ax1.axhspan(r_start - 0.5, r_end + 0.5, color=ss_col, alpha=ss_alpha)
else:
    ax1.text(0.5, 0.5, 'No RDCs computed.\nRerun with minimize_energy=True.',
             ha='center', va='center', transform=ax1.transAxes)

# ── Right: Schematic comparing information content ────────────────────────
ax2 = axes[2]
ax2.set_xlim(0, 10); ax2.set_ylim(0, 10); ax2.axis('off')

def box(ax, x, y, w, h, color, title, lines):
    from matplotlib.patches import FancyBboxPatch
    rect = FancyBboxPatch((x, y), w, h, boxstyle="round,pad=0.15",
                           linewidth=1.5, edgecolor='black',
                           facecolor=color, alpha=0.18)
    ax.add_patch(rect)
    ax.text(x + w/2, y + h - 0.35, title, ha='center', va='top',
            fontsize=10, fontweight='bold')
    for i, line in enumerate(lines):
        ax.text(x + 0.2, y + h - 0.85 - i * 0.6, f'• {line}',
                fontsize=8.5, va='top')

box(ax2, 0.3, 5.5, 4.4, 4.0, '#e74c3c', '📍 NOEs (Local)',
    ['Distance < 5–6 Å', 'Through space (r⁻⁶)', 'Secondary structure', 'Short-range only', 'Many restraints'])
box(ax2, 5.3, 5.5, 4.4, 4.0, '#2980b9', '🧭 RDCs (Global)',
    ['All residues', 'Bond orientation vs tensor', 'Tertiary fold topology', 'Long-range by design', 'One per N–H'])

# Arrow showing complementarity
ax2.annotate('', xy=(5.1, 7.5), xytext=(4.9, 7.5),
             arrowprops={"arrowstyle": '<->', "color": 'black', "lw": 2.0})
ax2.text(5.0, 8.1, 'Complementary!', ha='center', fontsize=9,
         fontweight='bold', color='#27ae60')

ax2.text(5.0, 4.8, 'Combined = uniquely determined fold', ha='center',
         fontsize=9, fontstyle='italic', color='#555')
ax2.text(5.0, 1.4,
     'Bewley & Clore (2000):\n'
     'RDC constraints refined HIV-1\n'
     'protease dimer to 0.4 Å RMSD\n'
     'without any additional NOEs.',
     ha='center', fontsize=8, color='#333',
     bbox={"boxstyle": 'round,pad=0.4', "fc": '#f9f9f9', "ec": '#bbb', "lw": 1})
ax2.set_title('Information Content\nNOE vs RDC', fontsize=11)

fig.suptitle(
    'NOEs probe LOCAL structure (contacts < 5 Å)\n'
    'RDCs probe GLOBAL orientation (all bond vectors vs alignment tensor)',
    fontsize=12, fontweight='bold'
)
plt.tight_layout()
plt.show()
# ============================================================
# Section 4: RDC vs NOE — Direct Complementarity Visualisation
# ============================================================
from synth_pdb.contact import compute_contact_map
from synth_pdb.nmr import calculate_synthetic_noes

# Compute RDCs (current slider values) and NOEs
Da_final, R_final = 10.0, 0.1
rdcs_final = calculate_rdcs(structure, Da=Da_final, R=R_final)
noes       = calculate_synthetic_noes(structure, cutoff=5.0)
ca_dists   = compute_contact_map(structure, method='ca', power=None)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# ── Left: NOE Contact Map (distance matrix) ──────────────────────────────
ax0 = axes[0]
im0 = ax0.imshow(ca_dists, cmap='viridis_r', interpolation='nearest',
                 vmin=0, vmax=30)
plt.colorbar(im0, ax=ax0, label='Cα–Cα distance (Å)', shrink=0.80)
ax0.set_title('NOE / Contact Map\n(local: Cα–Cα distances)', fontsize=11)
ax0.set_xlabel('Residue Index'); ax0.set_ylabel('Residue Index')

# Overlay NOE pairs as scatter points
noe_i = [n['index_1'] - 1 for n in noes if n['index_1'] != n['index_2']]
noe_j = [n['index_2'] - 1 for n in noes if n['index_1'] != n['index_2']]
ax0.scatter(noe_j, noe_i, c='white', s=2, alpha=0.4, label=f'{len(noes)} NOE pairs')
ax0.legend(fontsize=8, loc='upper right')

# ── Middle: RDC Bar Chart (global information) ────────────────────────────
ax1 = axes[1]
if rdcs_final:
    r_ids  = sorted(rdcs_final.keys())
    r_vals = [rdcs_final[r] for r in r_ids]
    c_vals = ['#e74c3c' if v >= 0 else '#2980b9' for v in r_vals]
    ax1.barh(r_ids, r_vals, color=c_vals, alpha=0.85,
             edgecolor='white', linewidth=0.4)
    ax1.axvline(0, color='black', lw=0.8)
    ax1.axvline( 2*Da_final, color='#c0392b', ls='--', lw=1.2,
                label=f'+2Dₐ = {2*Da_final:.0f} Hz')
    ax1.axvline(-Da_final,   color='#1a5276', ls='--', lw=1.2,
                label=f'−Dₐ = {-Da_final:.0f} Hz')
    ax1.set_xlabel('¹D_{NH} (Hz)', fontsize=11)
    ax1.set_ylabel('Residue Number', fontsize=11)
    ax1.set_title(f'RDC Per Residue\n(global: bond orientations, Dₐ={Da_final} Hz)', fontsize=11)
    ax1.legend(fontsize=8)
    ax1.grid(axis='x', ls='--', alpha=0.35)
    # Shade SS regions on the y-axis
    for r_start, r_end, _, ss_col, ss_alpha in SS_REGIONS:
        ax1.axhspan(r_start - 0.5, r_end + 0.5, color=ss_col, alpha=ss_alpha)
else:
    ax1.text(0.5, 0.5, 'No RDCs computed.\nRerun with minimize_energy=True.',
             ha='center', va='center', transform=ax1.transAxes)

# ── Right: Schematic comparing information content ────────────────────────
ax2 = axes[2]
ax2.set_xlim(0, 10); ax2.set_ylim(0, 10); ax2.axis('off')

def box(ax, x, y, w, h, color, title, lines):
    from matplotlib.patches import FancyBboxPatch
    rect = FancyBboxPatch((x, y), w, h, boxstyle="round,pad=0.15",
                           linewidth=1.5, edgecolor='black',
                           facecolor=color, alpha=0.18)
    ax.add_patch(rect)
    ax.text(x + w/2, y + h - 0.35, title, ha='center', va='top',
            fontsize=10, fontweight='bold')
    for i, line in enumerate(lines):
        ax.text(x + 0.2, y + h - 0.85 - i * 0.6, f'• {line}',
                fontsize=8.5, va='top')

box(ax2, 0.3, 5.5, 4.4, 4.0, '#e74c3c', '📍 NOEs (Local)',
    ['Distance < 5–6 Å', 'Through space (r⁻⁶)', 'Secondary structure', 'Short-range only', 'Many restraints'])
box(ax2, 5.3, 5.5, 4.4, 4.0, '#2980b9', '🧭 RDCs (Global)',
    ['All residues', 'Bond orientation vs tensor', 'Tertiary fold topology', 'Long-range by design', 'One per N–H'])

# Arrow showing complementarity
ax2.annotate('', xy=(5.1, 7.5), xytext=(4.9, 7.5),
             arrowprops={"arrowstyle": '<->', "color": 'black', "lw": 2.0})
ax2.text(5.0, 8.1, 'Complementary!', ha='center', fontsize=9,
         fontweight='bold', color='#27ae60')

ax2.text(5.0, 4.8, 'Combined = uniquely determined fold', ha='center',
         fontsize=9, fontstyle='italic', color='#555')
ax2.text(5.0, 1.4,
     'Bewley & Clore (2000):\n'
     'RDC constraints refined HIV-1\n'
     'protease dimer to 0.4 Å RMSD\n'
     'without any additional NOEs.',
     ha='center', fontsize=8, color='#333',
     bbox={"boxstyle": 'round,pad=0.4', "fc": '#f9f9f9', "ec": '#bbb', "lw": 1})
ax2.set_title('Information Content\nNOE vs RDC', fontsize=11)

fig.suptitle(
    'NOEs probe LOCAL structure (contacts < 5 Å)\n'
    'RDCs probe GLOBAL orientation (all bond vectors vs alignment tensor)',
    fontsize=12, fontweight='bold'
)
plt.tight_layout()
plt.show()

📈 Section 5: The RDC Distribution — Reading the Histogram¶

The histogram of RDC values for a folded protein has a characteristic shape that reveals tensor quality and structure information:

Bimodal distribution (two humps): typical of a folded protein with many residues in regular secondary structure — suggests a well-defined, non-random orientation distribution.
Single Gaussian: may indicate a partly disordered or random-coil protein.
Wide spread reaching the limits: suggests a large, well-defined alignment (good signal).

This cell also shows the sine-squared weighting that governs the distribution: since there are more bond vector orientations near the equator (θ ≈ 90°) than near the poles (θ ≈ 0°, 180°), and equatorial orientations give negative RDCs, perfectly isotropic proteins would be skewed towards negative values.

In [ ]:

Copied!





# ============================================================
# Section 5: RDC Distribution Histogram
# ============================================================

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# Compute for a few different Da values
da_vals = [5.0, 10.0, 20.0]
hist_colors = ['#3498db', '#e74c3c', '#2ecc71']

ax = axes[0]
for Da_h, col in zip(da_vals, hist_colors):
    rdcs_h = calculate_rdcs(structure, Da=Da_h, R=0.1)
    if rdcs_h:
        vals_h = list(rdcs_h.values())
        ax.hist(vals_h, bins=14, color=col, alpha=0.65,
                edgecolor='white', linewidth=0.5,
                label=f'Dₐ = {Da_h:.0f} Hz  (range: [{min(vals_h):.1f}, {max(vals_h):.1f}] Hz)')

ax.axvline(0, color='black', lw=0.8, alpha=0.5)
ax.set_xlabel('¹D_{NH} (Hz)', fontsize=12)
ax.set_ylabel('Residue Count', fontsize=12)
ax.set_title('RDC Distribution for Three Dₐ Values\n(R = 0.10, same structure)', fontsize=11)
ax.legend(fontsize=9)
ax.grid(axis='y', ls='--', alpha=0.35)

# ── Right: The sine² angular weighting ────────────────────────────────────
# Explains why purely random orientations are skewed towards negative RDC
ax2 = axes[1]
theta_lin    = np.linspace(0, np.pi, 500)
sinusoid_wgt = np.sin(theta_lin)  # Solid angle element: dΩ ∝ sin(θ)dθ
rdc_lin      = rdc_formula(np.rad2deg(theta_lin), 0, 10.0, 0.0)

# Weighted histogram of RDC values for isotropic orientation distribution
# Each theta contributes weight sin(theta) → probability element on sphere
n_samples = 50000
# Sample theta uniformly in cos(theta) — correct for isotropic distribution
cos_theta  = np.random.uniform(-1, 1, n_samples)
phi_rand   = np.random.uniform(0, 360, n_samples)
rdc_random = rdc_formula(np.rad2deg(np.arccos(cos_theta)), phi_rand, 10.0, 0.0)

ax2.hist(rdc_random, bins=60, color='#9b59b6', alpha=0.7,
         edgecolor='white', linewidth=0.3, density=True)
ax2.axvline(0, color='black', lw=0.8, alpha=0.5)
ax2.set_xlabel('¹D_{NH} (Hz)', fontsize=12)
ax2.set_ylabel('Probability Density', fontsize=12)
ax2.set_title(
    'Expected RDC Distribution for RANDOM Bond Orientations\n'
    '(isotropic sampling, Dₐ=10, R=0) — bimodal at ±2Dₐ and 0',
    fontsize=10
)
ax2.grid(axis='y', ls='--', alpha=0.35)

# Annotation
ax2.annotate('Poles: D = +2Dₐ\n(few vectors, strong coupling)',
             xy=(20, 0.005), xytext=(12, 0.025),
             arrowprops={"arrowstyle": '->', "color": 'black', "lw": 1},
             fontsize=8)
ax2.annotate('Equator: D = −Dₐ\n(many vectors, negative coupling)',
             xy=(-10, 0.025), xytext=(-8.5, 0.055),
             arrowprops={"arrowstyle": '->', "color": 'black', "lw": 1},
             fontsize=8)

plt.tight_layout()
plt.show()

print("\n📊 Reading RDC histograms:")
print("  • Peaks near ±2Dₐ and −Dₐ → random/isotropic orientation mixture")
print("  • peaks shifted from these exact positions → non-random alignment (folded structure)")
print("  • Very narrow distribution → poorly aligned sample (small Da) or disordered protein")
# ============================================================
# Section 5: RDC Distribution Histogram
# ============================================================

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# Compute for a few different Da values
da_vals = [5.0, 10.0, 20.0]
hist_colors = ['#3498db', '#e74c3c', '#2ecc71']

ax = axes[0]
for Da_h, col in zip(da_vals, hist_colors):
    rdcs_h = calculate_rdcs(structure, Da=Da_h, R=0.1)
    if rdcs_h:
        vals_h = list(rdcs_h.values())
        ax.hist(vals_h, bins=14, color=col, alpha=0.65,
                edgecolor='white', linewidth=0.5,
                label=f'Dₐ = {Da_h:.0f} Hz  (range: [{min(vals_h):.1f}, {max(vals_h):.1f}] Hz)')

ax.axvline(0, color='black', lw=0.8, alpha=0.5)
ax.set_xlabel('¹D_{NH} (Hz)', fontsize=12)
ax.set_ylabel('Residue Count', fontsize=12)
ax.set_title('RDC Distribution for Three Dₐ Values\n(R = 0.10, same structure)', fontsize=11)
ax.legend(fontsize=9)
ax.grid(axis='y', ls='--', alpha=0.35)

# ── Right: The sine² angular weighting ────────────────────────────────────
# Explains why purely random orientations are skewed towards negative RDC
ax2 = axes[1]
theta_lin    = np.linspace(0, np.pi, 500)
sinusoid_wgt = np.sin(theta_lin)  # Solid angle element: dΩ ∝ sin(θ)dθ
rdc_lin      = rdc_formula(np.rad2deg(theta_lin), 0, 10.0, 0.0)

# Weighted histogram of RDC values for isotropic orientation distribution
# Each theta contributes weight sin(theta) → probability element on sphere
n_samples = 50000
# Sample theta uniformly in cos(theta) — correct for isotropic distribution
cos_theta  = np.random.uniform(-1, 1, n_samples)
phi_rand   = np.random.uniform(0, 360, n_samples)
rdc_random = rdc_formula(np.rad2deg(np.arccos(cos_theta)), phi_rand, 10.0, 0.0)

ax2.hist(rdc_random, bins=60, color='#9b59b6', alpha=0.7,
         edgecolor='white', linewidth=0.3, density=True)
ax2.axvline(0, color='black', lw=0.8, alpha=0.5)
ax2.set_xlabel('¹D_{NH} (Hz)', fontsize=12)
ax2.set_ylabel('Probability Density', fontsize=12)
ax2.set_title(
    'Expected RDC Distribution for RANDOM Bond Orientations\n'
    '(isotropic sampling, Dₐ=10, R=0) — bimodal at ±2Dₐ and 0',
    fontsize=10
)
ax2.grid(axis='y', ls='--', alpha=0.35)

# Annotation
ax2.annotate('Poles: D = +2Dₐ\n(few vectors, strong coupling)',
             xy=(20, 0.005), xytext=(12, 0.025),
             arrowprops={"arrowstyle": '->', "color": 'black', "lw": 1},
             fontsize=8)
ax2.annotate('Equator: D = −Dₐ\n(many vectors, negative coupling)',
             xy=(-10, 0.025), xytext=(-8.5, 0.055),
             arrowprops={"arrowstyle": '->', "color": 'black', "lw": 1},
             fontsize=8)

plt.tight_layout()
plt.show()

print("\n📊 Reading RDC histograms:")
print("  • Peaks near ±2Dₐ and −Dₐ → random/isotropic orientation mixture")
print("  • peaks shifted from these exact positions → non-random alignment (folded structure)")
print("  • Very narrow distribution → poorly aligned sample (small Da) or disordered protein")

📋 Summary and Key Takeaways¶

Concept	In Plain Language
Alignment tensor	A mathematical description of how the protein prefers to orient in the magnetic field
Dₐ (axial component)	How strongly the protein is aligned — controls the scale of all RDC values
R (rhombicity)	Whether the alignment is symmetric around one axis (R=0) or has an extra twist (R>0)
RDC value sign	Positive = N–H mostly parallel to Z-axis; Negative = N–H mostly perpendicular
Magic angle	θ ≈ 54.7° → zero RDC regardless of Da (when R=0)
Proline	No backbone N–H → always absent from RDC data
RDC vs NOE	RDCs are global (orientation relative to a common frame); NOEs are local (short distances)

📚 Key References¶

Tjandra, N. & Bax, A. (1997). Science 278, 1111. DOI: 10.1126/science.278.5340.1111 — Original RDC measurement in liquid crystals
Prestegard, J.H. et al. (2000). Q Rev Biophys 33, 371. DOI: 10.1017/S0033583500003656 — Comprehensive RDC review
Bax, A. & Grishaev, A. (2005). Curr Opin Struct Biol 15, 563. DOI: 10.1016/j.sbi.2005.08.006 — Multiple media strategies
Bewley, C.A. & Clore, G.M. (2000). J Am Chem Soc 122, 6009. — RDC-only structure refinement
Hansen, M.R. et al. (2000). J Biomol NMR 14, 85. — Pf1 phage alignment medium

🔗 Next Steps¶

virtual_nmr_spectrometer.ipynb — see how RDCs fit into the full NMR observable pipeline alongside relaxation, HSQC, and J-couplings
--output-rdcs CLI flag — generate RDC CSVs for your own sequences: synth-pdb --sequence ACDEF --output-rdcs rdcs.csv --rdc-da 12 --rdc-r 0.15
synth_pdb.rdc API — call calculate_rdcs(structure, Da=10, R=0.1) directly in Python