Contact Analysis Plans

from warp_md import HbondPlan

plan = HbondPlan(
    donors,
    acceptors,
    dist_cutoff=3.5
)

# With angle criterion
plan = HbondPlan(
    donors,
    acceptors,
    dist_cutoff=3.5
).with_hydrogens(
    hydrogens,
    angle_cutoff_deg=120.0
)

import warp_md as wmd
import matplotlib.pyplot as plt

system = wmd.System("protein.pdb")
traj = wmd.open_trajectory("md.nc", system)

# Protein backbone hydrogen bonds
donors = system.select("protein and name N")
acceptors = system.select("protein and name O")
hydrogens = system.select("protein and name H")

plan = wmd.HbondPlan(
    donors,
    acceptors,
    dist_cutoff=3.5
).with_hydrogens(
    hydrogens,
    angle_cutoff_deg=120.0
)

hbond_counts = plan.run(traj, system)

# Plot
plt.plot(hbond_counts)
plt.xlabel("Frame")
plt.ylabel("Number of H-bonds")
plt.title("Protein Backbone Hydrogen Bonds")
plt.grid(alpha=0.3)
plt.show()

print(f"Average H-bonds: {hbond_counts.mean():.1f}")
print(f"Std deviation: {hbond_counts.std():.1f}")
print(f"Min: {hbond_counts.min()}, Max: {hbond_counts.max()}")

# Specific interaction analysis
protein_donors = system.select("protein and (name N or name O)")
protein_acceptors = system.select("protein and (name O or name N)")

ligand_donors = system.select("resname LIG and (name N or name O)")
ligand_acceptors = system.select("resname LIG and (name O or name N)")

# Protein donates to ligand
plan1 = wmd.HbondPlan(protein_donors, ligand_acceptors, dist_cutoff=3.5)
hbonds_p_to_l = plan1.run(traj, system)

# Ligand donates to protein
plan2 = wmd.HbondPlan(ligand_donors, protein_acceptors, dist_cutoff=3.5)
hbonds_l_to_p = plan2.run(traj, system)

total_hbonds = hbonds_p_to_l + hbonds_l_to_p

print(f"Protein→Ligand H-bonds: {hbonds_p_to_l.mean():.1f}")
print(f"Ligand→Protein H-bonds: {hbonds_l_to_p.mean():.1f}")
print(f"Total: {total_hbonds.mean():.1f}")

plan = NativeContactsPlan(
    selection_a,
    selection_b,
    reference_mode="first",
    cutoff=4.5,
    pbc="none"
)

# Monitor protein folding/unfolding
protein = system.select("protein and name CA")

# Intra-protein contacts (same selection)
plan = wmd.NativeContactsPlan(
    protein,
    protein,
    reference_mode="first",  # Reference = folded state
    cutoff=8.0
)

native_fraction = plan.run(traj, system)

import matplotlib.pyplot as plt
plt.plot(native_fraction)
plt.xlabel("Frame")
plt.ylabel("Native Contact Fraction")
plt.ylim(0, 1)
plt.axhline(0.7, color='r', linestyle='--', label='Folded threshold')
plt.legend()
plt.title("Protein Folding Progress")
plt.show()

if native_fraction[-1] > 0.7:
    print("Protein remains folded")
else:
    print("Protein has partially unfolded")

# Protein-protein interface
chainA = system.select("protein and chain A")
chainB = system.select("protein and chain B")

plan = wmd.NativeContactsPlan(
    chainA,
    chainB,
    reference_mode="first",
    cutoff=5.0
)

interface_contacts = plan.run(traj, system)

# Identify dissociation events
bound_threshold = 0.5
bound_frames = interface_contacts > bound_threshold

import numpy as np
if bound_frames.all():
    print("Complex remains bound throughout simulation")
else:
    first_unbind = np.where(~bound_frames)[0][0]
    print(f"Complex dissociates at frame {first_unbind}")

plan = MindistPlan(
    selection_a,
    selection_b,
    pbc="orthorhombic"
)

# Minimum distance between protein and membrane
protein = system.select("protein")
membrane = system.select("resname POPC")

plan = wmd.MindistPlan(protein, membrane, pbc="orthorhombic")
min_dist = plan.run(traj, system)

import matplotlib.pyplot as plt
plt.plot(min_dist)
plt.xlabel("Frame")
plt.ylabel("Minimum Distance (Å)")
plt.axhline(3.0, color='r', linestyle='--', label='Contact threshold')
plt.legend()
plt.title("Protein-Membrane Distance")
plt.show()

# Check for binding
contact_frames = np.sum(min_dist < 3.0)
print(f"Protein in contact for {contact_frames}/{len(min_dist)} frames")
print(f"Contact fraction: {contact_frames/len(min_dist)*100:.1f}%")

# Monitor specific ion pair
lys_nz = system.select("resid 42 and name NZ")  # Lysine
glu_oe = system.select("resid 87 and name OE1 OE2")  # Glutamate

plan = wmd.MindistPlan(lys_nz, glu_oe, pbc="none")
salt_bridge_dist = plan.run(traj, system)

# Classify as salt bridge if < 4.0 Å
salt_bridge_formed = salt_bridge_dist < 4.0

import numpy as np
occupancy = np.sum(salt_bridge_formed) / len(salt_bridge_formed)
print(f"Salt bridge occupancy: {occupancy*100:.1f}%")

# Find average distance when formed
formed_dists = salt_bridge_dist[salt_bridge_formed]
if len(formed_dists) > 0:
    print(f"Average distance when formed: {formed_dists.mean():.2f} Å")

import warp_md as wmd
import numpy as np
import matplotlib.pyplot as plt

system = wmd.System("protein.pdb")
traj = wmd.open_trajectory("md.xtc", system)

# Get C-alpha atoms
ca_atoms = system.select("name CA")
n_residues = len(ca_atoms.indices)

# Compute average contact map
contact_map = np.zeros((n_residues, n_residues))
cutoff = 8.0  # Å

for i in range(n_residues):
    for j in range(i+1, n_residues):
        sel_i = system.select_indices([ca_atoms.indices[i]])
        sel_j = system.select_indices([ca_atoms.indices[j]])
        
        plan = wmd.MindistPlan(sel_i, sel_j, pbc="none")
        distances = plan.run(traj, system)
        
        # Fraction of frames in contact
        contact_fraction = np.sum(distances < cutoff) / len(distances)
        contact_map[i, j] = contact_fraction
        contact_map[j, i] = contact_fraction

# Plot
plt.figure(figsize=(10, 8))
im = plt.imshow(contact_map, cmap='YlOrRd', vmin=0, vmax=1)
plt.colorbar(im, label='Contact Probability')
plt.xlabel("Residue Index")
plt.ylabel("Residue Index")
plt.title(f"Contact Map (cutoff = {cutoff} Å)")
plt.tight_layout()
plt.show()