Polymer Analysis Plans

from warp_md import EndToEndPlan

plan = EndToEndPlan(selection)

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

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

# Select polymer backbone
backbone = system.select("name CA or name C or name N")

plan = wmd.EndToEndPlan(backbone)
ree = plan.run(traj, system)

# If multiple chains
n_chains = ree.shape[1]
for i in range(n_chains):
    plt.plot(ree[:, i], label=f"Chain {i+1}", alpha=0.7)

plt.xlabel("Frame")
plt.ylabel("End-to-End Distance (Å)")
plt.legend()
plt.title("Polymer Extension")
plt.show()

print(f"Average R_ee: {ree.mean(axis=0)}")
print(f"Std R_ee: {ree.std(axis=0)}")

plan = ChainRgPlan(selection, mass_weighted=True)

# Compare Rg and end-to-end distance
backbone = system.select("polymer")

ree_plan = wmd.EndToEndPlan(backbone)
rg_plan = wmd.ChainRgPlan(backbone, mass_weighted=True)

ree = ree_plan.run(traj, system)
rg = rg_plan.run(traj, system)

# Plot relationship
import matplotlib.pyplot as plt
plt.scatter(ree.flatten(), rg.flatten(), alpha=0.3)
plt.xlabel("End-to-End Distance (Å)")
plt.ylabel("Radius of Gyration (Å)")
plt.title("Polymer Conformation")

# Theoretical line for ideal chain: Rg = Ree/sqrt(6)
ree_range = np.linspace(ree.min(), ree.max(), 100)
plt.plot(ree_range, ree_range / np.sqrt(6), 'r--', label='Ideal chain')
plt.legend()
plt.show()

print(f"Rg/Ree ratio: {(rg.mean() / ree.mean()):.3f}")
print(f"Expected for ideal chain: {1/np.sqrt(6):.3f}")

plan = PersistenceLengthPlan(
    selection,
    bond_length=3.8
)

# Calculate persistence length
backbone = system.select("name CA")

plan = wmd.PersistenceLengthPlan(
    backbone,
    bond_length=3.8  # C-alpha spacing
)

result = plan.run(traj, system)

lp = result["persistence_length"]
corr = result["bond_correlations"]

print(f"Persistence length: {lp:.2f} Å")

# Plot correlation decay
import matplotlib.pyplot as plt
import numpy as np

separation = np.arange(len(corr)) * 3.8
plt.semilogy(separation, corr, 'o-')
plt.xlabel("Contour separation (Å)")
plt.ylabel("Bond correlation")
plt.axhline(np.exp(-1), color='r', linestyle='--', label='1/e')
plt.axvline(lp, color='r', linestyle='--', label=f'Lp = {lp:.1f} Å')
plt.legend()
plt.title("Persistence Length Analysis")
plt.show()

# Calculate Kuhn length from persistence length and contour length
result = plan.run(traj, system)
lp = result["persistence_length"]

# Kuhn length b = 2 * Lp
kuhn_length = 2 * lp

# Contour length (total chain length)
backbone_indices = list(backbone.indices)
n_bonds = len(backbone_indices) - 1
contour_length = n_bonds * 3.8  # bond_length

# Number of Kuhn segments
n_kuhn = contour_length / kuhn_length

print(f"Kuhn length: {kuhn_length:.2f} Å")
print(f"Contour length: {contour_length:.2f} Å")
print(f"Number of Kuhn segments: {n_kuhn:.1f}")

# Expected Rg for Gaussian chain
rg_expected = kuhn_length * np.sqrt(n_kuhn / 6)
print(f"Expected Rg (Gaussian chain): {rg_expected:.2f} Å")

import warp_md as wmd
import numpy as np

def classify_polymer(system, traj, selection):
    """Classify polymer behavior from scaling analysis."""
    
    # Calculate Rg and end-to-end distance
    rg_plan = wmd.ChainRgPlan(selection)
    ree_plan = wmd.EndToEndPlan(selection)
    
    rg = rg_plan.run(traj, system).mean()
    ree = ree_plan.run(traj, system).mean()
    
    # Calculate persistence length
    lp_plan = wmd.PersistenceLengthPlan(selection, bond_length=3.8)
    result = lp_plan.run(traj, system)
    lp = result["persistence_length"]
    
    # Estimate contour length
    n_atoms = len(selection.indices)
    contour_length = (n_atoms - 1) * 3.8
    
    # Flexibility ratio
    flexibility = contour_length / lp
    
    # Shape factor
    shape = rg / ree
    
    print("=== Polymer Classification ===")
    print(f"Contour length: {contour_length:.1f} Å")
    print(f"Persistence length: {lp:.1f} Å")
    print(f"Flexibility ratio (L/Lp): {flexibility:.2f}")
    print(f"Shape factor (Rg/Ree): {shape:.3f}")
    print()
    
    if flexibility < 1:
        regime = "Rod-like (rigid)"
    elif flexibility < 10:
        regime = "Semi-flexible"
    else:
        regime = "Flexible coil"
    
    print(f"Regime: {regime}")
    
    if 0.38 < shape < 0.42:
        conformation = "Gaussian coil"
    elif shape < 0.38:
        conformation = "Collapsed/Compact"
    else:
        conformation = "Extended"
    
    print(f"Conformation: {conformation}")
    
    return {
        "lp": lp,
        "contour_length": contour_length,
        "flexibility": flexibility,
        "shape_factor": shape,
        "regime": regime,
        "conformation": conformation
    }

# Usage
backbone = system.select("name CA")
result = classify_polymer(system, traj, backbone)