Skip to main content

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_selection
Selection
required
Cation atoms (e.g., Na+, K+)
anion_selection
Selection
required
Anion atoms (e.g., Cl-)
charges
dict
required
Dictionary mapping atom indices to charges: {atom_idx: charge}
temperature
float
default:"298.15"
Temperature in Kelvin

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
)
selection
Selection
required
All atoms to include in dipole calculation
temperature
float
default:"298.15"
Temperature in Kelvin
volume
float
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")

Formula

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
)
selection
Selection
required
Molecules to analyze
reference_vector
tuple
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
)
cation_selection
Selection
required
Cation atoms
anion_selection
Selection
required
Anion atoms
bins
int
default:"200"
Number of histogram bins
r_max
float
default:"10.0"
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
)
selection
Selection
required
Atoms to include in structure factor calculation
q_max
float
default:"20.0"
Maximum wavevector magnitude in Å⁻¹
q_bins
int
default:"200"
Number of q bins
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
)
selection
Selection
required
Atoms for equipartition analysis
temperature
float
default:"300.0"
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")

Build docs developers (and LLMs) love