Bio Applied Molecular Modeling

Applied Molecular Modeling

When to Use

• Explaining or prototyping force-field energy terms (bond, angle, dihedral, van der Waals, Coulomb) before running GROMACS/OpenMM/AMBER
• Comparing energy-minimization algorithms (steepest descent vs. conjugate gradient) on a toy potential before a real minimization step
• Choosing between QM, MM/force-field, and coarse-grained levels of theory for a modeling problem
• Building or QC-ing a homology model (template search → alignment → build → refine → Ramachandran/QMEAN validation)
• Interpreting MD trajectory QC metrics (RMSD, RMSF, radius of gyration, H-bond counts) or AutoDock Vina docking scores

Version Compatibility

numpy ≥1.24, pandas ≥2.0, matplotlib ≥3.7, rdkit ≥2023.09. Concepts apply to GROMACS ≥2023.x, AutoDock Vina ≥1.2.x, AMBER/CHARMM force fields — this skill covers the underlying math and result interpretation, not running the external MD/docking engines themselves.

Prerequisites

pip install numpy pandas matplotlib rdkit. Familiarity with classical mechanics (potential energy, gradients) helps. Related structure skills: bio-structural-biology-structure-io, bio-chemoinformatics-molecular-descriptors.

Levels of theory: QM (electrons explicit, small systems) → MM/force fields (parametrized potentials, full proteins) → coarse-grained (sacrifices atomic detail for timescale). Always a tradeoff: accuracy vs. computational cost. QM (ab initio, DFT, or semi-empirical AM1/PM3/PM7) is needed for bond breaking/forming, charge distributions, or force-field parametrization.

Force Fields

$$E_{\text{total}} = E_{\text{bonds}} + E_{\text{angles}} + E_{\text{dihedrals}} + E_{\text{electrostatic}} + E_{\text{vdW}}$$

Bond: $\tfrac12 k_b (r-r_0)^2$. Angle: $\tfrac12 k_\theta(\theta-\theta_0)^2$. Dihedral: $\sum_n V_n[1+\cos(n\phi-\gamma)]$. Never mix parameters from different force fields (AMBER ff14SB, CHARMM/CGenFF, OPLS-AA, GROMOS, MARTINI coarse-grained).

Goal: compute and plot each force-field energy term for a molecule. Approach: implement the harmonic/LJ/Coulomb formulas directly with NumPy, then plot to sanity-check equilibrium positions and well depths.

import numpy as np
import matplotlib.pyplot as plt

def lennard_jones(r, epsilon=1.0, sigma=1.0):
    """Non-bonded van der Waals potential: V(r) = 4*eps*((sigma/r)^12 - (sigma/r)^6)."""
    return 4 * epsilon * ((sigma / r) ** 12 - (sigma / r) ** 6)

def coulomb(r, q1, q2, k_coul=332.0637):
    """Electrostatic energy (kcal/mol) for charges q1, q2 (e) separated by r (Angstrom)."""
    return k_coul * q1 * q2 / r

r = np.linspace(0.9, 3.0, 500)
r_min = 2 ** (1 / 6)  # equilibrium distance in units of sigma
V = lennard_jones(r)

fig, ax = plt.subplots()
ax.plot(r, V)
ax.axvline(r_min, color="r", linestyle="--", label=f"r_min = {r_min:.3f} sigma")
ax.set_ylim(-1.5, 3)
ax.set_xlabel("distance r (sigma)"); ax.set_ylabel("V(r) (epsilon)")
ax.legend(); plt.tight_layout()

Energy Minimization

Remove steric clashes before simulation. Steepest descent: robust, slow near the minimum — good initial relaxation. Conjugate gradient: uses prior steps, faster near the minimum. L-BFGS: quasi-Newton, fastest for smooth surfaces. GROMACS convention: steepest descent until max force < 1000 kJ/mol/nm, then optionally conjugate gradient.

Goal: minimize a rough energy surface and confirm convergence to the known minimum. Approach: use the Rosenbrock function as a stand-in energy surface (steep in some directions, shallow in others, like a real protein energy landscape) and run gradient descent with an explicit, correct analytic gradient.

import numpy as np

def rosenbrock(x, y, a=1, b=100):
    """Toy energy surface: steep valley, analogous to a rugged protein energy landscape."""
    return (a - x) ** 2 + b * (y - x ** 2) ** 2

def rosenbrock_grad(x, y, a=1, b=100):
    """Analytic gradient of rosenbrock(); needed for any gradient-based minimizer."""
    dx = -2 * (a - x) + b * 2 * (y - x ** 2) * (-2 * x)
    dy = b * 2 * (y - x ** 2)
    return np.array([dx, dy])

def steepest_descent(x0, y0, lr=0.001, steps=3000):
    """Steepest-descent minimization; returns the (x, y) path taken."""
    path = [(x0, y0)]
    x, y = x0, y0
    for _ in range(steps):
        grad = rosenbrock_grad(x, y)
        x -= lr * grad[0]
        y -= lr * grad[1]
        path.append((x, y))
    return np.array(path)

path = steepest_descent(-1.5, 2.0, lr=0.001, steps=3000)
final_x, final_y = path[-1]
assert rosenbrock(final_x, final_y) < rosenbrock(-1.5, 2.0)  # energy decreased
print(f"final=({final_x:.4f},{final_y:.4f}) energy={rosenbrock(final_x, final_y):.6f} (true min: (1,1)=0)")

Homology Modeling

Proteins with similar sequences adopt similar structures.

Target sequence -> Template search (BLAST/HHpred/SWISS-MODEL)
-> Target-template alignment (critical step)
-> Model building (copy coords, build loops, add side chains)
-> Refinement (energy minimization, MD relaxation)
-> Quality check (Ramachandran, QMEAN, ProSA)

A good model should have >90% of residues in Ramachandran-favored regions and <0.5% outliers.

MD Trajectory QC and Docking Scores

Goal: sanity-check an MD run (energy conservation, equilibration) and parse docking results. Approach: verify a Verlet integrator conserves energy on a known analytic system, then parse AutoDock Vina's tabular log into a DataFrame for scoring/filtering.

import numpy as np
import pandas as pd

def verlet_harmonic(k=1.0, m=1.0, x0=1.0, v0=0.0, dt=0.05, n_steps=500):
    """Velocity-Verlet integration of a 1D harmonic oscillator; used to check energy conservation."""
    x, v = np.zeros(n_steps), np.zeros(n_steps)
    x[0], v[0] = x0, v0
    a = -k * x[0] / m
    for i in range(1, n_steps):
        x[i] = x[i - 1] + v[i - 1] * dt + 0.5 * a * dt ** 2
        a_new = -k * x[i] / m
        v[i] = v[i - 1] + 0.5 * (a + a_new) * dt
        a = a_new
    return x, v

x, v = verlet_harmonic(dt=0.1, n_steps=500)
total_energy = 0.5 * v ** 2 + 0.5 * x ** 2
assert abs(total_energy[-1] - total_energy[0]) < 1e-3  # energy conserved -> integrator/timestep OK

def parse_vina_output(text):
    """Parse an AutoDock Vina log's mode table into a DataFrame (mode, affinity, RMSD bounds)."""
    rows = []
    for line in text.strip().split("\n"):
        line = line.strip()
        if line and line[0].isdigit():
            mode, affinity, rmsd_lb, rmsd_ub = line.split()
            rows.append({
                "mode": int(mode),
                "affinity_kcal_mol": float(affinity),
                "rmsd_lb": float(rmsd_lb),
                "rmsd_ub": float(rmsd_ub),
            })
    return pd.DataFrame(rows)

Pitfalls

• Never mix force-field parameter sets (e.g. AMBER charges with CHARMM bonded terms) — energies become meaningless.
• A minimized/relaxed structure is not a sampled ensemble — don't draw thermodynamic conclusions from a single minimum.
• Homology models are only as good as the template: <30% sequence identity to the template is a red flag; always check Ramachandran/QMEAN before trusting a model.
• Vina "affinity" is a scoring-function estimate, not a measured Kd/Ki — use it for ranking within one target, not absolute affinity claims.
• RMSD/RMSF only mean "equilibrated" if computed after excluding the initial relaxation window; check the running average, not raw noisy values.

See Also

• bio-structural-biology-structure-io
• bio-structural-biology-geometric-analysis
• bio-chemoinformatics-molecular-descriptors
• molecular-dynamics