Skip to main content

​
PcaPlan

Perform principal component analysis on trajectory coordinates.

​
Constructor

from warp_md import PcaPlan

plan = PcaPlan(
    selection,
    n_vecs=2,
    fit=True,
    ref_selection=None
)
​
selection
Selection
required
Atoms to include in PCA calculation
​
n_vecs
int
default:"2"
Number of principal components to compute
​
fit
bool
default:"True"
Align structures before PCA
​
ref_selection
Selection
Atoms to use for alignment (if different from PCA selection)

​
run() Method

Returns: Tuple of (projections, (eigenvalues, eigenvectors))
  • projections: 2D array of shape (n_vecs, n_frames)
  • eigenvalues: 1D array of length n_vecs
  • eigenvectors: 2D array of shape (n_vecs, 3*n_atoms)

​
Example

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)

# PCA on C-alpha atoms
ca_atoms = system.select("name CA")
plan = wmd.PcaPlan(ca_atoms, n_vecs=3, fit=True)
proj, (evals, evecs) = plan.run(traj, system)

# Plot first two principal components
plt.figure(figsize=(8, 6))
plt.scatter(proj[0], proj[1], c=range(len(proj[0])), cmap='viridis')
plt.xlabel(f"PC1 ({evals[0]/evals.sum()*100:.1f}%)")
plt.ylabel(f"PC2 ({evals[1]/evals.sum()*100:.1f}%)")
plt.colorbar(label="Frame")
plt.title("PCA Projection")
plt.show()

# Variance explained
var_explained = evals / evals.sum() * 100
print(f"Variance explained by first 3 PCs: {var_explained.sum():.1f}%")
Eigenvalues represent the variance captured by each principal component.

​
AnalyzeModesPlan

Analyze and visualize principal components or normal modes.

​
Constructor

plan = AnalyzeModesPlan(
    selection,
    eigenvectors,
    eigenvalues,
    mode_index=0
)
​
selection
Selection
required
Atoms used in mode calculation
​
eigenvectors
ndarray
required
Eigenvectors from PCA or normal mode analysis
​
eigenvalues
ndarray
required
Corresponding eigenvalues
​
mode_index
int
default:"0"
Which mode to analyze (0-indexed)

​
Example

# First perform PCA
ca_atoms = system.select("name CA")
pca_plan = wmd.PcaPlan(ca_atoms, n_vecs=5)
proj, (evals, evecs) = pca_plan.run(traj, system)

# Analyze first principal component
mode_plan = wmd.AnalyzeModesPlan(
    ca_atoms,
    evecs,
    evals,
    mode_index=0
)

# Generate structures along PC1
mode_data = mode_plan.run(traj, system)

​
ProjectionPlan

Project trajectory onto pre-computed eigenvectors.

​
Constructor

plan = ProjectionPlan(
    selection,
    eigenvectors,
    eigenvalues,
    average_coords=None
)
​
selection
Selection
required
Atoms to project
​
eigenvectors
ndarray
required
Eigenvectors to project onto (from previous PCA)
​
eigenvalues
ndarray
required
Corresponding eigenvalues
​
average_coords
ndarray
Reference coordinates for centering (if not provided, uses trajectory average)

​
run() Method

Returns: 2D array of shape (n_vecs, n_frames) with projection coefficients

​
Example

# Train PCA on equilibrated portion
system = wmd.System("system.pdb")
traj_train = wmd.open_trajectory("equilibration.xtc", system)
traj_prod = wmd.open_trajectory("production.xtc", system)

ca_atoms = system.select("name CA")

# Compute PCA on training trajectory
pca_plan = wmd.PcaPlan(ca_atoms, n_vecs=5)
proj_train, (evals, evecs) = pca_plan.run(traj_train, system)

# Project production trajectory onto same PC space
proj_plan = wmd.ProjectionPlan(ca_atoms, evecs, evals)
proj_prod = proj_plan.run(traj_prod, system)

# Compare distributions
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].hist(proj_train[0], bins=30, alpha=0.7, label="Training")
axes[0].hist(proj_prod[0], bins=30, alpha=0.7, label="Production")
axes[0].set_xlabel("PC1")
axes[0].legend()

axes[1].scatter(proj_train[0], proj_train[1], alpha=0.5, label="Training")
axes[1].scatter(proj_prod[0], proj_prod[1], alpha=0.5, label="Production")
axes[1].set_xlabel("PC1")
axes[1].set_ylabel("PC2")
axes[1].legend()
plt.show()

​
Python Helper Functions

from warp_md.analysis import pca, projection

# High-level PCA interface
proj, (evals, evecs) = pca(
    traj,
    system,
    mask="name CA",
    n_vecs=3,
    fit=True
)

# Project onto existing modes
proj_new = projection(
    traj_new,
    system,
    mask="name CA",
    eigenvectors=evecs,
    eigenvalues=evals
)
PCA is powerful for identifying collective motions and essential dynamics. Use the first few PCs (typically 2-5) to capture the dominant conformational changes.

​
MatrixPlan

Calculate covariance or correlation matrix from trajectory coordinates. Constructor: 
from warp_md import MatrixPlan

plan = MatrixPlan(
    selection,           # Selection: atoms for matrix
    matrix_type="covar", # str: "covar" or "correl"
    mass_weighted=False, # bool: mass weighting
    fit=True,           # bool: align before calculation
    device="auto"        # str: execution device
)
​
selection
Selection
required
Atoms to include in the covariance/correlation matrix
​
matrix_type
str
default:"covar"
Matrix type: "covar" (covariance) or "correl" (correlation)
​
mass_weighted
bool
default:"False"
Whether to use mass-weighted coordinates
​
fit
bool
default:"True"
Whether to align structures before calculation
Returns: ndarray[3N, 3N] — Covariance or correlation matrix (N = number of atoms) Example: 
from warp_md import System, Trajectory, MatrixPlan
import numpy as np

system = System.from_pdb("protein.pdb")
traj = Trajectory.open_xtc("traj.xtc", system)
ca_atoms = system.select("name CA")

# Calculate covariance matrix
plan = MatrixPlan(ca_atoms, matrix_type="covar", fit=True)
covar_matrix = plan.run(traj, system)

# Eigenvalue decomposition for essential dynamics
eigenvalues, eigenvectors = np.linalg.eigh(covar_matrix)
print(f"Top 3 eigenvalues: {eigenvalues[-3:][::-1]}")

Build docs developers (and LLMs) love

Get started for freeTalk to us