Computational Physics Guide

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: 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

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

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

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

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)

# 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