## Verlet algorithm Verlet algorithm is a method with accuracy of $O(t^4)$, this can be shown by [[Hamilton's principle]], but it is easier using Taylor expansion. For a position $\vec{r}$ (write as $r$ for simplification), we have $r(t+\delta t) = r(t) + \dot{r}(t)\delta t + \ddot{r}(t)\cdot \frac{1}{2} \delta t^2 + \frac{1}{3!}\dddot{r} \delta t^3 +O(\delta t^4)$ and $r(t-\delta t) = r(t) - \dot{r}(t)\delta t + \ddot{r}(t)\cdot \frac{1}{2} \delta t^2 - \frac{1}{3!}\dddot{r} \delta t^3 +O(\delta t^4)$ The summation gives $r(t+\delta t) + r(t-\delta t) = 2r(t)+\ddot{r}(t) \delta t^2 + 2O(\delta t^4)$ $\frac{\mathrm{d}^2 r(t)}{\mathrm{d}t^2} = \frac{r(t+\delta t) - 2r(t) + r(t-\delta t)}{\delta t^2}$ Since $\frac{\mathrm{d}^2r(t)}{\mathrm{d}t^2} = a = \frac{F(t)}{m} = -\nabla V \cdot \frac{1}{m}$, we have $r(t+\delta t) = 2r(t)-r(t-\delta t) + \frac{F(t)}{m}\delta t^2$ This is ***Verlet algorithm***, written in $q$, we have $q^{k+1} = 2q^k-q^{k-1} + \frac{F(t)}{m} \Delta t^2$^position-verlet >[!note] >Always remember that force $F$, acceleration $a$, potential $V$ (and $\dot{p}$) are interchangeable! If we (or say, simulation) begin by a position and a velocity, i.e., use $v$ and $a$ in the Taylor expansion, we have $r(t+\delta t) = r(t)+v(t)\delta t + \frac{1}{2}a(t)\delta t^2$ Do a **half step expansion for $v$**, we have $v\left( t+\delta t \cdot \frac{1}{2} \right) = v(t) +\frac{1}{2} a(t)\delta t$ ^1152db >[!notice] >This is an expansion for $v$! Not $r$! Since $a(t)$ here is just $\frac{\partial v}{\partial t} = \frac{F(t)}{m}$, it can be written as $v\left( t+\frac{1}{2}\delta t \right) = v(t) + \frac{F(t)}{m}$ Then update $F$ to $F(t+\delta t)$, and write expression for $v(t+\delta t)$, we get $\begin{aligned} v(t+\delta t) &=v\left( t+\frac{1}{2}\delta t \right) + \frac{\delta t}{2} \frac{\delta v\left( t+\frac{1}{2}\delta t \right)}{\delta t}\\ & = v\left( t+\frac{1}{2}\delta t \right) + \frac{\delta t}{2} \frac{F(k+1)}{m} \end{aligned}$ ^velocity-verlet >[!notice] >In above equations, $\frac{\delta t}{2} \frac{\delta v\left( t+\frac{1}{2}\delta t \right)}{\delta t}$ was replace by the updated force, i.e., $k+1$. > >This velocity Verlet is actually a *predictor-corrector* algorithm. The first expansion of $v\left( t+\delta t \cdot \frac{1}{2} \right)$ is a predictor, then this velocity was represented by force, and at the final step got replace, i.e., corrector. Comment on velocity Verlet: - fast - not accurate for long step - require little memory - good short-term energy conservation - small long-term energy drift Problems of position Verlet: - Lyapunov (in)stability: **acceptable error** decrease exponentially with simulation time ![[Drawing 2023-07-26 20.35.49.excalidraw.svg]] - Require frequent update of forces to keep accuracy - Energy drift, $E$ goes up with time due to numerical integration artifacts >[!notice] >Energy drift is a common phenomenon in MD, Verlet itself is time reversal, so $E$ conservation *should be* preserved. >[!note] >In the slide it says "The time stepping scheme results in a non-physical, limited sampling of motions with frequency close to frequency of velocity updates" -> time step not match the min time step of the system! >[!info] >See more at https://en.wikipedia.org/wiki/Verlet_integration