M10 Angular momentum is conserved or not

Conservation of Angular Momentum

Keonwook Kang, Seunghwa Ryu and Wei Cai

Jul 06 , 2009

---

Linear and angular momentum conservations

From classical dynamics, if a system is isolated such that there are no external force (F_ext) and torque (T_ext) acting on it, its linear and angular momenta are conserved since

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 \begin{array}{rcl} \dot\mathbf{P} & = & \mathbf{F}_{\mathrm{ext}} = 0 \\ \dot\mathbf{L} & = & \mathbf{T}_{\mathrm{ext}} = 0, \end{array} }

where linear momentum P and angular momentum L are obtained 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 \begin{array}{rcl} \mathbf{P} & \equiv & \sum_i m_i\mathbf{v}_ i \\ \mathbf{L} & \equiv & \sum_i \mathbf{r}_i \times m_i\mathbf{v}_i \end{array} }

in multi-body system.

Since molecular dynamics simulations follows classical dynamics theory, we would naturally expect that linear and angular momenta of a system of NVE ensemble will be constant during the time integration. Or

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 \begin{array}{rcl} \sum_i m_i\mathbf{v}_ i = \mathrm{Const.} \\ \sum_i \mathbf{r}_i \times m_i\mathbf{v}_i = \mathrm{Const.} \end{array}. }

In this manual, we will check whether linear and angular momenta are indeed conserved in MD simulations with or without periodic boundary conditions.

Zeroing out linear/angular momentum

In molecular dynamics simulations, we usually prevent our system of interest from drifting. For this purpose, we subtract ceter-of-mass velocity from the velocity of each atom.

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{v}_{\mathrm{CM}}, }

where v_CM is defined 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}_{\mathrm{CM}} \equiv \frac{\sum_i m_i \mathbf{v}_i}{\sum_i m_i} } .

In this case, the linear momentum becomes zero.

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{P} = \sum_i m_i\mathbf{v}_i - \mathbf{v}_{\mathrm{CM}}\sum_i m_i= 0 }

This is a special case of linear momentum conservation. And it accords well with periodic boundary conditions (PBC). You can think of a simulation system that loses some amount of linear momentum that leaves through a boundary gains exactly the same amount of linear momentum from the opposite boundary due to PBCs.

Similarly, we can additionally subtract velocity components, $\Omega \times \mathbf {r} _{i}$ shown below, contributing rotation motion so that the whole system can not rotate.

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.

However, it is questionable that a system under periodic boundary conditions can indeed maintain constant angular momentum. Let's check the conservations of both linear and angular momentum in test MD simulations below.

---

MD simulations of NVE ensemble

Bulk Si

For the test simulation, a diamond cubic structure of 216 silicon atoms are equilbriated under periodic boundary conditions for 10,000 steps using NVE ensemble. Linear and angular momenta are initially zeroed out as explained above.

	Figure 1
(a)	(b)	(c)

During 10,000 steps, 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). What is interesting is that the angular momentum is fluctuating around zero. Even if the angular momentum is not initially zeroed out, its behavior is very much similar in terms of fluctuation maginutude and frequency, though it is not shown here. 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.

Si cluster

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.

---

Si NW

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.

In MD++ there is a variable zerorot, depending on whose value components of the angular momentum can be selletively zero. If zerorot = "x", x-component of angular momentum will be zero. If zerorot = "all", all three components of the angular momentum will be zero. The default value of zerorot is "none". The other possible values are "y" and "z".

In the next test simulation, a Si NW is aligned along the z-direction and we expect L_z = 0.

	Figure 3
(a)	(b)	(c)

In figure 3(b) and 3(c), we can see that all three components of linear momentum and z-component of angular momentum are conserved during equilibriation.

MD++ script for this simulation is given below.

# -*-shell-script-*-
# Check linear momentum and angular momentum conservation 

#*******************************************
# Definition of procedures
#*******************************************
proc initmd { n } {
#MD++ setnolog
MD++ setoverwrite
MD++ dirname = "runs/si-test-$n"
MD++ NNM = 300
}
#end of proc initmd

proc create_crystal { N } {
#--------------------------------------------
# Create Perfect Lattice Configuration
#
MD++ latticestructure = diamond-cubic
MD++ latticeconst = 5.431 #(A) for Si
MD++ latticesize  = \[ 1 0 0  $N \
                       0 1 0  $N \
                       0 0 1  $N \]
MD++ makecrystal finalcnfile = perf.cn writecn
}

#--------------------------------------------
proc setup_md { } { MD++ {
equilsteps = 0  timestep = 0.001 # (ps)
atommass = 28.0855 # (g/mol) Si
DOUBLE_T = 0  equilsteps = 0
saveprop = 1 savepropfreq = 10 openpropfile #run
savecn = 0 savecnfreq = 10000
plotfreq = 20
vt2 = 2e28  #1e28 2e28 5e28
wallmass = 2e3     # atommass * NP = 14380
boxdamp = 1e-3     # optimal damping for 216 atoms and wallmass 1e-3
randseed = 12345 srand48
#MD++  srand48bytime
saveH # Use current H as reference (H0), needed for specifying stress
fixboxvec  = [ 0 1 1
               1 0 1
               1 1 0 ]
stress = [ 0 0 0
           0 0 0
           0 0 0 ]
output_fmt = "curstep EPOT Tinst VIRIAL(0) VIRIAL(4) DVIRIALDexx(0) DVIRIALDexx(4)"
writeall = 1
} }

#*******************************************
# Main program starts here
#*******************************************
set T 300
initmd NW
create_crystal 5
MD++ input = \[ 3 2 0 0 5 0 \] makecylinder

MD++ input = \[1 1 1\] changeH_keepR
MD++ input = \[2 2 1\] changeH_keepR

setup_md
MD++ totalsteps = 10000
MD++ T_OBJ = $T DOUBLE_T = 1 timestep = 0.001

MD++ initvelocity_type = "Gaussian" zerorot = "z" initvelocity
MD++ finalcnfile = init.cn writecn
MD++ finalcnfile = init.cfg writeatomeyecfg
  
MD++ ensemble_type = "NVE" integrator_type = "VVerlet"
MD++ {output_fmt = "curstep EPOT KATOM Tinst HELMP \
              TSTRESSinMPa_xx TSTRESSinMPa_yy TSTRESSinMPa_zz \
              TSTRESSinMPa_xy TSTRESSinMPa_xz TSTRESSinMPa_yz \
              H_11 H_22 H_33 H_12 H_13 H_21 H_23 H_31 H_32 \
              [MD++_Get P_com(0)] [MD++_Get P_com(1)] [MD++_Get P_com(2)] \
              [MD++_Get L_com(0)] [MD++_Get L_com(1)] [MD++_Get L_com(2)]" }
MD++ run

MD++ finalcnfile = equil.cn writecn
MD++ finalcnfile = equil.cfg writeatomeyecfg
MD++ quit

The MD++ variable zerorot can be used with the MD++ command initvelocity as shown above. Or it can be used with another command zerorotation. Internally, initvelocity calls the function zero_rotation() depending on the value of zerorot.

---

MD simulations of NVT ensemble

Upto here we have considered MD simulations of NVE ensemble. People might wonder whether the linear and the angular momenta will be conserved in case of canonical or NVT ensemble. In MD++, Nosé-Hoover thermostat is used for NVT ensemble. Nosé's Hamiltonian is expressed in time scaling variable s 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 \mathcal{H} = \frac{1}{2} \sum_i \frac{|\mathbf{p}_i|^2}{m_i s^2} + U(\{\mathbf{r}_i\}) + \frac{p_s^2}{2Q} + \frac{g\ln s}{\beta}} ,

where r_i and p_i are scaled coordinate and momentum of atoms. And they are related with real coordinate r_iFailed 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 '} and momenutm p_iFailed 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 '} 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 \begin{array}{rcl} \mathbf{r}_i & = & \mathbf{r}'_ i \\ \mathbf{p}_i & = & s\mathbf{p}'_i \end{array}} .

The real time step 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 \Delta t'} are related with the scaled time step 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 \Delta t} 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 \Delta t ' = \Delta t / s } .

Q is thermal mass, g is number of degree of freedom, and p_s is the conjugate momentum of s. Then, the time derivative of the linear and the angular momenta in real coordinate and momentum can be 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 \begin{array}{rcl} \frac{d}{dt'}\sum_i\mathbf{p}'_i & = & \frac{1}{s}\frac{d}{dt'}\sum_i\mathbf{p}_i - \frac{1}{s^2}\frac{ds}{dt'}\sum_i\mathbf{p}_i = \frac{d}{dt}\sum_i\mathbf{p}_i - \frac{1}{s}\frac{ds}{dt}\sum_i\mathbf{p}_i \\ \frac{d}{dt'}\sum_i\mathbf{r}'_i\times\mathbf{p}'_i & = & \frac{1}{s}\frac{d}{dt'}\sum_i\mathbf{r}_i \times \mathbf{p}_i - \frac{1}{s^2}\frac{ds}{dt'}\sum_i\mathbf{r}_i \times \mathbf{p}_i = \frac{d}{dt}\sum_i\mathbf{r}_i \times \mathbf{p}_i - \frac{1}{s}\frac{ds}{dt}\sum_i\mathbf{r}_i \times \mathbf{p}_i \end{array}} .

The 1st terms on the right hand side are zeros, since the linear and the angular momentum are conserved in scaled formulation. However, s is varing as a function of time and 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 ds/dt \neq 0} . Thus the linear and angular momenta in real varaiable formulation would not be constant in time, unless the following is satisfied.

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 \begin{array}{rcl} \sum_i\mathbf{p}_i & = & 0 \\ \sum_i \mathbf{r}_i \times \mathbf{p}_i & = & 0 \end{array}} .

In other words, only if the linear and angular momenta are intially zeros, those quantites are conserved in NVT ensemble. Besides, it is known that Nosé thermostat will generate the trajectory of canonical ensemble if the total momentum is zero.

Figure 4 shows the linear and the angular momenta as functions of integration step in MD simulation of NVT ensemble when both linear and angular momenta are initially set to be zeros. You can see that both momenta are conserved.

	Figure 4
(a)	(b)	(c)

On the other hand, figure 5(c) shows how the angular momentum behaves in NVT ensemble simulation when it is not zero initially.

	Figure 5
(a)	(b)	(c)

According to figure 5(c), angular momentum exponentially decays in the early time and fluctuates around some constant value later. The equations of motion of Nosé-Hoover thermostat are given 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 \begin{array}{rcl} \frac{d \mathbf{p}'_i}{dt'} &=& - \frac{\partial U(\left\{ \mathbf{r}'_i \right\})}{\partial \mathbf{r}'_i} - \zeta \mathbf{p}'_i \\ \frac{d \mathbf{r}'_i}{dt'} &=& \frac{\mathbf{p}'_i}{m_i} \\ \dot{\zeta}' &=& \frac{gk_B}{Q} \left( T - T_\mathsf{obj} \right) \end{array}} .

From the 1st equation, we can get

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{r}'_i \times \frac{d \mathbf{p}'_i}{dt'} = - \mathbf{r}'_i \times \frac{\partial U}{\partial \mathbf{r}'_i} - \zeta \mathbf{r}'_i \times \mathbf{p}'_i } .

Rearranging time derivative term gives

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 \frac{d}{dt'} \left(\mathbf{r}'_i \times \mathbf{p}'_i\right) - \frac{d\mathbf{r}'_i} {dt'} \times \mathbf{p}'_i= - \mathbf{r}'_i \times \frac{\partial U}{\partial \mathbf{r}'_i} - \zeta \mathbf{r}'_i \times \mathbf{p}'_i } .

By substituting the 2nd equation for the 2nd term in the left hand side, the above equation becomes

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 \frac{d}{dt'} \left(\mathbf{r}'_i \times \mathbf{p}'_i\right) - \frac{\mathbf{p}'_i} {m_i} \times \mathbf{p}'_i= - \mathbf{r}'_i \times \frac{\partial U}{\partial \mathbf{r}'_i} - \zeta \mathbf{r}'_i \times \mathbf{p}'_i } ,

and the 2nd term in the left hand side becomes zero. If we take summation over i of remaining terms, we get

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 \sum_i \frac{d}{dt'} \left(\mathbf{r}'_i \times \mathbf{p}'_i\right) = - \sum_i \mathbf{r}'_i \times \frac{\partial U}{\partial \mathbf{r}'_i} - \zeta \sum_i \mathbf{r}'_i \times \mathbf{p}'_i }

Since the 1st term in the right-hand side is the summation of the internal torque, it is zero. Finally, the equation for the angular momentum becomes

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 \frac{d}{dt'} \mathbf{L}' = - \zeta \mathbf{L}'} .

This relation explains the exponential behavior at the early time. Later, when the system reaches equilibrium state, temperature will fluctuate around the objective temperature and 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 d\zeta/dt'} keeps changing its sign to compensate the difference, T-T_OBJ. Thus, 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 \zeta} will fluctuate around zero in the long run and that explains the flat region of angular momentum curve. Linear momentum would show similiar behavior if its initial value is not zero.