🧠

Molecular Dynamics Simulations

Simulating the motion of atoms and molecules by integrating Newton's equations. From classical force fields to enhanced sampling and QM/MM hybrid methods — the computational microscope for biomolecular systems.

1. Classical Molecular Dynamics

Molecular dynamics solves Newton's equations of motion for a system of N interacting particles. For each particle i with mass $m_i$ at position $\mathbf{r}_i$:

$$m_i \frac{d^2\mathbf{r}_i}{dt^2} = \mathbf{F}_i = -\nabla_i U(\mathbf{r}_1, \mathbf{r}_2, \ldots, \mathbf{r}_N)$$

Velocity Verlet Integration

The velocity Verlet algorithm is the most widely used integrator due to its symplectic nature, time-reversibility, and excellent energy conservation:

$$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\Delta t + \frac{1}{2}\frac{\mathbf{F}(t)}{m}\Delta t^2$$

$$\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{\Delta t}{2m}\left[\mathbf{F}(t) + \mathbf{F}(t + \Delta t)\right]$$

Typical timesteps are 1-2 fs for atomistic simulations. The SHAKE/LINCS algorithms constrain bond lengths involving hydrogen, allowing timesteps up to 2 fs. Hydrogen mass repartitioning (HMR) can extend this to 4 fs.

Derivation: Newton's Equations and the Verlet Integration Algorithm

We derive the velocity Verlet algorithm from Taylor expansions and show why it conserves energy over long trajectories.

Step 1: Start from Newton's second law

For each atom i in the system, the force is the negative gradient of the potential energy surface:

$$m_i \ddot{\mathbf{r}}_i = \mathbf{F}_i = -\nabla_i V(\mathbf{r}_1, \ldots, \mathbf{r}_N)$$

This gives 3N coupled second-order ODEs that must be integrated numerically.

Step 2: Taylor expand positions forward and backward in time

$$\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \dot{\mathbf{r}}(t)\Delta t + \frac{1}{2}\ddot{\mathbf{r}}(t)\Delta t^2 + \frac{1}{6}\dddot{\mathbf{r}}(t)\Delta t^3 + \mathcal{O}(\Delta t^4)$$

$$\mathbf{r}(t - \Delta t) = \mathbf{r}(t) - \dot{\mathbf{r}}(t)\Delta t + \frac{1}{2}\ddot{\mathbf{r}}(t)\Delta t^2 - \frac{1}{6}\dddot{\mathbf{r}}(t)\Delta t^3 + \mathcal{O}(\Delta t^4)$$

Step 3: Add the two expansions (original Verlet)

The odd-order terms cancel, giving the position Verlet algorithm with local error $\mathcal{O}(\Delta t^4)$:

$$\mathbf{r}(t + \Delta t) = 2\mathbf{r}(t) - \mathbf{r}(t - \Delta t) + \frac{\mathbf{F}(t)}{m}\Delta t^2$$

Step 4: Reformulate as velocity Verlet

The velocity Verlet variant explicitly tracks velocities, improving numerical stability and enabling kinetic energy calculations:

$$\mathbf{r}(t+\Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\Delta t + \frac{\mathbf{F}(t)}{2m}\Delta t^2$$

$$\mathbf{v}(t+\Delta t) = \mathbf{v}(t) + \frac{\Delta t}{2m}\bigl[\mathbf{F}(t) + \mathbf{F}(t+\Delta t)\bigr]$$

Step 5: Symplecticity and energy conservation

The Verlet integrator is symplectic: it preserves the phase-space volume (Liouville's theorem) and conserves a shadow Hamiltonian $\tilde{H} = H + \mathcal{O}(\Delta t^2)$. This means total energy fluctuates but does not drift over time, making it ideal for long MD trajectories spanning nanoseconds to microseconds.

Derivation: Equipartition Theorem and Temperature in MD

We derive the relationship between kinetic energy and temperature that allows MD simulations to connect microscopic velocities to macroscopic thermodynamic temperature.

Step 1: Classical equipartition theorem

In classical statistical mechanics, each quadratic degree of freedom in the Hamiltonian contributes $\frac{1}{2}k_BT$ to the average energy. For a quadratic term $ax_i^2$ in the Hamiltonian:

$$\langle ax_i^2 \rangle = \frac{\int ax_i^2 e^{-ax_i^2/k_BT} dx_i}{\int e^{-ax_i^2/k_BT} dx_i} = \frac{1}{2}k_BT$$

Step 2: Apply to the kinetic energy

Each atom has 3 translational degrees of freedom. The kinetic energy contribution from each velocity component is $\frac{1}{2}mv_x^2$, so:

$$\left\langle \frac{1}{2}mv_{ix}^2 \right\rangle = \frac{1}{2}k_BT$$

Step 3: Sum over all atoms and components

For N atoms, the total average kinetic energy sums over all 3N velocity components:

$$\langle KE \rangle = \sum_{i=1}^{N} \frac{1}{2}m_i\langle v_i^2 \rangle = \frac{3N}{2}k_BT$$

Step 4: Define the instantaneous temperature

In MD, the instantaneous temperature is computed from the kinetic energy at each timestep, accounting for constraints that remove degrees of freedom:

$$T = \frac{2\langle KE \rangle}{N_f k_B} = \frac{1}{N_f k_B}\sum_{i=1}^{N} m_i v_i^2$$

where $N_f = 3N - N_c$ is the number of degrees of freedom after removing $N_c$ constraints (e.g., 3 for center-of-mass translation, 3 for rotation, plus any bond constraints).

2. Molecular Mechanics Force Fields

The potential energy function (force field) is the sum of bonded and non-bonded interactions:

$$U = U_{\text{bonded}} + U_{\text{non-bonded}}$$

Bonded Terms

$$U_{\text{bonds}} = \sum_{\text{bonds}} k_b(r - r_0)^2$$

$$U_{\text{angles}} = \sum_{\text{angles}} k_\theta(\theta - \theta_0)^2$$

$$U_{\text{dihedrals}} = \sum_{\text{dihedrals}} k_\phi[1 + \cos(n\phi - \delta)]$$

Non-bonded Terms

$$U_{\text{LJ}} = \sum_{i<j} 4\varepsilon_{ij}\left[\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12} - \left(\frac{\sigma_{ij}}{r_{ij}}\right)^6\right]$$

$$U_{\text{Coulomb}} = \sum_{i<j} \frac{q_i q_j}{4\pi\varepsilon_0 r_{ij}}$$

Derivation: Lennard-Jones Potential — Equilibrium Distance and Well Depth

The Lennard-Jones 12-6 potential models van der Waals interactions. We derive its key features from the functional form.

Step 1: Write the LJ potential

The 12-6 form combines short-range Pauli repulsion ($r^{-12}$) with long-range London dispersion attraction ($r^{-6}$):

$$V(r) = 4\varepsilon\left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$$

where $\varepsilon$ is the depth of the potential well and $\sigma$ is the distance at which $V = 0$.

Step 2: Find the equilibrium distance (potential minimum)

Set $dV/dr = 0$:

$$\frac{dV}{dr} = 4\varepsilon\left[-12\frac{\sigma^{12}}{r^{13}} + 6\frac{\sigma^{6}}{r^{7}}\right] = 0$$

$$12\sigma^{12}/r^{13} = 6\sigma^6/r^7 \implies r^6 = 2\sigma^6 \implies r_{\min} = 2^{1/6}\sigma \approx 1.122\sigma$$

Step 3: Verify the well depth

Evaluating $V(r_{\min})$ at $r = 2^{1/6}\sigma$:

$$V(r_{\min}) = 4\varepsilon\left[\frac{1}{4} - \frac{1}{2}\right] = 4\varepsilon\left(-\frac{1}{4}\right) = -\varepsilon$$

Confirming that $\varepsilon$ is indeed the well depth.

Step 4: Derive the force from the potential

The force between two particles is:

$$F(r) = -\frac{dV}{dr} = \frac{4\varepsilon}{r}\left[12\left(\frac{\sigma}{r}\right)^{12} - 6\left(\frac{\sigma}{r}\right)^{6}\right] = \frac{24\varepsilon}{r}\left[2\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6}\right]$$

Repulsive for $r < r_{\min}$ and attractive for $r > r_{\min}$. For argon: $\sigma = 3.4$ A, $\varepsilon/k_B = 120$ K; for carbon: $\sigma = 3.55$ A, $\varepsilon/k_B = 28$ K.

AMBER/CHARMM

All-atom protein force fields with extensive parameterization for biomolecular simulations

OPLS-AA

Optimized for liquid-state thermodynamic properties and protein-ligand binding

MARTINI/SIRAH

Coarse-grained models enabling microsecond timescales for large biomolecular assemblies

3. Enhanced Sampling Methods

Standard MD is limited by the timescale problem — many biologically relevant processes occur on microsecond to second timescales, far beyond typical MD reach. Enhanced sampling methods overcome energy barriers to accelerate conformational exploration.

Umbrella Sampling

Applies biasing potentials along a reaction coordinate to sample high-energy regions. The free energy (PMF) is recovered via the Weighted Histogram Analysis Method (WHAM):

$$w_i(\xi) = \frac{1}{2}k_i(\xi - \xi_i^0)^2$$

$$G(\xi) = -k_BT \ln P_{\text{unbiased}}(\xi)$$

Metadynamics

Deposits Gaussian hills along collective variables to discourage revisiting sampled states. The bias potential builds up to fill free energy basins:

$$V_{\text{bias}}(\mathbf{s}, t) = \sum_{t'<t} w \exp\left(-\frac{|\mathbf{s} - \mathbf{s}(t')|^2}{2\sigma^2}\right)$$

Well-tempered metadynamics ensures convergence

Replica Exchange (REMD)

Runs multiple replicas at different temperatures. Periodic exchange attempts between adjacent temperatures follow the Metropolis criterion:

$$P_{\text{accept}} = \min\left(1, e^{(\beta_i - \beta_j)(U_i - U_j)}\right)$$

Accelerated MD (aMD)

Flattens the energy landscape by adding a boost potential when the potential falls below a threshold:

$$\Delta V(\mathbf{r}) = \frac{(E - V(\mathbf{r}))^2}{\alpha + (E - V(\mathbf{r}))}$$

GaMD provides Gaussian-shaped boost for reweighting

Derivation: Langevin Dynamics from the Fluctuation-Dissipation Theorem

Langevin dynamics adds stochastic and frictional forces to Newton's equations, implicitly modeling solvent effects. We derive the noise amplitude from the fluctuation-dissipation theorem.

Step 1: Write the Langevin equation of motion

Add friction (proportional to velocity) and a random force to Newton's equation:

$$m\ddot{x} = F(x) - \gamma\dot{x} + \eta(t)$$

where $\gamma$ is the friction coefficient and $\eta(t)$ is a stochastic (random) force representing thermal kicks from solvent molecules.

Step 2: Specify the statistical properties of the noise

The random force is Gaussian white noise with zero mean and delta-correlated in time:

$$\langle \eta(t) \rangle = 0, \quad \langle \eta(t)\eta(t') \rangle = 2A\,\delta(t - t')$$

where A is the noise amplitude we need to determine.

Step 3: Apply the fluctuation-dissipation theorem

At thermal equilibrium, the energy dissipated by friction must equal the energy input from the random force. The equipartition theorem requires $\langle \frac{1}{2}m\dot{x}^2 \rangle = \frac{1}{2}k_BT$. For a free particle (F = 0), solving the Langevin equation in the long-time limit gives:

$$\langle \dot{x}^2 \rangle = \frac{A}{m\gamma}$$

Step 4: Determine the noise amplitude

Equating with equipartition: $\frac{A}{m\gamma} = \frac{k_BT}{m}$, giving the Einstein relation:

$$A = \gamma k_BT \implies \langle \eta(t)\eta(t') \rangle = 2\gamma k_BT\,\delta(t - t')$$

Step 5: Write the final Langevin equation

Substituting $\eta(t) = \sqrt{2\gamma k_BT}\,\xi(t)$ where $\xi(t)$ is unit white noise:

$$m\ddot{x} = F(x) - \gamma\dot{x} + \sqrt{2\gamma k_BT}\,\xi(t)$$

This is the standard form used in biomolecular simulations. The friction coefficient $\gamma$ controls the coupling to the implicit solvent: low $\gamma$ recovers Newtonian dynamics; high $\gamma$ gives Brownian (overdamped) dynamics where inertia is negligible.

4. QM/MM Methods & Free Energy Calculations

Quantum Mechanics/Molecular Mechanics (QM/MM) methods treat a small chemically important region quantum mechanically while the remainder uses classical force fields. The total energy is:

$$E_{\text{total}} = E_{\text{QM}} + E_{\text{MM}} + E_{\text{QM/MM}}$$

Free energy perturbation (FEP) calculates free energy differences between states by alchemically transforming one system into another:

$$\Delta G_{A \to B} = -k_BT \ln \left\langle \exp\left(-\frac{U_B - U_A}{k_BT}\right)\right\rangle_A$$

Thermodynamic integration (TI) provides an alternative by integrating over a coupling parameter$\lambda$ that smoothly transforms the Hamiltonian:

$$\Delta G = \int_0^1 \left\langle \frac{\partial U(\lambda)}{\partial \lambda}\right\rangle_\lambda d\lambda$$

Derivation: Free Energy Perturbation (Zwanzig Equation)

We derive the FEP formula from the fundamental thermodynamic identity relating free energy to the partition function.

Step 1: Express free energies in terms of partition functions

The Helmholtz free energy of states A and B is related to their configurational partition functions:

$$G_A = -k_BT\ln Z_A, \quad G_B = -k_BT\ln Z_B$$

where $Z_A = \int e^{-U_A(\mathbf{r})/k_BT}d\mathbf{r}$ integrates over all configurations.

Step 2: Write the free energy difference

$$\Delta G = G_B - G_A = -k_BT\ln\frac{Z_B}{Z_A}$$

Step 3: Express the ratio as an ensemble average over state A

Multiply and divide $Z_B$ by $e^{-U_A/k_BT}$:

$$\frac{Z_B}{Z_A} = \frac{\int e^{-U_B/k_BT}d\mathbf{r}}{Z_A} = \frac{\int e^{-(U_B - U_A)/k_BT}\cdot e^{-U_A/k_BT}d\mathbf{r}}{Z_A}$$

The right side is the canonical ensemble average over state A of $e^{-\Delta U/k_BT}$:

$$\frac{Z_B}{Z_A} = \left\langle e^{-(U_B - U_A)/k_BT}\right\rangle_A$$

Step 4: Arrive at the Zwanzig (FEP) equation

Substituting back:

$$\Delta G_{A\to B} = -k_BT\ln\left\langle \exp\left(-\frac{U_B - U_A}{k_BT}\right)\right\rangle_A$$

Step 5: Practical considerations — convergence and staging

FEP converges poorly when states A and B differ significantly (the exponential average is dominated by rare low-energy configurations). In practice, the transformation is broken into small steps using $\lambda$-windows: $U(\lambda) = (1-\lambda)U_A + \lambda U_B$, and the total $\Delta G$ is accumulated over many small perturbation steps. This is the basis of alchemical free energy calculations used in drug design (binding affinity prediction).

Python Simulation: 2D Lennard-Jones MD

This simulation performs a 2D molecular dynamics simulation of Lennard-Jones particles using the velocity Verlet integrator with periodic boundary conditions, a Berendsen thermostat, and computes the radial distribution function g(r) to characterize the liquid structure.

2D Lennard-Jones Molecular Dynamics

Python

Velocity Verlet MD with periodic boundaries, Berendsen thermostat, and radial distribution function

script.py177 lines

import numpy as np

# 2D Lennard-Jones Molecular Dynamics with Velocity Verlet
# ==========================================================

print("=" * 65)
print("2D LENNARD-JONES MOLECULAR DYNAMICS SIMULATION")
print("Velocity Verlet Integrator")
print("=" * 65)

# Parameters
N = 16              # Number of particles (4x4 grid)
sigma = 1.0         # LJ sigma (reduced units)
epsilon = 1.0       # LJ epsilon (reduced units)
mass = 1.0          # Particle mass (reduced units)
dt = 0.002          # Time step (reduced units)
n_steps = 3000      # Total steps
T_target = 0.8      # Target temperature (reduced units)
L = 6.0             # Box size (reduced units)
r_cut = 2.5 * sigma # Cutoff radius

print(f"\nParticles: {N}")
print(f"Box size: {L} x {L} (reduced units)")
print(f"Temperature target: {T_target}")
print(f"Time step: {dt}")
print(f"Total steps: {n_steps}")
print(f"Cutoff: {r_cut}")

# Initialize positions on a grid
n_side = int(np.ceil(np.sqrt(N)))
spacing = L / n_side
positions = np.zeros((N, 2))
idx = 0
for i in range(n_side):
    for j in range(n_side):
        if idx < N:
            positions[idx] = [(i + 0.5) * spacing, (j + 0.5) * spacing]
            idx += 1

# Initialize velocities (Maxwell-Boltzmann)
np.random.seed(42)
velocities = np.random.randn(N, 2) * np.sqrt(T_target / mass)
# Remove center of mass velocity
velocities -= np.mean(velocities, axis=0)

def compute_forces(pos, L, sigma, epsilon, r_cut):
    """Compute LJ forces with minimum image convention."""
    N = len(pos)
    forces = np.zeros_like(pos)
    potential = 0.0
    virial = 0.0

for i in range(N):
        for j in range(i + 1, N):
            dr = pos[j] - pos[i]
            # Minimum image convention
            dr -= L * np.round(dr / L)
            r2 = np.sum(dr**2)
            r = np.sqrt(r2)

if r < r_cut and r > 0.5:
                sr6 = (sigma / r)**6
                sr12 = sr6 * sr6

# Force magnitude: F = 4*eps*(12*sig^12/r^13 - 6*sig^6/r^7)
                f_mag = 4.0 * epsilon * (12.0 * sr12 - 6.0 * sr6) / r2
                f_vec = f_mag * dr

forces[i] -= f_vec
                forces[j] += f_vec

# Potential (shifted at cutoff)
                V_cut = 4.0 * epsilon * ((sigma/r_cut)**12 - (sigma/r_cut)**6)
                potential += 4.0 * epsilon * (sr12 - sr6) - V_cut

# Virial for pressure
                virial += np.sum(f_vec * dr)

return forces, potential, virial

def kinetic_energy(vel, mass):
    """Compute kinetic energy."""
    return 0.5 * mass * np.sum(vel**2)

def temperature(KE, N):
    """Compute temperature from kinetic energy (2D)."""
    return 2.0 * KE / (2.0 * N - 2.0)  # 2D, remove 2 COM dof

# Initial force computation
forces, PE, virial = compute_forces(positions, L, sigma, epsilon, r_cut)
KE = kinetic_energy(velocities, mass)
T = temperature(KE, N)

print("\n--- Simulation Progress ---")
print(f"{'Step':>6} {'Time':>8} {'KE':>10} {'PE':>10} {'Total E':>10} {'Temp':>8} {'Pressure':>10}")
print("-" * 72)

energies = []
temperatures = []

for step in range(n_steps + 1):
    # Record
    total_E = KE + PE
    pressure = N * T / (L * L) + virial / (2.0 * L * L)

if step % 300 == 0:
        print(f"{step:6d} {step*dt:8.3f} {KE:10.4f} {PE:10.4f} {total_E:10.4f} {T:8.4f} {pressure:10.4f}")
        energies.append(total_E)
        temperatures.append(T)

if step == n_steps:
        break

# Velocity Verlet step 1: update positions
    positions += velocities * dt + 0.5 * forces / mass * dt**2

# Apply periodic boundary conditions
    positions = positions % L

# Compute new forces
    old_forces = forces.copy()
    forces, PE, virial = compute_forces(positions, L, sigma, epsilon, r_cut)

# Velocity Verlet step 2: update velocities
    velocities += 0.5 * (old_forces + forces) / mass * dt

# Berendsen thermostat (weak coupling)
    KE = kinetic_energy(velocities, mass)
    T = temperature(KE, N)
    if T > 0.01:
        tau = 0.5  # coupling constant
        scale = np.sqrt(1 + dt / tau * (T_target / T - 1))
        velocities *= scale
        KE = kinetic_energy(velocities, mass)
        T = temperature(KE, N)

# Analysis
print("\n--- Simulation Summary ---")
print(f"Final temperature: {T:.4f} (target: {T_target})")
print(f"Final total energy: {KE + PE:.4f}")
print(f"Energy conservation (std): {np.std(energies):.6f}")
print(f"Mean temperature: {np.mean(temperatures):.4f}")

# Radial distribution function
print("\n--- Radial Distribution Function g(r) ---")
n_bins = 20
r_max = L / 2.0
dr_bin = r_max / n_bins
g_r = np.zeros(n_bins)

for i in range(N):
    for j in range(i + 1, N):
        dr_vec = positions[j] - positions[i]
        dr_vec -= L * np.round(dr_vec / L)
        r = np.sqrt(np.sum(dr_vec**2))
        bin_idx = int(r / dr_bin)
        if bin_idx < n_bins:
            g_r[bin_idx] += 2

# Normalize
rho = N / (L * L)
for b in range(n_bins):
    r_inner = b * dr_bin
    r_outer = (b + 1) * dr_bin
    shell_area = np.pi * (r_outer**2 - r_inner**2)
    if shell_area > 0:
        g_r[b] /= (N * rho * shell_area)

print(f"{'r':>8} {'g(r)':>10}")
print("-" * 22)
for b in range(n_bins):
    r_mid = (b + 0.5) * dr_bin
    bar = '#' * int(g_r[b] * 10)
    print(f"{r_mid:8.3f} {g_r[b]:10.4f} {bar}")

print("\n[Simulation complete]")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Computation: LJ Force Calculator

This Fortran program computes all pairwise Lennard-Jones forces for a system of N particles using minimum image convention. It verifies Newton's third law and calculates thermodynamic properties including the virial pressure.

Lennard-Jones Force Computation

Fortran

Pairwise force calculation with minimum image convention and thermodynamic analysis

lj_forces.f90161 lines

program lj_force_computation
  ! Lennard-Jones Pair Potential and Force Computation for N Particles
  ! Computes all pairwise interactions with minimum image convention
  implicit none

integer, parameter :: dp = selected_real_kind(15, 307)
  integer, parameter :: N = 12      ! Number of particles
  real(dp), parameter :: sigma = 3.405_dp     ! Angstroms (Argon)
  real(dp), parameter :: epsilon = 0.01032_dp ! eV (Argon)
  real(dp), parameter :: L = 20.0_dp          ! Box size (Angstroms)
  real(dp), parameter :: r_cut = 8.5_dp       ! Cutoff (2.5*sigma)
  real(dp), parameter :: pi = 4.0_dp * atan(1.0_dp)

real(dp) :: x(N), y(N), z(N)
  real(dp) :: fx(N), fy(N), fz(N)
  real(dp) :: dx, dy, dz, r2, r, sr6, sr12
  real(dp) :: f_mag, V_total, V_pair, V_cut
  real(dp) :: virial, pressure, rho
  integer :: i, j, n_pairs

write(*,'(A)') '========================================================'
  write(*,'(A)') ' LENNARD-JONES FORCE COMPUTATION FOR N PARTICLES'
  write(*,'(A)') '========================================================'
  write(*,'(A)') ''
  write(*,'(A,I4)')       ' Number of particles: ', N
  write(*,'(A,F8.4,A)')   ' Sigma:    ', sigma, ' A'
  write(*,'(A,F10.6,A)')  ' Epsilon:  ', epsilon, ' eV'
  write(*,'(A,F8.2,A)')   ' Box size: ', L, ' A'
  write(*,'(A,F8.2,A)')   ' Cutoff:   ', r_cut, ' A'
  write(*,'(A)') ''

! V(r) = 4*eps*[(sig/r)^12 - (sig/r)^6]
  ! F(r) = 4*eps*[12*sig^12/r^13 - 6*sig^6/r^7]

! Initialize positions (simple cubic-like arrangement)
  call initialize_positions(x, y, z, N, L)

write(*,'(A)') '--- Initial Particle Positions ---'
  write(*,'(A6,A12,A12,A12)') 'Atom', 'x (A)', 'y (A)', 'z (A)'
  write(*,'(A)') '------------------------------------------'
  do i = 1, N
    write(*,'(I6,F12.4,F12.4,F12.4)') i, x(i), y(i), z(i)
  end do

! Compute forces and energy
  fx = 0.0_dp
  fy = 0.0_dp
  fz = 0.0_dp
  V_total = 0.0_dp
  virial = 0.0_dp
  n_pairs = 0

! Cutoff correction
  V_cut = 4.0_dp * epsilon * ((sigma/r_cut)**12 - (sigma/r_cut)**6)

write(*,'(A)') ''
  write(*,'(A)') '--- Pairwise Interactions ---'
  write(*,'(A4,A4,A10,A14,A14,A14)') &
    'i', 'j', 'r (A)', 'V(r) (eV)', 'F(r) (eV/A)', '|F_vec|'
  write(*,'(A)') '------------------------------------------------------------'

do i = 1, N
    do j = i + 1, N
      ! Minimum image convention
      dx = x(j) - x(i)
      dy = y(j) - y(i)
      dz = z(j) - z(i)
      dx = dx - L * nint(dx / L)
      dy = dy - L * nint(dy / L)
      dz = dz - L * nint(dz / L)

r2 = dx*dx + dy*dy + dz*dz
      r = sqrt(r2)

if (r < r_cut .and. r > 0.5_dp) then
        sr6 = (sigma / r)**6
        sr12 = sr6 * sr6

V_pair = 4.0_dp * epsilon * (sr12 - sr6) - V_cut
        f_mag = 4.0_dp * epsilon * (12.0_dp * sr12 - 6.0_dp * sr6) / r2

! Accumulate forces
        fx(i) = fx(i) - f_mag * dx
        fy(i) = fy(i) - f_mag * dy
        fz(i) = fz(i) - f_mag * dz
        fx(j) = fx(j) + f_mag * dx
        fy(j) = fy(j) + f_mag * dy
        fz(j) = fz(j) + f_mag * dz

V_total = V_total + V_pair
        virial = virial + f_mag * r2

n_pairs = n_pairs + 1
        if (n_pairs <= 30) then
          write(*,'(I4,I4,F10.4,F14.8,F14.8,F14.8)') &
            i, j, r, V_pair, f_mag * r, sqrt(f_mag*f_mag * r2)
        end if
      end if
    end do
  end do

if (n_pairs > 30) then
    write(*,'(A,I4,A)') ' ... (', n_pairs - 30, ' more pairs not shown)'
  end if

write(*,'(A)') ''
  write(*,'(A)') '--- Force Summary ---'
  write(*,'(A6,A14,A14,A14,A14)') &
    'Atom', 'Fx (eV/A)', 'Fy (eV/A)', 'Fz (eV/A)', '|F| (eV/A)'
  write(*,'(A)') '--------------------------------------------------------------'
  do i = 1, N
    write(*,'(I6,F14.8,F14.8,F14.8,F14.8)') &
      i, fx(i), fy(i), fz(i), sqrt(fx(i)**2 + fy(i)**2 + fz(i)**2)
  end do

! Verify Newton's third law
  write(*,'(A)') ''
  write(*,'(A,ES12.4)') ' Sum of Fx: ', sum(fx)
  write(*,'(A,ES12.4)') ' Sum of Fy: ', sum(fy)
  write(*,'(A,ES12.4)') ' Sum of Fz: ', sum(fz)
  write(*,'(A)') ' (Should be ~0 by Newton''s 3rd law)'

write(*,'(A)') ''
  write(*,'(A)') '--- Thermodynamic Properties ---'
  write(*,'(A,I6)')          ' Total pairs:        ', n_pairs
  write(*,'(A,F14.8,A)')     ' Total PE:           ', V_total, ' eV'
  write(*,'(A,F14.8,A)')     ' PE per particle:    ', V_total / N, ' eV'
  rho = real(N, dp) / (L**3)
  write(*,'(A,ES14.6,A)')    ' Number density:     ', rho, ' A^-3'
  write(*,'(A,F14.8,A)')     ' Virial:             ', virial, ' eV'

write(*,'(A)') ''
  write(*,'(A)') '[Computation complete]'

contains

subroutine initialize_positions(x, y, z, N, L)
    integer, intent(in) :: N
    real(dp), intent(in) :: L
    real(dp), intent(out) :: x(N), y(N), z(N)
    integer :: i, nx, ny, nz, idx
    real(dp) :: spacing

nx = ceiling(real(N, dp)**(1.0_dp/3.0_dp))
    spacing = L / real(nx, dp)

idx = 0
    do i = 1, nx
      do ny = 1, nx
        do nz = 1, nx
          idx = idx + 1
          if (idx > N) return
          x(idx) = (real(i, dp) - 0.5_dp) * spacing
          y(idx) = (real(ny, dp) - 0.5_dp) * spacing
          z(idx) = (real(nz, dp) - 0.5_dp) * spacing
        end do
      end do
    end do
  end subroutine

end program lj_force_computation

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

Video Lectures

Introduction to Molecular Dynamics

Comprehensive introduction to molecular dynamics simulation methodology, including force fields, integrators, thermostats, and barostats.

Enhanced Sampling Methods

Overview of enhanced sampling techniques including metadynamics, replica exchange, umbrella sampling, and their applications in biomolecular simulations.

Key Concepts Summary

Velocity Verlet

Symplectic, time-reversible integrator with excellent long-term energy conservation for MD

Force Fields

AMBER, CHARMM, OPLS parameterize bonded and non-bonded interactions for biomolecular systems

Enhanced Sampling

Metadynamics, REMD, and umbrella sampling overcome energy barriers to access rare events

QM/MM & Free Energy

Hybrid quantum/classical methods and FEP/TI for accurate thermodynamic calculations

← Experimental Techniques Back to Molecular Biology →