leonardodalinky/scider/physics-simulation

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}")

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()

# Read LAMMPS dump file
import numpy as np
def read_lammps_dump(filename):
    frames = []
    with open(filename) as f:
        while True:
            try:
                next(f)  # ITEM: TIMESTEP
                timestep = int(next(f).strip())
                next(f)  # ITEM: NUMBER OF ATOMS
                n_atoms = int(next(f).strip())
                for _ in range(3): next(f)  # box bounds lines
                next(f)  # ITEM: ATOMS header
                atoms = np.array([next(f).split() for _ in range(n_atoms)], dtype=float)
                frames.append({"timestep": timestep, "atoms": atoms})
            except StopIteration:
                break
    return frames

import numpy as np
def metropolis_mcmc(energy_func, x_init, n_steps, step_size, temperature=1.0):
    """Metropolis-Hastings MCMC sampler."""
    x = x_init.copy()
    E = energy_func(x)
    samples = [x.copy()]
    n_accepted = 0
    for _ in range(n_steps):
        x_prop = x + np.random.normal(0, step_size, size=x.shape)
        E_prop = energy_func(x_prop)
        delta_E = E_prop - E
        # Accept with probability min(1, exp(-ΔE/kT))
        if delta_E < 0 or np.random.rand() < np.exp(-delta_E / temperature):
            x = x_prop
            E = E_prop
            n_accepted += 1
        samples.append(x.copy())
    acceptance_rate = n_accepted / n_steps
    print(f"Acceptance rate: {acceptance_rate:.2%} (target: 20-50%)")
    return np.array(samples)
# Convergence diagnostics
def autocorrelation_time(samples):
    """Estimate integrated autocorrelation time."""
    n = len(samples)
    mean = samples.mean()
    var = samples.var()
    acf = np.array([np.mean((samples[:n-k] - mean) * (samples[k:] - mean)) / var
                    for k in range(min(100, n//4))])
    # Integrated autocorr time ≈ 1 + 2 * sum(ACF)
    return 1 + 2 * acf[1:][acf[1:] > 0].sum()

from scipy import signal
import numpy as np
# Generate example signal
fs = 1000  # Hz sampling rate
t = np.arange(0, 1, 1/fs)
clean = np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t)
noisy = clean + 2*np.random.randn(len(t))
# FFT analysis
freqs = np.fft.rfftfreq(len(noisy), 1/fs)
fft = np.abs(np.fft.rfft(noisy))
# Butterworth lowpass filter (< 100 Hz)
b, a = signal.butter(N=4, Wn=100, btype="low", fs=fs)
filtered = signal.sosfiltfilt(signal.butter(4, 100, "low", fs=fs, output="sos"), noisy)
# Power spectral density (Welch method)
freq_psd, psd = signal.welch(noisy, fs=fs, nperseg=256)
# Spectrogram
freq_spec, t_spec, Sxx = signal.spectrogram(noisy, fs=fs, nperseg=128)

from uncertainties import ufloat, umath
import uncertainties.unumpy as unp
# Define measurements with uncertainties
length = ufloat(10.5, 0.1)   # 10.5 ± 0.1 m
width = ufloat(3.2, 0.05)    # 3.2 ± 0.05 m
area = length * width        # automatic error propagation
print(f"Area = {area}")      # Area = 33.6 ± 0.6 m²
# Works with functions
import math
radius = ufloat(5.0, 0.2)
circumference = 2 * umath.pi * radius

from scipy.optimize import curve_fit
def model(x, a, b, c):
    return a * np.exp(-b * x) + c
x_data = np.linspace(0, 5, 50)
y_data = 3 * np.exp(-0.7 * x_data) + 0.5 + 0.1 * np.random.randn(50)
y_err = 0.1 * np.ones_like(y_data)
popt, pcov = curve_fit(model, x_data, y_data, sigma=y_err, absolute_sigma=True)
perr = np.sqrt(np.diag(pcov))  # 1-sigma uncertainties on parameters
# Reduced χ²
residuals = y_data - model(x_data, *popt)
chi2 = np.sum((residuals / y_err)**2)
n_dof = len(y_data) - len(popt)
chi2_reduced = chi2 / n_dof
print(f"Reduced χ² = {chi2_reduced:.2f}  (good fit ≈ 1.0)")

from pint import UnitRegistry
ureg = UnitRegistry()
Q = ureg.Quantity
# Define physical quantities
mass = Q(10.0, "kg")
velocity = Q(3.0, "m/s")
kinetic_energy = 0.5 * mass * velocity**2
print(kinetic_energy.to("joule"))         # 45.0 joule
print(kinetic_energy.to("eV"))            # convert to electron volts
# Physical constants
import scipy.constants as const
print(f"Boltzmann: {const.k} J/K")
print(f"Planck: {const.h} J·s")
print(f"Speed of light: {const.c} m/s")
print(f"Avogadro: {const.N_A} mol⁻¹")
# Natural units conversion (particle physics)
# In natural units: ℏ = c = 1, energies in eV
hbar_eV_s = const.hbar / const.e  # ℏ in eV·s

python <SKILL_BASE_DIR>/scripts/physics_eda.py simulation_output.csv --output physics_report.md