The algorithms desribed in this article were implemented in my work-in-progress TaichiMD simulation package. For a brief analysis on time integrators in molecular dynamics simulations, see my post in the Taichi forum.
Hessian in implicit time integrators
In some physical simulation algorithms, we are interested in derivative of the force with respect to its position, for example, a backward Euler integrator solved as a linear system:
\[\mathbf{x}_{t+1}=\mathbf{x}_t+\Delta t \mathbf{v}_{t+1}, \quad \mathbf{v}_{t+1}=\mathbf{v}_t+\Delta t M^{-1}\mathbf{f}(\mathbf{x}_{t+1})\]where the velocity after the time step can be solved after a 1st-order Taylor expansion of the force
\[\mathbf{v}_{t+1}=\mathbf{v}_t+\Delta t M^{-1}\mathbf{f}(\mathbf{x}_t+\Delta t \mathbf{v}_{t+1})=\mathbf{v}_t+\Delta t M^{-1}\Big[\mathbf{f}(\mathbf{x}_t)+\frac{\partial \mathbf{f}}{\partial \mathbf{x}}(\mathbf{x}_t)\Delta t\mathbf{v}_{t+1}\Big]\]which will give us a linear system as
\[\Big[\mathbf{I}-(\Delta t)^2M^{-1}\frac{\partial\mathbf{f}}{\partial\mathbf{x}}(\mathbf{x}_t)\Big]\mathbf{v}_{t+1}=\mathbf{v}_t+\Delta tM^{-1}\mathbf{f}(\mathbf{x}_t)\]Then \(\mathbf{v}_{t+1}\) can be obtained by solving this linear system.
For a particle-based simulation with a well-defined potential energy \(U(\mathbf{x})\), the force is defined as \(\mathbf{f}=-\partial U / \partial\mathbf{x}\), and the Hessian is defined as \(\mathbf{H}=-\partial^2U/\partial\mathbf{x}^2\). Therefore, the linear system becomes
\[\Big[\mathbf{I}+(\Delta t)^2M^{-1}\mathbf{H}(\mathbf{x}_t)\Big]\mathbf{v}_{t+1}=\mathbf{v}_t+\Delta tM^{-1}\mathbf{f}(\mathbf{x}_t)\]Calculating the Hessian for pairwise systems
Molecular simulation systems are predominantly pairwise. That is, the potential energy for the simulation system is calculated as the sum of the interaction energies between particle pairs (in molecular simulations:
\[U=\sum_i\sum_{j<i}u(r_{ij}^2)\]where \(r_i\) is the position of the i th particle (instead of \(x_i\) due to molecular simulation convention) and \(u(\cdot)\) is a function operating on the square magnitude of vector \(r_{ij}=r_j-r_i\). Notable examples of \(u(\cdot)\) include the Lennard-Jones potential for intermolecular interactions,
\[u(r_{ij})=4\epsilon\Big[(\frac{\sigma}{r_{ij}})^{12}-(\frac{\sigma}{r_{ij}})^6\Big]\]the Coulomb potential for electrostatic interactions,
\[u(r_{ij})=\frac{1}{4\pi\epsilon_0}\frac{q_iq_j}{r_{ij}}\]and the harmonic potential for chemical bonds (mass-spring system)
\[u(r_{ij})=\frac{k}{2}(r_{ij}-r_0)^2\]It would be useful to implement a general algorithm computing the Hessian for any pairwise potential, then new types of pairwise interactions can be conveniently implemented by only defining the potential energy and its derivatives with respect to the squared distance.
For consistency while working with matrix calculus, we assume a column vector notation \(r_i=(x_i, y_i, z_i)^T\). The partial derivatives for the squared distance are
\[\frac{\partial(r_{ij})^2}{\partial r_i}=\frac{\partial}{\partial r_i}(r_j-r_i)^T(r_j-r_i)=-2(r_j-r_i), \quad \frac{\partial(r_{ij})^2}{\partial r_j}=2(r_j-r_i)\]which gives us the first-order derivative of the potential energy with respect to the i th particle
\[\frac{\partial U}{\partial r_i}=\sum_{j\neq i}\frac{\partial u}{\partial r_i}=\sum_{j\neq i}u^{(1)}\frac{\partial (r_{ij})^2}{\partial r_i}=-2\sum_{j\neq i}u^{(1)}(r_j-r_i)\]Then an “element” in the Hessian is a 3x3 matrix can be obtained as the second-order derivative of the potential energy:
\[H_{ij}=\frac{\partial^2U}{\partial r_ir_j^T}=-2\frac{\partial}{\partial r_j^T}\sum_{k\neq i}u^{(1)}(r_k-r_i)\]This equations will give different results between diagonal and non-diagonal terms. For diagonal terms,
\[H_{ii}=2\sum_{k\neq i}[2u^{(2)}(r_k-r_i)(r_k-r_i)^T+u^{(1)}\mathbf{I}]\]For off-diagonal terms,
\[H_{ij}=H_{ji}=-2[2u^{(2)}(r_j-r_i)(r_j-r_i)^T+u^{(1)}\mathbf{I}]\]For a sanity check, the center of mass \(r_c=\sum_ir_i/n\) in a pairwise system should be stationary, which means the sum over each force term \(f_i\) and the sum over each Hessian term \(H_{ij}\) should be zero. This property is satisfied according to our derivations because the row and column sums of \(\mathbf{H}\) is zero.