Question

In molecular dynamics simulations using classical mechanics modeling, one is often faced with a large nonlinear ODE system of the form $$ M \mathbf{q}^{\prime \prime}=\mathbf{f}(\mathbf{q}), \text { where } \mathbf{f}(\mathbf{q})=-\nabla U(\mathbf{q}) $$ Here \(\mathbf{q}\) are generalized positions of atoms, \(M\) is a constant, diagonal, positive mass matrix, and \(U(\mathbf{q})\) is a scalar potential function. Also, \(\nabla U(\mathbf{q})=\left(\frac{\partial U}{\partial q_{1}}, \ldots, \frac{\partial U}{\partial q_{m}}\right)^{T} .\) A small (and somewhat nasty) instance of this is given by the Morse potential where \(\mathbf{q}=q(t)\) is scalar, \(U(q)=D(1-\) \(\left.e^{-S\left(q-q_{0}\right)}\right)^{2}\), and we use the constants \(D=90.5 \cdot 0.4814 \mathrm{e}-3, S=1.814, q_{0}=1.41\), and \(M=\) \(0.9953 .\) (a) Defining the velocities \(\mathbf{v}=\mathbf{q}^{\prime}\) and momenta \(\mathbf{p}=M \mathbf{v}\), the corresponding first-order \(\mathrm{ODE}\) system for \(\mathbf{q}\) and \(\mathbf{v}\) is given by $$ \begin{aligned} \mathbf{q}^{\prime} &=\mathbf{v} \\ M \mathbf{v}^{\prime} &=\mathbf{f}(\mathbf{q}) \end{aligned} $$ Show that the Hamiltonian function $$ H(\mathbf{q}, \mathbf{p})=\mathbf{p}^{T} M^{-1} \mathbf{p} / 2+U(\mathbf{q}) $$ is constant for all \(t>0\). (b) Use a library nonstiff RK code based on a \(4(5)\) embedded pair such as MATLAB's ode 45 to integrate this problem for the Morse potential on the interval \(0 \leq t \leq 2000\), starting from \(q(0)=1.4155, p(0)=\frac{1.545}{48.888} M .\) Using a tolerance tol \(=1 . \mathrm{e}-4\), the code should require a little more than 1000 times steps. Plot the obtained values for \(H(q(t), p(t))-H(q(0), p(0))\). Describe your observations.

Short answer

The given problem involves a molecular dynamics problem described by a nonlinear ODE system. To prove that the Hamiltonian function is constant for all \(t > 0\), we calculated its derivative with respect to time and showed that it is zero. We then computed the partial derivatives of the Hamiltonian function with respect to \(\mathbf{q}\) and \(\mathbf{p}\) and further computed their time derivatives. Then, we used the chain rule to compute the time derivative of the Hamiltonian function, which was found to be zero. To solve the problem numerically and plot the change in Hamiltonian function over time, we proposed a procedure, including defining the constants, Morse potential and its gradient, ODE system, and implementing the ODE solver using MATLAB's "ode45". After running the code, we expect to observe that the Hamiltonian function remains constant (fluctuating around zero) throughout the whole integration, confirming the analytical result obtained in the first part of the problem.

Step-by-step solution

Apply the chain rule to Hamiltonian function derivatives

To show that the Hamiltonian function is constant we must show that its derivative with respect to time, \(dH(t)/dt\), is zero. Calculate \(dH/dt\) using the chain rule for the given Hamiltonian function \(H(\mathbf{q}, \mathbf{p})=\mathbf{p}^{T} M^{-1} \mathbf{p} / 2+U(\mathbf{q})\) and show it is zero.

Set up the derivatives

Compute the partial derivatives of the Hamiltonian function with respect to \(\mathbf{q}\) and \(\mathbf{p}\): \(\frac{\partial H}{\partial \mathbf{q}} = \nabla U(\mathbf{q})\) and \(\frac{\partial H}{\partial \mathbf{p}} = M^{-1}\mathbf{p}\) Compute the time derivative of \(\mathbf{q}\) and \(\mathbf{v}\): \(\frac{d\mathbf{q}}{dt} = \mathbf{v}\) and \(\frac{d\mathbf{v}}{dt}=\frac{d}{dt}(M^{-1}\mathbf{p}) = M^{-1}\frac{d\mathbf{p}}{dt}\). Since \(\mathbf{p} = M \mathbf{v}\), we get \(\frac{d\mathbf{p}}{dt} = M\frac{d\mathbf{v}}{dt} = \mathbf{f}(\mathbf{q})\) using the given equation.

Compute the Hamiltonian time derivative

Compute the time derivative of the Hamiltonian function using the chain rule: \(\frac{dH}{dt} = \frac{\partial H}{\partial \mathbf{q}}\frac{d\mathbf{q}}{dt} + \frac{\partial H}{\partial \mathbf{p}}\frac{d\mathbf{p}}{dt} = \nabla U(\mathbf{q})\mathbf{v} + M^{-1}\mathbf{p}\mathbf{f}(\mathbf{q})\) Now substitute \(\mathbf{f}(\mathbf{q}) = -\nabla U(\mathbf{q})\): \(\frac{dH}{dt} = \nabla U(\mathbf{q})\mathbf{v} - M^{-1}\mathbf{p}\nabla U(\mathbf{q}) = 0\) Since \(dH/dt = 0\), the Hamiltonian function is constant for all \(t > 0\). #b_Solve numerically and plot# For this part, we must implement the code using MATLAB's "ode45". We provide the procedure instead of the actual code. After implementing and running the code, you should be able to observe the plotted result.

Procedure to solve the ODE system numerically

To solve the ODE system numerically, follow these steps: 1. Define the constants: \(D, S, q_0, M\), and the given initial conditions for \(q(0)\) and \(p(0)\). 2. Define the Morse potential \(U(q)\) and its gradient \(\nabla U(\mathbf{q})\) using the given Morse potential formula. 3. Define the ODE system \(d\mathbf{q}/dt = \mathbf{v}\) and \(d\mathbf{v}/dt = M^{-1}\mathbf{f}(\mathbf{q})\). 4. Implement the ODE solver using MATLAB's "ode45", with a time interval from \(0\) to \(2000\), and a tolerance of \(1 \cdot 10^{-4}\). 5. Calculate the Hamiltonian function \(H(\mathbf{q}, \mathbf{p})\) for each time step in the obtained solution. 6. Plot the change in Hamiltonian function over time, i.e., plot \(H(q(t), p(t))-H(q(0), p(0))\).

Expected observation:

From the plotted result, you should be able to observe that the Hamiltonian function remains constant (fluctuating around zero) throughout the whole integration, confirming the analytical result obtained in part (a).

Key concepts

Nonlinear ODE Systems

In molecular dynamics simulations, one commonly encounters nonlinear ordinary differential equation (ODE) systems. These systems can be written in the form \[ M \mathbf{q}'' = \mathbf{f}(\mathbf{q}) \] where \( M \) represents a mass matrix, and \( \mathbf{q} \) denotes the generalized positions of atoms. The function \( \mathbf{f}(\mathbf{q}) \) denotes forces acting on these atoms, which are derived from a potential energy function, \( U(\mathbf{q}) \). This leads to nonlinearities because the forces potentially depend on the squared terms or other non-linear features of the positions.

Nonlinear ODE systems can describe a wide range of physical phenomena, including celestial mechanics and quantum systems. Due to their complexity, numerical simulations are often employed to find solutions. This offers insight into the behavior of molecular systems under various potentials, such as the Morse potential used in this particular problem.

Hamiltonian Function

The Hamiltonian function plays a crucial role in characterizing the dynamics of a physical system. In molecular dynamics, it combines kinetic and potential energy to encapsulate the total energy of the system. \[ H(\mathbf{q}, \mathbf{p}) = \frac{\mathbf{p}^T M^{-1} \mathbf{p}}{2} + U(\mathbf{q}) \]Here, \( \mathbf{p} \) is the momentum, which is a product of the mass matrix \( M \) and velocity \( \mathbf{v} \). The potential \( U(\mathbf{q}) \) derives from interactions between atoms.

For stable and conservative systems, the Hamiltonian remains constant over time. This conservation is significant because it ensures energy is neither added to nor removed from the system due to external factors. In practice, confirming that \( \frac{dH}{dt} = 0 \) validates the numerical methods and models used in simulations. The Hamiltonian is not a mere theoretical construct but a tool that helps validate and verify the accuracy of numerical experiments.

Numerical Integration

Numerical integration is a method used to approximate the solutions of both single and multiple-variable functions where analytic solutions are hard to find. In the context of molecular dynamics, these numerical methods solve the ODE system describing particle motion over time. The importance lies in the ability to simulate the progression of a system from known initial conditions over a specified time interval.

Techniques of numerical integration involve discretizing the continuous problem into small time steps, which are manageable computationally and allow one to approximate values at each point. Ensuring high precision while reducing computational cost is vital. Proper numerical integration helps achieve accurate predictions of phenomena such as molecular interactions and energy conservation in simulations.

Runge-Kutta Method

The Runge-Kutta method is a cornerstone technique for solving ODEs numerically. Specifically, the fourth-order Runge-Kutta method, often denoted as RK4, provides an effective balance of accuracy and computational efficiency. This method estimates the value of a future point using a weighted average of slopes (i.e., derivatives). This series of calculations offers superior precision over simpler methods like Euler's method.

An advanced version, the embedded Runge-Kutta method (such as the 4(5) pair used in MATLAB's ode45), provides adaptive step sizes based on error estimation. This means that the step size for each iteration is determined dynamically, allowing for more accurate results without unnecessary computations. In the context of our molecular dynamics problem, utilizing such methods ensures robust and faithful integration, maintaining the integrity of properties like energy conservation throughout the simulation.