GistGridPlan
Compute Grid Inhomogeneous Solvation Theory (GIST) analysis using a pre-defined grid.
Constructor
from warp_md import GistGridPlan
plan = GistGridPlan(
grid_origin,
grid_spacing,
grid_dimensions,
solvent_selection,
reference_density=0.0334
)
(x, y, z) coordinates of grid origin in Å
Grid spacing in Å (typically 0.5)
(nx, ny, nz) number of grid points in each dimension
Water molecules to analyze (typically oxygen atoms)
Bulk water density in molecules/Ų (0.0334 for SPC/E at 298K)
run() Method
Returns: Dictionary containing GIST grids:
population: Average water count per voxelenergy: Solvation energy (kcal/mol)entropy: Entropy contribution- Additional thermodynamic properties
Example
import warp_md as wmd
import numpy as np
system = wmd.System("protein_water.prmtop")
traj = wmd.open_trajectory("md.nc", system)
# Select water oxygens
water_o = system.select("resname WAT and name O")
# Define grid around active site
origin = (10.0, 15.0, 20.0)
spacing = 0.5
dimensions = (40, 40, 40) # 20Å × 20Å × 20Å
plan = wmd.GistGridPlan(
origin,
spacing,
dimensions,
water_o,
reference_density=0.0334
)
gist_data = plan.run(traj, system)
# Find high-density water sites
population = gist_data["population"]
high_density = np.where(population > 2.0 * 0.0334)
print(f"Found {len(high_density[0])} conserved water sites")
GIST analysis is computationally intensive. Start with a small grid or reduced trajectory for testing.
GistDirectPlan
Direct GIST calculation without pre-defined grid (automatic grid placement).
Constructor
plan = GistDirectPlan(
solute_selection,
solvent_selection,
spacing=0.5,
padding=3.0
)
Protein or solute atoms to analyze around
Distance in Å to extend grid beyond solute
Example
# Automatic GIST around protein
protein = system.select("protein")
water = system.select("resname WAT and name O")
plan = wmd.GistDirectPlan(protein, water, spacing=0.5, padding=5.0)
gist_results = plan.run(traj, system)
# Export to dx format for visualization
from warp_md.analysis import gist
gist.write_dx("gist_energy.dx", gist_results["energy"])
DensityPlan
Compute 3D density maps for selected atoms.
Constructor
plan = DensityPlan(
selection,
grid_origin,
grid_spacing,
grid_dimensions,
normalize=True
)
Atoms to compute density for
Grid origin (x, y, z) in Å
Spacing between grid points in Å
Normalize density by number of frames
Example
# Water density around protein
protein = system.select("protein")
water_o = system.select("resname WAT and name O")
# Get protein bounding box
coords = system.get_coords()
protein_coords = coords[list(protein.indices)]
min_coord = protein_coords.min(axis=0) - 5.0
max_coord = protein_coords.max(axis=0) + 5.0
spacing = 0.5
dimensions = tuple(((max_coord - min_coord) / spacing).astype(int))
plan = wmd.DensityPlan(
water_o,
tuple(min_coord),
spacing,
dimensions,
normalize=True
)
density = plan.run(traj, system)
# Save as DX file for VMD/PyMOL
import numpy as np
np.save("water_density.npy", density)
WatershellPlan
Identify water molecules in the first hydration shell.
Constructor
plan = WatershellPlan(
target_selection,
water_selection,
cutoff=3.5,
pbc="orthorhombic"
)
Solute atoms to analyze around
Water molecules (typically residues, not atoms)
Distance cutoff in Å for hydration shell
pbc
str
default:"orthorhombic"
Periodic boundary conditions
run() Method
Returns: Time series of hydration shell water counts
Example
# Count first-shell waters around ligand
ligand = system.select("resname LIG")
water = system.select("resname WAT and name O")
plan = wmd.WatershellPlan(
ligand,
water,
cutoff=3.5,
pbc="orthorhombic"
)
hydration_counts = plan.run(traj, system)
import matplotlib.pyplot as plt
plt.plot(hydration_counts)
plt.xlabel("Frame")
plt.ylabel("Hydration shell waters")
plt.title("Ligand Hydration")
plt.show()
print(f"Average hydration: {hydration_counts.mean():.1f} waters")
print(f"Std deviation: {hydration_counts.std():.1f}")
Typical first hydration shell cutoffs:
- Ions: 2.5-3.5 Å
- Proteins: 3.5-4.0 Å
- Small molecules: 3.0-3.5 Å
WaterCountPlan
Count the number of water molecules (or other solvent) within a distance cutoff of a selection.
Constructor:
from warp_md import WaterCountPlan
plan = WaterCountPlan(
solute_selection, # Selection: solute atoms
solvent_selection, # Selection: solvent molecules
cutoff=3.5, # float: distance cutoff (Å)
device="auto" # str: execution device
)
Solute atoms around which to count solvent
Solvent molecules to count
Distance cutoff in Ångströms
ndarray[frames] — Number of solvent molecules within cutoff at each frame
Example:
from warp_md import System, Trajectory, WaterCountPlan
system = System.from_pdb("solvated.pdb")
traj = Trajectory.open_xtc("traj.xtc", system)
protein = system.select("protein")
water = system.select("resname SOL or resname WAT or resname TIP3")
plan = WaterCountPlan(protein, water, cutoff=3.5)
water_counts = plan.run(traj, system)
print(f"Average hydration shell: {water_counts.mean():.1f} ± {water_counts.std():.1f} waters")
VolmapPlan
Create volumetric density map on a 3D grid.
Constructor:
from warp_md import VolmapPlan
plan = VolmapPlan(
selection, # Selection: atoms for density map
origin=(0, 0, 0), # tuple: grid origin (Å)
delta=(0.5, 0.5, 0.5), # tuple: grid spacing (Å)
dimensions=(100, 100, 100), # tuple: grid dimensions
device="auto" # str: execution device
)
Atoms to include in density calculation
Grid origin coordinates in Ångströms
delta
tuple
default:"(0.5, 0.5, 0.5)"
Grid spacing in each dimension (Å)
dimensions
tuple
default:"(100, 100, 100)"
Number of grid points in each dimension
ndarray[dimensions] — 3D density map (normalized counts)
Example:
from warp_md import System, Trajectory, VolmapPlan
import numpy as np
system = System.from_pdb("solvated.pdb")
traj = Trajectory.open_xtc("traj.xtc", system)
water_oxygen = system.select("resname SOL and name OW")
# Create 1 Å resolution density map
plan = VolmapPlan(
water_oxygen,
origin=(-25, -25, -25),
delta=(1.0, 1.0, 1.0),
dimensions=(50, 50, 50)
)
density_map = plan.run(traj, system)
# Find high-density regions
threshold = np.percentile(density_map, 95)
high_density_voxels = np.where(density_map > threshold)
print(f"High-density voxels: {len(high_density_voxels[0])}")