# Bio Applied Molecular Modeling

> Compute force-field energy terms (bond/LJ/Coulomb), run energy minimization, and QC MD/homology models (RMSD, RMSF, Ramachandran) in NumPy. Use for force-field, minimization, or homology/docking validation.

- Skill: `pavel-kravchenko/bio-applied-molecular-modeling` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-molecular-modeling`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-molecular-modeling/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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.

```python
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.

```python
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.

```text
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.

```python
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`