Simulation Protocol¶

AdsPro uses a three-phase simulation strategy: a fast rigid-body screen (Phase 0), full Langevin MD relaxation (Phase 1), and free energy calculation via umbrella sampling (Phase 2, described in PMF & WHAM).

Langevin Dynamics¶

All MD in AdsPro uses Langevin dynamics at constant temperature T (NVT ensemble). The equation of motion for bead i with mass m_i is:

\[ m_i \ddot{\mathbf{r}}_i = \mathbf{F}_i^\text{ff} - \gamma m_i \dot{\mathbf{r}}_i + \boldsymbol{\xi}_i(t) \]

where: - \(\mathbf{F}_i^\text{ff}\) — deterministic force from the force field (Morse, DH, Gō, backbone) - \(-\gamma m_i \dot{\mathbf{r}}_i\) — Langevin damping (friction), \(\gamma\) = friction coefficient per unit mass - \(\boldsymbol{\xi}_i(t)\) — Gaussian random force (thermal noise)

The random force satisfies the fluctuation-dissipation theorem:

\[ \langle \xi_{i,\alpha}(t) \, \xi_{j,\beta}(t') \rangle = 2 \gamma m_i k_B T \, \delta_{ij} \, \delta_{\alpha\beta} \, \delta(t - t') \]

This guarantees that the simulation samples the correct canonical ensemble at temperature T, regardless of the damping coefficient γ.

The random force amplitude at each step is:

\[ |\boldsymbol{\xi}_i| \sim \sqrt{\frac{2 \gamma m_i k_B T}{\Delta t}} \cdot \mathcal{N}(0,1) \]

Damping Coefficient¶

The friction coefficient γ = 1 ps⁻¹ (mass-independent, applied uniformly). This value is appropriate for an implicit-solvent CG model: large enough to maintain temperature, small enough that the dynamics are not overdamped (Brownian limit).

Velocity-Verlet Integration¶

Each MD step uses a velocity-Verlet algorithm with Langevin forces applied as additional force terms:

Step 1 — Half-step velocity update with forces at time t: [ \mathbf{v}_i!\left(t + \tfrac{\Delta t}{2}\right) = \mathbf{v}_i(t) + \frac{\Delta t}{2} \frac{\mathbf{F}_i(t)}{m_i} ]

Step 2 — Full position update: [ \mathbf{r}_i(t + \Delta t) = \mathbf{r}_i(t) + \Delta t \, \mathbf{v}_i!\left(t + \tfrac{\Delta t}{2}\right) ]

Step 3 — Force re-evaluation at \(t + \Delta t\), including new Langevin noise and damping: [ \mathbf{F}_i(t+\Delta t) = \mathbf{F}_i^\text{ff}(t+\Delta t) - \gamma m_i \mathbf{v}_i!\left(t+\tfrac{\Delta t}{2}\right) + \boldsymbol{\xi}_i(t+\Delta t) ]

Step 4 — Second half-step velocity update: [ \mathbf{v}_i(t + \Delta t) = \mathbf{v}_i!\left(t + \tfrac{\Delta t}{2}\right) + \frac{\Delta t}{2} \frac{\mathbf{F}_i(t + \Delta t)}{m_i} ]

This scheme is time-reversible, symplectic for the conservative part, and numerically stable for time steps up to 0.02 ps in typical CG models.

Time Step¶

The default time step is Δt = 0.01 ps (10 fs). This resolves the fastest motions in the Gō model (backbone stretching, ~0.1 ps period) with at least 10 steps per oscillation. The maximum recommended Δt is 0.02 ps.

Phase 0: Rigid-Body Orientation Screening¶

Before expensive MD, AdsPro rapidly evaluates 100 protein orientations using rigid-body docking — no relaxation.

Orientation Generation¶

Orientations are generated by sampling the rotation group SO(3) uniformly. This requires three independent angles (Euler angles), but uniform sampling in Euler angles does not give uniform coverage of SO(3). AdsPro uses a factored scheme: sample azimuthal angle φ uniformly in [0, 2π), cosine of polar angle θ uniformly in [−1, 1], and rotation around the protein axis ψ uniformly in [0, 2π).

Scoring¶

For each rigid orientation, the protein is placed at placement_distance from the NP surface and the interaction energy is computed:

\[ E_\text{rigid}(\Omega) = E_\text{Morse}(\Omega) + E_\text{DH}(\Omega) \]

The top n_orientations by this score are passed to Phase 1.

Protected Orientations¶

Three orientations are always included regardless of Phase 0 ranking, to ensure important binding modes are never missed:

Affinity-patch orientation — the 15–20 residues with the most negative Morse ε are rotated to face the surface. This ensures the highest-affinity face is always explored.
Charge-complementary orientation — positively charged residues (ARG, LYS, HIS) are rotated to face a negatively charged surface (ζ < 0), or negative residues toward a positive surface (ζ > 0). This captures the electrostatically favored orientation.
Dipole orientation — the protein electric dipole moment vector μ is aligned pointing toward the surface center. For strongly polar proteins with |μ| > 50 Debye, this is often the lowest-energy orientation.

Phase 1: Langevin MD Relaxation¶

Each orientation selected from Phase 0 undergoes full Langevin MD to allow the protein to: - Translate toward or away from the NP - Rotate to optimize contacts - Internally deform (stretch native contacts, form new surface bonds)

Energy Minimization¶

Before MD, each orientation undergoes a short energy minimization using the steepest-descent algorithm:

\[ \mathbf{r}_i^{(n+1)} = \mathbf{r}_i^{(n)} - h \cdot \nabla_i E \]

with step size h adapted to give a maximum force below em_tol (default 500 kJ/mol/nm) within em_max_steps (default 5000). This removes any initial clashes.

COM Restraint¶

A harmonic COM distance restraint prevents the protein from drifting away from the surface during Phase 1:

\[ V_\text{COM}(r_\text{COM}) = k_\text{COM} \, (r_\text{COM} - r_\text{target})^2 \quad \text{if } r_\text{COM} > r_\text{max} \]

with \(k_\text{COM} = 5000\) kJ/(mol·nm²). This is a one-sided harmonic: the protein is free to move closer to the surface (improving contacts) but not to drift far away. com_restraint_max_distance: auto sets the ceiling to the initial COM distance for each orientation.

Convergence Detection¶

Phase 1 terminates when the total energy change over the last equilibrium_window steps is below equilibrium_threshold (default 0.01 kJ/mol). This indicates the protein has settled into a locally stable adsorbed conformation.

If convergence is not reached within max_steps (default 100,000), the simulation stops anyway and uses the final conformation.

Best Orientation Selection¶

After Phase 1, all orientations are ranked by their final Morse energy (most negative = strongest adsorption). The best orientation is selected as the starting point for Phase 2.

The orientation summary table printed during the run shows:

  Orient    E_morse     ΔSASA    ΔE_Gō    ΔE_bonds   Contacts   Tilt    Converged
    1        -48.25     +13.01   +676.96   +551.15       28      39.7°     No
    8        -91.84     +16.96   +677.81   +565.67       31      51.1°     No  ← BEST

Parallel Phase 0¶

Phase 0 orientation scoring is parallelized across n_workers CPU cores. With n_workers=4, 100 orientations are evaluated in parallel batches, reducing Phase 0 from ~400 sequential evaluations to ~100 wall-clock time units.

Phase 1 orientations run sequentially (not in parallel) to avoid memory contention with the surface bead arrays. Parallelizing Phase 1 is a planned future enhancement.