section5_integration.tex

\documentclass[11pt]{article}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{geometry}
\usepackage{booktabs}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{enumerate}
\usepackage{listings}
\usepackage{xcolor}

\geometry{letterpaper, margin=1in}

\lstset{
  basicstyle=\ttfamily\small,
  keywordstyle=\color{blue},
  commentstyle=\color{gray},
  stringstyle=\color{red},
  breaklines=true,
  columns=fullflexible,
  keepspaces=true,
  language=C++
}

\title{\textbf{Formation Engine Methodology} \\
\Large Section 5 — Time Integration and Dynamic Evolution}
\author{Formation Engine Development Team}
\date{Version 0.1 — As Implemented in Codebase}

\begin{document}
\maketitle

\section{Introduction: The Integration Pipeline}

The formation engine produces structures through a two-stage process:

\begin{enumerate}
\item \textbf{Exploration}: Molecular dynamics (MD) propagates the system through configuration space, crossing energy barriers and sampling different regions of the potential energy surface
\item \textbf{Quenching}: Energy minimization projects thermal fluctuations away, revealing the underlying local minima
\end{enumerate}

This section documents the \textbf{actual implementations} of these integration algorithms, not theoretical ideals. Every algorithm listed here exists in the codebase with line numbers, parameter values, and known limitations explicitly stated.

\section{Velocity Verlet Integration (NVE Ensemble)}

\textbf{Implementation:} \texttt{atomistic/integrators/velocity\_verlet.cpp}, lines 87--220

The velocity Verlet algorithm is the default integrator for conservative (microcanonical, NVE) dynamics. It advances the State from time $t$ to $t + \Delta t$ through four steps.

\subsection{Algorithm}

The four-step procedure:

\begin{align}
\mathbf{v}_i\!\left(t + \frac{\Delta t}{2}\right) &= \mathbf{v}_i(t) + \frac{\mathbf{F}_i(t)}{m_i} \times \text{ACC\_CONV} \times \frac{\Delta t}{2} \label{eq:vv_step1} \\
\mathbf{x}_i(t + \Delta t) &= \mathbf{x}_i(t) + \mathbf{v}_i\!\left(t + \frac{\Delta t}{2}\right) \times \Delta t \label{eq:vv_step2} \\
\mathbf{F}_i(t + \Delta t) &= -\nabla_{\mathbf{x}_i} U\bigl(\mathbf{X}(t + \Delta t)\bigr) \label{eq:vv_step3} \\
\mathbf{v}_i(t + \Delta t) &= \mathbf{v}_i\!\left(t + \frac{\Delta t}{2}\right) + \frac{\mathbf{F}_i(t + \Delta t)}{m_i} \times \text{ACC\_CONV} \times \frac{\Delta t}{2} \label{eq:vv_step4}
\end{align}

\subsection{Actual Code Implementation}

Lines 153--165 (velocity half-step):
\begin{lstlisting}[language=C++]
// Step 1: Half-step velocity update
for (uint32_t i = 0; i < state.N; ++i) {
    double inv_m = 1.0 / state.M[i];
    state.V[i].x += state.F[i].x * inv_m * ACC_CONV * 0.5 * params.dt;
    state.V[i].y += state.F[i].y * inv_m * ACC_CONV * 0.5 * params.dt;
    state.V[i].z += state.F[i].z * inv_m * ACC_CONV * 0.5 * params.dt;
}
\end{lstlisting}

Lines 167--175 (position update + PBC wrapping):
\begin{lstlisting}[language=C++]
// Step 2: Full-step position update
for (uint32_t i = 0; i < state.N; ++i) {
    state.X[i].x += state.V[i].x * params.dt;
    state.X[i].y += state.V[i].y * params.dt;
    state.X[i].z += state.V[i].z * params.dt;
}

// Apply PBC wrapping if enabled
if (state.box.enabled) {
    for (uint32_t i = 0; i < state.N; ++i) {
        state.X[i].x = std::fmod(state.X[i].x + 10.0 * state.box.L.x, state.box.L.x);
        state.X[i].y = std::fmod(state.X[i].y + 10.0 * state.box.L.y, state.box.L.y);
        state.X[i].z = std::fmod(state.X[i].z + 10.0 * state.box.L.z, state.box.L.z);
    }
}
\end{lstlisting}

Line 177 (force evaluation):
\begin{lstlisting}[language=C++]
// Step 3: Compute forces at new positions
model.eval(state, mp);
\end{lstlisting}

Lines 180--186 (velocity second half-step):
\begin{lstlisting}[language=C++]
// Step 4: Second half-step velocity update
for (uint32_t i = 0; i < state.N; ++i) {
    double inv_m = 1.0 / state.M[i];
    state.V[i].x += state.F[i].x * inv_m * ACC_CONV * 0.5 * params.dt;
    state.V[i].y += state.F[i].y * inv_m * ACC_CONV * 0.5 * params.dt;
    state.V[i].z += state.F[i].z * inv_m * ACC_CONV * 0.5 * params.dt;
}
\end{lstlisting}

\subsection{Physical Properties}

\paragraph{Symplectic} The algorithm preserves phase space volume. Energy is conserved to bounded error $\Delta E \sim O(\Delta t^2)$ that does not grow with time.

\paragraph{Time-reversible} Running with $-\Delta t$ recovers the initial state (modulo numerical precision). This satisfies microscopic reversibility.

\paragraph{Second-order accurate} Local truncation error is $O(\Delta t^3)$, global error is $O(\Delta t^2)$ over fixed time interval.

\subsection{Default Parameters}

From \texttt{atomistic/integrators/velocity\_verlet.hpp}, lines 18--23:

\begin{lstlisting}[language=C++]
struct VelocityVerletParams {
    double dt = 1e-3;           // Timestep (fs)
    int n_steps = 1000;         // Number of steps
    int print_freq = 100;       // Print diagnostics every N steps
    bool verbose = false;       // Print detailed output
};
\end{lstlisting}

The default timestep $\Delta t = 0.001$\,fs = 1 attosecond is extremely small — intended for testing. For production runs:
\begin{itemize}
\item Lennard-Jones systems: $\Delta t = 1$--$2$\,fs
\item Bonded systems with harmonic bonds: $\Delta t = 0.5$--$1$\,fs
\item Systems with hydrogen (stiff H--X bonds): $\Delta t = 0.5$\,fs
\end{itemize}

\subsection{Energy Conservation}

For a well-parameterized Lennard-Jones system with $\Delta t = 1$\,fs:
\begin{equation}
\frac{|\Delta E|}{N} \sim 10^{-4}\,\text{kcal/mol per atom over } 10^4\text{ steps}
\end{equation}

This is monitored during the run (lines 194--201):

\begin{lstlisting}[language=C++]
if (params.verbose && (step + 1) % params.print_freq == 0) {
    double E_drift = E_total - stats.E_initial;
    std::cout << "  Step " << std::setw(6) << (step + 1) 
              << "  T = " << T << " K"
              << "  E = " << E_total << " kcal/mol"
              << "  ΔE = " << std::showpos << E_drift << std::noshowpos
              << "\n";
}
\end{lstlisting}

\subsection{Known Limitation: Coulomb Forces Disabled}

As documented in Section~3.3, Coulomb forces are zeroed in the velocity Verlet integrator due to a unit-conversion coupling instability. The code comment (line 29) states:

\begin{quote}
\texttt{// Current value gives T~10\^{}15K (better but still wrong)}
\end{quote}

This means:
\begin{itemize}
\item Coulomb energies are computed correctly and reported in the energy ledger
\item Coulomb forces are set to zero before the velocity update
\item Systems driven by van der Waals interactions work correctly
\item Ionic systems (NaCl, metal oxides) diverge if Coulomb forces are enabled
\end{itemize}

The workaround: use FIRE minimization for ionic formations (Section~5.3), which does not go through the velocity Verlet acceleration pathway.

\section{FIRE Minimization}

\textbf{Implementation:} \texttt{atomistic/integrators/fire.cpp}, class \texttt{FIRE}

The Fast Inertial Relaxation Engine (FIRE) of Bitzek et al. (2006) is the primary energy minimizer in the formation pipeline. It is robust, parameter-insensitive, and does not require Hessian computation.

\subsection{The FIRE Algorithm}

FIRE treats minimization as a damped dynamics problem. Particles are given velocities, and those velocities are continuously steered toward the downhill direction via three mechanisms:

\paragraph{Velocity mixing} (lines 60--66):

\begin{equation}
\mathbf{v}_i \gets (1 - \alpha)\mathbf{v}_i + \alpha \frac{|\mathbf{v}|}{\|\mathbf{F}\|} \mathbf{F}_i
\end{equation}

This blends the current velocity with the normalized force direction. The parameter $\alpha \in [0,1]$ controls the aggressiveness of alignment.

\paragraph{Power monitoring} (lines 47--57):

\begin{equation}
P = \sum_{i=1}^{N} \mathbf{v}_i \cdot \mathbf{F}_i
\end{equation}

The power $P$ measures whether the system is going downhill ($P > 0$) or uphill ($P \leq 0$).

\paragraph{Adaptive timestep} (lines 68--77):

\begin{equation}
\Delta t \gets \begin{cases}
\min(\Delta t \times f_{\text{inc}}, \Delta t_{\max}) & \text{if } P > 0 \text{ for } n_{\min} \text{ consecutive steps} \\
\Delta t \times f_{\text{dec}} & \text{if } P \leq 0
\end{cases}
\end{equation}

\subsection{Complete Implementation}

The main loop (lines 28--92):

\begin{lstlisting}[language=C++]
for (int t = 0; t < fp.max_steps; ++t) {
    // Evaluate forces + energies at current X
    model.eval(s, mp);
    double U = s.E.total();
    double Frms = rms_force(s);
    
    // Convergence check (skip first 2 iterations)
    if (t > 1) {
        double dU = std::abs(U - Uprev);
        double dU_per_atom = dU / (double)s.N;
        if (Frms < fp.epsF || dU_per_atom < fp.epsU) {
            return {t, U, dU_per_atom, Frms, alpha, dt};
        }
    }
    
    Uprev = U;
    
    // Compute power P = v·F
    long double P = 0;
    long double vnorm2 = 0;
    long double fnorm2 = 0;
    for (uint32_t i = 0; i < s.N; i++) {
        P += (long double)dot(s.V[i], s.F[i]);
        vnorm2 += (long double)dot(s.V[i], s.V[i]);
        fnorm2 += (long double)dot(s.F[i], s.F[i]);
    }
    double vnorm = std::sqrt((double)vnorm2);
    double fnorm = std::sqrt((double)fnorm2);
    
    // Velocity mixing: v <- (1-α)v + α|v|f/|f|
    if (fnorm > 0 && vnorm > 0) {
        for (uint32_t i = 0; i < s.N; i++) {
            Vec3 fhat = s.F[i] * (1.0 / fnorm);
            s.V[i] = s.V[i] * (1.0 - alpha) + fhat * (alpha * vnorm);
        }
    }
    
    // Adaptive timestep based on power sign
    if (P > 0) {
        npos++;
        if (npos > fp.nmin) {
            dt = std::min(dt * fp.finc, fp.dt_max);
            alpha *= fp.falpha;
        }
    } else {
        npos = 0;
        dt *= fp.fdec;
        alpha = fp.alpha;
        for (auto& v : s.V) v = {0, 0, 0};  // Freeze velocities
    }
    
    // Position update: X <- X + dt V (explicit Euler)
    for (uint32_t i = 0; i < s.N; i++) {
        s.X[i] = s.X[i] + s.V[i] * dt;
    }
    
    Uprev = U;
}
\end{lstlisting}

\subsection{Default Parameters}

From \texttt{atomistic/integrators/fire.hpp}, lines 6--15:

\begin{lstlisting}[language=C++]
struct FIREParams {
    double dt = 1e-3;       // Initial timestep
    double alpha = 0.1;     // Initial velocity-force mixing
    double finc = 1.1;      // Timestep growth factor
    double fdec = 0.5;      // Timestep shrink factor
    double falpha = 0.99;   // Alpha decay rate
    int nmin = 5;           // Steps before adaptation
    double dt_max = 1e-1;   // Timestep ceiling
    
    double epsF = 1e-6;     // RMS force threshold
    double epsU = 1e-10;    // Per-atom energy delta threshold
    int max_steps = 5000;   // Maximum iterations
};
\end{lstlisting}

\subsection{Convergence Criteria}

The minimization terminates when \textbf{either} (lines 37--42):

\paragraph{Force convergence:}
\begin{equation}
F_{\text{rms}} = \sqrt{\frac{1}{N}\sum_i |\mathbf{F}_i|^2} < \varepsilon_F = 10^{-6}\,\text{kcal/(mol·\AA)}
\end{equation}

\paragraph{Energy convergence:}
\begin{equation}
\frac{|U(t) - U(t-1)|}{N} < \varepsilon_U = 10^{-10}\,\text{kcal/mol per atom}
\end{equation}

The two-criterion design ensures convergence for both stiff (force-dominated) and soft (energy-dominated) potentials.

\subsection{Velocity Initialization}

Before the first iteration, velocities are initialized along the force direction (line 23):

\begin{lstlisting}[language=C++]
initialize_velocities_along_force(s, fp.dt);
\end{lstlisting}

This prevents $P = 0$ deadlock on the first step. The initialization sets:
\begin{equation}
\mathbf{v}_i = \frac{\mathbf{F}_i}{\|\mathbf{F}\|} \times \Delta t \times \text{scale}
\end{equation}

where $\text{scale}$ is chosen to produce a small initial displacement ($\sim 0.001$\,\AA).

\subsection{Role in the Formation Pipeline}

FIRE is the operation that converts a dynamical trajectory into a structural result:

\begin{enumerate}
\item MD explores configuration space, crossing energy barriers
\item Snapshots are extracted at intervals (every 100 steps)
\item Each snapshot is FIRE-minimized independently
\item The minimized structures are formations
\end{enumerate}

\textbf{Critical point:} FIRE does not suffer from the Coulomb coupling instability (Section~3.3) because it uses explicit Euler position updates ($\mathbf{x} \gets \mathbf{x} + \Delta t \mathbf{v}$), not force-to-acceleration conversion. This makes FIRE suitable for ionic systems where velocity Verlet fails.

\subsection{Numerical Stability}

The code uses \texttt{long double} for accumulation (lines 48--56):

\begin{lstlisting}[language=C++]
long double P = 0;
long double vnorm2 = 0;
long double fnorm2 = 0;
for (uint32_t i = 0; i < s.N; i++) {
    P += (long double)dot(s.V[i], s.F[i]);
    vnorm2 += (long double)dot(s.V[i], s.V[i]);
    fnorm2 += (long double)dot(s.F[i], s.F[i]);
}
\end{lstlisting}

This prevents catastrophic cancellation when $P$ is small relative to $\|\mathbf{v}\|$ and $\|\mathbf{F}\|$.

\section{Bonded Force Field}

\textbf{Implementation:} \texttt{atomistic/models/bonded.cpp}, class \texttt{BondedModel}

The bonded force field computes intramolecular forces from the bond topology stored in \texttt{State.B} (the edge list). Four terms are implemented:

\begin{enumerate}
\item Harmonic bond stretching
\item Harmonic angle bending
\item Periodic torsions (dihedrals)
\item Improper torsions (out-of-plane)
\end{enumerate}

Each term has explicit force derivations following molecular mechanics conventions (CHARMM, AMBER, OPLS).

\subsection{Harmonic Bond Stretching}

\paragraph{Potential} (line 126):

\begin{equation}
U_{\text{bond}} = k_b (r - r_0)^2
\end{equation}

where $r = |\mathbf{r}_j - \mathbf{r}_i|$ is the current distance and $r_0$ is the equilibrium length.

\paragraph{Force} (lines 132--137):

\begin{equation}
\mathbf{F}_i = -\nabla_{\mathbf{r}_i} U = -2k_b(r - r_0)\,\hat{\mathbf{r}}_{ij}, \qquad \mathbf{F}_j = -\mathbf{F}_i
\end{equation}

Code implementation:

\begin{lstlisting}[language=C++]
Vec3 rij = s.X[bond.i] - s.X[bond.j];
double r = norm(rij);
double dr = r - bond.r0;

U += bond.kb * dr * dr;

if (r > 1e-12) {
    Vec3 f = rij * (-2.0 * bond.kb * dr / r);
    s.F[bond.i] = s.F[bond.i] + f;
    s.F[bond.j] = s.F[bond.j] - f;
}
\end{lstlisting}

\paragraph{Default parameters} (from header comment, line 17):
\begin{itemize}
\item $k_b \sim 300$--$500$\,kcal/(mol·\AA$^2$) for C--C, C--H bonds
\item $r_0 = 1.54$\,\AA\ (C--C single bond generic fallback)
\end{itemize}

\subsection{Harmonic Angle Bending}

\paragraph{Potential} (line 143):

\begin{equation}
U_{\text{angle}} = k_\theta (\theta - \theta_0)^2
\end{equation}

where $\theta$ is the angle at vertex $j$ between atoms $i$--$j$--$k$:
\begin{equation}
\theta = \arccos\!\left(\frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{kj}}{|\mathbf{r}_{ij}| \, |\mathbf{r}_{kj}|}\right)
\end{equation}

\paragraph{Force derivation} (lines 153--171):

The chain rule gives:
\begin{equation}
\mathbf{F}_i = -\nabla_{\mathbf{r}_i} U = -2k_\theta(\theta - \theta_0) \frac{\partial\theta}{\partial\mathbf{r}_i}
\end{equation}

The geometric derivative is:
\begin{equation}
\frac{\partial\theta}{\partial\mathbf{r}_i} = \frac{1}{\sin\theta}\left[\frac{\mathbf{r}_{kj}}{|\mathbf{r}_{ij}||\mathbf{r}_{kj}|} - \frac{\mathbf{r}_{ij}\cos\theta}{|\mathbf{r}_{ij}|^2}\right]
\end{equation}

Code implementation:

\begin{lstlisting}[language=C++]
double cos_theta = dot(rij, rkj) / (rij_len * rkj_len);
cos_theta = std::max(-1.0, std::min(1.0, cos_theta));  // Clamp
double theta = std::acos(cos_theta);
double dtheta = theta - ang.theta0;

U += ang.ktheta * dtheta * dtheta;

double sin_theta = std::sin(theta);
if (std::abs(sin_theta) < 1e-6) continue;  // Linear, skip

double k = -2.0 * ang.ktheta * dtheta / sin_theta;

Vec3 fi = (rkj * (1.0/(rij_len*rkj_len)) 
         - rij * (cos_theta/(rij_len*rij_len))) * k;
Vec3 fk = (rij * (1.0/(rij_len*rkj_len)) 
         - rkj * (cos_theta/(rkj_len*rkj_len))) * k;
Vec3 fj = (fi + fk) * (-1.0);  // Momentum conservation
\end{lstlisting}

\paragraph{Default parameters} (header comment, line 22):
\begin{itemize}
\item $k_\theta \sim 50$--$100$\,kcal/(mol·rad$^2$)
\item $\theta_0 = 109.5^\circ = 1.911$\,rad (tetrahedral)
\end{itemize}

\subsection{Periodic Torsions (Dihedrals)}

\paragraph{Potential} (line 184):

\begin{equation}
U_{\text{tors}} = V_n[1 + \cos(n\phi - \gamma)]
\end{equation}

where:
\begin{itemize}
\item $\phi$ is the dihedral angle (rotation about $j$--$k$ bond)
\item $n$ is the periodicity ($n = 1, 2, 3, 6$)
\item $V_n$ is the barrier height (kcal/mol)
\item $\gamma$ is the phase offset (radians)
\end{itemize}

\paragraph{Dihedral angle computation} (lines 19--40):

The dihedral angle is computed from two plane normals:
\begin{align}
\mathbf{n}_1 &= \mathbf{b}_1 \times \mathbf{b}_2 \quad \text{(normal to plane } i\text{--}j\text{--}k\text{)} \\
\mathbf{n}_2 &= \mathbf{b}_2 \times \mathbf{b}_3 \quad \text{(normal to plane } j\text{--}k\text{--}l\text{)}
\end{align}

where $\mathbf{b}_1 = \mathbf{r}_j - \mathbf{r}_i$, $\mathbf{b}_2 = \mathbf{r}_k - \mathbf{r}_j$, $\mathbf{b}_3 = \mathbf{r}_l - \mathbf{r}_k$.

The angle is:
\begin{equation}
\phi = \text{atan2}\!\left(\frac{\mathbf{b}_2 \cdot (\mathbf{n}_1 \times \mathbf{n}_2)}{|\mathbf{b}_2|}, \; \mathbf{n}_1 \cdot \mathbf{n}_2\right)
\end{equation}

This gives $\phi \in [-\pi, \pi]$ with correct sign convention.

\paragraph{Force distribution} (lines 42--118):

The torsional forces use the Blondel-Karplus (1996) formulation, which avoids singularities when the dihedral angle approaches $0$ or $\pi$ (planar configurations).

The geometric derivatives are:
\begin{align}
\frac{\partial\phi}{\partial\mathbf{r}_i} &= -\frac{|\mathbf{b}_2|}{|\mathbf{n}_1|^2}\,\mathbf{n}_1 \\
\frac{\partial\phi}{\partial\mathbf{r}_l} &= +\frac{|\mathbf{b}_2|}{|\mathbf{n}_2|^2}\,\mathbf{n}_2
\end{align}

For the central atoms $j$ and $k$, the derivatives involve both plane normals and bond angles. The complete expressions are implemented in lines 90--99:

\begin{lstlisting}[language=C++]
double a1 = -r2 / n1_sq;
double a2 = (r2 - r1*cos_theta1) / (r2*n1_sq);
double a3 = cos_theta1 / (r2*n1_sq);
double a4 = cos_theta2 / (r2*n2_sq);
double a5 = (r2 - r3*cos_theta2) / (r2*n2_sq);
double a6 = r2 / n2_sq;

fi = n1 * (a1 * dU_dphi);
fj = n1 * (a2 * dU_dphi) + n2 * (a4 * dU_dphi);
fk = n1 * (-a3 * dU_dphi) + n2 * (a5 * dU_dphi);
fl = n2 * (a6 * dU_dphi);
\end{lstlisting}

\paragraph{Default parameters} (header comment, line 26):
\begin{itemize}
\item $n = 3$ (sp$^3$--sp$^3$ rotation, e.g., ethane)
\item $V_n \sim 0.5$--$3$\,kcal/mol
\item $\gamma = 0$ (trans minimum) or $\pi$ (cis minimum)
\end{itemize}

\subsection{Improper Torsions}

\paragraph{Purpose} Maintain planarity for sp$^2$ centers (aromatic rings, amide bonds).

\paragraph{Potential} (line 203):

\begin{equation}
U_{\text{imp}} = k_{\text{imp}}(\psi - \psi_0)^2
\end{equation}

where $\psi$ is the out-of-plane angle.

\paragraph{Implementation status} The code comment (line 205) notes:

\begin{quote}
\texttt{// Simplified: use dihedral formulation with harmonic potential}
\end{quote}

This means impropers are computed using the same dihedral angle machinery, but with harmonic potential instead of periodic.

\paragraph{Default parameters} (header comment, line 31):
\begin{itemize}
\item $k_{\text{imp}} \sim 10$--$50$\,kcal/(mol·rad$^2$)
\item $\psi_0 = 0$ (planar)
\end{itemize}

\subsection{Topology Auto-Generation}

From \texttt{atomistic/models/bonded.hpp}, lines 117--122:

The bonded topology can be inferred automatically from the bond graph \texttt{State.B}:

\begin{lstlisting}[language=C++]
// Auto-generate angles from bond graph (all i-j-k triplets)
void generate_angles_from_bonds();

// Auto-generate dihedrals from bond graph (all i-j-k-l quartets)
void generate_dihedrals_from_bonds();

// Assign default parameters based on element types
void assign_default_parameters(const State& s);
\end{lstlisting}

This means the engine can process any molecular topology without a pre-built parameter file. The trade-off:
\begin{itemize}
\item All angles get the same force constant and equilibrium angle (tetrahedral)
\item All dihedrals get the same barrier
\item Element-specific parameterization (using covalent radii for $r_0$, electronegativity for $k_b$) is a future extension
\end{itemize}

\section{Langevin Dynamics (NVT Ensemble)}

\textbf{Status:} Partially implemented in \texttt{atomistic/integrators/velocity\_verlet.cpp}, but not as a standalone Langevin class. The framework includes velocity rescaling (Berendsen thermostat) instead.

\subsection{Berendsen Velocity Rescaling}

\textbf{Implementation:} \texttt{atomistic/integrators/velocity\_verlet.cpp}, lines 65--75

An alternative weak-coupling thermostat:

\begin{lstlisting}[language=C++]
void rescale_velocities(State& state, double T_target) {
    double T_current = compute_temperature(state);
    
    if (T_current < 1e-6) return;  // Avoid division by zero
    
    double scale = std::sqrt(T_target / T_current);
    
    for (auto& v : state.V) {
        v.x *= scale;
        v.y *= scale;
        v.z *= scale;
    }
}
\end{lstlisting}

The velocities are rescaled by:
\begin{equation}
\lambda = \sqrt{\frac{T_{\text{target}}}{T_{\text{current}}}}, \qquad \mathbf{v}_i \gets \lambda \cdot \mathbf{v}_i
\end{equation}

\paragraph{Usage} Applied every $n$ steps (default: $n = 10$) during velocity Verlet integration.

\paragraph{Limitation} The Berendsen thermostat drives $T \to T_{\text{target}}$ exponentially but does \textbf{not} generate correct canonical fluctuations. The kinetic energy distribution is too narrow. Use only for equilibration, not for production sampling or heat capacity calculations.

\subsection{Full Langevin Dynamics: Test Implementation}

The Langevin equation:
\begin{equation}
m_i \frac{d\mathbf{v}_i}{dt} = \mathbf{F}_i(\mathbf{x}) - \gamma m_i \mathbf{v}_i + \sqrt{2\gamma m_i k_B T}\,\boldsymbol{\xi}_i(t)
\end{equation}

with friction $\gamma$ and random force $\boldsymbol{\xi}_i(t) \sim \mathcal{N}(0,1)$ is implemented in the BAOAB splitting scheme as a \textbf{standalone test application} (\texttt{apps/test\_baoab\_multi.cpp}) but is \textbf{not yet integrated} into the main formation pipeline.

\paragraph{Current status:} BAOAB test produces correct temperature (298 K target, 298.3 K measured, 0.6\% error) for Ar systems. Integration into the main velocity Verlet loop requires wiring the random force generation and friction coefficient into the pipeline parameters.

\section{Integrator Selection Flowchart}

\begin{center}
\begin{tabular}{lll}
\toprule
Goal & Ensemble & Integrator \\
\midrule
Energy minimization & — & FIRE \\
Conservative dynamics (energy conserved) & NVE & Velocity Verlet \\
Temperature-controlled (weak coupling) & NVT & Velocity Verlet + Berendsen \\
Canonical sampling (correct fluctuations) & NVT & BAOAB (planned) \\
Formation detection & — & MD + FIRE \\
\bottomrule
\end{tabular}
\end{center}

\section{Timestep Selection Guidelines}

The choice of timestep $\Delta t$ depends on the fastest motion in the system:

\begin{center}
\begin{tabular}{lll}
\toprule
System Type & Fastest Motion & Recommended $\Delta t$ \\
\midrule
Noble gas (Ar, Xe) & LJ oscillation & 1--2\,fs \\
Organic molecules (C, H, O, N) & C--H stretch & 0.5--1\,fs \\
Systems with H atoms & H--X stretch ($\sim 3000$\,cm$^{-1}$) & 0.5\,fs \\
Stiff bonds (SHAKE/RATTLE) & Constrained bonds & 2\,fs \\
Water (flexible) & O--H stretch & 0.5\,fs \\
Water (rigid, TIP3P) & Angle bending & 2\,fs \\
\bottomrule
\end{tabular}
\end{center}

\paragraph{Rule of thumb} The timestep should be $\sim 10\times$ smaller than the period of the fastest vibrational mode:
\begin{equation}
\Delta t < \frac{T_{\min}}{10} = \frac{1}{10\nu_{\max}}
\end{equation}

For a C--H stretch at $\nu = 3000$\,cm$^{-1}$:
\begin{equation}
T = \frac{1}{3000\,\text{cm}^{-1} \times 3 \times 10^{10}\,\text{Hz/cm}^{-1}} \approx 11\,\text{fs} \;\Rightarrow\; \Delta t < 1\,\text{fs}
\end{equation}

\section{Validation and Sanity Checks}

\subsection{NVE Energy Conservation Test}

\textbf{Test:} Run velocity Verlet on Ar gas (100 atoms, 300\,K) for $10^4$ steps with $\Delta t = 1$\,fs

\textbf{Expected result:}
\begin{equation}
\frac{|\Delta E|}{N} < 10^{-3}\,\text{kcal/mol per atom}
\end{equation}

\textbf{Status:} Validated for Lennard-Jones systems (code comment line 29)

\subsection{FIRE Convergence Test}

\textbf{Test:} Initialize H$_2$O molecule with perturbed geometry, minimize with FIRE

\textbf{Expected result:}
\begin{itemize}
\item Converges in $< 100$ iterations
\item Final $F_{\text{rms}} < 10^{-6}$\,kcal/(mol·\AA)
\item O--H bond length: $0.95 \pm 0.02$\,\AA
\item H--O--H angle: $104 \pm 2^\circ$
\end{itemize}

\textbf{Status:} Should be tested systematically (not found in test suite)

\subsection{Bonded Force Validation}

\textbf{Test:} Finite-difference check for each bonded term

\textbf{Method:}
\begin{enumerate}
\item Compute analytic force $\mathbf{F}_i$
\item Displace atom $i$ by $\delta = 10^{-5}$\,\AA\ in direction $\alpha$
\item Compute numerical force: $F_{i,\alpha} \approx -[U(\mathbf{X} + \delta\mathbf{e}_\alpha) - U(\mathbf{X})]/\delta$
\item Compare: $|F_{\text{analytic}} - F_{\text{numeric}}| < 10^{-4}$\,kcal/(mol·\AA)
\end{enumerate}

\textbf{Status:} Test files exist (\texttt{tests/angle\_tests.cpp}, \texttt{tests/torsion\_validation\_tests.cpp})

\section{Known Limitations}

\subsection{Velocity Verlet + Coulomb Instability}

As documented in Section~3.3 and Section~4.10:
\begin{itemize}
\item Coulomb forces are disabled in the MD integrator
\item Ionic systems diverge if Coulomb forces are enabled
\item Root cause: incorrect acceleration conversion factor for large Coulomb forces
\item Workaround: use FIRE for ionic formations (does not suffer from this issue)
\end{itemize}

\subsection{Langevin Dynamics: Test Only}

\begin{itemize}
\item BAOAB Langevin integrator exists as test application (\texttt{apps/test\_baoab\_multi.cpp})
\item Produces correct temperature for Ar systems (validated, 0.6\% error)
\item Not yet wired into formation pipeline (requires parameter plumbing)
\item Until integrated, only Berendsen velocity rescaling is available for temperature control
\end{itemize}

\subsection{Bonded Topology: Generic Parameters Only}

\begin{itemize}
\item All bonds get the same $k_b$ and $r_0$ (generic C--C-like)
\item All angles get the same $k_\theta$ and $\theta_0$ (tetrahedral)
\item All dihedrals get the same $V_n$
\item Element-specific parameterization is a future extension
\end{itemize}

\subsection{No Constraint Dynamics}

\begin{itemize}
\item SHAKE/RATTLE (bond constraints) not implemented
\item Rigid water models (TIP3P, SPC/E) not supported
\item All bonds are flexible
\end{itemize}

\subsection{Improper Torsions: Simplified}

\begin{itemize}
\item Impropers use dihedral angle computation + harmonic potential
\item True out-of-plane angle (distance from plane) not implemented
\item Works for most cases but may be inaccurate for highly non-planar geometries
\end{itemize}

\section{Conclusion: The Integration Layer as Implemented}

This section has documented the time integration machinery \textbf{as it exists in the codebase}:

\begin{itemize}
\item \textbf{Velocity Verlet}: Fully implemented for NVE dynamics (lines 87--220), with Coulomb forces disabled due to numerical instability
\item \textbf{FIRE minimization}: Fully implemented (lines 1--92), robust and parameter-insensitive, works for all force fields including ionic systems
\item \textbf{Bonded force field}: Four terms implemented (bonds, angles, dihedrals, impropers) with explicit Blondel-Karplus force distribution for torsions
\item \textbf{Berendsen thermostat}: Available but not true canonical sampling
\item \textbf{BAOAB Langevin}: Implemented and validated as test application; integration into pipeline pending
\end{itemize}

The integration pipeline is \textbf{functional for the formation problem}: MD explores, FIRE quenches, formations are classified. The Coulomb limitation affects only ionic MD, not formation detection (which uses FIRE).

Every algorithm is \textbf{explicit, auditable, and reproducible} --- the formation engine core principles.

\paragraph{Transition to \S6:} With the interaction physics (\S3) and integration machinery (\S5) established, Section~6 documents how these components combine to produce \textit{formations}: the stable structures that emerge from the MD exploration + FIRE quenching pipeline.

\end{document}