BAGOGAB integrator for in-triangle points
When a simulation contains in-triangle points, baoab! does not run the textbook BAOAB Langevin splitting. It runs BAGOGAB — BAOAB with an extra symmetric pair of “G” sub-steps inserted between A and O. This page derives the Hamiltonian that the G step integrates, describes how the splitting maps to in-triangle motion on the membrane mesh, and explains how the velocity thermostat handles in-triangle points.
Motivation
A naive way to keep a particle stuck to a triangle mesh is to add an attractive point-to-membrane force pulling it onto the surface, but then the integrator has to track which triangle each particle currently belongs to and keep the attracting force pointed at the right one as the mesh deforms and the point moves — a bookkeeping problem that gets ugly fast.
A second option is to step the particle without holding it to the mesh and then enforce the in-triangle constraint with a SHAKE-style solve. Done correctly, this solve has to couple across the entire connected mesh: a constraint impulse on a point requires an equal-and-opposite impulse on the three vertices of its current triangle, which changes the motion of those vertices and can violate the constraints of every other in-triangle point on the same triangle — or on any neighbouring triangle, since the vertices are shared along edges. A future direction could be an approximate, local constraint solver that avoids this global iteration.
BAGOGAB sidesteps both. We construct a Hamiltonian (derived below) in which the in-triangle motion has its own kinetic term. The G step then exactly integrates pure geodesic drift on the frozen mesh — a cheap per-particle operation with no global solve and no attracting force — at the cost of an additional pair of half-steps inside the BAOAB splitting.
Setup: one triangle with one in-triangle point
Consider a single triangle in 3-space with vertices \(R_1, R_2, R_3\), vertex masses \(M_1, M_2, M_3\), and vertex momenta \(P_1, P_2, P_3\). A point particle of mass \(m\) lives on the triangle. Its 3D position \(r\) is determined by its barycentric coordinates \((\lambda_1, \lambda_2, \lambda_3)\) constrained by \(\lambda_1 + \lambda_2 + \lambda_3 = 1\):
\[ r = \lambda_1 R_1 + \lambda_2 R_2 + \lambda_3 R_3. \]
Take \((\lambda_2, \lambda_3)\) as the independent coordinates with \(\lambda_1 \equiv 1 - \lambda_2 - \lambda_3\). Define the edge matrix and metric:
\[ A = \begin{bmatrix} R_2 - R_1 & R_3 - R_1 \end{bmatrix}, \qquad G = A^{\mathsf T} A. \]
Then \(r = R_1 + A \lambda\) with \(\lambda = (\lambda_2, \lambda_3)\).
There is a potential energy \(PE(R_1, R_2, R_3, r)\). Write \(f = -\partial PE / \partial r\) for the force on the particle and \(F_i = -\partial PE / \partial R_i\) for the explicit vertex forces (the chain-rule contribution from \(r\) will appear separately below).
A non-physical but tractable Lagrangian
The “physically honest” particle velocity in 3-space is \(\dot r = \dot R_1 + \dot A\,\lambda + A\,\dot\lambda\). Using this in the kinetic energy couples the particle’s motion to the vertices’ motion in a way that is awkward to integrate symplectically. Instead, in the spirit of the Andersen barostat (Andersen 1980), drop the \(\dot R_1 + \dot A\,\lambda\) contribution and use only the in-plane particle velocity, ignoring vertex motion:
\[ v := A \dot\lambda = \dot\lambda_1 R_1 + \dot\lambda_2 R_2 + \dot\lambda_3 R_3, \qquad \dot\lambda_1 := -(\dot\lambda_2 + \dot\lambda_3). \]
The resulting Lagrangian is
\[ \mathcal{L} = \tfrac{1}{2} m\,(A\dot\lambda)^{\mathsf T} (A\dot\lambda) + \tfrac{1}{2} \sum_i M_i\,\dot R_i^{\,2} - PE(R_1, R_2, R_3, R_1 + A\lambda). \]
This kinetic term is not physically correct — it ignores the motion the particle inherits from the moving vertices — but, as in the Andersen barostat, integrating out the momenta at temperature \(kT\) gives the configurational marginal \(\rho \propto \sqrt{\det G}\,\exp(-PE/kT)\). The \(\sqrt{\det G}\) is exactly the triangle area Jacobian \(dA = \sqrt{\det G}\,d\lambda_2\,d\lambda_3\), so in physical \(r\)-space the equilibrium is the correct \(\exp(-PE/kT)\,dA\).
Hamiltonian
Conjugate momenta:
\[ p_\lambda = \partial \mathcal{L}/\partial \dot\lambda = m\,G\,\dot\lambda, \qquad P_i = \partial \mathcal{L}/\partial \dot R_i = M_i\,\dot R_i. \]
Legendre transform gives
\[ \mathcal{H} = \frac{1}{2 m}\,p_\lambda^{\mathsf T} G^{-1} p_\lambda + \sum_i \frac{1}{2 M_i}\,P_i^{\,2} + PE(R_1, R_2, R_3, R_1 + A\lambda). \]
Equations of motion
Differentiating \(\mathcal{H}\):
\[ \begin{aligned} \dot\lambda &= \partial \mathcal{H}/\partial p_\lambda &&= \tfrac{1}{m}\,G^{-1} p_\lambda, \\ \dot{p}_\lambda &= -\partial \mathcal{H}/\partial \lambda &&= A^{\mathsf T} f, \\ \dot R_i &= \partial \mathcal{H}/\partial P_i &&= \tfrac{1}{M_i}\,P_i, \\ \dot P_i &= -\partial \mathcal{H}/\partial R_i &&= F_i + \lambda_i\,f + m\,\dot\lambda_i\,v, \\ v &= A\dot\lambda &&= \tfrac{1}{m}\,A\,G^{-1} p_\lambda. \end{aligned} \]
Three contributions to \(\dot P_i\):
- \(F_i\) — explicit vertex-vertex potential gradient (e.g. neighbouring edge springs).
- \(\lambda_i\,f\) — chain-rule kick from \(r = R_1 + A\lambda\). The particle pushes on each vertex weighted by that vertex’s barycentric coordinate.
- \(m\,\dot\lambda_i\,v\) — the recoil term, from differentiating \(G^{-1}\) with respect to \(R_i\). It is the configurational force the vertices must feel to reproduce the equilibrium distribution of a real massive particle constrained to the surface — the “2D gas” surface tension reduction that pushes outward on the mesh and bounces off triangle boundaries.
Because \(\sum_i \dot\lambda_i = 0\) and \(\sum_i \lambda_i = 1\), summing \(\dot P_i\) over the three vertices gives \(\sum_i F_i + f\): the vertex springs cancel pairwise and the recoils sum to zero, so total momentum tracks just the net external force on the system.
Scaling up to the full MEDYAN system
The toy above is one triangle plus one in-triangle point. A real simulation is much bigger. The state is organised around beads — the 3D positions that the user-defined \(PE\) is a function of. Concretely a bead is one of:
- a filament node;
- a filament material direction unit vector (twistable filaments);
- a membrane vertex;
- a ball center;
- an in-triangle point, whose 3D position \(r = R_1 + A\lambda\) is derived from its current triangle’s vertices and its barycentric coords.
A bead is a 3D position, not a DOF. The independent DOFs are a subset or reparametrisation of the bead coordinates — three Cartesian coords for an unconstrained bead, but only two barycentric coords for an in-triangle point, and material direction unit vectors and filament-length constraints remove further DOFs.
Each bead carries a 3D velocity, a mass \(M\), and a friction \(\gamma\). The user-defined \(PE\) is a single function of every bead’s 3D position — filament stretching/bending/twisting, membrane bending/area, ball–object collisions, external fields, etc. — and the conservative force on a bead is \(F = -\partial PE / \partial(\text{its 3D position})\).
Further changes from the toy:
- The membrane is a triangle mesh. Instead of a single isolated triangle, a membrane is a collection of triangles glued along shared edges. An in-triangle point lives in one current triangle at a time and must be able to transfer between neighboring triangles as smoothly as possible as it moves.
- Langevin dynamics at temperature \(T\). To sample Langevin dynamics and control the diffusion coefficients we add an Ornstein–Uhlenbeck velocity update — the O step — to each bead independently.
Splitting strategy
Write the Hamiltonian as three additive pieces and a friction/noise piece:
B — conservative potential \(PE\) (a function of every bead’s 3D position). Generates a kick on every bead’s momentum: \(\dot P_b = -\partial PE/\partial R_b\) for an unconstrained Cartesian bead, and for each in-triangle point \(\dot p_\lambda = A^{\mathsf T} f\) together with the chain-rule contribution \(\lambda_i\,f\) to its three triangle vertices.
A — free-bead kinetic energy \(\sum_b P_b^{\,2}/(2 M_b)\), summed over all non-intri beads. Generates uniform drift of each bead at fixed velocity. Filament beads carry holonomic constraints (fixed link length, plus material-direction perpendicularity for twistable filaments); these are enforced by a SHAKE-style solver that, per filament, sets up a block-tridiagonal matrix system and Newton-iterates the constraint impulses on velocity so that the post-drift positions land on the constraint manifold. This is not a simple projection: the constraints along a single filament are coupled and have to be solved jointly. For each in-triangle point, \(\lambda\) and \(p_\lambda\) are canonical constants of this flow, but its in-plane 3D velocity \(v = A\,G^{-1}\,p_\lambda / m\) changes because the supporting triangle’s vertices (and hence \(A\)) have moved.
G — in-triangle kinetic energy \(\sum_{\mathrm{intri}} p_\lambda^{\mathsf T} G^{-1} p_\lambda / (2 m)\), summed over every in-triangle point. For each one: generates \(\dot\lambda = G^{-1} p_\lambda / m\) — straight-line geodesic motion of the particle along the frozen mesh — and a recoil \(\dot P_i = m\,\dot\lambda_i\,v\) on the three vertices of the current triangle.
When the geodesic reaches an edge of the current triangle, the particle’s 3D velocity is rotated along the edge to be in the plane of the new triangle, or reflected if that edge is a mesh border. In-triangle kinetic energy \(\tfrac12 m\,|v|^2\) is conserved across the transfer. This transfer is not \(C^1\)-smooth, and there are two independent sources of kink. First, the particle’s velocity direction jumps wherever the edge has a nonzero dihedral angle — at curved interior edges, and at every border reflection. Second, even on a perfectly flat mesh the set of vertices receiving the recoil \(m\,\dot\lambda_i\,v\) changes discretely at the crossing: the old apex stops being recoiled and the new apex starts, so each affected vertex’s \(\dot P_i\) has a step discontinuity and \(P_i\) has a kink. Kinks like these are known to cause long-time energy-conservation problems for symplectic integrators; see Bond and Leimkuhler (2007) for analysis and stabilized schemes in the hard-sphere setting.
O — Ornstein–Uhlenbeck thermostat (added on top of the splitting of \(\mathcal{H}\)). The same isotropic 3D rule applies to every bead’s velocity:
\[ u \leftarrow c_1\,u + c_2\,\xi, \qquad c_1 = \exp(-\gamma\,dt / M), \qquad c_2^{\,2} = \tfrac{k_B T}{M}\,(1 - c_1^{\,2}), \]
with \(\xi\) a fresh standard normal 3-vector, \(M\) the bead’s mass, and \(\gamma\) its friction coefficient (MEDYAN units pN·s/nm). \(c_2^{\,2}\) is fixed by fluctuation–dissipation at temperature \(T\).
For in-triangle points the canonical variable is \(p_\lambda\), with equilibrium covariance \(m\,k_B T\,G\), so the textbook OU step \(p_\lambda \leftarrow c_1\,p_\lambda + \sqrt{m\,k_B T(1 - c_1^{\,2})}\,G^{1/2}\xi_2\) would need the matrix square root \(G^{1/2}\). The generic rule on \(v\) is the same step in disguise: applying \(p_\lambda = m\,A^{\mathsf T} v\) to \(v^{\mathrm{new}} = c_1 v + c_2 \xi\) gives \(p_\lambda^{\mathrm{new}} = c_1\,p_\lambda + c_2\,m\,A^{\mathsf T}\xi\), with the new noise having covariance \(c_2^{\,2}\,m^2\,A^{\mathsf T} A = m\,k_B T(1 - c_1^{\,2})\,G\) — exactly matching the \(G^{1/2}\xi_2\) noise above. The \(A^{\mathsf T}\) projection does the work of \(G^{1/2}\) for free and discards the out-of-plane part of the 3D draw on its next application (the following G or A step).
Applying a symmetric Strang splitting one step at a time yields the sequence
\[ B(\tfrac{dt}{2})\;A(\tfrac{dt}{2})\;G(\tfrac{dt}{2})\;O(dt)\;G(\tfrac{dt}{2})\;A(\tfrac{dt}{2})\;B(\tfrac{dt}{2}) \]
i.e. BAGOGAB. With no in-triangle points, the G sub-steps are empty and the sequence collapses to plain BAOAB, which is why the public function is still called baoab!.