Computational Physics Guide

Apply computational methods to physics research, including molecular dynamics, Monte Carlo simulations, quantum computing, and numerical methods for solving physical systems.

Computational Methods Overview

Method	Application	Scale	Key Software
Molecular Dynamics (MD)	Atomic-scale dynamics, materials	Atoms-molecules	LAMMPS, GROMACS, NAMD
Density Functional Theory (DFT)	Electronic structure, quantum chemistry	Electrons	VASP, Gaussian, Quantum ESPRESSO
Monte Carlo (MC)	Statistical mechanics, phase transitions	Configurable	Custom, CASINO
Finite Element Method (FEM)	Continuum mechanics, electrostatics	Macroscopic	COMSOL, FEniCS, Abaqus
Finite Difference (FDTD)	Electrodynamics, wave propagation	Macroscopic	Meep, Lumerical
N-body Simulation	Gravitational dynamics, plasma	Stars/particles	GADGET, REBOUND
Lattice QCD	Quantum chromodynamics	Subatomic	MILC, openQCD

Molecular Dynamics

Basic MD Algorithm

import numpy as np

def lennard_jones(r, epsilon=1.0, sigma=1.0):
    """Lennard-Jones potential and force."""
    r6 = (sigma / r) ** 6
    r12 = r6 ** 2
    potential = 4 * epsilon * (r12 - r6)
    force = 24 * epsilon * (2 * r12 - r6) / r
    return potential, force

def velocity_verlet(positions, velocities, forces, masses, dt):
    """Velocity Verlet integration step."""
    # Half-step velocity update
    velocities += 0.5 * forces / masses * dt
    # Full-step position update
    positions += velocities * dt
    # Compute new forces
    new_forces = compute_forces(positions)
    # Complete velocity update
    velocities += 0.5 * new_forces / masses * dt
    return positions, velocities, new_forces

def md_simulation(n_atoms, n_steps, dt=0.001, temperature=1.0):
    """Simple NVE molecular dynamics simulation."""
    # Initialize positions on a grid
    positions = initialize_fcc_lattice(n_atoms, box_size=10.0)
    velocities = np.random.randn(n_atoms, 3) * np.sqrt(temperature)
    velocities -= velocities.mean(axis=0)  # Remove center of mass motion

    forces = compute_forces(positions)
    trajectory = []

    for step in range(n_steps):
        positions, velocities, forces = velocity_verlet(
            positions, velocities, forces,
            masses=np.ones(n_atoms), dt=dt
        )
        if step % 100 == 0:
            ke = 0.5 * np.sum(velocities**2)
            pe = compute_potential_energy(positions)
            print(f"Step {step}: KE={ke:.4f}, PE={pe:.4f}, Total={ke+pe:.4f}")
            trajectory.append(positions.copy())

    return trajectory

LAMMPS Input Script Example

# LAMMPS input: Lennard-Jones fluid simulation
units           lj
atom_style      atomic
boundary        p p p

# Create simulation box and atoms
lattice         fcc 0.8442
region          box block 0 10 0 10 0 10
create_box      1 box
create_atoms    1 box

# Set mass and interactions
mass            1 1.0
pair_style      lj/cut 2.5
pair_coeff      1 1 1.0 1.0 2.5

# Initialize velocities at T=1.0
velocity        all create 1.0 87287 dist gaussian

# Thermostat: Nose-Hoover NVT
fix             1 all nvt temp 1.0 1.0 0.1

# Output settings
thermo          100
thermo_style    custom step temp pe ke etotal press
dump            1 all custom 1000 trajectory.lammpstrj id x y z vx vy vz

# Run simulation
timestep        0.005
run             100000

Monte Carlo Methods

Metropolis Algorithm for Ising Model

import numpy as np

def ising_monte_carlo(L, temperature, n_steps):
    """2D Ising model simulation using Metropolis algorithm."""
    # Initialize random spin configuration
    spins = np.random.choice([-1, 1], size=(L, L))
    beta = 1.0 / temperature

    energies = []
    magnetizations = []

    for step in range(n_steps):
        for _ in range(L * L):  # One sweep = L^2 single spin flips
            # Choose random spin
            i, j = np.random.randint(0, L, size=2)

            # Calculate energy change for flipping spin (i,j)
            neighbors = (
                spins[(i+1)%L, j] + spins[(i-1)%L, j] +
                spins[i, (j+1)%L] + spins[i, (j-1)%L]
            )
            delta_E = 2 * spins[i, j] * neighbors

            # Metropolis acceptance criterion
            if delta_E <= 0 or np.random.random() < np.exp(-beta * delta_E):
                spins[i, j] *= -1

        # Measure observables
        if step % 10 == 0:
            E = -np.sum(spins * (np.roll(spins, 1, 0) + np.roll(spins, 1, 1)))
            M = np.abs(np.sum(spins))
            energies.append(E / L**2)
            magnetizations.append(M / L**2)

    return energies, magnetizations

# Run near the critical temperature (T_c ≈ 2.269 for 2D Ising)
E, M = ising_monte_carlo(L=32, temperature=2.269, n_steps=10000)
print(f"Mean energy: {np.mean(E[-100:]):.4f}")
print(f"Mean magnetization: {np.mean(M[-100:]):.4f}")

Density Functional Theory

Quantum ESPRESSO Workflow

# Step 1: Self-consistent field (SCF) calculation
cat > si_scf.in << 'EOF'
&CONTROL
  calculation = 'scf'
  prefix = 'silicon'
  outdir = './tmp/'
  pseudo_dir = './pseudo/'
/
&SYSTEM
  ibrav = 2
  celldm(1) = 10.26  ! Lattice constant in Bohr
  nat = 2
  ntyp = 1
  ecutwfc = 30.0     ! Kinetic energy cutoff (Ry)
  ecutrho = 300.0    ! Charge density cutoff (Ry)
/
&ELECTRONS
  conv_thr = 1.0d-8
/
ATOMIC_SPECIES
  Si 28.086 Si.pbe-n-rrkjus_psl.1.0.0.UPF
ATOMIC_POSITIONS crystal
  Si 0.00 0.00 0.00
  Si 0.25 0.25 0.25
K_POINTS automatic
  8 8 8 0 0 0
EOF

pw.x < si_scf.in > si_scf.out

# Step 2: Band structure calculation
# (requires nscf + bands post-processing)

Python Interface (ASE + GPAW)

from ase.build import bulk
from gpaw import GPAW, PW

# Create silicon crystal structure
si = bulk('Si', 'diamond', a=5.43)

# DFT calculation with GPAW
calc = GPAW(mode=PW(300),       # Plane-wave cutoff: 300 eV
            kpts=(8, 8, 8),      # k-point mesh
            xc='PBE',            # Exchange-correlation functional
            txt='si_gpaw.txt')   # Output file

si.calc = calc
energy = si.get_potential_energy()
print(f"Total energy: {energy:.4f} eV")
print(f"Energy per atom: {energy/len(si):.4f} eV")

# Equation of state (find equilibrium lattice constant)
from ase.eos import EquationOfState
volumes, energies = [], []
for a in np.linspace(5.3, 5.6, 10):
    si = bulk('Si', 'diamond', a=a)
    si.calc = GPAW(mode=PW(300), kpts=(8,8,8), xc='PBE', txt=None)
    volumes.append(si.get_volume())
    energies.append(si.get_potential_energy())

eos = EquationOfState(volumes, energies)
v0, e0, B = eos.fit()
print(f"Equilibrium volume: {v0:.2f} A^3, Bulk modulus: {B:.1f} GPa")

Numerical Methods

Solving ODEs (Runge-Kutta)

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

# Example: Damped harmonic oscillator
# m*x'' + gamma*x' + k*x = 0
def damped_oscillator(t, y, gamma=0.1, omega0=1.0):
    x, v = y
    dxdt = v
    dvdt = -2*gamma*v - omega0**2 * x
    return [dxdt, dvdt]

sol = solve_ivp(damped_oscillator, [0, 50], [1.0, 0.0],
                t_eval=np.linspace(0, 50, 1000),
                method='RK45', rtol=1e-10)

plt.plot(sol.t, sol.y[0])
plt.xlabel('Time')
plt.ylabel('Displacement')
plt.title('Damped Harmonic Oscillator')
plt.savefig('oscillator.pdf', dpi=300)

Solving PDEs (Finite Differences)

# 2D Heat equation: du/dt = alpha * (d2u/dx2 + d2u/dy2)
def heat_equation_2d(Nx, Ny, Nt, alpha=0.01, dt=0.001):
    dx = dy = 1.0 / max(Nx, Ny)
    u = np.zeros((Nx, Ny))
    u[Nx//4:3*Nx//4, Ny//4:3*Ny//4] = 1.0  # Initial hot region

    for t in range(Nt):
        u_new = u.copy()
        u_new[1:-1, 1:-1] = u[1:-1, 1:-1] + alpha * dt / dx**2 * (
            u[2:, 1:-1] + u[:-2, 1:-1] + u[1:-1, 2:] + u[1:-1, :-2]
            - 4 * u[1:-1, 1:-1]
        )
        u = u_new
    return u

HPC and Parallelization

Approach	Tool	Best For
Shared memory (threads)	OpenMP	Multi-core CPU parallelism
Distributed memory (MPI)	mpi4py, MPI	Multi-node cluster computing
GPU computing	CUDA, CuPy, JAX	Massively parallel computations
Workflow management	Snakemake, Nextflow	Complex simulation pipelines
Job scheduling	SLURM, PBS	HPC cluster job submission

Research Resources

Resource	Description
arXiv cond-mat	Condensed matter preprints
arXiv hep-lat	Lattice field theory preprints
Journal of Computational Physics	Top computational physics journal
Physical Review E	Statistical, nonlinear, soft matter
Computer Physics Communications	Methods + software papers
NIST databases	Physical constants, atomic data