Energy analysis plans are currently in development. This page documents the planned API based on the available Python analysis functions.
Energy Analysis
Warp-MD provides energy analysis capabilities through high-level Python functions. These use the Amber sander engine for energy calculations.
energy_analysis()
Compute potential energy components for each frame.
from warp_md.analysis import energy_analysis
result = energy_analysis(
traj,
system,
prmtop="system.prmtop",
frame_indices=None
)
Amber parameter/topology file
Specific frames to analyze (None = all frames)
Returns
Dictionary containing energy components (all in kcal/mol):
total: Total potential energybond: Bond energyangle: Angle energydihedral: Dihedral energyvdw: Van der Waals energyelec: Electrostatic energyvdw_1_4: 1-4 Van der Waalselec_1_4: 1-4 Electrostatic
Example
import warp_md as wmd
from warp_md.analysis import energy_analysis
import matplotlib.pyplot as plt
system = wmd.System("system.prmtop")
traj = wmd.open_trajectory("md.nc", system)
# Compute energies
energies = energy_analysis(
traj,
system,
prmtop="system.prmtop"
)
# Plot energy components
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes[0, 0].plot(energies["total"])
axes[0, 0].set_title("Total Energy")
axes[0, 0].set_ylabel("Energy (kcal/mol)")
axes[0, 1].plot(energies["vdw"], label="VDW")
axes[0, 1].plot(energies["elec"], label="Electrostatic")
axes[0, 1].set_title("Non-bonded Energies")
axes[0, 1].legend()
axes[1, 0].plot(energies["bond"], label="Bond")
axes[1, 0].plot(energies["angle"], label="Angle")
axes[1, 0].plot(energies["dihedral"], label="Dihedral")
axes[1, 0].set_title("Bonded Energies")
axes[1, 0].legend()
axes[1, 0].set_xlabel("Frame")
axes[1, 1].plot(energies["vdw_1_4"], label="VDW 1-4")
axes[1, 1].plot(energies["elec_1_4"], label="Elec 1-4")
axes[1, 1].set_title("1-4 Interactions")
axes[1, 1].legend()
axes[1, 1].set_xlabel("Frame")
plt.tight_layout()
plt.show()
print("\nEnergy Statistics (kcal/mol):")
for component, values in energies.items():
print(f"{component:12s}: {values.mean():10.2f} ± {values.std():8.2f}")
Energy Decomposition
ene_decomp()
Decompose interaction energy between two selections.
from warp_md.analysis import ene_decomp
result = ene_decomp(
traj,
system,
prmtop="system.prmtop",
selection_a="protein",
selection_b="resname LIG",
frame_indices=None
)
First group (e.g., protein)
Second group (e.g., ligand)
Returns
Dictionary with interaction energies:
vdw: Van der Waals interactionelec: Electrostatic interactiontotal: Total interaction energy
Example: Protein-Ligand Binding Energy
import warp_md as wmd
from warp_md.analysis import ene_decomp
import numpy as np
system = wmd.System("complex.prmtop")
traj = wmd.open_trajectory("md.nc", system)
# Protein-ligand interaction energy
interaction = ene_decomp(
traj,
system,
prmtop="complex.prmtop",
selection_a="protein",
selection_b="resname LIG"
)
vdw = interaction["vdw"]
elec = interaction["elec"]
total = interaction["total"]
print("Protein-Ligand Interaction Energy:")
print(f" VDW: {vdw.mean():8.2f} ± {vdw.std():6.2f} kcal/mol")
print(f" Elec: {elec.mean():8.2f} ± {elec.std():6.2f} kcal/mol")
print(f" Total: {total.mean():8.2f} ± {total.std():6.2f} kcal/mol")
# Plot
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
axes[0].plot(vdw, label="VDW", alpha=0.7)
axes[0].plot(elec, label="Electrostatic", alpha=0.7)
axes[0].plot(total, label="Total", linewidth=2)
axes[0].set_ylabel("Energy (kcal/mol)")
axes[0].legend()
axes[0].set_title("Protein-Ligand Interaction Energy")
axes[0].grid(alpha=0.3)
axes[1].hist(total, bins=50, alpha=0.7)
axes[1].set_xlabel("Total Interaction Energy (kcal/mol)")
axes[1].set_ylabel("Count")
axes[1].axvline(total.mean(), color='r', linestyle='--',
label=f'Mean = {total.mean():.1f}')
axes[1].legend()
plt.tight_layout()
plt.show()
Per-Residue Decomposition
import warp_md as wmd
from warp_md.analysis import ene_decomp
import numpy as np
import pandas as pd
system = wmd.System("complex.prmtop")
traj = wmd.open_trajectory("md.nc", system)
# Get protein residues
atom_table = system.atom_table()
protein_resids = np.unique([
atom_table["resid"][i]
for i in range(len(atom_table["resid"]))
if atom_table["resname"][i] not in ["WAT", "Na+", "Cl-"]
])
# Compute per-residue interaction with ligand
residue_energies = []
for resid in protein_resids:
interaction = ene_decomp(
traj,
system,
prmtop="complex.prmtop",
selection_a=f"resid {resid}",
selection_b="resname LIG"
)
residue_energies.append({
"resid": resid,
"resname": atom_table["resname"][resid-1], # Approximate
"vdw": interaction["vdw"].mean(),
"elec": interaction["elec"].mean(),
"total": interaction["total"].mean()
})
# Create DataFrame
df = pd.DataFrame(residue_energies)
df = df.sort_values("total")
print("\nTop 10 favorable interactions:")
print(df.head(10))
print("\nTop 10 unfavorable interactions:")
print(df.tail(10))
# Plot energy footprint
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 6))
plt.bar(df["resid"], df["total"],
color=['green' if x < 0 else 'red' for x in df["total"]])
plt.xlabel("Residue ID")
plt.ylabel("Interaction Energy (kcal/mol)")
plt.axhline(0, color='k', linestyle='-', linewidth=0.5)
plt.title("Per-Residue Ligand Interaction Energy")
plt.tight_layout()
plt.show()
Advanced: Free Energy Estimation
LIE (Linear Interaction Energy)
from warp_md.analysis import lie
# Requires separate simulations of:
# 1. Ligand in complex (bound state)
# 2. Ligand in water (free state)
system_complex = wmd.System("complex.prmtop")
traj_complex = wmd.open_trajectory("complex_md.nc", system_complex)
system_water = wmd.System("ligand_water.prmtop")
traj_water = wmd.open_trajectory("ligand_md.nc", system_water)
# Calculate LIE binding free energy
dG = lie(
traj_complex,
system_complex,
traj_water,
system_water,
ligand_mask="resname LIG",
alpha=0.18, # VDW scaling
beta=0.50 # Electrostatic scaling
)
print(f"LIE Binding Free Energy: {dG:.2f} kcal/mol")
For rigorous binding free energy calculations, use alchemical methods (TI/FEP) or MM-PBSA. Energy decomposition provides quick estimates and insights into interaction patterns.