physics-simulation

Physics Simulation

Overview

This skill covers computational physics workflows: from solving differential equations and analyzing simulation trajectories to signal processing and rigorous uncertainty quantification.

When to Use This Skill

• Analyzing or running physical simulations (MD, Monte Carlo, FEM)
• Solving ODEs or PDEs numerically
• Processing experimental physics data (spectra, time series, sensor data)
• Propagating measurement uncertainties
• Working with physical units and dimensional analysis

---

1. ODE Solvers with SciPy

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# Example: damped harmonic oscillator
# y'' + 2γy' + ω₀²y = 0  →  state vector [y, y']
def damped_oscillator(t, y, gamma, omega0):
    return [y[1], -2*gamma*y[1] - omega0**2 * y[0]]

# Solve
sol = solve_ivp(
    fun=damped_oscillator,
    t_span=(0, 20),
    y0=[1.0, 0.0],         # initial displacement, velocity
    args=(0.1, 2.0),       # gamma, omega0
    method="RK45",         # default, good for non-stiff
    t_eval=np.linspace(0, 20, 500),
    rtol=1e-8, atol=1e-10,
)

if not sol.success:
    print(f"Solver failed: {sol.message}")

Method Selection

Method	Use when	Notes
RK45	Non-stiff, smooth solutions	Default, good general purpose
RK23	Non-stiff, less accuracy needed	Faster than RK45
DOP853	Non-stiff, high accuracy required	8th order, fewer function evaluations
Radau	Stiff systems	Chemical kinetics, circuit simulation
BDF	Very stiff, large time spans	Implicit, variable order
LSODA	Unknown stiffness	Auto-switches between stiff/non-stiff

Stiff system detection: If RK45 takes extremely small steps (t_eval coverage is sparse) or max_step warning appears → switch to Radau or BDF.

---

2. Molecular Dynamics

Trajectory Analysis with MDAnalysis

import MDAnalysis as mda
from MDAnalysis.analysis import rms, distances

# Load trajectory
u = mda.Universe("topology.tpr", "trajectory.xtc")

# RMSD relative to first frame
backbone = u.select_atoms("backbone")
R = rms.RMSD(backbone, backbone, ref_frame=0)
R.run()
# R.results.rmsd[:, 2] = RMSD values over time

# Radius of gyration over trajectory
Rg = []
for ts in u.trajectory:
    Rg.append(u.atoms.radius_of_gyration())

# Hydrogen bond analysis
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis
hbonds = HydrogenBondAnalysis(u, "protein", "protein")
hbonds.run()