M10 Angular momentum is conserved or not

Manual 10 for MD++
Angular momentum is conserved or not?

Keonwook Kang and Wei Cai

Jul 06 , 2009

---

In molecular dynamics simulations, positions and velocities of atoms are updated following Newton's equations of motion. Naturally, we would expect that linear and angular momenta (of an isolated system such as NVE ensemble) are conserved during the time integration, or

${\begin{array}{rcl}\mathbf {P} &\equiv &\sum _{i}m_{i}\mathbf {v} _{i}=\mathrm {Const.} \\\mathbf {L} &\equiv &\sum _{i}\mathbf {r} _{i}\times m_{i}\mathbf {v} _{i}=\mathrm {Const.} \end{array}},$

where P is linear momentum and L is angular mumentum. Usually, we subtract ceter-of-mass velocity from the velocity of each atom so that the whole system can not drift and the linear momentume becomes zero.

${\begin{array}{rcl}\mathbf {v} _{i}&:=&\mathbf {v} _{i}-\mathbf {v} _{\mathrm {CM} }\\\mathbf {P} &=&\sum _{i}m_{i}\mathbf {v} _{i}-\mathbf {v} _{\mathrm {CM} }\sum _{i}m_{i}=0\end{array}},$

where v_CM is defined as

$\mathbf {v} _{\mathrm {CM} }\equiv {\frac {\sum _{i}m_{i}\mathbf {v} _{i}}{\sum _{i}m_{i}}}$.

Similarly, we can additionaly subtract velocity component contributing rotation so that the whole system can not rotate and the angular momentum becomes zero as

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathbf{v}_i := \mathbf{v}_i - \mathbf{\Omega}\times\mathbf{r}_i } .

The angular velocity Ω is obtained from

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathbf{\Omega} \equiv \mathbf{I}^{-1} \mathbf{L} }

and moment-of-inertia matrix I is expressed as

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle \mathbf{I} \equiv \sum_i m_i(\mathbf{E}_3|\mathbf{r}_i|^2 - \mathbf{r}_i\otimes\mathbf{r}_i)} ,

where E₃ is a 3 by 3 identity matrix.

---

Let's see whether MD simulation results conserve both linear and angular momenta. For the test simulation, a diamomd cubic structure of 216 silicon atoms are equilbriated under periodic boundary conditions for 10,000 steps using NVE ensemble.

	Figure 1
(a)	(b)	(c)

The linear momentum is conserved within the machine precision in figure 1(b), while the angular momentum is not as shown in figure 1(c). According to Allen and Tildesley, the reason is told as

If the system is invariant to rotation about an axis, then the corresponding angular momentum component is conserved. ... In the case of the periodic boundary conditions, ... none of them were spherically symmetrical; in fact it is impossible (in Euclidean space) to construct a spherically symmetric periodic system. Hence, ... total angular momentum is not conserved in most molecular dynamics simulations.

Basically, what it says is that in most MD simulations angular momentum is no longer conserved if periodic boudary conditions are used and atoms move across the simulation cell boundary.

Let's check the above statement via another test simulation, where a Si cluster of 216 atoms are equilbrated with the same initial velocity. In this simulation, the simulation box has margin surrounding the Si cluster such that the cluster no longer sees its periodic images.

	Figure 2
(a)	(b)	(c)

In figure 2(b) and 2(c), we can see that linear and angular momenta are conserved within the machine precision in this test simulation.

---

Since the angular momentum is not conserved in MD simulations under periodic boundary conditions, it is usually unnecessary to make angular momentum be zero when the velocities of atoms are initialized. However, when we perform simulations of nanowires or rods, it is desirable to prevent the rotation around wire axis. For example, if a nanowire is aligned along the z-direction, we want to let the z-component of angular momentum be zero, or L_z = 0.