ConductivityPlan
Compute ionic conductivity from mean-square displacement of ions.
Constructor
from warp_md import ConductivityPlan
plan = ConductivityPlan(
cation_selection,
anion_selection,
charges,
temperature=298.15
)
Cation atoms (e.g., Na+, K+)
Dictionary mapping atom indices to charges: {atom_idx: charge}
run() Method
Returns: Conductivity value in S/m (Siemens per meter)
Example
import warp_md as wmd
import numpy as np
system = wmd.System("electrolyte.prmtop")
traj = wmd.open_trajectory("md.nc", system)
# Select ions
sodium = system.select("resname Na+")
chloride = system.select("resname Cl-")
# Assign charges
charges = {}
for idx in sodium.indices:
charges[idx] = 1.0
for idx in chloride.indices:
charges[idx] = -1.0
# Calculate conductivity
plan = wmd.ConductivityPlan(
sodium,
chloride,
charges,
temperature=298.15
)
conductivity = plan.run(traj, system)
print(f"Ionic conductivity: {conductivity:.3f} S/m")
Conductivity is calculated using the Nernst-Einstein relation from ion mean-square displacements.
DielectricPlan
Compute static dielectric constant from dipole fluctuations.
Constructor
plan = DielectricPlan(
selection,
temperature=298.15,
volume=None
)
All atoms to include in dipole calculation
System volume in Ų (if None, extracted from trajectory box)
run() Method
Returns: Static dielectric constant ε₀ (dimensionless)
Example
# Calculate dielectric constant of water
water = system.select("resname WAT")
plan = wmd.DielectricPlan(
water,
temperature=298.15
)
epsilon = plan.run(traj, system)
print(f"Dielectric constant: {epsilon:.2f}")
print(f"Literature value for SPC/E water: ~71")
The dielectric constant is computed using:
ε = 1 + (4π/3kTV) ⟨|ΔM|²⟩
where M is the total system dipole moment.
DipoleAlignmentPlan
Analyze dipole moment alignment and orientational order.
Constructor
plan = DipoleAlignmentPlan(
selection,
reference_vector=None
)
Reference direction (x, y, z). If None, uses system dipole
run() Method
Returns: 2D array of shape (n_frames, n_metrics) containing:
- Column 0: Average alignment (⟨cos θ⟩)
- Column 1: Order parameter S = (3⟨cos²θ⟩ - 1)/2
- Column 2: Total dipole magnitude
Example
# Analyze water dipole alignment in electric field
water = system.select("resname WAT")
# Reference = z-axis (field direction)
plan = wmd.DipoleAlignmentPlan(
water,
reference_vector=(0, 0, 1)
)
alignment_data = plan.run(traj, system)
avg_alignment = alignment_data[:, 0]
order_param = alignment_data[:, 1]
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 1, figsize=(8, 6))
axes[0].plot(avg_alignment)
axes[0].set_ylabel("⟨cos θ⟩")
axes[0].axhline(0, color='k', linestyle='--', alpha=0.3)
axes[0].set_title("Dipole Alignment")
axes[1].plot(order_param)
axes[1].set_ylabel("Order Parameter S")
axes[1].set_xlabel("Frame")
axes[1].axhline(0, color='k', linestyle='--', alpha=0.3)
plt.tight_layout()
plt.show()
print(f"Mean alignment: {avg_alignment.mean():.3f}")
print(f"Mean order parameter: {order_param.mean():.3f}")
Order parameter interpretation:
- S = 1: Perfect alignment
- S = 0: Random orientation
- S = -0.5: Perpendicular orientation
Advanced Example: Dipole Correlation
# Calculate dipole-dipole correlation function
from warp_md.analysis import timecorr
water = system.select("resname WAT")
# Get dipole vectors
plan = wmd.DipoleAlignmentPlan(water)
dipole_data = plan.run(traj, system)
# Extract time series of total dipole
total_dipole = dipole_data[:, 2]
# Compute autocorrelation
corr = timecorr(total_dipole, max_lag=1000)
import matplotlib.pyplot as plt
plt.plot(corr)
plt.xlabel("Lag (frames)")
plt.ylabel("Dipole Autocorrelation")
plt.title("Dipole Relaxation")
plt.show()
IonPairCorrelationPlan
Calculate ion-pair radial distribution function for electrolyte systems.
Constructor:
from warp_md import IonPairCorrelationPlan
plan = IonPairCorrelationPlan(
cation_selection, # Selection: cation atoms
anion_selection, # Selection: anion atoms
bins=200, # int: number of histogram bins
r_max=10.0, # float: maximum distance (Å)
device="auto" # str: execution device
)
Maximum distance in Ångströms
Returns: (r, g_r) — tuple of distance bins and correlation function
StructureFactorPlan
Calculate static structure factor S(q) from atomic positions.
Constructor:
from warp_md import StructureFactorPlan
plan = StructureFactorPlan(
selection, # Selection: atoms for S(q)
q_max=20.0, # float: maximum q value (Å⁻¹)
q_bins=200, # int: number of q bins
device="auto" # str: execution device
)
Atoms to include in structure factor calculation
Maximum wavevector magnitude in Å⁻¹
Returns: (q, S_q) — tuple of wavevectors and structure factor
EquipartitionPlan
Verify energy equipartition theorem for temperature validation.
Constructor:
from warp_md import EquipartitionPlan
plan = EquipartitionPlan(
selection, # Selection: atoms to analyze
temperature=300.0, # float: expected temperature (K)
device="auto" # str: execution device
)
Atoms for equipartition analysis
Expected system temperature in Kelvin
Returns: ndarray[frames, 3] — Kinetic energy per degree of freedom for each frame
Example:
from warp_md import System, Trajectory, EquipartitionPlan
system = System.from_pdb("md.pdb")
traj = Trajectory.open_dcd("md.dcd", system)
all_atoms = system.select("all")
plan = EquipartitionPlan(all_atoms, temperature=300.0)
ke_per_dof = plan.run(traj, system)
# Expected: 0.5 * k_B * T per DOF
expected = 0.5 * 0.001987 * 300.0 # kcal/mol
print(f"Mean KE per DOF: {ke_per_dof.mean():.4f} kcal/mol")
print(f"Expected: {expected:.4f} kcal/mol")