Computational Physics and Simulations

Introduction to Computational Physics

Computational physics occupies the third pillar of scientific inquiry alongside theory and experiment—solving physical problems too complex for analytical treatment using numerical algorithms implemented on computers. The exponential growth of computing power described by Moore's Law (doubling roughly every 18 months, 1965-2015; continuing more slowly through specialised accelerators: GPUs, TPUs, neuromorphic chips) has transformed computational physics from solving simple ODEs to simulating billions of atoms, petabyte fluid dynamics simulations, and machine-learning-accelerated quantum chemistry calculations. Key domains: molecular dynamics (Newton's equations for many-body systems); Monte Carlo methods (statistical sampling of configuration space); finite element/difference methods (discretising PDEs on spatial grids); density functional theory (electronic structure of materials); lattice QCD (quantum chromodynamics on a discrete lattice); and climate/weather simulation.

Modern computational physics increasingly interfaces with machine learning: neural network potentials (NNPs) trained on DFT calculations accelerate molecular dynamics by 10^6× while retaining quantum mechanical accuracy; physics-informed neural networks (PINNs) embed physical constraint equations into neural network training; and reinforcement learning optimises experimental parameters in particle accelerators and plasma experiments. The convergence of HPC (high-performance computing), cloud computing, and AI accelerators is reshaping how physicists interact with computation—enabling citizen science distributed computing (Folding@home: protein folding simulations on donated GPUs reaching exaFLOP scale in COVID-19 response) and open-access simulation notebooks democratising computational methods.

Molecular Dynamics and Monte Carlo

Molecular Dynamics Simulation

Molecular dynamics (MD) numerically integrates Newton's equations of motion for N particles interacting via interatomic potentials: m_i * d^2r_i/dt^2 = F_i = -grad_i U(r_1...r_N). Integration: Verlet algorithm and its velocity-Verlet variant conserve energy to O(dt^4) with time steps ~1 fs for atomic systems; RESPA multiple time-step integration separates fast (bonded) and slow (nonbonded) force evaluations. Classical force fields (AMBER, CHARMM, OPLS for biomolecules; LAMMPS-compatible ReaxFF for reactive chemistry; EAM for metals) provide U(r). Periodic boundary conditions eliminate surface effects for bulk simulations. GROMACS, NAMD, LAMMPS, and OpenMM are leading MD packages—GROMACS benchmarks above 10 microsecond/day for 100,000-atom protein systems on modern GPU nodes. Long-timescale methods: metadynamics (adding history-dependent bias potential discouraging revisiting explored configurations), replica exchange MD, and adaptive sampling accelerate rare events (protein folding, membrane permeation, crystal nucleation) inaccessible to plain MD.

Monte Carlo Methods

Monte Carlo (MC) methods use random sampling to evaluate high-dimensional integrals and statistical mechanical averages. Metropolis algorithm: propose random configuration change; accept with probability min(1, exp(-delta E / k_B T))—generates Boltzmann-distributed configurations for statistical ensemble averages. MC is exact in the ergodic limit and avoids force calculation entirely—advantage for hard-sphere, lattice model, and quantum mechanics applications where forces are undefined or expensive. Quantum MC: variational MC (VMC) and diffusion MC (DMC) solve the many-body Schrödinger equation stochastically—DMC exactly solves ground state for bosons and provides fixed-node approximation for fermions (quantum Monte Carlo for solids, Nobel-adjacent methods in quantum chemistry). Wang-Landau algorithm and nested sampling directly compute density of states enabling thermodynamic free energies, phase transitions, and partition functions. Monte Carlo simulations of the Ising model (spin lattice) reproduce critical phenomena: spontaneous symmetry breaking, diverging correlation length, and universality at the critical point T_c = 2J/(k_B * ln(1+sqrt(2))).

Density Functional Theory

DFT Foundations and Applications

Density functional theory (DFT, Hohenberg-Kohn theorem 1964; Kohn-Sham equations 1965; Nobel Chemistry 1998) maps the interacting many-electron problem to a set of single-particle equations in an effective potential—replacing the 3N-dimensional wavefunction with the 3D electron density n(r) as the fundamental variable. Kohn-Sham equations: (-hbar^2/2m)nabla^2 + V_eff[n(r)]) psi_i = epsilon_i psi_i; V_eff = V_ext + V_Hartree + V_xc (exchange-correlation functional is the unknown approximated term). Exchange-correlation functionals: LDA (local density approximation); GGA (PBE—most widely used); hybrid functionals (B3LYP for chemistry, HSE06 for materials bandgaps); meta-GGA (SCAN); methods beyond DFT (DFT+U for strongly correlated; GW approximation for accurate band gaps; TDDFT for optical properties). VASP, Quantum ESPRESSO, CP2K, and Abinit implement plane-wave basis DFT for periodic solids—calculating band structures, phonons, elastic constants, defect energetics, surface adsorption, and reaction barriers. DFT-based high-throughput materials screening (Materials Project, AFLOW, OQMD databases containing >1 million computed material properties) and machine learning interatomic potentials trained on DFT accelerate materials discovery for batteries, catalysts, and semiconductors.

Finite Element and Grid Methods

Finite element method (FEM) discretises PDEs (electromagnetic, structural, thermal, fluid) on unstructured meshes conforming to complex geometries—solving weak variational forms of governing equations. COMSOL Multiphysics, ANSYS, and Abaqus implement FEM for coupled multi-physics simulations of MEMS devices, electromagnetic scattering, thermal-structural coupling in electronics, and biomechanical loading. Finite difference time domain (FDTD, Yee 1966): discretises Maxwell's equations on staggered E and H grids; updates fields in time steps—used for photonic crystal simulation, antenna design, radar cross section. Spectral methods: expanding solutions in basis functions (Fourier, Chebyshev) achieving exponential convergence for smooth solutions—used in seismology, climate modelling. Lattice Boltzmann method (LBM): mesoscale approach to Navier-Stokes fluid dynamics through streaming and collision of probability distributions on a lattice—efficiently handles complex boundaries for porous media flow, microfluidics, DNS turbulence. GPU acceleration (CUDA, HIP, OpenCL) provides 100-1000× speedup over CPU for structured-grid solvers—enabling real-time computational fluid dynamics for jet engine simulation and automotive aerodynamics in wind tunnel replacement workflows.

Examples and Applications

Example 1: Protein Folding Simulations

Protein folding—the process by which a polypeptide chain adopts its functional 3D structure from a random coil—is a grand challenge computational physics and biochemistry problem. MD simulations of small fast-folding proteins (villin headpiece, 35 residues, ~5 microsecond folding time) were first achieved by Shaw Research DE Shaw's Anton specialised MD supercomputers (2010)—millisecond timescales now accessible for ~100 residue proteins on Anton 2/3. Folding@home distributed computing platform aggregates hundreds of thousands of computers globally to simulate rare folding events—characterized COVID-19 spike protein dynamics enabling therapeutics design access. AlphaFold2 (DeepMind 2021) predicts protein structure from sequence with experimental accuracy (median TM-score >0.9 on CASP14)—not a physics simulation but a learned model trained on PDB structures; follows from 50 years of computational structural biology. AlphaFold-Multimer and AlphaFold3 predict protein-protein, protein-nucleic acid, and protein-ligand complexes—essential for drug discovery reducing experimental structure determination of therapeutic targets.

Example 2: Climate and Weather Simulation

Numerical weather prediction (NWP) and climate models discretise the atmosphere on 3D grids (modern NWP ~1-3 km resolution globally; climate models ~10-100 km), solving primitive equations (simplified Navier-Stokes + thermodynamics) with parameterisation of sub-grid processes (cumulus convection, boundary layer, microphysics, radiation). ECMWF (European Centre for Medium-Range Weather Forecasts) runs the world's leading NWP model—10-day forecasts skillful to day 7 (compared to day 3 in the 1980s) from Cray HPE Nippon supercomputer (~8 exaFLOP). CMIP6 (Coupled Model Intercomparison Project) coordinates climate projections from 50+ Earth System Models globally—underlying the IPCC AR6 assessment report scenarios. Challenges: parameterisation uncertainty (cloud feedbacks dominate uncertainty in equilibrium climate sensitivity range 2.5-4°C/doubling CO2); computational cost of explicit convection at km-scale globally; and data assimilation of 10^7 observations daily (satellite radiances, radiosondes, surface stations) using 4D-Var or ensemble Kalman filter variational methods constraining initial conditions.

Example 3: Lattice QCD

Lattice QCD (quantum chromodynamics) formulates the strong nuclear force on a discrete Euclidean spacetime lattice—enabling non-perturbative calculations of quark-gluon dynamics from first principles. The path integral is evaluated by Monte Carlo importance sampling of gauge field configurations generated by Hybrid Monte Carlo with fermion determinant; observables extracted from correlator functions. Lattice QCD calculates: hadron masses (proton mass to <1% accuracy from quark masses); weak decay constants (f_pi, f_K, f_B) needed to extract CKM matrix elements from experimental rates; hadronic vacuum polarisation contribution to anomalous magnetic moment of the muon (muon g-2 tension with experiment requires ~0.1% precision lattice calculations); and QCD phase diagram (deconfinement transition temperature ~155 MeV from lattice). Exa-scale computing (Frontier, Aurora >10^18 FLOP/s) dramatically accelerates continuum extrapolation and physical quark mass simulations. Lattice QCD validation of the Standard Model relies on multi-petaFLOP allocations at ALCF, NERSC, and INCITE leadership computing facilities worldwide.

Example 4: Neural Network Potentials

Machine learning interatomic potentials (MLIPs, NNPs)—neural networks mapping local chemical environment descriptors to atomic energies and forces—combine DFT accuracy with 10^5-10^6× computational speed of classical force fields. Behler-Parrinello neural network potentials (2007), Gaussian approximation potentials (GAP, 2010, Csányi group Cambridge), and deep potential molecular dynamics (DeePMD, 2018) are leading frameworks. Active learning: iterative cycles of MD simulation, uncertainty detection (by committee disagreement), DFT calculation of uncertain configurations, and retraining extend coverage across configuration space while minimising DFT calculations. MACE (equivariant graph neural network potential, 2022) and CHGNet (2023) achieve near-universal materials potentials trained on 1M+ DFT calculations—enabling MD simulations of phase transitions, battery electrolyte decomposition, zeolite synthesis, and amorphous glass structure with first-principles accuracy. Google DeepMind's GNoME universal potential (2023) predicted 2.2 million stable crystal structures—expanding known stable inorganic crystal dataset from 48,000 to over 400,000 entries—demonstrating AI-accelerated materials discovery at scale unprecedented in human-guided research.

Example 5: Particle Physics Event Generation

Large Hadron Collider (LHC) experiments (ATLAS, CMS, LHCb, ALICE) process ~10^9 proton-proton collisions per second, comparing experimental distributions against Monte Carlo event generator predictions for new physics searches. Event generators: PYTHIA (parton shower, hadronisation, underlying event); HERWIG (angular-ordered parton shower); SHERPA (multi-parton NLO merging); MadGraph5_aMC@NLO (automated matrix element generation to NLO accuracy)—incorporating PDF uncertainties, parton showers, and hadronisation. Full detector simulation using Geant4 (particle tracking through geometry with material interactions: electromagnetic processes, hadronic models, optical photons) generates synthetic data samples statistically matching real detector response—consumed at ~10 petabytes/year at WLCG (Worldwide LHC Computing Grid, 200,000 CPU cores). Deep learning event selection: convolutional networks for jet tagging (boosted W, top, Higgs jet identification improving 10-100× over traditional cut-based methods); graph neural networks treating collision events as particle graphs for end-to-end classification; generative adversarial networks (GANs) for fast detector simulation replacing Geant4 for high-statistics samples.

Example 6: N-Body Gravitational Simulations

N-body gravitational simulations follow the trajectories of N particles (stars, dark matter, gas) under mutual gravitational attraction—fundamental to understanding galaxy formation, cluster dynamics, and planetary systems. Direct N-body summation scales as O(N^2)—feasible to N~10^6 on GPU clusters (NBODY7 code); tree algorithms (Barnes-Hut octree, O(N log N)) and Fast Multipole Methods (O(N)) extend to N~10^9. IllustrisTNG (cosmological simulation, 200 Mpc box, ~2×10^10 particles including hydrodynamics and sub-grid galaxy formation physics) ran on Hazel Hen Stuttgart supercomputer—reproducing galaxy morphology, size-mass relations, and metallicity gradients of observed galaxies. Millennium Simulation (2005, 10^10 dark matter particles) established the Lambda-CDM paradigm connecting initial conditions (CMB anisotropies) to observed galaxy clustering statistics. Stellar dynamics: NBODY-X codes handle stellar evolution, binary interactions, tidal stripping, and post-Newtonian corrections for compact binary mergers in globular clusters—predicting gravitational wave source rates from dynamically assembled binary black holes.

Example 7: Quantum Circuit Simulation

Classical simulation of quantum circuits—essential for verification, benchmarking, and error-correcting code research—faces exponential scaling: exact state vector simulation of N qubits requires 2^N complex amplitudes (N=50 requires 64 PB RAM). Google's quantum supremacy claim (53-qubit Sycamore, 2019) was based on claimed infeasibility for classical supercomputers (3 million seconds predicted)—disputed by IBM's tensor network simulations completing the same task in ~2.5 days. Tensor network methods (matrix product states for 1D circuits, projected entangled pair states for 2D, MERA for critical systems) efficiently represent limited-entanglement quantum states. Clifford circuit simulation (stabiliser formalism): poly(N) classical resources for circuits involving only Clifford gates—enabling fault-tolerant threshold calculation for 10^7-qubit error correction codes. Schrödinger-Feynman hybrid algorithms partition circuits into subcircuits classically combined—extending exact simulation feasibility to ~60 qubits. Simulation tools: IBM Qiskit Aer (GPU-accelerated state vector), Google's qsim, PennyLane for differentiable quantum computing—used daily by quantum algorithm researchers verifying gate decompositions for hardware execution.

Example 8: Finite Element Structural Simulation

Structural finite element analysis (FEA) solves the equations of solid mechanics—stress, strain, displacement—in complex geometries under loading conditions, from tiny MEMS accelerometers to 100-metre wind turbine blades and aircraft fuselages. Static analysis: assembles global stiffness matrix K from element stiffness ke(E, nu, geometry); solves Ku=F for displacement vector u; derives stress fromstrain-displacement. Modal analysis: eigenvalue problem Ku = omega^2 Mu identifies natural frequencies and mode shapes—used to validate aircraft wing flutter margins and avoid resonant excitation from engine harmonics. Topology optimisation (density-based SIMP method): iteratively redistributes material within design domain to minimise compliance or mass while satisfying stress/displacement constraints—generating organic-looking lattice structures manufactured by additive manufacturing 30-60% lighter than solid parts. Crash simulation (explicit time integration, LS-DYNA): simulates automotive frontal crash in ~5 ms event with 10^7 elements at ~1 microsecond time steps—requires 100+ CPU cores for hours; validated against physical crash tests per NCAP safety ratings guiding vehicle crumple zone and airbag design.

Try it live

Everything above runs in your browser — open SPH Fluid and change the parameters while it is running. Nothing is installed, nothing is uploaded, the whole model lives in one tab.