Please upload one PDF file to Gradescope (Entry code: G64YK2).
For this homework you will implement code (on PrairieLearn) to implement the velocity Verlet and Langevin dynamics integration schemes used in molecular dynamics simulations. We’ll try this out with a simple model system: the one-dimensional harmonic oscillator.
A one-dimensional harmonic oscillator consists of a particle tethered to the origin by a spring. The particle has one position coordinate \(x\) and one momentum coordinate \(p\). The spring is represented by a harmonic potential \(U = k x^2/2\) having spring constant \(k\), and the kinetic energy of the particle is \(K = p^2/(2m)\). Hence, the Hamiltonian is \[ \mathcal{H} = \frac{p^2}{2m} + \frac{k}{2} x^2. \]
As part of your homework, you will need to compare your results to the analytical, exact solution. For this, use numpy arrays with position and momentum with 361 steps each that cover 10 oscillation periods of the harmonic oscillator. Use \(\varepsilon\) to represent the unit of energy, \(\ell\) to represent the unit of length, and \(m\) to represent the unit of mass; the unit of time is \(\tau = \sqrt{m \ell^2/\varepsilon}\). The particle has mass \(1.0\,m\), the value of the spring constant is \(k = 1.0\,\varepsilon/\ell^2\), and the initial particle coordinates are \(x(0) = 1.0\,\ell\) and \(p(0) = 0.0\,m \ell/\tau\). Complete the following analysis and prepare a report.
Simulate the motion of the particle in the microcanonical (NVE) ensemble for 10 periods of oscillation using the velocity Verlet algorithm with an appropriate timestep \(\Delta t\). Plot \(x(t)\) and \(p(t)\) versus \(t\) and the phase-space trajectory \(p(t)\) versus \(x(t)\). Include the analytical solutions for harmonic oscillation in both plots.
Repeat the simulation for a few values of \(\Delta t\) that are both less than and greater than your chosen timestep, and comment on how your results qualitatively change. To quantify this change, compute and plot the root mean-squared error in \(x(t)\) between your simulations and the analytical solution as a function of \(\Delta t\).
Simulate the motion of the particle using Langevin dynamics with \(T=1.0\,\varepsilon/k_{\rm B}\) and \(\gamma=0.1\,m/\tau\) for \(10^4\,\tau\) using an appropriate timestep. Hint: the timestep is much smaller than in Part:1 NVE. Use the BAOAB algorithm [“Robust and efficient configurational molecular sampling via Langevin dynamics,” Leimkuhler and Matthews, J. Chem. Phys. 138, 174102 (2013) doi:10.1063/1.4802990]. Record the values of \(x\) and \(p\) every \(10\,\tau\). Plot the points in phase space \((x,p)\) and compare to your NVE result. Compute the average kinetic energy \(\langle K \rangle\) and compare to the expected value. Compute and plot the marginal probability distributions of \(x\) and \(p\), and compare to the expected distributions.
BAOAB algorithm. From Leimkuhler and Matthews, the “BAOAB” algorithm refers to how we split the three updates in Langevin dynamics. In Langevin dynamics, we’re discretizing the stochastic differential equation \[\mathrm{d}q = M^{-1} p\; \mathrm{d}t\] \[\mathrm{d}p = -\nabla U(q) \mathrm{d}t - \gamma p\; \mathrm{d}t + \sigma M^{1/2} \mathrm{d}W\] for positions \(q\), momenta \(p\), masses \(M\), potential energy \(U(q)\), and a \(W=W(t)\) a vector of \(3N\) independent Wiener processes (Gaussian white noise, uncorrelated in time), and \(\gamma>0\). With \(\sigma = \sqrt{2\gamma k_\text{B}T}\), this samples the canonical distribution. The parameter \(\gamma\) introduces a frictional force, while the Wiener processes are random “kicks” from the bath to the particles. Essentially, the friction removes memory from the system while the kicks move the momenta towards a proper Maxwell-Boltzman distribution.
To propagate this forward in time, there are three update steps to consider:
Thus, the full update is A+B+O. For A, there is no update to \(p\), and for B and O, there is no update to \(q\). Each of these three equations can be solved “exactly” independently; for example, with A, we have \(q(t+\Delta t) = q(t) + M^{-1}p\Delta t\), and for B, we have \(p(t+\Delta t) = p(t) - \nabla U(q(t))\Delta t\). For the O term, it is more complicated, but has the solution \[p(t+\Delta t) = e^{-\gamma\Delta t}p(t) + \frac{\sigma}{\sqrt{2\gamma}}\sqrt{1-e^{-2\gamma\Delta t}} M^{1/2} R_t\] where \(R_t \sim \mathcal{N}(0,1)\) is a vector of uncorrelated Gaussian noise.
As you remember from our discussion of velocity Verlet, we can do sequentially with small \(\Delta t\), but symmetrically. The velocity Verlet, which has no O update, amounts to doing a B update for \(\Delta t/2\), an A update for \(\Delta t\), then a final B update for \(\Delta t/2\). This would be written as “BAB”. To do the Langevin dynamics, you’ll do BAOAB in order to move forward by a time step of \(\Delta t\): a \(\Delta t/2\) with B, \(\Delta t/2\) with A, \(\Delta t\) with O, \(\Delta t/2\) with A, and finally \(\Delta t/2\) with B.